%% if you are submitting an initial manuscript then you should have submission as an option here
%% if you are submitting a revised manuscript then you should have revision as an option here
%% otherwise options taken by the article class will be accepted
\documentclass[finalversion]{FPSAC2025}
\articlenumber{116}
%% but DO NOT pass any options (or change anything else anywhere) which alters page size / layout / font size etc

%% note that the class file already loads {amsmath, amsthm, amssymb}

\definecolor{cb-yellow}{RGB}{221,170,51}
\definecolor{cb-red} {RGB}{187,85,102}
\definecolor{cb-teal}{RGB}{0,153,136}
\definecolor{cb-blue} {RGB}{0,68,136}
\definecolor{cb-green}{RGB}{17,119,51}
\definecolor{cb-purple} {RGB}{170,68,153}
\definecolor{cb-palegrey} {RGB}{221,221,221}

\newcommand{\colS}{red}
\newcommand{\colT}{yellow}

\usepackage{tikz,pgfplots}
\usetikzlibrary{positioning}
\pgfplotsset{compat=1.18}

% \usepackage{hyperref}
% \hypersetup{
%   pdfauthor={Michael Borinsky, Chiara Meroni, Maximilian Wiesmann},
%   pdftitle={Asymptotic count of edge-bicolored graphs},
%   pdfsubject={},
%   pdfkeywords={}
% }

\usepackage{algorithm}
\usepackage[noend]{algpseudocode}

\newtheorem{theorem}{Theorem}[section]
%\newtheorem*{theorem*}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{problem}[theorem]{Problem}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{examplex}[theorem]{Example}
\newenvironment{example}
{\pushQED{\qed}\renewcommand{\qedsymbol}{$\diamond$}\examplex}
{\popQED\endexamplex}

\newcommand{\GG}{\mathcal{G}}
\newcommand{\GGA}{\mathcal{G}^\star}
\DeclareMathOperator{\SG}{{\mathbb S}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\dd}{\mathrm{d}}
\DeclareMathOperator{\crit}{crit}
\DeclareMathOperator{\Aut}{Aut}
\newcommand{\ol}{\Lambda}
\newcommand{\bigO}{\mathcal{O}}



%% define your title in the usual way
\title{Asymptotic count of edge-bicolored graphs}  %[Optional shorter title]

%% define your authors in the usual way
%% use \addressmark{1}, \addressmark{2} etc for the institutions, and use \thanks{} for contact details
%% for more than 3 names, use the Oxford comma: One Author, Two Author, Red Author, \and Blue Author
\author[M. Borinsky, C. Meroni, M. Wiesmann]{\texorpdfstring{Michael Borinsky,\thanks{\href{mailto:michael.borinsky@pitp.ca}{michael.borinsky@pitp.ca}. 
Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Colleges and Universities. Michael Borinsky was supported by Dr.\ Max Rössler, the Walter Haefner Foundation, and the ETH Zürich Foundation.}\addressmark{1,2} Chiara Meroni,\thanks{\href{mailto:chiara.meroni@eth-its.ethz.ch}{chiara.meroni@eth-its.ethz.ch}. Chiara Meroni was supported by Dr.\ Max Rössler, the Walter Haefner Foundation, and the ETH Zürich Foundation.}\addressmark{2} and Maximilian Wiesmann\thanks{\href{mailto:wiesmann@mis.mpg.de}{wiesmann@mis.mpg.de}.}\addressmark{3}}{Michael Borinsky, Chiara Meroni, and Maximilian Wiesmann}}

%% then use \addressmark to match authors to institutions here
\address{
\addressmark{1}Perimeter Institute, Waterloo, Canada.\\
\addressmark{2}ETH Institute for Theoretical Studies, Z\"urich, Switzerland. \\
\addressmark{3}Max Planck Institute for Mathematics in the Sciences, Leipzig, Germany.}

%% put the date of submission here
\received{\today}

%% leave this blank until submitting a revised version
%\revised{}

%% put your English abstract here, or comment this out if you don't have one yet
%% please don't use custom commands in your abstract / resume, as these will be displayed online
%% likewise for citations -- please don't use \cite, and instead write out your citation as something like (author year)
\abstract{
We show that specific exponential bivariate integrals serve as generating functions of labeled edge-bicolored graphs.
Based on this, we prove an asymptotic formula for the number of regular edge-bicolored graphs with arbitrary weights assigned to different vertex structures. The asymptotic behavior is governed by the critical points of a polynomial.
As an application, we discuss the Ising model on a random 4-regular graph and show how its phase transitions arise from our formula.}

%% put your French abstract here, or comment this out if you don't have one
% \resume{\lipsum[2]}

%% put your keywords here, or comment this out if you don't have them yet
\keywords{Edge-bicolored graphs, generating function, asymptotic behavior, bivariate exponential, critical points}

%% you can include your bibliography however you want, but using an external .bib file is STRONGLY RECOMMENDED and will make the editor's life much easier
%% regardless of how you do it, please use numerical citations; i.e., [xx, yy] in the text

%% this sample uses biblatex, which (among other things) takes care of URLs in a more flexible way than bibtex
%% but you can use bibtex if you want
\usepackage[backend=bibtex, giveninits]{biblatex}
\addbibresource{references.bib}
%% note the \printbibliography command at the end of the file which goes with these biblatex commands

\begin{document}

\maketitle
%% note that you DO NOT have to put your abstract here -- it is generated by \maketitle and the \abstract and \resume commands above

\section{Introduction}\label{sec:intro}

The (asymptotic) enumeration of graphs is a subject with a long history (see, e.g., \cite{MR3966531} and the references therein) and many variants, such as assumptions regarding the sparsity of the graph, or the presence of weights that decorate vertices and edges as well as the whole graph itself. 
In this work, which is an extended abstract of the paper \cite{BMW:EdgeBicolored} by the same authors, we focus on central objects in extremal graph theory, namely edge-bicolored graphs (see \cite[Chapter~V]{bollobas2004extremal}, \cite{CaiYe:EdgeBicolored}).
Each graph is weighted by the reciprocal of the order of its automorphism group and the product of an arbitrary set of parameters $\Lambda_{u,w}$ assigned to each bicolored incidence structure (encoded in the pair $(u,w)$) of a vertex. For instance, the graph in \Cref{fig:graph_intro} has weight $\frac{1}{8}\Lambda_{2,2}\Lambda_{0,4}$, since the left vertex has two red and two yellow half-edges and the right vertex has four yellow half-edges; the order of the automorphism group of this graph is $8$, which is explained in \Cref{sec:graphs}.
\begin{figure}[H]
    \centering
    \scalebox{0.8}{
    \input{FigColored/threebubble_ryy}}
    \caption{An edge-bicolored graph.}
    \label{fig:graph_intro}
\end{figure}

We will demonstrate that the number of these weighted graphs, with fixed Euler characteristic (difference between number of vertices and edges), gives the coefficients of the large-$z$ asymptotic expansion of the following family of bivariate integrals:
\begin{equation}
\label{eq:exponentialIntegral}
I(z) = \frac{z}{2\pi}
\int_D \exp\left( {z\, g(x,y)} \right) \, \dd x \dd y\,,
\end{equation}
where $D$ is a neighborhood of the origin, and $g$ is a polynomial fulfilling specific conditions that we will describe below.
Integrals as $I(z)$ arise naturally in applications. For instance, they appear in Bayesian statistics as \emph{marginal likelihood integrals} (e.g.,\ \cite[Section 1]{watanabe2009algebraic}), but they can also be interpreted as \emph{path integrals} associated to a zero-dimensional quantum system with two interacting fields (e.g.,~\cite[Section 2]{skinner2018quantum} or \cite{bessis1980quantum}). The setups might differ, however, in the integration domain $D$ and the assumption on $g$, leading to different asymptotic behaviors (see \cite{lin2011algebraic} for an asymptotic analysis in the realm of statistics).

After a formal definition of edge-bicolored graphs in \Cref{sec:graphs}, we draw a connection between them and the integral in \eqref{eq:exponentialIntegral}, using a bivariate version of the Laplace method (cf. \Cref{prop:laplace}). Then, we interpret the coefficients of the asymptotic expansion from a combinatorial point of view. 
Therefore, $I(z)$ plays the role of a \emph{generating function} of edge-bicolored graphs (cf. \Cref{thm:laplace-graphs}).
We derive an effective algorithm (cf. \Cref{alg:coeffs}) for the computation of these coefficients. The code is available at \cite{mathrepo}, implemented in \texttt{Julia}.
In \Cref{sec:asymptotics}, we prove an asymptotic formula for the weighted number of \emph{regular} edge-bicolored graphs of fixed degree, in the limit where the number of edges and vertices goes to infinity. Our main result, \Cref{thm:main} below, relates this asymptotic formula to the critical points of a polynomial $g(x,y)$ whose shape is governed by the vertex incidence structure of the graphs. 
\begin{theorem}
    \label{thm:main}
    Let $k\geq 3$ and consider a polynomial 
    $$g(x,y) = - \frac{x^2}{2} - \frac{y^2}{2} + \sum_{\substack{u+w = k}} \ol_{u,w} \frac{x^u y^w}{u!w!}.$$
    Let $\Psi$ be the set of its nonzero complex critical points with real coordinate ratio and smallest distance to the origin, and assume these points are non-degenerate (full-rank Hessian $H_g$). Then,
    \[
    \sum_{G \in \GG_{-n}^k} \frac{1}{ \left\lvert \Aut(G) \right\rvert} \prod \limits _{v\in V^G} \ol_{\deg(v)} \sim 
    \frac{1}{2\pi}
    \Gamma(n)
    \sum\limits_{(x,y) \in \Psi} \frac{(-g(x,y))^{-n}}{\sqrt{-\det H_g(x,y)}}, \qquad n\to \infty,
    \]
    where $\GG_{-n}^k$ is the set of all regular edge-bicolored graphs with vertex degree $k$ and Euler characteristic $-n$, and $\deg(v)\in \Z^2$ is the number of \colS{}/\colT{} half-edges in a vertex $v$ of $G\in \GG_{-n}^k$.
\end{theorem}
We showcase that, unlike the monochromatic case, which has previously been discussed in \cite[Chapter~3]{Borinsky:2018mdl}, only critical points satisfying some reality constraint (namely, in $\C\cdot \R^2$) contribute to the asymptotics.
We conclude with \Cref{sec:outlook}, with two conjectures that generalize our main result and are the subject of ongoing and future work.

\section{Edge-bicolored graphs}
\label{sec:graphs}

A \emph{graph} is an at most one-dimensional finite CW complex, i.e., it has finitely many vertices and undirected edges, with loops and multiple edges allowed.\footnote{In parts of combinatorics and computer science, graphs with loops and multiple edges are often defined as a pair of a set of vertices and a \emph{multiset} of edges. For our purposes, this necessitates the cumbersome definition of \emph{compensation factors} to write down appropriate generating functions (see, e.g. \cite[Equation 1.1]{janson1993birth}). To avoid this, we borrow the elementary notion of a CW complex from topology, which is also natural in the (quantum) physical context.} It is \emph{edge-bicolored} if each edge has one of two different colors. We will represent graphs using only discrete data.
A \emph{(set) partition} $P$ of a finite set $H$ is a set of non-empty and mutually disjoint subsets of $H$
whose union equals $H$. The elements of $P$ are called \emph{blocks}.
\begin{definition}
\label{def:bigraph}
Given two disjoint finite sets $S$ and $T$ of labels, an \emph{$[S,T]$-labeled edge-} {\parfillskip=0pt \par} \noindent \emph{bicolored graph} is a tuple $\Gamma = (V,E_S,E_T)$, where
the \emph{vertex set} $V$ is a partition of $S\sqcup T$, $E_S$ is a partition of $S$ into blocks of size 2, and $E_T$ is a partition of $T$ into blocks of size 2.
% \begin{enumerate}
% \item the \emph{vertex set} $V$ is a partition of $S\sqcup T$,
% \item $E_S$ is a partition of $S$ into blocks of size 2,
% \item $E_T$ is a partition of $T$ into blocks of size 2.
% \end{enumerate}
\end{definition}
\noindent
We think of the elements of $S$ and $T$ as \emph{half-edge labels} colored 
\colS{} and \colT{}, respectively.
These half-edges are bundled together in vertices via the partition $V$. 
The edge sets $E_S$ and $E_T$ pair the half-edges into edges of the respective color. Every edge-bicolored graph without isolated vertices can be represented by at least one $[S,T]$-labeled graph.
All graphs here are edge-bicolored, so we drop this adjective now. Half-edge labeled graphs are a common graph model in the study of Feynman graphs, see \cite[Section 5]{yeats2017combinatorial}.

\begin{example}
    Let $S = \{s_1,s_2,\dots,s_6\}$ and $T = \{t_1,t_2\}$. The partitions 
    \begin{gather*}
        V = \left\{ \{s_1,s_2,s_3,s_4\}, \{s_5,s_6,t_1,t_2\} \right\}, \\
        E_S = \left\{ \{s_1,s_2\}, \{s_3,s_4\}, \{s_5,s_6\} \right\}, \quad
        E_T = \left\{ \{t_1,t_2\} \right\},
    \end{gather*}
    form an $[S,T]$-labeled graph representing the graph depicted in \Cref{fig:extwotadpole2s_rrry}.
\end{example}
\begin{figure}[ht]
    \centering
    \scalebox{0.5}{\input{FigColored/twotadpole2s_rrry_horizontal}}
    \caption{An edge-bicolored graph with two connected components.}
\label{fig:extwotadpole2s_rrry}
\end{figure}

It follows immediately from our definition that there are many equivalent ways to assign labels to the same graph. In order to get rid of this redundancy, we define isomorphisms. An \emph{isomorphism} from an $[S_1,T_1]$-labeled graph $(V^{1},E_S^{1},E_T^{1})$ to an $[S_2,T_2]$-labeled graph $(V^{2},E_S^{2},E_T^{2})$
is a pair of bijections $j_S:S_1\rightarrow S_2$, $j_T:T_1\rightarrow T_2$ such that $j(V^{1}) = V^{2}, j(E_S^{1}) = E_S^{2}$, and $j(E_T^{1}) = E_T^{2}$
with $j$ being the map that $j_S$ and $j_T$ canonically induce on the subsets of $S$, $T$, and $S\sqcup T$.
An \emph{automorphism} of an $[S,T]$-labeled graph $\Gamma$ is an isomorphism to itself. Those form the group $\Aut(\Gamma)$. By the orbit stabilizer theorem, the isomorphism class of any $[\{1,\ldots,2s\},\{1,\ldots,2t\}]$-labeled graph $\Gamma$ has size
\begin{equation}\label{eq:size_iso_class}
\frac{(2s)!(2t)!}{|\Aut(\Gamma)|}.
\end{equation}

\begin{example}
    An $[S,T]$-labeled graph $\Gamma$ representing the graph depicted in \Cref{fig:extwotadpole2s_rrry} has automorphism group isomorphic to $\left(
    \SG_2 \times \SG_2 \rtimes \SG_2 \times \SG_2 \right) \times \SG_2$, where $\rtimes$ denotes the semidirect product of groups. It is a subgroup of $\SG_6 \times \SG_2$, where $\SG_6$ refers to permutations of the six {\colS} half-edges in $S$ and $\SG_2$ to the two {\colT} half-edges in $T$.
\end{example}

We write $\GG$ for the set of isomorphism classes of graphs. For each $G \in \GG$, we write $V^G$, $E_S^G$, $E_T^G$, $E^G=E_S^G\sqcup E_T^G$ and $\Aut(G)$ 
for the respective set/group of some $[S,T]$-labeled representative of $G$.
The \emph{Euler characteristic} of $G$ is defined by $\chi(G) = |V^G|-|E^G|$.
The \emph{bidegree} of a graph's vertex $v\in V^G$ is the pair of integers $\deg(v) = (u,w)$ where $u$ counts the number of half-edges in $v$ that lie in the \colS{}-colored set $S$ and $w$ the half-edges in the \colT{}-colored part $T$. The \emph{vertex degree} of $v$ is $|\deg(v)| = u+w$.

We define a family of polynomials $a_{s,t}$ indexed by integers $s,t \geq 0$ in a two-fold infinite set of variables $\lambda_{u,w}$ indexed by $u,w \geq 0$ with $u+w \geq 1$.  Consider the ring of polynomials in these variables $\mathcal R = \Q[\lambda_{0,1},\lambda_{1,0},\lambda_{1,1}, \lambda_{0,2}, \ldots]$.
The polynomials $a_{s,t}(\lambda) \in \mathcal R$ are defined by the generating function
\begin{equation}
\label{eq:coeffs_extraction_gen_funct}
    \sum_{s,t \geq 0} a_{s,t}(\lambda)\, x^s y^t =
\exp\left(\sum_{\substack{u,w \geq 0\\u+w \geq 1}} \lambda_{u,w} \frac{x^u y^w}{u!w!}\right) \in \mathcal R[[x,y]].
\end{equation}
For instance,
$a_{0,0}(\lambda) = 1$,
$a_{1,0}(\lambda) = \lambda_{1,0}$,
and 
$a_{2,0}(\lambda) = \frac12 ( \lambda_{2,0} + \lambda_{1,0}^2 )$.

\begin{proposition}
\label{prop:genfun2colgraphs}
The generating function for graphs with marked bidegrees is
\begin{gather*}
        \sum_{G\in \GG} \frac{\eta^{|E^G|} }{
\left\lvert \Aut(G) \right\rvert} 
\prod \limits _{v\in V^G}
\lambda_{\deg(v)}
= \sum_{s,t\geq 0} \eta^{s+t} \cdot  (2s-1)!!\cdot (2t-1)!!\cdot
a_{2s,2t}(\lambda) \in \mathcal R[[\eta]],
\end{gather*}
where $a_{s,t}$ is defined as in \eqref{eq:coeffs_extraction_gen_funct}.
\end{proposition}

This generating function can also be seen as the \emph{exponential generating function} of bicolored graphs with labeled half-edges in the language of \cite{flajolet2009analytic}. Our convention that involves a $1/|\Aut(G)|$ factor and the sum over all isomorphism classes of graphs is more in line with the category theoretical perspective from~\cite{joyal1981theorie} and comes with a lighter notation, as we do not need to retain the number of half-edges as a counting variable.

The proof of \Cref{prop:genfun2colgraphs} relies on the size of isomorphism classes in \eqref{eq:size_iso_class} and on the number of partitions of a set $S\sqcup T$ with prescribed numbers of elements in $S$ and $T$. We illustrate the result with an example.
\begin{example}
    \Cref{prop:genfun2colgraphs} provides a recipe to count our graphs for a given number of edges, grouping them according to their bidegrees. For instance, the coefficient of $\eta^1$ on the left-hand side counts graphs with one edge:
    \begin{align*}
    \sum_{\substack{G\in \GG, \\ |E^G|=1}} \frac{1}{\left\lvert \Aut(G) \right\rvert} \prod \limits _{v\in V^G} \lambda_{\deg(v)} &= \;\;\; \scalebox{.3}{\input{FigColored/tadpole_r}} \;\;+\;\;\: \scalebox{.3}{\input{FigColored/segment_r}} \:\;\;+\;\;\: \scalebox{.3}{\input{FigColored/tadpole_y}} \:\;\;+\;\; \scalebox{.3}{\input{FigColored/segment_y}} \\
    &= \frac{1}{2} \lambda_{2,0} + \frac{1}{2} \lambda_{1,0}^2 + \frac{1}{2} \lambda_{0,2} + \frac{1}{2} \lambda_{0,1}^2.
    \end{align*}
    Using the power series on the right-hand side of \Cref{prop:genfun2colgraphs}, this can be obtained simply as $a_{2,0}+a_{0,2}$, and by expanding the exponential in \eqref{eq:coeffs_extraction_gen_funct}, we get exactly the above expression. If for $|E^G|=1$ these two approaches may seem equally complicated, already for graphs with two edges it is clear that the use of the generating function speeds up the computation. In fact, there are seven (monochromatic) graphs with two edges: \input{figs/segment2comp}\, , \input{figs/tadpole2comp}\, , \input{figs/segmenttadpole_2comp}\, , \input{figs/tadpole2}\, , \input{figs/banana2}\, , \input{figs/segmenttadpole_1comp}\, , \input{figs/segments2}, which turn into $23$ edge-bicolored graphs. On the other hand, a simple expansion of the exponential function gives
    \begin{align*}
        \sum_{\substack{G\in \GG, \\ |E^G|=2}} &\frac{1}{\left\lvert \Aut(G) \right\rvert} \prod \limits _{v\in V^G} \lambda_{\deg(v)} = 3 a_{4,0} + a_{2,2} + 3 a_{0,4} \\
        =&\,\, \lambda_{0,1} \lambda_{1,0} \lambda_{1,1}+\tfrac{\lambda_{0,1}^4}{8}+\tfrac{3 \lambda_{0,1}^2 \lambda_{0,2}}{4}+\tfrac{\lambda_{0,1}^2 \lambda_{1,0}^2}{4}+\tfrac{\lambda_{0,1}^2 \lambda_{2,0}}{4}
        +\tfrac{\lambda_{0,1} \lambda_{0,3}}{2}
        +\tfrac{\lambda_{0,1} \lambda_{2,1}}{2}
        +\tfrac{3 \lambda_{0,2}^2}{8}+\tfrac{\lambda_{0,2} \lambda_{1,0}^2}{4}\\
        &+\tfrac{\lambda_{0,2} \lambda_{2,0}}{4}+\tfrac{\lambda_{0,4}}{8}+\tfrac{\lambda_{1,0}^4}{8}
        +\tfrac{3 \lambda_{1,0}^2 \lambda_{2,0}}{4}+\tfrac{\lambda_{1,0} \lambda_{1,2}}{2}+\tfrac{\lambda_{1,0} \lambda_{3,0}}{2}+\tfrac{\lambda_{1,1}^2}{2}+\tfrac{3 \lambda_{2,0}^2}{8}+\tfrac{\lambda_{2,2}}{4}+\tfrac{\lambda_{4,0}}{8}.\qedhere
    \end{align*}
\end{example}

\subsection{Bivariate Laplace method}\label{sec:laplace}
Next we want to interpret the generating function for edge-bicolored graphs as a large $z$ expansion of the bivariate integral \eqref{eq:exponentialIntegral}.
As anticipated in the introduction, we assume that the integral $I(z)$ exists for $z > 0$, the domain $D$ contains a neighborhood of the origin, and that the polynomial $g$ attains its \emph{unique} supremum in $D$ at the origin, around which it has a converging expansion of the form
\begin{equation}
\label{eq:gdef}
g(x,y) = - \frac{x^2}{2} - \frac{y^2}{2} + \sum_{\substack{u,w \geq 0\\ u+w \geq 3}} \ol_{u,w} \frac{x^u y^w}{u!w!}\, .
\end{equation}
The last condition ensures that $I(z)$ resembles a Gaussian integral when $x$ and $y$ are small. This allows to approximate the integral by a slightly perturbed Gaussian for large $z$ and relate its \emph{asymptotic expansion} to the polynomials $a_{s,t}$ defined in \eqref{eq:coeffs_extraction_gen_funct}.

Given functions $f(z), g(z), h(z)$, the notation $f(z) = g(z) + \bigO(h(z))$ means that the limit superior of $\left|\frac{f(z)-g(z)}{h(z)}\right|$ for large $z$ is finite.
The asymptotic expansion notation
$f(z) \sim \sum_{n \geq 0} g_n(z)$  means that  $f(z) - \sum_{n=0}^{R-1} g_n(z) \in \bigO(g_{R}(z))$
for all $R \geq 0$.

\begin{proposition}
\label{prop:laplace}
If $I(z)$, $g$, $D$ and $\ol_{u,w}$ are related as above, then $I(z) \sim \sum_{n\geq 0} A_n z^{-n}$, for large $z$, where $A_n$ is the coefficient of $z^{-n}$ in the formal power series
\[
    \sum_{s,t\geq 0} z^{-(s+t)} (2s-1)!! \cdot (2t-1)!! \cdot a_{2s,2t} (z \cdot \ol) \in \R[[z^{-1}]]\,,
\]
where $(2s-1)!! = (2s-1)(2s-3)\cdots 3\cdot 1$ and $a_{2s,2t} (z \cdot \ol)$ is the polynomial $a_{2s,2t}(\lambda)$ defined in \eqref{eq:coeffs_extraction_gen_funct},
with 
\begin{equation}
    \label{eq:substitution_lambda_Lambda}
    \lambda_{u,w} = \left\{ \begin{array}{ll}
        0 & \text{for}~ u,w\geq 0 ~\text{and}~ 1 \leq u+w < 3, \\
        z\,\Lambda_{u,w} & \text{for}~ u,w\geq 0 ~\text{and}~ u+w \geq 3 .\\
    \end{array} \right.
\end{equation}
\end{proposition}

\noindent
The proof of \Cref{prop:laplace} uses the classical \emph{Laplace method} for the asymptotic expansion of the integral $I(z)$. The result allows to compute the coefficients $A_n$ explicitly using \Cref{alg:coeffs}. This recursive algorithm is implemented in \texttt{Julia} and is available at \cite{mathrepo}.
\begin{algorithm}[h!]
    \caption{Computing $A_0,\ldots,A_n$}
    \label{alg:coeffs}
    \textsc{Input:} $n\in \Z$, $\Lambda_{u,w}$ for all $u,w \geq 0$, $u+w\leq 2n + 2$\\
    \textsc{Output:} $A_0,A_1,\ldots,A_n$ such that $I(z) \sim \sum A_m z^{-m}$.
    \begin{algorithmic}[1]
    \State $F_k (x,y) \gets \sum_{u+w=k+2} \Lambda_{u,w}\frac{x^u y^w}{u!w!}$ for $k \in \{1,\ldots, 2n\}$
    \State $Q_0(x,y) \gets 1$
    \For{$i = 1,\dots, 2n$}
        \State $Q_i(x,y) \gets \frac{1}{i} \sum_{k=1}^{i} F_{k}(x,y) Q_{i-k}(x,y)$
        \State $q^{(i)}_{s,t} \gets$ coefficients of $Q_{i}(x,y) = \sum_{s,t\geq 0}q^{(i)}_{s,t} x^s y^t$
    \EndFor
    \State \Return $A_i \gets \sum_{s,t \geq 0} (2s-1)!! (2t-1)!! \, q^{(2i)}_{2s,2t}$  for $i \in \{1,\ldots, n\}$
    \end{algorithmic}
\end{algorithm}


\begin{example}
\label{ex:Ising_model_deg4_part1}
Fix $D = [-1,1]^2$ and $g(x,y) = -\frac{x^2}{2} - \frac{y^2}{2} + \frac{x^4}{4!} + \lambda \frac{x^2 y^2}{4} + \lambda^2 \frac{y^4}{4!}$ with $\lambda$ some nonnegative real number. 
By \Cref{prop:laplace}, the associated integral $I(z)$ has an asymptotic expansion
$I(z) \sim \sum_{n \geq 0} A_n z^{-n}$.
Combining \Cref{prop:laplace} and \eqref{eq:coeffs_extraction_gen_funct}, we find that
    \begin{align*}
        A_0 &= 1, \\
        A_1 &= \frac{1}{8} + \frac{1}{4}\lambda + \frac{1}{8}\lambda^2, \\
        A_2 &= \frac{35}{384} + \frac{5}{32}\lambda + \frac{19}{64}\lambda^2 + \frac{5}{32}\lambda^3 + \frac{35}{384}\lambda^4,\\
        A_3 &= \frac{385}{3072} + \frac{105}{512}\lambda + \frac{1295}{3072}\lambda^2 + \frac{175}{256}\lambda^3 + \frac{1295}{3072}\lambda^4 + \frac{105}{512}\lambda^5 + \frac{385}{3072}\lambda^6.
    \qedhere
    \end{align*}
\end{example}

We are ready to endow the obtained analytic expressions for $I(z)$ with a combinatorial interpretation, putting together \Cref{prop:genfun2colgraphs,prop:laplace}. This process is inspired by quantum field theory in physics, where perturbative expansions of observables, which are combinatorially controlled via \emph{Feynman graphs} (e.g., edge-bicolored graphs), relate to \emph{path integrals} (e.g., $I(z)$). 

\begin{theorem} 
\label{thm:laplace-graphs}
If $I(z)$, $g$, $D$ and $\ol_{u,w}$ are related as above, then $I(z) \sim \sum_{n\geq 0} A_n z^{-n}$,
for large $z$, where $A_n$ is given by
\[
A_n = \sum_{G \in \GGA_{-n}} 
        \frac{1 }{
\left\lvert \Aut(G) \right\rvert} 
\prod \limits _{v\in V^G}
\ol_{\deg(v)},
\]
where we sum over the set $\GGA_{-n}$ of all edge-bicolored graphs with vertex degrees at least $3$ and Euler characteristic equal to $-n$.
\end{theorem}

Hidden in the statement is the connection between the polynomial $g$ and the graphs we are actually considering in the summation for $A_n$. Indeed, because of the product of $\Lambda_{\deg(v)}$, if the monomial $x^s y^t$ does not appear in $g$, then the graphs that will contribute to $A_n$ will not have any vertex of bidegree $(s,t)$. The coefficients of the nonzero monomials can be interpreted as weights on the vertices with corresponding bidegree.

\begin{example}\label{ex:Ising_model_deg4_part2}
Continuing \Cref{ex:Ising_model_deg4_part1}, let $c^{(k)}_n$ be the coefficient of $\lambda^k$ in $A_n$. By \Cref{thm:laplace-graphs}, $c^{(k)}_n$ counts automorphism-weighted graphs with Euler characteristic $-n$ and vertex degree four, such that $k_1$ vertices have exactly two yellow half-edges and $k_2$ vertices have four yellow half-edges, so that $k_1+k_2 = k$. We can view this explicitly for $n=2$, as follows. Among the $21$ (monochromatic) graphs with $\chi = -2$, there are only three $4$-regular graphs. These are \input{figs/twotadpole2s}, \input{figs/banana4}, \input{figs/threebubble}. Considering all bicolorings, we get
    \begin{align*}
        c_2^{(0)} &= \scalebox{.3}{\input{FigColored/twotadpole2s_rrrr}} \; +\;  \scalebox{.3}{\input{FigColored/banana4_rrrr}} \;+\; \scalebox{.3}{\input{FigColored/threebubble_rrr}} = \tfrac{1}{128} + \tfrac{1}{48} + \tfrac{1}{16} = \tfrac{35}{384}, \\
        c_2^{(1)} &= \scalebox{.3}{\input{FigColored/twotadpole2s_rrry}} \;+\; \scalebox{.3}{\input{FigColored/threebubble_rry}} = \tfrac{1}{32} + \tfrac{1}{8} = \tfrac{5}{32}, \\
        c_2^{(2)} &= \scalebox{.3}{\input{FigColored/twotadpole2s_rryy}}  + \scalebox{.3}{\input{FigColored/twotadpole2s_ryry}} +  \scalebox{.3}{\input{FigColored/banana4_rryy}} + \scalebox{.3}{\input{FigColored/threebubble_ryr}} + \scalebox{.3}{\input{FigColored/threebubble_yry}} = \tfrac{1}{64} + \tfrac{1}{32} + \tfrac{1}{8} + \tfrac{1}{16} + \tfrac{1}{16} = \tfrac{19}{64}, \\   
        c_2^{(3)} &= \scalebox{.3}{\input{FigColored/twotadpole2s_ryyy}} \;+\; \scalebox{.3}{\input{FigColored/threebubble_ryy}} = \tfrac{1}{32} + \tfrac{1}{8} = \tfrac{5}{32}, \\
        c_2^{(4)} &= \scalebox{.3}{\input{FigColored/twotadpole2s_yyyy}} \; +\;  \scalebox{.3}{\input{FigColored/banana4_yyyy}} \;+\; \scalebox{.3}{\input{FigColored/threebubble_yyy}} = \tfrac{1}{128} + \tfrac{1}{48} + \tfrac{1}{16} = \tfrac{35}{384}. \qedhere
    \end{align*}
\end{example}

\begin{remark}[Ising model]
\label{rmk:ising}
Our examples are motivated from the physical Ising model.
The partition function of the \emph{Ising model} on a monochromatic graph $G$ is defined by
$$
Z(G,\lambda) = \sum_{\substack{\gamma \subset G\\\gamma \text{ Eulerian}}} \lambda^{|E(\gamma)|},
$$
where we sum over all Eulerian subgraphs $\gamma$ of $G$, i.e., subgraphs with all vertex degrees even. See, e.g., \cite{MR2989454} for a treatment of the Ising model on graphs. 
A pair $(G,\gamma)$ of a monochromatic graph $G$ and an Eulerian subgraph $\gamma \subset G$ is equivalent to an edge-bicolored graph in which an even number of \colT{} edges belongs to each vertex.

Notice that we effectively designed the polynomial $g(x,y)$ from \Cref{ex:Ising_model_deg4_part1}, {\parfillskip=0pt \par} \noindent and equivalently the coefficients $\Lambda_{u,w}$, such that the coefficient of $\lambda^k$ in $A_n$ is the {\parfillskip=0pt \par} \noindent  automorphism-weighted number of $4$-regular graphs with $k$ \colT{} edges where an even number of \colT{} edges belong to each vertex.
Hence, with $A_n$ as defined in \Cref{ex:Ising_model_deg4_part1}, we find that $A_n = \sum_{\substack{G}} \frac{Z(G,\lambda)}{|\Aut G|}$, where we sum over all monochromatic graphs $G$ that are $4$-regular 
and which have Euler characteristic $-n$. We can thus 
interpret $A_n$ as the partition function of the critical Ising 
model of a \emph{random} 4-regular monochromatic graph of fixed Euler characteristic. 
Here, random means that each monochromatic graph $G$ is sampled with probability $1/|\Aut G|$.
\end{remark}


\section{Asymptotics and critical points}
\label{sec:asymptotics}

In this section, we study the asymptotic behavior of the coefficients 
$A_n$ from \Cref{thm:laplace-graphs}, for large $n = -\chi(G)$. 
Here, we will restrict ourselves to 
\emph{regular} edge-bicolored graphs, 
meaning that each vertex has a fixed degree $k \geq 3.$
For fixed coefficients $\ol_{u,w}$ given for $u,w\geq 0$ with $u+w=k$, we study the weighted sum over graphs
\[
A_n = \sum_{G \in \GG_{-n}^k} 
        \frac{1}{
\left\lvert \Aut(G) \right\rvert} 
\prod \limits _{v\in V^G}
\ol_{\deg(v)},
\]
where $\GG_{-n}^k$ is the set of all regular (edge-bicolored) graphs with vertex degree $k$
and Euler characteristic $-n=|V^G|-|E^G|$.
As for each $k$-regular graph $G$ we have $k|V^G| = 2|E^G|$,
all graphs in $\GG_{-n}^k$ have $\frac{2n}{k-2}$ vertices and 
$\frac{n k}{k-2}$ edges. The negative Euler characteristic is also sometimes known as the \emph{excess} of the graph.


It is common belief in physics that the asymptotic behavior of $A_n$ depends on the critical points of $g(x,y)$ (see, e.g.,~\cite{le2012large}). This will be formalized in the next result. Notice that it is well-known, also in applied mathematics, that identifying the critical point that contributes most to the asymptotics is a complicated \emph{connection problem} \cite{berry1990hyperasymptotics}.

We write $\operatorname{crit}_D f$ for the set of critical points of $f$ restricted to the domain $D$. Let
\begin{equation}\label{eq:defPsi}
\Psi = \left\{ (x,y) \in \crit_{\C\cdot\R^2}g\!\setminus\! \{\mathbf{0}\} \,:\, \|(x,y)\|\leq \|(x',y')\| \, \forall (x',y') \in \crit_{\C\cdot\R^2}g\!\setminus\! \{\mathbf{0}\} \right\},
\end{equation}
where $\C\cdot\R^2$ is the set of complex points $(x,y)$ whose ratio (if well-defined) is real.
We call points in $\Psi$ \emph{non-degenerate} if the Hessian (the matrix of second derivatives) of $g$ has full rank. We are now ready to present the main result, which we already stated in \Cref{thm:main} and reformulate as follows.

\begin{theorem}\label{thm:crit_g}
    Let $\frac{z}{2\pi}\int_D \exp\left( {z\, g(x,y)} \right) \, \dd x \dd y\, \sim \sum_{n\geq 0} A_n z^{-n}$ and let $\Psi$ be defined as in \eqref{eq:defPsi}.
    Assume that $g(x,y) + \frac{x^2}{2} + \frac{y^2}{2}$ is homogeneous and that all extrema in $\Psi$ are non-degenerate. Then,
    \begin{equation}
    \label{eq:niceFormula}
    A_n \sim 
    \frac{1}{2\pi}
    \Gamma(n)
    \sum\limits_{(x,y) \in \Psi} \frac{(-g(x,y))^{-n}}{\sqrt{-\det H_g(x,y)}} \qquad \text{ as } \quad n \rightarrow \infty.
    \end{equation}
\end{theorem}

\begin{example}\label{ex:Ising_numerical_An}
    To continue the running example of the Ising model (cf. \Cref{ex:Ising_model_deg4_part1} and \ref{ex:Ising_model_deg4_part2}), let $g(x,y) = -\frac{x^2}{2} - \frac{y^2}{2} + \frac{x^4}{4!} + \lambda \frac{x^2 y^2}{4} + \lambda^2 \frac{y^4}{4!}$.
    Consider the system of critical equations for $g$
    \begin{equation}\label{eq:system_crit_pts}
        {\color{cb-green} \frac{\partial g}{\partial x}(x,y)=0}, \quad
        {\color{cb-purple} \frac{\partial g}{\partial y}(x,y)=0},
    \end{equation}
    shown in \Cref{fig:sing_various_lambdas_3D}, and its complex non-trivial solutions, for $\lambda>0$:
    \[
     (\pm \sqrt{6}, 0), \left(0,\pm \frac{\sqrt{6}}{\lambda}\right), \left( \pm \sqrt{\frac{9-3\lambda}{4\lambda}}, \pm \frac{\sqrt{9 \lambda -3}}{2\lambda} \right).
    \]
    \begin{figure}[ht]
        \centering
        \includegraphics[width=0.35\textwidth]{FigColored/gsurf.png}\\
        \input{FigColored/Ising_sing_g}
        \caption{The system \eqref{eq:system_crit_pts} for the function $g$ from \Cref{ex:Ising_numerical_An}: $\frac{\partial g}{\partial x}$ in green, $\frac{\partial g}{\partial y}$ in purple. Above: All values of $\lambda \in (\frac{1}{5},4)$ on the vertical axis. At each level $\lambda=\text{const.}$ the black points are the critical points of $g$.
        Below: Sections of the $3$D figure, for $\lambda = \frac{1}{4}, \frac{1}{2}, 4$. The non-trivial critical points closest to the origin are marked with squares. Notice the different scaling in the $y$-axis, for the sake of clarity.
        }
        \label{fig:sing_various_lambdas_3D}
    \end{figure}
    Among these solutions, some are real for every $\lambda>0$. 
    The last type of singular points is real or has a real ratio if and only if $\lambda \in [\frac{1}{3}, 3]$.
    Comparing their distances to the origin gives the set $\Psi$, as pictured in \Cref{fig:sing_various_lambdas_3D}.
    We can then use \Cref{thm:crit_g} to find
    $A_n \sim c \:\Gamma(n) \alpha^n$,
    where $c = c(\lambda)$ and $\alpha = \alpha(\lambda)$ are piecewise defined as
    \begin{center}
        {\def\arraystretch{1.8}
        \begin{tabular}{c|c|c|c}
             & $\Psi$ & $\alpha(\lambda)$ & $c(\lambda)$ \\
            \hline
           $0<\lambda<\tfrac{1}{3}$ & $(\pm\sqrt{6}, 0 )$ & $\frac{2}{3}$ & $\frac{1}{\pi} \sqrt{\frac{1}{2-6\lambda}}$ \\
           $\tfrac{1}{3}<\lambda<3$ & $\left(\pm \sqrt{\tfrac{9-3\lambda}{4\lambda}}, \pm \tfrac{\sqrt{9 \lambda -3}}{2\lambda} \right)$ & $\frac{-16\lambda^2}{3\lambda^2-18\lambda+3}$ & $\frac{1}{\pi} \sqrt{\frac{8\lambda}{-3\lambda^2+10\lambda -3}}$ \\
           $\lambda>3$ & $\left(0, \pm\tfrac{\sqrt{6}}{\lambda} \right)$ & $\frac{2\lambda^2}{3}$ & $\frac{1}{\pi} \sqrt{\frac{\lambda}{2\lambda -6}}$ \\
        \end{tabular}
        }
    \end{center}
    The function $\alpha$ is continuous, it is not $C^1$-differentiable at $\lambda=\frac{1}{3}$, and it is $C^1$- but not $C^2$-differentiable at $\lambda=3$. On the other hand, the limits of $c(\lambda)$ at $\frac{1}{3}$, $3$ go to infinity from both sides.
    This can be observed in \Cref{fig:Ising_deg4_alpha}.
    The points $\lambda=\frac13$ and $\lambda=3$, where the functions $\alpha(\lambda)$ and $c(\lambda)$ are non-analytic, 
    are \emph{phase transition} points.
    Phase transitions are of pivotal interest in statistical physics. Here, we find the phase transitions of the Ising model on a random 4-regular graph. In each of the three regions for the parameter $\lambda$, the statistical system is expected to exhibit intrinsically different behaviors. We postpone the study of such behaviors to future work.
    \begin{figure}[!ht]
        \centering
        \includegraphics[width = 0.45\textwidth]{FigColored/2dising.png}
        \quad
        \includegraphics[width = 0.45\textwidth]{FigColored/plot_c_Ising_transitions.png} 
        \caption{The behavior of $\alpha(\lambda)$ and $c(\lambda)$ in the Ising model from \Cref{ex:Ising_numerical_An}. 
        The phase transitions at $\lambda=\frac{1}{3}, 3$ can be detected in both quantities.}
        \label{fig:Ising_deg4_alpha}
    \end{figure}
\end{example}
Notice the significant difference between the univariate, monochromatic case and the bivariate, bicolored case. In the univariate case \cite[Chapter~3]{Borinsky:2018mdl}, the critical points of $g$ closest to the origin contribute to the asymptotic behavior of $A_n$ regardless of their reality. In contrast, the bivariate case imposes a restrictive condition: critical points must lie within $\C \cdot \R^2$. The difference seems to be an artifact of $\C \cdot \R = \C$.


\section{Outlook}\label{sec:outlook}
Although \Cref{thm:crit_g} assumes $g(x,y) + \frac{x^2}{2} + \frac{y^2}{2}$ to be homogeneous, numerical computations suggest that this condition is not necessary (see, e.g., \cite[Example 5.8]{BMW:EdgeBicolored}).
We conjecture that \Cref{thm:crit_g} is also valid in the inhomogeneous case, i.e., for graphs that are not necessarily regular. 
\begin{conjecture}
    Let $\frac{z}{2\pi}\int_D \exp\left( {z\, g(x,y)} \right) \, \dd x \dd y\, \sim \sum_{n\geq 0} A_n z^{-n}$
    and let $\Psi$ be defined as in \eqref{eq:defPsi}. Assume that all points in $\Psi$ are non-degenerate. Then, 
    $$
    A_n \sim 
    \frac{1}{2\pi}
    \Gamma(n)
    \sum\limits_{(x,y) \in \Psi} \frac{(-g(x,y))^{-n}}{\sqrt{-\det H_g(x,y)}}
    \qquad \text{ as } \quad n \rightarrow \infty.
    $$
\end{conjecture}

We remark that the only assumption that is left, namely non-degeneracy, could be removed via technical computation. Indeed, one would only have to adapt the asymptotic expression to higher-order Gaussian integrals, based on the lowest degree in the Taylor expansion of $g$ around points in $\Psi$, similarly to, e.g., \cite[Theorem VIII.10]{flajolet2009analytic}.

A second natural direction of research, which is the topic of ongoing work, is to generalize the results to a multicolor setting, where $g(x) = g(x_1,\ldots, x_m)$. In this framework, many more things can happen, such as the critical locus of $g$ having positive dimension. Therefore, it is necessary to appropriately generalize the concept of non-degeneracy, potentially in the direction of Morse--Bott theory, as well as to refine the notion of the real ratio of points in $\Psi$. 
% The very first step, however, can be stated in the case of regular multicolored graphs, where $g$ has only finitely many critical points. 
We pose the following challenge.

\begin{problem}
    Let 
    \[
    \left( \frac{z}{2\pi}\right)^{\frac{m}{2}}\int_D \exp\left( {z\, g(x)} \right) \, \dd x_1 \cdots \dd x_m\, \sim \sum_{n\geq 0} A_n z^{-n}.
    \]
    Find the set $\Psi \subset \crit g$ and a factor $c_m$ which possibly depends on $m$ such that
    \[
    A_n \sim 
    c_m \,
    \Gamma(n)
    \sum\limits_{x \in \Psi} \frac{(-g(x))^{-n}}{\sqrt{-\det H_g(x)}}
    \qquad \text{ as } \quad n \rightarrow \infty.
    \]
\end{problem}
% We can make the following conjecture.
% \begin{conjecture}
%     Let 
%     \[
%     \left( \frac{z}{2\pi}\right)^{\frac{m}{2}}\int_D \exp\left( {z\, g(x)} \right) \, \dd x_1 \cdots \dd x_m\, \sim \sum_{n\geq 0} A_n z^{-n}
%     \]
%     and 
%     \[
%     \Psi = \left\{ x \in \crit_{\C\cdot\R^m}g\!\setminus\! \{\mathbf{0}\} \,:\, \|x\|\leq \|x'\| \, \forall x' \in \crit_{\C\cdot\R^m}g\!\setminus\! \{\mathbf{0}\} \right\}.
%     \]
%     Assume that $g(x) + \sum_{i=1}^m \frac{x_i^2}{2}$ is homogeneous and that $\Psi$ consists of isolated non-degenerate (i.e., with full-rank Hessian) points. Then, 
%     $$
%     A_n \sim 
%     (2\pi)^{-\frac{m}{2}} \,
%     \Gamma(n)
%     \sum\limits_{x \in \Psi} \frac{(-g(x))^{-n}}{\sqrt{-\det H_g(x)}}
%     \qquad \text{ as } \quad n \rightarrow \infty.
%     $$
% \end{conjecture}

% \acknowledgements{TBD}

\printbibliography



\end{document}
