Nalini Anantharaman

SOME RELATIONS BETWEEN THE SPECTRA OF SIMPLE AND NON-BACKTRACKING RANDOM WALKS

Keywords:

We establish some relations between the spectra of simple and non-backtracking random walks on non-regular graphs, generalizing some well-known facts for regular graphs. Our two main results are 1) a quantitative relation between the mixing rates of the simple random walk and of the non-backtracking random walk 2) a variant of the "Ihara determinant formula" which expresses the characteristic polynomial of the adjacency matrix, or of the laplacian, as the determinant of a certain non-backtracking random walk with holomorphic weights.

Results

It has been noted by many authors that non-backtracking random walks on graphs looking locally like trees are simpler, from a combinatorial point of view, than usual random walks. This is for instance an ingredient of the proof of Alon's conjecture on random regular graphs by Friedman [START_REF] Friedman | A proof of Alon's second eigenvalue conjecture and related problems[END_REF] or, more recently, by Bordenave [START_REF] Bordenave | A new proof of friedman's second eigenvalue theorem and its extension to random lifts[END_REF]. The study of the spectrum of non-backtracking random walks is also at the heart of the "community detection" problem and of the solution to the "spectral redemption conjecture" for various models of random graphs [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular Ramanujan graphs[END_REF]. Non-backtracking random walks are known to mix faster than the usual ones [START_REF] Noga Alon | Non-backtracking random walks mix faster[END_REF]. In [START_REF] Lubetzky | Cutoff on all Ramanujan graphs[END_REF], the non-backtracking random walk is used to prove a cut-off phenomenon for the usual random walk on regular graphs. In geometric group theory, the "cogrowth" is directly related to the leading eigenvalue of non-backtracking random walks : this is discussed in [START_REF] Ortner | Non-backtracking random walks and cogrowth of graphs[END_REF], where a possibility to use this to extend the notion of co-growth to non-regular graphs is suggested. The non-backtracking random walk may also be used as an analog of "classical dynamics" in the field of quantum chaos on discrete graphs : see the papers by Smilansky [START_REF] Smilansky | Discrete graphs -a paradigm model for quantum chaos[END_REF][START_REF] Smilansky | Quantum chaos on discrete graphs[END_REF], which were a source of inspiration for this work. This note is a contribution to the study of the relation between the spectra of simple and non-backtracking random walks, for non-regular graphs. It originates in the work [START_REF] Anantharaman | Poisson kernel expansions for schrödinger operators on trees[END_REF], where we used non-backtracking random walks to study the quantum ergodicity problem for eigenfunctions of Schrödinger operators on non-regular expander graphs.

Let G = (V, E) be a graph without multiple edges and self-loops. We assume that the degree D(x) of a vertex x is bounded above and below : 2 ≤ D(x) ≤ D. In the results about spectral gaps, we actually have to assume that D(x) ≥ 3. We are interested in relating the spectrum of various operators, such as the laplacian, the adjacency matrix, and certain weighted non-backtracking random walks. Such relations are well-known for regular graphs (i.e. those for which D(x) is constant), and the goal of this note is to partially extend what is known to non-regular graphs. Our two main results are 1) a relation between the mixing rates of the simple random walk and of the non-backtracking random walk 2) a variant of the "Ihara determinant formula" [START_REF] Ihara | On discrete subgroups of the two by two projective linear group over p-adic fields[END_REF][START_REF] Ihara | Discrete subgroups of PL(2, k℘)[END_REF] which expresses the characteristic polynomial of the adjacency matrix, or of the laplacian, as the determinant of a certain non-backtracking transfer matrix with holomorphic coefficients.

We will write y ∼ x to mean that y is a neighbour of x in G.

The first operator we are interested in is the adjacency matrix A. It acts on C V by the formula Af (x) = y∼x f (y), and is self-adjoint on ℓ 2 (V, u) if V is endowed with the uniform measure u.

We are also interested in the spectrum of the laplacian, defined by :

P : C V -→ C V P f (x) = 1 D(x) y∼x f (y). (1)
If we endow the set of vertices V with the measure π(x) = D(x), then P is self-adjoint on ℓ 2 (V, π). Note that this implies x P f (x)π(x) = x f (x)π(x). The space ℓ 2 o (V, π) of functions orthogonal to constants in ℓ 2 (V, π) is preserved by P .

Let B be the set of oriented edges, endowed with the uniform measure U (each edge has weight 1). Denoting Q(x) = D(x) -1, the "transfer operator" is defined by

S : ℓ 2 (B, U ) -→ ℓ 2 (B, U ) Sf (e) = 1 Q(o(e)) e ′ ;e f (e ′ ) (2)
where e ′ ; e means that t(e ′ ) = o(e) and e is not the reverse of e ′ . The operator S is stochastic, it is the generator of the non-backtracking random walk. It is not self-adjoint, but we have e Sf (e) = e f (e). This is equivalent to saying that S * : ℓ 2 (B, U ) -→ ℓ 2 (B, U ) is also stochastic. When G is finite, this implies that S preserves the space ℓ 2 o (B, U ) of functions orthogonal to constants in ℓ 2 (B, U ). We will finally use the non-stochastic operator

B : C B -→ C B Bf (e) = e ′ ;e f (e ′ ) (3)
If the graph G is finite and (q + 1)-regular, then A and P (resp. B and S) are the same operator up to a homothety, and there is an explicit relation between the characteristic polynomials of A and B :

(4) det(I |B| -uB) = (1 -u 2 ) r-1 det((1 + u 2 q)I |V | -uA)
where r = |E| -|V | + 1 is the rank of the fundamental group. This is the contents of the Ihara determinant formula [START_REF] Ihara | On discrete subgroups of the two by two projective linear group over p-adic fields[END_REF][START_REF] Ihara | Discrete subgroups of PL(2, k℘)[END_REF], generalised in stages by Hashimoto, Bass and Kotani-Sunada [START_REF] Hashimoto | Zeta functions of finite graphs and representations of p-adic groups[END_REF][START_REF] Hashimoto | Selberg-Ihara's zeta function for p-adic discrete groups[END_REF][START_REF] Hashimoto | On zeta and L-functions of finite graphs[END_REF][START_REF] Hashimoto | Artin L-functions of finite graphs and their applications[END_REF][START_REF] Hashimoto | Artin type L-functions and the density theorem for prime cycles on finite graphs[END_REF][START_REF] Bass | The Ihara-Selberg zeta function of a tree lattice[END_REF][START_REF] Kotani | Zeta functions of finite graphs[END_REF]. For finite non-regular graphs the relation reads

(5) det(I |B| -uB) = (1 -u 2 ) r-1 det(I |V | -uA + u 2 Q)
where Q is the diagonal matrix with components Q(x) = D(x) -1. Note that the righthand side in ( 5) is not directly related to the characteristic polynomial of A. The identity (5) relates eigenvalues (and eigenvectors) of S to solutions of (

I |V | -uA + u 2 Q)v = 0,
which are not eigenvectors of A. Our goal is twofold :

• for finite graphs, compare the mixing rates of S and P . For regular graphs, there is an exact relation between eigenvalues of S and P , which implies that the spectral gap of P on ℓ 2 o (V, π) is explicitly related to the spectral gap of S on ℓ 2 o (B, U ). However since S is not self-adjoint, the knowledge of its spectral gap is not sufficient to determine its mixing rate, one needs to control the angles beween eigenvectors. This analysis, done in [START_REF] Anantharaman | Quantum ergodicity on large regular graphs[END_REF] and with more details in [START_REF] Lubetzky | Cutoff on all Ramanujan graphs[END_REF], uses the fact that the eigenvectors of S are explicitly related to those of P . Such explicit relations are not available for non-regular graphs and we have to find a more general method; • extend formula (4) to non-regular graphs in a way different of ( 5), by finding an identity involving the characteristic polynomial of A on the right-hand side.

The result about spectral gaps also holds for infinite graphs. We recover, in a less geometric but more quantitative way, the result of Ortner and Woess [START_REF] Ortner | Non-backtracking random walks and cogrowth of graphs[END_REF] saying that P has a spectral gap on ℓ 2 (V, π) if and only if the spectral radius of S on ℓ 2 (B, U ) is strictly less than 1.

We do not discuss the "spectral gap" of A as this is not as properly defined as for P (the top eigenvalue and eigenvector of A are not explicit in general).

1.1. Spectral gap and mixing rate for P and S. Let us first assume that G is finite, connected and non-bipartite. This is equivalent to assuming that 1 is a simple eigenvalue of P 2 . In other words, the spectrum of

P 2 in ℓ 2 o (V, π) is contained in [0, 1 -β],
for some β > 0 which measures the mixing rate of the simple RW on G.

Our first result is the following :

Theorem 1.1.
Assume that G is finite and that D(x) ≥ 3 for all x. Assume that the spectrum of

P 2 on ℓ 2 o (V, π) is contained in [0, 1 -β]. Then the spectrum of S * 2 S 2 on ℓ 2 o (B, U ) is contained in [0, 1 -c(D, β)],
where c(D, β) depends only on D and β, and is positive if β is so.

Note that S ℓ 2 o (B,U )-→ℓ 2 o (B,U ) = 1, as Sf = f as soon as f is a function on B
that is constant on edges having the same terminus. However, our theorem says that S 2

ℓ 2 o (B,U )-→ℓ 2 o (B,U ) ≤ (1 -c(D, β)) 1/2 . The value of c(D, β) is given in (16). Corollary 1.2. For all n ≥ 1, S n ℓ 2 o (B,U )-→ℓ 2 o (B,U ) ≤ (1 -c(D, β)) ⌊n/4⌋ .
This gives the rate of mixing of the non-backtracking RW. The converse is easier : in the course of the proof, we will also see that if the spectrum of

S * 2 S 2 on ℓ 2 o (B, U ) is contained in [0, 1 -c], then the spectrum of P 2 on ℓ 2 o (V, π) is contained in [0, 1 -D -2 c] (Remark 2.2).
We were primarily interested in finite graphs in view of the application to quantum ergodicity [START_REF] Anantharaman | Poisson kernel expansions for schrödinger operators on trees[END_REF], but the result also holds for infinite graphs : Theorem 1.3. Assume that G is infinite and that D(x) ≥ 3 for all x.

(i) If the spectrum of S * 2 S 2 on ℓ 2 (B, U ) is contained in [0, 1c], then the spectrum of

P 2 on ℓ 2 (V, π) is contained in [0, 1 -D -2 c]. (ii) If the spectrum of P 2 on ℓ 2 (V, π) is contained in [0, 1 -β], then the spectrum of S * 2 S 2 on ℓ 2 (B, U ) is contained in [0, 1 -c(D, β)],
where c(D, β) is given by (16).

Remark 1.4. It is well-known that G is amenable iff the spectral radius of P is 1 (see [START_REF] Kesten | Full Banach mean values on countable groups[END_REF][START_REF] Dodziuk | Combinatorial Laplacians and isoperimetric inequality[END_REF][START_REF] Dodziuk | Difference equations, isoperimetric inequality and transience of certain random walks[END_REF]). Thus Theorem 1.3 says that G is amenable iff S 2 ℓ 2 (B,U )-→ℓ 2 (B,U ) = 1. It was proven before by Ortner and Woess that G is amenable iff the spectral radius of S is 1 [START_REF] Ortner | Non-backtracking random walks and cogrowth of graphs[END_REF]. This can be recovered by our methods : indeed, one direction results from our Theorem 1.3; in the other direction, one can follow the same lines as in our Remark 2.2 to show that if

S * n S n ℓ 2 (B,U )-→ℓ 2 (B,U ) < 1 then P 2(n-1) ℓ 2 (V,π)-→ℓ 2 (V,π) < 1 (with an explicit bound).
Our method is very down-to-earth and gives, by basic manipulations, a quantitative relation between the spectral gap of P and S 2 ℓ 2 (B,U )-→ℓ 2 (B,U ) . The method in [START_REF] Ortner | Non-backtracking random walks and cogrowth of graphs[END_REF] is less direct and more geometric : it starts from the general fact that G is amenable iff SOLG (the symmetrized oriented line graph) is amenable. And then it is shown that SOLG is amenable iff the spectral radius of S is 1.

1.2. Determinant relation. We now assume that G is finite.

Let T = (V (T ), E(T )) be the universal cover of G : T is a tree, and there exists a subgroup Γ of of automorphism group of T , acting without fixed points on V (T ), such that G = Γ\T . Let à be the adjacency matrix of T . The Green function on T will be denoted by G(x, y; z) = δ x , ( Ãz) -1 δ y ℓ 2 (V (T ))

for z ∈ C \ R. Given v, w ∈ T with v ∼ w, we denote by T (v|w) the tree obtained by removing from T the branch emanating from v that passes through w. We define the restriction H (v|w) (x, y) = H(x, y) if v, w ∈ T (v|w) and zero otherwise. We then denote G (v|w) (•, •; z) the corresponding Green function.

Given z ∈ C \ R, v ∈ V , w a neighbour of v we denote G z (v) = G(ṽ, ṽ; z) and ζ z (w, v) = -G (ṽ| w) (ṽ, ṽ; z) .

where (ṽ, w) is a lift of the edge (v, w) in T . This definition does not depend on the choice of the lifts. If e = (w, v) ∈ B, we also use the notation

G z (e) = G( w, ṽ; z), ζ z (e) = ζ z (w, v).
Note that G z (e) is invariant under edge-reversal, whereas ζ z (e) is not. In the formula below, the function ζ z on B acts on C B as a multiplication operator.

Theorem 1.5. For all z ∈ C \ R, (6) e∈E
(-G z (e)) • det Ä (ζ z ) -1 I |B| -B ä = det Ä zI |V | -A ä • x∈V (-G z (x))
Remark 1.6. In the case of a (q +1)-regular graph, ζ z is a constant function, which solves the quadratic equation ( 7)

z = qζ z + 1 ζ z See Lemma 3.1 below. We also have G z (x) = ζ z (ζ z ) 2 -1 and G z (e) = (ζ z ) 2 (ζ z ) 2 -1
for all x and all e. It can be checked that Theorem 1.5 reduces to (4) by setting u = ζ z . It is, however, different from (5) for non-regular graphs. Although extensions of the Ihara formula (4) have been studied by many authors, the variant (6) seems to be new.

Note that (6) holds for any functions ζ z that are solutions of the system of algebraic equations appearing in Lemma 3.1. These are by no means unique : for instance, in the regular case, there are 2 solutions to equation [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular Ramanujan graphs[END_REF]. It is nice, however, to know an explicit solution of this system that can be expressed in terms of Green functions.

Remark 1.7. The theorem generalizes to the case where A is replaced by a discrete "Schrödinger operator" of the form A p +W : C V -→ C V , (A p +W )f (x) = y∼x p(x, y)f (y)+ W (x)f (x) where W is a real-valued function on V , and p is such that p(x, y) = p(y, x) ∈ R and p(x, y) = 0 iff x ∼ y. The definitions of the Green functions G z and ζ z should be modified (in the obvious manner) to incorporate the weights p and the potential W . The definition of B should be modified to p(e)

B p f (e) =
• det

Ä (ζ z ) -1 I |B| -B p ä = det Ä zI |V | -A p -W ä • x∈V (-G z (x)).
This remark, in particular, allows to cover the case of det(zI |V | -P ), noting that P is conjugate to A p with p(x, y) = (D(x)D(y)) -1/2 . upon work supported by the Agence Nationale de la Recherche under grant No.ANR-13-BS01-0007-01, . I am very thankful to Mostafa Sabri for his careful reading and numerous useful comments on the manuscript.

Proof of Theorem 1.1

We start by noting that if

f ∈ ℓ 2 o (V, π), then (8) 1 2 x∈V 1 D(x) y,y ′ ∼x |f (y) -f (y ′ )| 2 ≥ β f 2 ℓ 2 (V,π) .
This just comes from the identity

1 2 x∈V 1 D(x) y,y ′ ∼x |f (y) -f (y ′ )| 2 = x D(x)|f (x)| 2 - x D(x)|P f (x)| 2 = f, (I -P 2 )f ℓ 2 (V,π)
To prove Proposition 1.1, we use the following decomposition of the space of functions on B :

(9) ℓ 2 o (B, U ) = O(ℓ 2 o (V, π)) ⊕ T (ℓ 2 o (V, π)) ⊕ (O(ℓ 2 o (V, π)) ⊥ ∩ T (ℓ 2 o (V, π)) ⊥ ). The space O(ℓ 2 o (V, π)) is the image of ℓ 2 o (V, π) under the map Of (e) = f (o(e)). Thus O(ℓ 2 o (V, π))
is the space of functions (orthogonal to constants) such that f (e) depends only on the origin of e. Similarly, T (ℓ 2 o (V, π)) is the space of functions such that f (e) depends only on the terminus of e. It is the image of

ℓ 2 o (V, π) under the map T f (e) = f (t(e)). If G is non-bipartite, we have O(ℓ 2 o (V, π)) ∩ T (ℓ 2 o (V, π)) = {0}. Note that each space O(ℓ 2 o (V, π)), T (ℓ 2 o (V, π)) is orthogonal to 1l in ℓ 2 (B, U
), but that the two spaces are NOT orthogonal to each other. The space Remark 2.1. In the infinite case, the proof of Theorem 1.3 will be similar, using the decomposition

O(ℓ 2 o (V, π)) ⊥ ∩ T (ℓ 2 o (V, π)) ⊥ is of dimension 2|E| -2|V | + 1 = r -1,
ℓ 2 (B, U ) = O(ℓ 2 (V, π)) ⊕ T (ℓ 2 (V, π)) ⊕ (O(ℓ 2 (V, π)) ⊥ ∩ T (ℓ 2 (V, π)) ⊥ ).
To start the proof we use the Dirichlet identity for f ∈ ℓ 2 (B, U ):

f, (I -S * 2 S 2 )f ℓ 2 (B,U ) = 1 2 e,e ′ |f (e) -f (e ′ )| 2 S * 2 S 2 (e, e ′ ).
Let us decompose f according to [START_REF] Dodziuk | Difference equations, isoperimetric inequality and transience of certain random walks[END_REF] :

f = F + G + H where F ∈ O(ℓ 2 o (V, π)), G ∈ T (ℓ 2 o (V, π)), H ∈ (O(ℓ 2 o (V, π)) ⊥ ∩ T (ℓ 2 o (V, π)) ⊥
). We are first going to prove that [START_REF] Hashimoto | Zeta functions of finite graphs and representations of p-adic groups[END_REF] f,

(I -S * 2 S 2 )f ℓ 2 (B,U ) ≥ Q -4 β F + H 2 ℓ 2 (B,U )
where Q = D -1 (and, recall, D is an upper bound on the degree).

In order to have S * 2 S 2 (e, e ′ ) > 0, there must exist e 1 , e ′ 1 , e 2 ∈ B such that e ; e 1 ; e 2 and e ′ ; e ′ 1 ; e 2 . Counting the number of possibilities, we see that S * 2 S 2 (e, e ′ ) ≥

(D(t(e)) -2)Q -4 ≥ Q -4 if t(e) = t(e ′ ).
Here we use the assumption that D(t(e)) ≥ 3. Thus,

f, (I -S * 2 S 2 )f ℓ 2 (B,U ) ≥ Q -4 2 e,e ′ :t(e)=t(e ′ ) |f (e) -f (e ′ )| 2 ≥ Q -4 2 e,e ′ :t(e)=t(e ′ ) 1 D(t(e)) |f (e) -f (e ′ )| 2 = Q -4 2 e,e ′ :t(e)=t(e ′ ) 1 D(t(e)) |(F + H)(e) -(F + H)(e ′ )| 2 .
Let us fix a vertex x ∈ V . Using the fact that F depends only on the origin, and that H satisfies ( 10),

e,e ′ ,t(e)=t(e ′ )=x

|(F + H)(e) -(F + H)(e ′ )| 2 = y,y ′ ∼x |F (y) -F (y ′ )| 2 + e,

F (e)H(e)

Summing now over x, and using the fact that F and H are orthogonal,

e,e ′ ,t(e)=t(e ′ )

1 D(t(e)) |(F + H)(e) -(F + H)(e ′ )| 2 = x D(x) -1 y,y ′ ∼x |F (y) -F (y ′ )| 2
|F (y) -F (y ′ )| 2 + 2 e |H(e)| 2 ≥ x D(x) -1 y,y ′ ∼x |F (y) -F (y ′ )| 2 + 2 H 2 ℓ 2 (B,U ) ≥ 2β F 2 + 2 H 2 ≥ 2β F + H 2 .
On the last line, we have used [START_REF] Dodziuk | Combinatorial Laplacians and isoperimetric inequality[END_REF]. This concludes the proof of [START_REF] Hashimoto | Zeta functions of finite graphs and representations of p-adic groups[END_REF]. Now, let G ∈ T (ℓ 2 o (V, π). Then again, by looking at what it means to have S * 2 S 2 (e, e ′ ) > 0, we see that if y, y ′ are two vertices such that dist(y, y ′ ) = 2 (in other words, y and y ′ have a common neighbour x), then we can find edges e, e ′ such that t(e) = y, t(e ′ ) = y ′ and S * 2 S 2 (e, e ′ ) ≥ Q -4 . Indeed, we may choose e, e ′ such that t(e) = y and o(e) = x, t(e ′ ) = y ′ and o(e ′ ) = x, and

S * 2 S 2 (e, e ′ ) ≥ (D(x) -2)Q -4 ≥ Q -4 . Thus G, (I -S * 2 S 2 )G ℓ 2 (B,U ) = 1 2 e,e ′ |G(e) -G(e ′ )| 2 S * 2 S 2 (e, e ′ ) ≥ Q -4 2 x y,y ′ ∼x |G(y) -G(y ′ )| 2 ≥ Q -4 β G 2 .
Remark 2.2. We can also write (using the fact that S * 2 S 2 is stochastic)

G, (I -S * 2 S 2 )G ℓ 2 (B,U ) = 1 2 y,y ′ :d(y,y ′ )=2 |G(y) -G(y ′ )| 2
e,e ′ :t(e)=y,t(e ′ )=y ′ S * 2 S 2 (e, e ′ )

≤ D 2 1 2 x∈V 1 D(x) y,y ′ :y∼x∼y ′ |G(y) -G(y ′ )| 2 = D 2 G, (I -P 2 )G ℓ 2 (V,π)
and this proves part (i) of Theorem 1.3.

Let A > 1 (to be chosen later, depending on β and D). Let f = F + G + H as before. Assume first that A F + H ≥ G . Then by the triangular inequality f ≤ (1 + A) F + H . In addition, as we have seen, [START_REF] Hashimoto | Artin L-functions of finite graphs and their applications[END_REF] f, (

I -S * 2 S 2 )f ℓ 2 (B,U ) ≥ Q -4 β F + H 2 ℓ 2 (B,U ) ≥ Q -4 β(1 + A) -2 f 2 .
Otherwise, A F + H ≤ G , and f ≤ (1 + A -1 ) G . Noting that the operator norm of I -S * 2 S 2 is less than 1, we write for all f = F + G + H,

(15) f, (I -S * 2 S 2 )f ℓ 2 (B,U ) ≥ G, (I -S * 2 S 2 )G ℓ 2 (B,U ) -2A -1 G 2 -A -2 G 2 ≥ (Q -4 β -3A -1 ) G 2 ≥ (Q -4 β -3A -1 ) (1 + A -1 ) 2 f 2 .
Choosing A such that A -1 = Q -4 β/6, and gathering ( 15) and ( 14) we get the result with ( 16)

c(D, β) = min Ç Q -4 β 2(1 + Q -4 β/6) 2 , Q -4 β (1 + 6Q 4 /β) 2 å .

Proof of the determinant relation

The relations in the next lemma follow from the resolvent identity, and are proven (for instance) in [START_REF] Anantharaman | Poisson kernel expansions for schrödinger operators on trees[END_REF]. For a vertex v of T , N v stands for the set of neighbouring vertices.

Lemma 3.1. For any v ∈ V (T ), z = E + iη ∈ C + , if we let 2m z (v) = -1 G(v,v;z) , we have z = u∼v ζ z (v, u) + 2m z (v) and z = u∈Nv\{w} ζ z (v, u) + 1 ζ z (w, v) .
For any non-backtracking path (v 0 , . . . , v k ) in T ,

G(v 0 , v k ; z) = -k-1 j=0 ζ z (v j+1 , v j ) 2m z (v k ) = -k-1 j=0 ζ z (v j , v j+1 ) 2m z (v 0 ) . (17)
Also, for any w ∼ v, we have

(18) ζ z (w, v) = m(w) z m(v) z ζ z (v, w) , 1 ζ z (w, v) -ζ z (v, w) = 2m z (v) , Remark 3.2.
We can note that (17) may be written as

((ζ z ) -1 I |B| -B) -1 (e, e ′ ) = δ x=y + +∞ k=0 ζ z (e ′ )(ζ z B) k (e, e ′ ) = -2m z (x)(A -z) -1 (x, y)
for all e, e ′ ∈ B and x = o(e ′ ), y = t(e).

In the case of regular graphs, this is formula (2.4) in [START_REF] Ortner | Non-backtracking random walks and cogrowth of graphs[END_REF], where it is attributed to Grigorchuk, with various proofs published by Woess, Szwarc [START_REF] Woess | Cogrowth of groups and simple random walks[END_REF][START_REF] Szwarc | A short proof of the Grigorchuk-Cohen cogrowth theorem[END_REF], Northshield, Bartholdi [START_REF] Northshield | Cogrowth of regular graphs[END_REF][START_REF] Bartholdi | Counting paths in graphs[END_REF].

3.1. Operator relations. In this section z ∈ C + is fixed, so we write ζ(x, y) instead of ζ z (x, y), m(x) instead of m z (x). If e = (x, y) ∈ B, we write m 1 (e) = m(x) and m 2 (e) = m(y). A function on B defines a multiplication operator on C B (i.e. an operator which is diagonal in the canonical basis). We use the same notation for a function and the associated operator.

Let us introduce the notation

Pf (x) = 1 D(x) y∼x f (x, y).
This is a projector on the space of functions depending only on the origin, which may be identified with ℓ 2 (V, π), isometrically embedded into ℓ 2 (B, U ) by the map ψ → O(ψ) defined in the previous section. Let L = D(2m 1 ) -1 P. Let

Hg(x) = y,y∼x 1 2m(y) (ζ(y, x)g(y, x) -g(x, y))
Theorem 1.5 is based on the following exact relation :

Proposition 1. H • (ζ -1 I -B) = (A -zI) • L Proof. Let φ = Lf and g = -(ζ -1 I -B)f . The latter relation implies that for any y ∼ x, φ(x) = (2m(x)) -1 Ç f (x, y) + g(y, x) + f (y, x) ζ(y, x) å .
We then calculate which is the desired relation.

Aφ(x) = y,y∼x φ(y) = y,y∼x 1 2m(y)
We note that Hg = -D P((2m 2 ) -1 g) + P((2m 2 ) -1 ι(ζg)) where ιf (x, y) = f (y, x) is the edge reversal involution. So H itself is of the form H = DP • K where K = (2m 2 ) -1 (ιζ -I). We have proven that This yields the announced relation.

e

  ′ ;e p(e ′ )f (e ′ ) and (6) becomes e∈E (-G z (e))

  where r is the rank of the fundamental group of G. By definition, it is the space of functions f : B -→ C such that, for all x ∈ V ,[START_REF] Friedman | A proof of Alon's second eigenvalue conjecture and related problems[END_REF] e,o(e)=x f (e) = 0 and e,t(e)=x f (e) = 0.

Çff

  (y, x) + g(x, y) + f (x, y) ζ(x, y) y, x)2m(x)φ(x)ζ(y, x)f (x, y)ζ(y, x)g(y, x)) = y,y∼x ζ(x, y)φ(x) + y,y∼x 1 2m(y) (-ζ(y, x)f (x, y)ζ(y, x)g(y, x)) = (z -2m(x))φ(x) + y,y∼x 1 2m(y) (-ζ(y, x)f (x, y)ζ(y, x)g(y, x)) . (x, y)ζ(y, x)g(y, x) + g(x, y) å = (z -2m(x))φ(x) + (y, x)g(y, x) + g(x, y)) = (z -2m(x))φ(x) + 2m(x)φ(x) + (y, x)g(y, x) + g(x, y)) = zφ(x) + y,y∼x1 2m(y)(-ζ(y, x)g(y, x) + g(x, y))

  DP • K • (ζ -1 I -B) = (Az) • (2m 1 ) -1 DP.This is equivalent to the two relationsP • K • (ζ -1 I -B) • P = D -1 (Az) • (2m 1 ) -1 DP and P • K • (ζ -1 I -B) • (I -P) = 0. The latter implies that K • (ζ -1 I -B) sends KerP to itself.We use the decomposition C B = ImP ⊕ KerP. The two relations above tell us thatdet[K • (ζ -1 I -B)] = det[(Az) • (2m 1 ) -1 ] × det[K • (ζ -1 I -B)] KerP-→KerP . But f ∈ KerP is equivalent to Bf = -ιf . Thus if f ∈ KerP, K • (ζ -1 I -B)f (x, y) = (2m(y)) -1 f (x, y) Ç ζ(y, x) -1 ζ(x, y) å = -f (x, y).

  Finally we obtaindet K det(ζ -1 I-B) = det(A-z) x∈V (2m(x)) -1 (-1) dim KerP = det(A-z) x∈V (-G(x)) (-1) dim KerP = det(Az) x∈V (-G(x))(-1) |B|-|V | = det(z -A) x∈V (-G(x)) since |B| = 2|E| is even.To prove the determinant relation, there remains to compute det K. But K is diagonal by blocks of size 2 in the canonical basis of C B . More precisely, denoting by δ e the element of C B that takes the value 1 on e and 0 elsewhere,Kδ e = -12m(y) δ e + ζ(x, y) 2m(x) δ ê if e = (x, y). Thus det K = e={x,y}∈E (2m(x)) -1 (2m(y)) -1 (1ζ(x, y)ζ(y, x)) = e={x,y}∈E (-G(x, y)).

Acknowledgements :

The author is supported by the Labex IRMIA and USIAS of Université de Strasbourg, and by Institut Universitaire de France. This material is based


domain: math phys chim nlin spi

text:
Markus Schmuck
email: m.schmuck@hw.ac.uk

Homogenization of solid-electrolyte composite intercalation cathodes for energy storage systems

Keywords: batteries, homogenization, current-voltage curves, lithium intercalation, phase separation, Butler-Volmer reactions

We investigate well-accepted formulations describing charge transport in composite cathodes of batteries. Our upscaling of carefully selected microscopic equations shows three main features:

(i) a novel set of six equations equipped with nine effective parameters which systematically couple the microscale to the macroscale;

(ii) the coupling of transport and flow equations allows to account for three scales: pore scale, Darcy scale, and macroscale;

(iii) the upscaled equations take phase separation during Li-intercalation into account as well as specific particle configurations.

The wide range of applications and interest in energy storage devices make these results a promising tool to study the influence of the microstructure on current-voltage characteristics and to optimize cathode designs.

Introduction

We look at basic and widely accepted equations describing charge transport and reactions in composite cathode materials of energy storage systems such as Li-batteries. The main goal of this article is to systematically and reliably derive effective macroscopic equations that account for characteristic material properties such as geometric features of interfaces and material phases as well as electric and ionic conductivities and reaction rates. In particular, the subsequent mathematical modelling shall help to gain a deeper understanding into how material properties influence experimentally measured current-voltage characteristics. With the introduction of LiFePO 4 as a promising cathode material in 1997, strong improvements in charging and discharging rates as well as in safety have helped to make phosphates a promising battery material [START_REF] Hautier | Phosphates as lithium-ion battery cathodes: An evaluation based on high-throughput ab initio calculations[END_REF]. A crucial role for the charging and discharging rate plays the so-called Li-intercalation in which Li-ions enter into crystalline phosphate particles of the composite cathode. During this process, the phosphate material can phase separate into a Li-rich phase and a Li-poor phase for which a first modelling seems to be [START_REF] Han | Electrochemical modeling of intercalation processes with phase field models[END_REF]. Previous formulations describe the intercalation by a shrinking core using radial diffusion [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF]. Crystalline intercalation compounds show a topotactic transition, i.e., the compound may change its structure due to gain or loss of material (which is Li here). Recent research suggests that a change in structure/phase separation can be suppressed in LiFePO 4 nanoparticles, e.g. [START_REF] Bazant | Theory of chemical kinetics and charge transfer based on nonequilibrium thermodynamics[END_REF]. Suppression can occur in equilibrium [START_REF] Meethong | Size-dependent lithium miscibility gap in nanoscale li1-x fepo4[END_REF], but when a current flows [START_REF] Malik | Kinetics of non-equilibrium lithium incorporation in LiFePO 4[END_REF], non-idealities are enhanced and hence the suppression seems even more likely. This leads immediately to the important and interesting question what influence can the mircro-geometry have on the appearance of phase separation in these composite cathodes and on the current-voltage characteristics. For instance, it is possible that different parts of the composite cathode are in a different state and hence phase separation can be expected to be a very local property over time. In order to better understand the key properties of composite cathodes for optimizing design and safety, systematic mathematical and thermodynamic modelling will play a crucial role towards a reliable, effective macroscopic description of intercalation kinetics of Li + -ions into solid particles such as FePO 4 .

Existing porous electrode theory for batteries relies on intuitive, formal, or ad hoc averaging strategies based on integrating the equations directly over a representative volume element [START_REF] Bazant | Theory of chemical kinetics and charge transfer based on nonequilibrium thermodynamics[END_REF][START_REF] García | Microstructural modeling and design of rechargeable lithium-ion batteries[END_REF][START_REF] Thomas-Alyea | Electrochemical Systems[END_REF]. This kind of intuitive averaging seems to take dominant transport characteristics such as porosity and surface area into account but generally neglects more subtle geometric features such as pore shape and properties that go beyond porosity and tortuosity characterisations. Currently, it is still an open question how to rigorously justify the volume averaging approach. Often, these methods are refined with physically derived corrections to transport coefficients. Famous porous media corrections to transport parameters are the so-called Bruggeman relation [START_REF] Bruggeman | Berechnung verschiener physikalischer Konstanten von heterogenen Substanzen I. Dielektrizitaätskonstanten und Leitfähigkeiten der Mischköerper aus isotropen Substanzen[END_REF], i.e., the diffusion coefficient D p of the pore phase is obtained by multiplying the free space diffusion parameter D f by the factor p 3/2 where p is the porosity. More refined relations also take the tortuosity into account [START_REF] Van Brakel | ANALYSIS OF DIFFUSION IN MACROPOROUS MEDIA IN TERMS OF A POROSITY , A TORTUOSITY AND A CONSTRICTIVITY FACTOR XA RSA + XB / ( MA CrB ˜[END_REF], i.e., D p = 1/τ D f . The parameter τ is the tortuosity defined as the quotient τ := La Ls where L a is the length of the actual transport path through the pore and L s is the shortest distance connecting entrance and exit of a reference pore. There is increasing interest in the context of ionic transport to systematically derive upscaled/homogenized equations [START_REF] Allaire | Homogenization of the linearized ionic transport equations in rigid periodic porous media[END_REF][START_REF] Kekenes-Huskey | Predicting the influence of longrange molecular interactions on macroscopic-scale diffusion by homogenization of the smoluchowski equation[END_REF][START_REF] Ray | Multiscale Modeling of Colloid and Fluid Dynamics in Porous Media Including an Evolving Microstructure[END_REF][START_REF] Schmuck | Modeling and deriving porous media {S}tokes-{P}oisson-{N}ernst-{P}lanck equations by a multi-scale approach[END_REF][START_REF] Schmuck | First error bounds for the porous media approximation of the Poisson-Nernst-Planck equations[END_REF][START_REF] Schmuck | Homogenization of the poisson-nernst-planck equations for ion transport in charged porous media[END_REF]. For fuel cell and battery systems recent upscaling results are [START_REF] Gully | Effective transport properties of porous electrochemical materials -a homogenization approach[END_REF][START_REF] Lai | Mathematical modeling of porous battery electrodesRevisit of Newman's model[END_REF][START_REF] Richardson | Multiscale modelling and analysis of lithium-ion battery charge and discharge[END_REF][START_REF] Schmuck | Effective macroscopic equations for species transport and reactions in porous catalyst layers[END_REF] for instance. A summary of work closely related to the results presented in this article is given in Tab. 1. Due to increasing interest and number of publications, we do not claim this table to be complete.

Besides a wide interest in thermodynamics and materials science, the highly heterogeneous and multiphase character of composite cathodes requires a systematic and reliable derivation of effective macroscopic transport equations taking physically observed principles as well as geometric features of the microscale into account while still allowing for feasible and efficient computational schemes. The homogenization of charge transport equations in cathodes is the main goal of this article. The basic idea in homogenization is to account for at least one smaller length scale L, the so-called microscale y = x L ∈ Y ⊂ R d , which takes geometric features such as pore shape or changing material properties such as electric permittivity/conductivity into account. The parameter L is the macroscale (e.g. length of the composite cathode) such that = L represents the material's heterogeneity. The microscale is identified by a periodic reference cell (periodic homogenization). Homogenized equations can then be obtained by applying a so-called asymptotic Table 1: Summary of derivations of effective macroscopic equations for battery systems with Li-intercalation. The fact that the majority of work relies on volume averaging can be in part motivated by the fact that this intuitive method is more formal as it does not require rigorous existence and uniqueness results in the derivation like the asymptotic two-scale expansion applied here. Due to an increasingly growing literature on Li-battery modelling over the last decade, we do not claim this list to be complete.

Article

Volume averaging Homogenization Doyle et al. 1993, [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF] • concentrated solution • Nernst-Planck transport • Butler-Volmer reactions • "fast diffusion" intercalation • no flow Bazant 2013, [START_REF] Bazant | Theory of chemical kinetics and charge transfer based on nonequilibrium thermodynamics[END_REF] • "phase field" intercalation Then well-posedness theory allows us to identify the equations for u 0 and u 1 which represent the relevant terms in the upscaled/homogenized problem. This is a well-accepted methodology by now [START_REF] Chechkin | Homogenization: Methods and Applications[END_REF][START_REF] Cioranescu | An Introduction to Homogenization[END_REF][START_REF] Hornung | Homogenization and Porous Media[END_REF][START_REF] Mei | Homogenization Methods for Multiscale Mechanics[END_REF]. We describe the composite cathode of lithium ion batteries together with available charge transport equations in Section 2. Then, we formulate the microscopic composite cathode model which we will subsequently upscale/homogenize and which fully resolves the different phases in Section 2.2. Based on well accepted charge transport models, we present the homogenized/upscaled charge transport equations for composite cathode materials in Section 3. A derivation based on the asymptotic two-scale expansion follows then in Section 3.1.

Microscopic composite cathode formulation

We consider a composite cathode which consists of solid intercalation particles Y s ⊂ Y := [0,1] d , d ∈ N dimension of space, and a polymer solution Y p ⊂ Y which together define a characteristic reference cell Y := Y p ∪ Y s , see Fig. 1 (Left). Moreover, the particles Y s are connected via thin carbon black pathways allowing for electron transport. We state the basic equations that account for charge transport in a composite as defined in the reference cell Y . We refer the reader not interested in an overview on available physical models directly to Section 2.2 where we formulate the microscopic starting point of our upscaling/homogenization derivation.

Basic physical processes and models

Before we look at the upscaling/homogenization of the composite cathode, we introduce models available in the literature to describe charge transport.

(S) Li-Intercalation compound/particles (solid matrix):

The intercalation particles show two different transport characteristics, i.e, electron transport and Li-intercalation into an electrochemical solid (subject to constraints such as charge and mass conservation, see [START_REF] García | Microstructural modeling and design of rechargeable lithium-ion batteries[END_REF] for instance.) The general atomic transport/diffusion mechanisms through solids are either by interstitial diffusion (defects in lattice where atoms assume a normally unoccupied site), vacancy diffusion (driven by thermal energy of atomic vibrations), or via impurities (interdiffusion driven by concentration gradients). Interstitial diffusion is generally faster than vacancy diffusion due to the higher number of interstitial sites and the weaker bonding of interstitials to neighbouring atoms. Since Li resides in interstitial sites, crystalline intercalation hosts such as FePO 4 are favourable for battery applications. We focus here on topotactic compounds, that means, hosts/particles that may change its structure due to gain or loss of material. As a consequence, the local Li concentration can serve as a field variable which differs between a Li-rich and a Li-poor phase separated by a diffuse interface.

(i) Li-interacalation: The thermodynamics of phase separation can be systematically described by phase field models [START_REF] Bazant | Theory of chemical kinetics and charge transfer based on nonequilibrium thermodynamics[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF][START_REF] Elliott | Numerical computation of advection and diffusion on evolving diffuse interfaces[END_REF][START_REF] Han | Electrochemical modeling of intercalation processes with phase field models[END_REF][START_REF] Pego | Front migration in the nonlinear cahn-hilliard equation[END_REF] which rely on a thermodynamic free energy functional

F r (c s ) = n v Ds f r (c s ) + λ(∇c s ) 2 dx, (2.1)
where f is the homogeneous free energy of a regular solution, that is,

f r (c s ) = ωc s (1 -c s ) + k B T (c s log c s + (1 -c s )log(1 -c s )) . (2.
2)

The parameter ω in (2.2) is an effective interaction energy between a Li interstitial site and its neighbouring (Li-free) sites. The parameter λ, which accounts for the strength of penalisation of concentration gradients, is proportional to the interfacial width and assumes in the case of LiFePO 4 as intercalation material the following form,

λ = a 2 ω/6, (2.3)
where ω =59meV is the nearest neighbour interaction energy and a = 2.5 Å the nearest neighbour distance. Using also n v = 8.396 × 10 28 [/m 3 ] then leads to λ = n v λ = 2.48 × 10 -11 J/m as in [START_REF] Han | Electrochemical modeling of intercalation processes with phase field models[END_REF]. Taking a mass conserving gradient flow, i.e., computing the variational derivative with respect to test functions v ∈ H 1 (D s ), we obtain up to boundary terms

∂ t c s = -div M∇ ∇ H 1 cs F r (c s ) , (2.4)
where the derivative ∇ H 1 c represents the Gâteaux derivative. That means, for all v ∈ H 1 (D s ) we compute

δ v F r (c s ) δ v c s = lim θ→0 F r (c s + θv) -F r (c s ) θ = ((f r (c s ) -λ∆c s ),v) = ∇ L 2 cs F r (c s ),v , (2.5)
where the last equality identifies ∇ L 2 cs F r (c s ) by Riesz's theorem. In the same way we can apply Riesz's theorem with respect to the semi-inner product M∇u,∇v

H 1
for u,v ∈ H 1 (D s ) and a symmetric positive definite tensor M to define ∇ H 1 cs F r (c s ) under homogeneous Neumann boundary conditions as follows

∇ H 1 cs F r (c s ),v = M∇∇ L 2 cs F r (c s ),∇v = -div M∇ L 2 cs F r (c s ) ,v . (2.6)
Putting (2.12)-(2.6) together, leads with f r (c s ) = ω(1 -2c s ) + k B T log cs 1-cs to the phase field equation for intercalation compounds

∂ t c s = div M∇ f r (c s ) -λ∆c s , (2.7)
where the mobility tensor M is often chosen as a constant M s since generally it is difficult to trace back the irreversible behaviour to the underlying frictional processes after coarse graining. A state equation with concentration dependent mobility is suggested in [START_REF] Han | Electrochemical modeling of intercalation processes with phase field models[END_REF], i.e., M = c s (1 -c s )M 0 where M 0 is a constant.

(ii) Electron transport: We describe charge transport in the solid matrix by the current

i s = -σ s ∇ψ s , (2.8)
where σ s is the conductivity of the solid intercalation phase D s and ψ s is the electric potential in D s . Equation (2.8) is widely used in the literature, e.g. [START_REF] Thomas-Alyea | Electrochemical Systems[END_REF].

(E) Electrolyte/separator: A large part of the electrochemical literature is dedicated to transport equations for binary electrolytes. These denote solutions containing a single binary salt of the form X 1 X 2 plus a solvent (generally polymer in Li-batteries) where X + 1 and X - 2 represent a cation and an anion, respectively. Electrolytes X 1 X 2 that dissociate into an equal number of cations ν + and anions ν -are called symmetric. Subsequently, we present four different formulations (E1)-(E4) for binary electrolytes.

(E1) Dilute solution theory: charge transport by Poisson-Nernst-Planck (PNP). There is a wide range of different modelling approaches as well as analytical studies of charge transport formulations such as [START_REF] Eisenberg | Energy variational analysis of ions in water and channels: Field theory for primitive models of complex ionic fluids[END_REF][START_REF] Jerome | Analytical approaches to charge transport in a moving medium[END_REF][START_REF] Li | Solutions to a reduced poisson-nernst-planck system and determination of reaction rates[END_REF][START_REF] Markowich | Semiconductor equations[END_REF][START_REF] Roubicek | Incompressible ionized non-newtonian fluid mixtures[END_REF][START_REF] Schmuck | Analysis of the NavierStokesNernstPlanckPoisson System[END_REF][START_REF] Schmuck | Homogenization of the poisson-nernst-planck equations for ion transport in charged porous media[END_REF]. Recently, the PNP system has been modified to describe electrolytes beyond the dilute approximation in [START_REF] Dreyer | Overcoming the shortcomings of the Nernst-Planck model[END_REF], see also Remark 2.1 below. The widely used and well-accepted Nernst-Planck formulation allows us to take the densities c i of each species i in the solution into account. Since we consider a polymer electrolyte as a solution, we have a neutral (z 0 = 0) density c 0 for the polymer solvent, the cation density c + = c 1 of lithium ions X + 1 , and the anion density c -= c 2 of negative salt ion X - 2 . The following classical and well-accepted description is based on the PNP system, (PNP)

∂ t c i + u • ∇c i = div (D i ∇c i + ez i c i M i ∇ψ) , for i = 0,1,2 , -div(ε∇ψ) = e 2 i=0 z i c i , (2.9)
where M i , D i , and z i are the diffusion, mobility, and charge number of species i, respectively. The parameter ε is the electric permittivity of the electrolyte and e denotes the elementary charge.

The fluid velocity u is the solution of the incompressible Navier-Stokes (iNS) equations extended by a Coulomb force term on the right-hand side [START_REF] Jerome | Analytical approaches to charge transport in a moving medium[END_REF][START_REF] Roubicek | Incompressible ionized non-newtonian fluid mixtures[END_REF][START_REF] Schmuck | Analysis of the NavierStokesNernstPlanckPoisson System[END_REF], i.e., (iNS)

∂ t u + u • ∇u -∆u + ∇p = -e 2 i=0 z i c i ∇ψ , divu = 0.
(2.10)

In the context of the dilute solution theory, see Remark 2.1 below for more details, one does not use the flux of the solvent, i.e., i = 0 in (2.9), but instead the solvent is solely described by the momentum equation (2.10). For simplicity, we will subsequently denote the cation density by c + , the anion density by c -.

(E2) Binary electrolyte [START_REF] Thomas-Alyea | Electrochemical Systems[END_REF]. Since bulk electroneutrality is a widely confirmed characteristic property of solutions, we make the frequent assumption of an electrically neutral solution, i.e., i z i c i = 0.

(2.11)

We note that the neutrality assumption (2.11) is violated in interfacial regions (e.g. near electrodes) where thin electrical double layers (EDL) appear. Under electroneutrality, the concentration of the electrolyte is consistently defined by

c := c + ν + = c - ν - .
(2.12)

Moreover, we assume constant mobilities and diffusion coefficients. Applying the fluxes j i , i = 1,2, from the Nernst-Planck setting (2.9) leads to

∂ t c + u • ∇c = D∆c, (z + M + -z -M -)F div(c∇ψ) + (D + -D -)∆c = 0, (2.13)
where

D = z + M + D --z -M -D + z + M + -z -M - . Remark 2.1 (Validity)
The dilute solution theory is recommended as a general formulation since it is widely used in the electrochemical literature and at the same time it accounts for the essential physical processes involved without very complicated/less accessible ingredients [48, p.274]. But one should be aware that the fluxes defined in (2.9) 1 break down because electro-migration and diffusion are not consistently defined with respect to an average velocity u of the fluid. A way out of this is a so-called dilute solution where one does not apply (2.9) to the solvent, i.e., i = 0, but solely to the ionic species and where one attributes the fluid velocity u to the velocity of the solvent.

Recently, the authors of [START_REF] Dreyer | Overcoming the shortcomings of the Nernst-Planck model[END_REF] have carefully specified the fluid velocity u in order to overcome the above noted difficulty.

(E3) Concentrated solution. We state a charge transport model for a polymer electrolyte as proposed in [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF] for composite cathodes. This transport formulation relies on experimentally measurable quantities such as the transference number. We still assume that the polymer electrolyte is a solution of a binary salt and a solvent. For a Li-polymer insertion cell the solution is for instance a binary salt LiX and the solvent a polymer. Taking the polymer as a reference species, the following material balance on the lithium salt electrolyte c = c i /ν i ,

∂ t c = div D(c) 1 - d(lnc 0 ) d(lnc) ∇c - i p • ∇t 0 + (c) z + ν + F , (2.14)
where

t j = z 2 j M j c j i z 2
i M i c i is the experimentally obtained transference number which accounts for the current associated with species j. In [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF], the following fit to experimental data is employed t 0 + = 0.0107907 + 1.48837 × 10 -4 c.

(2.15)

Following again the concentrated solution theory [START_REF] Thomas-Alyea | Electrochemical Systems[END_REF] applied to lithium-polymer insertion cells in [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF], the current i p in the polymer phase admits the following representation

i p := -κ(c)∇ψ p - κ(c)RT F 1 + ∂ ln f ± ∂ ln c s + nν + + t 0 + (c) z + ν + ∇ln c, (2.16)
where the electric potential ψ p in the polymer phase is measured with a reference electrode, κ = F i z 2 i M i c i is the conductivity of the solution not accounting for concentration gradients. The variable f ± is the activity coefficient of the salt and is only accounted for if data are available [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF].

(E4) Onsager-type equations away from equilibrium. The subsequent charge transport formulation has been introduced in [START_REF] Efendiev | Upscaling of an isothermal li-ion battery model via the homogenization theory[END_REF]. The concentration of lithium c + in the electrolyte phase is governed by the following system

∂ t c + = div k + 11 (c + )∇c + + k + 12 ∇ψ = 0, 0 = div k + 21 (c + )∇c + + k + 22 ∇ψ , (2.17)
where the coefficients are defined for the ion conductivity κ + and ion diffusion coefficient

D + by k + 11 (c + ) = D + + RT F 2 (t + ) 2 κ + c + , k + 12 = κ + t + F , k + 21 (c + ) = RT F κ + t + c + , k + 22 = κ + .
(2.18)

Remark 2.2 (Comments to formulations (E1)-(E4))

The formulations (E1)-(E2) are widely used in theory and computational modelling whereas (E3) is generally used in experimental science in view of particular parameters that represent directly measurable quantities. The transport equations (E4) are interesting as they resemble in their stationary form the classical Onsager relations [START_REF] Onsager | Reciprocal relations in irreversible processes I[END_REF][START_REF] Onsager | Reciprocal relations in irreversible processes II[END_REF].

(I) Solid-electrolyte interphase: The interfacial transition of Li from the polymer solution into the solid insertion compound is governed by the Butler-Volmer kinetics [START_REF] Bard | Electrochemical Methods: Fundamentals and Applications[END_REF][START_REF] Bedeaux | Nonlinear coupled equations for electrochemical cells as developed by the general equation for nonequilibrium reversible-irreversible coupling[END_REF]. Following Doyle et al. [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF], the following reaction

Li + -Θ p + Θ s + e - Li -Θ s + Θ p , (2.19)
governs the Li insertion into the solid matrix. The variables Θ p and Θ s stand for the available sites in the polymer and the solid matrix (interstitial). For an overpotential η := ψ s -ψ between the solid (ψ s ) and the bulk electrolyte (ψ) and the voltage

U = U θ 2 -U θ ref + RT /F (βc s + ζ), the insertion process (2.19) induces the current i = F k fb (c pm -c) αc c αa c s exp α a F RT (η -U ) -(c sm -c s )exp - α c F RT (η -U ) , (2.20)
where c sm and c pm are the maximal densities in the solid and the polymer phase, respectively. The variable ζ is a constant expressing corrections to activity. We will set these parameters according to [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF], i.e., c sm = 29000mol/m 3 , c pm = 3920mol/m 3 . Additionally, the transfer coefficients are α a = α c 0.5, and the reaction rate at the particle-electrolyte interface is k fb = 10 -10 m 4 /(mol• s).

     i = -s i nF j p • n on ∂D p , i = -j s • n on ∂D s , i = -σ s ∇φ s on ∂D s , (2.21)
where n is the outward normal of the corresponding domains and the number of electrons is one, that is, n = 1 = s i , where s i is the stochiometric parameter.

Microscopic equations for composite cathodes

Based on the previous Section 2.1, we can now state the microscopic system of equations describing charge transport in a periodically extended composite cathode, see Fig. 2. We distinguish two different microscales, that is, a reference pore scale of characteristic length l p and a so-called Darcy scale of characteristic length l D . Herewith, mass transport is governed by the heterogeneity p = lp L of pores where l p is the characteristic size of a reference particle configuration and L is the macroscopic length of the composite cathode. The microscopic formulation fully resolves the highly heterogeneous geometry of the composite cathode in which we have to account for three essential transport mechanisms: The following equations fully describe the microscopic kinetic processes in the composite cathode and hence are indexed by the pore heterogeneity p . Moreover, we allow for a further length scale l D > l p , called Darcy scale, see Fig. 1 Right. This then leads to a so-called Darcy's heterogeneity D := l D L . We note that the subsequent results are valid for scenarios where the pore scale p differs from the Darcy scale D , i.e., where p < D . Hence, the scaling p < D enables us to account for a third scale via the momentum equation. Hence, fluid flow is considered on a perforated domain

D D p := z∈Z d D (Y p D + z) ∩ D
which is defined by a slightly larger representative volume element Y p D accounting for the Darcy scale l D . As the boundary conditions on ∂D do not have a direct influence in the upscaling inside the composite cathode D, we do not state here how D interacts with external forces and systems as this can still be made precise at a later stage.

Due to confined material structure, one can expect a fast relaxation time for the flow as well as a small characteristic flow field. Therefore, the flow is governed by a small Reynolds number Re =

ρU ref L µ
(where ρ is the density of the fluid, U ref is a reference velocity, µ the dynamic viscosity, and L a characteristic length of the system) and hence instead of using the full time-dependent Navier-Stokes equations to define u D , we can apply a quasi-stationary formulation based on Stokes equations, i.e.,

-2 D µ∆u D + ∇p D = f p := -e(z + c p + -z -c p -)∇ψ p p in D D p , divu D = 0 in D D p , (2.22)
where -e(z + c p + -z -c p -)∇ψ p p is a Coulomb driving force [START_REF] Jerome | Analytical approaches to charge transport in a moving medium[END_REF][START_REF] Roubicek | Incompressible ionized non-newtonian fluid mixtures[END_REF][START_REF] Schmuck | Analysis of the NavierStokesNernstPlanckPoisson System[END_REF]. In the polymer electrolyte (binder), we have

     ∂ t c p + + u p • ∇c p + = div k B T M + ∇c p + + ez + c p + ∇ψ p p in D p p , ∂ t c p -+ u p • ∇c p -= div k B T M -∇c p --ez -c p -∇ψ p p in D p p , -div ε p ∇ψ p p = e(z + c p + -z -c p -) in D p p , (2.23)
where the parameter ε p denotes the electric permittivity of the polymer and hence differs from the pore heterogeneity p . The variable c = c + denotes the density of Li + -ions. The mean field equations describing interstitial and electron transport inside a solid crystal read in the pore scale as follows

∂ t c p s = div M s ∇ f r (c p s ) -λ∆c p s in D p s , -div σ s ∇φ p s = 0 in D p s , (2.24)
where the electric conductivity of the solid intercalation compound is σ s . We note that the Poisson equations (2.33) 3 and (2.24) 2 can be written as a single equation for the electric potential φ p by

-div(ε(x/ p )∇ψ p ) = e(z + c p + z -c p -)χ D p p in D , (2.25)
valid on the whole domain D for the combined permeability and conductivity parameter

ε(x/ ) := ε p χ Yp (x/ p ) + σ s χ Ys (x/ p ), (2.26)
and with the characteristic function χ D p p (x) = 1 if x ∈ D p p . In (2.26), we employ the microscopic variable y ∈ Y defined on the reference cell Y , see Fig. 1 for a sample polymer-particle configuration of a composite cathode.

In order to complete the microscopic formulation, we still need to state the interfacial processes taking place on the polymer/binder-solid interface I p , i.e.,

                   p R p (c p + ,c p s ,η p ,U p ) = RT F M + ∇c p + + ez + c p + ∇ψ p -u D c p + • n ps on I p , 0 = RT F M -∇c p --ez -c p -∇ψ p -u D c p -• n ps on I p , p R (c p + ,c p s ,η p ,U p ) = -M s ∇ f r (c p s ) -λ∆c p s • n sp on I p , p R p (c p ,c p s ,η p ,U p ) = σ s ∇ψ s • n sp on I p , σ s ∇ψ p p • n sp = ε p ∇ψ p s • n sp on I p , u D = 0 on I p , (2.27)
where n sp = -n ps are normals pointing outward from the solid into the polymer or from the polymer into the solid phase, respectively. The current R p = R p (c p ,c p s ,η p ,U p ) is defined by the Butler-Volmer kinetics (2.20) extended here towards the highly heterogeneous composite cathode setting indicated by the pore heterogeneity p , that is,

R p = i ps c p s c p sm exp α sm α a F RT (η p -U p ) - (c p sm -c p s ) c p sm exp -α sm α c F RT (η p -U p ) , (2.28)
where i ps = F k fb (c p pm -c p ) αc (c p ) αa plays the role of an exchange current density and c p sm = exp(α sm ) brings the term in the angled brackets into dimensionless form. For the limits lim x→I p s and lim x→I p p , which denote passing from the solid phase and the polymer phase to the interface I p , respectively, we can identify the overpotential in (2.28) by

η p = lim x→I p s ψ p -lim x→I p p ψ p .
(2.29) Equation (2.28) additionally depends on the maximum capacity of ions in the polymer phase c p pm = 3920mol/m 3 as well as in the solid crystal c p sm = 29000mol/m 3 . We also recall from equation (2.20) that a fitting to experimental data [START_REF] Doyle | Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell[END_REF][START_REF] West | Modeling of porous insertion electrodes with liquid electrolyte[END_REF] denotes the dimensionless Debye length and σ := σs σ ref the dimensionless electric conductivity. We set the reference conductivity to the bulk conductivity in the electrolyte phase, i.e., σ ref = εpD λ 2 . For notational convenience, we assume M = M + = M -, and hence D = D + = D -, and we will subsequently drop the tilde from dimensionless variables. With Einstein's relation we then have D = k B T M . Moreover, we note that the convection-diffusion-electromigration equation reads with the above scalings as follows

D = u D U ref , c p ± = c p ± c ref , ψ = F ψ RT , t = t τ , x = x L , η = 2c ref RT L U ref µ , Pe = LU ref D , (2.31
L 2 τ D ∂c p ± ∂ t + Pe ũ D ∇c p ± = div z ± c p ± ∇ ψ .
As noted earlier in this section, in the composite cathode one expects a small Re = Pe/Pr but a large Prandtl number Pr := ν D , e.g. in solute transport in groundwater flow D = 10 -5 cm 2 /s, ν = 10 -2 cm 2 /s, and hence Pr = 10 3 , where ν = µ ρ is the the kinematic viscosity and ρ the density. Therefore, we allow for Pe = O(1) [START_REF] Mei | Homogenization Methods for Multiscale Mechanics[END_REF].

Herewith and after choosing τ = τ u , and setting τs τu = 1 for simplicity, we obtain the following dimensionless microscopic formulation

(micro (SE)) -2 D ∆u D + ∇p D = f p := -η(z + c p + -z -c p -)∇ψ p in D D p , divu D = 0 in D D p . (2.32)
(micro (E1))

Pe ∂ t c p + + u p • ∇c p + = div ∇c p + + z + c p + ∇ψ p in D p p , Pe ∂ t c p -+ u p • ∇c p -= div ∇c p --z -c p -∇ψ p in D p p , (2.33)
(micro (S))

∂ t c p s = div ∇ f r (c p s ) -λ∆c p s in D p s , (2.34)
(micro (P)) -div λ D ∇ψ p p = (z

+ c p + + z -c p -) in D p p , -div σ∇ψ p s = 0 in D p s , (2.35)
(micro (I))

                     p β + R p (c p + ,c p s ,η p ,U p ) = ∇c p + + z + c p + ∇ψ p -u D c p + • n ps on I p , 0 = ∇c p --z -c p -∇ψ p -u D c p -• n ps on I p , p β s R (c p + ,c p s ,η p ,U p ) = -∇ f r (c p s ) -λ∆c p s • n sp on I p , p β ψ R p (c p + ,c p s ,η p ,U p ) = σ∇ψ s • n sp on I p , σs σ ref ∇ψ p p • n sp = εp εs ∇ψ p s • n sp on I p , u D = 0 on I p . ( 2

Main results: Effective macroscopic equations

Our main result is the reliable and systematic derivation of a novel set of effective macroscopic composite cathode equations that take the crucial chemical and transport processes described in Section 2.2 into account: (a) ionic transport; (b) interfacial reactions; (c) intercalation dynamics under a possible phase transformation; (d) electron transport.

The subsequent upscaling result relies on the following property. With the help of the asymptotic two-scale expansion method, we derive the following Main result: (Upscaled composite cathode equations) The microscopic composite cathode formulation based on the equations (2.32), (2.33), (2.34), and (2.36) admits under local thermodynamic equilibrium (Definition 3.1) the following upscaled/homogenized equations

Definition 3.1 (Local Thermodynamic Equilibrium (LTE)) The chemical potentials µ ± (C ± ,Ψ) = logC ± ± Ψ and µ s (C s ) = f r (C s ) -λ q div M∇C s are said to be in local thermodynamic equilibrium if for every k ∈ N, 1 ≤ k ≤ d, ∂µ ι (C ι ,Ψ) ∂x k = 0 if it appears in the cell problem in Y × D , ∂µ ι (C ι ,Ψ) ∂x k if it
                             V(x,t) = -κ µ (f 0 -∇P ) in D , divV = 0, in D , p∂ t C + + PeV • ∇C + = div D∇C + + z + C + M∇Ψ p -β + R(C + ,C s ,η 0 ,U 0 ) in D , p∂ t C -+ PeV • ∇C -= div D∇C --z -C -M∇Ψ p in D , -div(ε∇Ψ p ) = p(z + C + + z -C -) in D , -div Σ∇Ψ s = β ψ R(C + ,C s ,η 0 ,U 0 ) in D , q∂ t C s = div Ms ∇f r (C s ) + λ q div Ms ∇C s + β s R(C + ,C s ,η 0 ,U 0 ) in D , (3.39
)

where p = |Y 1 | |Y | , q = 1 -p, f 0 = -η (C + -C -) M∇Ψ p , β + = i 0 L|Λ| eD , β ψ = F RT ipsL|Λ| σ ref , β s = i 0 L|Λ| eMs , η 0 = Ψ s - Ψ p , U 0 = 2.17 + RT F -0.000558C s I p + 8.

, and the material's correction tensors

κ = {κ ij } d i,j=1 , M = {m ij } d i,j=1 , D = d ij d i,j=1 , ε = {ε ij } d i,j=1 , Σ = {σ ij } d i,j=1
, and Ms = m s ij d i,j=1 are defined by

κ ij := - 1 µ|Y D | Y p D w j i (y)dy, ε ik = 1 |Y | Yp λ D δ ik -δ ij ∂ξ k ψ (y) ∂y j dy, m ik := 1 |Y | Yp δ ij (y) δ ik -δ ij ∂ξ k ψ (y) ∂y j dy, m s ik = 1 |Y | Ys δ ij (y) δ ik -δ ij ∂ξ k φ (y) ∂y j dy, d ik := 1 |Y | Yp d ij (y) δ ik -δ ij ∂ξ k ± (y) ∂y j dy, σ ik = 1 |Y | Ys σ δ ik -δ ij ∂ξ k ψ (y)
∂y j dy.

(3.40)

The correctors w j i , ξ k ψ , ξ k φ , and ξ k solve the reference cell problems (3.48), (3.63) 2 , (3.72), and (3.63) 1 , respectively.

Before we derive these effective macroscopic formulation, we provide values for the parameters used in general Li batteries, see Tab. 3.1 Derivation of the upscaled cathode equations (3.39).

As the upscaling involves several equations, we will start with the Stokes equations driven by a Coulomb force f p . Then we consider the PNP system, and finally we look at the phase field equation describing the Li-intercalation process.

Homogenizing the generalized Stokes equations (2.32) 1 : Our derivation gives an extension of the classical formal homogenization result [START_REF] Hornung | Homogenization and Porous Media[END_REF] towards Coulomb driving forces [START_REF] Jerome | Analytical approaches to charge transport in a moving medium[END_REF][START_REF] Roubicek | Incompressible ionized non-newtonian fluid mixtures[END_REF][START_REF] Schmuck | Analysis of the NavierStokesNernstPlanckPoisson System[END_REF][START_REF] Schmuck | Modeling and deriving porous media {S}tokes-{P}oisson-{N}ernst-{P}lanck equations by a multi-scale approach[END_REF].

Step 1: (Asymptotic expansion) We make the ansatz

u D (x) = U(x,y D ) + D u 1 (x,y D ) + D 2 u 2 (x,y D ) + ••• , p D (x) = P (x,y D ) + D p 1 (x,y D ) + D 2 p 2 (x,y D ) + ••• . (3.41)
The microscale We insert (3.41) into (2.32). This leads to the following problem where we only state the leading order terms,

y D = x/ D ∈ Y D is Y D -periodic
-D -1 ∇ y P (x,y D ) + D 0 µ∆ y D U(x,y D ) -∇ y D p 1 (x,y D ) + ∇ x P (x,y D ) + D 1 (...) + h.o.t. = f p = -e z + c p + -z -c p -∇φ p p , in Y p D , (3.42)
D -1 div y D U(x,y D ) + D 0 div y D u 1 (x,y D ) + div x U(x,y D ) + D 1 (...) + h.o.t. = 0, in Y p D , (3.43)
D 0 U(x,y D ) + D u 1 (x,y D ) + h.o.t. = 0 on I 0 . (3.44)
The variable I 0 := ∂Y p D ∩ ∂Y p D denotes the solid-pore interface on the level of the referenc cell Y D . Higher order terms are referred to by h.o.t..

Step 2: (Comparing terms of equal order) The lowest order problem in Step 1 is obtained by comparing coefficients of the terms of order D -1 , i.e., we obtain the equations

∇ y D P (x,y D ) = 0, in Y p D , div y D U(x,y D ) = 0, in Y p D . (3.45)
This allows us to conclude that the leading order pressure component P is independent of the microscale y D but not for leading order fluid velocity U.

At the next order O( D 0 ), we obtain

-µ∆ y D U(x,y D ) + ∇ y D p 1 (x,y D ) + ∇ x P (x,y D ) = f p = -e z + c p + -z -c p -∇φ p p , in Y p D , div y D u 1 (x,y D ) + div x U(x,y D ) = 0, in Y p D . (3.46)
Step 3: (Ansatz for first order terms) This step is not entirely standard and requires physical intuition. Making the usual ansatz for the pressure p 1 in the asymptotic expansion and treating the the fluid velocity U instead of u 1 in the same way, which is motivated by its dependence on the microscale y D , i.e.,

p 1 (x,y D ) = d j=1 f p j - ∂P (x) ∂x j π j (y D ), U(x,y D ) = 1 µ d j=1 f p j - ∂P (x) ∂x j w j (y D ), (3.47)
then problem (3.46) 1 turns into the following

(Cell Problem)      -∆ y D w j (y D ) + ∇ y D π j (y D ) = e j , in Y p D , div y D w j (y D ) = 0, in Y p D , w j (y D ) = 0, on I 0 , (3.48)
via ansatz (3.47) 2 as well as into Darcy's law after integrating (3.47) 2 over the pore space Y p D , i.e.,

V(x,t) = |Y p D | |Y D | M Y p D (U(x,y D )) = - 1 µ κ(f p -∇ x P (x)) , (3.49)
where the hydraulic conductivity/permeability tensor κ := {κ ij } d i,j=1 is defined by the volume average

κ ij = - 1 µ 1 |Y D | Y p D w j i (y D )dy D . (3.50)
It leaves to show that the effective Darcy velocity V is divergence free, that means, that it satisfies (2.32) 2 . To this end, we can make use of equation (3.46) 2 , which leads after integration over the pore space Y p D and with the help of Gauss' theorem to

div x V(x) = ∂Y D n • u 1 (x,y D )do(y D ) - I 0 n • u 1 (x,y D )do(y D ) = 0, (3.51)
where we used the Y D -periodicity and the homogeneous boundary condition (3.44) at order O( D ). Finally, we note that once we have the upscaling result for the transport equations (2.33), which are formulated with respect to the pore scale l p , then we can also determine the upscaled form of the Coulomb force term f p by f 0 , i.e., by formally passing p → 0.

2

Homogenizing the PNP system: The following derivation relies on homogenization results for the PNP system obtained in the context of microfluidics and soil mechanics such as [START_REF] Schmuck | Modeling and deriving porous media {S}tokes-{P}oisson-{N}ernst-{P}lanck equations by a multi-scale approach[END_REF][START_REF] Schmuck | First error bounds for the porous media approximation of the Poisson-Nernst-Planck equations[END_REF][START_REF] Schmuck | Homogenization of the poisson-nernst-planck equations for ion transport in charged porous media[END_REF]. The subsequent statements extend the derviation in [START_REF] Schmuck | First error bounds for the porous media approximation of the Poisson-Nernst-Planck equations[END_REF] to include fluid flow.

Step 1: (Asymptotic expansion) We consider the asymptotic two-scale expansions

c ± p (t,x) = C ± (t,x,y) + p c ± 1 (t,x,y) + 2 p c ± 2 (t,x,y) + h.o.t., ψ p (t,x) = Ψ(t,x,y) + p ψ 1 (t,x,y) + 2 p ψ 2 (t,x,y) + h.o.t.. (3.52)
For notational convenience, we introduce for smooth enough functions v p the following operators,

A s v p (x) := - N i,j=1 ∂ ∂x i ε(x/ p )δ ij ∂v p ∂x j = p -2 A 0 + p -1 A 1 + A 2 v p (x,x/ p ), B ± s v p (x) := - N i,j=1 ∂ ∂x i δ ij ∂v p ∂x j = p -2 B ± 0 + p -1 B ± 1 + B ± 2 v p (x,x/ p ), B 3 s v p (x) := - N i,j=1 ∂ ∂x i v p δ ij ∂v p ∂x j = p -2 B 3 0 + p -1 B 3 1 + B 3 2 v p (x,x/ p ), (3.53)
where A i are related to the Poisson equation and are defined by,

A 0 := - d i,j=1 ∂ ∂y i ε(y)δ ij ∂ ∂y j , A 1 := - d i,j=1 ∂ ∂x i ε(y)δ ij ∂ ∂y j - d i,j=1 ∂ ∂y i ε(y)δ ij ∂ ∂x j , A 2 := - d i,j=1 ∂ ∂x i ε(y)δ ij ∂ ∂x j , (3.54)
The operators B ± i and B 3 i related to the Nernst-Planck equations are defined for c ± by,

B ± 0 := - d i,j=1 ∂ ∂y i δ ij ∂ ∂y j , B ± 1 := - d i,j=1 ∂ ∂x i δ ij ∂ ∂y j - d i,j=1 ∂ ∂y i δ ij ∂ ∂x j , B ± 2 := - d i,j=1 ∂ ∂x i δ ij ∂ ∂x j , B 3 0 := -z ± d i,j=1 ∂ ∂y i δ ij ∂Ψ ∂y j , B 3 1 := -z ± d i,j=1 ∂ ∂x i δ ij ∂Ψ ∂y j -z ± d i,j=1 ∂ ∂y i δ ij ∂Ψ ∂x j -z ± d i,j=1 ∂ ∂y i δ ij ∂ψ 1 ∂y j , B 3 2 := - d i,j=1 ∂ ∂x i δ ij ∂Ψ ∂x j -z ± d i,j=1 ∂ ∂x i δ ij ∂ψ 1 ∂y j -z ± d i,j=1 ∂ ∂y i δ ij ∂ψ 1 ∂x j -z ± d i,j=1 ∂ ∂y i δ ij ∂ψ 2 ∂y j . (3.55)
The definitions (3.53), (3.54), and (3.55) allow to immediately obtain the following sequence of problems by equating terms of equal power in p , i.e.,

O( p -2 ) :          B ± 0 C ± + B 3 0 C ± = 0 in Y p , C ± is Y p -periodic, A 0 Ψ = 0 in Y , Ψ is Y -periodic, (3.56) O( p -1 ) :          B ± 0 + B 3 0 c ± 1 + Pediv y (u D C ± ) = -B ± 1 + B 3 1 C ± in Y p , c ± 1 is Y p -periodic, A 0 ψ 1 = -A 1 Ψ in Y , ψ 1 is Y -periodic, (3.57) O(1) :                B ± 0 + B 3 0 c ± 2 + Pediv y u D c ± 1 = -B ± 1 + B 3 1 c ± 1 -B ± 2 + B 3 2 C ± -∂ t C ± -Pediv x (u D C ± ) in Y p , c ± 2 is Y p -periodic, A 0 ψ 3 2 = (C + -C -) -A 1 ψ 1 -A 2 Ψ in Y , ψ 2 is Y -periodic.
     -∆ y C + -div y (C + ∇ y Ψ) = 0 in Y p , -∆ y C --div y (C -∇ y Ψ) = 0 in Y p , -div y (ε(y)∇ y Ψ) = 0 in Y . ( 3
     -∆ y c ± 1 -div y (z ± C ± ∇ y ψ 1 ) + Pediv y (u D C ± ) = div y z ± c ± 1 ∇ y Ψ + div x ∇ y C ± +div x (z ± C ± ∇ y Ψ) + div y ∇ x C ± + div y (z ± C ± ∇ x Ψ) in Y p , -div y (ε(y)∇ y ψ 1 ) = div x (ε(y)∇ y Ψ) + div y (ε(y)∇ x Ψ) in Y , (3.60)
where the convective terms Pediv y (u D C ± ) disappear since C ± and u D are both independent of y. The independence of C ± and Ψ from y together with the linearity of (3.60) motivates to make the following ansatz,

c ± 1 (t,x,y) = - d j=1 ξ j ± (y)
∂C ± (t,x) ∂x j , and ψ 1 (t,x,y) = -

d j=1 ξ j ψ (y) ∂Ψ(t,x) ∂x j ,. (3.61)
With the local thermodynamic equilibrium property, i.e.,

0 = ∂ ∂x j logC ± ± z ± ∂ ∂x j Ψ, (3.62)
we can simplify (3.60) to, The last problem (3.58) is again well-posed by Lax-Milgram's theorem or the Fredholm alternative. A necessary and sufficient condition for the existence and uniqueness of solutions to (3.63) is the requirement that corresponding right-hand sides are zero after integrating over Y p . Hence, we immediately obtain, porous media PNP:

             -N i,j=1 ∂ ∂y i δ ij ∂ ∂y j ξ j ± (y) -y i = -N i,j=1 ∂ ∂y i δ ij ∂ ∂y j ξ j ψ (y) -y i in Y p , ξ j ± (y) is Y p -periodic and M Yp (ξ j ± ) = 0, -N i,j=1 ∂ ∂y i ε(y)δ ij ∂ ∂y j ξ j ψ (y) -y i = 0 in Y , ξ j ψ (y) is Y -periodic and M Y (ξ j ψ ) = 0.
           |Y p |∂ t C ± + Pe |Y p |div x (u D C ± ) -d i,j,k=1 ∂ ∂x i δ ik -δ ij ∂ξ k ± ∂y j ∂C ± ∂x k ,1 Yp -d i,j,k=1 ∂ ∂x i z ± C ± δ ik -δ ij ∂ξ k ψ ∂y j ∂Ψ ∂x k ,1 Yp = 0, -d i,j,k=1 ∂ ∂x i ε(y) δ ik -δ ij ∂ξ k ψ ∂y j ∂Ψ ∂x k ,1 Y = |Y p |(C + -C -) . (3.64)
We note that u D is not upscaled/homogenized here since this allows us to account for multiple scales such as the pore scale l p , the Darcy scale l D , and the macroscale L. For details regarding the verifaction of the cell problems with Lax-Milgram's theorem or the Fredholm alternative [START_REF] Cioranescu | An Introduction to Homogenization[END_REF], we refer the interested reader also to [START_REF] Schmuck | New porous medium Poisson-Nernst-Planck equations for strongly oscillating electric potentials[END_REF] where these tools are applied in a related problem by a modified asymptotic expansion.

It leaves to homogenize the right-hand side of the Stokes equations (2.32), which is scaled with respect to the pore scale l p . With the corrector ξ k ψ we immediately obtain at the level O( 0 p ) that the effective Coulomb force term is

f 0 = -η (C + -C -) M∇Ψ p since at O( 0 p ) we have with (3.61) -η (C + -C -)(∇ x Ψ p + ∇ y ψ 1 ) = -η (C + -C -) d i,j,k=1 δ ik -δ ij ∂ξ k ψ ∂y j ∂Ψ p ∂x k . 2
Homogenizing the phase field intercalation equation Step 1: (Asymptotic expansion) The derivation of the phase field equations follows the essential steps presented in [START_REF] Schmuck | Effective macroscopic interfacial transport equations in strongly heterogeneous environments for general homogeneous free energies[END_REF] and which we repeat here for the reader's convenience. The heterogeneity parameter = p is the pore scale as for the ionic transport in the polymer solution.

In the context of asymptotic two-scale expansions, the Laplace operator ∆ can be written as follows,

A 0 = -d i,j=1 ∂ ∂y i δ ij ∂ ∂y j , A 1 = -d i,j=1 ∂ ∂x i δ ij ∂ ∂y j + ∂ ∂y i δ ij ∂ ∂x j , A 2 = -d i,j=1 ∂ ∂x j δ ij ∂ ∂x j , (3.65)
such that we can identify A := -2 A 0 + -1 A 1 + A 2 = ∆. We account for the multiscale nature of strongly heterogeneous environments [START_REF] Schmuck | First error bounds for the porous media approximation of the Poisson-Nernst-Planck equations[END_REF][START_REF] Hornung | Homogenization and Porous Media[END_REF][START_REF] Mei | Homogenization Methods for Multiscale Mechanics[END_REF] by the following ansatz η ≈ η 0 (x,y,t) + η 1 (x,y,t) + 2 η 2 (x,y,t) + ... , for η ∈ {w,φ} .

(3.66)

Before we can insert (3.66) into the microscopic formulation (2.34), we approximate the derivative of the nonlinear homogeneous free energy density f r by a Taylor expansion around the homogenized order parameter φ 0 , that means,

f r (φ ) ≈ f r (φ 0 ) + f r (φ 0 )(φ -φ 0 ) + 1 2 f r (φ 0 )(φ -φ 0 ) 2 + O (φ -φ 0 ) 3 , (3.67)
where φ 0 denotes the leading-order term in (3.66). Substituting (3.66) and (3.67) into (2.34), which we split into two second order problems as suggested in [START_REF] Schmuck | Effective macroscopic equations for species transport and reactions in porous catalyst layers[END_REF], and using (3.65) provides the following sequence of problems

O( -2 ) :          A 0 [λw 0 + f r (φ 0 )] = 0 in Y s , under no flux b.c. and with w 0 is Y s -periodic, A 0 φ 0 = 0 in Y s , with n sp • ∇φ 0 = 0 on ∂Y w s := ∂Y s ∩ ∂Y p and φ 0 Y s -periodic, (3.68)
O( -1 ) :

         A 0 [λw 1 + f r (φ 0 )φ 1 ] = -A 1 [λw 0 + f r (φ 0 )] in Y s , under no flux b.c. and w 1 Y s -periodic, A 0 φ 1 = -A 1 φ 0 in Y s , with n sp • ∇φ 1 = 0 on ∂Y s ∩ ∂Y p and φ 1 Y s -periodic, (3.69) O( 0 ) :                    A 0 λw 2 + 1 2 f r (φ 0 )φ 2 1 + f r (φ 0 )φ 2 = -(A 2 [λw 0 + f r (φ 0 )] + A 1 [λw 1 -f r (φ 0 )φ 1 ]) -∂ t A -1 2 w 0 in Y s , under no flux b.c. and with w 2 Y s -periodic, A 0 φ 2 = -A 2 φ 0 -A 1 φ 1 + w 0 in Y s , with n sp • ∇φ 2 = g on ∂Y s ∩ ∂Y p and φ 2 Y s -periodic. (3.70)
Step 2: (Analysis of the sequence of problems) The first problem (3.68) immediately implies that the leading-orders φ 0 and w 0 are independent of the microscale y [START_REF] Schmuck | First error bounds for the porous media approximation of the Poisson-Nernst-Planck equations[END_REF][START_REF] Mei | Homogenization Methods for Multiscale Mechanics[END_REF]. This suggests the following ansatz for w 1 and φ 1 , i.e., ∂x k ∂x i , where we applied the chain rule ∂ 2 fr(φ 0 ) ∂x k ∂x j = f r (φ 0 ) ∂φ 0 ∂x k ∂φ 0 ∂x j + f r (φ 0 ) ∂ 2 φ 0 ∂x k ∂x j and integration by parts. Adding now the term -A 2 f (φ 0 ) and using (3.72) finally gives -A 1 [f r (φ 0 )φ 1 ] -A 2 f r (φ 0 ) = div Ms ∇f r (φ 0 ) .

w 1 (x,y,t) = - d k=1 ξ k w (y) ∂w 0 ∂x k (x,t), φ 1 (x,y,t) = - d k=1 ξ k φ (y) ∂φ 0 ∂x k (x,t) = φ 1 (x,y,t). ( 3
ξ k φ : -d i,j=1 ∂ ∂y i δ ik -δ ij ∂ξ k φ ∂y j = -div e k -∇ y ξ k φ = 0 in Y s , (3.72
These considerations finally lead to the following effective equation for φ 0 , i.e., q ∂φ 0 ∂t = div Ms ∇f r (φ 0 ) + λ q div Ms ∇ div Ms ∇φ 0 -1 λ g0 .

(3.78)

In the case where the homogeneous free energy F r is the classical double-well potential F r (φ) = 1/4(φ 2 -1) 2 , one can verify the solvability of (3.78) along with the arguments in [START_REF] Novick-Cohen | On cahn-hilliard type equations[END_REF]. 2

Discussion and conclusions

We have given a short summary of widely used charge transport formulations in Li-ion battery systems. Based on well-accepted microscopic composite cathode equations which take into account ionic transport in the polymer electrolyte, electron transport as well as Li-intercalation undergoing a possible phase transformation in the solid phase, and Butler-Volmer reactions across the solidelectrolyte interface, we derived a novel upscaled/homogenized formulation. To the best of our knowledge, this seems to be the first two-scale asymptotic homogenization of composite cathodes for Li-batteries that also accounts for thermodynamically relevant phase transformations of crystalline host materials such as FePO 4 . For a literature overview on recent results available, we refer the interested reader to Tab. 1.

The promising feature of being able to account for three different scales, i.e., pore scale l p , Darcy's scale l D , and macroscale L, with the two-scale asymptotic expansion method by exploiting coupling of different physical processes, provides novel modelling opportunities not only in the context of batteries, but also in fuel cell systems, transport in porous media, and oil recovery for instance. We expect that the derived equations serve as a convenient computational model that allows to apply well-known numerical strategies available for homogeneous domains without ending up with high-dimensional discrete equations. Since the homogenization method is, compared to volume averaging, still less popular in applied sciences and engineering, we believe that the presented upscaling ideas and mathematical tools will become increasingly interesting in applications due to its reliable and systematic basis relying on rigorous convergence methods such as two-scale convergence [START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF] or Gamma convergence, e.g. [START_REF] Maso | An introduction to Γ-convergence[END_REF].

Figure 1 :

 1 Figure 1: Left: Reference cell Y consisting of a polymer solution Y p together with the solid intercalation compound Y s , i.e., Y = Y p ∪ Y s . The cell is electrically conducting due to the thin carbon black pathways connecting Y s . Right: Reference cell Y D , here for instance composed of an elementary pore cell Y such as depicted on the left-hand side. The cell Y D takes the influence of obstacles to momentum transport into account. The union of the Y p phases in Y D are then denoted by Y p D and similarly Y s D represents the union of the Y s phases in Y D .

Figure 2 :

 2 Figure 2: Left: A 2D periodic composite cathode D = D p ∪ D s of characteristic length L with solid intercalation particles D s surrounded by a polymer electrolyte solution. This black lines are electric conducting pathways in the form carbon black. This cell motivates a layer by layer assembly. We note that this figure can serve as a reference cell that represents a characteristic configuration of intercalation particles electrically connected in various ways via carbon black, and surrounded by polymer electrolyte. Right: A 3D sample reference cell of a possible cathode composed of FePO 4 particles connected via carbon black particles (thin lines). Obviously, more complex cells can be used.

( 1 )( 2 )

 12 Diffusive and electro-migrational transport in the polymer solution D p p := z∈Z d p (Y p + z) ∩ D ; Interstitial diffusion inside a solid intercalation compound D p s := z∈Z d p (Y s + z) ∩ D (described by a mean field phase field equation); and (3) Reaction kinetics on the interface I p := ∂D p p ∩ ∂D p s given by the Butler-Volmer equations (2.20) which are imposed as boundary conditions.

  ) where τ is either the diffusion time τ D := L 2 D , the interfacial transport time τ s := L 2 Ms , or the convection time τ u = L U ref . U ref is a reference velocity and c ref a reference salt concentration. We also use the identity k B T e = RT F and the electric potential ψ stands for the electric potentials of the solid phase ψ s and of the liquid phase ψ p . The parameter λ D := pkB T 2e 2 c ref 1/2

2 D.

 2 .36) where the parameters coupling the Butler-Volmer reactions to the different charge transport equations are defined by β + = ipsL eD , β s = ipsL eM , and β ψ = F RT ipsL σ ref for σ ref = εpD λ A summary of physical and dimensionless parameters are given in Tab. 2.

2 .Table 2 :

 22 Summary of available material parameters to model composite cathodes.ParameterValue range (a) Name L ≈ 100 -1000µm thickness of cathode ≈ 1µm radius of cathode particles/size of reference cellλ D < 1nm (d) Debye length D = D ± ≈ 1.5 -4.5 × 10 -10 m 2 /s (b) diffusion coeff. in D p p D s = k B T M s ≈ 5.0 × 10 -13 diffusion coeff. in D p s σ ≈1.0 × 10 4 S/m conductivity in D p s α a , α c 0.5 transfer c sm ≈ 29000mol/m 3 max. concentration in D p s c pm ≈ 3920mol/m 3 max. concentration in D p p k f b ≈ 10 -10 m 4 /(mol • s) reaction rate at I p i ps ≈ 12.6A/m 2 exchange current density |Λ| ≈ 5 -15 × 10 3 m 2 /kg (c) specific interfacial/surface area (BET) (a) Doyle et al. [14] & Garcia et al. [19], (b) Capiglia et al. [9], (c) Choy et al. [11],(d) Rademaker et al.[START_REF] Rademaker | On the deviation of electro-neutrality in li-ion battery electrolytes[END_REF]

  as in Section 2.2 where we introduced the microscopic formulation of the composite cathode based on a reference FePO 4 -particle configuration Y and a reference flow/filtration cell Y D .

  54) and (3.55) we can rewrite (3.56) in the more intuitive form,

( 3 .

 3 63) System (3.63) defines the reference cell problems for the upscaled Poisson-Nernst-Planck equations describing binary symmetric electrolytes in porous media. Thanks to the simplification of (3.60) via (3.62), the well-posedness of (3.63) follows by Lax-Milgram's theorem.

  .71) Inserting (3.71) into (3.69) 2 provides an equation for the correctors ξ k w and ξ k φ . The resulting equation for ξ k φ can be immediately written for 1 ≤ k ≤ d as,

a 1 χ ∂Y w 1 s+ a 2 χ ∂Y w 2 s∂ 2

 1122 ) for ξ k φ (y) Y -periodic with M Ys (ξ k φ ) = 0 and associated boundary condition n • ∇ξ k φ + e k = 0 on ∂Y w s := ∂Y s ∩ ∂Y p . To study (3.69) 1 , we first rewrite B 0 [f r (φ 0 )φ 1 ] and B 1 f r (φ 0 ) as follows A 0 [f r (φ 0 )φ 1 ] =the same for w 1 and w 0 and using (3.71) leads then to r (φ 0 ) ∂x k , (3.74) in Y s . Since the chemical potential µ(φ) = δE(φ) δφ is locally in thermodynamic equilibrium, we finally obtain the reference cell problem for ξ k w , 1 ≤ k ≤ d, for given ξ k φ w (y) is Y -periodic with boundary condition d i,j=1 n i δ ij ∂ξ k w ∂y j -δ ik + δ ik -∂Y w s = ∂Y s ∩ ∂Y p and M Ys (ξ k w ) = 0. Step 3: (Second order problem: upscaled equation) The last problem (3.70) then gives the upscaled equations by a solvability constraint, i.e., the Fredholm alternative. Hence, the solvability of (3.70) 2 is achieved by setting -0 ∂x i ∂x k = |Y s |w 0 + g0 . (3.76) where g0 := -γ C h ∂Y w s do(y) and ∂Y w s = ∂Y w 1 s ∪ ∂Y w 2 s . (3.76) can be written compactly as -∆ Ms φ 0 := -div Ms ∇φ 0 = qw 0 + g0 by (3.40) 4 . Applying the same ideas to equation (3.70) 1 leads to Ys-λ(A 2 w 0 + A 1 w 1 ) -1 λ A 1 [f r (φ 0 )φ 1 ] -1 λ A 2 f r (φ 0 ) -∂ t A -1 2 w 0 dy = 0,(3.77)where the first two terms can be rewritten with (3.75) by Ys -(A 2 w 0 + A 1 w 1 ) dy = div Ms ∇w 0 .The third integrand in (3.77) becomes -A 1 [f r (φ 0 )φ 1 ] = -fr(φ 0 )

  .59) Equation (3.59) 3 is classical in elliptic homogenization theory[START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Chechkin | Homogenization: Methods and Applications[END_REF][START_REF] Cioranescu | An Introduction to Homogenization[END_REF][START_REF] Hornung | Homogenization and Porous Media[END_REF]. Hence, we immediately can conclude that Ψ(t,x) is independent of the microscale y ∈ Y . Entering with this result into (3.59) 1 and (3.59) 2 further implies that also C ± (t,x) and Ψ(t,x) are independent of the microscale y ∈ Y p .The next problem in the sequence, i.e., (3.57), has the following explicit form,

Acknowledgements

We acknowledge financial support from EPSRC Grant No EP/P011713/1.