\documentclass[submission]{FPSAC2024}
\usepackage[utf8]{inputenc}


%
\usepackage{graphicx, enumitem}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{amsthm, xcolor}
\usepackage{comment, placeins,hyperref}

%\usepackage[margin=7mm]{caption}

\theoremstyle{theorem}
\newtheorem{theorem}{Theorem}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{conjecture}[theorem]{Conjecture}
\numberwithin{theorem}{section}

\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{notation}[theorem]{Notation}
\newtheorem{method}[theorem]{Method}
\newtheorem{example}[theorem]{Example}


\newcommand{\singeq}{\hat{\sim}}
%
%
\usepackage[backend=bibtex]{biblatex}
\addbibresource{biblio.bib}






%

\title[On the scaling of random Tamari intervals]{On the scaling of  random Tamari intervals and Schnyder woods of random triangulations (with an asymptotic D-finite trick)}

\author{Guillaume Chapuy%
\thanks{Université Paris Cité, CNRS, IRIF, F-75013, Paris, France
Email:~{\tt guillaume.chapuy@irif.fr}.
}
}

\date{\today}

\abstract{ 
	We consider a Tamari interval of size $n$ (i.e., a pair of Tamari-comparable Dyck paths) chosen uniformly at random. We show that the typical height of points on the paths scales as $n^{3/4}$, with an explicit limit law.
By the Bernardi-Bonichon bijection, 
this also applies to 
Schnyder trees of random plane triangulations.
	The exact solution of the model is based on polynomial equations with one and two catalytic variables. 
	To deduce convergence in law, we use a simple analytic method based on D-finiteness, which is essentially automatic.
%	We are not sure to have seen this simple trick before. 

%	The same method applies (with the same limit law) to the parking trees of Contat-Curien-Hénard. 
	%It would be interesting to study the universality of this convergence  for decomposition trees associated to positive Bousquet-Mélou--Jehanne equations.
}


%
%



\begin{document}


\maketitle





\section{Introduction and main results}


\begin{figure}[h]
	\centering
	\vspace*{2mm}
	\includegraphics[width=0.45\linewidth]{defTamari.pdf}
	\vspace*{-2mm}
	\hspace{1mm}
	\hfill
	\includegraphics[width=0.26\linewidth]{intervalle.png}
	\hfill
	\includegraphics[width=0.26\linewidth]{ratio.png}
	\caption{Left: The covering relation defining the Tamari lattice.
	%\caption{
	Right: a  uniform random Tamari interval $(P_n,Q_n)$ of size $n=65536$  (blue) generated with a Python code generously  provided by Wenjie Fang, and the plot of $Q_n/P_n$ (purple).}
	\label{fig:defTamari}
	%\label{fig:randomPicture}
\end{figure}

For $n\geq 0$, a \emph{Dyck path of size $n$} is a lattice path made of $n$ up-steps and $n$-down steps, starting (and ending) at height $0$, and whose height stays always nonnegative. %See Figure~\ref{fig:defTamari}.
%Dyck paths of size $n$ are in immediate bijection with well-formed parenthesis words of length $2n$, and they are counted by the Catalan number
%$$
%Cat(n):=\frac{1}{n+1} {2n \choose n}.
%$$
The set $\mathcal{D}_n$ of Dyck paths of size $n$ 
is
%can be 
endowed with the 
%several interesting (partial) orders. Maybe the most interesting and natural one is the 
\emph{Tamari partial order}, whose covering relation is
described as follows.
%nothing but the edge-flip once interpreted through the classical bijection between Dyck paths and triangulations of a polygon. The Tamari partial order is a lattice, whose covering relation can be described as follows:
Let $\mathfrak{p}$ be a Dyck path, and let $d$ be a down-step in $\mathfrak{p}$, followed by an up-step. Let $\mathfrak{e}$ be the shortest excursion following $d$ in $\mathfrak{p}$ (an excursion is a path staying higher than its starting point except for its last point). Then the Dyck path $\mathfrak{q}$ obtained from $\mathfrak{p}$ by exchanging $d$ and $\mathfrak{e}$ is declared \emph{larger} than $\mathfrak{p}$ for the Tamari order. The reflexive transitive closure of this relation defines the Tamari lattice. See Figure~\ref{fig:defTamari}.

The Tamari lattice plays an important role in many facets of algebraic combinatorics and discrete mathematics, in relation with polytopes (the associahedron), flip graphs, hyperbolic geometry, or mixing times of Markov chains.
%To name a few: its Hasse-diagram is the 1-skeleton of a famous polytope called the associahedron, see e.g.~\cite{TamariFestschrift,PilaudSantosZiegler}; determining its diameter has been a longstanding problem solved only recently~\cite{diameterPournin}, despite mysterious connections to hyperbolic geometry~\cite{diameterTarjan}; and, quite excitingly, determining  the mixing time on this graph is still an open problem despite recent partial results~\cite{McShineTetali,Eppstein2022}.

This paper is motivated by another famous connection, between Tamari intervals and planar maps.
%yet another 
%the connection, related to \emph{Tamari intervals}. Following the standard combinatorics terminology, we define  
An \emph{interval} of size $n$ in the Tamari lattice is a pair $(\mathfrak{p},\mathfrak{q}) \in (\mathcal{D}_n)^2$ such that $\mathfrak{p}\leq \mathfrak{q}$ (for the Tamari partial order).
%\footnote{Note that the interval is only the \emph{pair} of comparable paths, not all elements in between. We follow the standard terminology of the field.}. 
We let $\mathcal{I}_n$ be the set of Tamari intervals of size $n$.
In a famous paper, Chapoton proved that the number of elements of $\mathcal{I}_n$ was given by 
%\begin{align}\label{eq:Chapoton}
%|\mathcal{I}_n| = 
$\frac{2}{n(n+1)}{4n+1 \choose n-1}$,
%\end{align}
which is also the number of rooted planar triangulations
%\footnote{defined as embedded planar graphs with only faces of degree $3$. These objects are much more complex than the triangulations~\emph{of a single polygon} Dyck paths are in bijection with.}
of size $n$~\cite{Tutte:censusTriangulations}. An elegant, and deep, direct bijective proof has been given by Bernardi and Bonichon in~\cite{BernardiBonichon}.
Since Chapoton's discovery, the analogy between Tamari intervals and planar maps has been much developped, from the existence of refined product formulas~\cite{BMCPR} strongly resembling the ones appearing in enumerative geometry, to numerous coincidences between the enumeration of special families of intervals and planar maps (e.g.,~\cite{Fang}).
These phenomena are still very partially understood and are the subject of active research.

%One the one hand, the enumeration of Tamari intervals can be refined (or weighted), giving rise to product formulas~\cite{BMCPR} surprisingly similar to the ones counting planar constellations~\cite{BMS}, suggesting the existence of an underlying object maybe as rich as the GUE 1-matrix model, with a potentially integrable structure, which is yet to be discovered.
%On the other hand, many different variants of Tamari intervals have been introduced and enumerated. Their enumeration exhibits a universality phenomenon very similar to the one existing in the world of planar maps, and in fact many of these variants are equi-enumerated with some well-known family of planar maps. The understanding of this phenomenon is still very partial, see~\cite{Fang} and references therein.
%More recently, non-lattice variants of the Tamari partial order have been introduced for which the similarity to maps is even more stunning~\cite{BMChapoton}.


%\subsection{Main results}

Although large random planar maps have been intensely studied in the last decades (see e.g.~\cite{ChassaingSchaeffer,
%BDFG:distances, 
LeGallsurvey}), it seems that the behaviour of random Tamari intervals has not been studied, and we are not sure to know why. It is however natural to ask this question:\smallskip\\
%\begin{quote}
	\smallskip{\hspace{3cm} \it What does a large uniformly random Tamari interval look like?} 
%\end{quote}
\\
\noindent  In this paper we give a first answer to this question. 
If $P\in \mathcal{D}_n$ and $i\in [0..2n]$, we write $P(i)$ for the height of the point of $P$ lying at abscissa $i$.
We show:
%are not aware of any results in this direction. We show:
%We will see that quasi-vanilla methods in analytic combinatorics are enough to show that the correct scaling order for the paths is $n^{3/4}$, at least for typical points. More precisely, we show:
\begin{theorem}[Main result]\label{thm:tamari}
	Let $(P,Q)$ be a Tamari interval of size $n$, chosen uniformly at random in $\mathcal{I}_n$.
	Let $I\in[0..2n]$ be an integer chosen uniformly at random.
	%, and let $Q_n(I)$ be the height of the point of the upper path $Q$ lying at abscissa $I$. 
	Then we have the convergence in law
	\begin{align}\label{eq:limitlaw}
	\frac{Q_n(I)}{n^{3/4}} \longrightarrow Z 	=(XY)^{1/4}
	%,
	%$$
	%when $n$ goes to infinity. Here the random variable $Z$ is equal to 
	%$$
	\end{align}
	when $n$ goes to infinity, where $X\sim \beta(\frac{1}{3},\frac{1}{6})$ and $Y\sim \Gamma(\frac{2}{3},\frac{1}{2})$ are independent random variables. In fact, we have the convergence of all moments: for integer $k\geq 0$,
	\begin{align}\label{eq:mainMoment}
		\mathbf{E}\left[ \left(\frac{Q_n(I)}{n^{3/4}}\right)^k \right]\longrightarrow \frac{\sqrt{3} \cdot2^{-\frac{k}{4}-1}}{\sqrt{\pi}} \frac{\Gamma(\frac{1}{4}k+\frac{1}{3})\Gamma(\frac{1}{4}k+\frac{2}{3})}{\Gamma(\frac{1}{4}k+\frac{1}{2})}.% = \frac{2\Gamma(3g)}{3\Gamma(2g)}.
	\end{align}
	For the lower path we have similarly, with again the convergence of all moments,
	\begin{align}\label{eq:limitlawlow}
		\frac{P_n(I)}{n^{3/4}} \longrightarrow \frac{Z}{3}.
	\end{align}
\end{theorem}
We recall that the random variables $\beta(a,b)$ and $\Gamma(\ell, \theta)$ have respective densities
$\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}x^{a-1}(1-x)^{b-1}$ on $(0,1)$
and 
$
\frac{x^{k-1}e^{-x/\theta}}{\Gamma(k)\theta^k}	
$ on $\mathbb{R}_+$. Their respective $k$-th moments are
$\frac{\Gamma(a+b)\Gamma(a+k)}{\Gamma(a)\Gamma(a+b+k)}$, and  $\frac{\theta^k\Gamma(l+k)}{\Gamma(\ell)}$, so it is direct to check that~\eqref{eq:mainMoment} with $k$ substituted by $4k$ is indeed equal to the $k$-th moment of $XY$. 
Interestingly, the random variable $Z^4$ already appears in a (seemingly unrelated) physics context, see~\cite{physics} or OEIS:A064352.
%This reference also gives a closed form for the density (but does not identify it as a product of a Gamma and a Beta variable)\footnote{I thank Philippe Biane for help identifying the variable from its moments, and Thomas Budzinski for suggesting to look for them in the OEIS, which is how I found the reference~\cite{physics}, see the sequence OEIS:A064352.}.

%\medskip
%
%One can also address the case of the lower path, at the price of more technical, and maybe less standard, methods. We obtain the following theorem.
%%, maybe surprising, fact, which roughly speaking says that the lower path is almost everywhere and with high probability half of the upper path.
%\begin{theorem}\label{thm:tamariLow}
%	Let $(P,Q)$ be a Tamari interval of size $n$, chosen uniformly at random in $\mathcal{I}_n$. 	Let $I\in[0..2n]$ be an integer chosen uniformly at random, and let $P_n(I)$ 
%	%and $Q_n(I)$ 
%	be the height of the point of the lower path $P$
%	%and upper path $Q$, respectively, 
%	lying at abscissa $I$. Then we have the convergence in law
%	\begin{align}\label{eq:limitlawLow}
%		\frac{P_n(I)}{n^{3/4}} \longrightarrow \frac{1}{3}Z.
%	%,
%	%$$
%	%when $n$ goes to infinity. Here the random variable $Z$ is equal to 
%	%$$
%	\end{align}
%	when $n$ goes to infinity, where $Z$ is an in theorem~\ref{thm:tamari}. Again, we have the convergence of all moments.
%\end{theorem}
%%\noindent We suspect that experts in the Bernardi-Bonichon bijection might be able to understand, at least heuristically (and maybe even prove) Theorem~\ref{thm:tamariLow} directly from that bijection, by interpreting one of the two heights as the value of a random walk whose length is proportional to the other height (and applying the central limit theorem) -- but this goes far beyond our knowledge.
%%

\par
In view of Theorem~\ref{thm:tamari} and simulations (Figure~\ref{fig:defTamari})
%and~\ref{thm:tamariLow}, 
it is natural to suspect that
%\begin{align}\label{eq:conjThird}
	$P_n(I)$ is close to $\frac{1}{3} Q_n(I)$ with high probability.
%\end{align}
%with a little-o in probability.
%, i.e. that, roughly speaking, the lower path is almost a third of the lower path, asymptotically almost surely and almost everywhere. See Figure~\ref{fig:randomPicture} which shows a simulation for $n=10^5$, which supports this conjecture.
%
%\begin{remark}
Unfortunately our methods based on functional equations make it hard to track the joint law of  $(P_n(I),Q_n(I))$. However, they can handle the joint law of $(\tilde{P}_n(J),\tilde{Q}_n(J))$ where $J$ is uniform on $[1,n]$ and where for a path $P$ we write  $\tilde{P}(j)$ for the height of the $j$-th up-step of $P$.
	%The corresponding equation is close to the one of Theorem~\ref{thm:tamariLow}, and the resolution is close to the one of Section~\ref{sec:solutionLow}, with the difference that the catalytic equation to solve in the analogue of Section~\ref{sec:eliminating2} is quadratic instead of linear. In this way one can prove that 
	In the long version of this paper we will prove with similar methods\footnote{One can write a functional equation for the generating function of intervals $(P,Q)$ of size $n$ with a marked $j\in [n]$ with control on the parameter $\tilde{Q}(j)-3\tilde{P}(j)$, very similarly to Equation~\eqref{eq:catalyticTamariLower}. The resolution of this equation is very similar to what is done in Section~\ref{sec:solutionLow}, with the difference that the equation with catalytic variable to solve in the second step is now quadratic -- and can be solved with the quadratic method. At the time of writing, we have not performed the full asymptotics of moments, but from the solution it is immediate to show that $\mathbf{E}[(\tilde{Q}_n(J)-3\tilde{P}_n(J))^2]=O(n)$, which is enough to conclude~\eqref{eq:conjThirdTilde} by the Chebyshev inequality. For the impatient reader, we have already included all details in~\cite{maple}.}  that indeed, in probability, 
\begin{align}\label{eq:conjThirdTilde}
\tilde{P}_n(J) = \frac{1}{3} \tilde{Q}_n(J) + o(\tilde{Q}_n(J)).
\end{align}
	To prevent confusion, we mention that the \emph{individual} convergence of $\tilde{P}_n(J)/n^{3/4}$ and $\tilde{Q}_n(J)/n^{3/4}$  to $Z$ and $Z/3$ follow from Theorem~\ref{thm:tamari}, using a coupling between $I$ and $J$.

%	We have chosen not to include this result here (apart from this remark).
%All computer calculations have been done but details are still being written, and they will be addressed in the forthcoming longer version of this paper. In passing, we note that the (individual) convergence of $\tilde{P}_n(J)$ and $\tilde{Q}_n(J)$ to $Z$ and $Z/3$ follow from our results, as in the proof of Corollary~\ref{cor:schnyder} below. Note also that~\eqref{eq:conjThird} and~\eqref{eq:conjThirdTilde} are a priori incomparable statements, but we could deduce one from the other if we knew that the paths have a \emph{continuous} continuum limit in scale $n^{3/4}$. 
%\end{remark}

%\begin{figure}
%	\hspace{1cm}
%	\includegraphics[width=0.4\linewidth]{intervalle.png}
%	\hfill
%	\includegraphics[width=0.4\linewidth]{ratio.png}
%	\hspace{1cm}
%	\caption{Left: A  uniform random Tamari interval $(P_n,Q_n)$ of size $n=65536$ generated with a python code generously  provided by Wenjie Fang. Right: plot of $Q_n/P_n$.}\label{fig:randomPicture}
%\end{figure}




As we said, Bernardi and Bonichon~\cite{BernardiBonichon} provided an explicit bijection between intervals in $\mathcal{I}_n$ and rooted plane triangulations of size $n$. Such a triangulation can always be equipped with a \emph{canonical Schnyder wood}, which is a partition of its internal edges into three trees (say red, blue, green) with certain conditions. See Figure~\ref{fig:Schnyder}. 
%In the Bernardi-Bonichon bijection, if the interval $(P,Q)$ is in correspondence with the triangulation $T$, then $P$ encodes the blue tree in the canonical Schnyder wood of $T$. 
%Therefore we have the following (the full proof is omitted and requires to argue that the height of a random vertex and a random edge in a plane tree are closely distributed).
It is then not too hard to deduce the following from~\eqref{eq:limitlawlow}:
\begin{corollary}\label{cor:schnyder}
	Let $T_n$ be a rooted plane triangulation of size $n$, chosen uniformly at random, and let $(T^{(1)}_n, T^{(2)}_n, T^{(3)}_n)$ be its canonical Schnyder wood, that is to say, the one associated to its minimal orientation in the sense of~\cite{BernardiBonichon}. Let $V$ be a uniform random internal vertex of $T_n$ and let $H_n^{(i)}$ be the height of the vertex $V$ in the tree $T^{(i)}_n$. Then, for any $i\in [3]$ we have
	$$
	\frac{H_n^{(i)}}{n^{3/4}} \longrightarrow \frac{1}{3} Z.
	$$
%	where $Z$ is in as in Theorem~\eqref{thm:tamari}.
\end{corollary}
%\begin{proof}
%	Theorem~\ref{thm:tamariLow} asserts that $P_n(I)/n^{3/4}\rightarrow \frac{1}{2}Z^{1/4}$, and this is true also if we choose $I$ as a uniform variable on $[2n]$ instead of $[0,2n]$ (since the total variation distance between these two laws is $o(1)$ and since $i=0$ does not contributes to the height).
%	Since the bijection between Dyck paths and plane trees sends abscissas in $[2n]$ to corners of the tree, the bijection of~\cite{BernardiBonichon} would directly imply the corollary for $i=1$ if we had defined $H_n^{(1)}$ as the height of the vertex incident to a uniformly chosen \emph{corner} in the tree $T^{(1)}_n$ (rather than a uniform random \emph{internal vertex}).
%
%
%	However, on any plane tree, one can couple the uniform measure on corners and the uniform measure on internal vertices, such that the chosen corner and vertex are at distance at most 1 (explicitly, pick a uniformly random internal vertex, and consider its unique parent edge, this edge is uniformly distributed. Then flip a coin to choose one of the two halves of this edge, and choose the corner standing to the right of this half-edge around the vertex it is incident to. This corner is uniformly random).
%
%	This implies the result for $i=1$, and by symmetry for $i\in [3]$. 
%\end{proof}

\begin{figure}
\centering	\includegraphics[width=0.6\linewidth]{BBbij.png}
	\caption{A rooted planar triangulation equipped with its minimal Schnyder-wood, and its image (a Tamari interval) by the Bernardi-Bonichon bijection. The  lower path is nothing but the contour function of the blue tree. Figure taken from~\cite{BernardiBonichon}. 
}\label{fig:Schnyder}
\end{figure}



\medskip


The proof of each half of Theorem~\ref{thm:tamari} (namely \eqref{eq:limitlaw} and~\eqref{eq:limitlawlow}) 
%and~\ref{thm:tamariLow} 
has two parts: the first one consists in solving "exactly" the model, by obtaining an algebraic equation for the generating function $f(t,s)$ of intervals having a marked abscissa, with control on the size~$n$, and on the lower or upper height. 
In each case this requires to solve an equation with catalytic variables. The second part is to deduce the asymptotic of moments from there, which is a problem of analytic combinatorics in 2 variables, for which we need to find good tools. We describe below a simple method that will do the trick. 


\subsection{A method to prove moment convergence from algebraic equations}
\label{sec:trick}


Assume that we have a combinatorial class $\mathcal{A}$ equipped with a size function $|\cdot|$ and an integer statistic $\mu$, and consider the random variable $X_n=\mu(A_n)$ where $A_n$ is an object of size $n$ in $\mathcal{A}$ chosen uniformly at random. Consider the generating function
\begin{align}\label{eq:deffts}
f(t,s) = \sum_{n\geq 0} \sum_{p \in \mathbb{Z}}a_{n,p} t^n s^p
= \sum_{n\geq 0} \sum_{p \in \mathbb{Z}} a_n \mathbf{P}[X_n=p] t^n s^p,
\end{align}
with $a_{n,p} = |\{\alpha \in \mathcal{A}, |\alpha|=n, \mu(\alpha)=p\}|$, and  $a_{n} = |\{\alpha \in \mathcal{A}, |\alpha|=n\}|$,
%In this paper, $\mathcal{A}$ will be the class of Tamari intervals $(P,Q)$ equipped with a marked abscissa $i$, and $\mu(A_n)$ will be either $Q_n(i)$ or $P_n(i)-2Q_n(i)$. 
and assume that the generating function $f(t,s)$ 
%converges in a neighbourhood of $(0,0)$, and that it
is algebraic.
%, i.e. there exists a nontrivial polynomial $A$ such that 
%$$
%A(f(t,s),t,s)=0.
%$$
%
%In this case, $f(t,1)$ and all its derivatives 
For $k\geq 0$ we consider the generating function of factorial moments
$$f^{(k)}\equiv f^{(k)}(t) :=\left.\left(\left(\frac{d}{ds}\right)^k f(t,s)\right)\right|_{s=1}
% are also algebraic.
%, and standard methods in analytic combinatorics~\cite{FO, FS} enable one to obtain, in principle, the asymptotic of coefficients of $((d/ds)^kf)(t,1)$ for any $k\geq 0$. Since 
%Note that
%
%f^{(k)} 
%
=
\sum_{n\geq 0} a_n \mathbf{E}[ (X_n)_{(k)} ]t^n,
$$
with $(X_n)_{(k)}:= X_n (X_n-1)\dots (X_n-k+1)$.
To study the convergence of the random variable $X_n$, a standard way called the method of moments consists in studying the asymptotics for fixed $k$ of its $k$-th moment $\mathbf{E}[ (X_n)^k ]$ (or factorial moment $\mathbf{E}[ (X_n)_{(k)} ]$). In our setting, this is equivalent to studying the asymptotics of the coefficient $[t^n]f^{(k)}(t)$, for fixed $k$. 
Now, since all the functions $f^{(k)}(t)$ are algebraic, they are amenable to singularity analysis in the sense of~\cite{FS}. Therefore, to study the asymptotics of their coefficients it is enough to determine the nature of the dominant singularity(ies) of $f^{(k)}$ for each $k\geq 0$.
%as we will recall now.

%Let us now recall the classical transfer theorems, slightly reformulated\footnote{We will consider only the case of a unique singularity of minimal modulus. See Remark~\ref{rem:general}}. Note that we allow $c=0$ in the statement, and also $\alpha\in \mathbb{N}$ provided $c=0$. 
%\begin{proposition}[Transfer theorem for algebraic functions~\cite{FO,FS}]\label{prop:transfer}
%	Let $g(t)=\sum_{n\geq 0} g_n t^n$ be an algebraic function, with radius of convergence $\rho>0$, and assume that $t=\rho$ is the unique singularity of $g$ of minimal modulus. Assume that $g(t)$ has a singular expansion near $t=\rho$ of the form
%\begin{align}
%	g(t) \singeq c (1-t/\rho)^\alpha
%	\ \ , \ \ t\rightarrow \rho,
%	\label{eq:singg}
%\end{align}
%	with $c\in \mathbb{R}$ and $\alpha \in \mathbb{R}$, with $c=0$ if $\alpha \in \mathbb{N}$.
%	Then the coefficients $g_n=[t^n]g(t)$ satisfy the asymptotics, when $n\rightarrow \infty$,
%	$$
%	\frac{g_n}{n^{-1-\alpha}\rho^{-n}}\rightarrow \frac{c}{\Gamma(-\alpha)}  .
%	$$
%\end{proposition}
%Therefore we only need to determine the expansion~\eqref{eq:singg} for each of the functions $g=f^{(k)}$.

Here comes the main trick: since the function $f(t,s)$ is algebraic, it is also D-finite, i.e. the coefficients of its expansion in any variable satisfy a linear equation with polynomial coefficients (see e.g.~\cite{Comtet}).
%or~\cite{Stanley}).
We will apply this to the coefficients of the expansion\footnote{The function $f(t,1+r)$ is algebraic, therefore it is $D$-finite in the variable $r$. No notion of convergence is required to say this. Of course, one has to be careful about which branch of this function one considers when performing actual calculations.} in the variable $r$ such that $s=r+1$. These coefficients are, up to a factorial factor, the functions $f^{(k)}$.
%%, i.e.
%$
%[r^k] f(t,1+r) = \frac{1}{k!} f^{(k)}(t).
%$
%
%To go further, we observe that, because algebraic functions are D-finite, if the function $f(t,s)$ is algebraic,  its derivatives with respect to $s$ satisfy a linear recurrence relation with polynomial coefficients. This can be used to compute the moments recursively. While this idea is very versatile and can be adapted to various situations, we give a lemma that covers the generic case encountered in this paper.
%Therefore, we can write a recurrence formula of the form:
%\begin{align}\label{eq:Dfinite}
%	\sum_{d=0}^L   P_{d}(t,k) \frac{f^{(k-d)}(t)}{(k-d)!}=0 
%\ \ , \ \  
%k \geq L,
%\end{align}
%with $L\geq 1$ and $P_0, \dots, P_L$ being bivariate polynomials over $\mathbb{Q}$.
%In practice, it will be useful to allow ourselves to work under some change of variable (of the variable $t$), and for this reason we will more generally assume that our algebraic function $f$ satisfies~\eqref{eq:Dfinite} with $P_0, \dots, P_L$ being polynomials in $k$ whose coefficients are algebraic functions 
%having no singularity on the closed disk of radius $\rho$, except maybe at $t=\rho$
%\footnote{In our examples, they will be polynomials in $k$ and $(1-4t)^{1/4}$, with $\rho=\frac{1}{4}$.} 
%of $t$ (this will be the case in our main applications of the method).
It follows that one can compute the $f^{(k)}$ by induction, with a recurrence of the form
\begin{align}\label{eq:Dfinite2}
	f^{(k)}(t)= \sum_{d=1}^{L} h_d(t,k)  f^{(k-d)}(t) 
\ \ , \ \  
k \geq L,
\end{align}
for some $L>0$, where for each $d\in[L]$, $h_d(t,k)$ is a rational function in $k$ whose coefficients are algebraic functions of $t$ (we could assume that the $h_k$ are rational in $t$, but for applications this weaker asumption is convenient: it will enable us to work under some algebraic change of variables). This leads us to:
\begin{method}[D-finite trick for moment pumping]
	Given a bivariate  algebraic function $f$, obtain a linear equation for its derivatives $f^{(k)}$ of the form~\eqref{eq:Dfinite2} using standard computer algebra tools. Then,  use it to determine the asymptotic of $f^{(k)}$ near $t=\rho$ by induction on $k$. Deduce the asymptotics of $[t^n]f^{(k)}$ using the transfer theorem.
\end{method}

The idea of this method is quite general, but of course one has to be careful to carry the analytical details in the induction.
%(check that no unwanted singularity appears, for example; or be careful about multiple dominant singularities). 
We give below a simple framework of application   whose proof is essentially immediate.
For a function $g$ and $\alpha, c \in \mathbb{R}$, $\rho>0$ we write
$
g(t)\singeq c (1-t/\rho)^\alpha
$
if $g(t)$, as an analytic function, has no singularity on the closed circle of radius $\rho$ except maybe at $t=\rho$, and if $g(t)$ has a Puiseux expansion of the form
$
g(t) = P(t) + c (1-t/\rho)^\alpha + o((1-t/\rho)^\alpha)
$
near $t=\rho$, where $P$ is a polynomial.

%For a real number $x$, we let $\{x\}$ be the smallest \emph{nonnegative} integer which is strictly larger than $x$ (for example $\{0\}=\{0.5\}=1$, $\{-3\}=\{-0.2\}=0$).
%Note that $\{x\}=(\lfloor x\rfloor +1) \vee 0$
%and $\{x\}-1=(\lfloor x\rfloor) \vee (-1)$

\begin{theorem}[D-finite trick for moment pumping, an instance]\label{thm:trick}
	Let $f(t,s)$ be a generating function of the form~\eqref{eq:deffts}, and assume that $f(t,s)$ is an algebraic function. Then $f$ is D-finite in the variable $s-1$, and it satisfies an equation of the form~\eqref{eq:Dfinite2}, with $L\geq 1$ and where for $d\in [L]$, $h_d(t,k)$ is a rational function of $k$ whose coefficients are algebraic functions of $t$. Assume that:


	\noindent	{\bf (i)}  There is $\beta>0$ such that for $d\in [L]$ ,
	%$h_d(t,k)$ has no singularity on the closed circle of radius $\rho$, except maybe at $t=\rho$, and that $h_d(t,k)=a_d(k) (1-t/\rho)^{-\beta d}+ o((1-t/\rho)^{-\beta d})$ near this point, with possibly $a_d(k)=0$.
	$h_d(t,k)(1-t/\rho)^{\beta d} \rightarrow a_d(k) $ for $t\rightarrow \rho$, with possibly $a_d(k)=0$. Moreover, $h_d(k)$ has no singularity other than $\rho$ on the circle of radius $\rho$. 

	\noindent{\bf (ii)} "Initial conditions": there is $\alpha \in \mathbb{R}\setminus\mathbb{N}$ and numbers $c_\ell$, with $c_0>0$, such that 
	%the main singularity of $f^{(\ell)}$ is given by 
	$f^{(\ell)}\singeq c_\ell (1-t/\rho)^{\alpha-\beta\ell}$
	for all $0\leq \ell\leq \ell_0$, with $\ell_0:=L+ \max(\lfloor\alpha/\beta\rfloor,-1)$.
	Moreover, $c_{\ell}=0$ if $\alpha-\beta\ell \in \mathbb{N}$.

	Then one has 
	$f^{(k)}\singeq c_k (1-t/\rho)^{\alpha-\beta k}$
	for any $k \geq 0 $, where $c_k$ is given by the recurrence %formula
	\begin{align}\label{eq:ck}
	c_{k}= \sum_{d=1}^{L} a_d(k)  c_{k-d} \ \  ,\ \  k> \ell_0,
\end{align}
	where the values of $c_k$ for $k\leq \ell_0$ are given by the initial conditions.
Moreover, for each $k\geq 0$,
\begin{align}\label{eq:trickmoment}
	\frac{\mathbf{E}[(X_n)^k]}{n^{\beta k}} \rightarrow \frac{c_k\Gamma(-\alpha)}{c_0\Gamma(\beta k-\alpha)}.
\end{align}
\end{theorem}
The proof is ommitted in this extended abstract. It consists in showing that $f^{(k)}\singeq c_k (1-t/\rho)^{\alpha-\beta k}$ for any $k \geq 0 $, which in fact follows from a direction induction, and applying the standard theorems for algebraic functions~\cite{FS}. The reason why the recurrence starts only at $\ell_0$ is because we want $\alpha-\beta \ell$ to be negative for all $\ell$ such that $c_\ell$ appears in~\eqref{eq:ck}, so that the symbol $\singeq$ is in fact a true equivalent, in order for the induction to work.

%\begin{proof}
%	We will show by induction that $f^{(\ell)}$ has a unique singularity at $\rho$, with $f^{(\ell)}=c_\ell (1-t/\rho)^{\alpha-\beta\ell} + o((1-t/\rho)^{\alpha-\beta\ell})$, for every $\ell>\frac{\alpha}{\beta}$.
%
%	For $\frac{\alpha}{\beta}<\ell \leq \ell_0$, this is true by the initial conditions (note that $\alpha-\beta\ell<0$ for these values of $\ell$, so the symbol $\singeq$ is a true equivalent).
%	
%	Now let $k > \ell_0$ and assume the induction hypothesis for $\frac{\alpha}{\beta}<\ell<k$. Note that we have $k-d>\frac{\alpha}{\beta}$ for all $d\in [L]$, we can thus use the induction hypothesis for all the functions $f^{(k-d)}$ appearing in~\eqref{eq:Dfinite2}. Using~\eqref{eq:Dfinite2}, the hypothesis on the singularity of $h_d$, and induction, it then immediately follows that
%	$$
%	f^{(k)}=c_k (1-t/\rho)^{\alpha-\beta k} + o((1-t/\rho)^{\alpha-\beta k})
%	$$
%	with $c_k$ given by~\eqref{eq:ck}, and moreover that $f^{(k)}$ has a unique dominant singularity at $\rho$.
%This concludes the induction, and shows that we have $f^{(\ell)}\singeq c_\ell (1-t/\rho)^{\alpha-\beta\ell}$ for all $\ell\geq0$.
%
%	By the transfer theorem, this implies that $[t^n]f^{(\ell)} \sim \frac{c_\ell}{\Gamma(\beta\ell-\alpha)}n^{\beta\ell-\alpha-1}\rho^{-n}$.
%	Since $\mathbf{E}[(X_n)_{(k)}]=\frac{[t^n]f^{(k)}}{[t^n]f^{(0)}}$, the asymptotics~\eqref{eq:trickmoment} follows for $\mathbf{E}[(X_n)_{(k)}]$ by the transfer theorem. The same asymptotics holds for $\mathbf{E}[(X_n)^{k}]$ by the triangular relation between moments and factorial moments.
%\end{proof}
%
%
%
\begin{remark}\label{rem:general}
	We do not try to provide minimal hypotheses in this theorem:
	%our goal is to emphasize the method, that seems applicable to various situations, including some D-finite non algebraic functions. T
	the only requirement for the method to work is that dominant singularities of $f^{(k)}$ are "not too hard" to track by induction from~\eqref{eq:Dfinite2}.
Also, it is conceivable to use this technique in more than bivariate examples,
%	Note also that there is  no direct objection to trying to apply the method for more than bi-variate examples,
%	i.e. for the asymptotic of joint moments of variables of the form $\mathbf{E}X_{n,1}^{k_1}\dots X_{n,j}^{k_j}$ for fixed $k_1,\dots,k_j$ when $n\rightarrow \infty$,
provided the corresponding multivariate generating function is algebraic.% (see also Section~\ref{sec:Rayleigh}).
%\end{remark}
%\begin{remark}

	We insist that Theorem~\ref{thm:trick} if applicable, is essentially automatic. Indeed, computer algebra softwares (e.g. the Maple package GFUN~\cite{gfun}) are able to provide a recurrence of the form~\eqref{eq:Dfinite2}, and to check the initial conditions, automatically from an algebraic equation for $f$.
	Apart from this, the major interest of this method is that it allows a wide variety of limit laws (including non Gaussian) as our main applications show.
	%Also, since all the $f^{(k)}$ are algebraic with a computable algebraic equation, it is automatic to check the initial condition required in the lemma for $0\leq \ell \leq \ell_0$. All this will be illustrated by the proofs of Theorems~\ref{thm:tamari}-\ref{thm:tamariLow}.

To conclude this discussion, let us recall that using D-finiteness to compute coefficients of algebraic functions in the univariate case is a well-known trick going back at least to Comtet~\cite{Comtet}.
%and is extensively used in the context of computer algebra. Here 
We are only recycling this idea in the context of  bivariate asymptotics.
\end{remark}

	\medskip

{\bf Example.} As a simple application, we invite our reader to rediscover the well known limit law for the height $H_n$ of a uniform random point on a uniform random Dyck path of size $n$.
Using standard path decompositions one easily obtains a quadratic equation for the corresponding generating function $f(t,s)$, and from there it is immediate, with Maple, to check that Theorem~\ref{thm:trick} applies (with $-\alpha=\beta=\tfrac{1}{2}$). The recurrence~\eqref{eq:ck} becomes $c_k=\frac{1}{4} k (k-1) c_{k-2}$, and one thus directly shows that $H_n/\sqrt{n}$ converges to a Rayleigh law, of $k$-th moment $\Gamma(\tfrac{k}{2}+1)$. 
%
%
%B%y standard path decompositions the corresponding generating functions satisfies
%
%%\begin{example}\label{eq:Dyck}
%	As a simple (and somewhat artificial) application, let us rediscover the well-known limit law for the height $H_n$ of a uniform random point on a uniform random Dyck path of size $n$. Using standard methods (e.g. last passage decompositions), the corresponding generating function is found to be
%$$
%	f(t,s)=\frac{E(t)^2}{1-st E(t)^2}\ \ , \ \ E(t)=1+tE(t)^2.
%$$
%Using a computer algebra system, it is immediate to eliminate $E(t)$ to find an algebraic equation for $f(t,s)$. %We then find automatically a linear recurrence equation for its derivatives at $s=1$.
%	and the Maple package gfun~\cite{gfun} gives us an equation of the form~\eqref{eq:Dfinite2} with $L=2$.
%\begin{align}\label{eq:rayleigh}
%	\frac{f^{(k)}(t)}{k!} = -\frac{f^{(k-1)}}{(k-1)!}  +\frac{t}{1-4t}\frac{f^{(k-2)}}{(k-2)!}=0 \ \ , k\geq 2.
%\end{align}
	%it is then immediate to prove by induction that $m_{k}(t)$ has a dominant singularity a $t=\frac{1}{4}$ of the form $m_k(t)\sim 2^k (1-4t)^{k/2}$ (note that the middle term in~\eqref{eq:rayleigh} do not contribute at the dominant order). Using transfer theorems, one deduces that $[t^n]m_k(t)\sim TODO$, which implies that the height $H_n$ of a point taken uniformly at random on a Dyck path of size $n$ satisfies
%	This recurrence (multiplied by $k!$) fits in the framework of the theorem with $L=2$, $a_1=0$, $a_2=\frac{1}{4}k(k-1)$, and $\beta=1/2$. It is then easy to check explicitely the initial conditions, which hold with $\alpha=-\frac{1}{2}$, thus $\ell_0=1$, and with $c_0=2$ and $c_1=1$.
%	It is then a direct check that Theorem~\ref{thm:trick} applies with $\alpha=\beta=1/2, \ell_0=1,$ and the recurrence~\eqref{eq:ck} becomes $c_k=\frac{1}{4} k (k-1) c_{k-2}$.
	%and initial conditions $c_0=2$ and $c_1=1$.
	%(note that the middle-term in~\eqref{eq:rayleigh} does not contribute, in the notation of the lemma we have $a_1=0$), 
%	(so that $c_k=2^{1-k}k!$)
%Once initial conditions are checked this recurrence is obviously solved, and~\eqref{eq:trickmoment} immediately gives that $\frac{\mathbf{E}[(H_n)^k]}{n^{k/2}}$ converges to
%	%$\frac{2^{-k}k!\Gamma(\frac{1}{2})}{\Gamma(\frac{k}{2}+\frac{1}{2})}=
%	$\Gamma\left(\frac{k}{2}+1\right),$
%	This direction implies that $H_n/\sqrt{n}$ converges in distribution to a random variable,
%	which we recognize as the $k$-th moment of a Rayleigh law of parameter $\sigma=1/\sqrt{2}$.
	%(with density $\frac{x}{\sigma^2}e^{-x^2/(2\sigma^2)}$ on $(0,\infty)$). 
	All calculations and checks are available in~\cite{maple}!


%Note that in principle, the method could work for more marked points. Consider a uniform random Dyck path of size $n$, and let $X_n,Y_n$, be the heights of two marked abscissas.

%	In passing, we remark that we can also apply the theorem with a value of $\beta$ larger than $1/2$. In this case get $a_d\equiv 0$, and $c_k=0$ for $k>0$. We conclude that all positive of moments of $H_n/n^{\beta}$ go to zero -- which is true but weaker than the statement we obtained with the optimal value $\beta=1/2$.
%\end{example}



%\subsection{Comments, questions, and plan of the paper}




%When we brought Theorem~\ref{thm:tamari} to the attention of Nicolas Curien, he mentionned a paper in preparation describing the scaling limit, in scale $n^{3/4}$, of the  parking trees of~\cite{parkingTrees}. Critical parking tree models are in the same universality class as Tamari intervals at least from the enumerative point of view, with asymptotic enumeration exponent $n^{-5/2}$ and a  description with an equation with one-catalytic variable. As far as we understand, the methods of that paper will apply similarly to Tamari intervals and describe the scaling limit of the (whole) upper path $Q$. We do not know if these methods will also provide alternative proofs of our results. We suspect however that, all calculations made, the methods of our paper should enable to obtain the limit law of a uniform random point in critical random  parking trees (but have not tried to do so).
%
\medskip


{\bf Plan of the paper.}
%The paper is organized as follows.
In Section~\ref{sec:tamari} we solve the exact counting problem which underlies the first half of Theorem~\ref{thm:tamari} (Equation~\eqref{eq:limitlaw} by studying the classical  equation with one catalytic variable for Tamari intervals, which we enrich by an extra variable marking the height. See Theorem~\ref{thm:exactUp}.
In Section~\ref{sec:tamariLow} we solve the exact counting problem which underlies the second half (Equation~\eqref{eq:limitlawlow})
%Theorem~\ref{thm:tamariLow} which deals with the combined height of the upper and lower path. See
see Theorem~\ref{thm:exactDown}. Our proof is more technical, as we only manage to write an equation with \emph{two} catalytic variables for this problem (enriched, once again, by an extra variable for the height). 
%It turns out that this equation with two catalytic variables is not too hard to solve, as it can be viewed as two nested equations with one catalytic variable each, which can be solved by standard methods.
In Section~\ref{sec:asymptotics}, we sketch the proof of Theorem~\ref{thm:tamari}, which given these results directly follows from 
%these exact calculations and 
%by performing the asymptotics of moments from the algebraic equations obtained in Theorems~\ref{thm:exactUp} and~\ref{thm:exactDown}. The proofs are  immediate verifications given
Theorem~\ref{thm:trick} and computer algebra calculations.

%Finally in Section~\ref{sec:other}, we prove Theorems~\ref{thm:parking} and~\ref{thm:tutte} by applying the same method as in Section~\ref{sec:tamari} to the equations for these models.
%
%\smallskip

%We will try in this paper to focus on the methods, and leave the actual 
All the calculations supporting our results (including results whose proofs are not fully presented in this abstract) are available in the accompanying Maple worksheet~\cite{maple}.
%Whether these calculations could be done or verified by hand will not be our concern.

\section{Upper path: exact solution}
\label{sec:tamari}

\subsection{The classical equation and its solution, after \cite{Chapoton, BMFPR}}
\label{sec:decomp}

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{decomp1.pdf}
	\caption{The classical decomposition of Tamari intervals. To the left, an interval of size $n+1$, where $v_1,v_2$ are the first contacts of the lower and upper path, respectively. The decomposition gives rises, to the right, to two Tamari intervals of total size $n$, the first of which has a marked contact, called here $\tilde{v}$. This construction is bijective.}\label{fig:decomp1}
\end{figure}

%Let $(P,Q)$ be a Tamari interval. 
Following~\cite{BMFPR}, we call \emph{contact} of a path a vertex of that path lying on the $x$-axis. We now present a recursive decomposition of Tamari intervals based on contacts of the lower path, following~\cite{Chapoton, BMFPR}.
%See Figure~\ref{fig:decomp1} as a support for this discussion.

%This decomposition is the starting point to prove Theorems~\ref{thm:tamari}-\ref{thm:tamariLow}.

Let $(P,Q)\in \mathcal{I}_{n+1}$ be a Tamari interval, for $n\geq 0$, and let $v_1$ and $v_2$ be respectively the leftmost contact of $P$ and $Q$ different from the origin.
By deleting the first up step of $P$ and $Q$, and the downstep of $P$ and $Q$ preceding respectively $v_1$ and $v_2$, one obtains two paths that can be naturally seen as the concatenation of two pairs of paths $(P_1,Q_1)$ and $(P_2,Q_2)$ as on Figure~\ref{fig:decomp1}.
%
%Then $v_2$ is also a contact of $P$ and it is preceded by a down-step in both $P$ and $Q$. We can thus split paths as $P=uP_0dP_2$, $Q=uQ_1dQ_2$, where we write paths as words with letters $u/d$ representing up/down steps respectively, and where $P_2,Q_2$ start at the contact $v_2$. Moreover, $v_1$ is a contact of $uP_0d$ (possibly equal to $v_2$), so we can write $uP_0d=uP'_1dP''_1$, with $P''_1$ possibly empty, and $P_1''$ starting at $v_1$.
%We now define the path $P_1:=P_1'P_1''$, which is naturally equipped with a marked contact~$\tilde{v}$. See Figure~\ref{fig:decomp1} again.
It is proved in~\cite{Chapoton, BMFPR} that $(P_1,Q_1)$ and $(P_2,Q_2)$ are two Tamari intervals. Moreover, this operation is a \emph{bijection} (!) between intervals of size $n+1$ and pairs of intervals of total size $n$, such that the first interval of the pair has a distinguished contact on its lower path -- inherited from the vertex $v_1$.
To translate this recursive construction into an equation for enumeration, one needs to introduce a two-parameter generating function:% of Tamari intervals:
$$
F(x)\equiv F(t,x) := \sum_{n\geq 0} t^n \sum_{(P,Q)\in \mathcal{I}_n} x^{\mathrm{contact}(P)}
$$
where $\mathrm{contact}(P)$ is the number of contacts of the lower path $P$.


Note that, if a path $\tilde{P}_1$ has $k$ contacts, there are $k$ possible ways to mark a contact in $\tilde{P}_1$, thus decomposing it as $\tilde{P}_1=\tilde{P}_1'\tilde{P}_1''$. When going through these $k$ choices, the number of \emph{non-final} contacts of the path $\tilde{P}_0:=u\tilde{P}_1'd\tilde{P}''_1$ goes through the values $1, 2, \dots, k$, see Figure~\ref{fig:operator2}(up).
At the level of generating functions, 
%the operation of marking a contact in an interval $(\tilde{P}_1, \tilde{Q}_1)$ and transforming into an interval of the form $(\tilde{P},\tilde{Q})$ with this construction is taken into account by the operator:
the corresponding operator is% given by
\begin{align}\label{eq:operator}
	&\Delta_x: x^k\hspace{3mm} \longmapsto x^k+\dots +x^{1} = x \frac{x^k-1}{x-1} \ \ , \ \ 
	\Delta_x: f(x) \longmapsto x \frac{f(x)-1}{x-1}.
\end{align}
%The reason why we choose to count only non-final contacts is to avoid counting the final contact twice when concatenating intervals (see Figure~\ref{fig:decomp1} again).
%
%Note that, in the construction above, the number of contacts of $P$ is equal to the number of contacts of $P_2$ plus the number of contacts of $P_1$ that follow its marked contact. 
%In general, if a path $P_1$ has $k$ contacts, there are $k$ ways to mark a contact. 
%
This observation being made, the recursive decomposition immediately translates into the following functional equation~\cite{Chapoton, BMFPR}
\begin{align}\label{eq:catalyticTamariSimple}
	F(x) = x+xtF(x)\frac{F(x)-F(1)}{x-1}.
\end{align}
Indeed, the first term accounts for the empty path (of size $n=0$), and in the second term the factor $F(x)$ is the contribution of the interval $(P_2,Q_2)$, and the factor $x\frac{F(x)-F(1)}{x-1}=\Delta_x F(x)$ is the contribution of the interval $(P_1, Q_1)$.

%\begin{figure}
%	\centering
%	\includegraphics[width=\linewidth]{operator.pdf}
%%	\caption{How the divided difference operator appears. On the left, the power of $x$ marks all contacts, while on the right it only marks contacts which are not the last one.}\label{fig:operator}
%\end{figure}

%Equation~\eqref{eq:catalyticTamariSimple} is a typical example of a \emph{functional equation with one-catalytic variable}, which is the name given in the combinatorics literature\footnote{We have been told that a part of the community is shifting towards the fancier terminology "discrete differential equation".%, but let's take our time and be old-fashioned.
%} to polynomial equations involving a bivariate series $F(t,x)$ as well as its specialisation at $x=1$. Such equations are not easy to solve \emph{a priori}, but fortunately a full theory has been developed in~\cite{BMJ}, which makes their resolution essentially automatic. In the present case, 
The solution can be written especially nicely with a rational parametrization.
\begin{proposition}[{\cite[Thm. 10]{BMFPR}}]
	The functions $F(x)$ and $F(1)$ are given by
	%have the explicit rational parametrisations:
\begin{align}\label{eq:solutionTamariSimple}
	F(x)=\frac{1+u}{(1+zu)(1-z)^3}(1-2z-z^2u)
	\  \  \ , \ \ \ F(1)=\frac{1-2z}{(1-z)^3},
\end{align}
with 
\begin{align}\label{eq:z} %\label{eq:u}
t=z(1-z)^3  \ \ , \ \   
	x=\frac{1+u}{(1+zu)^2}.
\end{align}
\end{proposition}
%Of course, this explicit expression of $F(1)$ directly implies the formula~\eqref{eq:Chapoton}, and this is historically the main reason why the equations of this section were written and solved in the first place.


\subsection{The enriched equation and its solution}

%One can apply the same decomposition as above to Tamari intervals in which the upper path $Q$ carries a marked point, and track where this marked point appears through the decomposition.
%: one only has to track recursively the position of the marked point into the different components of the decomposition.
We now study intervals in which the upper path carries a marked point. We let
%. We introduce the generating function
$$
H(x,s)\equiv H(t,x,s) := \sum_{n\geq 0} t^n \sum_{(P,Q)\in \mathcal{I}_n} x^{\mathrm{contact}(P)} \sum_{i=0}^{2n} s^{Q(i)}.
$$
%where by $Q(i)$ we mean the height of the point of the path $Q$ lying at abscissa $i$.
It is clear that the height of the marked point can be tracked in the decomposition above. This leads to the functional equation:
\begin{proposition}\label{prop:exactUp}
	The generating function $H(x,s)$
	%of Tamari intervals with a marked abscissa where $s$ marks the upper height 
	is solution of the equation:
\begin{align}\label{eq:catalyticTamariUpper}
	H(x,s) =F(x)+sxt\frac{H(x,s)-H(1,s)}{x-1}F(x)+xt\frac{F(x)-F(1)}{x-1}H(x,s).
\end{align}
\end{proposition}
\begin{proof}
	This follows directly from the combinatorial decomposition of Figure~\ref{fig:decomp1}, applied to intervals with a mark abscissa.
	%, see Figure~\ref{fig:decomp11}.
	The first term accounts for the case where the marked abscissa (and height) is zero. The second term accounts for the case where the marked abscissa appears before vertex $v_2$. 
	%In this case, the height of the marked vertex of the upper path is shifted by one in the decomposition, henve the factor of $s$.
 Through the decomposition the corresponding vertex of the upper path becomes a marked vertex of the path $Q_1$, and its height is shifted by~$1$, hence a contribution of $s\Delta_x H(s,x)$ for the interval $(P_1,Q_1)$,
while the interval $(P_2,Q_2)$ has contribution $F(x)$ as before. 
	The third term accounts for the case where the marked abscissa appears at $v_2$ of after, in which case there is no shift in the height and the second interval has contribution $H(s,x)$, while the first has contribution $\Delta_x F(x)$ as before.
	%, in which case the corresponding vertex becomes a marked vertex of the path $Q_2$, with no shift in the height, hence a contribution of $\Delta_xF(x)$ for $(P_1,Q_1)$ and $H(x,s)$ for $(P_2,Q_2)$.
\end{proof}

%\begin{figure}
%	\centering
%	\includegraphics[width=\linewidth]{decomp11.pdf}
%	\caption{The decomposition of Figure~\ref{fig:decomp1}, now with a marked abscissa. The three cases correspond to the three terms in Equation~\eqref{eq:catalyticTamariUpper}. In cases 2 and 3, the upper point at this abscissa becomes an upper point of the path $Q_1$ and $Q_2$ respectively, and in case 2 the height is shifted by one.}\label{fig:decomp11}
%\end{figure}
%Our most skeptical readers might convince themselves (again) that we haven't forgotten any terms in the equation by noting that $H(x,1)=(2\frac{xd}{dx}+1)F(x)$, since there are $2n+1$ possible abscissas to mark in an interval of size $n$, and that Equation~\eqref{eq:catalyticTamariUpper} for $s=1$ can be obtained by applying the operator $2\frac{xd}{dx}+1$ to~\eqref{eq:catalyticTamariSimple}.

%\medskip
%
%Note that, since $F(x)$ and $F(1)$ are known, this equation for $H(x)\equiv H(x,s)$ is again a polynomial equation with catalytic variable $x$, in which the variable $s$ only plays the role of a parameter. In some sense, this equation is simpler than the previous one since it is \emph{linear} in $H$. Linear equations with one catalytic variable are easy to solve, at least in principle, via the \emph{kernel method}, which we now apply.

Equation~\eqref{eq:catalyticTamariUpper} is easily solved via the kernel method (see e.g.~\cite[p.~508]{FS}).
First, let us work under the change of variables $(t,x)\leftrightarrow(u,z)$ of~\eqref{eq:z}, so we can consider $F(1)$ and $F(x)$ known.
%in kernel form as follows:
We can then write~\eqref{eq:catalyticTamariUpper} in "kernel" form
%\begin{align}\label{eq:kernelForm}
	$	K(x,s) H(x,s) = R.H.S.$
	%F(x)-\frac{sxtF(x)H(1)}{(x-1)}
%\end{align}
%with 
%\begin{align}\label{eq:kernel}
%	K(x,s)=\left(-1+sxt\frac{F(x)}{x-1}+xt\frac{F(x)-F(1)}{x-1}\right).
%\end{align}
%it is easily checked that there is a unique formal power series $X_0(t)$
%\in \mathbb{Q}[s][[t]]$ 
%that is such that 
%$
%K(X_0,s)=0.
%$
%Indeed, multiplying by $x-1$, this equation has the form 
%$X_0(t)=1+t\text{Polynomial}(X_0(t),t,s)$ which enables one to compute the coefficients of $X_0(t)$ recursively. 
%We will now apply the kernel method. The most obvious way to do it consists in substituting $x$ for $X_0(t)$ and solving for $H(1)$,
%This enables one to obtain explicitly the value of $H(1)$ by substituting $x=X$ in~\eqref{eq:catalyticTamariUpper} and solving. We get in this way an explicit expression of $H(1)$ in terms of $X,s$, and $F(X)$, which is satisfying since the function $F(x)$ is known from the previous section.
%%but we find more convenient to work directly with the variables $(u,z)$ rather than $x,t$.
where the kernel $K(x,s)$ is an explicit rational function of $z$ and $u$. One easily checks that there is a unique series $U=U(z)$ cancelling the kernel, given by
%Expressing the kernel in terms of these variables we have (see~\cite{maple} for the calculation)
%$$
%\tilde{K}\equiv \tilde{K}(z,u,s) := K(x,s,t)\Big|_{t=z(1-z)^3\atop x=\frac{1+u}{(1-zu)^2}}=\frac{P(u,z,s)}{(uz+1)u(uz^2+2z-1)},
%$$
%with 
%$$
%P(u,z,s)=s(1+u)^2z(uz^2+2z-1)-uz^3+3uz^2-3zu+u.
%$$
%Since $P(u,z,s)=u-sz + O(z^2+u^2+zu)$, there is a unique series $U\equiv U(z,s) \in \mathbb{Q}[s][[z]]$ such that $P(U,z,s)=0$. Its expansion starts as 
%$$U=sz + O(z^2).$$ 
%Moreover, since $P$ has degree $1$ in $s$, we can also write $s$ as a rational function of $U$ and $z$,
\begin{align}\label{eq:srat}
	s={U(1-z)^3}/[{z(1+U)^2(1-Uz^2-2z)}].
\end{align}
Solving the R.H.S. for $H(1)$ we obtain (see~\cite{maple} for full calculations and checks)
\begin{theorem}\label{thm:exactUp}
	The series $H(1)\equiv H(t,s,1)$ of Tamari intervals with a marked abscissa, where $t$ marks the size, and $s$ the upper height at this abscissa, has the following rational parametrisation:
	\begin{align}\label{eq:Hsexplicit}
		H(1)= {(1-2z-Uz^2)^2(1+U)}/{(1-z)^6}
\end{align}
	with the change of variables $(t,s)\leftrightarrow (z,U)$ given by $t=z(1-z)^3$ and by~\eqref{eq:srat}.
\end{theorem}
%\begin{proof}
%	We substitute $t=z(1-z)^3$ and $x=\frac{1+u}{(1-zu)^2}$ in~\eqref{eq:kernelForm}, and we further substitute
%	$u=U(z,s)$, which cancels the left-hand side. From the right-hand side, and using the explicit expression~\eqref{eq:solutionTamariSimple} of $F(1)$ in terms of $z$  we thus get an expression for $H(1)$ as a rational function of $z$, $s$, and $U=U(z,s)$, in which $s$ can be substituted by~\eqref{eq:srat}, thus giving an expression for $H(1)$ as a rational function of $U$ and $z$. The only thing to be done is to perform the explicit calculation. See~\cite{maple}.
%\end{proof}
%

\section{Lower path: exact solution}
\label{sec:tamariLow}

We will now apply the same decomposition as in the previous sections, but keep track of the height of points on the lower path $P$. In order to do this, we will have to treat differently the contacts of $P$ which appear before or after the marked abscissa, which will force us to work with \emph{two} catalytic variables.
We write $\mathrm{contact}_{< i}(P), \mathrm{contact}_{\geq i}(P)$ for the number of contacts of $P$ strictly before, or weakly after, abscissa $i$, respectively.
We introduce the generating function
$$
G(x,y)\equiv G(t,x,y,w) := \sum_{n\geq 0} t^n \sum_{(P,Q)\in \mathcal{I}_n}  \sum_{i=0}^{2n}
%s^{Q(i)} 
w^{P(i)}x^{\mathrm{contact}_{< i}(P)}y^{\mathrm{contact}_{\geq i}(P)}.
$$

\subsection{The enriched equation}


%In order to prove our wanted asymptotic theorem, we will only need to work under the specialisation $s=w, r=w^{-2}$, but for the sake of clarity in writing the next equation, we prefer to work with the two variables $r$ and $s$. 
%We write $H_1(x):=H(x,s)_{|s=1}$, this series has been computed in the previous section.
%We have
\begin{proposition}
	The generating function $G(x,y)\equiv G(t,x,y,w)$ of Tamari intervals with a marked abscissa where $w$ marks the lower height is solution of the equation:
	\begin{align}
	G(x,y) &=
	F(y)
+
	txw\frac{G(1,y)-G(1,1)}{y-1}
	F(y)
+
tx \frac{F(y)-yF(1)}{y-1} F(y) \nonumber \\
	&+
	t\frac{x^2}{y} \frac{G(x,y)-\frac{y}{x}F(x)-G(1,y)+yF(1)}{x-1} F(y) 
		+tx\frac{F(x)-F(1)}{x-1} G(x,y).\label{eq:catalyticTamariLower}
\end{align}
\end{proposition}
%\begin{figure}
%	\centering
%	\includegraphics[width=\linewidth]{decomp12.pdf}
%	\caption{The decomposition of Figure~\ref{fig:decomp1}, with a marked abscissa, but we now track the height of points of both paths. There are now five cases, which correspond to the four terms in Equation~\eqref{eq:catalyticTamariLower}.}\label{fig:decomp12}
%\end{figure}

\begin{proof}[Idea of the proof]
Given an interval $(P,Q)$ of size $n+1$ with a marked abscissa, we apply again the decomposition of Figure~\ref{fig:decomp1}. 
%	As before we let $v_1$, $v_2$, be the first non-initial contacts of the paths $P$ and $Q$, respectively, and we also let $i_1, i_2$ be their abscissa. We let $i\in [0,2n]$ be the marked abscissa.
	We let $i_1,i_2$ be the abscissa of the first nonzero contact of $P$ and $Q$ respectively, and $i$ the marked abscissa. We will distinguish fives cases depending on the fact that $i$ belongs to $\{0\}, [1,i_1-1], \{i_1\}, [i_1+1,i_2-1],[i_2,2n]$. These correspond (from left to right) to the five terms in~\eqref{eq:catalyticTamariLower}.

	In this extended abstract, we will only address the second case, which illustrates in a simple way why we  need two catalytic variables.
%
%	We distinguish five cases (we will treat the third one last), see Figure~\ref{fig:decomp12},
%\begin{itemize}
%	\item Case 1: $i=0$. Then we are only  considering a Tamari interval, where all contacts are marked by the variable $y$. The corresponding contribution if $F(y)$.
%	\item Case 2: 
		If $i \in [1,i_1-1]$, in the decomposition of Figure~\ref{fig:decomp1}, the corresponding vertex of $P$ becomes a marked vertex of $P_1$ with a shift of $1$ in the height, hence a contribution of $w$. 
%		Moreover, the marked contact of the path $P_1$ after applying the construction appears after these marked vertices. Therefore, 
		Moreover, configurations in this case are obtained by applying the construction  of Figure~\ref{fig:operator2}(top) but restricting it to contacts appearing after the marked abscissa, see Figure~\ref{fig:operator2}(center). Therefore, an interval $(\tilde{P_1}, \tilde{Q_1})$ having a marked abscissa, with respectively $k$ and $\ell$ contacts strictly before, and weakly after this abscissa (thus having a contribution of $x^ky^\ell$ in $G(x,y)$) gives rise to $\ell$ intervals contributing to this case, with a contribution of 
		$
		x^1 y^{\ell-1} + x^1 y^{\ell-2} + \dots + x^1 y^0 = x \frac{y^\ell-1}{y-1}.
		$
In total, the contribution for the first interval is thus $x \frac{G(1,y)-G(1,1)}{y-1}$. The contribution of the second interval $(P_2,Q_2)$ is just $F(y)$, since all corresponding contacts appear after the marked abscissa. In total, this 
		%cases gives a contribution of 
	%	$t x w \frac{G(1,y)-G(1,1)}{y-1} F(y)$, which is
		gives the second term in~\eqref{eq:catalyticTamariLower}.

		We omit other cases, but we point out that the fourth case requires a similar discussion where now the catalytic variable $x$ plays the main role, see Figure~\ref{fig:operator2}(bottom). Some care is needed  since the last vertex of the path $\tilde{Q}_1$ cannot be marked in this case, hence the slightly more complicated numerator in this term.
%	\item  Case 4: $i \in (i_1,i_2)$. In the decomposition, the vertex of $P$ becomes a vertex of $P_1$ with no shift in height. 
%		Moreover, these configurations are obtained by applying the construction  of Figure~\ref{fig:operator} but restricting it to contacts appearing before the marked abscissa, see Figure~\ref{fig:operator2}(down). More precisely, an interval $(\tilde{P_1}, \tilde{Q_1})$ having a marked abscissa which is not the final one, with respectively $k$ and $\ell$ contacts stricly before, and weakly after, this abscissa (thus having a contribution of $x^ky^\ell$ in $G(x,y)$) gives rise to $k$ intervals contributing to this case, with a contribution of 
%		$$
%		x^{k+1} y^{\ell} + x^{k-1} y^{\ell-1} + \dots + x^2 y^{\ell-1} = \frac{x^2}{y}\frac{x^k-1}{x-1} y^\ell.
%		$$
%In total, the contribution for the first interval is thus 
%		$$t\frac{x^2}{y} \frac{\tilde{G}(x,y)-\tilde{G}(1,y)}{x-1},$$
%		where $\tilde{G}$ is the same generating function as $G$, but in which the marked abscissa is not the last one.
%We have
%$$\tilde{G}(x,y) = G(x,y) - \frac{y}{x}F(x),
%$$
%since an interval with the last abscissa marked is the same as an unmarked interval (and the corresponding vertex has height zero), and since in such a case all contacts are counted with weight $x$ except the last one which receives a weight $y$.
%
%		Moreover, the contribution of the second interval $(P_2,Q_2)$ is just $F(y)$, since all corresponding contacts appear after the marked abscissa. In total, this case gives a contribution of 
%		$t\frac{x^2}{y} \frac{\tilde{G}(x,y)-\tilde{G}(1,y)}{x-1} F(y)$, wh
%		which is the fourth term in~\eqref{eq:catalyticTamariLower}.
%	\item  Case 5: $i\geq i_2$. In this case, the marked vertex becomes a marked vertex of $Q_2$ and there is not shift in  height. The contribution of $(P_2,Q_2)$ is simply $G(x,y)$, while the contribution of the first interval $(P_1,Q_1)$ is $xt \frac{F(x)-F(1)}{x-1}$ as in the classical case, since all its contacts appear before the marked abscissa. This gives the last term in~\eqref{eq:catalyticTamariLower}.
%	\item Case 3: $i=i_1$ and $i_1\not\in\{0,i_2\}$. Such configurations are made of the concatenation of two intervals.
%		The second one has no constraint, hence a contribution of $F(y)$. The first one has several constraints. First, its upper path has only two contacts. Arguing as in the classical decomposition, the generating function for intervals where the upper path has two contacts is 
%$$ ty\frac{F(y)-F(1)}{y-1}.
%$$
%		Second, in the first interval, the first nonzero contact of the lower path should be different from the one of the upper path (since we want $i_1\neq i_2$ in the end). We thus have to substract the corresponding contribution, and in total we find that the contribution of the first interval is 
%		$$
%		ty\frac{F(y)-F(1)}{y-1} - ty F(1).
%		$$
%		Putting things together, the  generating function Case 3 is 
%		$$
%		\frac{x}{y} \left(ty\frac{F(y)-F(1)}{y-1} -tyF(1)\right) F(y),
%		$$
%		where the prefactor $\frac{x}{y}$ accounts for the initial contact. 
%This is the third term in~\eqref{eq:catalyticTamariLower}.
%\end{itemize}
\end{proof}
%\noindent Note: We have tested each term of~\eqref{eq:catalyticTamariLower} up to order $n=8$, 
%against exhaustive generation.
%using an exhaustive generation of Tamari lattices made by independent means.
	%(we say this just in case the reader is worried by the many occasions to miss a detail in the five cases of the proof), and in fact we have tested the five cases independently.

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{operator2.pdf}
	\caption{
		Top: How the divided difference operator appears. On the left, the power of $x$ marks all contacts, while on the right it only marks contacts which are not the last one.
	Center and Bottom:	Refinement to distinguish the case where the marked abscissa is before (upper figure) or after (lower figure) the contact $v_1$.}\label{fig:operator2}
\end{figure}

%For the rest of this section we will work with the function 
%$$
%G(x,y)\equiv G(x,y,w)=D(x,y,w^{-2},w) = \sum_{n\geq 0} t^n \sum_{(P,Q)\in \mathcal{I}_n}  \sum_{i=0}^{2n} w^{P(i)-2Q(i)}x^{\mathrm{contact}_{\leq i}(P)}y^{\mathrm{contact}_{>i}(P)}.
%$$
%From the last proposition, it is solution of the equation:
%\begin{align}\label{eq:catalyticTamariLowerW}
%	G(x,y) =
%	F(y)+
%	w^{-1}
%	tx\frac{G(1,y)-G(1,1)}{y-1}
%	F(y)
%+
%	txw
%	\frac{G(x,y)-G(1,y)}{x-1}
%	F(y) 
%		+xt\frac{F(x)-F(1)}{x-1} G(x,y).
%\end{align}
%We now have to solve this equation.

\subsection{Solution}
\label{sec:solutionLow}

Although equations with two catalytic variables are notoriously difficult, Equation~\eqref{eq:catalyticTamariLower} is of a very particular kind, 
%First, it is \emph{linear} in the main unknown $G(x,y)$. Second, 
as it involves $G(x,y)$, $G(1,y)$, and $G(1,1)$, but not $G(x,1)$. This will enable us to treat~\eqref{eq:catalyticTamariLower} as two nested equations, each having only one catalytic variable.
In what follows we sketch the resolution, see~\cite{maple} for more details.
%: we will first eliminate the variable $x$, then $y$, using standard methods, in two steps.
%The reader interested in more difficult cases of equations with two catalytic variables might consult, for example, the reference~\cite{BBMR} for entry points into this fastly growing literature.

{\noindent \bf First step: eliminating the variable $x$ (or $u$).}
%We start by writing Equation~\eqref{eq:catalyticTamariLower} in kernel form
%\begin{align}\label{eq:kernelForm2}
%	K_2(x,y) G(x,y)=&
%	F(y)
%	+
%	txw\frac{G(1,y)-G(1,1)}{y-1}
%	F(y)\nonumber \\
%	&+
%tx \frac{F(y)-yF(1)}{y-1} F(y) 
%	+
%	t\frac{x^2}{y} \frac{-\frac{y}{x}F(x)-G(1,y)+yF(1)}{x-1} F(y) .
%	%tx\frac{G(1,y)-G(1,1)}{y-1}F(y)w^{-1}+F(y)-tx\frac{G(1,y)F(y)w}{x-1},
%\end{align}
%with
%\begin{align}\label{eq:kernel2}
%	K_2(x,y) =\frac{tx^2F(y)}{y(x-1)}-1+xt\frac{F(x)-F(1)}{x-1}.
%\end{align}
We will work under the change of variables~\eqref{eq:z}, with a new variable $v$ which is to $y$ what $u$ is to $x$, namely
\begin{align}\label{eq:zuv}
t=z(1-z)^3 \ \  , \ \   
	x=\frac{1+u}{(1+zu)^2} \ \ , \ \ 
	y=\frac{1+v}{(1+zv)^2}.
\end{align}
We write respectively $\tilde{G}(u,v), \tilde{G}_1(v), \tilde{G}_{11}$ for the quantities $G(x,y),G(1,y),G(1,1)$ expressed in the variables $z,u,v$ after the substitutions~\eqref{eq:zuv}.


%Its numerator is a polynomial $\tilde{k}_2(z,u,v,w)$ which, computations made~\cite{maple}, is given by 
%\begin{align*}
%	\tilde{k}_2(z,u,v,w)=
%	((vz^2+z^2)u+vz^2+2z-1)((vz^2-z^2+3z-1)u+vz^2+z).
%\end{align*}
Making the substitutions~\eqref{eq:zuv} and using the known expressions of $F(x)$, $F(y)$, Equation~\eqref{eq:catalyticTamariLower} takes the form
$$
\tilde K_2 (z,u,v) \tilde{G}(u,v) = \tilde{L}_2(z,u,v,\tilde{G}_1(v),\tilde G_{11})
$$
for some rational functions $\tilde{K}_2, \tilde{L}_2$ that can be written explicitly~\cite{maple}.
%(the numerator of $\tilde{L}_2$ has 228 terms,).
%We can write Equation~\eqref{eq:catalyticTamariLower} in kernel form by factoring out $\tilde{G}$. The kernel $\tilde{K}_2$ is an explicit rational function of $z,u,v$. 
One easily checks that the kernel $\tilde{K}_2$ has a unique root $U_0\equiv U_0(z,v)$ which is a power series in $z$.
%, given by the root of the second factor. Explicitly, it is given by with $U_0=\frac{z(1+vz)}{1-3z+z^2-vz^2}=z+O(z^2)$.
Substituting $u=U_0$ in the last equation, we cancel the left-hand side,
hence we also cancel the right-hand side.
We are thus left with the following polynomial equation:
%satisfied by 
%$z,v,w, \tilde{G}_1(v)$ and $\tilde{G}_{11}$,
\begin{align}\label{eq:B}
\tilde{L}_2(z,U_0(z,v),v,\tilde{G}_1(v),\tilde G_{11})=0.
\end{align}
%from which we can eliminate the variable $u$. We thus obtain a polynomial equation satisfied by 
%$z,v,w, \tilde{G}_1(v)$ and $\tilde{G}_{11}$. 
The numerator of this equation has 91 terms, but it has only degree one in $\tilde{G}_1(v)$.
At this stage, we have eliminated the unknown $\tilde{G}(u,v)$ and the  variable $u$.

\smallskip
{\noindent \bf Second step: eliminating the variable $y$ (or $v$).}
%\label{sec:eliminating2}
%At this stage, we have shown that the series
%$\tilde{G}_1(v)$ and $\tilde{G}_{11}$ satisfy the explicit polynomial equation~\eqref{eq:B}.
%of the form
%\begin{align}\label{eq:tutte1}
%	B(\tilde G_1(v),\tilde{G}_{11},v,z,w)=0,
%\end{align}
%with $B(g_1,g_{11},v,z,w)$ of degree $1$ in $g_1$. 
%But 
The last equation, \eqref{eq:B}, is nothing but an equation with one catalytic variable, which is now the variable $v$ (or $x$)!
Since it is linear in $\tilde{G}_1(v)$, we can just use the kernel method again: one first checks that there is a unique series $V_0(z)$ cancelling the kernel, thus giving two equations: the linear and constant coefficient in $\tilde{G}_1(v)$ in~\eqref{eq:B}, which both vanish when $v=V_0$.
%, and using polynomial elimination one 
%
%We write the numerator of~\eqref{eq:B} in the form
%$$
%a(\tilde{G}_{11},v,z,w)\tilde{G}_1(v)+
%b(\tilde{G}_{11},v,z,w)=0,
%$$
%%In principle we can
%%We will solve this equation by using the Bousquet-Mélou--Jehanne method~\cite{BMJ} -- with proper justification since our equation is not a priori covered by the general statements in this reference. 
%%
%where in fact $a(\tilde{G}_{11},v,z,w)=a(v,z,w)$ does not depend on $\tilde{G}_{11}$. A direct check shows that the equation $a(V_0,z,w)=0$ as a unique solution $V_0(z)$ which is a formal power series in $z$, with $V_0=zw+O(z^2)$.
%Substituting $v=V_0$ in the equation, we finally obtain the equation
%$$
%b(\tilde{G}_{11},V_0(z),z,w)=0,
%$$
%which shows that $\tilde{G}_{11}$ is algebraic.
%Explicitly, we can eliminate the variable $v$ from the equations $a(v,z,w)$ and $b(\tilde{G}_{11},v,z,w)$, and 
Eliminating $V_0$, we obtain a polynomial equation for the function $\tilde{G}_{11}$
%
%
%we can eliminate the variable $v$ between the two polynomial $a(\tilde{G}_{11},v,z,w)$ and $b(\tilde{G}_{11},v,z,w)$.
%
%%\begin{lemma}
%	There exists a power series $V_0\equiv V_0(z,w)$, with $V_0=\frac{1}{w}z +O(z)$ which is solution of the equation
%\begin{align}\label{eq:tutte2}
%	B'_1(\tilde{G}_1(V_0),\tilde{G}_{11},V_0,z,w),
%\end{align}
%where $B'_1$ is the derivative with respect to the first variable. 
%\end{lemma}
%\begin{proof}
%The proof relies on an explicit check~\cite{maple} from the coefficients of the equation. 
%	Writing $B = a \tilde{G}_1(v)^2 + b \tilde{G}_1(v) + c$ with $a,b,c$ explicit polynomials in $\tilde{G}_{11},z,v,w$, the equation considered is $B'_1 = 2a \tilde{G}_1(v) + b.$
%	Now, we can check~\cite{maple} that
%	$$a =w(w-1)zv-z^2(w-1)+\dots \ \ , \ \ b=w(w-1)v^2-(w-1)(\tilde{G}_{11}w+w+1)vz+2\tilde{G}_{11}(w-1)z^2+\dots$$ where the neglected terms are at least cubic in variables $z,v$. 
%	Moreover, one checks from the combinatorial equation that  $\tilde{G}_1(v)=1+v+O(z)$ and $\tilde{G}_{11}=1+O(z)$.
%All this shows that the equation $2a\tilde{G}_1(v)+b=0$ is of the form:
%	\begin{align}\label{eq:taylor}
%	w(w-1)v^2-(w-1)vz = E(v,z),
%	\end{align}
%	where $E(v,z)$ is a formal power series (well defined, but unknown, since its expression involves $\tilde{G}_1(v)$ and $\tilde{G}_{11}$) having no sub-cubic monomial. This equation has a formal power series solution $v=v_0$, whose coefficients can be computed inductively, starting with the germ $V_0(z) = \frac{1}{w} z +\dots$.
%%\end{proof}
%
%We can now apply the general strategy of~\cite{BMJ}. If we substitute $v=V_0$ in~\eqref{eq:tutte1}, we obtain the following system of \emph{three} equations:
%\begin{align}\label{eq:tutteSys}
%\left\{
%	\begin{array}{rl}
%	B   (\tilde G_1(V_0), \tilde{G}_{11},V_0,z,w)   &=0\\
%	B'_1(\tilde G_1(V_0), \tilde{G}_{11},V_0,z,w)&=0\\
%	B'_3(\tilde G_1(V_0), \tilde{G}_{11},V_0,z,w)&=0\\
%\end{array}\right.
%\end{align}
%Indeed, the first two equations are~\eqref{eq:tutte1},~\eqref{eq:tutte2}, and the third one is implied by the first two, see~\cite{BMJ}.
%One can now eliminate the two variables $\tilde{G}_1(V_0)$ and $V_0$ from these three equations.
%%All these computations can be performed with a computer algebra system without conceptual difficulty. 
%It turns out that this equation is not even so big, and computations made~\cite{maple} we finally obtain:
which is not even so big. We obtain:
\begin{theorem}\label{thm:exactDown}
	%The generating function $G(x,y)$ is an algebraic function. Moreover, t
	The function $G(1,1)$ after the change of variables $z\leftrightarrow t$ given by~\eqref{eq:z} satisfies the polynomial equation $C(G(1,1),z,w)=0$ with
	\begin{align}\label{eq:exactDown}
		C(h,z,w)&= 
wz(-1+z)^9h^3+(-1+z)^6(2w^2z^2-w^2z+2z^2+w-z)h^2+4wz^2-4wz+w\nonumber \\
			-(-1&+z)^3(w^2z^3-3w^2z^2-2wz^3+w^2z-2wz^2+z^3+5wz-3z^2-2w+z)h.
	\end{align}
\end{theorem}
%\begin{proof}
%	We have already given justification to arrive at the system~\eqref{eq:tutteSys}, let us now explain how to perform the polynomial elimination in it, since our computer algebra software (Maple) was not able to do it in one call.
%
%	We first eliminate $\tilde{G}_1(V_0)$ from the first two equations of~\eqref{eq:tutteSys}, which results in a polynomial equation which has three irreducible factors. Two of these factors do not involve the quantity $\tilde{G}_{11}$, they are thus polynomial equations in $V_0,z,w$ and one can check that the series $V_0$ is \emph{not} a solution of them (to see this, it suffices to check the first term of the Puiseux expansion of the $V_0$-root of each factor near $z=0$, which do not coincide with that of the actual series $V_0=\frac{1}{w}z+\dots$). We thus only keep the third factor, which give us a polynomial equation of the form $E_1(\tilde{G}_{11},V_0,z,w)=0$. As it turns out, this polynomial has degree $1$ in $\tilde{G}_{11}$, we can thus solve it explicitly and we obtain that 
%	$\tilde{G}_{11}=R(V_0,z,w)$, with
%$$
%	R(v,z,w)
%=
%	\frac{v^2wz^2+2vwz^2+vwz-2vz^2+2wz^2-vw+vz-wz-2z^2+z}{wz(vz^3-3vz^2+z^3+3vz-3z^2-v+3z-1)}.
%$$
%	Of course at this stage $V_0$ is still unknown, so we are not done.
%
%	At this stage, we have two independent ways to conclude. The first one would be to eliminate $\tilde{G}_1(V_0)$ between the third and second equation, and perform a similar analysis, to finally obtain a second equation relating $\tilde{G}_{11}$ and $V_0$. This approch gives rise to relatively heavy computations.
%
%	The following way to proceed turned out to be more efficient. This approach does \emph{not} use the third equation of~\eqref{eq:tutteSys} (since this equation is implied by the other two, this is not as surprising as it could be) and it is closer in spirit to the historical applications of the quadratic method.
%	When we substitute the last obtained expression of $\tilde{G}_{11}$ into the two first equations of~\eqref{eq:tutteSys}, we see that these two equations factor nicely (up to trivial factors of the form $w-1, (1+v),$, or $z$), into the two equations
%$$
%	S(V_0,z,w) T(\tilde{G}_1(V_0),V_0,z,w)^2=0  \ \ , \ \ 	S(V_0,z,w) T(\tilde{G}_1(V_0),V_0,z,w)=0,
%$$
%	where $S$ and $T$ are explicit polynomials, in particular
%\begin{align*}
%	S(v,z,w)
%	&=
%	v^4w^2z^4+v^4wz^4+2v^3w^2z^4+3v^3w^2z^3+2v^2w^2z^4-2v^3w^2z^2+4v^3wz^3-2v^3z^4+3v^2w^2z^3
%	\\&
%	-v^3wz^2+v^3z^3-2v^2z^4+4vw^2z^3-3v^2w^2z+5v^2wz^2-3v^2z^3-4vw^2z^2+v^2w^2-2v^2wz+2v^2z^2
%	\\&
%	+vw^2z-4vz^3+2vwz-vw+vz-2z^2+z.
%\end{align*}
%	Now one can check, by looking at the expansion near $z=0$, that $T(\tilde{G}_1(V_0),V_0,z,w)$ is nonzero (again this only requires to know the first few orders of $V_0$ and $\tilde{G}_1(V_0)$, which we get from the original equation).
%	It follows that our series $V_0$ satisfies the equation 
%	$$S(V_0,z,w) =0.
%	$$
%	This equation, together with the above expression of $\tilde{G}_{11}$, gives us a compact parametrization of $\tilde{G}_{11}$. 
%
%	Finally, one can eliminate the variable $V_0$ from the two equations giving this parametrization, to obtain the algebraic equation given in the statement. All computations are done in~\cite{maple}.
%\end{proof}

\section{Asymptotics of moments}
\label{sec:asymptotics}


The two parts of Theorem~\ref{thm:tamari} are direct applications of Theorem~\ref{thm:trick}, up to computer algebra calculations done in~\cite{maple}.
In both cases we have $\rho=\frac{27}{256}$, corresponding to $z=\frac{1}{4}$.

%\subsection{Proof of Theorem~\ref{thm:tamari}}

In the case of the upper path, we start from~\eqref{eq:Hsexplicit} in Theorem~\ref{thm:exactUp}, which tells us that the function $f(t,s)=H(1,s)$ is algebraic.
%, therefore $H(x,r+1)$ is algebraic. Therefore (see e.g.~\cite{Comtet} again), there exists a recurrence formula of the form~\eqref{eq:Dfinite2} to compute its derivatives
%$$
%f^{(k)}\equiv f^{(k)}(t) :=\left.\left(\left(\frac{d}{ds}\right)^k H(1,s)\right)\right|_{s=1}.
%$$
%Given the form of our parametrization, we will prefer to work under the variable $z$ given by~\eqref{eq:z}, and in fact it is even more convenient to introduce the variable $\delta=(1-4z)^{1/2}$.
With GFUN~\cite{gfun}, we directly obtain~\cite{maple} a recurrence formula for the derivatives at $s=1$
%that the main singularity of $H(1,1)$ is at $t=\rho:=\frac{27}{256}$, which corresponds to $z=\frac{1}{4}$, i.e. $\delta=0$. In what follows we will still write $f^{(k)}$ for the function $f^{(k)}(t)$ expressed in the variable $\delta$.
%Using the package gfun, we immediately find~\cite{maple} that the $f^{(k)}$ satisfy an equation 
of the form~\eqref{eq:Dfinite2} with $L=6$,
%namely 
%\begin{align*}
%	f^{(k)}(t)= \sum_{d=1}^{L} h_d(t,k)  f^{(k-d)}(t) 
%\ \ , \ \  
%k \geq L,
%\end{align*}
where the $h_d(t,k)$ for $d=1..6$ are Laurent polynomials in $\delta=\sqrt{1-4z}$ and rational functions of $k$.
It is automatic to check that the  hypotheses of Theorem~\ref{thm:trick} holds with $\beta=\frac{3}{4}$ and $\alpha=\frac{1}{2}$. One can explicitly check the values of the corresponding constants $a_d(k)$, which are nothing but the coefficients of $\delta^{-3d}$ in $h_d$, up to a scaling factor. They are given~\cite{maple} by
\begin{align}\label{eq:explicitAi}
	a_1(k),\dots,a_6(k) = 0, \frac{\sqrt{6}}{96}(3k-4)(3k-8), 0, 0, 0, 0,
\end{align}
so the main recurrence formula~\eqref{eq:ck} becomes 
$$
c_k=\frac{\sqrt{6}(3k-4)(3k-8)}{96}c_{k-2} \ \ , \ \  k> 6.
$$
The initial conditions require
%Now it only remains to check the initial conditions (ii). From the explicit Equation~\eqref{eq:catalyticTamariUpper}, the dominant singulariy of $f^{(0)}$, and indeed of each given $f^{(i)}$, can be computed automatically. 
%We find that $f^{(0)}=\frac{64}{(\delta+3)^3}$, with singular expansion $f^{(0)}\singeq c_0 (1-t/\rho)^\alpha$ with $\alpha=\frac{1}{2}$ and $c_0=-\frac{32\sqrt{6}}{27}$. 
%
%Given $\alpha, \beta$ we compute $\ell_0=L-0=6$, and to check the initial conditions we need
to estimate the main singularity of $f^{(1)}, \dots, f^{(6)}$, which is easily done automatically,
%This is done automatically and one gets $f^{(\ell)}\singeq c_\ell (1-t/\rho)^{\alpha-\beta\ell}$ for $\ell \geq 6$, with 
%given by
%$$
%c_0,\dots,c_6=
%-{\frac {32}{27}}\,\sqrt {6},
%{\frac {16}{27}}\,{3}^{3/4}\sqrt [4]{2},
%{\frac {8}{27}},
%{\frac {5}{54}}\,\sqrt [4]{3}{2}^{3/4},
%{\frac {8}{81}}\,\sqrt {2}\sqrt {3},
%{\frac {385}{2592}}\,{3}^{3/4}\sqrt [4]{2},
%{\frac {70}{81}}.
%$$
%
%
%We have now verified al hypotheses of Theorem~\ref{thm:trick}.
%
%with explicit initial conditions above, and i
and it is then a direct check that the solution is given by 
$$
c_k=
{ {16}/{(27\pi)}}\, {\Gamma  \left( \tfrac{k}{2}+\tfrac{1}{3} \right) \Gamma 
 \left( \tfrac{k}{2}-\tfrac{1}{3} \right) \sqrt{2} \cdot {4}^{-k}{6}^{3/4\,k}}%{\pi }}
.
$$
for all $k\geq 0$.
%Applying Theorem~\ref{thm:trick}, we immediately get:
%\begin{align*}
%		\mathbf{E}\left[ \left(\frac{Q_n(I)}{n^{3/4}}\right)^k \right]\longrightarrow 
%	{\frac {\Gamma  \left( \frac{k}{2}+\frac{1}{3} \right) \Gamma  \left( \frac{k}{2}-\frac{1}{3} \right) {4}^{-k}{6}^{3/4\,k}}{\sqrt {3\pi }\Gamma  \left( \frac{3k}{4}-\frac{1}{2}\right) }}
%=\frac{\sqrt{3} \cdot2^{-\frac{k}{4}-1}}{\sqrt{\pi}} \frac{\Gamma(\frac{1}{4}k+\frac{1}{3})\Gamma(\frac{1}{4}k+\frac{2}{3})}{\Gamma(\frac{1}{4}k+\frac{1}{2})},
%	\end{align*}
%where the equality follows from  the duplication formula for the Gamma function (see also~\cite{maple}).
Applying Theorem~\ref{thm:trick}, we directly obtain \eqref{eq:mainMoment}, and~\eqref{eq:limitlaw} follows for example from the Carleman criterion. See~\cite{maple} for full calculations.

The proof of the second half (i.e. \eqref{eq:limitlawlow}) follows similarly from Theorem~\ref{thm:exactDown}, this time with $\alpha=\frac{1}{2}, \beta=\frac{3}{4}$, and $L=9$. See~\cite{maple} again.

%\subsection{Proof of Theorem~\ref{thm:tamariLow}}
%
%The proof is similar, we will be more succinct (all calculations are in~\cite{maple}).
%We start with Equation~\eqref{eq:exactDown} in Theorem~\ref{thm:exactDown}. 
%We use this equation to compute a recurrence formula for the functions
%$$
%f^{(k)}\equiv f^{(k)}(t) :=\left.\left(\left(\frac{d}{dw}\right)^k G(1,1)\right)\right|_{w=1}.
%$$
%%Given the form of our parametrization, we will prefer to work under the variable $z$ given by~\eqref{eq:z}, and in fact it is even more convenient to introduce the variable $\delta=(1-4z)^{1/2}$.
%%It is easily cheched~\cite{maple} that 
%The main singularity of $f^{(0)}$ (which is the same as in the previous section) is at $t=\rho:=\frac{27}{256}$, which corresponds to $z=\frac{1}{4}$, i.e. $\delta=0$.
%
%Using the package gfun, we immediately find~\cite{maple} that the $f^{(k)}$ satisfy an equation of the form~\eqref{eq:Dfinite2} with $L=9$,
%where again the $h_d(t,k)$ for $d=1..p$ rational function in $\delta$ with coefficients which are rational functions of $k$.
%The poles of $h_1,\dots,h_9$ at $\delta=0$ are respectively of order $0, 6, 6, 10, 10, 10, 10, 10, 10$.
%This implies that hypothesis (i) of Theorem~\ref{thm:trick} holds with $\beta=\frac{3}{4}$.
%One can explicitly check the values of the corresponding constants $a_d(k)$, which are nothing but top coefficients of $h_d$ in $\delta^{-1}$ up to a scaling factor. They are given~\cite{maple} by
%\begin{align}\label{eq:explicitAi2}
%	a_1(k),\dots,a_9(k) = 0, \frac{\sqrt{6}(3k-4)(3k-8)}{864}, 0, 0, 0, 0, 0, 0, 0.
%\end{align}
%
%Now it only remains to check the initial conditions (ii). From the explicit Equation~\eqref{eq:catalyticTamariUpper}, the dominant singulariy of $f^{(0)}$, and indeed of each given $f^{(i)}$, can be computed automatically. 
%The function  $f^{(0)}=\frac{64}{(\delta+3)^3}$ is the same as in the previous section, with singular expansion $f^{(0)}\singeq c_0 (1-t/\rho)^\alpha$ with $\alpha=\frac{1}{2}$ and $c_0=-\frac{32\sqrt{6}}{27}$. 
%
%We now have $\ell_0=L+0=9$, and to check the initial conditions we need to estimate the main singularity of $f^{(1)}, \dots, f^{(9)}$. This is done automatically and one gets $f^{(\ell)}\singeq c_\ell (1-t/\rho)^{\alpha-\beta\ell}$ for $\ell \leq 9$, with 
%$c_\ell$ given by
%\begin{center}${\footnotesize
%\displaystyle c_1,...,c_9=
%\frac{16\cdot 3^\frac{3}{4}2^\frac{1}{4}}{81},
%\frac{8}{243},
%\frac{5\cdot 2^\frac{3}{4}3^\frac{1}{4}}{1458},
%\frac{8\sqrt{6}}{6561},
%\frac{385\cdot 3^\frac{3}{4}2^\frac{1}{4}}{629856},
%\frac{70}{59049},
%\frac{85085\cdot 3^\frac{1}{4}2^\frac{3}{4}}{181398528},
%\frac{700\sqrt{6}}{1594323},
%\frac{37182145\cdot3^\frac{3}{4}2^\frac{1}{4}}{78364164096}.}
%$
%\end{center}
%
%
%We have now verified al hypotheses of Theorem~\ref{thm:trick}. The main recurrence formula becomes 
%$$
%c_k=\frac{\sqrt{6}(3k-4)(3k-8)}{864}c_{k-2} \ \ , \ \  k> 9,
%$$
%with initial conditions above. 
%It direct check to obtain the solution of this recurrence and, this concludes the proof of Theorem~\ref{thm:tamariLow} in the same way as we did for Theorem~\ref{thm:tamari}.
%%\begin{align*}
%%	\mathbf{E}\left[ \left(\frac{Q_n(I)-2P_n(I)}{n^{1/2}}\right)^k \right]\longrightarrow 
%%2\,{\frac { \left( \Gamma  \left( 1/2+1/2\,k \right)  \right) ^{2}{2}^
%%{-1/2\,k}}{ \left( k-1 \right) \sqrt {\pi }\Gamma  \left( 1/2\,k-1/2
%%	\right) }} \mathbf{1}_{k\text{ even}}.
%%	\end{align*}
%%which is equivalent to Theorem~\ref{thm:tamariLow} after using the duplication formula for the Gamma function, see~\cite{maple}.
%%



%\section{Other models (DELETE?)}
%\subsection{Tutte exploration tree}
%
%\subsection{Parking trees}


%\subsection*{Acknowledgements}
%
% The author acknowledges funding from the grants  ANR-19-CE48-0011 ``Combiné'' and ANR-18-CE40-0033 ``Dimers''.
%
%\bibliographystyle{plain}
%\bibliography{biblio}

 \printbibliography

%\appendix



\end{document}






In order to prove the convergence of moments, we need to study, for any fixed $p\geq 0$, the asymptotics of the coefficients of the derivative:
$$
H_p \equiv H_p(t):= \left.\frac{d^p}{d^ps}  H(1)\right|_{s=1}.
$$
In order to do this, we will study the behaviour of each of these series, now viewed as analytic functions of $t$, close to their dominant singularity.
In fact, we will prefer to work again in the variable $z$. Moreover, in order to study the successive derivatives $H_p$ at $s=1$, we will make use of the very general fact that algebraic functions are D-finite, meaning that their coefficients as power series are solution of linear recurrence equations. Applying this general fact (to derivatives at the point $s=1$ with $z$ a spectator variable) will give us a very convenient, and essentially automatic, way to study the successive derivatives $H_p$ and their singular behaviour\footnote{we have not sure to have already seen this technique used in this way. That might be the main technical innovation of the present paper.}

The generating function $F(1)$ has a unique singularity for $t=27/256$, which corresponds to the point $z=1/4$. In the sake of uniformization, we will introduce the variable 
$\delta$ related to $z$ by 
\begin{align}\label{eq:delta}
	z=\frac{1}{4}-\delta^2.
\end{align}

The equation~\eqref{eq:Hsexplicit} implies that $H(1)$ is an algebraic function of $z$ and $s$. Its minimal polynomial can be explictly computed easily~\cite{maple}. This gives us an explicit polynomial such that
\begin{align}\label{eq:algH}
B(H(1),z,s)=1.
\end{align}
This polynomial has degree 3 in $H(1)$ and (only) 42 terms.


I follows by induction, differentiating~\eqref{eq:algH}, that for any $p\geq 0$, the derivative
$$
\left(\frac{d}{ds}\right)^p H(1)
$$
can be written as a rational function in $z,s,H(1)$ with denominator $B'_1(H(1),z,s)^p$.


%\begin{align*}\scriptsize
%h^3s^3z^{10}-9h^3s^3z^9+36h^3s^3z^8-84h^3s^3z^7+126h^3s^3z^6+4h^2s^2z^8-126h^3s^3z^5-26h^2s^2z^7+84h^3s^3z^4+73h^2s^2z^6-36h^3s^3z^3-116h^2s^2z^5-hs^2z^6+9h^3s^3z^2+115h^2s^2z^4+6hs^2z^5+2hsz^6-h^3s^3z-74h^2s^2z^3-15hs^2z^4-4hsz^5-hz^6+31h^2s^2z^2+20hs^2z^3-hsz^4+6hz^5-8h^2s^2z-15hs^2z^2+7hsz^3-15hz^4+h^2s^2+6hs^2z-5hsz^2+20hz^3-hs^2+hsz-15hz^2+6hz+4z^2-h-4z+1
%\end{align*}
