text:
# Untitled Document

SPhT/02-093

Census of Planar Maps: From the One-Matrix Model Solution to a Combinatorial Proof

J. Bouttier1 bouttier@spht.saclay.cea.fr, P. Di Francesco2 philippe@spht.saclay.cea.fr and E. Guitter3 guitter@spht.saclay.cea.fr

Service de Physique Théorique, CEA/DSM/SPhT

Unité de recherche associée au CNRS

CEA/Saclay

91191 Gif sur Yvette Cedex, France

We consider the problem of enumeration of planar maps and revisit its one-matrix model solution in the light of recent combinatorial techniques involving conjugated trees. We adapt and generalize these techniques so as to give an alternative and purely combinatorial solution to the problem of counting arbitrary planar maps with prescribed vertex degrees.

07/02

1. Introduction

Enumeration of planar maps has been a classical subject of combinatorics originally motivated by the famous four-color problem. Major advances in this field were obtained in the 60’s by W. Tutte in his famous “Census” papers \[1-4\], giving many explicit enumerations for various classes of planar maps. Fifteen years later, the same problem became popular among physicists in the context of the perturbative expansion of SU($`N`$) gauge field theory . Indeed, at large $`N`$, the dominant Feynman diagrams correspond precisely to planar maps. Explicit enumeration formulas extending previous results were derived by matrix integral techniques , a tool which proved very powerful for such problems. Even more recently, planar maps were used in physics as tessellations of random surfaces, in the context of both discretized 2D quantum gravity and fluid membranes (see e.g.  and ).

Planar maps are formally defined as proper embeddings of graphs into the two-dimensional Riemann sphere, considered up to continuous deformations. A map is characterized by a number of vertices, edges and faces. The general question we address here is the enumeration all such maps with prescribed vertex degrees, i.e, for each $`k`$, a fixed number $`n_k`$ of $`k`$-valent vertices. Turning to generating functions, we can instead consider the partition function of all maps with a weight $`g_k`$ per $`k`$-valent vertex. It proves also useful to introduce an additional weight $`t`$ per edge so that a term of given degree in $`t`$ accounts for a finite number of maps.

The effective counting can be performed in a mechanical way by solving the generic one-matrix model with standard techniques. This approach provides directly the net result, but its powerful nature may seem magical as it hides all the combinatorial aspects of the problem, and some of its calculations must be understood in a purely formal way. Recently, a new purely combinatorial method for studying planar maps was developed, using so-called conjugacy classes of trees, allowing for a new derivation of Tutte’s results as well as some generalizations (see  for a review). In comparison to our present goal, this approach was restricted to the enumeration of Eulerian planar maps, i.e. maps with only vertices of even valences.

In this paper, we give a combinatorial solution of the general problem of enumeration of planar maps with arbitrary valences, generalizing the method of conjugated trees and elucidating the combinatorial structure of the one-matrix model solution. The paper is organized as follows. In Section 2, we introduce the quantities of interest and recall their expression as provided by the one-matrix model (corresponding derivations are detailed in appendix A). In Section 3, we present our combinatorial proof. We first define appropriate decorated trees whose enumeration makes the connection with the one-matrix model solution. We then establish bijections between classes of such trees and planar maps. Particular cases are discussed in Section 4 and conclusions are gathered in Section 5.

2. Solution via matrix integral

2.1. Definitions

Let us now come to the precise definition of the various generating functions that we wish to compute.

Planar maps may display a large number of internal symmetries making the counting quite subtle. Indeed, in both matrix integral and combinatorial approaches, maps are naturally counted with an inverse symmetry weight. A convenient way to avoid this problem is to consider instead rooted maps, with a marked oriented edge which lifts up ambiguities due to symmetries. We will denote by $`E(t;\{g_k\}_{k1})`$ the generating function for rooted planar maps.

Beside these “closed diagrams”, it will prove convenient to also consider one-leg diagrams, namely planar maps with a distinguished univalent vertex (the endpoint of the external leg). For these diagrams, no symmetry problem arises and we denote by $`\mathrm{\Gamma }_1(t;\{g_k\})`$ their generating function with no weight for the distinguished vertex.

Finally, we will also consider the generating function $`\mathrm{\Gamma }_2(t;\{g_k\})`$ for two-leg diagrams with two unweighted distinguished (as, say, incoming and outcoming) univalent vertices adjacent to the same face.

Fig. 1: Typical planar maps contributing respectively to the generating functions $`E`$, $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$, namely (a) a rooted planar map with one regular univalent, one bivalent, three trivalent and two tetravalent vertices; (b) a one-leg diagram with one regular univalent, one bivalent, one trivalent, three tetravalent and one pentavalent vertices; (c) a two-leg diagram with one regular univalent, five trivalent and two tetravalent vertices. In the first case, the rooted edge is marked by an arrow. In the two latter, the distinguished univalent vertices attached to the external legs are marked by a cross. In (c), the two legs must lie in the same face and are distinguished as in- and out- coming via an orientation.

Examples of maps contributing to $`E`$, $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ are displayed in Fig. 1. Our convention to represent maps in the plane is to choose as the external face that on the right of the marked edge for closed diagrams, and that adjacent to the distinguished vertices for both one- and two-leg diagrams.

When cutting the marked edge in a rooted closed map, we get either two disconnected one-leg diagrams, or a two-leg diagram with legs in the same face, excluding that formed by a single oriented edge. This leads to the first general relation

$$E=\frac{\mathrm{\Gamma }_1^2+\mathrm{\Gamma }_2t}{t}$$

2.2. One-matrix model results

Formally, all quantities of interest can be obtained as the large $`N`$ limit of average values in a $`N\times N`$ one-matrix model with partition function

$$Z=𝑑M\mathrm{exp}N\mathrm{Tr}\left(\frac{M^2}{2t}+\underset{k1}{}\frac{g_k}{k}M^k\right)$$

where $`dM`$ denotes the standard Haar (SU($`N`$) invariant) measure over hermitian matrices (see e.g.  and  for details).

As shown in appendix A, the solution of the problem involves two functions $`S(t;\{g_k\})`$ and $`R(t;\{g_k\})`$ characterizing the eigenvalue distribution at large $`N`$. Introducing the potential

$$V(x)=\underset{k1}{}\frac{g_k}{k}x^k$$

and the formal parametrization

$$Q(z)=z+S+R/z$$

in terms of a dummy variable $`z`$, the functions $`S`$ and $`R`$ are implicitly determined by

$$\begin{array}{cc}\hfill S& =tV^{}(Q)|_{z^0}\hfill \\ \hfill R& =t+tV^{}(Q)|_{z^1}\hfill \end{array}$$

where $`V^{}(Q)|_{z^m}`$ denotes the coefficient of $`z^m`$ in $`V^{}(Q)`$ when viewed as a Laurent series in $`z`$. The correct determination is fixed by $`S(0;\{g_k\})=R(0;\{g_k\})=0`$, allowing for an expansion as a power series in $`t`$.

The generating functions for maps are then given in terms of $`S`$ and $`R`$ through

$$\begin{array}{cc}\hfill \mathrm{\Gamma }_1& =SV^{}(Q)|_{z^2}\hfill \\ \hfill \mathrm{\Gamma }_2& =RV^{}(Q)|_{z^3}\left(V^{}(Q)|_{z^2}\right)^2\hfill \end{array}$$

and $`E`$ follows from Eq. (2.1).

In principle, the above equations solve the enumeration problem. In practice, the generic term of the series can be computed explicitly in simple cases such as cubic maps ($`g_k=0`$ when $`k3`$), or Eulerian maps ($`g_k=0`$ for odd $`k`$) as discussed in Section 4 below. As an illustration, we list the first few terms in Eq. (2.1) by expanding the potential (2.1) up to the quartic term

$$\begin{array}{cc}\hfill S& =tg_1+tg_2S+tg_3(S^2+2R)+tg_4(S^3+6RS)+\mathrm{}\hfill \\ \hfill R& =t+tg_2R+tg_3(2RS)+tg_4(3S^2R+3R^2)+\mathrm{}\hfill \end{array}$$

which can be solved order by order in $`t`$ with the initial conditions $`R=t+O(t^2)`$, $`S=tg_1+O(t^2)`$.

In the next Section, we will recover precisely the same set of equations through a purely combinatorial approach.

3. Combinatorial Solution

We now come to the core of the paper, namely the combinatorial rederivation of the above equations (2.1) and (2.1). We first show that the functions $`S`$ and $`R`$ are indeed generating functions of decorated “blossom trees” generalizing those introduced for the counting of Eulerian maps in Ref. . These trees are to be closed into graphs to recover planar maps. As in Ref. , bijections can be established between planar maps and properly defined conjugacy classes of these trees, which we will describe precisely. We will show the equivalence to planar maps both in the case of one-leg and two-leg diagrams, allowing to recover all the results of the matrix model described above.

3.1. Rooted blossom trees

Fig. 2: A typical S-tree (S) and a typical R-tree (R). The leaves are represented with empty arrows, buds correspond to filled black arrows, while regular vertices are solid black dots. The leaves carry a charge $`+1`$ while the buds carry a charge $`1`$. The S-tree has a total charge $`0`$ ($`7`$ leaves and $`7`$ buds), while the R-tree has a total charge $`1`$ ($`5`$ leaves and $`4`$ buds). Both have the property that any descendent subtree not reduced to a bud has total charge $`0`$ or $`1`$.

Let us now define rooted S- and R- blossom trees. These are (finite) rooted planar trees with three species of endpoints that we will refer to as leaves, buds, and regular univalent vertices. We assign a “charge” $`q=+1`$ (resp. $`q=1`$) to leaves (resp. buds), while regular univalent vertices remain neutral ($`q=0`$). These trees are called S- (resp. R-) trees iff: (i) their total charge is $`0`$ (resp. $`1`$) and (ii) any descendent subtree not made of a single bud has total charge $`0`$ or $`1`$.

As an example, the smallest S-tree is made of a single regular univalent vertex attached to the root while the smallest R-tree is made of a single leaf attached to the root. Typical S- and R- trees are represented in Fig. 2.

Fig. 3: Illustration of the recursive generation of a rooted S-tree via the enumeration of all possible vertices of total charge $`0`$ attached to the root for terms up to $`g_4`$ in the potential $`V`$. This parallels the first line of the recursion (2.1) by viewing the labels $`S`$ and $`R`$ as generating functions for S- and R-trees.

Fig. 4: The same illustration as in Fig. 3 for a rooted R-tree (with now vertices of total charge $`1`$), to be paralleled with the second line of the recursion (2.1) for terms up to $`g_4`$ in the potential $`V`$.

Obviously, any descendent subtree of an S- or R-tree not reduced to a bud is itself an S- or an R-tree. This recursive property allows to interpret the functions $`S`$ (resp. $`R`$) of Eq. (2.1) as the generating functions for rooted S-trees (resp. R-trees) with a weight $`g_k`$ per $`k`$-valent vertex ($`k1`$) and a weight $`t`$ per edge not leading to a bud. First, we note that the dummy variable $`z`$ in Eq. (2.1) may be thought of as a fugacity per unit of charge. Viewing the three terms in $`Q(z)=z+S+R/z`$ as associated respectively to a bud ($`q=1`$), an S-tree ($`q=0`$) and an R-tree ($`q=+1`$), the power of $`1/z`$ measures the total charge of any composite object. This allows to interpret $`V^{}(Q)|_{z^m}`$ as generating all possible vertices with descendents being buds, S-trees and R-trees, with a total charge $`q=m`$ and with weight $`g_k`$ per $`k`$-valent vertex (i.e. with $`(k1)`$ descendents). Expressions (2.1) for S-trees and R-trees follow by enumeration of all possible configurations around the vertex attached to the root, with a total charge of $`0`$ and $`1`$ respectively. The factor $`t`$ in front of $`V^{}(Q)`$ accounts for the weight of the root edge, while the additional $`t`$ term in $`R`$ stands for the germ of the recursion. These recursions are depicted in Figs. 3 and 4, for the first few terms involving $`1`$, $`2`$, $`3`$ and $`4`$-valent vertices.

3.2. Conjugated trees

Starting with the above S- and R-trees, we may now wipe out the root by replacing it by a standard (unmarked) leaf, thus increasing the total charge by one unit. This defines unrooted S- and R-trees with total charge $`1`$ and $`2`$ respectively. Two rooted S- (resp. R-) trees are said to belong to the same conjugacy class iff they lead to the same unrooted tree.

The case of S-trees is simple in that reciprocally, given any unrooted S-tree, we obtain all the corresponding rooted S-trees by picking any of its leaves and choosing it as the root. This leads to a $`n`$-to-one correspondence between rooted and unrooted S-trees with $`n`$ leaves (including the former root). The proof goes as follows. Recall first that cutting any edge not leading to a bud in a rooted S-tree cuts out a descendent subtree of charge $`0`$ or $`1`$. Consequently, cutting the same edge in the corresponding unrooted S-tree separates the tree into a piece of charge $`0`$ or $`1`$ and a complementary piece of charge $`1`$ or $`0`$ respectively, as the total charge is $`1`$. This property, valid for any edge not leading to a bud, completely characterizes the unrooted S-trees. Indeed, starting from an unrooted tree with this property, we see that it has total charge $`1`$, and that therefore by replacing an arbitrary leaf by a root, we end up with a rooted tree of total charge $`0`$ whose descendent subtrees not reduced to a bud have charge $`0`$ or $`1`$, the characterization of rooted S-trees. As the choice of leaf taken as a root is arbitrary among the $`n`$ leaves, this moreover shows that the conjugacy class is made of $`n`$ rooted S-trees, by finally noting that the unrooted S-trees have no accidental rotational symmetry since the numbers of buds ($`n1`$) and leaves ($`n`$) are coprime. This completes the abovementioned $`n`$-to-one correspondence. A last remark is in order: cutting any edge not leading to a bud of an unrooted S-tree separates the tree into two rooted trees, one being a rooted S-tree and the other a rooted R-tree. Indeed these trees have respective charges $`0`$ and $`1`$ and any of their descendent subtrees not made of a single bud also have charge $`0`$ or $`1`$ as pieces of the unrooted S-tree.

Fig. 5: The core of a typical unrooted R-tree. We have indicated by double-lines the edges “of type $`0`$-$`2`$” to be cut in the core construction process. The circled pieces are the maximal S-subtrees (of charge $`0`$) attached to the core (of charge $`2`$), represented in thick lines. In this particular example, the core has three leaves and one bud.

The case of R-trees is more subtle. Indeed, it is no longer true that picking as a root any leaf of an unrooted R-tree leads to a rooted R-tree. In order to characterize the subset of admissible leaves (leading effectively to rooted R-trees), we first define the core of an unrooted R-tree by the following procedure: (i) we mark all the inner edges which separate the tree into two pieces of respective charge $`0`$ and $`2`$, (ii) we remove all these edges and (iii) we define the core as the only connected component with total charge $`2`$. This is illustrated in Fig. 5. To show the existence and uniqueness of the core, we consider an arbitrary admissible leaf and the corresponding rooted R-tree. Due to charge constraints, an edge is marked at step (i) iff the descendent subtree originating from it is an S-tree. Such a subtree is then amputated from its own S-subtrees, splitting eventually into connected components all of charge $`0`$. The connected component containing the root is the only one not obtained in this way, hence it has charge $`2`$. Note also that all trees attached to the core are rooted S-trees attached by their root. We call them maximal S-subtrees. As our procedure was defined independently of the admissible root at hand, it follows that all admissible leaves belong to the core. Conversely, replacing an arbitrary leaf in the core by a root, we get a rooted tree such that any descendent subtree not made of a single bud either is contained in an S-subtree, or originates from an edge in the core not marked at step (i): in both cases its total charge is $`0`$ or $`1`$ as wanted. To conclude, the conjugacy class associated with an unrooted R-tree is made of all trees rooted at leaves of the core. As opposed to the S-case, care must be taken with a possible two-fold rotational symmetry (half turn) which is the only one compatible with the fact that there are two more leaves than buds in an unrooted R-tree.

Fig. 6: The core of a typical rooted R-tree as obtained from the semi-developed representation using only the vertices of Fig. 4. Replacing the “leaves” marked with a label $`S`$ by arbitrary rooted S-trees, we generate all rooted R-trees sharing the same rooted core.

To conclude this section, let us note that a simple way to generate the core of R-trees is to use a semi-developed representation of rooted R-trees, namely by only using the vertices of Fig. 4 (second line of Eq. (2.1)) and keeping S as a label, as illustrated in Fig. 6.

3.3. Enumeration of one-leg diagrams:

Let us now show that one-leg diagrams are in one-to-one correspondence with unrooted S-trees.

Fig. 7: A sample unrooted S-tree is closed into a one-leg diagram in a unique way, by iteratively matching each bud to the closest available leaf in counterclockwise direction. The resulting edges are indicated in thin solid lines, and form a system of non-intersecting arches around to the tree. We have also indicated for each bud and leaf the corresponding depth, namely that of the arch in the arch system, starting with depth $`1`$ for arches adjacent to the external face. The unmatched vertex receives the depth $`0`$ and serves as the distinguished endpoint of the one-leg diagram.

Starting from an unrooted S-tree, we build a map by matching buds and leaves as in Ref. , by connecting iteratively each bud to the closest available leaf in counterclockwise direction as shown in Fig. 7, in such a way that the resulting graph is planar (i.e. with no intersection of edges). Such connected bud-leaf pairs are replaced by regular edges. As the total charge of $`1`$ counts the number of leaves minus that of buds, the matching procedure leaves exactly one unmatched leaf which we replace by a distinguished univalent vertex. The net result is a one-leg diagram.

Fig. 8: A sample one-leg diagram is cut into an unrooted S-tree by visiting edges in counterclockwise direction around the diagram and iteratively cutting those which do not disconnect the remaining diagram, until a tree is obtained. The cut edges are indicated by a double-line, they are to be replaced by a bud, followed by a leaf. We have also indicated the depth of the corresponding cut edges, which will translate into the depths of buds and leaves of the S-tree. These depths correspond to the number of visits to the distinguished endpoint before the edge is cut.

Conversely, starting from a one-leg diagram, we recover an unrooted S-tree by an inverse algorithm similar to that of Ref. . Starting from the distinguished endpoint, we successively visit all edges adjacent to the external face in counterclockwise direction as shown in Fig. 8. At each step, the edge is cut iff the cutting does not disconnect the diagram. In the cutting procedure, the first half of the edge is replaced by a bud and the second half by a leaf. After one turn, the external face has been merged with all its adjacent faces. We then iterate the procedure until all faces have been merged. Replacing the distinguished endpoint by a leaf, the resulting connected diagram is an unrooted blossom tree with buds, leaves and regular vertices, and with total charge $`1`$.

Fig. 9: The tree obtained by cutting a one-leg diagram is an unrooted S-tree. Any of its inner edges, $`e`$, separates the tree into two pieces $`T_1`$ and $`T_2`$, say with the distinguished endpoint in $`T_1`$. We have represented the cut bud-leaf pairs formerly connecting $`T_1`$ to $`T_2`$. The cutting procedure has successively cut these pairs in counterclockwise order around the edge $`e`$. In our example, there is one more bud-leaf pair on the left than on the right of $`e`$ (odd number of passings). We could also have the same number of bud-leaf pairs on both sides (even number of passings). Collecting the total charge in $`T_1`$ and $`T_2`$ we see that in the case of the figure $`T_1`$ has charge $`0`$ while $`T_2`$ has charge $`1`$. In the case of an even number of passings, $`T_1`$ would have charge $`1`$ and $`T_2`$ charge $`0`$. This is the charge characterization of an unrooted S-tree.

Let us now show that this tree actually is an S-tree by considering one of its inner edges (connecting regular vertices), which we denote by $`e`$ (see Fig. 9). This edge separates the tree into two pieces, say $`T_1`$ and $`T_2`$, which contain all the vertices of the original diagram. We choose for $`T_1`$ the piece containing the distinguished endpoint. In the cutting procedure, $`e`$ is not visited until all the other edges of the original diagram connecting a vertex of $`T_1`$ to one of $`T_2`$ have been visited once and cut into a bud-leaf pair. More precisely, these edges are replaced by a bud in $`T_1`$ and a leaf in $`T_2`$ when passing from $`T_1`$ to $`T_2`$ and vice versa. As the passings from $`T_1`$ to $`T_2`$ and from $`T_2`$ to $`T_1`$ alternate, these buds and leaves contribute by a net charge of $`0`$ in both $`T_1`$ and $`T_2`$ if the number of passings is even, while they contribute by a net charge of $`1`$ in $`T_1`$ and $`+1`$ in $`T_2`$ if the number is odd. Adding the charge $`+1`$ of the distinguished vertex changed into a leaf and noting that all other edges not connecting $`T_1`$ to $`T_2`$ have neutral contributions, we end up with two pieces of respective charges $`0`$ and $`1`$, the characterization of unrooted S-trees.

To further establish the bijection, it remains to prove that the two algorithms described above are inverse of each other. This is best seen by introducing the notion of depth of buds and leaves within an unrooted S-tree. Starting from a one-leg diagram and applying the above cutting algorithm, we may associate to each cut edge a natural “depth” $`d=1,2,3,\mathrm{}`$ by the number of visits to the distinguished vertex before this edge is cut. The same depth is attached to the corresponding bud and leaf, and finally the extra leaf replacing the distinguished vertex is given the depth $`0`$. We have indicated the depths of buds and leaves in Figs. 7 and 8 for illustration. Remarkably, the depth is a notion intrinsic to the resulting unrooted S-tree. Indeed, in the bud-leaf matching procedure of an unrooted S-tree, the created edges form a system of arches around the tree and the depth of the buds and leaves is nothing but that of the corresponding arches, starting with depth $`1`$ for external arches, and moreover associating the depth $`0`$ to the unmatched leaf. This system of arches may also be constructed now as leaf-bud pairs of decreasing depth as follows. Starting from the depth $`0`$ leaf we proceed clockwise around the tree and connect the first encountered leaf-bud pair of maximal depth, say $`k`$, and continue to connect leaves and buds of depth $`k\mathrm{}+1`$ after $`\mathrm{}`$ visits to the depth $`0`$ leaf. This alternative closing procedure is the exact inverse of the cutting algorithm for one-leg diagrams.

As a direct consequence of the above bijection, the enumeration of one-leg diagrams is equivalent to that of unrooted S-trees. The latter is very simply performed by noting that there is one more leaf than bud in such a tree. The corresponding generating function is therefore equal to the difference between that of unrooted S-trees with a marked leaf and that of unrooted S-trees with a marked bud. The generating function of unrooted S-trees with a marked leaf is nothing but $`S`$ as each of these leaves may be taken as a root. On the other hand, the generating function for unrooted S-trees with a marked bud is computed by noting that these trees are in one-to-one correspondence with rooted trees of total charge $`2`$ whose descendent subtrees are buds, S- or R-trees, as is easily seen by replacing the marked bud by a neutral root. The generating function for unrooted S-trees with a marked bud therefore reads $`V^{}(Q)|_{z^2}`$, leading to $`\mathrm{\Gamma }_1=SV^{}(Q)|_{z^2}`$ as in (2.1). QED.

3.4. Enumeration of two-leg diagrams:

Fig. 10: A typical two-leg diagram (a) is cut into an unrooted R-tree (b) by first applying the cutting algorithm (cut edges are indicated by parallel thin lines in (a)) and then replacing the in- and out- coming endpoints by leaves. We have indicated the depths of the leaves and buds in the R-tree (b), as well as maximal S-subtrees (circled pieces) and the core to which they are attached (thickened edges). Note that the outgoing (depth $`0`$) leaf is in the core, while the incoming one is not, as it lies in a maximal S-subtree. This gives an example of an unrooted R-tree in the set R1, i.e. with exactly one depth $`0`$ leaf in the core.

Let us now look for a similar equivalence in the case of two-leg diagrams. We will show that two-leg diagrams are in one-to-one correspondence with unrooted R-trees with a marked leaf of depth $`0`$ in the core.

We apply the same cutting procedure as before to two-leg diagrams (see Fig. 10 (a)) starting from the incoming distinguished endpoint and finally replacing the two distinguished endpoints by leaves. The resulting tree (see Fig. 10 (b)) has total charge $`2`$. Let us now show that: (i) it is an unrooted R-tree and (ii) the originally outcoming distinguished endpoint is a leaf belonging to the core. This is easily seen by noting that in the cutting procedure, the outcoming endpoint is treated as a regular univalent vertex. Therefore, viewing the original two-leg diagram as a one-leg diagram with an extra marked regular univalent vertex, we may apply the result of Section 3.3 showing that the resulting tree is an unrooted S-tree with this same marked regular univalent vertex. Cutting the edge leading to this vertex leaves us with two pieces, one rooted S-tree and one rooted R-tree. Obviously, the S-tree is the piece made of the univalent vertex attached to a root. Therefore, the other piece is a rooted R-tree. This proves that picking the outcoming endpoint as a root leads to an R-tree showing both (i) and (ii) as the admissible roots are the leaves of the core.

Up to this point, we have never used the fact that the two endpoints of the diagram both lie in the external face. In fact, the above construction shows that $`R`$ can be understood as the generating function $`\mathrm{\Gamma }_{1,1}`$ of two-leg diagrams whose legs do not have to lie in the same face. This is also proved in appendix A in the matrix language. The requirement that the two endpoints belong to the external face implies a further restriction on the unrooted R-tree obtained through the cutting procedure, namely that the leaf replacing the outcoming endpoint has depth $`0`$, as illustrated in Fig. 10 (b). As before, we indeed have an intrinsic notion of depth for unrooted R-trees, defined through the same bud-leaf counterclockwise matching algorithm. For unrooted R-trees, it now leads to two arch systems separated by two unmatched leaves. These two leaves are assigned a depth $`0`$ while all the other buds or leaves are assigned the depth of the corresponding arch. For two-leg diagrams, the external legs naturally separate the edges to be cut into two independent arch system, showing that, when replaced by leaves, the incoming and outcoming endpoints both have depth $`0`$ (see Fig. 10 (b)). Finally, the inverse of the above cutting procedure is easily identified as before with the suitable alternative matching algorithm connecting leaves and buds clockwise in decreasing depth order.

The enumeration of two-leg diagrams is finally reduced to that of unrooted R-trees with a marked leaf of depth $`0`$ in the core. As a preliminary remark, note that there are exactly two leaves of depth $`0`$ in an unrooted R-tree, which may or may not belong to the core. This suggests to classify the unrooted R-trees into three subsets R0, R1 and R2 according to the number $`0`$, $`1`$ or $`2`$ of depth $`0`$ leaves in the core. For instance, the unrooted R-tree of Fig. 10 (b) belongs to R1, as the incoming (depth $`0`$) leaf lies in a maximal S-subtree (circled in the figure), while the outcoming one is in the core (represented by thickened edges in the figure). Denoting by $`R_0`$, $`R_1`$ and $`R_2`$ the corresponding generating functions, we have

$$\mathrm{\Gamma }_2=R_1+2R_2=2\stackrel{~}{R}(R_1+2R_0)$$

where $`\stackrel{~}{R}R_0+R_1+R_2`$ is the generating function for unrooted R-trees.

$`\stackrel{~}{R}`$ can be computed as before by noting that there are exactly two more leaves than buds in the core. Hence $`2\stackrel{~}{R}`$ is the difference between the generating function for unrooted R-trees with a marked leaf in the core, and that for unrooted R-trees with a marked bud in the core. The former is nothing but $`R`$ according to Section 3.2. On the other hand, the generating function for unrooted R-trees with a marked bud in the core is computed by noting that these trees are in one-to-one correspondence with rooted trees of total charge $`3`$ whose descendent subtrees are buds, S- or R-trees. This is seen again by replacing the marked bud by a neutral root and checking that all proper descendent subtrees not reduced to a bud have charge $`0`$ or $`1`$. Indeed such a descendent subtree either is itself a descendent of an S-subtree, hence has charge $`0`$ or $`1`$, or originates from an edge of the core, hence has charge $`1`$. The generating function for unrooted R-trees with a marked bud in the core therefore reads $`V^{}(Q)|_{z^3}`$, leading to $`2\stackrel{~}{R}=RV^{}(Q)|_{z^3}`$. Remark that the accidental two-fold symmetry of unrooted R-trees is properly accounted for in $`\stackrel{~}{R}`$ in which symmetric contributions (which contribute only once to $`R`$ and $`V^{}(Q)|_{z^3}`$) receive a weight $`1/2`$.

The remaining term in Eq. (3.1) can be computed with the result $`(R_1+2R_0)=(V^{}(Q)|_{z^2})^2`$. The proof is slightly tedious and is detailed in appendix B.

This finally allows to express Eq. (3.1) in the form of Eq. (2.1). QED.

4. Discussion

4.1. Cubic maps

The first non-trivial case of planar maps is that of cubic maps, namely maps whose regular vertices are all trivalent (closed planar cubic maps are dual to triangulations of the sphere). These correspond in our notations to taking $`g_k=g\delta _{k,3}`$, in which case $`V^{}(Q)=gQ^2`$, and the fundamental relations (2.1) read

$$\begin{array}{cc}\hfill S& =gt(S^2+2R)\hfill \\ \hfill R& =t+2gtRS\hfill \end{array}$$

while the one- and two-leg diagram generating functions read

$$\begin{array}{cc}\hfill \mathrm{\Gamma }_1& =SgR^2\hfill \\ \hfill \mathrm{\Gamma }_2& =Rg^2R^4\hfill \end{array}$$

Note that S- and R-trees are now binary trees, as they originate from cubic maps. In this particular case, we also have the following additional relation

$$\mathrm{\Gamma }_1=gt(\mathrm{\Gamma }_2+\mathrm{\Gamma }_1^2)$$

obtained by cutting out from any one-leg diagram the leg and the trivalent vertex connected to it: the remaining part is made either of a connected two-leg diagram (term in $`\mathrm{\Gamma }_2`$) or of two disconnected one-leg diagrams (contributing to a term in $`\mathrm{\Gamma }_1^2`$). Alternatively, the relation (4.1) is also easily proved by use of Eqs. (4.1) and (4.1).

Fig. 11: In the cubic map case, the core of a typical unrooted R-tree is a chain with a leaf at each end and attached maximal S-subtrees on either side.

Another drastic simplification is that the core of any unrooted R-tree is a simple chain joining two leaves, with no buds and with attached maximal S-trees (see Fig. 11). The absence of buds is a consequence of $`V^{}(Q)|_{z^3}=0`$ in this case, leading moreover to $`2\stackrel{~}{R}=R`$, as there are two admissible leaves that can serve as roots.

The actual numbers of S- and R-trees as well as one- and two-leg diagrams with fixed number of edges are obtained by first eliminating $`R`$ from Eq. (4.1), and then applying the Lagrange inversion formula. Introducing the rescaled functions and variables $`\sigma =gtS`$, $`\rho =(gt)^2R`$, $`\theta =g^2t^3`$, we get

$$\begin{array}{cc}\hfill 2\theta & =\sigma (1\sigma )(12\sigma )\hfill \\ \hfill \rho & =\frac{\sigma (1\sigma )}{2}\hfill \end{array}$$

Denoting by $`\phi (\sigma )=\sigma (1\sigma )(12\sigma )`$, we may now express the series expansion in $`\theta `$ of any function $`h(\sigma )`$ (with $`h(0)=0`$) through the Lagrange inversion formula, a direct consequence of the Cauchy formula:

$$\begin{array}{cc}\hfill h(\sigma )& =\frac{ds}{2i\pi }\frac{h(s)\phi ^{}(s)}{\phi (s)\phi (\sigma )}\hfill \\ & =\underset{n=1}{\overset{\mathrm{}}{}}\frac{(2\theta )^n}{n}\frac{ds}{2i\pi }\frac{h^{}(s)}{\phi (s)^n}\hfill \end{array}$$

where we have used the first line of (4.1) and integrated by parts. Upon picking respectively $`h(x)=x`$ and $`h(x)=x(1x)/2`$, we get the series coefficients $`\sigma _n`$ and $`\rho _n`$ for $`S=(gt)^1_{n=1}^{\mathrm{}}\sigma _n\theta ^n`$ and $`R=(gt)^2_{n=1}^{\mathrm{}}\rho _n\theta ^n`$ as contour integrals, easily calculated as $`\sigma _1=2`$, $`\rho _1=1`$ and for $`n2`$:

$$\begin{array}{cc}\hfill \sigma _n& =\frac{2^n}{n}\frac{ds}{2i\pi s^n}\frac{1}{(1s)^n(12s)^n}\hfill \\ & =\frac{2^n}{n}\underset{k=0}{\overset{n}{}}2^k\left(\genfrac{}{}{0pt}{}{2nk2}{n1}\right)\left(\genfrac{}{}{0pt}{}{n+k1}{n1}\right)\hfill \\ & =\frac{2^{2n1}}{n!}(n+1)(n+3)\mathrm{}(3n5)(3n3)\hfill \\ \hfill \rho _n& =\frac{2^{n1}}{n}\frac{ds}{2i\pi s^n}\frac{1}{(1s)^n(12s)^{n1}}\hfill \\ & =\frac{2^{n1}}{n}\underset{k=0}{\overset{n}{}}2^k\left(\genfrac{}{}{0pt}{}{2nk2}{n1}\right)\left(\genfrac{}{}{0pt}{}{n+k2}{n2}\right)\hfill \\ & =\frac{2^{2n2}}{n!}n(n+2)(n+4)\mathrm{}(3n6)(3n4)\hfill \end{array}$$

for the numbers $`\sigma _n,\rho _n`$ of rooted S- and R-trees with $`n`$ leaves.

Similarly, we get expressions for the numbers of one- (resp. two-) leg diagrams with $`3n1`$ (resp. $`3n2`$ edges), denoted by $`\gamma _n^{(1)}`$ (resp. $`\gamma _n^{(2)}`$) by expressing the relevant generating functions in terms of $`\sigma `$ only, with the result $`\gamma _1^{(1)}=1`$, $`\gamma _1^{(2)}=1`$ and for $`n2`$:

$$\begin{array}{cc}\hfill \gamma _n^{(1)}& =\frac{2^{2n1}}{(n+1)!}(n+1)(n+3)\mathrm{}(3n5)(3n3)\hfill \\ \hfill \gamma _n^{(2)}& =\frac{3}{n+2}\frac{2^{2n2}}{n!}n(n+2)(n+4)\mathrm{}(3n6)(3n4)\hfill \end{array}$$

As a direct consequence of the S-tree conjugacy, we have $`\gamma _n^{(1)}=\sigma _n/(n+1)`$, as a corresponding unrooted S-tree has exactly $`n+1`$ leaves which all may serve as a root in the conjugacy class. We note a slightly less trivial relation between the numbers of two-leg diagrams and R-trees: $`\gamma _n^{(2)}=3\rho _n/(n+2)`$ which still awaits a good combinatorial interpretation.

Finally, we get the number $`e_n`$ of rooted cubic maps with $`3n`$ edges by substituting Eq. (4.1) into Eq. (2.1), leading to $`E=\mathrm{\Gamma }_1/(gt^2)1`$, and therefore

$$e_n=\gamma _{n+1}^{(1)}=\frac{2^{2n+1}}{(n+2)!}(n+2)(n+4)(n+6)\mathrm{}(3n2)(3n)$$

for all $`n1`$.

Fig. 12: The local environment in an unrooted S-tree in the cubic case around (a) a leaf, (b) a bud, (c) a vertex with no bud, and (d) an edge not leading to a bud.

Introducing the rescaled one-leg diagram generating function $`G_1=gt\mathrm{\Gamma }_1`$, we easily get the following differential equations

$$\begin{array}{cc}\hfill \left(\theta \frac{d}{d\theta }+1\right)G_1& =\sigma \hfill \\ \hfill \theta \frac{d}{d\theta }G_1& =\frac{\rho ^2}{\theta }\hfill \\ \hfill \left(\theta \frac{d}{d\theta }1\right)G_1& =\frac{\rho \sigma ^2}{\theta }\hfill \\ \hfill \left(3\theta \frac{d}{d\theta }1\right)G_1& =\frac{\rho \sigma }{\theta }\hfill \end{array}$$

obtained by marking in an unrooted S-tree respectively a leaf, a bud, a trivalent vertex with no bud and an edge, in respective numbers $`n+1`$, $`n`$, $`n1`$, $`3n1`$. The combinatorial origin of these relations is depicted in Fig. 12.

4.2. Eulerian maps

The case of Eulerian maps, namely of diagrams with only regular vertices of even valences was solved by Tutte , and takes a much simpler form than the general case. Eulerian maps may be recovered here by demanding that $`V`$ be an even function, namely that all $`g_{2k1}=0`$, $`k=1,2,\mathrm{}`$ The matrix model solution (see appendix A) shows in that case that $`S=0`$ as a consequence of the symmetry $`xx`$ of the planar eigenvalue density. From the present combinatorial point of view, it is easy to see that the recursion relations of Figs. 3 and 4 lead to infinite rooted trees unless we decide that there are no rooted S-trees at all. In other words, our finiteness requirement forces us to take the trivial generating function $`S=0`$.

With this simplification, we are left with only rooted R-trees, which become the fundamental combinatorial objects, now characterized by the property that (i) their total charge is $`1`$ and (ii) all their descendent subtrees not reduced to a bud have charge $`1`$ as well (thus are themselves R-trees). These trees are those which have been considered in Ref.  in the original connection between Eulerian maps and conjugated trees. Moreover, in the absence of S-trees, it is easy to see that the core of an unrooted R-tree is actually the whole tree itself, with all the leaves as admissible roots. The notion of conjugacy therefore becomes much simpler and coincides with that used in Ref. .

5. Conclusion

In this paper, we have elucidated the combinatorial structure of the one-matrix model solution to the problem of enumerating planar maps with prescribed vertex degrees. Our construction generalizes that of Ref.  for Eulerian maps using two types of trees instead of one and a refined notion of conjugacy of trees. The need for these two types of trees is inherent to the absence of (twofold) symmetry in arbitrary planar graphs as opposed to Eulerian ones (face bicolorability) which can be traced back to the absence of $`\text{ZZ}_2`$ symmetry in the matrix model.

Here we have considered arbitrary graphs with no restrictions on reducibility. Irreducible maps were enumerated in some particular cases by Tutte    and a connection with conjugated trees was obtained recently in Ref. . In the matrix integral language, the construction of ($`p`$-particle) irreducible diagrams can be performed as already shown in  through appropriate renormalizations. It does not involve the introduction of new fundamental objects and we therefore believe that the corresponding enumeration can be understood also in terms of S- and R- trees.

Finally, we have many other solvable matrix models at hand and we may expect that, for some of them, the underlying combinatorics can be understood directly in terms of appropriate conjugated trees. This is already the case for particular two-matrix models, describing the enumeration of bipartite regular graphs, recently reformulated in terms of conjugated trees \[12-15\]. The hope is that it may also apply for instance to the case of interacting two matrix models such as that describing the two-dimensional Ising model  or hard particle model  on random lattices. This also suggests to look for a more general picture relating the tree structure to the (Toda-like) integrable structure of matrix models.

Appendix A. Planar maps from the one matrix model

In this appendix, we recall the derivation of the planar limit of the integral (2.1). For the sake of technical simplicity, we must take for the potential $`V(x)`$ a truncated polynomial form namely

$$V(x)=\underset{k=1}{\overset{I}{}}g_k\frac{x^k}{k}$$

The result is then trivially extended so as to involve arbitrarily many $`g_k`$’s.

The planar limit of the integral (2.1) is evaluated in a succession of standard steps as follows (see e.g. Refs.  and  for details). The first step consists in reducing the integral (2.1) to one over the eigenvalues of $`M`$. This is readily done by changing variables $`M(m,U)`$, where $`M=UmU^{}`$, and $`m`$ is a real diagonal matrix, while $`UU(N)/U(1)^N`$ is a unitary matrix defined up to the right multiplication by an arbitrary diagonal matrix with entries of the form $`e^{i\varphi _j}`$ with real phases $`\varphi _j`$, $`j=1,2,\mathrm{},N`$. The Jacobian of this change of variables reads $`J=\mathrm{\Delta }(m)^2`$, where $`\mathrm{\Delta }(m)=det(m_i^{j1})|_{1i,jN}`$ is the Vandermonde determinant of $`m`$. As the integrand of (2.1) only depends on $`m`$, the angular variables $`U`$ may be integrated out, with the result

$$Z_N(t,\{g_k\})=𝑑m\mathrm{\Delta }(m)^2e^{N\mathrm{Tr}\left(\frac{m^2}{2t}V(m)\right)}$$

up to an overall unimportant multiplicative constant (the volume of the angular variables’ locus).

The second step consists in evaluating the large $`N`$ limit of (A.1) by putting this integral into a more familiar form, and evaluating it via a saddle-point method. Indeed, introducing the functional

$$S(m)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left(\frac{m_i^2}{2t}V(m_i)\right)\frac{1}{N^2}\underset{1ijN}{}\mathrm{Log}|m_im_j|$$

the integral (A.1) takes the form $`Z_N=𝑑me^{N^2S(m)}`$. In the limit of large $`N`$, this integral is dominated by the solution $`m`$ to the saddle-point equations $`S/m_i=0`$, namely

$$\frac{m_i}{t}V^{}(m_i)=\frac{1}{N}\underset{\genfrac{}{}{0pt}{}{j=1}{ji}}{\overset{N}{}}\frac{1}{m_im_j}$$

In particular, we get the planar (genus zero) free energy of the model as the limit

$$f(t,\{g_k\})\underset{N\mathrm{}}{lim}\frac{1}{N^2}\mathrm{Log}Z_N(t,\{g_k\})=S(m)$$

evaluated at the solution $`m`$ to (A.1). To further compute this solution, let us introduce the resolvent $`\omega _N(\zeta )`$

$$\omega _N(\zeta )=\frac{1}{N}\underset{i=1}{\overset{N}{}}\frac{1}{\zeta m_i}$$

evaluated on the solution $`m`$ to (A.1). This function is related to the eigenvalue density $`\rho _N(\zeta )`$ defined as

$$\rho _N(\zeta )=\frac{1}{N}\underset{i=1}{\overset{N}{}}\delta (\zeta m_i)$$

also evaluated on the solution $`m`$ to (A.1). More precisely we have $`\rho _N(\zeta )=(\omega (\zeta +i0)\omega (\zeta i0))/(2i\pi )`$, so the density of eigenvalues is concentrated on the (real) singularities of $`\omega _N`$. By the definition (A.1), it is easy to see that $`\omega _N(\zeta )1/\zeta `$ when $`\zeta `$ is large, expressing the normalization of the density $`_{\mathrm{IR}}𝑑x\rho _N(x)=1`$. Let us now multiply both sides of (A.1) by $`1/(\zeta m_i)`$ and sum over $`i=1,2,\mathrm{},N`$. We find the following simple quadratic differential equation for $`\omega _N(\zeta )`$

$$\frac{1}{N}\frac{d\omega _N(\zeta )}{d\zeta }=\omega _N(\zeta )^2\left(\frac{\zeta }{t}V^{}(\zeta )\right)\omega _N(\zeta )+P_N(\zeta )$$

where $`P_N(\zeta )`$ is defined by

$$P_N(\zeta )=\frac{1}{t}\frac{1}{N}\underset{i=1}{\overset{N}{}}\frac{V^{}(m_i)V^{}(\zeta )}{m_i\zeta }$$

Note that, as $`V^{}`$ is a polynomial of degree $`I1`$, $`P_N`$ is itself a polynomial of degree $`I2`$. In the large $`N`$ limit, the left hand side of (A.1) is negligible, and we may therefore obtain the large $`N`$ resolvent $`\omega (\zeta )lim_N\mathrm{}\omega _N(\zeta )`$ as the solution of a quadratic equation, in the form

$$\omega (\zeta )=\frac{1}{2}\left(\frac{\zeta }{t}V^{}(\zeta )\pm \sqrt{(\frac{\zeta }{t}V^{}(\zeta ))^24P(\zeta )}\right)$$

where $`P`$ denotes the large $`N`$ limit of $`P_N`$ (it is also a polynomial of degree $`I2`$), and the sign $`\pm `$ is selected by requiring that the large $`\zeta `$ asymptotics of $`\omega (\zeta )`$ read $`\omega (\zeta )1/\zeta `$.

In a third step, the planar resolvent $`\omega (\zeta )`$ is further completely fixed by making the standard “one-cut” hypothesis, namely that $`\omega (\zeta )`$ only has a square root singularity of the form $`\sqrt{(\zeta a)(\zeta b)}`$ where $`a`$ and $`b`$ are real, with say $`a<b`$. This expresses the fact that the planar eigenvalue density $`\rho (\zeta )lim_N\mathrm{}\rho _N(\zeta )`$ is assumed to have a compact support made of a single real interval $`[a,b]`$. Picking say $`g_I`$ to be positive, the sign in (A.1) is fixed to be $`+`$, and the resolvent takes the form

$$\omega (\zeta )=\frac{1}{2}\left(\frac{\zeta }{t}V^{}(\zeta )+G(\zeta )\sqrt{(\zeta a)(\zeta b)}\right)$$

where $`G(\zeta )`$ is a polynomial of degree $`I2`$. Writing that $`\omega (\zeta )1/\zeta `$ at large $`\zeta `$ now determines all coefficients of $`G`$ and both $`a`$ and $`b`$. More precisely, we have

$$G(\zeta )=\left[\frac{V^{}(\zeta )\frac{\zeta }{t}}{\sqrt{(\zeta a)(\zeta b)}}\right]_+$$

where we must expand the bracket as a power series of $`\zeta `$ at $`\mathrm{}`$, and the subscript $`+`$ indicates that we must retain only the polynomial part of this Laurent series. Denoting moreover by $`H(\zeta )=(V^{}(\zeta )\frac{\zeta }{t})/\sqrt{(\zeta a)(\zeta b)}`$ the Laurent series in the bracket, we may express directly

$$\omega (\zeta )=\frac{1}{2}[H(\zeta )]_{}\sqrt{(\zeta a)(\zeta b)}$$

where the subscript $``$ indicates that we only retain negative powers of $`\zeta `$ in the Laurent series. Using the Cauchy formula around the infinity in the complex plane, we may express the coefficients of $`[H(\zeta )]_{}=_{m1}\zeta ^mH(\zeta )|_{\zeta ^m}`$ as

$$H(\zeta )|_{\zeta ^m}=\frac{1}{2i\pi }w^{m1}𝑑w\frac{V^{}(w)\frac{w}{t}}{\sqrt{(wa)(wb)}}$$

This integral may be drastically simplified if we change variables to $`z`$:

$$w=Q(z)z+S+\frac{R}{z}$$

and pick $`R`$ and $`S`$ so that the square root disappears. Indeed,

$$(wa)(wb)=z^2+(2Sab)(z+\frac{R}{z})+(Sa)(Sb)+2R+\frac{R^2}{z^2}$$

is a perfect square iff

$$\begin{array}{cc}\hfill S& =\frac{a+b}{2}\hfill \\ \hfill R& =\left(\frac{ba}{4}\right)^2\hfill \end{array}$$

in which case $`\sqrt{(wa)(wb)}=zR/z`$ around the infinity. The change of variables $`wz`$ then leads to

$$H(\zeta )|_{\zeta ^m}=\frac{1}{2i\pi }\frac{dz}{z}\left(V^{}(Q(z))\frac{Q(z)}{t}\right)Q(z)^{m1}$$

For the first few values of $`m=1,2,3,4`$, this gives

$$\begin{array}{cc}\hfill H(\zeta )|_{\zeta ^1}& =V^{}(Q)|_{z^0}\frac{S}{t}\hfill \\ \hfill H(\zeta )|_{\zeta ^2}& =2\left(V^{}(Q)|_{z^1}\frac{R}{t}\right)+SH(\zeta )|_{\zeta ^1}\hfill \\ \hfill H(\zeta )|_{\zeta ^3}& =2V^{}(Q)|_{z^2}+2SH(\zeta )|_{\zeta ^2}+(2RS^2)H(\zeta )|_{\zeta ^1}\hfill \\ \hfill H(\zeta )|_{\zeta ^4}& =2V^{}(Q)|_{z^3}+3SH(\zeta )|_{\zeta ^3}+3(RS^2)H(\zeta )|_{\zeta ^2}+(S^33RS)H(\zeta )|_{\zeta ^1}\hfill \end{array}$$

In the second, third and fourth lines of (A.1), we have used the symmetry $`zR/z`$ of the change of variables to express that $`f(Q)|_{z^m}=R^mf(Q)|_{z^m}`$ for any Laurent series $`f`$, and all $`m0`$, in order to reexpress the result in terms mostly of $`V^{}(Q)`$. We now simply have to substitute the above expressions for the coefficients of $`H`$ into the formula (A.1) for $`\omega (\zeta )`$, and impose that $`\omega (\zeta )=1/\zeta +O(1/\zeta ^2)`$ at infinity. This gives respectively for the order $`\zeta ^0`$ and $`\zeta ^1`$ terms the two equations

$$\begin{array}{cc}\hfill \frac{S}{t}V^{}(Q)|_{z^0}& =0\hfill \\ \hfill \frac{R}{t}V^{}(Q)|_{z^1}& =1\hfill \end{array}$$

which are nothing but (2.1).

The planar resolvent $`\omega (\zeta )`$ may also be used as the generating function for the quantities $`\theta _k=lim_N\mathrm{}\frac{1}{N}\mathrm{Tr}(M^k)`$ that enumerate the (possibly disconnected) planar maps with one special $`k`$-valent vertex and arbitrary regular $`m`$-valent ones weighted by $`g_m`$, $`m=1,2,3,\mathrm{}`$ and with a weight $`t`$ per edge. Indeed, the definition of $`\omega (\zeta )`$ entails that

$$\omega (\zeta )=\underset{k=0}{\overset{\mathrm{}}{}}\zeta ^{k1}\underset{N\mathrm{}}{lim}\frac{1}{N}\mathrm{Tr}(m^k)=\underset{k=0}{\overset{\mathrm{}}{}}\zeta ^{k1}\theta _k$$

where $`m`$ stands for the solution to the saddle-point equations (A.1). In particular, let us use the direct expression (A.1) to compute the coefficients $`\theta _1`$ and $`\theta _2`$ corresponding to one and two-leg diagrams respectively. Substituting (A.1) into (A.1), we get

$$\begin{array}{cc}\hfill \theta _1& =SV^{}(Q)|_{z^2}\hfill \\ \hfill \theta _2& =R+S^2V^{}(Q)|_{z^3}2SV^{}(Q)|_{z^2}\hfill \end{array}$$

Finally, we note that the one-leg diagrams counted by $`\theta _1`$ are always connected hence $`\theta _1=\mathrm{\Gamma }_1`$, while the two-leg ones may be disconnected, and we must write $`\mathrm{\Gamma }_2=\theta _2\theta _1^2`$. This finally leads to Eqs. (2.1).

A first remark is in order. The quantity $`R`$ may be directly expressed as a connected matrix average, namely

$$R=\theta _{1,1}\theta _1^2$$

where

$$\theta _{1,1}=\underset{N\mathrm{}}{lim}\frac{1}{N^2}\mathrm{Tr}(M)\mathrm{Tr}(M)$$

where the average is computed by attaching two external legs to connected planar maps and those two legs do not have to lie in the same face (as they do in $`\mathrm{\Gamma }_2`$). To prove this, the simplest way is to use the linear term in the potential to generate insertions of Tr$`(M)`$ by differentiating with respect to $`g_1`$. For instance, the connected average $`\mathrm{\Gamma }_{1,1}\theta _{1,1}\theta _1^2`$ is easily obtained by differentiating $`\theta _1`$ with respect to $`g_1`$:

$$\begin{array}{cc}\hfill \mathrm{\Gamma }_{1,1}& =\frac{\theta _1}{g_1}\hfill \\ & =\frac{S}{g_1}\frac{V^{}(Q)|_{z^2}}{g_1}\hfill \\ & =(1V^{\prime \prime }(Q)|_{z^2})\frac{S}{g_1}V^{\prime \prime }(Q)|_{z^1}\frac{R}{g_1}\hfill \end{array}$$

Let us now use the second equation of (A.1) to express the $`1`$ in the first factor as $`1=R/tV^{}(Q)|_{z^1}`$, and upon differentiating the first line of (A.1) with respect to $`g_1`$ and multiplying by $`R`$, we may eliminate $`(R/t)S/g_1V^{\prime \prime }(Q)|_{z^1}R/g_1=R+RV^{\prime \prime }(Q)|_{z^0}S/g_1`$, to get

$$\mathrm{\Gamma }_{1,1}=R+(RV^{\prime \prime }(Q)|_{z^0}V^{}(Q)|_{z^1}V^{\prime \prime }(Q)|_{z^2})\frac{S}{g_1}$$

Finally, it is easy to see that the prefactor of $`S/g_1`$ vanishes, as it is nothing but the residue of a total derivative: $`\frac{d}{dz}\left(zV^{}(Q)\right)|_{z^1}=0`$. We finally get the desired result $`\mathrm{\Gamma }_{1,1}=R`$.

A last remark is in order. When all $`g_{2k1}=0`$, $`k=1,2,3,\mathrm{}`$, the partition function (2.1) counts possibly disconnected Eulerian maps, namely with only regular vertices of even valence. In our present formulation, we see that the complete potential $`\frac{x^2}{2t}V(x)`$ is an even function of $`x`$, which simplifies the results considerably. Indeed, we may immediately infer that the planar density of eigenvalues is also even, so that $`a+b=0`$, henceforth $`S=0`$.

Appendix B. Proof of the relation $`R_1+2R_0=(V^{}(Q)|_{z^2})^2`$

Fig. 13: A schematic representation of an Sk-tree (on the left) and of an $`\overline{\mathrm{S}}_l`$-tree (on the right). We represent only the buds and leaves which play a role in the evaluation of the depth of the marked (by a cross) leaf or bud. We also indicate the closing arches in dashed lines.

Let us introduce refined definitions of S-trees as follows. We define an Sk\- (resp. an $`\overline{\mathrm{S}}_l`$\- ) tree as an unrooted S-tree with a marked leaf (resp. a marked bud) of depth $`k`$ (resp. $`l`$). An illustration of these definitions is given in Fig. 13. We denote by $`S_k`$ (resp. $`\overline{S}_l`$) the corresponding generating functions. As the number of leaves and buds of given depth $`k1`$ are the same in any unrooted S-tree, we have $`S_k=\overline{S}_k`$ for all $`k1`$. For $`k=0`$ we have $`\overline{S}_0=0`$ while $`S_0=\mathrm{\Gamma }_1`$ as there is a unique depth $`0`$ leaf in any unrooted S-tree. Moreover we have $`S=_{k0}S_k`$ as the choice of root for an S-tree is that of a leaf of arbitrary depth. Using the first line of (2.1) (for which we have already given a combinatorial proof in Section 3.3), we may express

$$V^{}(Q)|_{z^2}=S\mathrm{\Gamma }_1=\underset{k1}{}S_k=\underset{k1}{}\overline{S}_k$$

and therefore

$$(V^{}(Q)|_{z^2})^2=\underset{k1}{}S_k\overline{S}_k+2\underset{k>l1}{}S_k\overline{S}_l$$

Let us now identify the right hand side of this equation with $`R_1+2R_0`$.

The quantity $`S_k\overline{S}_l`$ counts all trees obtained by gluing into an edge the marked leaf of an Sk-tree with the marked bud of an $`\overline{\mathrm{S}}_l`$-tree. The resulting tree is easily seen to be an unrooted R-tree (of total charge $`2`$) with a particular marked edge. This edge separates the unrooted R-tree into a rooted S-subtree (the former Sk-tree with its marked leaf replaced by a root) and a part of charge $`2`$ (obtained by removing the marked bud of the $`\overline{\mathrm{S}}_l`$-tree). Let us show that the S-subtree is maximal hence the marked edge connects it to the core of the unrooted R-tree. Indeed, the maximal S-tree containing this S-subtree is connected to the core through a $`0`$-$`2`$ edge. This edge separates the $`\overline{\mathrm{S}}_l`$-tree into a piece of charge $`2`$ and one of charge $`1`$, which in an unrooted S-tree is possible only if the part of charge $`1`$ is reduced to a bud. This identifies the above $`0`$-$`2`$ edge with the marked one, and proves that the S-subtree is maximal. Note moreover that the core of the R-tree is actually obtained by gluing the cores of all the R-subtrees of the $`\overline{\mathrm{S}}_l`$-tree originating from the vertex to which the marked bud is attached. The quantity $`S_k\overline{S}_l`$ counts unrooted R-trees with a singularized maximal S-subtree and particular depth restrictions.

Fig. 14: A schematic view of the unrooted R-tree resulting from the connection of an Sk-tree with an $`\overline{\mathrm{S}}_l`$-tree via their marked leaf and bud. The positions of the two unmatched leaves in the closing procedure depend on whether (a) $`k=l`$ (one on each side) (b) $`k>l`$ (both on the left) or (c) $`k<l`$ (both on the right).

Now it is easy to check that if $`k=l`$, the bud-leaf matching procedure leaves us with one unmatched leaf on each side (see Fig. 14). On the other hand, for $`kl`$, the two unmatched leaves are on the same side, which is that of the singularized maximal S-subtree iff $`k>l`$. The right hand side of (B.1) therefore counts the number of unrooted R-trees with a singularized maximal S-subtree and with a marked depth $`0`$ leaf on it. This also counts the unrooted R-trees twice if their two depth $`0`$ leaves lie outside the core (i.e. in either the same or two distinct maximal S-subtrees), and once if exactly one depth $`0`$ leaf lies outside the core (i.e. is in a maximal S-subtree). This is nothing but $`2R_0+R_1`$. QED.

References

relax W. Tutte, A Census of planar triangulations Canad. Jour. of Math. 14 (1962) 21-38. relax W. Tutte, A Census of Hamiltonian polygons Canad. Jour. of Math. 14 (1962) 402-417. relax W. Tutte, A Census of slicings Canad. Jour. of Math. 14 (1962) 708-722. relax W. Tutte, A Census of Planar Maps, Canad. Jour. of Math. 15 (1963) 249-271. relax G. ’t Hooft, A planar diagram theory for strong interactions, Nucl. Phys. B72 (1974) 461-473. relax E. Brézin, C. Itzykson, G. Parisi and J.-B. Zuber, Planar Diagrams, Comm. Math. Phys. 59 (1978) 35-51. relax P. Di Francesco, P. Ginsparg and J. Zinn–Justin, 2D Gravity and Random Matrices, Physics Reports 254 (1995) 1-131. relax B. Eynard, Random Matrices, Saclay Lecture Notes (2000), available at http://www-spht.cea.fr/lectures\_notes.shtml relax G. Schaeffer, Conjugaison d’arbres et cartes combinatoires aléatoires PhD Thesis, Université Bordeaux I (1998). relax G. Schaeffer, Bijective census and random generation of Eulerian planar maps, Electronic Journal of Combinatorics, vol. 4 (1997) R20. relax D. Poulhalon and G. Schaeffer, A bijection for loopless triangulations of a polygon with interior points, proceedings of the conference FPSAC’02, Melbourne (2002), available at http://www.loria.fr/$``$schaeffe/ relax M. Bousquet-Mélou and G. Schaeffer, Enumeration of planar constellations, Adv. in Applied Math., 24 (2000) 337-368. relax D. Poulhalon and G. Schaeffer, A note on bipartite Eulerian planar maps, preprint (2002), available at http://www.loria.fr/$``$schaeffe/ relax P. Di Francesco, B. Eynard and E. Guitter, Coloring Random Triangulations, Nucl. Phys. B516 \[FS\] (1998) 543-587. relax J. Bouttier, P. Di Francesco and E. Guitter, Counting colored Random Triangulations, Saclay preprint SPhT/02-075, cond-mat/0206452, to appear in Nucl. Phys. B. relax see for instance D. Boulatov and V. Kazakov, The Ising model on a random planar lattice: the structure of the phase transition and the exact critical exponents, Phys. Lett. B186 (1987) 379-384. relax J. Bouttier, P. Di Francesco and E. Guitter, Critical and tricritical hard objects on bicolourable random lattices: exact solutions J. Phys. A: Math. Gen. 35 (2002) 3821-3854.


format: text

text:
# Isotopic and spin selectivity of H2 adsorbed in bundles of carbon nanotubes

## 1 Introduction

The adsorption of gases within single-walled, carbon nanotubes (SWNT) has recently attracted broad attention among physicists, chemists, materials scientists and engineers . Experiments have shown that one can create ordered arrays of nanotubes , which form a close-packed bundle (or rope) of SWNT’s. Ideal bundles of nanotubes consist of very long strands of nearly parallel tubes held together at equilibrium separation by intermolecular forces. The small diameter and large aspect ratio of the nanotubes make them interesting systems for gas adsorption. In such a bundle, physisorption may occur at four distinct sites: i) inside the tubes (endohedral adsorption); ii) in the interstitial channels (IC) between three contiguous nanotubes (exohedral adsorption); iii) in the grooves between adjacent nanotubes on the external surface of the bundle, and iv) elsewhere on the outside surface of the bundle. In this paper, we assess the ability of nanotube bundles to preferentially adsorb specific isotope and spin species.

Small atoms and molecules, in contrast to larger molecules, have been predicted to fit well in the ICs of a nanotube bundle . For such species, the IC’s environment provides a large number of neighboring C atoms at nearly optimal distance from the adatom, so the interstitial binding energy is larger than is found in other known environments . The grooves also present strongly attractive environments for both small and large molecules. The binding energies in these sites have been determined experimentally to be between 50 % and 100 % greater than the values on graphite .

Molecular hydrogen inside SWNTs is a particularly appealing system to study, since its adsorption in ICs and grooves is relevant to both gas storage and isotope separation . “Selectivity” is the term used to describe the separation, or selective adsorption, of one species relative to the other species of a mixture. In a mixture of two components, the selectivity of component 1 relative to that of 2 is:

$$S=(x_1/x_2)/(y_1/y_2)$$
(1)
with $`x_i(y_i)`$ the pore (bulk) molar fractions. Since two isotopes have similar sizes, shapes and interaction potentials, the separation of isotopic mixtures is a difficult and energy-expensive process, requiring special experimental techniques, such as cryogenic distillation, diffusion separation, laser isotope separation or microwave molecular separation . However, most of these processes have low selectivity for separating $`H_2`$ isotopes. Recently, a novel separation technique called quantum sieving was predicted to be particularly efficient for nanotubes . Quantum sieving separates lighter molecules from heavier ones by selectively adsorbing heavier molecules. Such selective adsorption can be explained by the higher zero-point energy of the light species, which makes their adsorption relatively unfavorable. Quantum sieving can be implemented when the adsorbate is effectively confined to a one-dimensional (1D) channel or a 2D surface or a small cavity (0D). The selectivities observed for adsorption of $`H_2`$ isotopes on common substrates such as graphite, zeolites, and alumina are low - typically in the range of 1.1 to 3 . However, differences in zero-point energies of the adsorbed species are expected to be particularly large when molecules are confined in very narrow pores. Wang et al. have shown that nanopores with diameters wider than 7 $`\stackrel{̊}{A}`$ exhibit weak selectivity while smaller nanopores were predicted to exhibit large selectivities . Carbon nanotubes typically have diameters larger than 7 $`\stackrel{̊}{A}`$, but the interstitial channels are smaller, so that they have a pore size and solid-fluid potential which can effectively sieve mixtures of $`H_2`$ and $`D_2`$. Indeed, Path Integral Monte Carlo calculations at variable pressures  yielded large isotope selectivity in the IC. Our calculations are performed in the limit of low coverage (virtually zero pressure), using what we believe to be an improved gas-solid interaction potential.

The solid-fluid potential is usually modelled as a pair-wise sum of Lennard-Jones interactions between adatoms and carbon atoms . The potential used in our study is the one recently developed by Kostov et al. . In addition to terms present in previous calculations, Kostov’s potential includes also the interactions of the $`H_2`$ quadrupole moment with the local electrostatic field of C nanotubes (found from ab initio calculations ) and interactions of the $`H_2`$ static multipole moments with the image charges induced on the surrounding nanotubes. The model yields a strongly confining potential for $`H_2`$ and $`D_2`$, of depth about 1800 K ($``$ 160 meV), which alters the rotational spectrum and induces a large difference between isotopic zero point energies. As we will show, the calculated difference between $`H_2`$ and $`D_2`$ binding energies is consistent with recent isosteric heat data of Wilson et al. . This large difference, when exponentiated in an Arrhenius expression at low temperature, yields huge isotope selectivities in the IC and groove, far exceeding that on common sorbents (e.g. graphite and zeolites).

The rotational hindrance of $`H_2`$ molecules adsorbed in carbon nanotube bundles is a controversial issue. There exist experiments showing no alteration of the rotational spectrum  and other experiments (with different samples) which manifest large changes compared to the free rotation energy spectrum . The calculations of Kostov et al. indicate that interstitially adsorbed $`H_2`$ has a significantly hindered rotational motion . The existence of a large rotational barrier leads to large splittings and shifts (relative to an orientationally invariant potential) of the $`H_2`$ and $`D_2`$ rotational energy levels (see Fig.1). In this work we show that a consequence of this hindered rotational motion is a significant spin selectivity in the ICs and grooves.

Recently, Hathorn et al. showed that quantization of the restricted rotational motion of $`H_2`$ and $`D_2`$ confined in SWNTs contributes significantly to quantum sieving . Their study, however, is substantially different from the present work, since they consider $`H_2`$ adsorbed within a single nanotube and approximate the rotational potential as that of the molecule on the axis of the nanotube. In addition, they employed a modified Lennard-Jones potential with $`ϵ`$ and $`\sigma `$ parameters far different from the potential parameters assumed in this work. They find significantly hindered rotational spectrum for the adsorbed molecule for small nanotubes. Their results are qualitatively consistent with our results although we are considering quite different sites.

The outline of this paper is as follows. In the next section, we calculate the ortho-para selectivity as a function of temperature (T), taking into account the hindered rotational motion. In Section III, we derive the specific heat due to the rotational motion of the adsorbed molecules. Section IV presents results for the isotope selectivity of mixtures of $`H_2`$ and $`D_2`$. Section V presents calculations of the isosteric heats of $`H_2`$ and $`D_2`$ and compares these with the experimental results of Wilson et al. . Section VI summarizes and discusses this work.

## 2 Ortho-para selectivity

In order to compute the selectivity of component 1 (e.g. ortho) relative to component 2 (e.g. para) at low coverage, we use the equilibrium condition relating the chemical potential of each component of the 3D vapor(outside the nanotube bundles) to that of the 1D gas adsorbed in the IC:

$$\mu _i^{vapor}=\mu _i^{IC}$$
(2)
where i stands for each component, ortho and para. We employ statistical mechanics to find the chemical potential $`\mu =(F/N)_{T,V}`$, where $`F=k_BT\mathrm{ln}Q`$ is the Helmholtz free energy, $`Q`$ is the canonical partition function and $`k_B`$ is Boltzmann’s constant. In the low density limit, the partition function involves energy contributions from (hindered) rotation, vibration and from translation of the center of mass (CM) in the presence of a transverse confining potential.

In the IC or groove, the molecules are assumed to move freely along the axial direction, so the fluid is described by a constant(mean) 1D density, $`\rho `$. The effect of corrugation on the motion is sensitive to assumptions about registry of adjacent tubes  and is neglected here. The CM transverse oscillation is described by the one-particle partition function, $`q^{cm}`$. The internal degrees of freedom consist of the intra-molecular vibrations and hindered rotations, described by $`q^{vib}`$ and $`q^{rot}`$, respectively. Therefore:

$$F_i^{IC}=Nk_BT\mathrm{ln}[\frac{e}{\rho _i\lambda }q_i^{cm}q_i^{vib}(g_iq_i^{rot})\mathrm{exp}(C_i/(k_BT))]$$
(3)
$$\mu _i^{IC}=C_ik_BT\mathrm{ln}(\frac{e}{\rho _i\lambda })k_BT\mathrm{ln}q_i^{cm}k_BT\mathrm{ln}q_i^{vib}k_BT\mathrm{ln}(g_iq_i^{rot})$$
(4)
where i stands for para and ortho, $`\lambda =\sqrt{2\pi \mathrm{}^2/(mk_BT)}`$, $`g_i`$ is the nuclear spin degeneracy and $`C_i`$ is the interaction energy experienced by molecules due to the confining environment.

The partition function of the coexisting gas outside the nanotube bundle has contributions from the 3D translation motion and from internal degrees of freedom (vibration and free rotation of the molecules). The free energy and the chemical potential are then:

$$F_i^{vapor}=Nk_BT\mathrm{ln}[\frac{e}{n_i\lambda ^3}q_i^{vib}(g_iq_i^{free})]$$
(5)
$$\mu _i^{vapor}=k_BT\mathrm{ln}\frac{e}{n_i\lambda ^3}k_BT\mathrm{ln}q_i^{vib}k_BT\mathrm{ln}(g_iq_i^{free})$$
(6)
where $`n_i`$ is the 3D density of ortho or para vapor outside the nanotube bundle.

Experimental data of Williams et al and calculations of Kostov et al indicate that the intra-molecular vibrations are essentially the same in the IC as in free space . The CM oscillations depend only on the molecular mass, so they give the same contribution to the ortho and para partition functions. Finally, the only relevant contributions to the spin selectivity come from the even and odd rotational partition functions, denoted $`q_+`$ and $`q_{}`$, respectively:

$$q_+=\underset{j=even}{}\underset{m=j}{\overset{j}{}}g_{jm}\mathrm{exp}(\beta ϵ_{jm})$$
(7)
where $`\beta =1/(k_BT)`$, $`g_{jm}`$ is the rotational degeneracy (not including the nuclear degeneracy), and $`ϵ_{jm}`$ are the rotational levels. For the odd partition function ($`q_{}`$), the summation in $`Eq.7`$ is done only over the odd rotational levels. In free space, both $`H_2`$ and $`D_2`$ can be modelled as free quantum rotors with $`g_{jm}=j(j+1)`$ and $`ϵ_{jm}=j(j+1)B`$. Here j is the rotational quantum number and $`B=\mathrm{}^2/(2I)`$, the rotational constant of $`H_2`$ or $`D_2`$. However, the different character of the nuclear spin (fermionic for $`H_2`$ and bosonic for $`D_2`$) imposes different conditions on the total wavefunction (antisymmetric for $`H_2`$ and symmetric for $`D_2`$), so that the nuclear degeneracies corresponding to even and odd rotational states are different: $`g_+^{H_2}=1`$, $`g_{}^{H_2}=3`$ and $`g_+^{D_2}=6`$, $`g_{}^{D_2}=3`$. By convention, the name ortho is given to the species with larger statistical weight, so the even states are called para for $`H_2`$ and ortho for $`D_2`$. Thus, at high T, the ortho and para molecules exist in the ratio of 3:1 in the case of $`H_2`$ and 6:3 in the case of $`D_2`$.

The spin selectivity in the low coverage (pressure) limit involves a ratio of densities, which can be found by imposing the equilibrium condition ($`Eq.2`$):

$$S=\frac{(\rho _{}/\rho _+)_{IC}}{(n_{}/n_+)_{vapor}}$$
(8)
where $``$ and $`+`$ subscripts correspond to the odd and even rotational states, respectively (ortho and para for $`H_2`$, para and ortho for $`D_2`$). All similar factors mentioned previously cancel out in $`S`$ and, in the end, the low coverage selectivity depends on just the even and odd partition functions, $`q_+`$ and $`q_{}`$:

$$S=\frac{(q_{}^{IC}/q_+^{IC})}{(q_{}^{vapor}/q_+^{vapor})}$$
(9)
A qualitative understanding of the spin selectivity can be achieved by examining the rotational spectrum. The ortho-para ratio in the IC, in the limit of high T, is the same as in the gas phase. But since the energy spectrum changes in the IC, the overall partition functions differ, with particularly large effects at low T. Fig.1 contrasts the rotational spectra of $`H_2`$ in the IC (as computed by Kostov) and in free space. In the IC, the reference energy is the energy the molecule would have if it does not rotate ($`C_{H_2}`$=-1275 K). Notice that the ground (para) state is shifted down from this value by $``$ 200 K and the first excited state (ortho $`j=1,m_j=0`$) is also shifted down ($``$ 60 $`K`$), close to the ground state. The $`D_2`$ rotational spectrum undergoes a similar alteration in the IC, but with a smaller shift of the ground (ortho) state (about 130 K) with respect to a different reference energy ($`C_{D_2}`$=-1445 K), and a slightly larger excitation energy to the ($`j=1,m_j=0`$) para state (about 85 K). The shift in the first excited state (ortho for $`H_2`$, para for $`D_2`$) towards the ground state (para for $`H_2`$, ortho for $`D_2`$) causes the ortho to para ratio of $`H_2`$ molecules, or para to ortho ratio of $`D_2`$ molecules, to increase substantially in the IC with respect to that in free space. In the groove site, these ratios are still larger than in free space, but smaller, however, than in the IC, due to the lesser confinement between the two nanotubes in the groove site. At low T, all but the lowest states yield negligible contributions to $`q^{rot}`$, so that the selectivity depends only on the values of the lowest excitation energies in the IC and free space.

$$S_o=\frac{\mathrm{exp}[\beta (ϵ_{10}ϵ_{00})]+2\mathrm{exp}[\beta (ϵ_{11}ϵ_{00})]}{3\mathrm{exp}[2\beta B]}$$
(10)
where the free space separation is twice the rotational constant B=85 K for $`H_2`$.

The results for the ortho-para selectivity are shown in Fig.2. We find that the spin selectivity is larger than 1 for an extended range of temperatures, up to 100 K in the IC and $``$ 75 K in grooves. This is consistent with the difference between the excitation energies in the IC (groove) and in free space which enters $`Eq.10`$. As discussed above, the alterations of the rotational spectrum favor the $`H_2`$ ortho species, whereas for $`D_2`$, the para molecules are preferentially adsorbed. Since the $`H_2`$ molecules have a larger zero point motion than $`D_2`$ molecules, they experience more rotational hindrance and their spectrum is more altered in the corresponding energy scale. Since the reference energy does not count for the spin selectivity, this exhibits larger values for $`H_2`$ than for $`D_2`$.

## 3 Rotational specific heat

In order to compute the thermodynamic properties of the adsorbate, one works with the composite equilibrium partition function:

$$q_{rot}^{H_2}=S(2S+1)q_++(S+1)(2S+1)q_{}$$
(11)
$$q_{rot}^{D_2}=(S+1)(2S+1)q_++S(2S+1)q_{}$$
(12)
where $`S=1/2`$ for $`H_2`$ and $`S=1`$ for $`D_2`$. We will discuss the $`H_2`$ case and point out the differences in the case of $`D_2`$. In many circumstances, $`H_2`$ is not in thermal equilibrium as regards the relative magnitudes of para and ortho components, because the probability of flipping the nuclear spin is very small (the lifetime in free space is $``$ one year). Therefore, the transition probability of a molecule from one nuclear spin state to another is negligible during a specific heat experiment. Consequently, the sample may be viewed as a non-equilibrium mixture of two independent species, which give additive contributions to the heat capacity.

$$C_{noneq}=f_pC_p+f_oC_o$$
(13)
Here $`f_p`$ and $`f_o`$ are the para and ortho molar fractions. In an equilibrium mixture these depend on temperature, whereas in a non-equilibrium mixture their relative concentrations are determined by the initial conditions. In our calculations, we let the para and ortho molar fractions have their high T values: $`f_p=1/4`$ and $`f_o=3/4`$. The para and ortho specific heats (in the low density limit) are:

$$C_{p,o}=N_{p,o}k_B\frac{d}{dT}(T^2\frac{d}{dT}\mathrm{ln}q_{+,})$$
(14)
At low T, only the lowest levels contribute to the specific heat. Let us consider only the first two levels: $`ϵ_{00}`$ and $`ϵ_{20}`$ for the para species, and $`ϵ_{10}`$ and $`ϵ_{11}`$ for the ortho species. Applying the formulae above, the para and ortho specific heats are:

$$C_p=Nk_B\frac{(ϵ_{20}ϵ_{00})e^{\beta (ϵ_{20}ϵ_{00})}}{(1+e^{\beta (ϵ_{20}ϵ_{00})})^2}$$
(15)
$$C_o=Nk_B\frac{2(ϵ_{11}ϵ_{10})e^{\beta (ϵ_{11}ϵ_{10})}}{(1+2e^{\beta (ϵ_{11}ϵ_{10})})^2}$$
(16)
Thus, at low T, the para and ortho heats in the IC depend only on the excitation energies $`(ϵ_{20}ϵ_{00})`$ and $`(ϵ_{11}ϵ_{10})`$, respectively. Fig. 3 (a,b) shows a comparison between the para, ortho and non-equilibrium specific heats in free space and in the IC. Notice that, at low T, in the IC the main contribution to the net specific heat comes from $`C_o`$, whereas in free space it comes from $`C_p`$. This is related to the finite excitation energy ($`ϵ_{11}ϵ_{10}`$) between (j=1) ortho states in IC, which are degenerate in free space (see fig.1). It can be shown that the peak in the heat capacity occurs at $`T\mathrm{\Delta }/(3k_B)`$, where $`\mathrm{\Delta }`$ is the first excitation energy present in a specific problem. Therefore, in free space the para peak is at about 80 K, while in the IC the ortho peak is at about 75 K.

In some environments, however, the ortho-para conversion may occur quickly (e.g. due to magnetic impurities). In this case, the net rotational specific heat is obtained from the composite partition function ($`Eqs.11,12`$):

$$C_{eq}=Nk_B\frac{d}{dT}(T^2\frac{d}{dT}\mathrm{ln}q_{rot})$$
(17)
At low T, the equilibrium specific heat depends only on the $`(ϵ_{10}ϵ_{00})`$ excitation energy:

$$C_{eq}=Nk_B\frac{(ϵ_{10}ϵ_{00})e^{\beta (ϵ_{10}ϵ_{00})}}{(1+e^{\beta (ϵ_{10}ϵ_{00})})^2}$$
(18)
The equilibrium specific heats of $`H_2`$ and $`D_2`$ in the IC are shown in Fig.4. In the case of $`H_2`$, there is a substantial ortho-para conversion peak at about 20 K, followed by a gentle rise above 75 K, corresponding to ortho-ortho excitation. However, in the case of $`D_2`$, the ortho-para conversion appears only as a small bump at about 25 K, since the molar weights of para and ortho species are now different. At high T, the rotational specific heat per molecule goes to $`k_B`$, the classical result for a free rotor.

## 4 Isotope slectivity

Let us now consider a mixture of $`H_2`$ and $`D_2`$. The equilibrium condition between the 1D gas adsorbed in the ICs or grooves and the coexisting vapor outside the nanotube bundles is given again by $`Eq.2`$: $`\mu _i^{IC}=\mu _i^{vapor}`$, where now $`i`$ stands for $`H_2`$ or $`D_2`$. The partition functions and the chemical potentials of each species consist of the same factors, but their meanings are different. For example, the rotational partition function is given now by the composite partition function ($`Eqs.11,12`$) of each isotope. Since the $`H_2`$ and $`D_2`$ species have different masses, their zero point energies will be different, and these will give a nontrivial contribution to the selectivity (as will differences in the rotational energies). The transverse CM oscillations of the molecule in the IC or grooveare treated as excitations of a 2D harmonic oscillator. The corresponding partition function is:

$$q_{vib}=[\underset{n=0}{\overset{\mathrm{}}{}}e^{\beta \mathrm{}\omega (n+1/2)}]^2=\frac{e^{\beta \mathrm{}\omega }}{(1e^{\beta \mathrm{}\omega })^2}$$
(19)
where $`\omega `$ is the frequency of CM oscillations, found from the force constant. The rotational hindrance is expressed in terms of the rotational energy shifts, $`ϵ_{mj}`$ and the additive constants (C) in the expansion of the rotational part of the interaction potential . This additive (reference) energy is also diferent for $`D_2`$ and $`H_2`$. The rotational hindrance leads to a shift of the binding energy of each species: $`E_b=(C+\mathrm{}\omega +ϵ_{00})`$.

The equilibrium condition for each species ($`i=D_2,H_2`$) yields the ratio of IC $`(\rho _i)`$ and vapor ($`n_i`$) densities for each species, which is then used to find the isotope selectivity.

$$\frac{\rho _i}{n_i}=\lambda _i^2\frac{\mathrm{exp}(\beta \mathrm{}\omega _i)}{(1\mathrm{exp}(\beta \mathrm{}\omega _i))^2}\frac{\mathrm{exp}(\beta C)q_{rot}^{IC}}{q_{rot}^{free}}$$
(20)
The isotope selectivity is:

$$S_0=\frac{\rho _{D_2}/\rho _{H_2}}{n_{D_2}/n_{H_2}}$$
(21)
$`Eq.20`$ implies that at low T the main contribution to selectivity comes from the difference between the binding energies

$$S_0\frac{1}{2}\mathrm{exp}[\beta (C_{D_2}+\mathrm{}\omega _{D_2}+ϵ_{00}^{D_2}C_{H_2}\mathrm{}\omega _{H_2}ϵ_{00}^{H_2})]=\frac{1}{2}exp[\beta (E_b^{D_2}E_b^{H_2})]$$
(22)
Fig.5 displays the results of our calculations for isotopic selectivity in the ICs and grooves of an (18,0) nanotube bundle, and on graphite (for comparison). All three surfaces favor $`D_2`$. However, the isotope selectivity in the IC and groove is much larger than that on graphite and within the tubes . Even at 300 K the selectivity in the IC and groove is substantial, due to the large difference between the $`D_2`$ and $`H_2`$ binding energies. This large difference comes primarily from the zero point energy (but the difference between the rotational energies also contributes). Table 1 presents our results for the different contributions to the the binding energies. The difference between $`D_2`$ and $`H_2`$ binding energies is $``$ 400 K in the IC and $``$ 170 K in the groove, much larger than the difference of $``$ 35 K on graphite . When exponentiated, these differences lead to the very large selectivities found in the IC and groove. We may compare our results with previous calculations. Challa et al. have studied isotope selectivity at variable pressure in the interstices of (10,10) nanotube bundle, using Ideal Adsorbed Solution Theory(IAST) and Path Integral Grand Canonical Monte Carlo simulations (PI-GCMC) . At low pressures, the IAST and PI-GCMC results converge to the analytic zero pressure results. In their calculation, only the transverse zero-point energy is considered, i.e. rotational motion is ignored. The selectivity calculated by them at 20 K and very low pressure ($`10^{13}`$ atm) is about 700. (10,10) and (18,0) nanotubes have similar diameters (13.6 and 13.8, respectively), so the interstices of a closed packed nanotube bundle are similar, too ($``$ 6 $`\stackrel{̊}{A}`$ diameter). Recently, Hathorn et al.  have estimated the quantum effects due to the restriction of the rotational motion of $`H_2`$ and $`D_2`$ adsorbed within a single nanotube. Their results reveal significantly enhanced rotational separation factor(selectivity), which at $`20`$ K is of order $`10^3`$. The zero pressure isotope selectivity calculated by us at 20 K, in the interstices of bundles of (18,0) carbon nanotubes is of order $`10^9`$. Our results differ to such a large degree since we considered both zero point and rotational contributions to selectivity and we used a slightly deeper, solid-fluid potential. The stronger the confinement, the larger the difference between the isotope zero point energies and between the rotational energies. Previous calculations inside (5,5) nanotubes at 20 K yielded much lower isotope selectivities . At low pressure, the computed selectivity of tritium to $`H_2`$ is 23, and that of tritium to $`D_2`$ is 1.7; hence $`S_0(D_2/H_2)=13`$. Thus, even if the inner channels of (5,5) nanotubes have about the same radius (3.3 $`\stackrel{̊}{A}`$) as the ICs of (18,0) nanotube bundles (3 $`\stackrel{̊}{A}`$), the computed selectivity in the IC is much larger due to the different arrangement of the carbon atoms in the IC, which yields greater confinement.

## 5 Isosteric heat

The isosteric heat is defined as:

$$Q_{st}=[\frac{d(lnP)}{d\beta }]_N$$
(23)
where P is the pressure of the vapor outside the nanotube bundles. In our model, at low density, the dependence of P on T can be found from the equilibrium condition $`(Eq.2)`$ using the ideal gas formula: $`n=\beta P`$. Finally, the isosteric heat is:

$$Q_{st}=2k_BT\mathrm{}\omega C+(\frac{d}{d\beta })[ln\frac{(q_{rot}^{IC}/q_{rot}^{free})}{(1e^{\beta \mathrm{}\omega })^2}]$$
(24)
At low T, the excitations can be neglected and the term involving the derivative with respect to $`\beta `$ goes to $`ϵ_{00}(ϵ_{10}ϵ_{00})e^{\beta (ϵ_{10}ϵ_{00})}+2Be^{2\beta B}`$, where $`ϵ_{00}`$ and $`ϵ_{10}`$ are the first two rotational energy levels. Therefore, at low T, the isosteric heat is a measure of the binding energy:

$$Q_{st}=E_b+2k_BT$$
(25)
since $`E_b=(C+\mathrm{}\omega +ϵ_{00})`$.

The isosteric heat can be found experimentally by taking the difference between two nearby isotherms. Recently, there have been reported experiments of adsorption of $`D_2`$ and $`H_2`$ in carbon nanotube bundles , which show isosteric heats in nanotube bundles to be a factor of 1.5 (for $`H_2`$) to 1.8 (for $`D_2`$) larger than those on graphite. This means that the $`D_2`$ and $`H_2`$ binding energies are almost twice as large as those on graphite. Moreover, the difference between isotope isosteric heats, which again is related to the difference between isotope binding energies, was found to be about 200 K, much greater than the difference of about 35 K on graphite . This is a mark of the greater confinement of molecules adsorbed in nanotube bundles than on the surface of graphite; the confinement leads to the enhanced separability of isotopes found in the previous section.

We performed calculations of the isotope isosteric heats in the low density limit ($`Eq.24`$). Table 2 and Figs. 6(a),(b) compare the results of the calculated isosteric heats and the low coverage experimental data at T=85 K. The experimental procedure of Wilson et al does not reveal where in the nanotube bundle the adsorption occurs. Since the experimental isosteric heats are much larger than on graphite, one assumes that the adsorption environment should be either the IC or the groove site. Our calculations of the isosteric heats in both IC and groove are lower than the experimental results. The calculated difference between $`D_2`$ and $`H_2`$ isosteric heats (374 K) in the IC is much larger than the one found experimentally, but the difference (137 K) in the groove is closer to the experimental difference (200 K). In addition, the calculated binding energies(see Table I) for $`D_2`$ and $`H_2`$ in the grooves are consistent with the experimental results and show that the adsorption on the grooves is much more likely than on the IC at low T. Assuming the accuracy of our approximations, we believe that our results are compatible with an adsorption in the groove, rather than the IC. This is also experimentally plausible, since the external grooves are more accessible for adsorption than the ICs. Indeed, the ICs may be completely blocked, preventing any adsorption there.

## 6 Summary and conclusions

We have investigated the spin (ortho versus para) and the isotope ($`D_2`$ versus $`H_2`$) selectivity at low coverage in the IC and groove channel of (18,0) nanotube bundles. The rotational hindrance of molecules in the IC and groove site induces shifts and splittings in the rotational spectrum, enabling the nanotube bundles to preferentially adsorb $`H_2`$ ortho molecules and $`D_2`$ para molecules. Our calculations show substantial spin selectivity for temperatures up to $``$ 100 K in the IC and up to $``$ 75 K in the groove. At low temperatures, the alterations of the rotational spectrum in the IC induces new features in the ortho heat capacity, as the ortho-ortho peak. The $`H_2`$ equilibrium heat capacity exhibits a distinctive ortho-para conversion peak at $``$ 20 K. The different features of non-equilibrium and equilibrium heat capacities may be a way to check the existence of the ortho-para conversion of $`H_2`$ and $`D_2`$ adsorbed in nanotube bundles. In the case of the isotope selectivity in nanopores, zero-point motion favors the heavier isotope , and the rotational motion enhances the preferential adsorption of the heavier isotope . We find substantial isotope selectivity even at temperatures as high as 300 K in the IC and groove. At 20 K, our calculation at low coverage (consequently zero pressure) yield isotope selectivities of order $`10^9`$ in the IC, orders of magnitude larger than the low pressure limit found in the calculations of Wang et al. and Challa et al.. Our much larger isotope selectivity is a consequence of the stronger confinement in the IC, which yields larger zero point and rotational energies. The difference between the $`D_2`$ and $`H_2`$ isosteric heats in the groove ($``$ 135 K) is close to the value found experimentally ($``$ 200 K). We conclude that carbon nanotube bundles are ideal surfaces for spin and isotope separation, providing a potential technology of quantum sieving.

We acknowledge stimulating discussions with J. K. Johnson, M. M. Calbi, T. Wilson, M. J. Bojan, O. Vilches, P. Sokol, D. Narehood and D. Stojkovic. This research has been supported by Petroleum Research Foundation of the American Chemical Society and Air Products and Chemicals, Inc.

Figure captions

Fig. $`1`$. $`H_2`$ rotational spectrum in free space and in IC. The number of lines in a level represents the degeneracy of that energy level. o and p mean ortho and para, respectively.

Fig. $`2`$. $`D_2`$ and $`H_2`$ spin selectivity in IC and groove channel as a function of T

Fig. $`3`$. Ortho (dashed line), para (dotted line) and non-equilibrium (full line) specific heat (without spin equilibration) of $`H_2`$ molecules (a) in IC, (b) in free space

Fig. $`4`$. The equilibrium rotational specific heat of $`H_2`$ (full line) and of $`D_2`$ (dashed line)

Fig. $`5`$. Isotope selectivity in the IC (full line), groove (dashed line) and on graphite (dotted line).

Fig. $`6`$. The calculated $`H_2`$ (dotted line), $`D_2`$ (dashed line) and the difference $`D_2H_2`$ (full line) isosteric heats (a) in the IC and (b) in the groove. The experimental values at 85 K are shown in symbols: circle for $`D_2`$, square for $`H_2`$ and triangle for the difference $`D_2H_2`$.