Matt Barker
email: m.barker17@imperial.ac.uk

Pierre Degond

Ralf Martin
email: r.martin@imperial.ac.uk

AND Mirabelle Mu Ûls

A MEAN FIELD GAME MODEL OF FIRM-LEVEL INNOVATION

Keywords: AMS subject classifications. 35Q89, 91B69, 91A16, 49N80 mean field games, knowledge spillovers, innovation model

Knowledge spillovers occur when a firm researches a new technology and that technology is adapted or adopted by another firm, resulting in a social value of the technology that is larger than the initially predicted private value. As a result, firms systematically under-invest in research compared with the socially optimal investment strategy. Understanding the level of underinvestment, as well as policies to correct it, is an area of active economic research. In this paper, we develop a new model of spillovers, taking inspiration from the available microeconomic data. We prove existence and uniqueness of solutions to the model, and we conduct some initial simulations to understand how indirect spillovers contribute to the productivity of a sector.

1. Introduction. When a business invests in research and development (R&D), such strategy only takes into account how a potential innovation may increase the investing company's private value. However, other businesses may utilise innovations made by the original investing company to increase their own profits. This is known in economic literature as the knowledge spillover effect. By only considering its private return, businesses systematically undervalue their own innovations and hence underinvest in R&D, compared with the socially optimal investment level. To counteract the under-investment, governments introduce R&D subsidy policies for certain sectors of the economy. In order to effectively allocate such subsidies, it is therefore important to understand the extent of under-investment and how it varies between sectors.

To understand the spillover effect we develop a mean field game (MFG) model of firms distributed heterogeneously between sectors and according to their productivity level, taking into account their microscopic behaviour. From a microeconomic perspective, the size of knowledge spillovers can be inferred from the network of patent citations [START_REF] Fung | Are knowledge spillovers driving the convergence of productivity among firms?[END_REF]. When an industrial technology is developed, it often gets patented. As part of the patent any previous technology that has been used must be cited. This results in a network of patent citations, where each citation can be used as a proxy for a spillover from one technology to another, so spillover sizes can be evaluated [START_REF] Dechezleprêtre | Knowledge spillovers from clean and dirty technologies: a patent citation analysis[END_REF].

In the model we develop, sectors are connected by a graph that is informed by and can be calibrated to the microeconomic network of patent citation data.

A first model of knowledge spillovers, by Cohen and Levinthal [START_REF] Cohen | Innovation and learning: the two faces of R & D[END_REF], considered the stock of knowledge of a firm to depend on the amount of investment in R&D of that firm and the total amount of investment by all other firms, through a mean field-type interaction. Only an initial analysis of the model was conducted in [START_REF] Cohen | Innovation and learning: the two faces of R & D[END_REF]. A later model, acknowledged in Section 13.2 of [START_REF] Acemoglu | Introduction to Modern Economic Growth[END_REF], started from a macroscopic perspec-tive, hence only the aggregate knowledge of the entire economy was considered and spillovers were assumed to increase the aggregate uniformly. This does not explain how spillovers heterogeneously affect firms. Similar models have also been used to study entrepreneurship and intellectual property rights, such as in [START_REF] Acs | Patents, knowledge spillovers, and entrepreneurship[END_REF]. There has been particularly extensive research of knowledge spillovers in cross-country models.

In such models, a country's own output is aggregated and the knowledge level increases at a rate that depends on the leading country's knowledge level. However, simplifying assumptions are made that may affect their accuracy, such as in [START_REF] Eeckhout | Knowledge spillovers and inequality[END_REF], where interactions take place in a discrete time setting, or [START_REF] Aghion | The effect of financial development on convergence: Theory and evidence[END_REF][START_REF] Howitt | implementation and stagnation: a Schumpeterian theory of convergence clubs[END_REF] where the interactions between firms was described only through the evolution of aggregate quantities. In this paper, we use an MFG model to both increase the complexity of the description of firms, and to link firm-level evolution directly to microeconomic data for spillover sizes.

There have been several other papers focussing on MFG-type models of knowledge spillovers, see [START_REF] Burger | On a Boltzmann mean field model for knowledge growth[END_REF][START_REF] Lucas | Knowledge growth and the allocation of time[END_REF]. The Boltzmann model studied in both papers doeas not consider how innovation among firms evolves, nor did it incorporate the microeconomic data related to patent citations in its formulation. As a result, the model studied in this paper can give greater insight into firm-level dynamics.

In this paper, we analyse a stationary MFG model describing the spillover effect.

The MFG model describes the long-term behaviour of firms with full anticipation of the future. MFGs were described mathematically by Lasry and Lions [START_REF] Lasry | Jeux à champ moyen. I --Le cas stationnaire[END_REF][START_REF] Lasry | Mean field games[END_REF], and simultaneously by Huang, Caines and Malhamé [START_REF] Huang | Distributed multi-agent decision-making with partial observations: asymptotic Nash equilibria[END_REF] and they build on the work of Aumann and related authors on anonymous games [START_REF] Aumann | Markets with a continuum of traders[END_REF][START_REF] Schmeidler | Equilibrium points of nonatomic games[END_REF]. The novelties of the system we develop are, first, that the distribution dependence enters into the drift term rather than in the cost functional and, second, that we are considering more than one population of agents. Therefore our MFG model can be classed as a multipopulation MFG with a non-separable Hamiltonian. There has been some work in both multi-population MFGs (see [START_REF] Cirant | Multi-population mean field sames systems with Neumann boundary conditions[END_REF]) and MFG models with non-separable Hamiltonians (see [START_REF] Ambrose | Strong solutions for time-dependent mean field games with non-separable Hamiltonians[END_REF][START_REF] Ferreira | Existence of weak solutions to stationary mean-field games through variational inequalities[END_REF]). However, we are aware of no literature for models that display both characteristics, so although our model is one dimensional, its interest reaches beyond this setting. As a result of the novelty of our model, the techniques we use to prove existence and uniqueness are also novel. However, they rely heavily on the ability to write a stationary Fokker-Planck equation in the form of an exponential. This characteristic has previously been used in [START_REF] Barker | Comparing the best-reply strategy and meanfield games: The stationary case[END_REF] to prove existence and uniqueness in MFG and BRS models in a slightly different framework.

The paper is organised as follows. In Section 2, we develop the spillover model by describing firm behaviour at a microscopic level and formally deriving the mean field limit. In section 3, we describe the MFG problem and prove existence of solutions to it. We also show uniqueness of such solutions holds, provided the coupling strength between sectors is small enough. In Section 4, we provide some deeper insights into the effects of the modelling parameters, through numerical simulations. The first simulations show how parameters describing effects unrelated to spillovers (for example the discount factor, the noise level and the labour efficiency) change the MFG model.

Our second group of simulations demonstrate the effect of the spillover network on the model. The spillover network is a sector-level network that aggregates the patent citation network. We show that the effect of a spillover on any sector is a result of all paths to that sector in the associated network, and not just the immediate connections between sectors, which is contrary to the current economic state of the art. Finally, in Section 5, we briefly discuss future research prospects for the model, including how we intend to apply the model to economic questions relevant to R&D subsidy policy.

This manuscript is for review purposes only.

2. Model development.

The microscopic model.

Firms. Assume there are L sectors within the economy, and in sector there are N firms. We assume firm i in group has s , i,j links with firm j in group , where s , i,j is a random variable, taking a value s ∈ N with probability 1 N p( , , s). The i th firm in sector has a productivity level Z ,i ∈ Ω = (0, z), which increases as a result of employing labour h ,i or due to knowledge spillovers from firms that they are linked with. The productivity dynamics are also affected by noise with strength σ ∈ (0, ∞). As a result, Z ,i evolves according to the following SDE

dZ ,i (t) =   (h ,i (t)) γ + 1 N L =1 N j=1 s , i,j Z ,j (t)   dt + σdB ,i (t) (2.1a) L (Z ,i (0)) = m 0 , (2.1b)
where N = L =1 N , B ,i is an independent Brownian motion with reflection at boundaries 0 and z, and γ ∈ (0, 1) represents the inefficiency in converting one unit of labour to one unit of knowledge. In the initial condition (2.1b), L (Z ,i (0)) denotes the law of the random variable Z ,i (0) and m 0 is an initial distribution, which may be different for each sector. We assume firms produce a quantity of differentiated good at a rate q ,i according to the production function

(2.2) q ,i = Z ,i .
Each firm sells their product at a market-determined price r ,i and maximises their profit subject to the other firms' decisions. Each agent's profit functional is given by

(2.3) J ,i (h) = E ∞ 0 (r ,i (t)q ,i (t) -wh ,i (t)) e -ρt dt , where h = (h ,i ) N i=1 L =1
. The wage, w, and the discount rate, ρ, are given constants.

Consumers. Assume there is a representative consumer with preferences given by Q

= 1 N L =1 N i=1 q α ,i
, the Dixit-Stiglitz constant elasticity of substitution (CES) form, and with average income Y . The value α ∈ (0, 1) is related to the elasticity of substitution. The demand for each variety can be found by maximising Q under the budget constraint that average expenditure is equal to average income, i.e.

1 N L =1 N i=1 r ,i q ,i = Y . This gives (2.4) q ,i = Br 1 α-1 ,i , B = Y R α 1-α , R = 1 N L =1 N i=1 r α α-1 ,i α-1 α
.

For the purposes of the firm-level optimisation problem, we assume R is fixed, in that it can't be changed by any individual firm -this becomes true as N → ∞ for each . For the mathematical analysis we assume B to be a fixed constant, which simplifies matters and enhances the model's relevant features. Later in the numerical simulations B will be determined as the solution to a fixed point problem, which endogenises the price formation through the interaction between firms and consumers.

This manuscript is for review purposes only.

Firm profits revisited. Now, the profit functional (2.3) can be rewritten as

J ,i (h) = E ∞ 0 (Z ,i (t)) α B α-1
-wh ,i (t) e -ρt dt , using the consumer behaviour (2.4) and the production function (2.2).

2.2. Mean field limit. When there are large numbers of firms in each sector, the microscopic model developed in Section 2.1 can become intractable. Instead, we assume the number of firms in each sector, N , goes to infinity while N N → A for some A ∈ (0, 1), which represents the proportion of firms in sector . In order to derive the limiting mean field model, we first define the empirical distributions for each sector = 1, . . . , L by m N = 1 N N i=1 δ Z ,i , where δ Z ,i is a Dirac delta at the point Z ,i . We can then rewrite the dynamics (2.1) using m N as

dZ ,i (t) =   (h ,i (t)) γ + L =1 (N ) 2 N Ω z dm N (z , t) 1 N N j=1 s , i,j   dt + σdB ,i (t) L (Z ,i (0)) = m 0 ,
Assuming m N has a limit, m , as N → ∞ then, in the limiting model, a representative firm in sector evolves according to the SDE

dZ h,m (t) = (h (t)) γ + L =1 A p( , ) Ω z dm (z , t) dt + σdB (t) (2.5a) L Z h,m (0) = m 0 , (2.5b)
where, by the law of large numbers, p( , ) = ∞ s=0 p( , , s). The corresponding profit functional is

(2.6) J (h; m) = E   ∞ 0   Z h,m (t) α B α-1 -wh (t)   e -ρt dt   .
If all firms act in the same way as the representative firm, then the distribution of firms with respect to productivity level is given by a system of L Fokker-Planck equations

∂ t m = -∂ z (h ) γ + L =1 A p( , ) Ω z dm (z , t) m + σ 2 2 ∂ 2 zz m (2.7a) -(h ) γ + L =1 A p( , ) Ω z dm (z , t) m + σ 2 2 ∂ z m z=0,z = 0 (2.7b) m (z, 0) = m 0 (z) . (2.7c)
3. The MFG model. on Ω, such that for any other sequence of controls h and every

J (h * , m * ) ≥ J (h , m * ) (3.1a) and m * = L Z h * ,m * . (3.1b)
Such a distribution is called an MFG equilibrium.

To find an MFG equilibrium we first describe the Hamilton-Jacobi-Bellman (HJB) PDE related to the optimisation part of the problem (3.1a). Then we couple the HJB PDE to the Fokker-Planck PDE (2.7) to solve the consistency part (3.1b).

We start by defining L Hamiltonians

H : Ω × H 1 (Ω) L × R → R, for = 1, . . . L as (3.2) H (z, m, λ) = sup h≥0 h γ + L =1 A p( , ) Ω z m (z ) dz λ + z α B α-1 -wh = (1 -γ) γ w γ 1-γ max(0, λ) 1 1-γ + λ L =1 A p( , ) Ω z m (z ) dz + z α B α-1 ,
where z ∈ Ω is productivity, m = (m ) L =1 is a distribution of firms in each sector and λ is an adjoint variable. The optimal control is given by h * = γ w max(0, λ)

1 1-γ , for = 1, . . . , L.
Then we define the running profit V (z, t), for = 1, . . . , L, by

(3.3) V (z, t) = sup h E ∞ t (Z (s)) α B α-1 -wh (z) e -ρ(s-t) ds Z (t) = z ,
where Z (s) follows (2.5). If we let the equilibrium distribution be given by m (for = 1, . . . , L), then the MFG PDE system is stationary and given by

V ∈ H 1 (Ω) (3.4a) m ∈ H 1 (Ω) (3.4b) - σ 2 2 V + ρV -H (z, m, V ) = 0 (3.4c) - σ 2 2 m + (∂ λ H (z, m, V ) m ) = 0 (3.4d) V | z=0,z = 0 (3.4e) - σ 2 2 m + ∂ λ H (z, m, 0) m z=0,z = 0 (3.4f) Ω m (z) dz = 1 . (3.4g)
It can be shown, using either the dynamic programming principle (c.f [START_REF] Tembine | Mean-field-type games[END_REF]) or the stochastic maximum principle (c.f. [START_REF] Carmona | Probabilistic Theory of Mean Field Games with Applications I[END_REF]), that V (z), as defined by (3.3), satisfies the HJB equation (3.4c), (3.4e). The Fokker-Planck system (3.4d), (3.4f), (3.4g) comes from the distribution in the previous section (2.7) and the consistency condition (3.1b). This manuscript is for review purposes only.

Theorem 3.3. There exists a solution (m, V ) ∈ C 2 (Ω) ∩ C 1 Ω 2L to (3.4).

Furthermore, if

L =1 A p( , ) is small enough for every = 1, . . . L, then the solution is unique.

Proof outline. As noted in the introduction, this proof is based on the proof of existence and uniqueness in [START_REF] Barker | Comparing the best-reply strategy and meanfield games: The stationary case[END_REF]. The proof presented here has some technical differences compared with the one in [START_REF] Barker | Comparing the best-reply strategy and meanfield games: The stationary case[END_REF], hence it is reproduced in full. However, it follows a similar framework and so we do not claim the proof to be new. First, for k ∈ [0, ∞) we introduce an auxiliary system of PDEs defined by

V k ∈ H 1 (Ω) (3.5a) - σ 2 2 V k + ρV k -H k z, V k = 0 (3.5b) V k z=0,z = 0 , (3.5c) m k ∈ H 1 (Ω) (3.6a) - σ 2 2 m k + γ w max(0, V k ) γ 1-γ + k m k = 0 (3.6b) - σ 2 2 m k + km k z=0,z = 0 (3.6c) Ω m k (z) dz = 1 , (3.6d) where H k (z, λ) = (1 -γ) γ w γ 1-γ (max (0, λ)) 1 1-γ + kλ + z α B α-1 .
We use a modified version of upper and lower solutions (c.f. [START_REF] Schmitt | Boundary value problems for quasilinear second order elliptic equations[END_REF]) to prove existence and uniqueness of a weak solution V k to (3.5) for any k ∈ [0, ∞), and use elliptic regularity the-

ory to show V k ∈ C 2 (Ω) ∩ C 1 Ω . Next we define m k = 1 mk 1 mk , where mk = e 2 σ 2 kz+ z 0 ( γ w max(0,(V k ) )) γ 1-γ dy and mk 1 = Ω mk dz, for k ∈ [0, ∞). We prove that m k ∈ C 2 (Ω) ∩ C 1 Ω
and that m k is the unique solution of (3.6). Finally, we de- Solutions to the auxiliary HJB PDE.

fine a map Φ : [0, ∞) L → [0, ∞) L by Φ (k) = L =1 A p( , ) Ω zm k (z) dz for = 1
Theorem 3.4. There exists a unique solution V k ∈ C 2,τ Ω to the auxiliary HJB PDE (3.5) for any k ∈ [0, ∞) and some τ ∈ (0, 1), where C 2,τ Ω is the set of C 2 functions on Ω whose second derivative is Hölder continuous with exponent τ .

Furthermore, 0 ≤ V k ≤ zα ρB α-1
Proof. The existence part of the proof uses the theory of upper and lower solutions, specifically Theorem 4.3. in [START_REF] Loc | Bernstein-Nagumo conditions and solutions to nonlinear differential inequalities[END_REF], and follows along similar lines to the proof of Proposition 3.12 in [START_REF] Barker | Comparing the best-reply strategy and meanfield games: The stationary case[END_REF]. This shows that a solution V k ∈ W 1,p (Ω) to the auxiliary HJB PDE exists, for some p ≥ 1, provided the following hold true:

This manuscript is for review purposes only.

There exist constants

V ≤ V such that ρ V -z α B α-1 ≤ 0 ≤ ρ V -z α B α-1 , for every z ∈ Ω. 2. There exist constants a k ∈ R and b k > 0 such that ρu -(1 -γ) γ w γ 1-γ (max(0, λ)) 1 1-γ -kλ - z α B α-1 ≤ a k + b k |λ| p ,
for every z ∈ Ω, u ∈ V, V and every λ ∈ R.

If these two properties hold, then V ≤ V k ≤ V . The first assertion is true by taking V = 0 and V = zα ρB α-1 , which also gives the required bounds for V k . The second assertion is true with

b k = k + (1 -γ) γ w γ 1-γ , a k = 2z α B α-1 + b k , and p = 2 1-γ , as then ρu -(1 -γ) γ w γ 1-γ (max(0, λ)) 1 1-γ -kλ - z α B α-1 ≤ ρ|u| + k + (1 -γ) γ w γ 1-γ max 1, |λ| 1 1-γ + zα B α-1 ≤ a k + b k |λ| p . Now, since Ω is bounded and p > 2, V k ∈ H 1 (Ω). To show V k ∈ C 2,τ Ω , take any solution V k to (3.5) and define f = 2 σ 2 σ 2 2 -ρ V k + k V k + (1 -γ) γ w γ 1-γ max(0, V k 1 1-γ + z α B α-1 .
Then V k is a solution of -u +u = f , where f ∈ L 2 (Ω). So, from the elliptic regularity result of Proposition 7.2. p.404 in [START_REF] Taylor | Partial Differential Equations I[END_REF],

V k ∈ H 2 (Ω). Therefore V k ∈ H 1 (Ω)
, and so f ∈ H 1 (Ω) because α ∈ (0, 1). So, from the elliptic regularity result of Proposition 7.4. p.407 in [START_REF] Taylor | Partial Differential Equations I[END_REF], V k ∈ H 3 (Ω). Then, by the Sobolev inequality (c.f. Theorem 6 p.270 in [START_REF] Evans | Partial Differential Equations[END_REF])

V k ∈ C 2,τ Ω .
To prove uniqueness we use the strong maximum principle and Hopf's lemma, as stated in [START_REF] Evans | Partial Differential Equations[END_REF] Section 6.4.2. pp. 330-333. Suppose, for some k ∈ [0, ∞), there are

two solutions V 1 , V 2 ∈ C 2 (Ω) ∩ C 1 ( Ω) to (3.5) and V 1 = V 2 . If we define u = V 1 -V 2 ,
then u must attain its maximum at some point z * ∈ Ω. Suppose at this point u > 0.

Note that if this were not the case, we could consider the minimum, as either its maximum or its minimum must be non-zero. The argument for the minimum is the same as the one for the maximum, so it is omitted. First suppose z * ∈ Ω. Since this

is the maximal point, u (z * ) = 0, so V 1 (z * ) = V 2 (z * ).
Hence, there exists an open, connected and bounded region U such that U ⊂ Ω, z * ∈ U and

- σ 2 2 u = -ρu + ku + (1 -γ) γ w γ 1-γ max (0, V 1 ) 1 1-γ -max (0, V 2 ) 1 1-γ ≤ 0 ,
for every z ∈ U . So, by the strong maximum principle, u is constant in U . In particular, using (3.5b), u(z * ) = 0. But this is a contradiction. The only other case is z * ∈ ∂Ω and u(z) < u(z * ) for every z ∈ Ω. Then, ∂u ∂ν z * > 0 by Hopf's Lemma, but by (3.5c), ∂u ∂ν = ∂V1 ∂ν -∂V2 ∂ν = 0. This again leads to a contradiction. Therefore

V 1 = V 2 and solutions to (3.5) are unique for every k ∈ [0, ∞). Proposition 3.5. Fix k, k 1 , k 2 ∈ [0, ∞).
Then, the unique classical solution to the auxiliary HJB PDE (3.5), as found in Theorem 3.4, satisfies the following properties:

1. V k is an increasing function on Ω i.e. V k ≥ 0
This manuscript is for review purposes only.

2. V k > 0 for all z ∈ Ω 3. V k ∞ = sup z∈Ω V k (z) ≤ zα (1-γ)B α-1 1-γ w γ γ 4. V k1 -V k2 ∞ ≤ 1 ρ zα (1-γ)B α-1 1-γ w γ γ |k 1 -k 2 |
5. V k is strictly increasing with respect to k

6. V k1 -V k2 ∞ ≤ 4z σ 2 zα (1-γ)B α-1 1-γ w γ γ |k 1 -k 2 | 7. V k (0) > 0 > V k (z).
Proof. Property (1): Suppose, for a contradiction, there exists z ∈ Ω such that V k (z) < 0. First, by the boundary condition (3.5c), z ∈ Ω. So, by the boundary conditions and continuity of V k , there exists

z 0 , z 1 ∈ Ω with z 0 < z 1 , V k (z 0 ) = V k (z 1 ) = 0 and V k (z) ≤ 0 for all z ∈ (z 0 , z 1 ). Suppose that z 0 , z 1 ∈ Ω. Then V k (z 0 ) ≤ 0 ≤ V k (z 1 ) by construction of z 0 , z 1 and differentiability of V k . Furthermore, V k (z 0 ) > V k (z 1 ) because V k < 0 in (z 0 , z 1 ). So, using (3.5b) 0 = - σ 2 2 V k (z 1 ) -V k (z 0 ) + ρ V k (z 1 ) -V k (z 0 ) - 1 B α-1 (z α 1 -z α 0 ) < 0 .
This is a contradiction, so z 0 = 0 or z 1 = z. Assume z 0 = 0, we will again prove a contradiction (the other two cases of z 1 = z and both z 0 = 0, z 1 = z follow along similar arguments so their proofs are omitted). Since V k (0) = V k (z 1 ) = 0 and

V k (z) < 0 for all z ∈ (0, z 1 ) then, by continuity of V k , we can find 1 , δ 1 ∈ (0, z1 2 ) such that V k (z) ≤ 0 for all z ∈ (0, 1 ] and V k (z) ≥ 0 for all z ∈ [z 1 -δ 1 , z 1 )
. Furthermore, V k is strictly decreasing on (z 0 , z 1 ). So, using these two facts and continuity of V k there exists δ ∈ (0, δ 1 ] and ∈ (0, 1 ] such that

1. V k ( ) = V k (z 1 -δ) = min V k ( 1 ), V k (z 1 -δ 1 ) 2. V k ( ) > V k (z 1 -δ) 3. V k ( ) ≤ 0 ≤ V k (z 1 -δ). Then -σ 2 2 V k (z 1 -δ) -V k ( ) + ρ V k (z 1 -δ) -V k ( ) -(z1-δ) α -α B α-1 < 0,
which contradicts the fact that V k is a solution to (3.5). Therefore, V k ≥ 0 in Ω.

Property (2): From Property (1), we know V k ≥ 0. Now suppose, for a contradiction, there exists

z * ∈ Ω such that V k (z * ) = 0. Then V k (z * ) = 0, since it is a minimum of V k . So, by (3.5b), V k (z * ) = z α ρB α-1 and, since V k (z * ) < d dz z α ρB α-1 , there exists z 0 , z 1 ∈ Ω with z 0 < z * < z 1 such that 1. V k (z 0 ) = V k (z 1 ) 2. V k (z 0 ) > z α 0 ρB α-1 and V k (z 1 ) < z α 1 ρB α-1 3. V k (z 0 ) ≤ 0 ≤ V k (z 1 ). Then -σ 2 2 V k (z 1 ) -V k (z 0 ) + ρ V k (z 1 ) -V k (z 0 ) -1 B α-1 (z α 1 -z α 0 ) < 0,
which is a contradiction of (3.5b). Therefore, V k (z) > 0 for all z ∈ Ω.

Property (3): Since V k is continuous on Ω, V k ≥ 0 and V k (0) = V k (z) = 0, then V k must have a maximum that it attains at some point z * ∈ Ω. Furthermore, since V k is continuously differentiable in Ω, then V k (z * ) = 0. So,
This manuscript is for review purposes only.

using the bound on V k found in Theorem 3.4

0 ≤ V k (z) ≤ V k (z * ) = w γ 1-γ (1 -γ) γ γ 1-γ ρV k (z * ) -k V k (z * ) - (z * ) α B α-1 1-γ ≤ zα (1 -γ)B α-1 1-γ w γ γ . Property (4): Take k 1 , k 2 ∈ [0, ∞) such that k 1 < k 2 . First we show V k2 - V k1 ≥ 0, then we show V k2 -V k1 ≤ (V k 1 ) ∞ ρ (k 2 -k 1 )
, and we can conclude using Property (3). Let u 1 = V k2 -V k1 and assume, for a contradiction, there exists z ∈ Ω such that u 1 (z) < 0. Then, u 1 attains a minimum at z * ∈ Ω and u 1 (z * ) < 0. First suppose z * ∈ Ω, then u 1 (z * ) = 0 and from (3.5b)

- σ 2 2 u 1 (z * ) = -ρu 1 (z * ) + (k 2 -k 1 ) V k1 (z * ) > 0 ,
since V k1 ≥ 0 and u 1 < 0. Then, by continuity of u 1 , there exists an open bounded, connected Ω ⊂ Ω such that z * ∈ Ω and u 1 < 0 for all z ∈ Ω . So, by the strong maximum principle, u 1 is constant in Ω . In particular, u 1 = 0, which contradicts

u 1 < 0 for all z ∈ Ω . So, z * ∈ ∂Ω and u 1 (z) < u 1 (z * ) for all z ∈ Ω. However, from
Hopf's lemma u 1 (z * ) = 0, which contradicts (3.5c). So, we conclude that u 1 ≥ 0.

Now let u 2 = V k2 -V k1 -, with = (V k 1 ) ∞ ρ (k 2 -k 1 ) < ∞.
We assume, for a contradiction, there exists z ∈ Ω such that u 2 (z) > 0. Then u 2 attains a maximum at z * ∈ Ω and u 2 (z * ) > 0. First suppose z * ∈ Ω, then u 2 (z * ) = 0 and from (3.5b)

- σ 2 2 u 2 (z * ) = -ρu 2 (z * ) + (k 2 -k 1 ) V k1 (z * ) -ρ < 0 , since u 2 > 0 and ρ ≥ (k 1 -k 2 ) V k1 (z * ).
Then, by continuity of u 2 , there exists an open bounded, connected Ω ⊂ Ω such that z * ∈ Ω and u 2 > 0 for all z ∈ Ω .

So, by the strong maximum principle, u 2 is constant in Ω . In particular, u 2 = 0, which contradicts u 2 > 0 for all z ∈ Ω . So, z * ∈ ∂Ω and u 2 (z) > u 2 (z * ) for all z ∈ Ω. However, from Hopf's lemma u 2 (z * ) = 0, which contradicts (3.5c). So, we can conclude that u 2 ≤ 0.

Property (5): The proof of Property (4) shows V k is increasing with respect to k. Now suppose, for a contradiction, there exists

z * ∈ Ω such that k 1 < k 2 but V k1 (z * ) = V k2 (z * ). First, assume z * ∈ Ω and define u = V k2 -V k1 . Then, u(z * ) = 0, u (z * ) = 0 and u (z * ) ≥ 0, since z * is a minimum of u. Furthermore, from
Property (2), V k1 (z * ) > 0. Therefore, using (3.5b), we get the contradiction

0 = - σ 2 2 u + (k 1 -k 2 ) V k1 < 0 .
Hence, z * ∈ ∂Ω, so u (z * ) = 0 ,using (3.5c). But, u (z * ) = 0 by Hopf's lemma, which is a contradiction. So, V k is strictly increasing with respect to k.

Property (6): Fix k 1 , k 2 ∈ [0, ∞). Let u = V k1 -V k2 . Then, u satisfies σ 2 2 u = ρu-k 1 u +(k 2 -k 1 ) V k2 -(1-γ) γ w γ 1-γ V k1 1 1-γ -V k2 1 1-γ .
This manuscript is for review purposes only.

Suppose for z ∈ Ω, u (z) ≥ 0. Then, since u (0) = 0, there exists z 0 ∈ [0, z] such that u (y) ≥ 0 for all y ∈ [z 0 , z] and u (z 0 ) = 0. Therefore

0 ≤ u (z) = z z0 u (y) dy ≤ 2z σ 2 ρ||u|| ∞ + |k 2 -k 1 | V k2 ∞ ≤ 4z σ 2 zα (1 -γ)B α-1 1-γ w γ γ |k 2 -k 1 | .
We can similarly show that u (z) ≥ -4z

σ 2 zα (1-γ)B α-1 1-γ w γ γ |k 2 -k 1 | if u (z) ≤ 0.
Hence, V k is Lipschitz continuous with respect to k with the required constant.

Property [START_REF] Burger | On a Boltzmann mean field model for knowledge growth[END_REF]: First, we will show V k (0) > 0 and then we will show V k (z) < 0. Both steps use a similar method. Note that V k (0) ≥ 0 and

V k (z) ≤ 0, because V k (0) = V k (z) = 0 and V k (z) > 0 for all z ∈ (0, z). So suppose, for a contradiction, that V k (0) = 0. Then, since V k ∈ C 2,τ Ω , we can use continuity of V k , V k , V k and (3.5b), (3.5c) to show V k (0) = 0.
We can also use continuity of V k to show that for every C > 0 there exists 1 > 0 such that

z ∈ (0, 1 ) =⇒ V k (z) < C. Therefore, for any z ∈ (0, 1 ) (3.7) V k (z) = z 0 y 0 V k (y ) dy dy -z V k (0) -V k (0) < C 2 z 2 .
Since V k > 0 in Ω, there exists 2 > 0 such that V k increases on (0, 2 ). Therefore, V k ≥ 0 on (0, 2 ). Take C = 2 ρB α-1 and = min( 1 , 2 , 1). Then, from (3.7)

(3.8) V k (z) < z 2 ρB α-1 ≤ z α ρB α-1 ,
for all z ∈ (0, ). But, by rearranging (3.5b) and using V k , V k ≥ 0, we get

V k (z) ≥ z α ρB α-1 , which contradicts (3.8). Hence, V k (0) > 0. Now suppose, for a contradiction, that V k (z) = 0. Then, since V k ∈ C 2,τ Ω , we can use continuity of V k , V k , V k and (3.5b), (3.5c) to show V k (z) = zα ρB α-1 .
We can also use continuity of V k to show that for every C > 0 there exists 1 > 0 such that

z ∈ (z -1 , z) =⇒ V k (z) > -C. Therefore, for any z ∈ (z -1 , z) (3.9) V k (z) = z z z y V k (y) dy+V k (z)-(z -z) V k (z) > zα ρB α-1 - C 2 (z-z) 2 . Since V k > 0 in Ω, there exists 2 > 0 such that V k decreases on (z -2 , z).
Therefore, V k ≤ 0 on (z -2 , z). Take C > 0 such that

kC ρ + C 2 + C 1 1-γ ρ (1 -γ) γ w γ 1-γ ≤ α zα-1 ρB α-1 ,
and = min( 1 , 2 , 1). Then, from (3.9)

(3.10) V k (z) > zα ρB α-1 - C 2 (z -z) 2 ≥ zα ρB α-1 - C 2 (z -z) ,
This manuscript is for review purposes only.

for all z ∈ (z -, z). But, from (3.5b) and using

V k ≤ C(z -z), V k ≤ 0, we get V k (z) = 1 ρ σ 2 2 V k (z) + k V k (z) + (1 -γ) γ w γ 1-γ V k 1 1-γ + z α B α-1 ≤ kC ρ (z -z) + C 1 1-γ ρ (1 -γ) γ w γ 1-γ (z -z) 1 1-γ + z α ρB α-1 ≤ z α ρB α-1 + kC ρ + C 1 1-γ ρ (1 -γ) γ w γ 1-γ (z -z) .
So, using (3.10), we get α zα-1

ρB α-1 (z -z) < kC ρ + C 1 1-γ ρ (1 -γ) γ w γ 1-γ (z -z), which contradicts the definition of C. Hence, V k (z) < 0.
The auxiliary Fokker-Planck equation.

Definition 3.6. Fix k ∈ [0, ∞) and let V k ∈ C 2 (Ω) ∩ C 1 Ω
denote the unique solution to (3.5). Then, we define the function m k : Ω → (0, ∞) by

(3.11a) mk = e 2 σ 2 kz+ z 0 ( γ w (V k ) ) γ 1-γ dy (3.11b) mk 1 = Ω mk dz (3.11c) m k = 1 mk 1 mk . Proposition 3.7. For every k ∈ [0, ∞), m k ∈ C 2 (Ω) ∩ C 1 Ω where m k is defined by (3.6).
Proof. First, note that m k is well defined because V k ≥ 0 and V k is uniformly bounded. Hence, there exists

C ∈ (1, ∞) such that mk (z) ∈ [1, C] and mk 1 ∈ [z, C z], so m k (z) ∈ 1 C z , C z . Furthermore, m k ∈ C Ω because V k ∈ C 1 Ω . Now, if m k ∈ C 2 (Ω) ∩ C 1 Ω , then its derivatives would be (3.12a) m k = 2 σ 2 k + γ w V k γ 1-γ m k (3.12b) m k = 2 σ 2 k + γ w V k γ 1-γ m k + 2γ 2 σ 2 w(1 -γ) γ w V k 2γ-1 1-γ V k m k . But, since V k ∈ C 1 Ω and m k ∈ C Ω , then 2 σ 2 k + γ w V k γ 1-γ m k is well- defined and continuous for all z ∈ Ω. Hence, m k ∈ C 1 Ω . Then, m k , V k and V k are continuous in Ω and from Proposition 3.5 V k > 0 in Ω. Hence, m k is well-defined in Ω, m k ∈ C(Ω) and m k ∈ C 2 (Ω) ∩ C 1 Ω .
This manuscript is for review purposes only. Proof. Take m k defined in Definition 3.6. Then, m k ∈ C 2 (Ω) ∩ C 1 Ω by Proposition 3.7. Furthermore, from (3.12), m k satisfies (3.6b), (3.6c). Finally, by construction, m k satisfies (3.6d). Therefore, a solution to the auxiliary Fokker-Planck equation (3.6) exists, it is given by m k , and m k ∈ C 2 (Ω) ∩ C 1 Ω . To prove uniqueness we follow the same proof in [START_REF] Barker | Comparing the best-reply strategy and meanfield games: The stationary case[END_REF]. For brevity we only outline the argument here.

First, with mk defined as in (3.11), we can use regularity of mk from Proposition 3.7 to show (3.6) is equivalent to

m k , m k mk ∈ H 1 (Ω) (3.13a) mk m k mk = 0 (3.13b) mk m k mk z=0,z = 0 , Ω m k dz = 1 . (3.13c)
Then, by multiplying (3.13b) by m k mk , integrating over Ω and using integration by parts, the system (3.13) is equivalent to

(3.14) m k ∈ H 1 (Ω) , there exists Z > 0 such that m k = 1 Z mk , Ω m k dz = 1 .
From the previous results in this section, we have shown there exists a unique solution to (3.14) given by m k from Definition 3.6. Hence, existence and uniqueness of the auxiliary Fokker-Planck PDE follows from the equivalence between (3.6) and (3.14).

The fixed point problem. Proposition 3.10. The function Φ defined in Definition 3.9 is bounded. Furthermore, defining P as the L × L matrix with entries P , = p( , ) and A as the column vector (A 1 , . . . , A L ) T , then

0 ≤ Φ(k) 1 ≤ z P A 1 ,
where the 1-norm • 1 is defined as This manuscript is for review purposes only.

Φ (k) ≤ z L =1 A p( , ) z 0 m k (z) dz = z L =1 A p( , ). Therefore, 0 ≤ L =1 Φ (k) ≤ z L =1 L =1 A p( , ) = z P A 1 .
Theorem 3.12. The function Φ : [0, ζ] L → [0, ζ] L defined in Definition 3.9 is Lipschitz in the 1-norm on R L . The Lipschitz constant is given by C max =1,...,L A P , where P = L =1 p( , ) and C depends on P A 1 , but not explicitly on P or A.

Proof. First, fix k ∈ [0, ζ].
From Property (5) of Proposition 3.5, the continuity of V k , V k , V k with respect to z in Ω, and equations (3.5b), (3.5c), we find

V k (0) = 2ρ σ 2 V k (0) ≥ 2ρ σ 2 V 0 (0) = V 0 (0) > 0 ,
with the middle inequality an equality if and only if k = 0. Similarly, V k (z) ≤

V ζ (z) < 0 with the first inequality an equality if and only if k = 0. Moreover, V k is continuous with respect to k due to (3.5b) and continuity of V k , V k with respect to k, which was proven in Proposition 3.5. Therefore, there exists 1 , 2 ∈ (0, 1)

and

C 1 , C 2 > 0, independent of k, such that V k (z) = z 0 V k (y) dy ≥ z 0 V 0 (y) dy ≥ C 1 z , if z ∈ [0, 1 ] V k (z) = - z z V k (y) dy ≥ - z z V ζ (y) dy ≥ C 2 z , if z ∈ [z -2 , z] .
Furthermore, by continuity of V k with respect to k and compactness of [0, ζ], there exists

C 3 > 0 such that inf k∈[0,ζ] V k (z) ≥ C 3 if z ∈ [ 1 , z -2 ]. Note that C j for j = 1, 2, 3 are all independent of k ∈ [0, ζ]. Therefore, if γ ≤ 1 2 , for any k 1 , k 2 ∈ [0, ζ]: (3.15) z 0 min V k1 (z), V k2 (z) 2γ-1 1-γ dz ≤ 1 0 (C 1 z) 2γ-1 1-γ dz + z-2 1 C 2γ-1 1-γ 3 dz + z z-2 (C 2 (z -z)) 2γ-1 1-γ dz ≤ 1 -γ γ (C 2γ-1 1-γ 1 + C 2γ-1 1-γ 2 ) + C 2γ-1 1-γ 3 z , while, using Proposition 3.5, if γ ≥ 1 2 (3.16) z 0 max V k1 (z), V k2 (z) 2γ-1 1-γ dz ≤ z zα (1 -γ)B α-1 2γ-1 w γ γ(2γ-1) 1-γ
. Now, with the definition of mk in (3.11), for any

k 1 , k 2 ∈ [0, ζ] we have | mk1 -mk2 | = e 2 σ 2 k1z+ z 0 ( γ w (V k 1 ) (y)) γ 1-γ dy -e 2 σ 2 k2z+ z 0 ( γ w (V k 2 ) (y)) γ 1-γ dy = 2 σ 2 k1z+ z 0 ( γ w (V k 1 ) (y)) γ 1-γ dy k2z+ z 0 ( γ w (V k 2 ) (y)) γ 1-γ dy e 2 σ 2 u du .
This manuscript is for review purposes only.

Then, using the uniform bound on V k (y) with respect to k given by Proposition 3.5, we get

| mk1 -mk2 | ≤ 2 C1 σ 2 (k 1 -k 2 )z + z 0 γ w V k1 (y) γ 1-γ - γ w V k2 (y) γ 1-γ dy ≤ 2 C1 σ 2 |k 1 -k 2 |z + γ w γ 1-γ z 0 (V k 1 ) (y) (V k 2 ) (y)
| mk1 -mk2 | ≤ 2 C1 σ 2 |k 1 -k 2 |z + γ w γ 1-γ γ 1 -γ || V k1 -V k2 || ∞ z 0 max V k1 (y) 2γ-1 1-γ , V k2 (y) 2γ-1 1-γ dy ≤ 2 C1 σ 2 z + C2 |k 1 -k 2 | , where C2 = 4z σ 2 z (1-γ)B α-1 1-γ γ w γ 2 1-γ C 2γ-1 1-γ 1 + C 2γ-1 1-γ 2 + γ 1-γ C 2γ-1 1-γ 3 z , if γ < 1 2 . While C2 = γ 1-γ 4z 2 σ 2 wz α γ(1-γ)B α-1 γ , if γ ≥ 1 2 . Note that for any k ∈ [0, ζ], mk
z 0 z(m k1 -m k2 ) dz ≤ 1 mk1 1 z 0 z( mk1 -mk2 ) dz + z 0 z mk2 mk1 1 mk2 1 dz mk1 1 -mk2 1 ≤2z z 0 | mk1 -mk2 | dz ≤ 4 C1 z σ 2 z 0 (z + C2 ) dz|k 1 -k 2 | = 2 C1 z2 (z + 2 C2 ) σ 2 |k 1 -k 2 | := C |k 1 -k 2 | . Now take k (1) , k (2) ∈ [0, ζ] L . Define P = L =1 p( , )
, then recalling the definition of Φ given in Definition 3.9 and using (3.19)

(3.20) Φ(k (1) ) -Φ(k (2) ) 1 = L =1 L =1 A p( , ) z 0 z(m k (1) -m k (2) ) dz ≤ C L =1 A P k (1) -k (2) ≤ C max =1,...,L A P k (1) -k (2) 1 ,
which concludes the proof.

This manuscript is for review purposes only.

Theorem 3.13. For any given data, there exists a solution to the innovation MFG (3.4). Furthermore, if P A 1 is fixed, this solution is unique provided A P < 1 C for every = 1, . . . , L.

Proof. From Proposition 3.10 and Theorem 3.12, the function Φ :

[0, ζ] L → [0, ζ] L
is a continuous function from a convex compact subset of R L to itself. Therefore, by Brouwer's fixed point theorem, Φ has a fixed point. Furthermore, Theorem 3.12 shows that Φ is a Lipschitz function in • 1 . The Lipschitz constant is given by C max =1,...,L A P , where C depends on P A 1 but not directly on P or A .

Therefore, for fixed P A 1 , Φ is a contraction map provided A P < 1 C for every = 1, . . . , L, and in this case the fixed point is unique. Remark 3.14. In practical terms we can guarantee the condition A P < 1 C holds for every = 1, . . . , L provided L is large enough. This is because

L =1 A = 1. So,
for fixed P A 1 , when L is sufficiently large we can take A to be sufficiently small so that A P < 1 C 4. Numerical simulations.

4.1. Consumers. In the previous analysis, we assumed that consumers play a passive role in the model. In particular, the constant B has been fixed. However, in doing so we have not modelled the active nature of consumers in determining the price index R. To include this when implementing our numerical methods we return to (2.4) and we normalise economic output to Y = 1. Then, by rearranging (2.4) and using the production function q ,i = Z ,i , we get

B = 1 N L =1 N i=1 Z α ,i 1 α-1 . So, as
the number of firms in each sector goes to infinity, B =

L =1 A Ω z α m (z) dz 1 α-1 .
Note that this now needs to be solved as a fixed point, as m itself depends on B.

Simulations.

We computed simulations with synthetic data, using the numerical method outlined in Appendix A. From an economics perspective it is important to understand how the model affects the sector-level productivity. The purpose of the simulations is to provide initial insights into the role of the modelling parameters and of the network configuration.

Parameter effects. The MFG depends on the parameters σ, w, α, γ and ρ.

Recall that σ > 0 is the strength of noise in an individual's dynamics, w > 0 is the wage paid to employees, α ∈ (0, 1) is a parameter in the consumer optimisation problem which ensures convexity, and γ ∈ (0, 1) is the returns to labour i.e. the inefficiency in converting one unit of labour to one unit of knowledge, it also ensures convexity of the firm-level optimisation problem.

This manuscript is for review purposes only.

In order to separate the parameter effects from any effects caused by the sector network, we ran simulations with just a single sector. We fixed z = 2, A = A 1 = 1 and P = 0.1, where z is the maximum productivity level, A 1 is the proportion of firms in sector 1 and P is the strength of connection from sector 1 to itself. For baseline values, This manuscript is for review purposes only.

we took σ = 1, w = 1, ρ = 1, γ = 0.5 and α = 0.5. For each simulation, we varied one parameter while keeping all others at the baseline level. Figures 1a and1b show that the relationship between α and the distribution of firms is a complex one. There is some α * ∈ (0, 1) where the average productivity reaches a maximum, while on (0, α * ] average productivity is monotonically increasing, and on [α * , 1) average productivity is monotonically decreasing. Note that, for fixed productivity level and firm distribution, a firm's revenue is

r q = Z α 1 N L =1 A Ω z α m (z) dz -1
, which consists of a term that increases with respect to α multiplied by a term that decreases with respect to α. This results in a competing effect between α and a firm's revenue, which in turn affects a firm's return on investment, and therefore its level of investment in labour. Since labour investment has an increasing effect on average productivity, the competing terms in the revenue equation directly correspond to the behaviour exhibited in figure 1b.

Figures 1c and1d shows the effect of γ on the sector-level productivity. Figure 1d shows that as γ increases, the average productivity decreases. Since γ relates to the inefficiency of converting one unit of labour to one unit of productive work, it seems counter-intuitive at first that average productivity would be a decreasing function of γ. Recall that the optimal level of employment is given by h * = γ w max(0, V )

1 1-γ ,
which increases productivity at a rate (h * ) γ . Then, h * is increasing with respect to γ for fixed V if and only if V ≥ w γ e γ-1 γ and (h * ) γ is increasing if and only if

V ≥ w γ e γ-1
. Hence, the effect of γ on the average productivity depends on V and how it changes with respect to γ.

The effects of ρ and σ on the average productivity, shown in Figures 1e,1f and 2a, 2b respectively, show the same trend: average productivity decreases as each parameter increases. The size of ρ is the extent to which a firm discounts future profits. As ρ increases, firms care less about the future state of the system and so they are less willing to invest in labour; it is an investment whose effect is only on the future value of productivity. This results in reduced average productivity in the long run, which can be seen in Figure 1f. As σ increases, the randomness in productivity evolution of each firm increases. So, the impact of labour on productivity decreases with increasing σ, and this is reflected in Figure 2b.

Finally, Figures 2c and2d shows that average productivity also decreases with increasing wage, w. The wage rate increases the cost of labour. So, we can directly see that as the wage increases, the optimal level of employment, and hence the average productivity, decreases.

Spillover size effects. The sector-level network, encoded by the vertex weights

A for sector , and the edge weights p( , ) for a transfer of knowledge from sector to sector , is called the spillover network as it describes how knowledge and productivity spills over from one sector to another. A path in the spillover network is called a spillover path, or just spillover if there is no ambiguity. A path of length 1 from sector to sector is called a direct spillover, a path of length 2 or greater from sector to sector is called an indirect spillover, and in both cases sector is called the receiving sector and sector is the originating sector.

In almost all economic literature, only direct spillovers have been modelled and we are aware of no models that pay attention to the effect indirect spillovers have on economic productivity. In this subsection, we begin investigating how the productivity of a sector is affected by the structure of the spillover network, and in particular the effect of indirect spillovers on productivity. To undertake this investigation, we

This manuscript is for review purposes only. the distribution of firms. In the second simulations, we randomly generated spillover networks in models with three sectors and used the collected data to hypothesise a relationship between the average productivity of a sector and the size of spillovers (direct and indirect) it received. In the final simulations, we tested our hypothesis on more randomly generated spillover networks, this time for models with 10 sectors, which more closely resembles the number of sectors in the real economy. We showed that the hypothesis developed accurately describes the relationship between the spillover network and the average productivity of firms, moreover there was a 20% reduction in error when direct and indirect spillovers were taken into account, compared with when only direct spillovers were considered. Therefore, our conclusion from this preliminary investigation is that indirect spillovers have a significant effect on economic productivity in our model and they should not be ignored.

The networks in Figure 3 length 2, and in network 3 it has two direct spillovers. So, the difference in productivity in sector C between network 2 and network 1 will show the effect of an indirect spillover compared with having no spillover, and the difference between networks 3 and 2 will

show the effect of an indirect spillover compared with a direct spillover. The differences in density of sector C are plotted in Figure 4a. From the plots, it can be seen that the density of firms is larger at high productivity levels in network 2 compared with network 1 and the density is lower at low productivity levels. This means that the indirect spillover from sector A to sector C has a positive effect on sector C, skewing the distribution towards higher productivity levels. The same behaviour can be seen when we compare sector C in network 3 to network 2, however the effect is an order of magnitude larger. Therefore, although an indirect spillover path has some positive effect compared with no path at all, the effect is less strong than a direct spillover path.

In Figure 4b, sector D of networks four to six were modelled. For sector D: in network 4 there is one indirect spillover with path length 2; network 5 has one indirect spillover with path length 2 and one with path length 3; finally network 6 has an infinite number of indirect spillovers, one for every path length. We have plotted the difference in density of sector D between network 5 and network 4 and between networks 6 and 5. The difference between network 5 and network 4 shows the effect of an indirect spillover of length 3, while the difference between network 6 and network 5 shows the effect of indirect spillovers of all lengths greater than 3. For the difference between network 5 and network 4, the same qualitative result as the difference between network 2 and network 1, in Figure 4a, is observed. This suggests that having spillover paths of greater length do have positive impacts on productivity, but with reduced impact for increased path lengths. Interestingly, sector D in network 6 is less productive than sector D in sector 5. Further investigation showed that if B is fixed, rather than the solution of a fixed point problem, then the effect that more paths result in greater productivity returns, see figure 4c. The reason for this is not immediately obvious and warrants further study. Since the observed change is very small, it can't be ruled out that this result is an artefact from simplifications in the model.

In the second set of simulations, we took a closer look at how the spillover network structure affects the average productivity within each sector. Recall that if, given a network, we know the value of the fixed point, k * , of the function Φ defined in Definition 3.9. Then the average productivity in sector is Ω zm (z) dz = Ω zm k * (z) dz.

So, to understand the relationship between average productivity and the network, we first need to understand the relationship between Ω zm k (z) dz and k , for any k ≥ 0. Then, we also need to understand the relationship between k * and the L × L matrix S with entries defined by S , = A p( , ), because

k * = S Ω zm k * (z) dz L =1 .
In Figure 5 Fig. 3: Sector-level networks for simulations in Figures 4a and4b of three vertices, with A = 1/3 for all . We created a random network between the vertices by choosing a connection probability p, and making a directed edge between vertices with probability p. We then weighted each directed edge with a random weight, chosen from a uniform distribution on [0, 1]. We repeated this 100 times for each connection probability, and recorded both the size of direct spillovers to each sector and the value of the fixed point of Φ. Figures 6a and6b shows a scatter plot of k * , the th co-ordinate of the fixed point of Φ, against the sum of direct spillover strengths L =1 S , . In the simulations with a high connection probability, Figure 6a, there is a strong linear relationship between k * and L =1 S , . However, with low connection probabilities, Figure 6b, the simulations tend to follow one of two weaker linear relationships with the row sum.

To understand the relationships further, we can look at the equation that k

* ∈ [0, ∞) implicitly satisfies: k * = L =1 S , Ω zm k * dz,
(k) = Ω zm k dz (4.2) k * = f (0) L =1 S , .
We can see this linear relationship between k * and L =1 S , in Figures 6c and6d, where we have taken the simulated points in Figure 6b, and split the data into those points which have only direct spillovers and those that have indirect spillovers as well.

In Figure 6c, where sectors with only direct spillover paths are considered, the linear relationship described by (4.2) can be clearly seen.

To understand how the value of k * depends on the matrix S in the case of indirect spillovers, we can return to the definition of the spillover size and f (k). If we assume that f is approximately linear for sectors with indirect spillovers, i.e.

f (k) = f 0 + f 1 k,
This manuscript is for review purposes only.

k * = (S(f 0 1 + f 1 k * )) ,
where 1 is the vector of length L with ones in every entry. Using the identity (I +

f 1 S) -1 = ∞ n=0 f n 1 S n , we can rearrange (4.3) (4.4) k * = f 0 ∞ n=0 f n 1 S n+1 1 ,
which gives a way to estimate the value k * directly from the initial data. Therefore, combining estimates (4.1) and (4.4), we can estimate the value of average productivity from the matrix S by (4.5)

Ω zm (z) dz = z - b 0 f 0 L =1 ∞ n=0 f n 1 (S n+1 ) , b1 + b 2 .
This manuscript is for review purposes only. This manuscript is for review purposes only.

The relationship suggests that the average productivity depends on S n for every n i.e. on indirect spillovers of every path length. Moreover, if f 1 is small enough, the effect of a spillover path is decreasing by an order of magnitude for every increase in path length, which agrees with our initial simulations of networks 1-6.

In order to verify the hypothesis, in the final simulations we ran a regression to estimate the parameters f 0 , f 1 , b 0 , b 1 , b 2 and provide evidence that approximation (4.5) is accurate. We performed 1000 simulations on networks of ten vertices, with connection probability chosen randomly and uniformly distributed in [0, 1], with connection strength chosen randomly and uniformly in [0, 3], and with sector sizes A also randomly chosen. We ran a nonlinear regression, of the form (4.4), on sectors with indirect spillovers, to obtain optimal values of f 0 and f 1 . Then, using the optimal values of f 0 and f 1 we ran a second nonlinear regression, of the form (4.5), to find the optimal values of b 0 , b 1 and b 2 . Table 1 gives estimates for the parameters f i and b i .

We found that average productivity does behave approximately according to (4.5),

with table 1 suggesting a statistically significant result. Visually, this can be seen in Figure 7b, where we plotted (4.5) using the optimal values of f i and b i . We also computed estimates for the model This is reconfirmed by the 20% reduction in R-squared error when indirect spillover paths are included in the model. Therefore, indirect spillover paths can not be ignored as a factor determining a sector's productivity.

Conclusion and future research.

We have developed an MFG model of firm-level innovation from a microscopic formulation. The model can be calibrated to fit economic data of spillovers, so its economic validity can be verified. We have been able to prove existence of solutions and, under a smallness assumption on the data, uniqueness. We have investigated numerically how the modelling parameters and the spillover network affects the sector-level productivity, through the development of a simple algorithm that takes advantage of the structure of the proof of existence and uniqueness.

In future work, we hope to compare the MFG model with the socially optimal behaviour, as described by the mean field optimal control problem. We will also use patent-level data to calibrate and test the two models for their accuracy. We hope the comparison between the social optimum and the competitive equilibrium will suggest a method for implementing socially optimal subsidy policies for R&D.

Appendix A. Numerical Methods. The numerical method we designed to solve (3.4) is informed by the structure of the proof of existence and uniqueness. The method of proof relies on the contraction mapping theorem to find a fixed point of the map Φ, defined in Definition 3.9. We are also required to solve a fixed point problem to find the value of the parameter B. In light of this, our numerical method proceeds as follows, after choosing an initial guess k 0 ∈ [0, ∞) L , B 0 ∈ [0, ∞) and tolerances

δ 1 , δ 2 .
This manuscript is for review purposes only. , we compute the next iteration, V k i ,B i n+1 , using a Newton-Raphson method:

V k i ,B i n+1 = V k i ,B i n -dF V k i ,B i n -1 F V k i ,B i n
.

(e) Continue iteratively until F V

k i ,B i n 1
≤ δ 1 and define V i, = V This manuscript is for review purposes only.

3. 1 .

 1 Problem formulation. The MFG problem is related to the search for Nash equilibria in the optimisation of the profit functional (2.6) while agents evolve according to the dynamics (2.5). Definition 3.1. The MFG problem is to find a pair (h * , m * ), where h * = (h * ) L =1 is a sequence of controls and m * = (m * ) L =1 is a sequence of probability distributions This manuscript is for review purposes only.

3. 2 .

 2 Existence and uniqueness of solutions to the MFG. Definition 3.2. A solution to the innovation MFG model (3.4) is defined to be a tuple (m, V ) = (m 1 , . . . , m L , V 1 , . . . , V L ) such that m : Ω → (0, ∞), V : Ω → R satisfy (3.4) in the weak sense for each = 1, . . . , L.

  , . . . , L, and using the Brouwer fixed point theorem we prove there exists k ∈ [0, ∞) L such that Φ k = k. We use the contraction mapping theorem to prove uniqueness under certain smallness assumptions for the data. Then it follows, by replacing k with Φ k in (3.5) and(3.6), that m k, V k = m k1 , . . . , m kL , V k1 , . . . , V kL is a (unique) solution to (3.4) with the required regularity.

Theorem 3 . 8 .

 38 There exists a unique solution m k ∈ C 2 (Ω) ∩ C 1 Ω to the auxiliary Fokker-Planck PDE (3.6) for any k ∈ [0, ∞).

Definition 3 . 9 .

 39 Fix k = (k ) L =1 ∈ [0, ∞) L . For = 1, . . . , L, let V k be the unique solution to the auxiliary HJB PDE (3.5) with constant k , and let m k be the unique solution to the auxiliary Fokker-Planck PDE (3.6) with constant k . Then we define the function Φ : [0, ∞) L → [0, ∞) L by k (z) dz , = 1, . . . , L .

x 1 =L

 1 =1 |x | for any x ∈ R L . Remark 3.11. Due to this proposition, we can define ζ = z P A 1 and consider only the restriction of Φ to [0, ζ] L , which we will still denote by Φ for convenience. Proof. Take = 1, . . . , L. Then Φ (k) = L =1 A p( , ) z 0 zm k (z) dz ≥ 0, since p( , ) ≥ 0 and m k ≥ 0. Similarly, since m k is a probability distribution,

dz ≥ 1 ,

 1 as V k ≥ 0. So, for any k 1 , k 2 ∈ [0, ζ], using (3.17) and (3.18), we have(3.19)

Theorems 3 .

 3 4 and 3.8 proved existence and uniqueness of solutions to equations (3.5) and (3.6) respectively for any k ∈ [0, ζ]. Now, if k * is a fixed point of Φ then (m * , V * ) := m k * , V k * L =1 is a solution to (3.4), which can be seen by replacing k * with Φ (k * ) in (3.5), (3.6) for every = 1, . . . , L. Conversely, if (m * , V * ) is a solution to (3.5), (3.6), then clearly, by defining k * co-ordinate wise as k * = L =1 A p( , ) z 0 zm (z) dz, k * ∈ [0, ζ] L is a fixed point of Φ. Furthermore, by uniqueness of (3.5), (3.6), m k * , V k * = (m * , V * ). So, existence and uniqueness of solutions to the innovation MFG (3.4) is equivalent to existence and uniqueness of fixed points of Φ. Hence, there exists a solution to the innovation MFG. Furthermore, this solution is unique, provided A P < 1 C for every = 1, . . . , L.

Fig. 1 :

 1 Fig. 1: Simulations of MFG with varying α, γ, and ρ

( a )Fig. 2 :

 a2 Fig. 2: Simulations of MFG with varying σ and w

  provide insight into how indirect spillovers affect the distribution of firms, in comparison to direct spillovers. In network 1, sector C has one direct spillover, in network 2 it has one direct spillover and one indirect spillover of This manuscript is for review purposes only.

  , we have plotted Ω zm k (z) dz against k. The relationship appears to approximately follow (4.1)Ω zm k (z) dz = z -b 0 k b1 + b 2 ,for some b 0 , b 1 , b 2 > 0, as can be seen by the second line in Figure5. To understand the relationship between the fixed point of Φ and the matrix S, we considered networks This manuscript is for review purposes only.

( a )= 1 Fig. 4 : 3 then

 a143 Fig. 4: Simulations of MFG comparing distribution of firms in sectors C and D with respect to productivity in networks 1, 2 and 3, and networks 4, 5 and 6 from Figure 3

Fig. 5 : 1 Fig. 6 :

 516 Fig. 5: Plot of average productivity against size of k in auxiliary Fokker-Planck equation (3.6) and plot of y = z -(z-1.4) k 2 +1 for comparison

,

  which assumes average productivity depend on direct spillovers only, and plotted the result in Figure7b. Comparing plots 7a and 7b shows that the model (4.5), which includes the effects of indirect spillovers, provides a more accurate estimate for average productivity than model (4.6), which only accounts for the effect of direct spillovers.

1 .

 1 Given k i ∈ [0, ∞) L and B i ∈ [0, ∞), solve (3.5b), (3.5c) using the following method, based on a Newton-Raphson method in a Banach space.(a) Define F (v) = -σ 2 2 v + ρv -kv -(1 -γ) -z α B α-1 . We want to find zeros of F (v). (b) We define dF (v)(u) = -σ 2 2 u + ρu -ku -γ w v γ 1-γ u , which is theFréchet derivative of F .

( c )

 c Denote by V k i ,B i 0 the initial guess for the th component of the solution to (3.5b), (3.5c) with k = k i and B = B i . (d) Given V k i ,B i n

k i ,B i n 2 .

 2 Given V i , compute the solution to (3.6b) using(3.11) and denote it bym i 3. Define k i+1 = Φ k i and B i+1 = L =1 A Ω z α m (z) dz 1 1-α 4. If k i+1 -k i 1 + |B i+1 -B i | ≤ δ 2 thenstop the iteration process and define the MFG solution (m, V ) = m i , V i . Otherwise return to Step 1. Variable Coefficient estimate Standard error t stat p value f

( a )

 a Plot of average productivity against right hand side of (4.5), with optimal values for fi, bi (b) Plot of average productivity against right hand side of (4.6), with optimal values for f0, bi

Fig. 7 :

 7 Fig. 7: Plots of average productivity against (4.5) with models for k given by considering direct and indirect spillovers (4.4) or only direct spillovers (4.6)

Table 1 :

 1 2080 0 Table of regression results related to linear regression (4.5)

This work was supported by the Natural Environment Research Council [grant number NE/L002515/1] (p.degondPD acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grants no. EP/M006883/1 and by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048. PD's affiliation from 01/01/2021 is "Institut de . MB is supported by the Natural Environment Research Council


Qianqian Wang

Jie Han
email: jiehan@xjtu.edu.cn

Society organization, not pathogenic viruses, is the fundamental cause of pandemics

Qianqian Wang 1 • Jie Han 1 • Hong Chang 2 • Chaoqi Wang 1 • Eric Lichtfouse 3,4 In cases of emergencies such as pandemics, it is common sense-and somewhat too easy-to design a culprit or enemy, e.g., a pathogenic virus, to be fought and eradicated, thus lessening the human responsibility of the issue. This is the classical strategy of the occidental world where the most diseases and other problems are mainly tackled by treating the symptoms, thus largely overlooking the profound causes [START_REF] Lichtfouse | Climate change, society issues and sustainable agriculture[END_REF][START_REF] Lichtfouse | Society issues, painkiller solutions, dependence and sustainable agriculture[END_REF]. There is now mounting evidence that societal and environmental factors have largely contributed to the coronavirus disease 2019 pandemic . For instance, since pathogens are occurring widely in nature, one might hypothesize that the pandemic would never have happened without animal-to-human transmission, which is controlled by society organization. Would the pandemic have been so deadly and fast-developing worldwide under low-population levels, density, and travel habits (Fig. 1)? Probably not. In order to prevent future pandemics, which will most likely happen because we have eradicated only the symptomatic virus, we should adopt a more oriental, holistic strategy that considers that the good functioning of any system depends on the good interaction of all its elements, summarized in the term 'harmony.' Practically, all societal and environmental factors controlling the pandemic should be studied. Here, we discuss viral transmission, imperfect protection, population density, and urbanization to unveil societal factors controlling pandemics, and to propose societal changes, such as an urban exodus, that would help prevent future pandemics.

Viruses travel fast

The World Health Organization (WHO) declared the novel coronavirus disease (COVID-19) a pandemic on 13 March 2020. As of 15 October 2021, nearly 240 million confirmed cases of COVID-19, including nearly five million deaths, have been reported around the globe (WHO 2021a). Humanto-human transmission through respiratory droplets, aerosols, contact surfaces or fomites is the primary cause for the transmission of COVID-19 [START_REF] Chen | Role of indoor aerosols for COVID-19 viral transmission: a review[END_REF][START_REF] Han | Microfiber pillow as a potential harbor and environmental medium transmitting respiratory pathogens during the COVID-19 pandemic[END_REF][START_REF] He | Electrostatic fine particles emitted from laser printers as potential vectors for airborne transmission of COVID-19[END_REF]WHO 2020a). In addition, fecal transmission and contact with virus-contaminated frozen foods and packages can spread the coronavirus (Han and He 2021;Han et al. 2021).

Infectious agents, hosts, and environmental transmission are the prerequisites for the spread of an infectious disease in the population. Possible modes of transmission for SARS-CoV-2 include airborne, contact, fecal-oral, and others (CDC 2021a; WHO 2020b). All these transmission routes are dependent on human behavior. Through the airborne route, viral transmission occurs primarily through respiratory droplets and aerosols generated during coughing or sneezing by infected individuals [START_REF] Li | Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia[END_REF][START_REF] Zhou | Modes of transmission of SARS-CoV-2 and evidence for preventive behavioral interventions[END_REF]. Aerosols generally have smaller particle sizes, less than 5 μm, than the most respiratory droplets and can travel longer distances. A recent study reported that SARS-CoV-2 could remain viable in aerosols for up to three hours at 21-23 °C, with a reduction in infectious titer from 10 3.5 to 10 2.7 TCID 50 per liter of air titer (TCID: median tissue culture infectious dose, [START_REF] Van Doremalen | Aerosol and surface stability of SARS-CoV-2 as compared with SARS-CoV-1[END_REF]. Similarly, SARS-CoV-2 maintained infectivity when suspended led to 20-40% reduction in virus intake by the wearer, while wearing a surgical mask reduced approximately 50% [START_REF] Ueki | Effectiveness of face masks in preventing airborne transmission of SARS-CoV-2[END_REF]). N95 particulate respirators showed the highest efficacy, enabling approximately 80-90% of reduction in virus intake at a social distance of 50 cm. The protection is twofold. For infected individuals, wearing a cotton or surgical mask also blocked more than 50% of the viruses from their respiratory activities [START_REF] Ueki | Effectiveness of face masks in preventing airborne transmission of SARS-CoV-2[END_REF].

Sometimes masks may not provide adequate means of protection. For instance, in crowded places people are very close to each other, and the situation is often worse in enclosed spaces with poor ventilation [START_REF] He | Electrostatic fine particles emitted from laser printers as potential vectors for airborne transmission of COVID-19[END_REF][START_REF] Wang | Unprotected mothers and infants breastfeeding in public amenities during the COVID-19 pandemic[END_REF]. Expiratory droplets and aerosols can escape from masks due to imperfect sealing or fitting. Even when properly worn, medical face masks cannot completely block the transmission of respiratory virus droplets and aerosols [START_REF] Ueki | Effectiveness of face masks in preventing airborne transmission of SARS-CoV-2[END_REF]. Indeed, aerosols containing influenza viruses can easily go around the face mask and into the user's nose and mouth [START_REF] He | Effect of particle size on the performance of an N95 filtering facepiece respirator and a surgical mask at various breathing conditions[END_REF]. Moreover, when wearers are talking, 53% of small respiratory droplets (0.5-20 μm) are leaked from the top, 15% from the side, and 9% from the bottom of standard surgical masks [START_REF] Cappa | Expiratory aerosol particle escape from surgical masks due to imperfect sealing[END_REF]. This is particularly relevant in situations where people wearing face masks are in close proximity to one another Fig. 1 The spread of infectious agents among humans is exacerbated by high population densities in modern societies. It is also impractical to maintain adequate social distances in certain environments, where wearing masks alone may not provide adequate protection for in aerosols for up to 16 h at room temperature [START_REF] Fears | Comparative dynamic aerosol efficiencies of three emergent coronaviruses and the unusual persistence of SARS-CoV-2 in aerosol suspensions[END_REF].

Transmission can also occur by viruses deposited on inanimate surfaces, which could infect people through their eyes, noses, and mouths after touching by hands. Moreover, SARS-CoV-2 viruses were found to show long persistence under low-temperature conditions, showing similar traits to SARS-CoV, which could survive up to two years at -20 °C (WHO 2020c). Researchers at the United States Army Medical Research Institute of Infectious Diseases (USAM-RIID) reported that SARS-CoV-2 survived on swine skin for 14 days at 4 °C [START_REF] Harbourt | Modeling the stability of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) on skin, currency, and clothing[END_REF]. There were also reports on re-emergent outbreaks of COVID-19 caused by virus-contaminated frozen foods and food packages in areas previously achieving local eradication of COVID-19 (Han et al. 2021). Overall, transmission is controlled not only by biological and physical factors, but also by human behavior.

Imperfect masks

Masks and face coverings are key measures to mitigate human-to-human transmission. Wearing face masks or cloth face coverings can considerably reduce the spread of respiratory pathogens. For healthy adults, wearing a cotton mask reducing the spread of respiratory pathogens. The mass availability of public transport and increased globalization with travel and trade activities may act synergistically and allow rapid spread of infectious agents in the world population (Fig. 2). Such risks may be exacerbated by the recent emergence of the B.1.617.2 variant, also known as the 'Delta variant,' which has triggered a new wave of massive infections of COVID-19 [START_REF] Kupferschmidt | Delta variant triggers new phase in the pandemic[END_REF]. As of 5 October 2021, the Delta variant had spread to 192 countries out of 194 WHO member countries (WHO 2021b), and became the dominant viral strain in the newly confirmed cases in the United States (CDC 2021b). A recent study showed that the respiratory viral load of people infected with the Delta variant was 1260 times higher than that of those infected with the original SARS-CoV-2 [START_REF] Li | Viral infection and transmission in a large, well-traced outbreak caused by the SARS-CoV-2 delta variant[END_REF]). This makes it even more difficult to rely upon mask-wearing to prevent infection in places with inadequate social distancing and poor ventilation. Overall, imperfect mask protection and the proximity of people in public spaces have favored viral transmission.

City hotspots

Frequent social interactions are crucial for viral respiratory infections to develop into an epidemic or pandemic. Previous studies showed that seasonal infectious human diseases such as avian influenza are closely related to population aggregation and movement [START_REF] Maya | Seasonal patterns in human A (H5N1) virus infection: analysis of global cases[END_REF]. This is a typical example where 'we knew,' but no change in society organization has been done to prevent future pandemics. High population density, frequent movement, and situations of close contact appear as major factors controlling the virus spread (Fig. 1). For instance, many studies have noted that population density, which is generally high in cities, strongly influences the likelihood of a disease outbreak [START_REF] Connolly | Extended urbanisation and the spatialities of infectious disease: demographic change, infrastructure and governance[END_REF]. High population density generally leads to closer contact and more interactions among residents, making cities the likely 'hotspots' during the early spread of emerging infectious diseases. Indeed, urban areas were initially hit much harder than suburban and rural areas in the COVID-19 pandemic. In the USA, cases and deaths were heavily concentrated in the metropolitan areas of New York, New Orleans, Boston, and Detroit during the early months of the pandemic (PRC 2021). The Lombardy region, which has the largest population and the highest population density in Italy, reported the highest incidence of coronavirus infection of any sub-national region in the country, whereas greener Mediterranean regions such as Sardinia, Calabria, and Basilica showed lower mortality rates [START_REF] Cortright | Cities and coronavirus: some thoughts[END_REF][START_REF] Lorusso | Paraneoplastic neurological syndromes: study of prevalence in a province of the Lombardy region[END_REF][START_REF] Roviello | Less COVID-19 deaths in southern and insular Italy explained by forest bathing, Mediterranean environment, and antiviral plant volatile organic compounds[END_REF]. Overall, highly populated areas showed higher transmission rates and cases, suggested that the future society should decrease population density in urban areas and moving people 'back' to the countryside.

No safe island

Public transportation provides the means of mass human movement, which is an essential utility in the modern society but also favors the spread of infectious diseases. These efficient and extensive public transport networks help pathogens move rapidly between cities and countries, even continents. Such phenomena have been observed in previous human disease outbreaks, e.g., Ebola (Sharifi and Khavarian 2020). In the COVID-19 pandemic, the coronavirus outbreak hit densely populated urban areas first and hardest; then, quickly spread into the suburbs due to the busy connections in public transport [START_REF] Tabernese | America's biggest cities were already losing their allure[END_REF]. Official statistics released in August 2020 suggested that 10% of new infections in Greece were linked to tourism (Roache 2021). [START_REF] Zhang | Exploring the roles of high-speed train, air and coach services in the spread of COVID-19 in China[END_REF] found that the spread of COVID-19 in China was closely correlated to the routes and frequencies of domestic air, train, and coach services. Viruses reached more rural and exurban areas along with the human movements. Eventually, there are no 'safe islands' during a pandemic in the modern, interconnected human society.

Deadly closed contact

Close contact was defined as being within six feet of someone infected with the virus for a cumulative total of 15 min or more over a 24-h period (CDC 2021a). Situations of close contact are common in urban environments such as in elevators, cinemas, churches, shopping centers, and transportation hubs. Under those circumstances, it is often difficult to maintain the six-foot social distance [START_REF] Wang | Unprotected mothers and infants breastfeeding in public amenities during the COVID-19 pandemic[END_REF], especially during peak hours (Fig. 2). During COVID-19, governments and public health authorities closed many public spaces to minimize non-essential contacts between people, and implemented symptom screening and the mandatory requirements of mask wearing in public transportation vehicles and facilities.

One of the exceptions, however, is religions activities [START_REF] Pennisi | The changing geography of COVID-19 in the U.S[END_REF]. For instance, the Kumbh Mela festival drew tens of millions of pilgrims across India to bathe in sacred rivers with celebrations and gathering for up to 48 days, which contributed to the surging numbers of COVID-19 infections between March and April 2021. On 7 May 2021, India registered 414,188 new COVID-19 cases, the highest daily figure in any country or territory since the outset of the COVID-19 (WHO 2021c). Italy, one of the catholic countries and the initial epicenter of the COVID-19 pandemic in Europe, also reported mass infections related to church activities [START_REF] Mogi | The influence of social and economic ties to the spread of COVID-19 in Europe[END_REF].

Human migrations

Historically, the emergence of pandemics and epidemics has coincided with the establishment of large urban settlements [START_REF] Antoni | The 1918 "Spanish Flu" in Spain[END_REF][START_REF] Ellwanger | Control and prevention of infectious diseases from a One Health perspective[END_REF]. The nineteenth century witnessed dramatic changes in the dynamic of urbanization. In Europe, two distinct phases of urbanization and human migration have occurred. The first phase, which was driven by industrialization and began in the 1850s, was characterized by the surging increase of urban population. Some industrial sectors became almost completely mechanized in the early mid-nineteenth century (Britannica 2021). Since more economic production was located in urban areas, people started to migrate from the countryside, e.g., farmers, to work near industrial plants, which was the beginning of the 'rural exodus' [START_REF] Mercer | Rural migration to urban settings: educational and welfare problems[END_REF]. Driven by better education, healthcare, infrastructure, and living conditions in cities, a great number of individuals continued to migrate to populated industrial areas [START_REF] Mervyn | Factors influencing migration and population movements-part 1. Future directions international strategic analysis paper[END_REF]. Between 1805 and 1911, the proportion of the French population living in cities increased from 25 to 44%, while in Germany the proportion of urban population increased from 30 to 60% (DHR 2021).

Since the 1960s, the industrialized nations in Europe have witnessed the second phase of mass human migration, known as 'suburbanization.' Many people started to migrate from rural areas and city centers to suburban areas. While European central cities sustained losses-to the order of more than three million residents-in reality suburbanization increased the total number of people living in the megacities. From the mid-1960s to the beginning of the 2000s, the metropolitan areas of Europe with populations over one million have gained about 24 million residents in total [START_REF] Cox | Examining sprawl in Europe and America[END_REF]. This trend continued into recent decades, although it has slowed down significantly. In 2021, there were 11.1 and 9.4 million of people living in the metropolitan areas of Paris and London, compared with only 7.4 and 8.2 million in 1960, respectively (WPR 2021).

A positive feedback loop exists between population growth and public transport services connecting metropolitan centers and suburban areas. On one hand, the mass availability of public transport provides a convenient means of daily commuting for people to travel between suburbs and city centers. On the other hand, the growing population living in city suburbs increases the demand for public transport services. Expressways, train stations, ferries, and airports are built to meet the demand of more people using public transport for daily commuting and inner-city traveling. Overall, a long-term and dramatic increase in the population living in urban and suburban areas has increased the risk of epidemics or pandemics.

A painful lesson

Recent evidence has indicated that SARS-CoV-2 might have been quietly circulating in communities before the first known outbreak (He and Han 2021). For example, in Italy, researchers found traces of SARS-CoV-2 in municipal wastewater sampled at Milan and Turin on December 18, 2019 (La [START_REF] Rosa | SARS-CoV-2 has been circulating in northern Italy since December 2019: Evidence from environmental monitoring[END_REF], suggesting that COVID-19 had already spread in their local population. Similarly, [START_REF] Fongaro | The presence of SARS-CoV-2 RNA in human sewage in Santa Catarina, Brazil, November 2019[END_REF] detected the coronavirus in sewage samples collected on 27 November 2019 in Brazil. An earlier study in Spain even detected SARS-CoV-2 in archived sewage samples collected in Barcelona in as early as March 2019 (Chavarria-Miró et al. 2020), several months before the first outbreak was reported in Wuhan, China [START_REF] Ma | Coronavirus: China's first confirmed Covid-19 case traced back to November 17. South China morning post[END_REF]. Regardless of the political debates over the geographic origin of COVID-19, the virus has spread to people all around the world in a quick and highly efficient manner.

These findings also raise the question that whether there will be other novel infectious agents, including zoonotic pathogens, circulating in the human population while they mutate, evolve, and eventually become epidemic or pandemic-causing agents. Although the answer is not immediately clear for SARS-CoV-2, evidence has been reported in several previous outbreaks of coronaviruses, directly supporting this hypothesis (He and Han 2021). There exists a strong link between anthropogenic activities and spillover of zoonotic pathogens (He et al. 2021). One evident example is that, in Peru, the destruction of thousands of acres of rainforests created new mosquito habitat and caused more than 120,000 cases of malaria in the late 1990s, compared with only less than 150 cases per year in the early 1990s [START_REF] Shah | The spread of new diseases and the climate connection[END_REF][START_REF] Vittor | The effect of deforestation on the human-biting rate of Anopheles darlingi, the primary vector of Falciparum malaria in the Peruvian Amazon[END_REF]. Climate change may also accelerate the transmission of infectious agents to and between humans [START_REF] Beyer | Shifts in global bat diversity suggest a possible role of climate change in the emergence of SARS-CoV-1 and SARS-CoV-2[END_REF]). Overall, it seems inevitable that humans will encounter other infectious disease outbreaks in a not distant future.

With two-thirds of the global population expected to live in cities by 2050, more people will be exposed to infectious agents in public spaces and through social interactions [START_REF] Fulmer | What in the world is infrastructure? PEI infrastructure investor[END_REF]. The global outbreak of COVID-19-along with the recent re-emergent outbreaks caused by its variants-has taught humans a painful lesson. The organization of modern human society must acknowledge these looming threats and adapt to the challenges and future uncertainties, through forward-thinking urban planning and by building more resilient cities, communities, and infrastructure. Since population density and movements are major controlling factors [START_REF] Tatem | Mapping population and pathogen movements[END_REF], de-centralizing the human society, therefore, will be a fundamental factor to attain a more robust societal structure when facing the threat of a fast-spreading pandemic.

There are realistic actions we can take in the near future to move toward this goal. Most human migrations are driven by resources. Young talents and labors seek work and opportunities, which are often more abundant in cities and industrialized regions. Spreading resources such as employments, education, and healthcare geographically, therefore, will fundamentally retard further migration of the population into large cities. Beijing, for instance, has been doing this for several years, and the population inflow has indeed slowed [START_REF] Du | Population of Beijing growing, aging[END_REF]). In addition, governments and municipalities may offer incentives, e.g., tax breaks and rent subsidies for enterprises and employers to relocate to smaller cities, towns, and countryside while at the same time, invest in building more infrastructure, schools, and hospitals to accommodate the population inflow. The less expensive and more natural living environments are major factors attracting people to move to and settle in smaller cities and towns. While the COVID-19 is not eradicated, people can make use of technologies to maintain productivity and avoid the risks of virus infection in public transport or work environments by working remotely from home. In the long term, we should build cities with more open space, spread-out infrastructure, and more even population densities in cities, suburbs, and surrounding towns. It is also helpful to provide more areas in buildings and public transport facilities to allow more social distances with means of using natural ventilation in air-conditioned enclosed spaces to lower the risk of transmission of human pathogens [START_REF] Wang | Unprotected mothers and infants breastfeeding in public amenities during the COVID-19 pandemic[END_REF]. Lastly, it is imperative to build more resilient waste management infrastructures to mitigate the possible transmission of infectious diseases in extreme weathers and natural disasters (Han and He 2021), which have become increasingly common due to anthropogenic activities and climate change.

Fig. 2

 2 Fig. 2 A metro station as in Paris, France, during the COVID-19 pandemic in February 2021, showing a large crowd of people waiting for the subway train. Although most people were wearing face masks, the efficacy of protection is questionable in the very crowded space with

Acknowledgements This work was funded by the Young Talent Support Plan of Xi'an Jiaotong University. The authors thank Ms. Xueting Yang for her assistance in collecting data and evidence in the early stage of preparing this manuscript.

Declarations

Conflict of interest The authors declare that they have no conflict of interest in this work.


Dominique Planson
email: dominique.planson@insa-lyon.fr

C Sonneville
email: bcamille.sonneville@insa-lyon.fr

Pascal Bevilacqua
email: cpascal.bevilacqua@insa-lyon.fr

L V Phung

B Asllani
email: ebesar.asllani@insa-lyon.fr

D Tournier

P Brosselard

V Phung

4H-SiC PiN diode protected by narrow field rings investigated by the micro-OBIC method

Keywords: 4H-SiC bipolar diode, high voltage device, electro-physical characterisation, periphery protection, OBIC

come

Introduction

Silicon Carbide (SiC) is an attractive semiconductor material for high power and high temperature applications. To reach the expected breakdown voltage, an effective junction termination is required to spread the electrical field profile towards edges mitigating the crowding effect. Among different techniques used for periphery protection [START_REF] Godignon | this conference[END_REF], Junction Termination Extension (JTE) is a well-known technique often chosen due to their relative design simplicity. However, the structure of the present paper rely exclusively on guard rings, a more challenging approach where small technological deviations on the rings position can completely negate the benefits of the periphery protection. To get the most from those rings, key parameters such as the distance between each ring and the number of rings must be carefully chosen. The efficiency of this periphery protection is optimal when the electric field at the edge of each ring is nearly equal for all rings. Optical Beam Induced Current (OBIC) technique is a non-destructive powerful technique to analyze the electric field distribution at the edge of the device periphery. This technique has been used to investigate silicon diodes protected by field rings [START_REF] Stengl | High-voltage planar junctions investigated by the OBIC method[END_REF], with spacing greater than 10 µm, as well as 4H-SiC PiN diodes protected by JTE [3], [4]. In this paper, micro-OBIC technique is applied to SiC high voltage bipolar diodes where the design of the guard rings is more tightly constrained to analyze the experimental behavior of the periphery protection with respect to the applied reverse voltage.

Experimental Setup

In the micro-OBIC experimental setup, a 349 nm UV pulsed laser is used to generate electron-hole pairs (EHPs). A set of properly adjusted optics (mirrors, microscope objective …) is inserted on the beam path to finally get a focused laser spot with a diameter of about 1-4 µm as shown on Fig. 1. The sample is placed on a motorized stage and the position is controlled with LabView (i.e. on X, Y, Z automatically and teta manually) of the focal point. The energy at the top surface of the device is high enough to generate electron-hole pairs (EHPs), so that an OBIC current could be measured. The micro-OBIC principle has already been described in [4].

The testbench allows us to realize a spatial mapping of OBIC signal, either on a (X, Y) plane or only on X (or Y) lines. For X or Y lines, the step between two displacements of the laser beam is 200 nm, Anode-Cathode biased. The DUT can be reverse biased by applying voltage up to 500V thanks to the SMU (Source Measure Unit) Keithley 237. The cross section of the DUT and its periphery is sketched on the Figure 2. The epilayer is 14.8 µm thick, n-type doped (6´10 15 cm -3 ). The main junction is realized by ion implantation, at the same time as the field rings, which were carefully designed by Caly Technologies. The diode is then protected by a passivation layer on the front side.

The measurements were realized on two different square-shaped bipolar diodes protected by field rings, namely device 1 and device 2, with and without passivation layer respectively. The metallized anode contact is square-shaped, with a side length of 200 µm x 300 µm, and a corner curvature radius of 200 µm as shown in Fig. 2-b.

This non-destructive technique allows to measure OBIC currents by scanning the outer edge of the biased diode, in order to assess the periphery protection efficiency. OBIC current increases as the electric field intensity increases in the structure. There is no OBIC current on the metallization, neither far from the periphery of the diode. These OBIC current profiles would come in handy to the designer to optimize the breakdown voltage. The diode is protected by 22 field rings as shown in Fig. 2, with increasing spacing between them (0.07 µm from one ring to the next). The distance between the main junction and the first ring is 1.2 µm. All the rings have the same width (1.5 µm).

Inside a vacuum probing chamber, on-wafer measurements of the device demonstrated a stable blocking voltage up to 2 kV as shown in Fig. 3.

Results and Analysis

By applying different reverse voltages (up to 500V) to the same diode and re-scanning the very same location, several OBIC currents can be compared to each other, as shown in Fig. 4. As the reverse voltage increases, the OBIC current increases for the same location, due to the increase of the electric field. One can also observe the extension of the space charge region as the reverse voltage increases and the role of outer field rings, witness of the peripheral protection good design. As the reverse voltage increase, the number of OBIC peaks increase too. From 20 µm up to 97 µm, each maximum is related to the presence of a ring as described in [START_REF] Godignon | this conference[END_REF] and, more precisely, is located on the outer edge of each ring junction. As the reverse voltage increases, the presence of the rings is more clearly revealed by µ-OBIC. Indeed, the number visible rings by OBIC increase with the reverse voltage. These measurements up to 500V, show the presence of very close rings (between 1.2 µm for the first ring and 2.67 µm for the last ring). One can count all the 22 field rings for a reverse voltage of 500V plus one more peak for the border of the main junction. That means that at 500V all the rings are revealed by the OBIC method. The amplitude of the maximum doesn't increase much from one voltage to another. The distance between the peaks is in total accordance with the layout of the rings as shown in the figure and in Table 1. Table 1 gives the position of each ring based on its left lateral edge that is oriented towards the main junction when viewed from the cross-section depicted in fig. 2a. The OBIC signal amplitude is given for each left lateral edge. Comparing the data from Table 1 and the OBIC profile shown in Fig. 4, each dip in the profile perfectly matches with the left lateral edge of each ring. All the data in the present table were extracted from diode #2 biased under 500 V. After measuring diode #1, some scratches appear at the surface of the device under test, due to the sweeping of the laser. In fact, this device was passivated with a polyimide layer, allowing a uniform spreading of the electric field as shown in Fig. 4. With respect to the beam penetration, this additional layer acts as a filter, and the OBIC signal is reduced by a factor of 100 approximately. Nonetheless, Figure 4 shows the benefits of polyimide layer at the top in terms of electric field distribution: at 500V the OBIC current is almost constant whereas Figure 5 shows a decrease in OBIC current value. With a voltage of 300V, one can observe all the rings. Fig. 4. OBIC currents at the periphery of the diode device #1 (with polyimide) for different reverse voltages (ranging from 0V up to 500V). The mask layout with the spacing between the rings is in grey color.

In order to study the periphery protection without any passivation layer, the polyimide was removed with a very high-density plasma. Then, the same experiment has been performed on diode without polyimide layer at the top. The applied reverse voltage was also limited to 500V. The results are shown in Fig. 5.

Compared to the previous measurements on both devices, performing measurements without the polyimide layer give OBIC profile that reflects the peripheral protection more regular. In both cases, the peaks are related to the distance between the rings, and that is true even for very small distance as low as 1.2 µm.

Fig. 5. OBIC currents at the periphery of the diode device #2 (without polyimide) for different reverse voltages (ranging from 0V up to 500V by 50V step). The mask layout with the spacing between the rings is in grey color.

Fig. 1 .

 1 Fig. 1. Schematic representation of the micro-OBIC testbench.

Fig. 2 .

 2 Fig. 2. a) Schematic cross section view of half vertical PiN diodes protected by 22 field rings, b) photograph of the diode with the 50 µm (diameter) gold ball bonding

Fig. 3

 3 Fig. 3 Reverse characteristics of both devices (diode #1 with polyimide passivation, diode #2 without any passivation)

Table 1 :

 1 Distance from the main junction for ring mask layout and respective OBIC signal values.

	Ring	#1	#2	#3	#4	#5	#6	#7	#8	#9	#10 #11
	number
	Distance	1.2	3.97	6.81	9.72	12.7	15.75 18.87	22.06	25.32 28.65 32.05
	from main
	junction
	(µm)
	OBIC signal 1.19 3.99 6.99 9.79 12.79 16	18.8	22.21	25.38 28.79 32.2

Summary

This paper presents OBIC measurements performed on field rings protected bipolar diodes at different voltages. This test-bench has a well improved spatial resolution compared to the previous OBIC test-bench. We estimated theoretically and experimentally the spatial resolution of the micro-OBIC to be around 1-4 μm. A better resolution is interesting to characterize peripheral protections, such as JTE rings, whose lateral extension is of the order of magnitude of the micron. JTE rings are directly observable in these experiments, even with narrow distances between each 22 rings. It would be then interesting to complete this study by measurements at higher reverse bias under vacuum and to be closed to the breakdown voltage of the device (1 700 V). It would be also interesting to study these same diodes by other complementary characterization methods such as micro-Raman spectroscopy, SEM/EDX and TOF-SIMS.