Nicolas Bonichon
email: bonichon@labri.fr.

Mireille Bousquet-Mélou

AND Paul Dorbec
email: dorbec@labri.fr.

Claire Pennarun
email: cpennaru@labri.fr.

ON THE NUMBER OF PLANAR EULERIAN ORIENTATIONS

Keywords: planar maps, Eulerian orientations, algebraic generating functions

The number of planar Eulerian maps with n edges is well-known to have a simple expression. But what is the number of planar Eulerian orientations with n edges? This problem appears to be difficult. To approach it, we define and count families of subsets and supersets of planar Eulerian orientations, indexed by an integer k, that converge to the set of all planar Eulerian orientations as k increases. The generating functions of our subsets can be characterized by systems of polynomial equations, and are thus algebraic. The generating functions of our supersets are characterized by polynomial systems involving divided differences, as often occurs in map enumeration. We prove that these series are algebraic as well. We obtain in this way lower and upper bounds on the growth rate of planar Eulerian orientations, which appears to be around 12.5.

Introduction

The enumeration of planar maps (graphs embedded on the sphere) has received a lot of attention since the sixties. Many remarkable counting results have been discovered, which were often illuminated later by beautiful bijective constructions. For instance, it has been known 1 since 1963 that the number of rooted planar Eulerian maps (i.e., planar maps in which every vertex has even degree) with n edges is [START_REF] Tutte | A census of planar maps[END_REF]:

m n = 3 • 2 n-1 (n + 1)(n + 2) 2n n . (1)
A bijective explanation involving plane trees can be found in [START_REF] Bousquet-Mélou | Enumeration of planar constellations[END_REF]. The associated generating function M (t) = n≥0 m n t n is known to be algebraic, that is, to satisfy a polynomial equation. More precisely: t 2 + 11t -1 -(8t 2 + 12t -1)M (t) + 16t 2 M (t) 2 = 0. Beyond their enumerative implications, bijections involving maps have been applied to encode, sample and draw maps efficiently [START_REF] Bonichon | Canonical decomposition of outerplanar maps and application to enumeration, coding and generation[END_REF][START_REF] Castelli Aleardi | Succinct representations of planar maps[END_REF][START_REF] Fusy | Transversal structures on triangulations: a combinatorial study and straight-line drawings[END_REF][START_REF] Poulalhon | Optimal coding and sampling of triangulations[END_REF]. More recently, they have played a key role in the study of large random planar maps, culminating with the existence of a universal scaling limit known as the Brownian map [START_REF] Gall | Random geometry on the sphere[END_REF].

Planar maps equipped with an additional structure (e.g. a spanning tree [START_REF] Mullin | On the enumeration of tree-rooted maps[END_REF], a proper colouring [START_REF] Tutte | Chromatic sums for rooted planar triangulations: the cases λ = 1 and λ = 2[END_REF][START_REF] Tutte | Map-colourings and differential equations[END_REF], an Ising or Potts configuration [START_REF] Baxter | Dichromatic polynomials and Potts models summed over rooted maps[END_REF][START_REF] Borot | Loop models on random maps via nested loops: case of domain symmetry breaking and application to the Potts model[END_REF][START_REF] Boulatov | The Ising model on a random planar lattice: the structure of the phase transition and the exact critical exponents[END_REF][START_REF] Bousquet-Mélou | The degree distribution of bipartite planar maps: applications to the Ising model[END_REF][START_REF] Bouttier | Blocked edges on Eulerian maps and mobiles: application to spanning trees, hard particles and the Ising model[END_REF][START_REF] Daul | q-States Potts model on a random planar lattice[END_REF][START_REF] Eynard | The Potts-q random matrix model: loop equations, critical exponents, and rational case[END_REF][START_REF] Kazakov | Ising model on a dynamical planar random lattice: exact solution[END_REF]...) are also much studied, both in combinatorics and in theoretical physics, where maps are considered as a model for two-dimensional quantum gravity [START_REF] Di Francesco | 2D gravity and random matrices[END_REF]. However, for many of these structures, we are still in the early days of the study, as even their enumeration remains elusive, not to mention bijections and asymptotic properties.

Recent progresses in this direction include the enumeration of planar maps weighted by their Tutte polynomial, or equivalently, maps equipped with a Potts configuration. The associated generating function P (t) is known to be differentially algebraic. That is, there exists a polynomial equation relating P (t) and its derivatives [START_REF] Bernardi | Counting colored planar maps: algebraicity results[END_REF][START_REF] Bernardi | Counting coloured planar maps: differential equations[END_REF]. The Tutte polynomial has many interesting specializations (in particular, it counts all structures cited above, like spanning trees and colourings) and several special cases had been solved earlier. One key tool in the solution is that the Tutte polynomial of a map can be computed inductively, by deleting and contracting edges.

Another solved example, which does not seem to belong to the Tutte/Potts realm, consists of maps (in fact, triangulations) equipped with certain orientations called Schnyder orientations. The results obtained there have analogies with those obtained for another class of orientations, called bipolar (which do belong to the Tutte realm). Indeed, for both classes of oriented maps:

• oriented maps are counted by simple numbers, which are also known to count other combinatorial objects (various lattice paths and permutations, among others); • there exist nice bijections explaining these equi-enumeration results [START_REF] Bonichon | A bijection between realizers of maximal plane graphs and pairs of non-crossing Dyck paths[END_REF][START_REF] Bonichon | Baxter permutations and plane bipolar orientations[END_REF][START_REF] Felsner | Bijections for Baxter families and related objects[END_REF][START_REF] Fusy | Bijective counting of plane bipolar orientations and Schnyder woods[END_REF];

• for a fixed map M , the set of Schnyder/bipolar orientations of M has a lattice structure [START_REF] Propp | Lattice structure for orientations of graphs[END_REF][START_REF] Felsner | Lattice structures from planar graphs[END_REF][START_REF] Ossona De Mendez | Orientations bipolaires[END_REF]. The above bijections, once specialized to maps equipped with their (unique) minimal orientation, coincide with attractive bijections designed earlier for (unoriented) maps [START_REF] Bernardi | Intervals in Catalan lattices and realizers of triangulations[END_REF][START_REF] Bonichon | Baxter permutations and plane bipolar orientations[END_REF]; • specializing the bijections further to maps that have only one Schnyder/bipolar orientation also yields interesting combinatorial results [START_REF] Bernardi | Intervals in Catalan lattices and realizers of triangulations[END_REF][START_REF] Bonichon | Baxter permutations and plane bipolar orientations[END_REF].

These observations led us to wonder about another natural class of orientations, namely those in which every vertex as equal in-and out-degree, known as Eulerian orientations (Figure 1). Clearly, a map needs to be Eulerian to admit an Eulerian orientation. The condition is in fact sufficient (such maps even admit an Eulerian circuit [START_REF] Hierholzer | Ueber die Möglichkeit, einen Linienzug ohne Wiederholung und ohne Unterbrechung zu umfahren[END_REF]). One analogy with the above two classes is that the set of Eulerian orientations of a given planar map can be equipped with a lattice structure [START_REF] Propp | Lattice structure for orientations of graphs[END_REF][START_REF] Felsner | Lattice structures from planar graphs[END_REF]. Moreover, Eulerian maps (equivalently, Eulerian maps equipped with their minimal Eulerian orientation) have rich combinatorial properties: not only are they counted by simple numbers (see [START_REF] Alm | Random self-avoiding walks on one-dimensional lattices[END_REF]), but they are equinumerous with several other families of objects, like certain trees [START_REF] Bousquet-Mélou | Enumeration of planar constellations[END_REF] and permutations [START_REF] Bóna | Exact enumeration of 1342-avoiding permutations: a close link with labeled trees and planar maps[END_REF][START_REF] Fusy | Bijective counting of involutive Baxter permutations[END_REF]. And they are often related to them by beautiful bijections.

Hence our plan to count Eulerian orientations with n edges. However, this appears to be a difficult problem. In fact, we even lack a way to compute the corresponding numbers in, say, polynomial time. This leads us to resort to approximation methods that are ubiquitous when studying hard counting problems, like the enumeration of self-avoiding walks [START_REF] Alm | Random self-avoiding walks on one-dimensional lattices[END_REF][START_REF] Fisher | Excluded-volume problem and the Ising model of ferromagnetism[END_REF][START_REF] Guttmann | Effect of confinement: polygons in strips, slabs and rectangles[END_REF][START_REF] Pönitz | Improved upper bounds for self-avoiding walks in Z d[END_REF], or polyominoes [START_REF] Barequet | Counting polyominoes on twisted cylinders[END_REF][START_REF] Jensen | Counting polyominoes: A parallel implementation for cluster computing[END_REF][START_REF] Klarner | A procedure for improving the upper bound for the number of n-ominoes[END_REF]: denoting by O the set of Eulerian orientations, we construct subsets and supersets of O, indexed by an integer parameter k, which converge to O as k increases. And we count the elements of these sets.

One difference between our study and those dealing with tricky objects on regular lattices (like the above mentioned self-avoiding walks and polyominoes) is worth noting. The subsets and supersets that are defined to approximate lattice objects often have a one-dimensional structure, and rational generating functions that can be obtained using a transfer matrix approach. A typical example is provided by self-avoiding walks confined to a strip of fixed width. But our subsets and supersets of orientations belong to the world of maps (or random lattices in the physics terminology), and have algebraic generating functions. More precisely, our subsets have a branching, tree-like structure, which yields a system of algebraic equations for their generating functions, and a universal asymptotic behaviour in λ n n -3/2 (for a growth rate λ depending on the index k). The generating functions of our supersets are more mysterious. They are bivariate series given by systems of equations involving divided differences of the form F (t; x) -F (t; 1)

x -1 , and we have to resort to a deep theorem in algebra, due to Popescu [START_REF] Popescu | General Néron desingularization and approximation[END_REF], to prove their algebraicity for all k (we also solve these systems for small values of k). We conjecture that their asymptotic behaviour is also universal, this time in λ n n -5/2 , as for planar maps (again, for varying λ).

Here is now an outline of the paper. In Section 2 we first present a simple recursive decomposition of (rooted) Eulerian orientations, based on the contraction of the root edge, and then a variant of this decomposition. Thanks to this variant, we can compute the number o n of Eulerian orientations having n edges for n ≤ 15 (Figure 2). By attaching two orientations at their root vertex, we see that the sequence (o n ) n≥0 is super-multiplicative:

o m+n ≥ o m o n .
This classically implies that the limit µ of o (

) 2
(see Fekete's Lemma in [59, p. 103]). We call µ the growth rate of Eulerian orientations. It is bounded from below by the growth rate 8 of Eulerian maps, and from above by the growth rate 16 of Eulerian maps equipped with an arbitrary orientation. Our data for n ≤ 15 suggest than µ is around 12.5 (Figure 2, right). Using differential approximants [35], Tony Guttmann predicts µ = 12.568 . . ., and an asymptotic behaviour o n ∼ cµ n n -γ with γ = 2. [START_REF]The on-line encyclopedia of integer sequences[END_REF]). Right: A plot of o n+1 /o n vs. 1/n, for n = 4, . . . , 14, suggests that the growth rate of Eulerian orientations, located at the intercept of the curve and the y-axis, is around 12.5.

Sections 3 and 4 deal with two families of subsets of Eulerian orientations. The first family uses our first recursive decomposition of Eulerian orientations, and should be considered as a warm up. The second family uses the variant of the standard decomposition of orientations. Its study is a bit more involved, but it gives better bounds on the growth constant. Both families are proved to have algebraic generating functions and a tree-like asymptotic behaviour in λ n n -3/2 . The next two sections deal with two families of supersets of Eulerian orientations. Both are proved to have algebraic generating functions, and we conjecture a map-like asymptotic behaviour in λ n n -5/2 . We solve our systems of equations explicitly for small values of k, and thus obtain Table 1. All calculations are supported by Maple sessions available on our web pages. We gather in Section 7 a few final comments and questions.

Table 1. Growth rates and cardinalities of subsets (L (k) and L (k) ) and supersets (U (k) and U (k) ) of Eulerian orientations. The table also records the degrees of the associated generating functions, which are systematically algebraic. The symbol refers to a numerical estimate. The other growth rates are algebraic numbers known exactly via their minimal polynomial. Let us mention that counting Eulerian orientations of 4-valent (rather than Eulerian) maps might be simpler: in this case, the number of Eulerian orientations is a specialization of the Tutte polynomial [START_REF] Welsh | The Tutte polynomial[END_REF], and in fact some results exist in the physics literature [START_REF] Kostov | Exact solution of the six-vertex model on a random lattice[END_REF][START_REF] Zinn-Justin | The six-vertex model on random lattices[END_REF]. In the final section, we discuss further this problem, which seems to deserve more attention.

Recursive decompositions of Eulerian orientations

In this section, after a few basic definitions, we recall the standard recursive decomposition of Eulerian maps based on the contraction of the root edge, which can be traced back to the early papers of Tutte (e.g. [START_REF] Tutte | On the enumeration of planar maps[END_REF]). We then adapt it to decompose Eulerian orientations. We also introduce a variant of the standard decomposition of Eulerian maps, based on a notion of prime maps, and adapt it again to Eulerian orientations. This variant is slightly more involved, but turns out to be more effective: it allows us to compute the numbers o n for larger values of n, and it also leads to better lower and upper bounds on these numbers (Sections 4 and 6).

Definitions

A planar map is a proper embedding of a connected planar graph in the oriented sphere, considered up to orientation preserving homeomorphism. Loops and multiple edges are allowed (Figure 1). The faces of a map are the connected components of its complement. The number of edges of a planar map M is denoted by e(M ). The degree of a vertex is the number of edges incident to it, counted with multiplicity (e.g., a loop counts twice). A corner is a sector delimited by two consecutive edges around a vertex; hence a vertex of degree k defines k corners. Alternatively, a corner can be described as an incidence between a vertex and a face.

For counting purposes it is convenient to consider rooted maps. A map is rooted by choosing a corner, called the root corner. The vertex and face that are incident at this corner are respectively the root vertex and the root face. The root edge is the edge that follows the root corner in counterclockwise order around the root vertex. In figures, we indicate the rooting by an arrow pointing to the root corner, and take the root face as the infinite face (Figure 1).

From now on, every map is planar and rooted, and these precisions will often be omitted. By convention, we include among rooted planar maps the atomic map having one vertex and no edge.

A map M is Eulerian if every vertex has even degree. In this case, we denote by dv(M ) the half degree of the root vertex. An Eulerian orientation is a map with oriented edges, in which the in-and out-degrees of every vertex are equal. Note that the underlying map must be Eulerian. We denote by M the set of Eulerian maps, and by O the set of Eulerian orientations.

Eulerian maps: standard decomposition

Consider an Eulerian map M , not reduced to the atomic map, and its root edge e. If e is a loop, then M is obtained from two smaller maps M 1 and M 2 by joining M 1 and M 2 at their root vertices and adding a loop surrounding M 1 (Figure 3, left). The maps M 1 and M 2 are themselves Eulerian (because the sum of vertex degrees in a map is even, so that one cannot have a single odd vertex in M 1 or M 2 ). We call this operation the merge of M 1 and M 2 .

If the root edge e is not a loop, then we contract it, which gives a smaller Eulerian map M . Note however that several maps give M after contracting their root edge. All such maps can be obtained from M as follows (see Figure 3, right): we split the root vertex v of M into two vertices v and v joined by an edge (which will be the root edge), and distribute the edges adjacent to v between v and v in such a way the degrees of v and v remain even. This operation is called a split of M , and more precisely an i-split if v has degree 2i in the larger map. Note that if v has degree 2d in M , then i can be chosen arbitrarily between 1 and d.

Let M (t; x) be the generating function of Eulerian maps, counted by edges (variable t) and by the half degree of the root vertex (variable x):

M (t; x) = M ∈M t e(M ) x dv(M ) = d≥0 x d M d (t),
where M d (t) denotes the edge generating function of Eulerian maps with root vertex degree 2d. The above construction translates into the following functional equation, which we explain below:

M (t; x) = 1 + txM (t; x) 2 + t d≥0 M d (t)(x + x 2 + • • • + x d ) = 1 + txM (t; x) 2 + t d≥0 M d (t) x d+1 -x x -1 = 1 + txM (t; x) 2 + tx x -1 (M (t; x) -M (t; 1)). (3)
On the first line, the term 1 accounts for the atomic map, the next term for maps obtained by merging two smaller maps, and the third term for maps obtained from a split.

Eulerian orientations: standard decomposition

Our recursive decomposition for Eulerian orientations is essentially the same as for Eulerian maps: if the root edge is a loop, we delete it and obtain two orientations, which are both Eulerian (in any oriented map, the sum over all vertices of in-degrees equals the sum of out-degrees, hence one cannot have a single vertex with distinct in-and out-degrees); otherwise we contract the root edge, which gives a smaller Eulerian orientation.

However, care must be taken when going in the opposite direction, that is, when constructing large orientations from smaller ones. The first type of orientations, obtained by a merge, do not raise any difficulty; one can orient the new root edge (the loop) in two different ways (Figure 4, left). But consider now an Eulerian orientation O , with root vertex v of degree at least 2i, and perform an i-split on O : is there a way to orient the new edge so as to obtain an orientation O that is still Eulerian? The answer is yes if and only if the numbers of in-and out-edges in the 4, right). In this case, we say that the i-split, performed on O , is legal. Note that the 1-split and the d-split are always legal, where 2d is the degree of the root vertex of O . The fact that not all splits are legal makes it difficult to write a single functional equation for the generating function of Eulerian orientations. However, we can write an infinite system of equations relating the generating functions of orientations with prescribed orientations at the root.

} 2i -1 v v v M 1 M 2 M M M Figure 3.
Let us be more precise. Given an Eulerian orientation O with root vertex v of degree 2d, the root word w(O) of O is a word of length 2d on the alphabet {0, 1} describing the orientation of the edges around v (in counterclockwise order, starting from the root corner): the k-th letter of w(O) is 0 (resp. 1) if the k-th edge around v is in-going (resp. out-going). Note that this word is always balanced, meaning that it contains as many 0's as 1's. We call a word w quasi-balanced if the number of 0's and 1's in w differ by ±1. The length (number of letters) of w is denoted by |w|, while the number of occurrences of the letter a in it is denoted by |w| a . We define the balance of w to be b(w) := ||w| 1 -|w| 0 |. The empty word is denoted by ε.

Now we can decide from the root word of O if the i-split of O is legal: this holds if and only if the last 2i -1 letters of w(O ) form a quasi-balanced word.

For w a word on {0, 1}, let O w (t) ≡ O w be the generating function of Eulerian orientations having w as root word, counted by their edge number. Clearly, O w = 0 if w is not balanced and O ε = 1 (accounting for the atomic map). Now if w is non-empty and balanced,

O w = t auāv=w O u O v + t u O uws . ( 4
)
This identity is illustrated in Figure 4. Here, a stands for any of the letters 0, 1, and the first sum runs over all factorisations of w of the form auāv, with ā := 1 -a. This sum counts orientations obtained by a merge. The second sum runs over all possible words u, and w s denotes the suffix of w of length |w| -1. This sum counts orientations obtained by a (legal) split of an orientation having root word uw s . Now the generating function O of Eulerian orientations is

w O w ,
where the sum runs over all (balanced) words w.

We do not know how to solve this system. But a map with n edges has a root word of length at most 2n, and hence we can use our system to compute the numbers o n for n small. We obtain in this way the first 11 values of Figure 2. In Sections 3 to 6, we define subsets and supersets of O that we can generate by just looking at the last 2k -1 letters of the root word (for k fixed). This allows us to write finitely many equations for the generating functions of these subsets and supersets. Solving them gives lower and upper bounds on the growth rate of Eulerian orientations. However, we obtain more precise bounds by using a variant of the standard decomposition of maps and orientations. We now present this variant.

} 2i -1 v v v O 1 O 2 O O O w(O 1 ) = u w(O 2 ) = v w(O) = auāv w(O ) = uw s , |w s | = 2i -1 w(O) = aw s , balanced

Prime decomposition of maps and orientations

A (non-atomic) map is said to be prime if the root vertex appears only once when walking around the root face. A planar map M can be seen as a sequence of prime maps M 1 , . . . , M (Figure 5). We say that the M i are the prime submaps of M , and denote Now take a prime Eulerian map M , and apply the standard decomposition of Section 2.2, illustrated in Figure 3: either M is an (arbitrary) Eulerian map M 1 surrounded by a loop, or M is obtained by an i-split in another Eulerian map M , provided the last prime submap of M (in counterclockwise order) has root degree at least 2i (otherwise, the resulting map would not be prime). Alternatively, if M = M 1 • • • M , we can obtain M by performing an i-split in the prime map M , and attaching the map M := M 1 • • • M -1 at the new vertex v created by this split (Figure 6).

M = M 1 • • • M . Note that if M is Eulerian, then each M i is Eulerian too. M 1 M 2 M 3 M M = M 1 M 2 M 3
} 2i -1 v M (prime) M + M 1 M v v M Figure 6
. Construction of a prime Eulerian map: add a loop around any Eulerian map, or split a prime Eulerian map, and attach an arbitrary Eulerian map at the end of its root edge. This alternative decomposition of Eulerian maps gives a system of two equations defining the generating function M (t; x) of Eulerian maps (still counted by edges and root vertex degree) and its counterpart M (t; x) for prime maps:

M (t; x) = 1 + M (t; x)M (t; x), M (t; x) = txM (t; x) + txM (t; 1) M (t; x) -M (t; 1)
x -1 .

In the first equation, the term M (t; x) accounts for the last prime submap attached at the root vertex (denoted M above). In the second equation, the divided difference (M (t; x) -M (t; 1))/(x -1) has the same explanation as in (3). This equation is easily recovered by eliminating M (t; x) from the above system.

This decomposition can also be applied to Eulerian orientations: an Eulerian orientation is a sequence of prime Eulerian orientations, and a prime orientation is either obtained by adding an oriented loop around another orientation, or by performing a legal split in a prime orientation, and attaching another orientation at the vertex v created by the split.

Thus, denoting again by O w the generating function of orientations with root word w, and by O w its counterpart for prime orientations, we now have

O w = O w = 0 if w is not balanced, O ε = 1, O ε =
0 and finally for w balanced and non-empty,

O w = uv=w O u O v , O w = tO wc + tO u O uws .
In the second equation, w c denotes the central factor or w of length |w| -2, and O = w O w is the generating function of all Eulerian orientations. Recall that w s is the suffix of w of length |w| -1.

Using these equations, we have been able to push further the enumeration of Eulerian orientations of small size, thus obtaining the values of Figure 2.

Subsets of Eulerian orientations, via the standard decomposition

In this section and the three following ones, we define certain subsets and supersets of Eulerian orientations, indexed by an integer k, which converge (monotonously) to the set O of all Eulerian orientations as k tends to infinity. Those sets are respectively denoted by L (k) (and L (k) ), U (k) (and U (k) ), as they give lower and upper bounds on the numbers o n and their growth rate. For each set, we give a system of functional equations defining its generating function: for the subsets L (k) and L (k) , these systems are algebraic, so that the associated generating functions are algebraic series. For the supersets U (k) and U (k) , the systems define bivariate series and involve divided differences as in [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF]. However, we prove that the resulting series are also algebraic.

Recall from Figure 4 that planar Eulerian orientations can be obtained recursively from the atomic map by either:

• the merge of two orientations O 1 , O 2 ∈ O (with the root loop oriented in either way),

• or a legal split on an orientation O ∈ O.

Definition 1. Let k ≥ 1. Let L (k)
be the set of planar orientations obtained recursively from the atomic map by either:

• the merge of two orientations O 1 , O 2 ∈ L (k) (with the root loop oriented in either way),

• or a legal i-split on an orientation O ∈ L (k) such that i ≤ k or i = dv(O ). In other words, the only allowed splits are the small splits (i ≤ k) and the maximal split

(i = dv(O )).
Obviously, all orientations of L (k) are Eulerian. Moreover, the sets L (k) form an increasing sequence since more and more (legal) splits are allowed as k grows. Finally, all Eulerian orientations of size n belong to L (n) (and even to L (n-2) ). Hence the limit of the sets L (k) is the set O of all Eulerian orientations.

Figure 7 shows a (random) orientation of L (1) .

Figure 7. An Eulerian orientation in L (1) , taken uniformly at random among those with 20 edges.

3.1. An algebraic system for L (k) In this section, k is a fixed integer.

Definition 2. A word w on {0, 1} is valid (for k) if there exists a balanced word of length 2k having w as a factor. Equivalently, the balance of w satisfies b(w) ≤ 2k -|w|. This holds automatically if |w| ≤ k.

Given a word w, it will be convenient to have notation for several words that differ from w by one or two letters. We have already defined w c , the central factor of w of length |w| -2, and w s , the suffix of w of length |w| -1. We similarly define w p as the prefix of w of length |w| -1.

Finally, if w is quasi-balanced, then ←w stands for the unique balanced word of the form aw, for a ∈ {0, 1}.

For any word w, we denote by L (k) w (t) the generating function of orientations of L (k) whose root word ends with w, counted by edges. In particular, the generating function counting all orientations of

L (k) is L (k) ε (t). We denote by K (k)
w (t) the generating function of orientations of L (k) having root word exactly w. In order to lighten notation, we often omit the dependence of our series in t and the superscript (k).

We now give equations defining the series L w (for |w| ≤ 2k -1) and K w (for |w| ≤ 2k). First, we note that K w = 0 if w is not balanced, and that K ε = 1. Now for w balanced of length between 2 and 2k, we have:

K w = t w=auāv K u K v + tL ws , (5)
where, as before, a is any of the letters 0, 1. This equation is analogous to (4): the first term counts orientations obtained from a merge, the second orientations obtained from a split. Now for L w , with w of length at most 2k -1, we have:

L w = 1 w=ε + 2tL ε L w + t w=uav L u K v + t w=auāv K u K v + t (L w -1 w=ε ) + t u=vw 2≤|u|≤2k u balanced (L us -K u ). (6)
This equation deserves some explanations. The first line counts the atomic map (if w = ε), and the orientations obtained by a merge. The second (resp. third, fourth) term of this line counts orientations such that no (resp. one, both) half-edge(s) of the root loop is/are involved in the suffix w of the root word. Equivalently, denoting by O 1 and O 2 the merged orientations, those three terms respectively correspond to

|w| ≤ 2 dv(O 2 ), 2 dv(O 2 ) < |w| < 2 + 2 dv(O 1 ) + 2 dv(O 2 ) and |w| = 2 + 2 dv(O 1 ) + 2 dv(O 2 ).
The second line counts orientations O obtained by a legal i-split in a smaller orientation O . The first term accounts for maximal splits (i = dv(O )), which, we recall, do not change the root word (note also that no split is possible on the atomic map). The second term counts orientations O obtained from a non-maximal split. The word u stands for the root word of O. The subtraction of K u comes from the condition that the split is not maximal. Proposition 3. Consider the collection of equations consisting of:

• Equation (5), written for all balanced words w of length between 2 and 2k, • Equation (6), written for all valid words w of length at most 2k -1.

In this collection, replace all trivial K-series by their value: K w = 0 when w is not balanced, K ε = 1. Let S 0 denote the resulting system. The number of series it involves is

f (k) = 2k + 2 k + 1 -1 + k-1 i=1 2i i . ( 7
)
The system S 0 defines uniquely these f (k) series. Its size can be (roughly) divided by two upon noticing that replacing all 0's by 1's, and vice-versa, in a word w, does not change the series L w nor K w .

Proof. To see that S 0 defines all the series that it involves, it suffices to note the factor t in the right-hand sides of ( 5) and ( 6), and to check that each series occurring in the right-hand side of some equation also occurs as the left-hand side of another. This is readily done, as any factor of a valid word is still valid.

Let us now count the equations of the system. The number of non-empty balanced words of length at most 2k is

k i=1 2i i .
Then, all words of length at most k are valid, while the number of valid words of length k + i,

for 1 ≤ i ≤ k -1, is k j=i k + i j .
(One can interpret j as the number of occurrences of 0 in the word.) Hence the number of equations in the system is

f (k) = k i=1 2i i + k i=0 2 i + k-1 i=1 k j=i k + i j .
The second sum evaluates to 2 k+1 -1. The third one is

k j=1 min(j,k-1) i=1 k + i j = k-1 j=1 j i=1 k + i j + k-1 i=1 k + i k = 1 + 2k + 1 k -2 k+1 + k k + 1 2k k -1 = 3k + 1 k + 1 2k k -2 k+1 .
The sums are evaluated using classical summation identities, or Gosper's algorithm [START_REF] Petkovšek | [END_REF]. The expression of f (k) given in the proposition then follows after elementary manipulations.

Remark 4. If w is such that 0w and 1w are both valid of length less than 2k, we can define L w by a simpler "forward" equation, without increasing the size of the system:

L w = K w + L 0w + L 1w . (8)
This is obviously smaller than (6), and possibly better suited to feed a computer algebra system. However, mixing equations of type ( 6) and ( 8) makes some proofs of Section 3.3 heavier.

Examples

3.2.1.

When k = 1, the system S 0 contains f (1) = 5 equations and reads

               K 01 = tK ε K ε + tL 1 , K 10 = tK ε K ε + tL 0 , L ε = 1 + 2tL ε L ε + t(L ε -1) + t(L 0 -K 10 + L 1 -K 01 ), L 0 = 2tL ε L 0 + tL ε K ε + tL 0 + t(L 0 -K 10 ), L 1 = 2tL ε L 1 + tL ε K ε + tL 1 + t(L 1 -K 01 ), (9)
with K ε = 1. Using the 0/1 symmetry, this system can be compacted into

     K 01 = t + tL 0 , L ε = 1 + 2tL ε L ε + t(L ε -1) + 2t(L 0 -K 01 ), L 0 = 2tL ε L 0 + tL ε + tL 0 + t(L 0 -K 01 ).
The variant mentioned in Remark 4 consists in replacing the second equation by

L ε = 1 + 2L 0 .
The reader may check that this is consistent with the above system.

Eliminating L 0 and K 01 gives a quadratic equation for the generating function

L ε = L (1) ε of Eulerian orientations in L (1) : 2tL 2 ε -L ε (1 -t) 2 -t 2 -2t + 1 = 0. (10)
We defer to Section 3.3 the study of the asymptotic behaviour of its coefficients.

3.2.2. When k = 2, the system S 0 contains f (2) = 21 equations, or 11 if we exploit the 0/1 symmetry:

                                                     K 10 = K 01 = t + tL 0 , K 1100 = tK 10 + tL 100 , K 1010 = t(K 10 + K 01 ) + tL 010 , K 0110 = tK 10 + tL 110 , L ε = 1 + 2tL ε L ε + t(L ε -1) + 2t(L 0 -K 10 + L 100 -K 1100 +L 010 -K 1010 + L 110 -K 0110 ), L 0 = L 1 = 2tL ε L 0 + tL ε + tL 0 + t(L 0 -K 10 + L 100 -K 1100 +L 010 -K 1010 + L 110 -K 0110 ), L 00 = L 11 = 2tL ε L 00 + tL 0 + tL 00 + t(L 100 -K 1100 ), L 10 = L 01 = 2tL ε L 10 + tL 1 + t + tL 10 + t(L 0 -K 10 + L 010 -K 1010 + L 110 -K 0110 ), L 100 = 2tL ε L 100 + tL 10 + tL 100 + t(L 100 -K 1100 ), L 010 = 2tL ε L 010 + t(L 01 + L ε K 10 ) + tL 010 + t(L 010 -K 1010 ), L 110 = 2tL ε L 110 + t(L 11 + L ε K 10 ) + tL 110 + t(L 110 -K 0110 ). (11)
The variant mentioned in Remark 4 consists in replacing the equations defining L ε , L 0 and L 10 by L ε = 1 + 2L 0 , L 0 = L 00 + L 10 and L 10 = K 10 + L 010 + L 110 respectively.

Eliminating all series but L ε gives a quartic equation for the generating function

L ε = L (2) ε
of Eulerian orientations in L (2) :

8t 3 L 4 ε -4t 2 (3t 3 + 4t 2 -6t + 3)L 3 ε + 2t(3t 5 -12t 4 -10t 3 + 14t 2 -10t + 3)L 2 ε + (t -1)(11t 5 -10t 4 -6t 3 -3t 2 -t + 1)L ε + (t -1)(5t 5 -4t 4 + 6t 3 -7t 2 + 5t -1) = 0. ( 12
)
We defer to Section 3.3 the study of the asymptotic behaviour of its coefficients.

Asymptotic analysis for subsets of Eulerian orientations

Here, we apply the theory of positive irreducible polynomial systems [30, Sec. VII.6] to prove the following asymptotic result. Proposition 5. For k ≥ 1, let ρ k denote the radius of convergence of the series L (k) ε , which counts orientations of L (k) . Then ρ k is the only singularity of L (k) ε of minimal modulus, and it is of the square root type: as t tends to ρ k from below, 1)) for non-zero constants α and β depending on k.

L (k) ε (t) = α -β 1 -t/ρ k (1 + o(
The number

(k)
n of orientations of size n in L (k) satisfies, as n tends to infinity:

(k) n ∼ cλ n k n -3/2 , where λ k = 1/ρ k and c = -β/Γ(-1/2).
Proof. We use the terminology of [START_REF] Flajolet | Analytic combinatorics[END_REF]Sec. VII.6.3]. Our first objective is to transform the system S 0 of Proposition 3 into a positive one. The obstructions to positivity come from the expression (6) of L w , and more precisely from the terms L ε -1 (when w = ε) and L us -K u , where u is balanced. These terms can be written L w -K← - w , where w = u s is quasi-balanced and ←w is the unique balanced word of the form aw, for a ∈ {0, 1}.

This leads us to define, for w quasi-balanced of length less than 2k, the series L + w := L w -K← - w . We will also need to define, for w balanced, L + w := L w -K w . These series have natural combinatorial interpretations in terms of orientations whose root word ends strictly with w (if w is balanced) or ←w (if w is quasi-balanced). Then we alter the original system S 0 as follows.

(i) For w balanced or quasi-balanced, we replace the equation ( 6) defining L w by an equation defining L + w :

L + w = 2tL ε L w + t w=uav (L u -K u )K v + tL + w + t u=vw,u =w |u|≤2k-1 u quasi-balanced L + u . (13)
To obtain it, either we get back to the explanation of ( 6) and remove from its right-hand side the terms that count orientations with root word exactly w (if w is balanced) or ←w (if w is quasi-balanced). Or we simply subtract from ( 6) Equation ( 5), written for w if w is balanced, for ←w if w is quasi-balanced. (ii) In the new system thus obtained, we replace every series K w such that w is not balanced by 0, every series L w such that w is balanced by K w +L + w , and every series L w such that w is quasi-balanced by K← - w + L + w . In particular, the series L u -K u occurring in (13) becomes L + u when u is balanced, L u otherwise. The only series L w that remain in the system are such that the balance of w is at least 2.

We thus obtain a positive system, denoted S 1 , defining the following series:

• K w , for w balanced of length between 2 and 2k, • L + w , for w balanced or quasi-balanced of length less than 2k, • L w , for w valid of length less than 2k and balance at least 2.

For instance, when k = 1, the system (9) becomes (after exploiting the 0/1 symmetry):

     K 01 = t + t(K 01 + L + 0 ), L + ε = 2t(1 + L + ε ) 2 + tL + ε + 2tL + 0 , L + 0 = 2t(1 + L + ε )(K 01 + L + 0 ) + tL + ε + tL + 0 .
Similarly, when k = 2, the system (11) is replaced by:

                                           K 10 = K 01 = t + t(K 10 + L + 0 ), K 1100 = tK 10 + t(K 1100 + L + 100 ), K 1010 = t(K 10 + K 01 ) + t(K 1010 + L + 010 ), K 0110 = tK 01 + t(K 0110 + L + 110 ), L + ε = 2t(1 + L + ε ) 2 + tL + ε + 2t(L + 0 + L + 100 + L + 010 + L + 110 ), L + 0 = L + 1 = 2t(1 + L + ε )(K 10 + L + 0 ) + tL + ε + tL + 0 + t(L + 100 + L + 010 + L + 110 ), L 00 = L 11 = 2t(1 + L + ε )L 00 + t(K 10 + L + 0 ) + tL 00 + tL + 100 , L + 10 = L + 01 = 2t(1 + L + ε )(K 10 + L + 10 ) + t(K 01 + L + 1 ) + tL + 10 + t(L + 110 + L + 010 ), L + 100 = 2t(1 + L + ε )(K 1100 + L + 100 ) + tL + 10 + tL + 100 , L + 010 = 2t(1 + L + ε )(K 1010 + L + 010 ) + t(L + 01 + L + ε K 10 ) + tL + 010 , L + 110 = 2t(1 + L + ε )(K 0110 + L + 110 ) + t(L 11 + L + ε K 10 ) + tL + 110 .
The second condition that we need is properness (again, in the sense of [30, Sec. VII.6.3]). But the system S 1 that we have just obtained is proper, thanks to the factor t occurring in the right-hand side of ( 5), ( 6) and [START_REF] Boulatov | The Ising model on a random planar lattice: the structure of the phase transition and the exact critical exponents[END_REF].

The next condition is aperiodicity. The coefficients of t 1 and t 2 in the series L + ε (t) are both non-zero. This implies that this series is aperiodic. Consequently, if we prove that the system S 1 is irreducible (which will be our final objective below), then it will be aperiodic [30, p. 483].

So let us finally prove that S 1 is irreducible. Recall that in such a polynomial system, a series F depends on a series G if G occurs in the right-hand side of the equation defining F . Irreducibility means that the digraph of dependences is strongly connected. Recall that S 1 involves two families of series: the series K w , for w balanced of length between 2 and 2k, and L w (or L + w ) for any valid w of length at most 2k -1. To lighten notation, for any w we denote by Lw the corresponding L-series, be it L w or L + w . Let us first prove that every series in S 1 depends on Lε = L + ε . By ( 6) and ( 13), this holds for every Lw . Now, each K w depends on at least one L-series (see ( 5)), and thus by transitivity on L + ε .

Conversely, let us prove that L + ε depends on all other series occurring in S 1 . • First, Equation ( 13) applied to w = ε shows that L + ε depends on all series Lu such that u is quasi-balanced.

• Let us now prove, by induction on the balance b(u), that L + ε depends on Lu for each valid word u of length at most 2k-1. We have already seen this for b(u) = 1. If b(u) = 0, then |u| ≤ 2k -2, and for any letter a the word w := ua is valid and quasi-balanced. The second term in [START_REF] Boulatov | The Ising model on a random planar lattice: the structure of the phase transition and the exact critical exponents[END_REF] shows that Lw depends on Lu . By transitivity, this implies that L + ε depends on Lu . We have thus set the initial cases of our induction, for balances 0 and 1. Now assume b(u) ≥ 2. There exists a letter a such that w := ua is valid and has balance b(u) -1. If b(w) = 1 (resp. b(w) > 1), the second (resp. third) term in the equation ( 13) (resp. ( 6)) defining Lw shows that Lw depends on Lu . By the induction hypothesis, L + ε depends on Lw , and thus by transitivity on Lu . • Finally, let u be balanced of length between 2 and 2k. Then w := u s is quasi-balanced.

The first term of the equation ( 13) defining Lw involves L w = K u + Lw , so that by transitivity, L + ε depends on K u . We have now checked all conditions of Theorem VII.6 of [30, p. 489]. Applying it gives our proposition.

Back to examples

We now return to the cases k = 1 and k = 2 studied in Section 3.2. We refer to [START_REF] Flajolet | Analytic combinatorics[END_REF]Sec. VII.7] for generalities on the singularities of algebraic series, and on the asymptotic behaviour of their coefficients. When k = 1, we have obtained for L ε the quadratic equation [START_REF] Bonichon | Baxter permutations and plane bipolar orientations[END_REF]. Its dominant coefficient only vanishes at t = 0, and its discriminant is ∆ 1 (t) := t 4 + 4t 3 + 22t 2 -12t + 1. The radius ρ 1 must be one of the roots of ∆ 1 . The only real positive roots are around 0.1032 and 0.3998. By solving (10) explicitly, we see that the smallest of these roots is indeed a singularity of L ε . Hence ρ 1 = 0.1032 . . . and the corresponding growth rate is λ 1 = 1/ρ 1 = 9.684 . . ., which improves on the lower bound 8 coming from Eulerian maps.

When k = 2, we have obtained for L ε the quartic equation [START_REF] Borot | Loop models on random maps via nested loops: case of domain symmetry breaking and application to the Potts model[END_REF]. Its dominant coefficient does not vanish away from 0, and its discriminant is

∆ 2 (t) := 64t
239969t 7 + 723151t 6 -1106764t 5 + 820832t 4 -316644t 3 + 65424t 2 -6780t + 268).
The only roots in (0, 1) are 0.0984 . . . and 0.2714 . . .. The radius ρ 2 must be the first one (the other would give a growth rate smaller than 8). Hence the corresponding growth rate is λ 2 = 1/ρ 2 = 10.16 . . ., which improves on the previous bound λ 1 .

We do not push our study to larger values of k, as we will obtain better bounds with the prime decomposition in the next section.

Subsets of Eulerian orientations, via the prime decomposition

In this section, we combine the restriction on allowed splits of the previous section with the prime decomposition of Section 2.4 to obtain a new family of subsets of Eulerian orientations. The results and proofs are similar to those of the previous section, and we give fewer details. The new subsets L (k) satisfy L (k) ⊂ L (k) (Proposition 8), hence they give better lower bounds on the growth rate µ than those obtained in the previous section. Moreover these bounds increase to µ as k increases (Proposition 7).

Recall from Section 2.4 that an Eulerian orientation is a sequence of prime Eulerian orientations, and that a prime (Eulerian) orientation can be obtained recursively from the atomic map by either:

• adding a loop, oriented in either way, around an orientation O 1 ,

• or performing a legal split on a prime orientation O ∈ O, followed by the concatenation of an arbitrary Eulerian orientation O at the new vertex created by the split (Figure 6). Definition 6. Let k ≥ 1. Let L (k) be the set of planar orientations obtained recursively from the atomic map by either:

• concatenating a sequence of prime orientations of L (k) ,

• or adding a loop, oriented in either way, around an orientation O 1 of L (k) ,

• or performing a legal i-split on a prime orientation O ∈ L (k) , with i = dv(O ) or i ≤ k, followed by the concatenation of an arbitrary orientation of L (k) at the new vertex created by the split.

Clearly, the sets L (k) increase to the set O of all Eulerian orientations as k increases, hence their growth rates λk form a non-decreasing sequence of lower bounds on µ. But we have in this case a stronger result. Proposition 7. For k ≥ 1, the sequence ( ¯ (k) n ) n≥0 that counts orientations of L (k) by their size is super-multiplicative. Consequently, the associated growth rate λk := lim

n ( ¯ (k) n ) 1/n = sup n ( ¯ (k) n ) 1/n ( 14
)
increases to µ as k tends to infinity.

Proof. By definition of L (k) , concatenating two orientations of L (k) at their root vertex gives a new element of L (k) , which implies super-multiplicativity and the identity (14) (by Fekete's Lemma [59, p. 103]). Now since L (k) converges to O, for any n, there exists k such that o n = (k) n (one can take k = n, or even k = n -2). Hence

o 1/n n = ( ¯ (k) n ) 1/n ≤ λk ≤ lim k λk ,
and it follows now from (2) that µ ≤ lim λk . Since λk ≤ µ, the proposition follows.

Proposition 8. For k ≥ 1, the subset of orientations L (k) includes the subset L (k) defined in Section 3.

Proof. We prove this by induction on the number of edges. The inclusion is obvious for orientations with no edge. Now let O ∈ L (k) , having at least one edge.

If O is the merge of two orientations O 1 and O 2 , then the induction hypothesis implies that O 1 and O 2 are in L (k) . The structure of L (k) implies that every prime sub-orientation of O 2 (attached at the root of O 2 ) also belongs to L (k) . Then O can be obtained as an orientation of L (k) by first adding a loop around O 1 (this is the second construction in Definition 6), then concatenating one by one the prime sub-orientations of O 2 (first construction in Definition 6).

Otherwise, O is obtained by a legal split in an orientation O formed of the prime suborientations P 1 , . . . , P . By the induction hypothesis, O , and its prime sub-orientations P 1 , . . . , P , belong to L (k) . Let us say that the split occurs in P i (this means that the sub-orientations P 1 , . . . , P i-1 are attached to the new created vertex v , while P i+1 , . . . , P remain attached to the original vertex v, the root vertex of O). Then the orientation O 1 obtained by deleting from O the sub-orientations P i+1 , . . . , P can be obtained by a legal split in the prime orientation P i , followed by the concatenation of P 1 , . . . , P i-1 at the new created vertex. This is the third construction in Definition 6, hence O 1 belongs to L (k) . It remains to concatenate P i+1 , . . . , P at the root (first construction in Definition 6), and we recover O as an element of L (k) . 4.1. An algebraic system for L (k) We now fix k ≥ 1. For w a word on {0, 1}, let L (k) w (t) ≡ L w denote the generating function of orientations of L (k) whose root word ends with w. Let K w be the generating function of those that have root word exactly w. Let L w and K w be the corresponding series for prime orientations of L (k) . We are especially interested in the series L ε that counts all orientations of L (k) .

If w is not balanced, K w = K w = 0, while if w = ε, K w = 1 and K w = 0. For w non-empty and balanced, of length at most 2k, we have

K w = w=uv K u K v , (15)
since an orientation of L (k) is a sequence of orientations of L (k) . Now the description of prime orientations of L (k) (Definition 6) gives

K w = tK wc + tL ε L ws . (16)
The first term corresponds to adding a loop, and the second to a legal i-split, where i ≤ k is the half-length of w. The factor L ε accounts for the orientation O attached at the end of the root edge. Now let w be a valid word of length at most 2k -1, and let us write equations for the series L w and L w . For L w , the sequential structure of orientations of L (k) gives

L w = 1 w=ε + L ε L w + w=uv,v =w L u K v . ( 17
)
The second (resp. third) term counts orientations in which the root word of the last prime component ends with w (resp. is shorter than w). Finally, for the series L w we obtain the following counterpart of (6):

L w = 2tL ε 1 w=ε + tL wp 1 w =ε + tK wc 1 w =ε balanced + tL ε      L w + u=vw 0<|u|≤2k u balanced (L us -K u )      . (18)
The first three terms count orientations in which the root edge is a loop, and the last one those obtained by a split. Proposition 9. Consider the collection of equations consisting of:

• Equation (15), written for all balanced words w of length between 2 and 2k -2,

• Equation [START_REF] Bousquet-Mélou | The degree distribution of bipartite planar maps: applications to the Ising model[END_REF], written for all balanced words w of length between 2 and 2k, • Equation [START_REF] Bouttier | Counting colored random triangulations[END_REF], written for all valid words w of length at most 2k -2,

• Equation (18), written for all valid words w of length at most 2k -1.

In this collection, replace all trivial Kand K -series by their value: K w = K w = 0 when w is not balanced, K ε = 1, K ε = 0. Let S 0 denote the resulting system. The number of series it involves is 2f (k) -2 2k k , where f (k) is given by [START_REF] Bernardi | Counting coloured planar maps: differential equations[END_REF]. Moreover, S 0 defines uniquely all these series. Its size can be (roughly) divided by two upon exploiting the 0/1 symmetry. Proof. To prove that all series are well defined by the system, we first check that every series occurring in the right-hand side of some equation is the left-hand side of another equation. Then we note that:

• the equations for prime orientations, namely ( 16) and [START_REF] Bouttier | Blocked edges on Eulerian maps and mobiles: application to spanning trees, hard particles and the Ising model[END_REF], have a factor t in their right-hand sides, • for the other two equations, ( 15) and ( 17), every non-trivial term in the right-hand side has a series of prime orientations as a factor. Now the number of equations: every series that was occurring in the system S 0 of Proposition 3 now has two copies (one with a prime, one without), except for the series K w , for w balanced of length 2k, and L w , for w quasi-balanced of length 2k -1, which have only one copy. Since there are 2k k balanced words of length 2k, and 2 2k-1 k = 2k k quasi-balanced words of length 2k -1, the result follows.

Remark 10. As in Remark 4, if w is such that 0w and 1w are both valid of length less than 2k -2 (resp. 2k -1), we can replace (17) (resp. ( 18)) by the simpler forward equation:

L w = K w + L 0w + L 1w (resp. L w = K w + L 0w + L 1w ).
This does not increase the size of the system.

Examples

4.2.1.

When k = 1, the system S 0 contains 2(f (1) -2) = 6 equations, or 4 of we exploit the 0/1 symmetry:

         K 10 = t + tL ε L 0 , L ε = 1 + L ε L ε , L ε = 2tL ε + tL ε (L ε + 2L 0 -2K 10 ), L 0 = tL ε + tL ε (L 0 + L 0 -K 10 ).
Eliminating all series but L ε gives a cubic equation for the generating function L ε ≡ L

(1) ε of Eulerian orientations in L (1) :

t 2 L 3 ε + t(t -4)L 2 ε + (2t + 1)L ε -1 = 0. (19)
4.2.2. When k = 2, the system S 0 contains 2(f (2) -6) = 30 equations, or 16 if we exploit the 0/1 symmetry:

                                                                   K 01 = K 10 = K 01 , K 10 = K 01 = t + tL ε L 1 , K 1100 = tK 10 + tL ε L 100 , K 1010 = tK 01 + tL ε L 010 , K 0110 = tL ε L 110 , L ε = 1 + L ε L ε , L 0 = L 1 = L ε L 0 , L 00 = L 11 = L ε L 00 , L 01 = L 10 = L ε L 01 , L ε = 2tL ε + tL ε (L ε + 2(L 0 -K 10 + L 100 -K 1100 + L 010 -K 1010 + L 110 -K 0110 )), L 0 = L 1 = tL ε + tL ε (L 0 + L 0 -K 10 + L 100 -K 1100 + L 010 -K 1010 + L 110 -K 0110 ), L 00 = tL 0 + tL ε (L 00 + L 100 -K 1100 ), L 10 = L 01 = tL 1 + t + tL ε (L 10 + L 1 -K 01 + L 010 -K 1010 + L 110 -K 0110 ), L 100 = tL 10 + tL ε (L 100 + L 100 -K 1100 ), L 010 = tL 01 + tL ε (L 010 + L 010 -K 1010 ), L 110 = tL 11 + tL ε (L 110 + L 110 -K 0110 ). (20)
Eliminating all series but L ε gives an equation of degree 6 for the generating function

L ε ≡ L (2) ε of Eulerian orientations in L (2) : 2t 5 L 6 ε -t 4 (t+8)L 5 ε -t 3 (3t 2 -16)L 4 ε +t 2 (2t+3)(2t-5)L 3 ε -t(2t 2 -7t-7)L 2 ε -(5t+1)L ε +1 = 0. (21)

Asymptotic analysis for subsets of Eulerian orientations (prime decomposition)

We now prove for the polynomial system of Proposition 9 an analogue of Proposition 5.

Proposition 11. For k ≥ 1, let ρk denote the radius of convergence of the series L (k) ε that counts orientations of L (k) . Then ρk is the only singularity of L (k) ε of minimal modulus, and it is of the square root type. Consequently, there exists a constant c such that the number ¯ (k) n of orientations of size n in L (k) satisfies, as n tends to infinity:

¯ (k) n ∼ c λn k n -3/2 , with λk = 1/ρ k .
Proof. Again, we apply the theory of positive irreducible polynomial systems [30, Sec. VII.6].

The system of Proposition 9 is not positive. To correct this, we replace the series L w (for w balanced) and L w (for w balanced or quasi-balanced) by their "positive" versions:

L + w := L w -K w , L + w := L w -K w (w balanced), L + w := L w -K ← - w (w quasi-balanced).
In particular, L ε is replaced by L + ε := L ε -1 and L + ε coincides with L ε . We alter the original system S 0 as follows:

(i) For w balanced, we replace the equation ( 17) defining L w by the difference between ( 17) and ( 15):

L + w = L + ε L w + L + w + w=uv,v =w L u K v . (22)
(ii) For w balanced or quasi-balanced, we replace the equation ( 18) defining L w by the difference between ( 18) and ( 16) (written for w if w is balanced, for ←w otherwise):

L + w = 2tL ε 1 w=ε + t(L wp -K wp )1 w =ε + tL ε       L + w + u=vw,u =w |u|≤2k-1 u quasi-balanced L + u       . ( 23
)
(iii) In the new system thus obtained, we replace every series K w such that w is not balanced by 0, every series L w (resp. L w ) such that w is balanced by K w + L + w (resp. K w + L + w ), and every series L w such that w is quasi-balanced by K ← - w + L + w . In particular, the series L wp -K wp occurring in [START_REF] Di Francesco | 2D gravity and random matrices[END_REF] becomes L + wp when w p is balanced, L wp otherwise. The series L w (resp. L w ) that remain in the system are such that w has balance at least 1 (resp. 2).

We thus obtain a positive system, denoted S 1 , defining the following series:

• K w , for w balanced of length between 2 and 2k -2, • K w , for w balanced of length between 2 and 2k, • L w for w valid of balance at least 1 and length at most 2k -2, • L w for w valid of balance at least 2 and length at most 2k -1, • L + w , for w balanced of length at most 2k -2, • L + w , for w balanced or quasi-balanced of length at most 2k -1.

For instance, when k = 1 we obtain the following system:

         K 10 = t + t(1 + L + ε )(K 10 + L + 0 ), L + ε = L + ε L ε + + L + ε , L + ε = 2t(1 + L + ε ) + t(1 + L + ε )(L ε + + 2L + 0 ), L + 0 = tL + ε + t(1 + L + ε )L + 0 .
Recall that the series we are interested in is L + ε . But then we can drop the first equation of the above system. This size reduction occurs for any value of k, and the positive system S 2 that we will study is finally obtained by performing one last change:

(iv) Delete the equations defining the series K w , for w of length 2k.

Observe that all the series involved in S 2 are well-defined by this system. This comes from the fact that all the series K u , for u of length 2k, that occurred in S 0 came from the term L us -K u of [START_REF] Bouttier | Blocked edges on Eulerian maps and mobiles: application to spanning trees, hard particles and the Ising model[END_REF], which now reads L + us .

Here is for instance the system obtained for k = 2, which has three equations less than [START_REF] Brown | On the existence of square roots in certain rings of power series[END_REF]:

                                                     K 01 = K 10 = K 01 , K 10 = K 01 = t + t(1 + L + ε )(K 01 + L + 1 ), L + ε = L + ε L + ε + L + ε , L 0 = (1 + L + ε )(K 10 + L + 0 ), L 00 = L 11 = (1 + L + ε )L 00 , L + 01 = L + ε (K 01 + L + 01 ) + L + 01 , L + ε = 2t(1 + L + ε ) + t(1 + L + ε )(L + ε + 2(L + 0 + L + 100 + L + 010 + L + 110 )), L + 0 = L + 1 = tL + ε + t(1 + L + ε )(L + 0 + L + 100 + L + 010 + L + 110 ), L 00 = tL 0 + t(1 + L + ε )(L 00 + L + 100 ), L + 10 = L + 01 = tL 1 + t(1 + L + ε )(L + 10 + L + 010 + L + 110 ), L + 100 = tL + 10 + t(1 + L + ε )L + 100 , L + 010 = tL + 01 + t(1 + L + ε )L + 010 , L + 110 = tL 11 + t(1 + L + ε )L + 110 .
Let us now discuss properness [30, p. 489]. The system S 2 that we have just obtained is not proper. However, the right-hand sides of the equations that define series with a prime (K , L and L + ) are multiples of t (see ( 16), ( 18) and ( 23)). In the remaining equations, that is ( 15), ( 17) (for w not balanced) and ( 22) (for w balanced), each term on the right-hand side involves a series with a prime: hence after one iteration of S 2 , one obtains a new system S 3 which is positive and proper.

Aperiodicity holds as in the previous section, and we are left with irreducibility. Note that proving irreducibility for S 2 or its iterated version S 3 is equivalent, so we focus on S 2 . As in the previous section, we denote by Lw the series L + w or L w , depending on whether w is balanced or not. Similarly, L w denotes L + w if w is balanced or quasi-balanced, and L w otherwise. Let us prove that all series depend on Lε . We first observe that this holds for every K -or Lor L -series (see ( 16), ( 17), ( 18), ( 22), ( 23)). We are left with the series K w : but it depends on K w (see [START_REF] Bousquet-Mélou | Enumeration of planar constellations[END_REF]), and hence on Lε .

Conversely, let us prove that Lε depends on every other series in the system. By [START_REF] Daul | q-States Potts model on a random planar lattice[END_REF], it depends on L ε . Then by [START_REF] Di Francesco | 2D gravity and random matrices[END_REF] applied to w = ε, it depends on every series L u , where u is quasibalanced. Going back and forth between the equations defining the Lseries and the L -series (see [START_REF] Bouttier | Counting colored random triangulations[END_REF], ( 18), ( 22), ( 23)), and using an induction on the balance, we then see that Lε depends on all series Lu (for |u| ≤ 2k -2) and all series L u (for |u| ≤ 2k -1). Then the first term of [START_REF] Daul | q-States Potts model on a random planar lattice[END_REF], written as L + ε (K w + L + w ), shows that Lε depends on all series K v with |v| ≤ 2k -2. It remains to prove that Lε depends on the K-series. Let u = au s be balanced of length at most 2k -2, and define w = u s ā. This word has balance 2. The second term of [START_REF] Bouttier | Blocked edges on Eulerian maps and mobiles: application to spanning trees, hard particles and the Ising model[END_REF] involves

L wp = L us = K u + L + us .
Hence L w depends on K u , and by transitivity, Lε depends on K u . This proves the irreducibility of the system and concludes the proof of the proposition.

Back to examples

We first return to the cases k = 1 and k = 2 studied in Section 4.2. When k = 1, we obtained the cubic equation [START_REF] Brak | Algebraic approximants: a new method of series analysis[END_REF] for L [START_REF] Alm | Random self-avoiding walks on one-dimensional lattices[END_REF] ε . The discriminant has three positive roots, which are 1, and (approximately) 0.094 and 15.9. The second one is the radius of convergence, and we obtain the lower bound λ1 10.603 on the growth rate of Eulerian orientations. This improves significantly on the growth rate λ 1 = 9.68 . . . obtained from the set L (1) .

For k = 2, we obtained the equation ( 21) satisfied by L

(2)

ε . The discriminant has two roots in (0, 1), which are approximately 0.0911 and 0.414. The first one is the radius of convergence, and we obtain the lower bound λ2 10.9759 on the growth rate of Eulerian orientations.

When k = 3, we find that L (3) ε satisfies an equation of degree 20 (see the Maple sessions available on our web pages). The dominant coefficient only vanishes at t = 8, and the discriminant has only one relevant root, around 0.089. This gives the lower bound λ3 11.2289 on the growth rate of Eulerian orientations.

When k = 4, we did not compute the equation satisfied by L (4)

ε , but we estimated λ4 from the first 30 coefficients of L (4) ε using quadratic approximants [START_REF] Brak | Algebraic approximants: a new method of series analysis[END_REF]. We predict λ4 11.41. This value has then been confirmed by Bruno Salvy using the Maple package NewtonGF [START_REF] Pivoteau | Algorithms for combinatorial structures: well-founded systems and Newton iterations[END_REF], with which he obtained 10 digits of λ4 . This package also allows us to compute more coefficients in L (4) ε . Moreover, Jean-Charles Faugère [26] has finally been able to determine the equation for L

Supersets of Eulerian orientations, via the standard decomposition

We now want to define, and count, supersets of Eulerian orientations. Their generating functions will be described by functional equations involving divided differences (as in ( 3)). The proof of their algebraicity is non-trivial, relying on a deep result from Artin's approximation theory (Theorem 16).

Recall that Eulerian orientations can be obtained recursively from the atomic map by either:

• the merge of two orientations O 1 , O 2 ∈ O (with the root loop oriented in either way),

• or a legal split on an orientation O ∈ O. We now define the sets U (k) . The idea is that we allow illegal i-splits, provided i is larger than k. Definition 12. Let k ≥ 1. Let U (k) be the set of planar orientations obtained recursively from the atomic map by either:

• the merge of two orientations O 1 , O 2 ∈ U (k) (with the root loop oriented in either way),

• or a legal i-split on a map O ∈ U (k) with i ≤ k (small split),

• or an arbitrary split on a map O ∈ U (k) with i > k (large split). If the split is legal, the root edge is oriented in the only way that makes the new orientation Eulerian. Otherwise, it is oriented away from the root vertex.

Observe that all Eulerian orientations belong to U (k) . Moreover, the sets U (k) form a decreasing sequence, as fewer illegal splits are performed as k grows. Finally, for k ≥ n (and even for k ≥ n -2), all orientations of size n in U (k) are Eulerian. Hence the limit of the sets U (k) is the set O of all Eulerian orientations.

Another important observation is that, if the root vertex of an orientation of U (k) has degree at most 2k, then the root word of this orientation is balanced.

Functional equations for U (k)

We now fix an integer k. For a word w on {0, 1}, let U (k) w (t; x) ≡ U w (x) denote the generating function of orientations of U (k) whose root word ends with w, counted by the edge number (variable t) and the half-degree of the root vertex (variable x). Let T (k) w (t) ≡ T w denote the generating function of orientations of U (k) having root word exactly w. We do not record in this series the root degree (which is the length of w). To lighten notation, we often denote simply by U w the edge generating function U w (1) ≡ U w (t, 1), and by U x w the refined generating function U w (x) ≡ U w (t; x).

Note that T w = 0 if w is not balanced and that T ε = 1. Now, for w balanced of length between 2 and 2k, we have

T w = t auāv=w T u T v + tU ws . ( 24
)
The first term counts orientations obtained by a merge. The second one counts those obtained by a split, which is necessarily small since we have assumed |w| ≤ 2k. Note the analogy with [START_REF] Bernardi | Intervals in Catalan lattices and realizers of triangulations[END_REF]. For w valid of length at most 2k -1, let us now prove the following identity:

U x w = 1 w=ε + 2txU x ε U x w + tx w=uav U x u x |v|/2 T v + tx |w|/2 w=auāv T u T v + t u=vw 2≤|u|≤2k u balanced x |u|/2 U us + tx x -1 U x w -x k U w - tx x -1 u=vw |u|≤2k-2 T u (x |u|/2 -x k ). (25)
The first line is similar to the first line of ( 6): it counts the atomic map and orientations obtained from a merge. The only difference is that we now record the root degree. On the second line, the first sum counts orientations obtained by a small split (with root word u). Let us explain the remaining terms, which count orientations obtained by a large split, legal or not, of an orientation O whose root word ends (necessarily) with w. Given an orientation O with root vertex degree 2d, with d > k, the generating function of orientations obtained from O by a large split is

t 1+e(O ) x k+1 + x k+2 + • • • + x d = t 1+e(O ) x d+1 -x k+1
x -1 .

Let us underline that we cannot apply a large split to an orientation O whose root word u satisfies |u| ≤ 2k. Hence the generating function of orientations obtained by a large split is

tx x -1     U x w - u=vw,|u|≤2k x |u|/2 T u   -x k   U w - u=vw,|u|≤2k T u     ,
which gives the last two terms of (25) (the terms T u with |u| = 2k do not contribute).

Remark 13. In the proof of (25), we have tried to follow the same steps as in the proof of (6). However, comparing [START_REF] Drmota | Universal exponents and tail estimates in the enumeration of planar maps[END_REF] and ( 25) suggests to replace (25) by a lighter equation:

U x w = x |w|/2 T w + 2txU x ε U x w + tx w=uav U x u x |v|/2 T v + t u=vw |u|≤2k-1 u quasi-balanced x (1+|u|)/2 U u + tx x -1 U x w -x k U w - tx x -1 u=vw |u|≤2k-2 T u (x |u|/2 -x k ). ( 26
)
Proposition 14. Consider the collection of equations consisting of:

• Equation (24), written for all balanced words w of length between 2 and 2k -2,

• Equation (26), written for all valid words w of length at most 2k -1.

In this collection, replace all trivial T -series by their value: T w = 0 when w is not balanced, T ε = 1. Let R 0 denote the resulting system. The number of series it involves is f (k) -2k k , where f (k) is given by [START_REF] Bernardi | Counting coloured planar maps: differential equations[END_REF]. Moreover, R 0 defines uniquely these series. Its size can be (roughly) divided by two upon exploiting the 0/1 symmetry.

The proof is similar to the proof of Proposition 3.

Remark 15. As in Remark 4, if w is such that 0w and 1w are both valid of length less than 2k, we can replace (26) by the simpler forward equation:

U x w = x |w|/2 T w + U x 0w + U x 1w .
This does not increase the size of the system.

Examples

5.2.1.

When k = 1, the system of Proposition 14 contains f (1) -2 = 3 equations. Upon exploiting the 0/1 symmetry, it reads:

U x ε = 1 + 2tx(U x ε ) 2 + 2txU 0 + tx x-1 (U x ε -xU ε ) + tx, U x 0 = 2txU x ε U x 0 + txU x ε + txU 0 + tx x-1 (U x 0 -xU 0 ). (27)
The first equation can be replaced by the forward equation U x ε = 1 + 2U x 0 . We explain in Section 5.4 how to solve this system. 5.2.2. When k = 2, the system of Proposition 14 contains f (2) -6 = 15 equations. Upon exploiting the 0/1 symmetry and (some) forward equations, it reads:

                                 T 10 = t + tU 0 , U x ε = 1 + 2U x 0 , U x 0 = U x 1 = 2txU x ε U x 0 + txU x ε + txU 0 + tx 2 (U 100 + U 010 + U 110 ) + tx x-1 (U x 0 -x 2 U 0 ) + tx 2 T 10 , U x 10 = U x 01 = xT 10 + U x 110 + U x 010 , = xT 10 + 2txU x ε U x 10 + txU x 1 + tx 2 (U 010 + U 110 ) + tx x-1 (U x 10 -x 2 U 10 ) + tx 2 T 10 , U x 00 = U x 11 = 2txU x ε U x 00 + txU x 0 + tx 2 U 100 + tx x-1 (U x 00 -x 2 U 00 ), U x 100 = 2txU x ε U x 100 + txU x 10 + tx 2 U 100 + tx x-1 (U x 100 -x 2 U 100 ), U x 010 = 2txU x ε U x 010 + tx (U x 01 + U x ε xT 10 ) + tx 2 U 010 + tx x-1 (U x 010 -x 2 U 010 ), U x 110 = 2txU x ε U x 110 + tx (U x 11 + U x ε xT 10 ) + tx 2 U 110 + tx x-1 (U x 110 -x 2 U 110 ). (28)
We explain in Section 5.4 below how to solve this system.

Algebraicity

Since the early work of Brown in the sixties on the quadratic method [START_REF] Brown | On the existence of square roots in certain rings of power series[END_REF], a lot has been known about equations involving divided differences of the form (F (t; x) -F (t; 1))/(x -1)). However, most of the literature deals with a single equation, not with a system [START_REF] Bousquet-Mélou | Polynomial equations with one catalytic variable, algebraic series and map enumeration[END_REF][START_REF] Goulden | Combinatorial enumeration[END_REF]. In order to prove that the series U (k) ε (t; x) that counts orientations of U (k) is algebraic, we use a deep theorem from Artin's approximation theory, due to Popescu [START_REF] Kurke | Die Approximationseigenschaft lokaler Ringe[END_REF][START_REF] Popescu | General Néron desingularization and approximation[END_REF]. The form we will need is given below. We recall that C[[z 1 , . . . , z ]] is the ring of formal power series in the variables z 1 , . . . , z , with complex coefficients, and that a series Z in this ring is algebraic if it satisfies a non-trivial polynomial equation Pol(z 1 , . . . , z , Z) = 0. Theorem 16 ([50], Thm. 1.4). Consider a polynomial system of n equations in + n variables over C, written as P i (z 1 , . . . , z , y 1 , . . . , y n ) = 0, for 1 ≤ i ≤ n. Let (d 1 , . . . , d n ) be a sequence of integers in {0, 1, . . . , }. Assume that there exists an n-tuple Y = (Y 1 , . . . , Y n ) of series in C[[z 1 , . . . , z ]] that satisfies the following conditions:

• the n-tuple Y solves this system, that is,

P i (z 1 , . . . , z , Y 1 , . . . , Y n ) = 0 for 1 ≤ i ≤ n,
• for 1 ≤ i ≤ n, the series Y i does not depend on the variables z j such that j > d i (if d i = , then there is no condition on the series Y i ). Then there exists an n-tuple (Z 1 , . . . , Z n ) of algebraic series in C[[z 1 , . . . , z ]] that solves the system and satisfies the same dependence conditions as Y.

In particular, if the system has a unique solution satisfying the dependence conditions, then this solution is algebraic.

An application. To our knowledge, this theorem has not been applied yet in a combinatorial context. So, before we use it to prove the algebraicity of U (k) ε (t; x), let us examine its application to a simple equation, namely [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF].

First, observe that the algebraicity of M (t; x) is not obvious. Clearly, if we could prove that M (t; 1) is algebraic, we would be done with M (t; x) as well, but why should M (t; 1) be algebraic? We can apply the above theorem as follows. Let us denote t = z 1 and x = 1 + z 2 (we shall explain later why we need to translate the variable x). We consider the system in z 1 , z 2 , y 1 and y 2 consisting of the following (single) equation:

z 2 y 2 = z 2 + z 1 z 2 (1 + z 2 )y 2 2 + z 1 (1 + z 2 )(y 2 -y 1 )
. Take d 1 = 1 and d 2 = 2. Then (3) shows that the pair (Y 1 , Y 2 ) := (M (t; 1), M (t; x)) solves the above equation. Moreover Y 1 = M (t; 1) is independent of z 2 = x -1, while Y 2 = M (t; x) depends on both variables z 1 and z 2 , in accordance with d 1 = 1 and d 2 = 2.

Let us now prove that there cannot be another solution of this system in the ring

C[[z 1 , z 2 ]] such that Y 1 is independent of z 2 .
First, setting z 2 = 0 in the equation shows that Y 1 must be the specialization of Y 2 at z 2 = 0. This, combined with the factor z 1 occurring in every non-initial term in the right-hand side, implies that the coefficient of z n 1 in Y 2 can be computed by induction of n, starting from the constant coefficient 1. Hence the uniqueness of (Y 1 , Y 2 ). The algebraicity of M (t; x) now follows from the above theorem.

Note that, if we had used z 2 = x instead of z 2 = x -1, we could not apply the last part of Theorem 16. The equation would read

(z 2 -1)y 2 = (z 2 -1) + z 1 (z 2 -1)z 2 y 2 2 + z 1 z 2 (y 2 -y 1
), but this equation has many solutions in the ring C[[z 1 , z 2 ]] of formal power series in z 1 = t and z 2 = x. For instance, one can take Y 1 = 0 and

Y 2 = 1 -x + tx -(1 -x + tx) 2 -4tx(1 -x) 2 2tx(1 -x) .
Theorem 16 tells us that at least one of these solutions is algebraic, but we need uniqueness to conclude that our solution is algebraic. The key point is that a series in C[[z 1 , z 2 ]] can always be specialized at z 2 = 0, but not at z 2 = 1.

We now apply Theorem 16 to the larger example of orientations of U (k) .

Proposition 17. For any k ≥ 1, the generating function

U (k) ε (t; x) that counts orientations of U (k) is algebraic.
Proof. Again, we take as variables z 1 = t and z 2 = x -1. For short, we denote z 2 by z. We consider the polynomial system consisting of the following equations, which mimic [START_REF] Drmota | Universal exponents and tail estimates in the enumeration of planar maps[END_REF] and (26). For w balanced of length between 2 and 2k -2,

A w = t auāv=w A u A v + tC ws ,
and for w valid of length at most 2k -1,

zB w = z(1 + z) |w|/2 A w + 2tz(1 + z)B ε B w + tz(1 + z) w=uav B u (1 + z) |v|/2 A v + tz u=vw |u|≤2k-1 u quasi-balanced (1 + z) (1+|u|)/2 C u + t(1 + z) B w -(1 + z) k C w -t(1 + z) u=vw |u|≤2k-2 A u ((1 + z) |u|/2 -(1 + z) k ), (29)
where A ε = 1. The variables A w , B w and C w play the role of the y i in Theorem 16. By construction, the series

A w := T w (t),
B w := U w (t; 1 + z), C w := U w (t; 1)

solve the system. Moreover, A w and C w do not depend on z 2 = z.

By Theorem 16, it suffices to prove that the system in A, B, C has a unique solution in C[[t, z]] satisfying these dependence relations to conclude that all our series T and U counting orientations are algebraic.

So assume that A w , B w and C w solve the system and satisfy the required dependences. Then by setting z = 0 in (29), we see that C w must be the specialization of B w at z = 0, for all w valid of length at most 2k -1. Then the form of the system implies that the coefficient of t n in all series can be computed by induction on n, the initial values being B ε = C ε = 1 + O(t) and A w = B w = C w = O(t) for w non-empty (recall that we have set A ε = 1). This proves the uniqueness of the solution (with the required dependences) and concludes the proof. [START_REF] Felsner | Lattice structures from planar graphs[END_REF]. In the first equation, replace U 0 by (U ε -1)/2 to obtain a single equation involving only U x ε = U ε (x) and U ε = U ε (1). For simplicity, we now drop the index ε. This equation reads:

Pol(U (x), U (1), t, x) = 0, with Pol(x 0 , x 1 , t, x) = (x -1)(-x 0 + 1 + 2txx 2 0 + tx(x 1 - 1
)) + tx(x 0 -xx 1 ) + tx(x -1). We apply Brown's quadratic method. Its principle is the following: if there exists a formal power series X ≡ X(t) such that Pol x0 (U (X), U (1), t, X) = 0,

then this series X must be a double root of the discriminant ∆(U (1), t, x) of Pol(x 0 , U (1), t, x) with respect to x 0 (the notation Pol x0 stands for the derivative of Pol with respect to its first variable). The proof is elementary (see [34, Sec. 2.9] or [START_REF] Bousquet-Mélou | Polynomial equations with one catalytic variable, algebraic series and map enumeration[END_REF]). Equation ( 30) reads

X = 1 + tX + 4tX(X -1)U (X),
and has a unique power series solution X(t), whose coefficients can be computed by induction from those of U (x) (we do not need to determine X, just to know that it exists). Thus X is a double root of ∆(U (1), t, x), and hence the discriminant in x of ∆ must vanish. This gives the following cubic equation for U (1) (see our Maple sessions):

64t 3 U (1) 3 + 2t(24t 2 -36t + 1)U (1) 2 + (-15t 3 + 9t 2 + 19t -1)U (1) + t 3 + 27t 2 -19t + 1 = 0. ( 31
)
The series U (1) has a unique positive singularity τ 1 , around 0.0765, which is a root of 216t 3 -81t 2 + 18t -1. This gives the upper bound µ 1 = 1/τ 1 = 13.0659 • • • on the growth rate of Eulerian orientations. Expanding the series near τ 1 (using for instance the Maple function algeqtoseries [START_REF] Salvy | Gfun: a Maple package for the manipulation of generating and holonomic functions in one variable[END_REF]) shows that it has a singularity in (1 -µ 1 t) 3/2 , as the generating function of many families of planar maps.

5.4.2.

The case k = 2. We now return to the system [START_REF] Felsner | Bijections for Baxter families and related objects[END_REF]. Observe that we can reduce it to a system of three equations defining the series U x ε , U x 10 and U x 100 :

     U x ε = 1 + 4txU x ε U x 0 + 2txU x ε + 2txU 0 + 2tx 2 (U 100 + U 10 ) + 2tx x-1 (U x 0 -x 2 U 0 ), U x 10 = xt(1 + U 0 ) + 2txU x ε U x 10 + txU x 0 + tx 2 U 10 + tx x-1 (U x 10 -x 2 U 10 ), U x 100 = 2txU x ε U x 100 + txU x 10 + tx 2 U 100 + tx x-1 (U x 100 -x 2 U 100 ), (32)
in which we still need to plug

U 0 = U ε -1 2 and U x 0 = U x ε -1 2 . ( 33
)
To solve this system, we could develop a matrix analogue of the quadratic method, where (30) would be replaced by the cancellation of the Jacobian of the system. However, we prefer a step by step approach here, among other reasons because our system is not generic (its Jacobian has a multiple root).

From now on, we lighten notation by denoting

A = U x ε , A 1 = U ε , B = U x 10 , B 1 = U 10 , C = U x
100 and C 1 = U 100 . We will determine three polynomial equations relating the onevariable series A 1 , B 1 and C 1 , and then eliminate B 1 and C 1 to obtain a polynomial equation satisfied by A 1 = U ε .

We now describe the various steps of our calculation, without giving the intermediate equations: we refer to our web pages for a Maple session where the calculations are performed.

The first equation of [START_REF] Fusy | Bijective counting of involutive Baxter permutations[END_REF], after injecting [START_REF] Fusy | Bijective counting of plane bipolar orientations and Schnyder woods[END_REF], involves only one x-dependent series, namely

A = U x ε = U ε (x).
Once the denominators are cleared out, the degree in A is 2, and we can apply the quadratic method of Section 5.4.1: the discriminant (in x) of a certain discriminant (in x 0 ) vanishes, and this gives a first equation between A 1 , B 1 and C 1 .

We then move to the second equation of [START_REF] Fusy | Bijective counting of involutive Baxter permutations[END_REF], which (after injecting [START_REF] Fusy | Bijective counting of plane bipolar orientations and Schnyder woods[END_REF]) involves two xdependent series, namely A and B. It is linear in the latter series, with coefficient:

1 -x + tx + 2tx(x -1)A. ( 34
)
This coefficient vanishes for a (unique) series in t, denoted X, satisfying

X = 1 + tX + 2tX(X -1)A 2 , with A 2 := U ε (X).
Replacing x by X in the second equation of [START_REF] Fusy | Bijective counting of involutive Baxter permutations[END_REF] gives another equation between X and A 2 , from which we compute

2t(A 1 -2B 1 )X 2 + (1 -t -2tA 1 )X = 1, (35)
A 2 = 2XB 1 X -1 -A 1 .
We now eliminate X and A 2 between the last two identities and the first equation of ( 32), specialized at x = X. This gives a second equation between our three main unknown series A 1 , B 1 and C 1 .

We finally consider the third equation of (32) (after injecting [START_REF] Fusy | Bijective counting of plane bipolar orientations and Schnyder woods[END_REF]), which now involves all three x-dependent series. It is linear in C, again with coefficient [START_REF] Goulden | Combinatorial enumeration[END_REF]. Setting x = X in this equation gives an expression of U 10 (X):

B 2 := U 10 (X) = XC 1 /(X -1).
We now get back to the second equation of [START_REF] Fusy | Bijective counting of involutive Baxter permutations[END_REF], differentiate it with respect to x and set x = X. Replacing B 2 and A 2 by the above expressions gives:

A 2 := ∂U ε ∂x (X) = 2 2C 1 t(2X -1)(X -1)A 1 -4Xt(2X -1)B 1 C 1 -B 1 t(X -1) -(t -1)(X -1)C 1 t(X -1) 2 (4C 1 X + X -1)
.

It remains to differentiate the first equation of [START_REF] Fusy | Bijective counting of involutive Baxter permutations[END_REF] with respect to x, specialize it at x = X, and plug the above values of A 2 , B 2 and A 2 to obtain one more equation between A 1 , B 1 , C 1 and X. Eliminating X thanks to (35) gives our third and last equation between A 1 , B 1 and C 1 .

From this system, we eliminate B 1 and C 1 , and obtain an equation of degree 27 for A 1 = U ε . Its dominant coefficient does not vanish away from 0, and its discriminant has three roots in [1/10, 1/16] (where we know that the radius must be found), respectively located around 0.07509, 0.07658 and 0.07727. Following numerically the branches that start from 1 at t = 0 shows that the radius of U ε is the second one, giving the upper bound µ 2 = 13.057 . . . on the growth rate µ of Eulerian orientations. From numerical estimates of the singular exponent, we predict that the series has again a "planar map" singularity in (1 -µ 2 t) 3/2 . This is known to hold for many series satisfying an equation with divided differences [START_REF] Drmota | Universal exponents and tail estimates in the enumeration of planar maps[END_REF]. This leads us to complete Proposition 17 as follows.

Conjecture 18. For every k, the algebraic series U (k) ε (t; 1) that counts orientations of U (k) has a unique dominant singularity τ k = 1/µ k which is of the planar map type: as t approaches τ k from below,

U (k) ε (t; 1) = c 0 + c 1 (1 -µ k t) + c 2 (1 -µ k t) 3/2 1 + o(1) with c 2 = 0.

Supersets of Eulerian orientations, via the prime decomposition

In this section, we combine the illegal large splits of the previous section with the prime decomposition of Section 2.4 to obtain a new family of supersets of Eulerian orientations. These new supersets U (k) satisfy U (k) ⊂ U (k) (Proposition 20), hence they give better bounds on the growth rate µ than those obtained from the standard decomposition. Many arguments are similar to those of the previous section, and we give fewer details.

Recall from Section 2.4 that an Eulerian orientation is a sequence of prime Eulerian orientations, and that a prime (Eulerian) orientation can be obtained recursively from the atomic map by either:

• adding a loop, oriented in either way, around an orientation O 1 ,

• or a legal split on a prime orientation O ∈ O, followed by the concatenation of an arbitrary Eulerian orientation O at the new vertex created by the split (Figure 6).

Definition 19. Let k ≥ 1. Let U (k) be the set of planar orientations obtained recursively from the atomic map by either:

• concatenating a sequence of prime orientations of U (k) ,

• or adding a loop, oriented in either way, around an orientation O 1 of U (k) ,

• or performing a legal i-split on a prime orientation O ∈ U (k) , with i ≤ k, followed by the concatenation of an arbitrary orientation O of U (k) at the new vertex created by the split (small split), • or performing an arbitrary i-split on a prime orientation O ∈ U (k) , with i > k, followed by the concatenation of an arbitrary orientation O of U (k) at the new vertex created by the split (large split). If the split is legal, then the new edge is given the only orientation that makes the root word balanced, otherwise the root edge is oriented away from the root vertex.

Again, the sets U (k) decrease to the set O of all Eulerian orientations as k increases, hence their growth rates μk form a non-increasing sequence of upper bounds on µ. We do not know if this sequence converges to µ. At any rate, the convergence appears to be rather slow, as shown by the estimates of μk in Table 1.

Proposition 20. For k ≥ 1, the superset of orientations U (k) is contained in the superset U (k) defined in Section 5.

Proof. We prove this by induction on the number of edges. The inclusion is obvious for orientations with no edges. Now let O ∈ U (k) , having at least one edge.

If O is prime and is obtained by adding a loop around a smaller orientation O 1 of U (k) (second construction in Definition 19), then O 1 belongs to U (k) by the induction hypothesis, and so does O (first construction in Definition 12).

Assume now that O is prime and is obtained by an i-split in a prime orientation O of U (k) , followed by the concatenation of an orientation O of U (k) at the new vertex (third or fourth construction in Definition 19). Then the orientation Õ obtained by concatenating O and O at their root belongs to U (k) (first construction in U (k) ) and hence to U (k) by the induction hypothesis. But then one can recover O by performing an i-split in Õ, which is allowed in Õ as it was allowed in O . This is the second construction in Definition 12, hence O is in U (k) .

Assume finally that O is obtained by concatenating a prime orientation P of U (k) and another orientation O 2 of U (k) (first construction in Definition 19). By the induction hypothesis, both P and O 2 are in U (k) . If the root edge of P is a loop, deleting it from P leaves an orientation O 1 which is in U (k) . Then we can reconstruct O by a merge of O 1 and O 2 as in the first construction of Definition 12. If the root edge of P is not a loop (Figure 8), then P was obtained by the third or fourth construction in Definition 19: allowed split in a prime orientation P of U (k) , followed by the concatenation of some O ∈ U (k) at the new vertex. Let Õ be obtained by concatenating O , P and O 2 (in counterclockwise order) at their roots. Then Õ is in U (k) , but also in U (k) by the induction hypothesis. Then O can be recovered by a split in Õ, which is allowed in Õ as it was allowed in P (the split may have been small in P and become large in Õ, because of the orientation O 2 , but the converse is not possible). This is the second construction in Definition 12, hence O is in U (k) .

O P O 2 O , P , O 2 ∈ U (k) O = P O 2 ∈ U (k) ∩ U (k) Õ = O P O 2 ∈ U (k) ∩ U (k) (ind. hyp.) + O 2 P { split in P in U (k) concatenation in U (k) concatenation in U (k) split in Õ in U (k) O O P O 2 O 2 Figure 8. Two constructions of the orientation O: (top) in U (k)
, via an i-split, (bottom) in U (k) , via a j-split, with j = i + dv(O 2 ).

Functional equations for U (k)

We now fix an integer k. For a word w on {0, 1}, let U (k) w (t; x) ≡ U w (x) denote the generating function of orientations of U (k) whose root word ends with w, counted by the edge number (variable t) and the half-degree of the root vertex (variable x). Let T (k) w (t) ≡ T w denote the generating function of orientations of U (k) having root word exactly w. We define analogous generating functions U w (x) and T w for prime orientations. As in the previous section, we often denote simply by U w (resp. U w ) the edge generating function U w (t; 1) (resp. U w (t; 1)), and by U x w (resp. U

x w ) the refined generating function U w (t; x) (resp. U w (t; x)). Note that T w = T w = 0 if w is not balanced and that T ε = 1, T ε = 0. For w balanced of length between 2 and 2k, we have both a sequential equation

T w = w=uv T u T v (36)
analogous to [START_REF] Bousquet-Mélou | Enumeration of planar constellations[END_REF], and an equation for prime orientations:

T w = tT wc + tU ε U ws , (37)
analogous to [START_REF] Bousquet-Mélou | The degree distribution of bipartite planar maps: applications to the Ising model[END_REF]. The factor U ε accounts for the orientation concatenated after a split. For w valid of length at most 2k -1, we have a sequential equation, analogous to (17) but taking care of the root degree:

U x w = 1 w=ε + U x ε U x w + w=uv,v =w U x u x |v|/2 T v . (38)
Finally, we have the following equation for prime orientations, which is the counterpart of ( 18) and involves ingredients of (25) for orientations obtained by a split:

U x w = 2txU x ε 1 w=ε + txU x wp 1 w =ε + tx |w|/2 T wc 1 w =ε balanced + tU ε      u=vw 2≤|u|≤2k u balanced x |u|/2 U us + x x -1 U x w -x k U w - x x -1 u=vw |u|≤2k-2 T u (x |u|/2 -x k )    . (39)
The first line counts orientations obtained by adding a loop, and the second those obtained by a split.

Remark 21. As in the previous section, we can use [START_REF] Hierholzer | Ueber die Möglichkeit, einen Linienzug ohne Wiederholung und ohne Unterbrechung zu umfahren[END_REF] to replace (39) by a slightly lighter equation:

U x w = x |w|/2 T w + 2txU x ε 1 w=ε + txU x wp 1 w =ε + tU ε      u=vw |u|≤2k-1 u quasi-balanced x (1+|u|)/2 U u + x x -1 U x w -x k U w - x x -1 u=vw |u|≤2k-2 T u (x |u|/2 -x k ))    . (40)
Proposition 22. Consider the collection of equations consisting of:

• Equation (36), written for all balanced words w of length between 2 and 2k -4, • Equation (37), written for all balanced words w of length between 2 and 2k -2,

• Equation [START_REF]The on-line encyclopedia of integer sequences[END_REF], written for all valid words w of length at most 2k -2,

• Equation (40), written for all valid words w of length at most 2k -1.

In this collection, replace all trivial Tand T -series by their value: T w = T w = 0 when w is not balanced, T ε = 1, T ε = 0. Let R 0 denote the resulting system. The number of series it involves is 2f (k) -3 2k k -2k-2 k-1 1 k>1 , where f (k) is given by [START_REF] Bernardi | Counting coloured planar maps: differential equations[END_REF]. Moreover, R 0 defines uniquely all these series. Its size can be (roughly) divided by two upon exploiting the 0/1 symmetry.

The proof is similar to the proofs of Propositions 3 and 9.

Remark 23. As always, we can alternatively write forward equations:

U x w = x |w|/2 T w + U x 0w + U x 1w , U x w = x |w|/2 T w + U x 0w + U x 1w .
6.2. Examples 6.2.1. When k = 1, the system of Proposition 22 contains 2f (1) -3 • 2 = 4 equations. Upon exploiting the 0/1 symmetry, it reads:

       U x ε = 1 + U x ε U x ε , U x ε = 2txU x ε + tU ε 2xU 0 + x x-1 (U x ε -xU ε ) , U x 0 = txU x ε + tU ε xU 0 + x x-1 (U x 0 -xU 0 ) . (41)
The second equation can be replaced by the forward equation U x ε = 2U

x 0 . We solve this system in Section 6.4. 6.2.2. When k = 2, the system of Proposition 22 contains 2f (2) -3 • 6 -2 = 22 equations. Upon exploiting the 0/1 symmetry and the forward equations, it reads:

                                                           T 01 = t + tU ε U 1 , U x ε = 1 + U x ε U x ε , U x 0 = U x ε U x 0 , U x 10 = U x 01 = U x ε U x 10 , U x 00 = U x 11 = U x ε U
x 00 , U (42) We explain in Section 6.4 below how to solve this system.

Algebraicity

The analogue of Proposition 17 holds for the supersets obtained via the prime decomposition. Proof. Again, the idea is to apply Theorem 16 to the system of Proposition 22, after writing x = 1 + z. The proof is roughly the same as that of Proposition 17: we define a polynomial system involving two variables, t and z, and six families of unknowns A w , A w , B w , B w , C w , C w . The equations they satisfy are those of Proposition 22, rewritten with

T w → A w , T w → A w , U x w → B w , U x w → B w , U w → C w , U w → C w .
In fact, the only series U w occurring in our system is U ε , so that the polynomial system we construct involves C ε , but no other C-series. The prescribed dependences are that the A, A , C and C series are independent of z. For instance, when k = 1 we convert (41) into:

     B ε = 1 + B ε B ε , B ε = 2t(1 + z)B ε + tC ε 2(1 + z)C 0 + 1+z z (B ε -(1 + z)C ε ) , B 0 = t(1 + z)B ε + tC ε (1 + z)C 0 + 1+z z (B 0 -(1 + z)C 0 )
. However, this is not sufficient, because this system does not imply that C ε is B ε at z = 0 (it does however imply that C ε is B ε at z = 0, and similarly for C 0 and B 0 ). Hence, in order to apply Popescu's theorem, we need to add to our collection of equations the case x = 1, w = ε of (38), namely C ε = 1 + C ε C ε . The rest of the argument mimics the proof of Proposition 17.

Back to examples

We first return to the system (41) obtained for k = 1. In the second equation, replace U x ε by 1/(1 -U x ε ), U ε by 1/(1 -U ε ) and U 0 by U ε /2. This gives a polynomial equation involving only U x ε and U ε , which can be solved by the quadratic method already used in Section 5.4. This gives for U ε a cubic equation. Getting back to U ε = 1/(1 -U ε ), we obtain for the generating function U

(1) ε of orientations in U (1) the same cubic equation [START_REF] Fusy | Transversal structures on triangulations: a combinatorial study and straight-line drawings[END_REF] as for orientations of U (1) . In fact, one can check that U (1) = U (1) . Of course, the upper bound on µ is μ1 = µ 1 = 13.0659 . . . Let us now solve the system (42) obtained for k = 2. It can be reduced to a system of three equations defining the series U ε , U 10 and U 100 :

         U x ε = 2txU x ε + tU ε xU ε + 2x 2 (U 100 + U 10 ) + x x-1 (U x ε -x 2 U ε ) , U
in which we inject

U ε = 1 1 -U ε , U x ε = 1 1 -U x ε and U x 10 = U x 10
1 -U x ε . ( 44
)
We lighten notation by denoting A = U [START_REF] Felsner | Bijections for Baxter families and related objects[END_REF]. Again, we refer to our web pages for the corresponding Maple session. The intermediate steps are as follows. We first apply the quadratic method to the first equation. We then turn to the second one. The equation satisfied by X is

X = 1 + t 1 -t -A 1 ,
(Note that it gives X explicitly in terms of A 1 , whereas we had a quadratic equation (35) in the previous case.) We then derive

A 2 := U ε (X) = 1 + t 1 -A 1 t -2B 1 (1 -A 1 )
.

We finally consider the third equation of (43) (after injecting (44)), and derive:

B 2 := U 10 (X) = (1 -A 2 )C 1 /t.
Then it follows from the second equation that: 3 .

A 2 := ∂U ε ∂x (X) = 2(1 -A 1 )(1 -t -A 1 ) 2 ((1 -t -A 1 )C 1 + 2(-1 + A 1 )B 2 1 + tB 1 ) (t -2B 1 (1 -A 1 ))
At the end, we obtain an equation of degree 28 for A 1 = U ε , and then for U ε . Its dominant coefficient does not vanish away from 0, and its discriminant has only one root in [1/10, 1/16] (where we know that the radius must be found), around 0.0766. This gives the upper bound μ2 = 13.047 . . . on the growth rate µ of Eulerian orientations. From numerical estimates of the singular exponent, we predict that the series has again a "planar map" singularity in (1 -μ2 t) 3/2 . For k = 3, 4 and 5, we have generated our systems of equations and computed the first 100 coefficients of U (k) ε (t; x). From this we get the estimates of the growth rates μk shown in Table 2. The singularity still appears to be in (1 -μk t) 3/2 . We conjecture that this holds for any k.

Conjecture 25. For every k, the algebraic series U (k) ε (t; 1) that counts orientations of U (k) has a unique dominant singularity τk = 1/μ k which is of the planar map type: as t approaches τk from below, U (k) ε (t; 1) = c 0 + c 1 (1 -µ k t) + c 2 (1 -μk t) 3/2 1 + o(1) for c 2 = 0.

Final comments

As mentioned at the end of the introduction, it might be easier to count Eulerian orientation of 4-valent maps. In such orientations, each vertex has two in-going and two outgoing edges, so that counting Eulerian orientations means solving the so-called ice-model on (random) 4-valent maps [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF]Chap. 8]. The number of Eulerian orientations of a 4-valent planar map is known to be a third of the number of proper 3-colourings of its dual [START_REF] Welsh | The Tutte polynomial[END_REF]. Thus counting these orientations is equivalent to counting 3-coloured planar quadrangulations. A number of enumerative results are already known about coloured maps. In particular, 3-coloured planar maps, and 3-coloured planar triangulations, have algebraic generating functions [START_REF] Bernardi | Counting colored planar maps: algebraicity results[END_REF][START_REF] Bouttier | Counting colored random triangulations[END_REF][START_REF] Tutte | A census of planar maps[END_REF][START_REF] Tutte | Chromatic sums revisited[END_REF]. More generally, qcoloured planar maps and triangulations have differentially algebraic generating functions. So this could be true for 3-coloured quadrangulations as well, and hence for Eulerian orientations of 4-valent maps. Several results on this problem appear in the physics literature, but there does not seem to be an explicit exact solution at the moment [START_REF] Kostov | Exact solution of the six-vertex model on a random lattice[END_REF][START_REF] Zinn-Justin | The six-vertex model on random lattices[END_REF].

To finish, let us recall two questions left open by this paper (beyond the enumeration of Eulerian orientations!).

• Do the growth rates μk of orientations of U (k) decrease to the growth rate µ of Eulerian orientations, or to a larger value? • In Sections 3.3 and 4.3, we have used a general result about positive irreducible systems of polynomial equations to prove that the generating functions of our subsets of Eulerian orientations have a square root singularity. Could one define a notion of positive irreducible system of polynomial equations with divided differences whose solutions would systematically exhibit a singularity in (1 -µt) 3/2 ? Hopefully this would apply to our supersets of orientations and prove Conjectures 18 and 25. A first step in this direction, applicable to a single equation, is achieved in [START_REF] Drmota | Universal exponents and tail estimates in the enumeration of planar maps[END_REF].

Figure 1 .

 1 Figure 1. A rooted Eulerian map and a rooted Eulerian orientation.

Figure 2 .

 2 Figure 2. Left: First values of o n , for n from 0 to 15 (entry A277493 of the OEIS [38]). Right: A plot of o n+1 /o n vs. 1/n, for n = 4, . . . , 14, suggests that the growth rate of Eulerian orientations, located at the intercept of the curve and the y-axis, is around 12.5.

Figure 4 .

 4 Figure 4. Construction of an Eulerian orientation: merge two Eulerian orientations (the loop can be oriented in either way), or split (legally) an Eulerian orientation. Observe how the root word changes.

Figure 5 .

 5 Figure 5. Decomposition of an Eulerian map M into prime Eulerian maps M 1 , M 2 , M 3 .

5. 4 .

 4 Back to examples 5.4.1. The case k = 1. Let us go back to System

1 =xT 10 + txU x 1 +

 11 txU x ε + tU ε xU 0 + x 2 (U 100 + U 010 + U 110 ) + x x-1 (U x 0 -x 2 U 0 ) + x 2 T 10 , tU ε x 2 U 010 + x 2 U 110 + x x-1 (U x 10 -x 2 U 10 ) + x 2 T 10 , U x 00 = txU x 0 + tU ε x 2 U 100 + x x-1 (U x 00 -x 2 U 00 ) , U x 100 = txU x 10 + tU ε x 2 U 100 + x x-1 (U x 100 -x 2 U 100 ) , U x 010 = txU x 01 + tU ε x 2 U 010 + x x-1 (U x 010 -x 2 U 010 ) , U x 110 = txU x 11 + tU ε x 2 U 110 + x x-1 (U x 110 -x 2 U 110 ) .

Proposition 24 .

 24 For any k ≥ 1, the generating function U

ε

  (t; x) that counts orientations of U (k) is algebraic.

x 10 =

 10 xt(1 + U ε U ε /2) + tx(U x ε -1)/2 + tU ε x 2 U 10 + x x-1 (U x 10 -x 2 U 10 ) , U x 100 = txU x 10 + tU ε x 2 U 100 + x x-1 (U x 100 -x 2 U 100 ) ,

x ε , A 1 =

 1 U ε , B = U x 10 , B 1 = U 10 , C = U x 100 and C 1 = U 100 , and we follow the steps used in Section 5.4.2 to solve System

  Construction of an Eulerian map with n edges: merge an ordered pair of Eulerian maps M 1 , M 2 with n 1 and n 2 edges (n 1 + n 2 = n -1) and add a loop, or make a split on an Eulerian map with n -1 edges. The new edge (here thicker) is the root edge of M . last 2i -1 edges incident to v in O differ by ±1 (edges are visited in counterclockwise order, starting from the root corner). The orientation of the root edge of O is then forced (Figure

  12 (t -1)(81t 21 + 1863t 20 + 11322t 19 + 38592t 18 + 101105t 17 + 226631t 16 + 393423t 15 + 532907t 14 + 665167t 13 + 719797t 12 + 454804t 11 + 355710t 10 + 360159t 9 -262135t 8 -

Acknowledgements. We are grateful to Alin Bostan, Jean-Charles Faugère, Tony Guttmann and Bruno Salvy for their help with the algebraic systems, various calculations and suggestions. MBM also thanks the organizers of the workshop "Approximation and Combinatorics" in 2015 in CIRM (Luminy, F) Herwig Hauser and Guillaume Rond, for very interesting discussions on the algebraic aspects of functional equations arising in map enumeration.

NB, MBM and CP were partially supported by the French "Agence Nationale de la Recherche", via grants JCJC EGOS ANR-12-JS02-002-01 (NB and CP) and Graal ANR-14-CE25-0014 (MBM).