# Phases and Duality
in Fundamental Kazakov-Migdal Model
on the Graph

So Matsuura<EMAIL_ADDRESS>Hiyoshi Departments of Physics, and Research and
Education Center for Natural Sciences, Keio University, Yokohama, Kanagawa
223-8521, Japan Kazutoshi Ohta<EMAIL_ADDRESS>Institute
for Mathematical Informatics, Meiji Gakuin University, Yokohama, Kanagawa
244-8539, Japan

Abstract

We examine the fundamental Kazakov-Migdal (FKM) model on a generic graph,
whose partition function is represented by the Ihara zeta function weighted by
unitary matrices. The FKM model becomes unstable in the critical strip of the
Ihara zeta function. We discover a duality between small and large couplings,
associated with the functional equation of the Ihara zeta function for regular
graphs. Although the duality is not precise for irregular graphs, we show that
the effective action in the large coupling region can be represented by a
summation of all possible Wilson loops on the graph similar to that in the
small coupling region. We estimate the phase structure of the FKM model both
in the small and large coupling regions by comparing it with the Gross-Witten-
Wadia (GWW) model. We further validate the theoretical analysis through
detailed numerical simulations.

##  1 Introduction

Gauge theory serves as a fundamental framework for describing the interactions
in nature. Specifically, the property of the asymptotic freedom in the non-
Abelian gauge theory, which describes QCD, imposes limitations on analysis
solely through perturbation theory. As a result, non-perturbative approaches,
especially lattice gauge theories [1], have emerged as crucial tools. A
characteristic feature of lattice gauge theory is that it divides continuous
space-time into discrete lattices (graphs) and approximates infinite degrees
of freedom with finite degrees of freedom. This not only non-perturbatively
defines the theory but also enables the application of powerful numerical
analysis techniques.

In our recent series of studies [2, 3], we have developed a discretized gauge
theory on graphs as a generalization of the so-called Kazakov-Migdal (KM)
model [4]. The KM model on the graph has a mathematically remarkable property
that its partition function is represented by the graph zeta functions, that
is, graph-theoretical analogues of the Riemann zeta function [5, 6, 7, 8, 9,
10, 11] (For comprehensive reviews and references, see also [12]). This
enables the generalized KM model to accurately enumerate Wilson loops on the
graph and to precisely evaluate the partition function in the large $N_{c}$
limit. While the generalized KM model displays fascinating characteristics,
its application as a model for induced QCD — the foundational motivation
behind the KM model — faces the same difficulties as the original KM model due
to the presence of an undesirable local $U(1)$ symmetry [13, 14, 15, 16, 17,
18].

To construct a model closer to QCD, we have replaced the scalar fields of the
generalized KM model in the adjoint representation with ones in the
fundamental representation. We call this model the fundamental Kazakov-Migdal
(FKM) model [19], as a generalization of the model introduced in [20]. The FKM
model does not have the extra $U(1)$ symmetry with retaining the property that
the partition function is represented by the graph zeta function. As a result,
the effective action of the FKM model includes the conventional Wilson action
of gauge theory as a limit of certain theoretical parameters, suggesting the
model’s potential to describe a realistic gauge theory. We have analytically
shown that the FKM model on the cycle graph undergoes the Gross-Witten-Wadia
(GWW) third-order phase transition [21, 22], and have found theoretical
evidences for the existence of non-trivial phase transitions for the general
graphs, although the order of these phase transition is not so clear.

In this paper, we investigate the phase structure of the FKM model whose
partition function is described by the Ihara zeta function, which is the
ancestor of the graph zeta function, in more detail. Further theoretical and
numerical analysis yields the following interesting results on the phase
structure of the FKM model:

First of all, the FKM model on a regular graph, where each vertex has the same
degree, exhibits a duality between the small and large coupling regions. This
duality is a consequence of the functional equation of the Ihara zeta function
on the regular graph. Secondly, we discover a dual transformation applicable
to the Ihara zeta function of any graph, including irregular graphs, as a
generalization of the functional equation of the Ihara zeta function of the
regular graph. While the exact duality does not hold for the FKM model on
irregular graphs, the effective action in the large coupling region still can
be determined by using the dual transformation. In particular, the phase
structure of the FKM model, even in the large coupling region, can be
elucidated through the asymptotic expansion of the effective action,
especially when the number of “flavors” is large enough for the GWW model to
become a good approximation of the FKM model. Analytical and numerical
investigations of the phase structure of the FKM model on various graphs
reveals that it is closely related to the structure of the graph.

The organization of this paper is as follows:

In the next section, we introduce the notation of the graph theory and the
Ihara graph zeta function used in this paper, and then briefly review the
basic properties of the FKM model. In Sec. 3, we explain the functional
equation of the Ihara zeta function, which is the origin of the duality of the
FKM model. In Sec. 4, we show that there is an unstable parameter region in
the FKM model corresponding to the distribution of the poles of the Ihara zeta
function and explain the duality of the FKM model. After discussing the exact
results of the cycle graph, we theoretically analyze the phase structure of
the FKM model on some specific examples of both regular and irregular graphs.
In Sec. 5, we give the results of the numerical simulation for the regular
graphs. The exact analytical results for the cycle graph is used for the check
of the simulation code. In Sec. 6, we show the results of the numerical
simulation for the irregular graphs. The last section is devoted to the
conclusion and discussions.

##  2 Review of the Fundamental Kazakov-Migdal Model

###  2.1 Notations in the graph theory

The graph $G$ is an object which consists of vertices and edges that connect
two vertices. We denote the set of the vertices and edges by $V$ and $E$,
respectively, with their cardinalities represented as $n_{V}$ and $n_{E}$. We
assume that the edges have directions. Therefore, the term “graph” within this
paper refers to a directed graph (digraph). A directed edge $e$ starting from
a vertex $v$ and ending at a vertex $v^{\prime}$ is represented as $\langle
v,v^{\prime}\rangle$. The two vertices $v$ and $v^{\prime}$ of $e=\langle
v,v^{\prime}\rangle$ are called the “source” $v=s(e)$ and the “target”
$v^{\prime}=t(e)$ of $e$, respectively. The degree of a vertex $v$, denoted as
$\deg v$, is the number of its neighboring vertices. It is useful to consider
both a directed edge $e=\langle v,v^{\prime}\rangle$ and its opposite
$e^{-1}=\langle v^{\prime},v\rangle$. We then enhance the set of the edges $E$
to
$E_{D}=\\{{\boldsymbol{e}}_{a}|a=1,\cdots,2n_{E}\\}=\\{e_{1},\cdots,e_{n_{E}},e_{1}^{-1},\cdots
e_{n_{E}}^{-1}\\}$, combining with the opposite directed edges. In the
following, we assume that there is at most one edge connecting two vertices,
and no edge connecting the same vertex. In other words, we consider only
simple digraphs in this paper.

A path (or walk) of length $l$ on the graph $G$ is a sequence of $l$ edges,
$P={\boldsymbol{e}}_{a_{1}}{\boldsymbol{e}}_{a_{2}}\cdots{\boldsymbol{e}}_{a_{l}}$
satisfying the conditions
$t({\boldsymbol{e}}_{a_{i}})=s({\boldsymbol{e}}_{a_{i+1}})$ ($i=1,\cdots,l$
with $a_{l+1}=a_{1}$). The inverse of the path $P$ is defined as
$P^{-1}={\boldsymbol{e}}_{a_{l}}^{-1}\cdots{\boldsymbol{e}}_{a_{1}}^{-1}$. A
path $C$ whose starting and ending vertices are the same is called a cycle. We
can define the product of two paths
$P_{1}={\boldsymbol{e}}_{a_{1}}\cdots{\boldsymbol{e}}_{a_{l_{1}}}$ and
$P={\boldsymbol{e}}_{a^{\prime}_{1}}\cdots{\boldsymbol{e}}_{a^{\prime}_{l_{2}}}$
when $t({\boldsymbol{e}}_{a_{l}})=s({\boldsymbol{e}}_{a^{\prime}_{1}})$ as
$P_{1}P_{2}={\boldsymbol{e}}_{a_{1}}\cdots{\boldsymbol{e}}_{a_{l_{1}}}{\boldsymbol{e}}_{a^{\prime}_{1}}\cdots{\boldsymbol{e}}_{a^{\prime}_{l_{2}}}$.
In particular, a power of a cycle $C^{n}$, which is constructed through this
product rule, is also a cycle. In this paper, we consider only connected
graphs, that is, graphs where any pair of vertices is connected by a path. We
also assume that the degrees of all vertices are more than or equal to two and
the graph is simple, that is, it does not contain more than one edge between
two vertices and does not contain any edge that connects a vertex to itself
(loop).

A part of a path of length $l$ satisfying
${\boldsymbol{e}}_{a_{j+1}}={\boldsymbol{e}}_{a_{j}}^{-1}$ ($j=1,\ldots,l-1$)
is called a backtracking. If a cycle of length $l$ satisfies
${\boldsymbol{e}}_{a_{l}}={\boldsymbol{e}}_{a_{1}}^{-1}$, this part is called
a tail of the cycle $C$. For a cycle $C$, the equivalence class $[C]$ is
defined by the set of cyclic rotations of $C$ as

$[C]=\\{{\boldsymbol{e}}_{a_{1}}{\boldsymbol{e}}_{a_{2}}\cdots{\boldsymbol{e}}_{a_{l}},\,{\boldsymbol{e}}_{a_{2}}\cdots{\boldsymbol{e}}_{a_{l}}{\boldsymbol{e}}_{a_{1}},\,\ldots\,,{\boldsymbol{e}}_{a_{l}}{\boldsymbol{e}}_{a_{1}}\cdots{\boldsymbol{e}}_{a_{l-1}}\\}\,.$

For the equivalence class of the cycle $[C]$, the backtracking and the tail
are identical, hence they are collectively called a bump.

A cycle $C$ is called reduced when it contains no bump. A cycle $C$ is called
primitive when $C$ does not satisfy $C\neq B^{r}$ for any cycle $B$ and $r\geq
2$. We denote the set of representatives of reduced cycles as
$[{\mathcal{P}}_{R}]\subset[{\mathcal{P}}]$. Since a primitive reduced cycle
$C$ has its inverse $C^{-1}$ also be a primitive reduced cycle of equal
length, the set $[{\mathcal{P}}_{R}]$ can be partitioned into two disjoint
unions; $[{\cal P}_{R}]=[\Pi_{+}]\sqcup[\Pi_{-}]$, where $[\Pi_{-}]$ consists
of the inverses of elements in $[\Pi_{+}]$. These elements in $[\Pi_{+}]$ are
referred to as (the representatives of) chiral primitive reduced cycles.

###  2.2 Ihara zeta functions

The Ihara zeta function of a graph $G$ is defined by a product over the
primitive reduced cycles on $G$ [5, 6, 7]

$\zeta_{G}(q)\equiv\prod_{C\in[{\cal P}_{R}]}\frac{1}{1-q^{|C|}}\,,$ (2.1)

which is an analog to the Euler product expression of the Riemann zeta
function. In the above expression, $|C|$ denotes the length of the cycle $C$.
One of the important properties of the Ihara zeta function is that it is
expressed as the reciprocal of a polynomial

$\zeta_{G}(q)=\bigl{(}1-q^{2}\bigr{)}^{-(n_{E}-n_{V})}\det\bigl{(}\boldsymbol{1}_{n_{V}}-qA+q^{2}Q\bigr{)}^{-1}\,,$
(2.2)

where a square matrix $A$ of size $n_{V}$ is called the adjacency matrix
defined by

$A_{vv^{\prime}}=\sum_{{\boldsymbol{e}}\in
E_{D}}\delta_{\langle{v},{v^{\prime}}\rangle,{\boldsymbol{e}}}\quad(v,v^{\prime}\in
V)\,,$ (2.3)

and $Q$ is a diagonal matrix given by

$Q\equiv\operatorname{diag}\left(t_{1},t_{2},\cdots,t_{n_{V}}\right)\,,$ (2.4)

with $t_{i}\equiv\deg v_{i}-1$. The Ihara zeta function (2.1) has another
expression introduced by Hashimoto [8, 9],

$\zeta_{G}(q)=\det\left(\boldsymbol{1}_{2n_{E}}-qW\right)^{-1}\,,$ (2.5)

where a square matrix $W$ of size $2n_{E}$ is called the edge adjacency matrix
whose elements are given by

$\displaystyle
W_{{\boldsymbol{e}}{\boldsymbol{e}}^{\prime}}=\begin{cases}1&{\rm if}\
t({\boldsymbol{e}})=s({\boldsymbol{e}}^{\prime})\ {\rm and}\
{\boldsymbol{e}}^{\prime-1}\neq{\boldsymbol{e}}\\\ 0&{\rm
others}\end{cases}\,,$ (2.6)

where ${\boldsymbol{e}},{\boldsymbol{e}}^{\prime}\in E_{D}$.

###  2.3 The FKM model

On the graph $G$, we can construct a gauge theory by putting the scalar field
$\Phi^{I}_{v}$ on the vertices and the gauge fields $U_{e}$ on the edges. The
edge variables $U_{e}$ are $N_{c}\times N_{c}$ unitary matrices and the scalar
fields $\Phi^{I}_{v}$ transform in the fundamental representation of
$U(N_{c})$, where the index $I=1,2,\cdots,N_{f}$ labels the “flavors”. Using
these fields, we define a model by the action,

$S=\sum_{v\in V}m_{v}^{2}\Phi_{v}^{\dagger I}\Phi_{vI}-q\sum_{e\in
E}\left(\Phi_{s(e)}^{\dagger\ \ I}\,U_{e}\,\Phi_{t(e)I}+\Phi_{t(e)}^{\dagger\
\ I}\,U_{e}^{\dagger}\,\Phi_{s(e)I}\right)\,,$ (2.7)

where $q$ is a coupling constant111 Although $q$ is called a coupling constant
in this paper, we should note that it becomes the inverse of the gauge
coupling constant when we compare the model with the GWW model. Therefore, the
small/large region of $q$ in this paper corresponds to the strong/weak
coupling region in the sense of the gauge coupling.  and $m_{v}$ are the
masses for the scalar fields. We call this model the fundamental Kazakov-
Migdal (FKM) model introduced in [19]. Note that the term proportional to $q$
in the action represents the neighborhood interaction of $\Phi_{v}^{I}$.

In contrast to the traditional Kazakov-Migdal model, where the scalar fields
are in the adjoint representation and the effective action of the adjoint
scalar fields rather than the edge variables is obtained by using the Harish-
Chandra-Itzykson-Zuber integral [4], we first integrate out the scalar fields
$\Phi_{v}^{I}$ in the partition function,

$Z_{G}=\int\prod_{I=1}^{N_{f}}\prod_{v\in
V}d\Phi_{vI}d\Phi^{{\dagger}vI}\prod_{e\in E}dU_{e}\,e^{-S}\ .$ (2.8)

Then, we obtain

$Z_{G}=\left(2\pi\right)^{N_{c}N_{f}n_{V}}\int\prod_{e\in
E}dU_{e}\,\det\left(m_{v}^{2}{\delta_{v}}^{v^{\prime}}{\bf
1}_{N_{c}}-q{(A_{U})_{v}}^{v^{\prime}}\right)^{-N_{f}}\ ,$ (2.9)

where $A_{U}$ is a matrix of size $N_{c}n_{V}$ which is essentially the
adjacency matrix (2.3) but is weighted by unitary matrices as

${(A_{U})_{v}}^{v^{\prime}}=\sum_{{\boldsymbol{e}}\in
E_{D}}U_{{\boldsymbol{e}}}\delta_{\langle
v,v^{\prime}\rangle,{\boldsymbol{e}}}\ .$ (2.10)

Note that $U_{e^{-1}}$ is regarded as $U_{e}^{{\dagger}}$ with respect to the
unitary matrix $U_{e}$ along the original direction.

The integrand of the partition function (2.9) can be expanded in terms of the
various kind of the Wilson loops along the cycles on $G$. This is a reason why
this model appears to have the potential to induce QCD [20], but it is
difficult to control the coupling constant of the Wilson loops systematically
because of the mismatch between the order of the mass parameter $m_{v}$ and
the net length of the Wilson loop (cycles).

This difficulty is resolved by tuning the mass parameter $m_{v}$, following
[19], as

$m_{v}^{2}=1+q^{2}(\deg v-1)=1+q^{2}t_{v}\,.$ (2.11)

Then, the partition function (2.9) can be rewritten through the unitary matrix
weighted Ihara zeta function [2, 3]222$Q$ has been redefined to
$Q=\operatorname{diag}(t_{1}\boldsymbol{1}_{N_{c}},t_{2}\boldsymbol{1}_{N_{c}},\cdots,t_{n_{V}}\boldsymbol{1}_{N_{c}})$,
but we believe there is no confusion to the matrix $Q$ used to express the
original Ihara zeta function (2.2). ,

$\zeta_{G}(q;U)=\bigl{(}1-q^{2}\bigr{)}^{-{N_{c}}(n_{E}-n_{V})}\det\bigl{(}\boldsymbol{1}_{{N_{c}}n_{V}}-qA_{U}+q^{2}Q\bigr{)}^{-1}\,,$
(2.12)

as

$\displaystyle Z_{G}$
$\displaystyle=\left({2\pi}\right)^{{N_{f}}{N_{c}}n_{V}}\left(1-q^{2}\right)^{{N_{f}}{N_{c}}(n_{E}-n_{V})}\int\prod_{e\in
E}dU_{e}\,\zeta_{G}(q;U)^{N_{f}}\ .$ (2.13)

Similar to the original Ihara zeta function (2.1), the matrix weighted Ihara
zeta function can be expressed in the Euler product form,

$\zeta_{G}(q;U)\equiv\prod_{C\in[{\cal
P}_{R}]}\frac{1}{\det\bigl{(}\boldsymbol{1}_{{N_{c}}}-q^{|C|}U_{C}\bigr{)}}\
.$ (2.14)

Expanding the product (2.14) by the Wilson loops $U_{C}$, we can exactly count
the length of the cycles in the powers of $q$ as

$\displaystyle Z_{G}$
$\displaystyle=\left({2\pi}\right)^{{N_{f}}{N_{c}}n_{V}}\left(1-q^{2}\right)^{{N_{f}}{N_{c}}(n_{E}-n_{V})}$
$\displaystyle\qquad\times\int\prod_{e\in
E}dU_{e}\,\prod_{C\in[\Pi_{+}]}\exp\left\\{N_{f}\operatorname{Tr}\left[\sum_{n=1}^{\infty}\frac{q^{n|C|}}{n}\left(U_{C}^{n}+U_{C}^{\dagger\,n}\right)\right]\right\\}\
.$ (2.15)

Note that we can generalize the parametrization (2.11) as

$m_{v}^{2}=1-q^{2}(1-u)^{2}+q^{2}(1-u)\deg v\,,$ (2.16)

so that the partition function is expressed by the matrix weighted Bartholdi
zeta function [10], which is one parameter extension of the Ihara zeta
function. However, we focus on the Ihara zeta function in order to simplify
the following discussion in this paper. The generalization of the analysis of
this paper to the Bartholdi zeta function will be discussed separately.

From the partition function (2.15), the effective action of the model can be
read off as

$\displaystyle S_{\rm eff}(q;U)$ $\displaystyle=-N_{f}\log\zeta_{G}(q;U)$
$\displaystyle=-\gamma
N_{c}\operatorname{Tr}\left[\sum_{C\in[\Pi_{+}]}\sum_{n=1}^{\infty}\frac{q^{n|C|}}{n}\left(U_{C}^{n}+U_{C}^{\dagger\,n}\right)\right]\,,$
(2.17)

where we have defined $\gamma\equiv N_{f}/N_{c}$ and ignored an irrelevant
constant term. This action is a generalization of the Wilson action of the
lattice gauge theory on the graph $G$ in the sense that it contains infinitely
many different Wilson loops controlled by the Ihara zeta function. So it is
interesting to investigate the dynamics and phases of this complicated model
including the higher order interactions.

We note that the unitary matrix weighted Ihara zeta function (2.14) also has
the Hashimoto expression [8] similar to (2.5),

$\displaystyle\zeta_{G}(q;U)$
$\displaystyle=\det\left(\boldsymbol{1}_{2N_{c}n_{E}}-qW_{U}\right)^{-1}\,,$
(2.18)

where the matrix $W_{U}$ is an extension of the edge adjacency matrix (2.6)
weighted by unitary matrices,

$\displaystyle(W_{U})_{{\boldsymbol{e}}{\boldsymbol{e}}^{\prime}}=\begin{cases}U_{{\boldsymbol{e}}}&{\rm
if}\ t({\boldsymbol{e}})=s({\boldsymbol{e}}^{\prime})\ {\rm and}\
{\boldsymbol{e}}^{\prime-1}\neq{\boldsymbol{e}}\\\ 0&{\rm
others}\end{cases}\,.$ (2.19)

The equivalence of both expressions (2.14) and (2.19) has been proved in [3]
by using the similar technique given by Bass [9].

###  2.4 Gauge fixing and the path representation

So far, we have considered unitary matrices on all the edges, but we can fix
the gauge as $U_{e}=1$ on a spanning tree (maximal tree) [23]. The spanning
tree $T$ of the graph $G$ is a tree subgraph of $G$ which does not contain any
closed cycle and includes all of the vertices. If we remove the spanning tree
$T$ from the graph $G$, it consists of the disconnected pieces of the edges.
After fixing the gauge, the path integral is taken over the remaining
$U_{a}$’s where the edge label $a$ runs over the edges on $G-T$. We call the
the number of the edges on $G-T$ as the rank $r$ of the graph $G$.

To perform the path integral over $U_{a}$ on $G-T$, it is useful to define the
path representation of the Ihara zeta function. To this end, we need to
introduce some additional notations: As we have explained, the index
$a=1,\cdots,r$ corresponds to the edges in $G-T$. Corresponding to the inverse
edges of them, we introduce the index $\bar{a}=r+1,\cdots,2r$ and we express
them collectively as ${\boldsymbol{a}}\equiv(a,\bar{a})$ which runs from $1$
to $2r$. From this definition, it is natural to set $U_{\bar{a}}\equiv
U_{a}^{\dagger}$. Since $T$ is a tree graph, such a reduced path
$P_{{\boldsymbol{a}}{\boldsymbol{b}}}$ that starts from $e_{\boldsymbol{a}}\in
G-T$ and ends at $s(e_{\boldsymbol{b}})$ is unique and contains only the edges
in $T$ except for $e_{\boldsymbol{a}}$. Using these notations, the path
representation of the Ihara zeta function is given by [24]

$\zeta_{G}(q;U)=\det\left(1-P_{U}\right)\,,$ (2.20)

where $P_{U}$ is a square matrix of size $2r$ whose elements are given by

$\left(P_{U}\right)_{{\boldsymbol{a}}{\boldsymbol{b}}}=\begin{cases}q^{|P_{{\boldsymbol{a}}{\boldsymbol{b}}}|}U_{{\boldsymbol{a}}}&({\boldsymbol{b}}\neq\bar{{\boldsymbol{a}}})\\\
0&({\boldsymbol{b}}=\bar{{\boldsymbol{a}}})\end{cases}\,.$ (2.21)

Comparing to the vertex representation (2.2) and the edge representation
(2.5), where the zeta function is described by the determinant of sparse
matrices, the path representation is efficient for numerical simulations
because the matrix $P$ is dense.

##  3 The Functional Equation of the Ihara Zeta Function and the Riemann
Hypothesis

In this section, we show that there is a dual expression of the Ihara zeta
function of any connected graph in general, followed by recalling the
functional equation of the Ihara zeta function of a regular graph. We will use
the results in this section to show a duality in the FKM model in the
following section.

###  3.1 Functional equation of the Ihara zeta function of a regular graph

Before discussing the functional equation of the Ihara zeta function, it is
instructive to recall the more familiar Riemann zeta function, since the Ihara
zeta function has many analogies to it. The Riemann zeta function is defined
by

$\zeta(s)=\sum_{n=1}^{\infty}\frac{1}{n^{s}}\ ,$ (3.1)

and it is shown that $\zeta(s)$ satisfies the following functional equation,

$\zeta(s)=2^{s}\pi^{s-1}\sin\left(\frac{\pi s}{2}\right)\Gamma(1-s)\zeta(1-s)\
.$ (3.2)

If we define the so-called the completed zeta function ($\xi$-function) as

$\xi(s)\equiv\frac{1}{2}\pi^{-\frac{s}{2}}s(s-1)\Gamma\left(\frac{s}{2}\right)\zeta(s)\,,$
(3.3)

which satisfies the symmetric functional equation,

$\xi(s)=\xi(1-s)\,.$ (3.4)

This reflective symmetry is a key to understand that non-trivial zeros lies at
least in the open strip of $0<\operatorname{Re}s<1$, which is called the
critical strip. Furthermore, the conjecture that there exists non-trivial
zeros only on the line along $\operatorname{Re}s=\frac{1}{2}$ (critical line)
is called the Riemann hypothesis.

The Ihara zeta function of a regular graph enjoys a similar relation to the
functional equation of the Riemann zeta function. The regular graph is such a
graph that all of the degree of all vertices are the same. If we set the
common degree by $\deg v=t+1$, the Ihara zeta function (2.2) reduces to

$\zeta_{G}(q)=\left(1-q^{2}\right)^{-n_{E}+n_{V}}\det\left((1+tq^{2})\boldsymbol{1}_{n_{V}}-qA\right)^{-1}\
.$ (3.5)

Similar to the Riemann zeta function, we can show that the completed Ihara
zeta function,

$\displaystyle\xi_{G}(q)$
$\displaystyle\equiv\left(1-q^{2}\right)^{n_{E}-\frac{n_{V}}{2}}\left(1-t^{2}q^{2}\right)^{\frac{n_{V}}{2}}\zeta_{G}(q)$
$\displaystyle=\left(1-q^{2}\right)^{\frac{n_{V}}{2}}\left(1-t^{2}q^{2}\right)^{\frac{n_{V}}{2}}\det\left((1+tq^{2})\boldsymbol{1}_{n_{V}}-qA\right)^{-1}\,,$
(3.6)

satisfies the symmetric functional equation,

$\xi_{G}(q)=(-1)^{n_{V}}\xi_{G}\left({1}/{tq}\right)\,.$ (3.7)

The proof is straightforward from the expression (3.6):

$\displaystyle\xi_{G}(1/tq)$
$\displaystyle=\left(1-\left(tq\right)^{-2}\right)^{\frac{n_{V}}{2}}\left(1-{q^{-2}}\right)^{\frac{n_{V}}{2}}\det\left(\left(1+(tq^{2})^{-1}\right)\boldsymbol{1}_{n_{V}}-(tq)^{-1}A\right)^{-1}$
$\displaystyle=\frac{(t^{2}q^{2}-1)^{\frac{n_{V}}{2}}(q^{2}-1)^{\frac{n_{V}}{2}}}{(tq^{2})^{n_{V}}}\times(tq^{2})^{n_{V}}\det\left((tq^{2}+1)\boldsymbol{1}_{n_{V}}-qA\right)^{-1}$
$\displaystyle=(-1)^{n_{V}}\xi_{G}\left({q}\right)\,.$ (3.8)

If we change the variable $q$ with

$q=t^{-s}\,,$ (3.9)

the functional equation (3.7) reduces to

$\xi_{G}(s)=(-1)^{n_{V}}\xi_{G}\left(1-s\right).$ (3.10)

This is nothing but an analog of the symmetric functional equation (3.4) of
the Riemann zeta function.

For the $(t+1)$-regular graph, the radius of the largest circle of convergence
$R_{G}$ with respect to $q$ is given by $1/t$ and all the poles exist in the
region $1/t\leq|q|\leq 1$, that is, $0\leq s\leq 1$ which is consistent with
the functional relation (For a proof, see Appendix A)333 In the Ihara zeta
function, we ara interested in the non-trivial poles instead of the non-
trivial zeros. . It is also interesting that the Ihara zeta function
$\zeta_{G}(q)$ of the $(t+1)$-regular graph has the non-trivial poles only on
the critical line ($\operatorname{Re}s=1/2$), namely it satisfies the Riemann
hypothesis, iff $G$ is Ramanujan, that is, the maximal eigenvalue
$\lambda_{\rm max}$ except for $\lambda=t+1$ of the adjacency matrix $A$
satisfies444 Since the adjacency matrix $A$ satisfies
$\sum_{v^{\prime}=1}^{n_{V}}A_{vv^{\prime}}=\deg v$, $A$ for a $(t+1)$-regular
graph has a trivial eigenvector $(1,\cdots,1)^{T}$ with the eigenvalue $t+1$.
This means that the Ihara zeta function of the $(t+1)$-regular graph has
always poles at $q=1$ and $q=1/t$.

$|\lambda_{\rm max}|\leq 2\sqrt{t}\,.$ (3.11)

In contrast to that the original Riemann hypothesis of the Riemann zeta
function is highly non-trivial, we can easily prove it as follows:

Since the Ihara zeta function of the $(t+1)$-regular graph is expressed as
(3.5), the non-trivial poles of $\zeta_{G}(q)$ is written down by using the
eigenvalue $\lambda$ of the adjacency matrix $A$ as $q_{\pm}\equiv
t^{-s_{\pm}}=\frac{\lambda\pm\sqrt{\lambda^{2}-4t}}{2t}$, which obey
$t^{-s_{+}}t^{-s_{-}}=t^{-1}$ in general. We assume that $\lambda$ is not
$t+1$ which is a trivial eigenvalue of $A$. When $\lambda^{2}-4t<0$, $s_{\pm}$
become complex conjugate with each other and these solutions satisfy
$t^{-s_{+}}t^{-s_{-}}=t^{-2\operatorname{Re}s_{\pm}}$, that means the non-
trivial poles lie along the line of $\operatorname{Re}s_{\pm}=1/2$ on
$s$-plane. $\Box$

(a) Tetrahedron $K_{4}$

(b) Double triangle $K_{4}-e$

Figure 1: Two graphs with four vertices. The double triangle $K_{4}-e$ (DT) is
obtained by removing one edge from the tetrahedron $K_{4}$ of the complete
graph. $K_{4}$ is the regular graph with the same degree 3 for each vertex,
but DT is not regular.

As a concrete example of a regular graph, let us consider the complete graph
on four vertices $K_{4}$, that is, the tetrahedron graph (see Fig. 1 (a)). The
Ihara zeta function is given by

$\zeta_{K_{4}}(q)=\frac{1}{(1-q)(1-2q)(1-q^{2})^{2}(1+q+2q^{2})^{3}}\,,$
(3.12)

and the convergence radius is $R_{K_{4}}=1/2$ since $t=2$ in this case. We see
that the eigenvalues of the adjacency matrix are $(3,-1,-1,-1)$, and the
maximal eigenvalue except for $t+1=3$ is $\lambda_{\rm max}=-1$ which
satisfies $|\lambda|\leq 2\sqrt{2}$. This means that the $K_{4}$ graph is
Ramanujan. The poles of the zeta function (3.12) exist at $q=\pm 1,-\frac{1\pm
i\sqrt{7}}{4},\frac{1}{2}$. Putting these poles on the $s$-plane as in Fig. 2
(a), we see that the non-trivial poles ($q=-\frac{1\pm i\sqrt{7}}{4}$ and
$\frac{1}{2}$) lie along the critical line $\operatorname{Re}s=1/2$, that is,
the Riemann hypothesis is satisfied.

(a) Tetrahedron $K_{4}$

(b) Double triangle $K_{4}-e$

Figure 2: The poles of the Ihara zeta function of the $K_{4}$ and $K_{4}-e$
(DT) graph on the $s$-plane. The $K_{4}$ case satisfies the Riemann hypothesis
while $K_{4}-e$ does not, but all the poles exist in the critical strip with
the boundary $0\leq\operatorname{Re}s\leq 1$.

###  3.2 Relations in the Ihara zeta function of an irregular graph

We next consider the Ihara zeta function of an irregular graph $G$. In this
case, the value of the convergence radius $R_{G}$ must be calculated
explicitly for each graph, but the lower and the upper bound of $R_{G}$ is
known as [25]

$\frac{1}{t_{\text{max}}}\leq R_{G}\leq\frac{1}{t_{\text{min}}},$ (3.13)

where $t_{\text{min}}=\min_{v\in V}(\deg v-1)$ and $t_{\text{max}}=\max_{v\in
V}(\deg v-1)$ (See Appendix A for a proof). Note that, if $G$ is the
$(t+1)$-regular graph, $R_{G}={1}/{t}$ is reproduced from this relation
because $t_{\rm max}=t_{\rm min}=t$. Since we can show that all the poles of
$\zeta_{G}(q)$ exist in the range $R_{G}\leq|q|\leq 1$ also for an irregular
graph (see Appendix A), it is reasonable to define a parameter $s$ as

$q\equiv R_{G}^{s}\,,$ (3.14)

corresponding to (3.9) for the regular graph. Then, all the poles again exist
in the critical strip with the boundary $0\leq\operatorname{Re}s\leq 1$ in
this parametrization as well as for regular graphs.

The most important difference from the regular graph is that there is no
functional relation for the graph zeta function of the irregular graph.
Nevertheless, we can show that the Ihara zeta function with $q$ replaced by
the inverse is again essentially the reciprocal of a polynomial of $q$ up to a
non-essential factor. This can be seen from the fact that the Ihara zeta
function satisfies the following identity in general

$\displaystyle\zeta_{G}\left(R_{G}/q\right)^{-1}$
$\displaystyle=\frac{\left((q/R_{G})^{2}-1\right)^{n_{E}-n_{V}}\prod_{v\in
V}(\deg v-1)}{(q/R_{G})^{2n_{E}}}$ $\displaystyle\hskip
56.9055pt\times\det\left({\bf
1}_{n_{V}}-q(R_{G}Q)^{-1}A+R_{G}^{-1}q^{2}(R_{G}Q)^{-1}\right)\,.$ (3.15)

The proof is straightforward as follows:

By changing the argument of the Ihara zeta function (2.2) $q\to 1/cq$ with a
parameter $c\neq 0$, we obtain

$\displaystyle\zeta_{G}(1/cq)^{-1}$
$\displaystyle=\bigl{(}1-(cq)^{-2}\bigr{)}^{n_{E}-n_{V}}\det\bigl{(}\boldsymbol{1}_{n_{V}}-(cq)^{-1}A+(cq)^{-2}Q\bigr{)}$
$\displaystyle=\bigl{(}1-(cq)^{-2}\bigr{)}^{n_{E}-n_{V}}\left(\det((cq)^{-2}Q)\right)\det\bigl{(}\boldsymbol{1}_{n_{V}}-cqQ^{-1}A+(cq)^{2}Q^{-1})\bigr{)}$
$\displaystyle=\frac{\left(c^{2}q^{2}-1\right)^{n_{E}-n_{V}}\prod_{v\in
V}(\deg v-1)}{(cq)^{2n_{E}N_{c}}}\det\left({\bf
1}_{n_{V}N_{c}}-cqQ^{-1}A+c^{2}q^{2}Q^{-1}\right)\,,$ (3.16)

which reduces to (3.15) by setting $c=1/R_{G}$ as a special case of the
identity. Note that $Q$ is invertible since we assume that the degrees are
more than or equal to two. $\square$

Note that, if the graph is $(t+1)$-regular, the identity (3.15) becomes

$\displaystyle\zeta_{G}\left(1/tq\right)^{-1}$
$\displaystyle=\frac{(1-t^{-2}q^{-2})^{n_{E}-\frac{n_{V}}{2}}(1-q^{-2})^{\frac{n_{V}}{2}}}{(1-q^{2})^{n_{E}-\frac{n_{V}}{2}}(1-t^{2}q^{2})^{\frac{n_{V}}{2}}}\zeta_{G}(q)^{-1}\,,$
(3.17)

which is equivalent to the functional relation (3.7).

The structure of the right-hand side of (3.15) is essentially the same as the
vertex representation of the Ihara zeta function (2.12). In this sense, (3.15)
is a generalization of the functional relation of the Ihara zeta function.
This transformation makes us possible to expand the graph zeta function by
$1/q$ and we will use it to show an approximate duality of the FKM model on
irregular graphs in the next section.

As an example of the irregular graph, let us see the double triangle (DT)
which is obtained by removing one edge from the tetrahedron $K_{4}$ (Fig. 1
(b)). For DT, the Ihara zeta function becomes

$\zeta_{\rm
DT}(q)=\frac{1}{\left(1-q^{4}\right)\left(1+q^{2}-2q^{3}\right)\left(1-q^{2}-2q^{3}\right)}\
.$ (3.18)

From this expression, we find $R_{\rm DT}=0.65729\cdots$, which satisfies
$1/t_{\rm max}<R_{\rm DT}<1/t_{\rm min}$ with $t_{\text{min}}=1$ and
$t_{\text{max}}=2$. We plot the poles in the $s$-plane in Fig. 2 (b). We see
that all the poles exist in the critical strip with the boundary
$0\leq\operatorname{Re}s\leq 1$ as expected, but they clearly do not satisfy
the Riemann hypothesis.

##  4 Duality of the FKM Model and the GWW Phase Transitions

###  4.1 Stability of the FKM model

We first would like to discuss the stability of the model on a general graph.

Since the action of the FKM model (2.7) is Gaussian with respect to
$\Phi_{v}^{I}$, the matrix $\Delta(q;U)\equiv{\bf
1}_{N_{c}n_{V}}-qA_{U}+q^{2}Q$ must be positive definite in order for the
system to be stable. Since we are interested in the stability around the
vacuum, namely, $U_{a}={\bf 1}_{N_{c}}$ ($a=1,\cdots,r$) after fixing the
gauge, we need only to examine the positive definiteness of the matrix of size
$n_{V}$,

$\Delta(q)\equiv\Delta(q;U=1)={\bf 1}_{n_{V}}-qA+q^{2}Q\,,$ (4.1)

which is called the deformed graph Laplacian.

Let us consider the behavior of the eigenvalues $\lambda_{i}(q)$
($i=1,\cdots,n_{V}$) of $\Delta(q)$ for $q>0$. The point is that the product
of the eigenvalues is essentially expressed by the inverse of the Ihara zeta
function,

$\det\Delta(q)=\prod_{i=1}^{n_{V}}\lambda_{i}=(1-q^{2})^{-n_{E}+n_{V}}\zeta_{G}(q)^{-1}\,.$
(4.2)

Let us first consider the region $0<q<R_{G}$. Since $\Delta(q)\to{\bf
1}_{n_{V}}$ in the limit of $q\to 0$, all the eigenvalues
$\\{\lambda_{i}(q)\\}$ ($i=1,\cdots,n_{V}$) satisfy $\lim_{q\to
0}\lambda_{i}(q)=1>0$. According to Kotani-Sunada’s theorem [25],
$\zeta_{G}(q)$ has a first pole at $q=R_{G}<1$ which has the smallest
magnitude of all poles. This means that the sign of one of the eigenvalues of
$\Delta(q)$ (we call it $\lambda_{1}(q)$) flips at $q=R_{G}$. Therefore, the
vacuum $U_{a}=\boldsymbol{1}_{N_{c}}$ is stable in $0<q<R_{G}$ while it
becomes unstable at $q=R_{G}$.

On the other hand, since $\lim_{q\to\infty}\Delta(q)=q^{2}Q$, the eigenvalues
$\lambda_{i}(q)$ $(i=1,\cdots,n_{V})$ asymptotically behave as
$\lim_{q\to\infty}\lambda_{i}(q)=q^{2}(\deg v-1)$ for ${}^{\exists}v\in V$.
Therefore, all the eigenvalues $\lambda_{i}(q)$ are positive for $q>1$ because
$\zeta_{G}(q)$ has no pole in $q>0$, that is, the vacuum
$U_{a}=\boldsymbol{1}_{N_{c}}$ is stable in $q>1$.

The matrix $\Delta(q)$ at $q=1$,

$\Delta(1)=\boldsymbol{1}_{n_{V}}-A+Q\,,$ (4.3)

is nothing but the graph Laplacian. As discussed in [26], the graph Laplacian
of a connected graph has only one zero mode, $\left(1,\cdots,1\right)^{T}$,
since it is expressed by the incidence matrix555 This expression of the
incidence matrix is possible since we assume that the graph is simple. ,

$L^{e}_{\ v}=\begin{cases}1&v=t(e)\\\ -1&v=s(e)\\\ 0&{\rm
others}\end{cases}\,,$ (4.4)

as

$\Delta=L^{T}L\,,$ (4.5)

and $\dim\ker L=1$ for a connected graph (See also [27, 28, 29]). Therefore,
one of the eigenvalues becomes zero at $q=1$, namely the FKM model becomes
unstable at $q=1$. Combining these results, we can conclude that the vacuum
$U_{a}={\bf 1}_{N_{c}}$ is stable in the regions $0<q<R_{G}$ and $q>1$, but it
is unstable at least in the vicinity of the boundaries of $R_{G}<q<1$.

Whether there are other stable regions in $R_{G}<q<1$ or not depends on the
positions of the real poles of $\zeta_{G}(q)$. If $\zeta_{G}(q)$ has no real
pole in $0<q<1$ except for $q=R_{G}$, the corresponding FKM model is unstable
in the whole region of $R_{G}<q<1$, namely in the critical strip, because
there is no chance to flip the sign of the eigenvalue $\lambda_{1}(q)$ other
than $q=R_{G}$ and $q=1$. Since all of the graphs we consider in the next
section have this property, the FKM models we numerically simulate in this
paper are unstable in the critical strip. However, this is not always the
case, as real poles can generally exist in the region $R_{G}<q<1$. For graphs
that have real poles in the range $0<q<1$ other than $q=R_{G}$, there can be
extra stable regions in the middle. Of course, even in that case, it remains
true that all unstable regions exist in the critical strip.

We note that, in the numerical simulations shown in Secs. 5 and 6, the system
shows instability when the value of $q$ is not in the critical strip but close
to its boundaries, that is, $q\lesssim R_{G}$ and $q\gtrsim 1$ ($s\gtrsim 1$
and $s\lesssim 0$, equivalently). This is due to the fluctuation of the
unitary matrices. In fact, the above discussion of the stability is based on
the assumption that all the unitary matrices are set to the unit matrix.
Although the stability of the model with non-trivial configurations of the
unitary matrices remains an unresolved issue, it is reasonable to assume that
these matrices can easily transition to an unstable region due to the
tunneling effect, particularly when the coupling constant approaches to the
critical strip. This is thought to be the reason why the numerical simulations
around the critical strip fail.

###  4.2 Duality of the FKM Model

Since the partition function of the FKM model is described by the unitary
matrix weighted Ihara zeta function, the functional equation of the Ihara zeta
function is carried over to the FKM model as a duality between the small and
large couplings. In fact, the proof of the transformation (3.15) is also
applicable to the unitary weighted Ihara zeta function and the following
identity for the matrix weighted Ihara zeta function holds666 $Q$ has been
again redefined to
$Q=\operatorname{diag}(t_{1}\boldsymbol{1}_{N_{c}},t_{2}\boldsymbol{1}_{N_{c}},\cdots,t_{n_{V}}\boldsymbol{1}_{N_{c}})$.

$\displaystyle\zeta_{G}\left({R_{G}}/{q};U\right)^{-1}$
$\displaystyle=\frac{\left((q/R_{G})^{2}-1\right)^{(n_{E}-n_{V})N_{c}}\prod_{v\in
V}(\deg v-1)^{N_{c}}}{(q/R_{G})^{2n_{E}N_{c}}}$ $\displaystyle\hskip
56.9055pt\times\det\left({\bf
1}_{n_{V}}-q(R_{G}Q)^{-1}A_{U}+R_{G}^{-1}q^{2}(R_{G}Q)^{-1}\right)\,.$ (4.6)

Similar to the usual Ihara zeta function, the relation (4.6) reduces to the
functional equation of the matrix weighted Ihara zeta function

$\displaystyle\zeta_{G}\left(1/tq;U\right)^{-1}$
$\displaystyle=\frac{(1-t^{-2}q^{-2})^{n_{E}-\frac{n_{V}}{2}}(1-q^{-2})^{\frac{n_{V}}{2}}}{(1-q^{2})^{n_{E}-\frac{n_{V}}{2}}(1-t^{2}q^{2})^{\frac{n_{V}}{2}}}\zeta_{G}(q;U)^{-1}\,,$
(4.7)

when $G$ is $(t+1)$-regular. Since the effective action of the FKM model is
also written in terms of the unitary matrix weighted Ihara zeta function as
(2.17), there exists the duality on the couplings between $q$ and $1/tq$ for
the FKM model on the regular graph,

$S_{\rm eff}(1/tq;U)=S_{\rm eff}(q;U)\,,$ (4.8)

up to an irrelevant constant shift. Therefore, when $G$ is regular, the
dynamics of the FKM model in the region $q>1$ is identical with that in
$0<q<R_{G}$ even if the series expansion of the effective action (2.17)
converges only in $0<q<R_{G}$.

On the other hand, when $G$ is irregular, the right-hand side of the
transformation (4.7) cannot be interpreted as a unitary matrix weighted Ihara
zeta function. Therefore, the dynamics of the FKM model in $q>1$ is not the
same with that in the region $0<q<R_{G}$. Nevertheless, the effective action
of the FKM model in $q>1$ can be expressed in a sum of all possible Wilson
loops on $G$. In fact, changing the argument of (4.6) from $q$ to $R_{G}/q$,
we obtain

$\zeta_{G}(q;U)^{-1}=\frac{\left(q^{-2}-1\right)^{(n_{E}-n_{V})N_{c}}\prod_{v\in
V}(\deg v-1)^{N_{c}}}{q^{-2n_{E}N_{c}}}\\\ \times\det\left({\bf
1}_{n_{V}}-q^{-1}Q^{-1}A_{U}+q^{-2}Q^{-1}\right)\,.$ (4.9)

This means that the effective action (2.17) of the FKM model can also be
written as

$S_{\rm eff}(q;U)=-\gamma N_{c}\log\zeta_{G}(q;U)=\gamma
N_{c}\operatorname{Tr}\log\left({\bf
1}_{n_{V}}-q^{-1}Q^{-1}A_{U}+q^{-2}Q^{-1}\right)\,,$ (4.10)

up to irrelevant constant terms in general. In particular, since the series
expansion by $q^{-1}$ of the right-hand side becomes well-defined for large
$q$, we can interpret again the effective action of the FKM model in a sum of
all possible Wilson loops on $G$ as announced. Of course, the detail of the
expansion depends on $Q$ which gives non-trivial weights to the expansion when
$G$ is irregular. Therefore, to see the difference of the behavior of the FKM
model in the region $q>1$ from that in the region $0<q<R_{G}^{-1}$, we need to
carry out numerical simulations. We will see it in the next section in
concrete examples.

###  4.3 Phase structure of the FKM model on the cycle graph

Before considering the phase structure of the model on a general graph, it is
instructive to recall the exact results of the FKM model on the cycle graph
whose effective action is given by

$S_{\rm eff}^{C_{n}}(q;U)=-\gamma
N_{c}\operatorname{Tr}\left[\sum_{m=1}^{\infty}\frac{q^{nm}}{m}\left(U^{m}+U^{-m}\right)\right]\,.$
(4.11)

In [19], the partition function (the free energy) of the model in $0<q<1$ of
this system is analytically evaluated. In this subsection, we will explain the
complete phase structure of the cycle graph by combining this exact analysis
with the duality discussed in the previous subsection.

Figure 3: A cycle graph with three vertices.

The cycle graph $C_{n}$ is a 2-regular graph with $n$ vertices $\\{v_{i}\\}$
and $n$ edges; $e_{i}=\langle v_{i},v_{i+1}\rangle$ for $i=1,\cdots,n-1$ and
$e_{n}=\langle v_{n},v_{1}\rangle$. (See e.g. $C_{3}$ in Fig. 3). Since a
maximal tree of $C_{n}$ can be set to $T=e_{2}\cdots e_{n}$, we can fix the
gauge as $U_{e_{2}}=\cdots=U_{e_{n}}=\boldsymbol{1}_{N_{c}}$. Therefore, the
rank of $C_{n}$ is $1$ and only the unitary matrix on $e_{1}$, $U\equiv
U_{e_{1}}$, survives after the gauge fixing. Then, the partition function is
expressed in terms of the eigenvalues $e^{i\theta}$ ($i=1,\cdots,N_{c}$) of
the unitary matrices $U$ as

$\displaystyle Z_{C_{n}}$ $\displaystyle={\cal
N}\int_{-\pi}^{\pi}\prod_{i=1}^{N_{c}}d\theta_{i}\,e^{\sum_{j\neq
k}\log\left|\sin\frac{\theta_{j}-\theta_{k}}{2}\right|-N_{f}\sum_{i}\log\left(1-2q^{n}\cos\theta_{i}+q^{2n}\right)}\,,$
(4.12)

where the first term $\sum_{j\neq
k}\log\left|\sin\frac{\theta_{j}-\theta_{k}}{2}\right|$ arises from the
Jacobian associated with diagonalization of $U$ and ${\cal N}$ is an
irrelevant constant. We note that this expression is valid only in $0<q<1$
because the expansion of $U$ in the effective action is defined only in this
region.

(a) $\alpha=0.01$

(b) $\alpha=0.08$

Figure 4:  The eigenvalue distributions of the FKM model on $C_{n}$ with
$\gamma=10$, that is, ${q^{*}}^{n}=\frac{1}{19}\simeq 0.053$. (a) For
$q<{q^{*}}$, the eigenvalue distribution is continuous on the cycle. (b) For
${q^{*}}<q$, the eigenvalue distribution is concentrated around $\theta=0$.

In the limit of $N_{c}\to\infty$ with fixing $\gamma=N_{f}/N_{c}$, the saddle
point approximation becomes exact and the eigenvalues $\theta_{i}$ can be
regarded as a continuous variable $\theta(x)$ $(x\in[0,1])$. As a result, the
density of the eigenvalues,

$\rho(\theta)\equiv\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}\delta(\theta-\theta_{i})\,,$
(4.13)

can be evaluated exactly as

$\rho(\theta)=\begin{cases}\displaystyle\frac{1}{2\pi}\left(1+2\gamma\frac{q^{n}\cos\theta-q^{2n}}{1-2q^{n}\cos\theta+q^{2n}}\right)\,,&\quad(\theta_{0}=\pi)\vspace{4mm}\\\
\displaystyle\frac{2(\gamma-1)q^{n}}{\pi}\frac{\cos\frac{\theta}{2}}{1-2q^{n}\cos\theta+q^{2n}}\sqrt{\sin^{2}\frac{\theta_{0}}{2}-\sin^{2}\frac{\theta}{2}}\,,&\quad(\theta_{0}<\pi)\end{cases}$
(4.14)

where $\theta_{0}$ is determined by

$\sin^{2}\frac{\theta_{0}}{2}=\frac{(1-q^{n})^{2}}{4q^{n}}\frac{2\gamma-1}{(\gamma-1)^{2}}\,,$
(4.15)

which expresses the boundary of the support of the density function, namely,
$\rho(\theta)\geq 0$ in the region $-\theta_{0}\leq\theta\leq\theta_{0}$ and
$\rho(\theta)=0$ outside. In other words, $\theta_{0}=\pi$ means the
eigenvalue distribution is continuous on the cycle $e^{i\theta}$, whereas
$\theta_{0}<\pi$ means that the eigenvalues are distributed over only a
portion of the circle around $\theta=0$ (See Fig. 4). Since $\theta_{0}=\pi$
is equivalent to

$(q^{*})^{n}=\frac{1}{2\gamma-1}\,,$ (4.16)

which gives the critical coupling $q={q^{*}}$. This transition indexed by the
eigenvalue distribution is the so-called GWW phase transition [21, 22]. For
more detail, see [19].

The free energy of the system is defined by

$\displaystyle F_{G}$
$\displaystyle\equiv-\lim_{N_{c}\to\infty}\frac{1}{N_{c}^{2}}\log Z_{G}\,.$
(4.17)

In order to evaluate the order of the phase transition, it is also convenient
to define the “internal energy”,

$E_{G}\equiv\gamma\frac{\partial F_{G}}{\partial\gamma}\,,$ (4.18)

the “specific heat”,

$C_{G}\equiv-\gamma^{2}\frac{\partial^{2}F_{G}}{\partial\gamma^{2}}\,,$ (4.19)

and the “derivative of the specific heat”,

$dC_{G}\equiv\gamma^{3}\frac{\partial^{3}F_{G}}{\partial\gamma^{3}}\,,$ (4.20)

by regarding the parameter $\gamma$ as inverse temperature.

From the eigenvalue density (4.14), the free energy of the FKM model on
$C_{n}$ can be evaluated as [19]

$F_{C_{n}}=\begin{cases}\gamma^{2}\log(1-q^{2n})&(0<q\leq{q^{*}})\\\
(2\gamma-1)\log(1-q^{n})+\frac{1}{2}\log(2\gamma
q^{n})+f_{1}(\gamma)&({q^{*}}<q<1)\end{cases}\,,$ (4.21)

where

$f_{1}(\gamma)=(\gamma-1)^{2}\log\left(1-\frac{1}{\gamma}\right)-2\left(\gamma-\frac{1}{2}\right)^{2}\log\left(1-\frac{1}{2\gamma}\right)$
(4.22)

is determined by the condition that the free energy is continuously connected
at the critical coupling $q={q^{*}}$. Substituting the expression (4.21) into
the definitions (4.18), (4.19) and (4.20), we obtain the internal energy,

$E_{C_{n}}=\begin{cases}\displaystyle
2\gamma^{2}\log(1-q^{2n})&(0<q\leq{q^{*}})\\\ \displaystyle
2\gamma\log(1-q^{n})+e_{1}(\gamma)&({q^{*}}<q<1)\end{cases}\,,$ (4.23)

where

$e_{1}(\gamma)\equiv
2\gamma\left(\gamma-1\right)\log\left(1-\frac{1}{\gamma}\right)-2\gamma(2\gamma-1)\log\left(1-\frac{1}{2\gamma}\right)\,,$
(4.24)

the specific heat,

$C_{C_{n}}=\begin{cases}\displaystyle-2\gamma^{2}\log(1-q^{2n})&(0<q\leq{q^{*}})\\\
\displaystyle-2\gamma^{2}\log\frac{4\gamma(\gamma-1)}{(2\gamma-1)^{2}}&({q^{*}}<q<1)\end{cases}\,,$
(4.25)

and the derivative of the specific heat,

$dC_{C_{n}}=\begin{cases}0&(0<q\leq{q^{*}})\\\
\displaystyle\frac{2\gamma^{3}}{2\gamma^{3}-3\gamma^{2}+\gamma}&({q^{*}}<q<1)\end{cases}\,,$
(4.26)

for $G=C_{n}$. It is easy to check that the internal energy (4.23) and the
specific heat (4.25) connect continuously at $q={q^{*}}$ while the derivative
of the specific heat does not. This means that the GWW phase transition of the
FKM model on $C_{n}$ is third-order as the same as the original GWW model777
In [19], it has been shown that the phase transition is third-order in the
sense that the first and second derivative of the free energy in terms of the
coupling constant $\mathit{q}$. It is interesting that the transition is still
third-order in terms of $\gamma$ even if $q$ and $\gamma$ are independent in
the FKM model. .

Although we have considered only the region $0<q<1$ so far, we can extend it
to $q>1$ by using the small/large coupling duality explained in Sec. 4.2.
Since the cycle graph is 2-regular graph, it obeys the functional equation,

$\zeta_{Cn}(1/q;U)=q^{-2nN_{c}}\zeta_{C_{n}}(q;U)\,.$ (4.27)

This means that the free energy in the region $q>1$ is given by

$F_{C_{n}}(1/q)=F_{C_{n}}(q)-2n\gamma\log q\,.$ (4.28)

In order to display the duality in a symmetric manner, it is convenient to
change the variable from $q$ to $s$ as (3.14). However, we cannot adopt
$q=R_{G}^{s}$ in this case because the cycle graph has exceptionally
$R_{G}=1$, in other words, the critical strip shrinks to the line associated
with $|q|=1$. We thus instead define

$s=-\log q\,,$ (4.29)

for $G=C_{n}$. In this parametrization, the original region $0<q<1$ and the
dual region $q>1$ correspond to $s>0$ and $s<0$, respectively. Corresponding
to (4.16), the critical coupling in terms of the parameter $s$ is

$s^{*}=\frac{1}{n}\log(2\gamma-1)\,.$ (4.30)

Therefore, we expect that the FKM model has two GWW phase transition points at
$s=\pm s^{*}$.

Taking into account the duality (4.28), the free energy for all region of $s$
becomes

$F_{C_{n}}=\begin{cases}\gamma^{2}\log(1-e^{-2ns})&(s^{*}\leq s)\\\
(2\gamma-1)\log(1-e^{-ns})+\frac{1}{2}\log(2\gamma
e^{-ns})+f_{1}(\gamma)&(0<s<s^{*})\\\
(2\gamma-1)\log(1-e^{ns})+\frac{1}{2}\log(2\gamma
e^{ns})+f_{1}(\gamma)-2n\gamma s&(-s^{*}<s<0)\\\
\gamma^{2}\log(1-e^{2ns})-2n\gamma s&(s\leq-s^{*})\end{cases}\,.$ (4.31)

Similarly, the internal energy, the specific heat and the derivative of the
specific heat become

$\displaystyle E_{C_{n}}$ $\displaystyle=\begin{cases}\displaystyle
2\gamma^{2}\log(1-e^{-2ns})&(s^{*}\leq s)\\\ \displaystyle
2\gamma\log(1-e^{-ns})+e(\gamma)&(0<s<s^{*})\\\ \displaystyle
2\gamma\log(1-e^{ns})+e(\gamma)-2n\gamma s&(-s^{*}<s<0)\\\ \displaystyle
2\gamma^{2}\log(1-e^{2ns})-2n\gamma s&(s\leq-s^{*})\end{cases}\,,$ (4.32)
$\displaystyle C_{C_{n}}$
$\displaystyle=\begin{cases}\displaystyle-2\gamma^{2}\log(1-e^{-2n|s|})&(s^{*}\leq|s|)\\\
\displaystyle-2\gamma^{2}\log\frac{4\gamma(\gamma-1)}{(2\gamma-1)^{2}}&(0<|s|<s^{*})\end{cases}\,,$
(4.33) $\displaystyle dC_{C_{n}}$
$\displaystyle=\begin{cases}0&(s^{*}\leq|s|)\\\
\displaystyle\frac{2\gamma^{3}}{2\gamma^{3}-3\gamma^{2}+\gamma}&(0<|s|<s^{*})\end{cases}\,,$
(4.34)

respectively. We see that the order of the phase transition at $s=-s^{*}$ is
also of the third-order as expected. We plot (4.32), (4.33) and (4.34) in Fig.
5 where we can see the order of the phase transitions and the duality clearly.

(a) $E_{C_{n}}$

(b) $C_{C_{n}}$

(c) $dC_{C_{n}}$

Figure 5:  The behavior of the observables of the FKM model on $C_{n}$ with
$n=3$, $\gamma=10$ and $N_{c}=16$. The vertical dashed lines show the
positions of the critical coupling $\pm s^{*}\simeq\pm 0.98$ in this
parametrization. (a) the internal energy, (b) the specific heat, (c) the
derivative of the specific heat.

###  4.4 GWW phase transitions in the FKM model for large $\gamma$

As shown in [19], the FKM model for $N_{c}\to\infty$ generally undergoes the
GWW phase transitions in the region $0<q<R_{G}$ when $\gamma$ exceeds a value
of ${\cal O}(1)$. Since the GWW phase transition occurs when the distribution
of eigenvalues of a Wilson loop (or a Polyakov loop) on $S^{1}$ breaks or
connects in the process of changing the coupling constant, the number of the
phase transitions in the FKM model on a general graph $G$ is essentially the
same as the rank of $G$, that is, the number of fundamental cycles of $G$
which generate the Wilson loops. However, as we will see shortly, the actual
number of the phase transitions depends on the structure of the graph since
some of the fundamental cycles simultaneously exhibit the GWW phase
transition.

In order to analyze the phase structure of the FKM model on a general graph,
it is instructive to look at the model on the cycle graph considered in the
previous section in the limit,

$\gamma\to\infty,\quad q\to 0,\quad\gamma q^{n}\equiv\lambda^{-1}\,:\,{\rm
fixed}\,.$ (4.35)

In this limit, the effective action (4.11) reduces to

$S_{\rm
eff}^{C_{n}}(q;U)\to-\frac{N_{c}}{\lambda}\operatorname{Tr}\left(U+U^{\dagger}\right)\,,$
(4.36)

which is nothing but the action of the GWW model. Not only the tree level
action, but also the limit of the exact solution of the eigenvalue density
(4.14) becomes that of the GWW model [19]. Since the phase transition of the
GWW model takes place at $\lambda=2$, the position of the phase transition of
the FKM model for large $\gamma$ is estimated as

$q^{*}\simeq\left(2\gamma\right)^{-1/n}\,,$ (4.37)

which is consistent with the exact result (4.16). In addition, the phase
transition of the GWW model is third-order, which is the also consistent with
the fact that the phase transition of the FKM model on $G=C_{n}$ is third-
order. From this observation, we see that we can estimate the nature of the
phase transition of the FKM model for large $\gamma$ when the effective action
can be approximated by the GWW model.

Whether the FKM model can be approximated by the GWW model or not depends on
the structure of the graph. To see it, we adopt as the set of the fundamental
cycles the one with the smallest length888 A choice of fundamental cycles is
not unique in general. ,

${\cal F}\equiv\\{[C_{a}]\,|\,a=1,\cdots,r,\
|C_{1}|\leq\cdots\leq|C_{r}|\\}\,,$ (4.38)

and suppose that ${\cal F}$ consists of cycles of lengths $l_{i}$
($i=1,\cdots,b_{G}$, $l_{1}<\cdots<l_{b_{G}}$). In other words, $b_{G}$ is the
number of the types of polygons in the fundamental cycles of the least
lengths. Correspondingly, we define $b_{G}$ subsets of ${\cal F}$,

${\cal F}_{i}\equiv\\{[C_{a}]\in{\cal F}\,|\,|C_{a}|=l_{i}\\},\quad
m_{i}\equiv|{\cal F}_{i}|\,,$ (4.39)

which is the set of $l_{i}$-gons in $\cal{F}$. In addition, we define a subset
of $[{\cal P}_{R}]$ whose length is $l_{i}$ but are not in ${\cal F}_{i}$ as

${\cal G}_{i}\equiv\left\\{[C]\in[{\cal P}_{R}]\,|\,|C|=l_{i},\
[C]\not\in{\cal F}_{i}\right\\}\,,$ (4.40)

whose elements are composite cycles by definition. In this setting, the
effective action of the FKM model in the region $0<q<R_{G}$ is expanded by $q$
as

$S_{\rm eff}(q;U)\\\ =-\gamma
N_{c}\operatorname{Tr}\left[\sum_{i=1}^{b_{G}}q^{l_{i}}\left(\sum_{C_{a}\in{\cal
F}_{i}}\left(U_{C_{a}}+U_{C_{a}}^{-1}\right)+\sum_{C^{\prime}\in{\cal
G}_{i}}\left(U_{C^{\prime}}+U_{C^{\prime}}^{-1}\right)\right)+{\cal
O}(q^{l_{b_{G}}+1})\right]\,.$ (4.41)

When ${\cal G}_{i}=\emptyset$ for all $i=1,\cdots,b_{G}$, the situation is the
easiest: For large $\gamma$ with $\gamma q^{l_{i}}={\cal O}(1)$, Wilson loops
of the length more than $l_{i}$ are decoupled from the system, and the
effective action in the limit of $\gamma\to\infty$ can be regarded as a
summation of independent GWW models. Therefore, we can expect that the GWW
phase transitions occur at

$q^{*}_{i}\simeq(2\gamma)^{1/l_{i}}\,\quad\text{for }i=1,\cdots,b_{G}.$ (4.42)

On the other hand, if ${\cal G}_{i}$ for some $i$ is not empty, the effective
action in the same limit is not simply a summation of the GWW models but also
includes non-trivial coupling terms of the fundamental Wilson loops. Because
of the existence of these coupling terms, we cannot use the result of the GWW
model. In this paper, we examine one of such models numerically in the next
section. A theoretical analysis of this kind is discussed separately.

As discussed in the previous section, the FKM model has a dual expression in
the region $q>1$ where the effective action can be expanded by $1/q$.
Therefore, as long as the phase transitions occur in the region $0<q<R_{G}$,
they are expected to occur also in the region $q>1$.

If the FKM model is defined on a regular graph, the region $q>1$ is just a
copy of the region $0<q<R_{G}$ because of the exact duality (4.8). However,
for an irregular graph, the phase structure in $q>1$ is different from that in
$0<q<R_{G}$ in general. In this case, we have to write down the
$1/q$-expansion of the effective action (4.10), which generally takes the
form,

$S_{\rm eff}(q;U)\\\ =-\gamma
N_{c}\operatorname{Tr}\left[\sum_{i=1}^{b_{G}}g_{i}(q)\left(\sum_{C_{a}\in{\cal
F}_{i}}\left(U_{C_{a}}+U_{C_{a}}^{-1}\right)+\sum_{C^{\prime}\in{\cal
G}_{i}}\left(U_{C^{\prime}}+U_{C^{\prime}}^{-1}\right)\right)+{\cal
O}(q^{-l_{b_{G}}-1})\right]\,,$ (4.43)

where $g_{i}(q)$ is a function of $q$ which is defined by $1/q$-expansion.
When all ${\cal G}_{i}$ are empty, we can expect that the phase transitions
occur at $\tilde{q}_{i}^{*}>1$ with

$\gamma g_{i}(\tilde{q}_{i}^{*})=\frac{1}{2}\,,$ (4.44)

as well as in the region $0<q<R_{G}$. On the other hand, if ${\cal G}_{i}$ is
not empty for some $i$, the situation is expected to be different from the GWW
model. In this case, we have to carry out the numerical computations to see
the phase structure.

In summary, the phase structure of the FKM model in large $\gamma$ depends on
whether ${\cal G}_{i}$ are empty or not and whether the graph $G$ is regular
or not. The cycle graph is a simplest example of the case ${\cal G}_{i}$ are
empty and the graph is regular. We will give concretely an example of a
regular graph with non-trivial ${\cal G}_{i}$ (the tetrahedron) and examples
of irregular graphs with empty ${\cal G}_{i}$ (the double triangle and the
triangle-square) in turn. We also consider the case that the number of the
phase transitions in $0<q<1$ and $q>1$ is different (the triple triangle).

Before analyzing concrete examples, we emphasize that, even if we can analyze
the FKM model by using the GWW model, the FKM model differs from the GWW model
because of the existence of the duality. For example, we cannot discuss the
phase transition at $q^{n}=2\gamma-1$, that is, the phase transition in the
dual region from the GWW model. This is because the theory is defined as a sum
of infinitely many Wilson loops.

###  4.5 Concrete examples

#### Tetrahedron

The tetrahedron ($K_{4}$) is an example of regular graphs as explained in Sec.
3.1 and depicted in Fig. 1 (a). The maximal convergence radius is determined
by the (common) degree of the vertices like $R_{K_{4}}=1/t=1/2$. In
particular, $q=R_{K_{4}}=1/2$ is the only nontrivial real pole of the Ihara
zeta function (3.18). The set of six directed edges is defined by

$E=\left\\{e_{1},\cdots,e_{6}\right\\}\equiv\left\\{\langle
v_{1},v_{2}\rangle,\langle v_{1},v_{3}\rangle,\langle
v_{1},v_{4}\rangle,\langle v_{2},v_{3}\rangle,\langle
v_{3},v_{4}\rangle,\langle v_{4},v_{2}\rangle\right\\}\,,$ (4.45)

where $V=\\{v_{1},v_{2},v_{3},v_{4}\\}$ stands for the set of the vertices
depicted in Fig. 1 (a).

The rank of $K_{4}$ is $r=3$ and the set ${\cal F}$ is not unique since it has
four triangles, but only three of them can be fundamental. We here choose

${\cal F}={\cal
F}_{1}=\\{[e_{1}e_{4}e_{2}^{-1}],[e_{2}e_{5}e_{3}^{-1}],[e_{3}e_{6}e_{1}^{-1}]\\}\equiv\\{[C_{1}],[C_{2}],[C_{3}]\\}.$
(4.46)

Then, the remaining triangle $e_{4}e_{5}e_{6}$ is regarded as a composite
cycle,

$e_{4}e_{5}e_{6}=C_{1}C_{2}C_{3}\,.$ (4.47)

Therefore, we see $b_{K_{4}}=1$ and

$l_{1}=3,\quad m_{1}=3,\quad{\cal G}_{1}=\\{[C_{1}C_{2}C_{3}]\\}\,.$ (4.48)

The $q$-expansion of the effective action is given by

$S_{\rm eff}(q;U)=-\gamma
N_{c}\operatorname{Tr}\Big{[}q^{3}\left(U_{1}+U_{2}+U_{3}+U_{1}U_{2}U_{3}+{\rm
h.c.}\right)+{\cal O}(q^{-4})\Big{]}\,,$ (4.49)

where the interaction term $U_{1}U_{2}U_{3}$ exists because ${\cal G}_{1}$ is
not empty, which prevents the theory from being associated with the GWW model.
We will investigate the nature of the phase transition numerically in the next
section. Since $K_{4}$ is regular, the FKM model on $K_{4}$ has the exact
duality between the regions $0<q<R_{K_{4}}$ and $q>1$. We will also check it
numerically in the next section.

#### Double triangle

As explained in Sec. 3.2 and depicted in Fig. 1 (b), the double triangle (DT)
is an example of irregular graphs and the maximal convergence radius is given
by $R_{\rm DT}=0.65729\cdots$. We see that $q=R_{\rm DT}$ is the only non-
trivial real pole of $\zeta_{\rm DT}(q)$. The set of five directed edges is
defined by

$E=\left\\{e_{1},\cdots,e_{5}\right\\}\equiv\left\\{\langle
v_{1},v_{2}\rangle,\langle v_{2},v_{3}\rangle,\langle
v_{3},v_{4}\rangle,\langle v_{4},v_{1}\rangle,\langle
v_{3},v_{1}\rangle\right\\}\,,$ (4.50)

where $V=\\{v_{1},v_{2},v_{3},v_{4}\\}$ stands for the set of the vertices
depicted in Fig. 1 (b).

The rank of DT is $r=2$ and the set of generators of the minimal lengths is

${\cal F}={\cal F}_{1}=\\{[e_{1}e_{2}e_{5}],[e_{5}^{-1}e_{3}e_{4}]\\}\,,$
(4.51)

which consists of two triangles. Since there is no other triangle in the
cycles of DT, we see $b_{\rm DT}=1$ and

$l_{1}=3,\quad m_{1}=2,\quad{\cal G}_{1}=\emptyset\,.$ (4.52)

Correspondingly, the $q$-expansion of the effective action in $0<q<R_{\rm DT}$
is given by

$S_{\rm eff}(q;U)=-\gamma
N_{c}\operatorname{Tr}\left[\sum_{a=1,2}q^{3}\left(U_{a}+U_{a}^{-1}\right)+{\cal
O}(q^{-4})\right]\,.$ (4.53)

Therefore, we expect that the FKM model on DT exhibits the GWW phase
transition only once at

$q^{*}\simeq(2\gamma)^{-\frac{1}{3}}\,,$ (4.54)

in the region $0<q<R_{G}$ for large $\gamma$.

On the other hand, for $q>1$, the $1/q$-expansion of the dual effective action
(4.10) on DT is given by

$S_{\rm eff}(q;U)=-\gamma
N_{c}\operatorname{Tr}\left[g(q)\sum_{a=1,2}\left(U_{a}+U_{a}^{-1}\right)+{\cal
O}(q^{-4})\right],$ (4.55)

up to constant terms, where

$g(q)=\frac{1}{4q^{3}}-\frac{3}{16q^{5}}+\frac{17}{64q^{7}}-\frac{63}{256q^{9}}+{\cal
O}(q^{-11})\,.$ (4.56)

As expected, it is not completely symmetric to the region $0<q<R_{\rm DT}$,
although there is still a symmetry between $U_{1}$ and $U_{2}$. The position
of the phase transition in the region $q>1$ is expected to be the solution of

$\gamma g(\tilde{q}^{*})=\frac{1}{2}\,,$ (4.57)

for large $\gamma$.

#### Triangle-square

Figure 6: A triangle-square (TS)

The triangle-square (TS) is an asymmetric irregular graph which contains two
cycles the different length 3 and 4 (See Fig. 6). The Ihara zeta function of
TS is given by

$\zeta_{\rm
TS}(q)=\frac{1}{(1-q^{2})(1+q+q^{2}-q^{4}-2q^{5})(1-q+q^{2}-2q^{3}+q^{4}-2q^{5})}.$
(4.58)

The maximal convergence radius $R_{\rm TS}$ of $\zeta_{\rm TS}(q)$ is read
from this zeta function as $R_{\rm TS}=0.70222\cdots$ and $q=R_{\rm TS}$ is
only the non-trivial real pole of $\zeta_{\rm TS}(q)$. The set of six directed
edges is defined by

$E=\left\\{e_{1},\cdots,e_{6}\right\\}\equiv\left\\{\langle
v_{1},v_{2}\rangle,\langle v_{2},v_{3}\rangle,\langle
v_{3},v_{4}\rangle,\langle v_{4},v_{5}\rangle,\langle
v_{5},v_{1}\rangle,\langle v_{3},v_{1}\rangle\right\\}\,,$ (4.59)

where $V=\\{v_{1},v_{2},v_{3},v_{4},v_{5}\\}$ stands for the set of the
vertices depicted in Fig. 6.

The rank of TS is again $r=2$ and the set of generators of the minimal lengths
is

${\cal F}=\\{[e_{1}e_{2}e_{6}],[e_{6}^{-1}e_{3}e_{4}e_{5}]\\}\,,$ (4.60)

which consists of a triangle and a square. In this case, ${\cal F}$ is the
disjoint union of

${\cal F}_{1}=\\{[e_{1}e_{2}e_{6}]\\}\quad{\rm and}\quad{\cal
F}_{2}=\\{[e_{6}^{-1}e_{3}e_{4}e_{5}]\\}\,.$ (4.61)

Since there is no other triangle and square in the cycles of TS, we see
$b_{\rm TS}=2$ and

$(l_{1},l_{2})=(3,4),\quad(m_{1},m_{2})=(1,1),\quad({\cal G}_{1},{\cal
G}_{2})=(\emptyset,\emptyset)\,.$ (4.62)

Correspondingly, the $q$-expansion of the effective action in $0<q<R_{\rm TS}$
is given by

$S_{\rm eff}(q;U)=-\gamma
N_{c}\operatorname{Tr}\Big{[}q^{3}\left(U_{1}+U_{1}^{-1}\right)+q^{4}\left(U_{2}+U_{2}^{-1}\right)+{\cal
O}(q^{-5})\Big{]}\,.$ (4.63)

Therefore, we expect that the FKM model on TS exhibits the GWW phase
transition at

$q_{1}^{*}\simeq(2\gamma)^{-\frac{1}{3}}\quad{\rm and}\quad
q_{2}^{*}\simeq(2\gamma)^{-\frac{1}{4}}\,,$ (4.64)

in the region $0<q<R_{G}$ for large $\gamma$.

On the other hand, the $1/q$-expansion of the dual effective action (4.10) in
$q>1$ is given by

$S_{\rm eff}(q;U)=-\gamma
N_{c}\operatorname{Tr}\left[\sum_{a=1,2}g_{a}(q)\left(U_{a}+U_{a}^{-1}\right)+{\cal
O}(q^{-5})\right]$ (4.65)

up to constant terms, where

$\begin{split}g_{1}(q)&=\frac{1}{4q^{3}}-\frac{3}{16q^{5}}+\frac{17}{64q^{7}}-\frac{39}{256q^{9}}+{\cal
O}(q^{-11}),\\\
g_{2}(q)&=\frac{1}{4q^{4}}-\frac{3}{16q^{6}}+\frac{17}{64q^{8}}-\frac{63}{256q^{10}}+{\cal
O}(q^{-12}).\end{split}$ (4.66)

Again, it is not completely symmetric to the region $0<q<R_{\rm TS}$ although
the structure is similar. From the effective action, the FKM model on TS is
expected to undergo phase transitions twice in the dual region $q>1$ as well
as in $0<q<R_{\rm TS}$, and the positions of the phase transitions
$\tilde{q}_{i}^{*}$ in $q>1$ are expected to be the solutions of

$\gamma g_{i}(\tilde{q}_{i}^{*})=\frac{1}{2}\,,$ (4.67)

for large $\gamma$.

#### Triple triangle

Figure 7: A triple triangle (TT)

The triple triangle (TT) is an irregular graph which contains three cycles
(See Fig. 7). The Ihara zeta function of TT is given by

$\zeta_{\rm
TT}(q)=\frac{1}{\left(2q^{2}+q+1\right)\left(3q^{3}+q-1\right)\left(2q^{4}+q^{3}+2q^{2}+q+1\right)\left(q+1\right)^{2}\left(q-1\right)^{3}}$
(4.68)

The maximal convergence radius $R_{\rm TT}$ of $\zeta_{\rm TT}(q)$ is read
from this zeta function as $R_{\rm TT}=0.53656\cdots$ and $q=R_{\rm TT}$ is
the only non-trivial real pole of $\zeta_{\rm TT}(q)$. The set of seven
directed edges is defined as

$E=\left\\{e_{1},\cdots,e_{7}\right\\}\equiv\left\\{\langle
v_{1},v_{2}\rangle,\langle v_{2},v_{3}\rangle,\langle
v_{3},v_{4}\rangle,\langle v_{4},v_{5}\rangle,\langle
v_{5},v_{1}\rangle,\langle v_{3},v_{1}\rangle,\langle
v_{3},v_{5}\rangle\right\\}\,,$ (4.69)

where $V=\\{v_{1},v_{2},v_{3},v_{4},v_{5}\\}$ expresses the set of the
vertices depicted in Fig. 7.

The rank of TT is again $r=3$ and the set of generators of the minimal lengths
is

${\cal F}={\cal
F}_{1}=\\{[e_{1}e_{2}e_{6}],[e_{3}e_{4}e_{7}^{-1}],[e_{5}e_{6}^{-1}e_{7}]\\}\,.$
(4.70)

Since there is no other triangle in the cycles of TT, we see $b_{\rm TT}=1$
and

$l_{1}=3,\quad m_{1}=3,\quad{\cal G}_{1}=\emptyset\,.$ (4.71)

Correspondingly, the $q$-expansion of the effective action in $0<q<R_{\rm TT}$
is given by

$S_{\rm eff}(q;U)=-\gamma
N_{c}\operatorname{Tr}\left[\sum_{a=1,2,3}q^{3}\left(U_{a}+U_{a}^{-1}\right)+{\cal
O}(q^{-4})\right]\,.$ (4.72)

Therefore, we expect that the FKM model on TS exhibits the GWW phase
transition only once at

$q^{*}\simeq(2\gamma)^{-\frac{1}{3}}\,,$ (4.73)

in the region $0<q<R_{G}$ for large $\gamma$.

Interestingly, for $q>1$, the $1/q$-expansion of the dual effective action
(4.10) on TT splits as follows

$\displaystyle S_{\rm eff}(q;U)=-\gamma
N_{c}\operatorname{Tr}\Bigl{[}g_{1}(q)\left(U_{1}+U_{2}+U_{1}^{-1}+U_{2}^{-1}\right)+g_{2}(q)\left(U_{3}+U_{3}^{-1}\right)+{\cal
O}(q^{-4})\Bigr{]}$ (4.74)

up to constant terms, where

$\begin{split}g_{1}(q)&=\frac{1}{6q^{3}}-\frac{11}{72q^{5}}+\frac{233}{864q^{7}}+{\cal
O}(q^{-9})\,,\\\
g_{2}(q)&=\frac{1}{12q^{3}}-\frac{17}{144q^{5}}+\frac{115}{576q^{7}}+{\cal
O}(q^{-9}).\end{split}$ (4.75)

Therefore, as opposed to the cases of DT and TS, the FKM model in $q>1$ on TT
is not only asymmetric with the region $0<q<R_{\rm TT}$, but also the phase
structure in large $\gamma$ differs from $q<1$. It is expected that the model
undergo phase transitions twice in the dual region $q>1$ unlike in $0<q<R_{\rm
TT}$, and the positions of the phase transitions $\tilde{q}_{i}^{*}$ ($i=1,2$)
in $q>1$ are expected to be the solutions of

$\gamma g_{i}(\tilde{q}_{i}^{*})=\frac{1}{2}\,,$ (4.76)

for large $\gamma$.

##  5 Numerical Results I: Regular Graphs

In this section, we numerically examine the phase structure of the FKM model
on the regular graphs to check the analytic investigation described so far.

###  5.1 Algorithm and parameter settings

The simulations are carried out by using the HMC algorithm. The action used in
the code is defined through the path representation of the unitary matrix
weighted Ihara zeta function (2.20) to optimize the calculation of the
Hamiltonian and the force.

The observable which we mainly measure is the specific heat, which shows a
cusp behavior at the phase transition point if the transition is of the third-
order as the standard of the GWW phase transition. We will plot the behavior
of the specific heat for each graph described in Sec. 4.5 toward the parameter
$s$ defined in (3.14) or (4.29) for $C_{n}$ in order to see the dual
expression discussed in Sec. 4.2 clearly. If the graph is regular, the plot of
any observable must be fold-symmetric around $s=1/2$, while the symmetry must
be broken if the graph is irregular.

The regions of the coupling constant $s$ to carry out the simulations are
$s>1$ and $s<0$ since the FKM model is unstable in the critical strip $0<s<1$
as discussed in Sec. 4.1. Although the simulations fail when $s$ is close to
the boundary of the critical strip, $s=0$ and $s=1$, it is natural as
discussed at the end of Sec. 4.1 and it does not affect to our purpose to see
the phase structure of the model.

Since the GWW phase transition takes place in large $N_{c}$, we have to use
sufficiently large $N_{c}$ to see the nature of the transition. On the other
hand, the cost of the computation becomes higher when we adopt larger value of
$N_{c}$. Therefore, we have to find an appropriate value of $N_{c}$. For this
purpose, we use the model on the cycle graph where the exact result is known
as shown in Sec. 4.3. As we will see shortly in Sec. 5.2, $N_{c}=16$ is
sufficient to consider $N_{c}$ to be large.

The parameters of the molecular dynamics (the number of divisions $N_{\tau}$
and the time step $D_{\tau}$) are tuned so that the acceptance ratio exceeds
$0.75$ with fixing the combination $N_{\tau}D_{\tau}$ to be $0.1$. We start
the simulation with $U_{a}=\boldsymbol{1}_{N_{c}}$ $(a=1,\cdots,r)$ and
discard the first $50000$ trajectories for thermalization. After that, we
compute the observables every ten trajectories and have gathered $20000$
samples for each set of parameters. The average and standard error of each
physical quantity are evaluated using the jackknife method with a bin size for
which the autocorrelation becomes sufficiently small. Correspondingly, the
number of independent samples is $1000$–$10000$ for each parameter.

###  5.2 Cycle graph

Figure 8:  The plots of the specific heat of the FKM model on $C_{3}$ with
varying $N_{c}=4,8,16,32$ and fixed $\gamma=128$. The solid line indicates the
exact result obtained in the large $N_{c}$ limit. The vertical dashed line
shows the phase transition point which is also evaluated at large $N_{c}$.

The cycle graph $C_{n}$ is a 2-regular graph with $n$ vertices. (See Fig. 3.)
In this case, we can choose only one representative unitary matrix $U$ after
fixing the gauge. First of all, we simulate and mesure the specific heat at
$\gamma=128$ while varying $N_{c}$. The plot of the specific heat near the
phase transition point is shown in Fig. 8. The plot indicates that the
behavior at $N_{c}=16$ and $N_{c}=32$ is almost the same. This means that the
analytical results in the large $N_{c}$ limit are well reproduced even at
$N_{c}=16$.999 On the cycle graph, the autocorrelation at smaller $N_{c}=4$ is
stronger than at $N_{c}=16$, and the errors of the measurements become large.
In this sense as well, $N_{c}=16$ is the best choice in the following
measurements.

Figure 9: The numerical results of the specific heat on the $C_{3}$ graph. We
have changed $\gamma$ as $\gamma=16,128,1024$ with fixing $N_{c}=16$. The
solid line and the vertical dashed line show the analytic result of the
specific heat and the phase transition points in the large $N_{c}$ limit,
respectively.

We then measure the specific heat with fixing $N_{c}=16$ while varying
$\gamma=16,128$ and $1024$, which is plotted in Fig. 9 on $s$-axis. Note that
we have removed some statistical outliers of the measurements near $s=0$,
namely the shrunk critical strip, where the numerical simulation becomes
unstable. The exact solution in the large $N_{c}$ limit is again drawn by the
solid line, which is valid at any value of finite $\gamma$. It shows that the
numerical results agree with the analytic results including the dependence of
$\gamma$. Similarly, it can be seen that the phase transition points of the
exact solutions, represented by the dashed lines in the same figure, also
perfectly match the results of the simulation. We can also see the exact
duality between the small ($s>1$) and large ($s<0$) coupling regions as
discussed in Sec. 4.3.

(a) $s=2.15848$

(b) $s=2.54084\simeq s^{*}$

(c) $s=2.84215$

Figure 10: The eigenvalue distributions of $U$ for $N_{c}=16$ and
$\gamma=1024$. (a) and (c) corresponds to the deconfinement and confinement
phase, respectively. (b) is the eigenvalue distribution near the phase
transition point.

Fig. 10 shows the eigenvalue distributions of $U$ at three typical parameters
of $s$. The eigenvalue distribution at $s\simeq s^{*}$ is depicted in Fig. 10
(b), and those of smaller and larger values of $s^{*}$ are depicted in (a) and
(c), respectively. This result indicates that this phase transition is the GWW
phase transition. In other words, it separates the deconfinement phase (a) and
the confinement phase (c).

The results of the above numerical calculations perfectly reproduce the
analytical results in the large $N_{c}$ limit, which means that our simulation
code can be trusted.

###  5.3 Tetrahedron

The tetrahedron, considering as a graph, is the complete graph $K_{4}$
depicted in Fig. 1 (a). Since $K_{4}$ is a 3-regular graph, the FKM model on
$K_{4}$ should have the duality between $q$ and $1/(2q)$ as discussed in Sec.
4.5. Thus, if we define $s$ by $q=2^{-s}$ as (3.9), the specific heat must
indicate the duality between $s$ and $1-s$ across the critical strip.

As discussed in Sec. 4.5, we expect that the phase transition occurs at the
same time for the three independent unitary matrices $U_{1}$, $U_{2}$, $U_{3}$
because of the symmetry of each cycle (face). Thus, the FKM model on $K_{4}$
expresses the phase transition only once, both in the small and large $q$
regions.

Figure 11: The specific heat of the FKM model on the tetrahedron $K_{4}$ with
varying $\gamma=128,1024,16384$ and fixing $N_{c}=16$. The solid line, drawn
to compare the results of the FKM and GWW models, represents three times the
specific heat of the GWW model for $\gamma=16384$.

The numerical result of the specific heat on $K_{4}$ is plotted in Fig. 11.
Exact duality is observed between the small ($s>1$) and large ($s<0$) coupling
regions, with one phase transition point in each region being perfectly
symmetrical. We also see that the behavior of the specific heat seems to be
quite different from that of the naive GWW approximation as expected from the
discussion in Sec. 4.5. Indeed, comparing the results of the simulation with
the solid line in Fig. 11), which is the specific heat of the GWW model for
the three independent cycles $U_{1}$, $U_{2}$, $U_{3}$, we see a large
discrepancies in the shapes. Concretely speaking, the transition points are
different and “pointed elbow” behaviors can be seen at the phase transition
points in the numerical simulation. This is because the interaction terms
$\operatorname{Tr}U_{1}U_{2}U_{3}+\rm{h.c.}$ in (4.49), which come from the
fourth Wilson loop on $K_{4}$, cannot be negligible even in the large $\gamma$
limit. We emphasize that the solid line in Fig. 11 is just a sum of the three
cycles without the interaction.

(a) $s=5.0$

(b) $s=5.11979$

Figure 12:  The eigenvalue distribution of $U_{1}$, $U_{2}$, $U_{3}$ and
$U_{1}U_{2}U_{3}$ in the FKM model on $K_{4}$ at $s=5.0$, the phase transition
point of the GWW model (a), and at $s=5.11979$, the value of $s$ closest to
the peak of the pointed elbow of the specific heat in Fig. 11 (b).

To check if the phase transition point of the FKM model is different from that
of the GWW model from another point of view, we compare the eigenvalue
distributions at different values of $s$ in Fig. 12. Ignoring the interaction
term allows the application of the GWW model’s results, predicting the phase
transition point at $s=5.0$. However, the eigenvalue distribution at this
value, shown in Fig. 12 (a), does not exhibit the characteristics of the GWW
phase transition. On the other hand,

Figure 13: The specific heat of the FKM model on $K_{4}$ near the phase
transition point with varying $N_{c}=4,8,16,32$ and fixing $\gamma=16384$.

$s=5.11979$ is the value in our data, which is closest to the peak of the
pointed elbow in Fig. 11, and the eigenvalue distribution at this value of
$s$, shown in Fig. 12 (b), displays the typical shape of the GWW phase
transition. This result confirms that the GWW phase transition occurs at the
pointed elbow of Fig. 11, which is shifted from the GWW model by the effect of
the interection term.

We next investigate in detail the shape of the “elbow” near the phase
transition. Fig. 13 shows an enlarged view of the behavior of the specific
heat near the phase transition at large $\gamma=16384$. In order to see the
change of the width of the elbow in taking the large $N_{c}$ limit, $N_{c}=32$
is included. As shown in the figure, the shape of the elbow does not change
when $N_{c}$ is enlarged. This implies that the behavior does not indicate a
jump in the specific heat but rather a discontinuous change in the
differential coefficient, similar to what is observed in the case of the
circle graph. Therefore, it is reasonable to assume that the order of the
phase transition of the FKM model on $K_{4}$ is again third-order.
Analytically exploring the significance of the pointed elbow and how the phase
transition point shifts due to the interaction term of three cycles is
intriguing. These aspects will be discussed separately.

##  6 Numerical Results II: Irregular Graphs

In this section, we numerically examine the phase structure of the FKM model
on the irregular graphs, whose zeta functions do not possess exact duality.
However, we find that the behavior of the specific heat still has a similar
structure between the small and large coupling regions.

###  6.1 Double triangle

Figure 14: The specific heat of the FKM model on DT with varying
$\gamma=128,1024,16384$ and fixing $N_{c}=16$. The solid lines in $s>1$ and
$s<0$ represent the predictions from the effective actions (4.53) and (4.55)
for $\gamma=16384$, respectively. The dashed line in $s<0$ region is just a
flip of the prediction of the GWW model in $s>1$.

The double triangle (DT) is depicted in Fig. 1 (b). As discussed in Sec. 3.2,
if we define $q=(R_{\rm DT})^{s}$ with $R_{\rm DT}=0.65729\cdots$, the
critical strip appears at $0<\operatorname{Re}s<1$. Since the model is
unstable in the critical strip, we need to avoid this region in the numerical
simulation.

The numerical result of the specific heat is shown in Fig. 14. The solid line
in the region $s>1$ expresses the specific heat derived from the effective
action (4.53) in the small $q$ region. It well fits the numerical results. At
first glance, one would see that the duality also holds in DT, but it is not
correct. Indeed, if we simply flip the prediction of the GWW model from the
small coupling region ($s>1$) to the large coupling region ($s<0$) as drawn in
dashed line, the numerical results in $s<0$ are slightly different from it.
The specific heat derived from the effective action (4.55) in the large $q$
region which is shown as the solid line in $s<0$ in Fig. 14. It well agrees
with the numerical results.

Thus, we can conclude that the duality does not exactly hold and the
asymptotic expansion in the dual side is sufficiently accurate to explain the
numerical result in the FKM model on DT.

###  6.2 Triangle-square

Figure 15: The specific heat of the FKM model on TS with varying
$\gamma=128,1024,16384$ and fixing $N_{c}=16$. The solid lines in $s>1$ and
$s<0$ represent the predictions from the effective actions (4.63) and (4.65)
for $\gamma=16384$, respectively. The dashed line in $s<0$ region is just a
flip of the prediction of the GWW model in $s>1$.

The triangle-square (TS) contains two fundamental cycles of different lengths
three and four (See Fig. 6). We plot the specific heat of the FKM model on TS
in Fig. 15. Due to the difference of the length of two fundamental cycles, two
phase transitions appear in small $q$ region ($s>1$). The numerical results in
$s>1$ well fit to the prediction of the sum of the GWW models of two different
coupling constants (4.63), which is depicted by a solid line, as expected.

The duality between the small and large $q$ regions is broken as well as the
DT case. Indeed, if we flip the prediction of the GWW in the small $q$ region
to the large $q$ region as drawn in dashed line, it slightly does not agree
with the numerical results in $s<0$. This slight difference comes from the
sub-leading terms of the effective action in the large $q$ region (4.65). The
numerical results rather agree with the specific heat predicted by (4.65)
drawn in a solid line in the figure.

(a) $s=7.35277$

(b) $s=9.79308$

Figure 16: The eigenvalue distributions of $U_{1}$ and $U_{2}$ at the phase
transition points of $N_{c}=16$ and $\gamma=16384$, which are associated with
the fundamental cycles of TS with length three and four, respectively. The
eigenvalue distribution of $U_{1}$ in (a) saturates the upper limit of the
figure since it is still too sharp and narrow in the deeper deconfinement
phase.

The exsistence of the two different phase transitions can be also seen from
the eigenvalue distributions of the fundamental cycles. Fig. 16 is the
eigenvalue distributions of $U_{1}$ and $U_{2}$, which associated with the
fundamental cycles of length three and four, respectively, at the two phase
transition points in $s>1$ at $N_{c}=16$ and $\gamma=16384$. We see that the
eigenvalue distribution of $U_{2}$ exhibits the GWW phase transition at
$s=7.35277$, while $U_{1}$ is still in the deconfinement phase. At the second
phase transition point of $s=9.79308$, the eigenvalue distribution of $U_{1}$
exhibits the transition. Therefore, there exists a mixed phase of the
deconfinement and confinement between $s=7.35277$ and $s=9.79308$.

###  6.3 Triple triangle

Figure 17: The specific heat of the FKM model on TT with varying
$\gamma=128,1024,16384,131072$ and fixing $N_{c}=16$. The solid lines in the
regions $s>1$ and $s<0$ are the predictions from the effective actions (4.72)
and (4.74), respectively. The dashed line in the region $s<0$ is just a
reflection of the solid line in $s>1$.  Figure 18: An enlarged view near the
phase transition points for $s<0$ shown in Fig. 17.

We plot the specific heat of the FKM model on triple triangle (TT) in Fig. 17.
As depicted in Fig. 7, TT consists of three triangles connected together.
Contrary to its apparent simplicity, the FKM model on TT exhibits an
interesting behavior.

In the small coupling region $(s>1)$, we see only one phase transition and the
numerical results well fit the prediction of the effective action (4.72) drawn
by the solid line in $s>1$. On the other hand, it is observed that there are
two phase transitions in the large coupling region $(s<0)$. This is expected
from the $1/q$ expansion of the effective action in the large coupling reagion
(4.74). The small and large coupling regions are obviously not dual.

The difference of the two transition points in $s<0$ is smaller than that of
TS shown in Fig. 15. This is due to the fact that, in the effective action of
TT, the phase transition points are separated not by the difference in the
powers of $q$ in the coupling constants of the two GWW models but just by the
difference in the coefficients, while, in the effective action of TS, they
separate by the difference in the powers of $q$. Since the coupling constants
of the GWW models that appear in the effective action are proportional to
$\gamma$, we need large $\gamma$ to clearly see the separation of phase
transition points in the large $q$ region of TT. This is why we have added the
results of $\gamma=131072$ in Fig. 15. To see the separation of the phase
transition points clearly, we zoom in around the two phase transition points
at $s<0$ in Fig. 18. From this figure, we observe that the numerical results
for $s<0$ do not align with the dashed line, which reflects the solid line for
$s>1$. Instead, they align with the solid line, representing the effective
action’s prediction in the large $q$ region as shown in equation (4.74).

(a) $s=6.67970$

(b) $s=-5.72085$

(c) $s=-5.34085$

Figure 19: The eigenvalue distributions of the FKM model on TT with $N_{c}=16$
and $\gamma=131072$. (a) The common phase transition point for all unitary
matrices in $s>1$. (b) The phase transition point of $U_{1}$ and $U_{2}$ in
$s<0$. (c) The phase transition point of $U_{3}$ in $s<0$.

Fig.19 shows the eigenvalue distributions of the fundamental Wilson loops
$U_{1}$, $U_{2}$ and $U_{3}$ with $N_{c}=16$ and $\gamma=131072$. In the small
coupling region at $s=6.67970$ shown in (a), all the fundamental Wilson loops
exhibits the GWW phase transition at the same time, while, in the large
coupling region, only $U_{3}$ behaves differently from $U_{1}$ and $U_{2}$,
resulting in two different phase transition points as shown in (b) and (c).
This is again expected from the effective actions (4.72) and (4.74). The
separated behavior in the region $s<0$ is due to the fact that the mid cycle
$C_{3}$ is asymmetric with $C_{1}$ and $C_{2}$. This asymmetry appears in the
sub-leading of the $q$-expansion in the small $q$ region, but in the large $q$
region it appears as a difference in the coefficient of leading of the
$1/q$-expansion. This is the reason why the phase transition point splits at
$s<0$. Therefore, although it cannot be seen at the current resolution, a tiny
separation of phase transition points in the small $q$ region is expected to
occur as well as in the large $q$ region.

##  7 Conclusion and Discussion

In this paper, we have investigated the phase structure and duality of the
fundamental Kazakov-Migdal (FKM) model on generic connected simple graphs. The
partition function of the FKM model is described by the graph zeta function
weighted by unitary matrices on the links [19]. In this paper, we have
concentrated on the case where the graph zeta function becomes the Ihara zeta
function $\zeta(q;U)$. The series expansion of the Ihara zeta function
guarantees that the effective action of the FKM model, $-\gamma
N_{c}\log\zeta(q;U)$, can be expressed as a summation of all possible Wilson
loops on the graph for small $q$. We found that the FKM model is unstable in
the critical strip of the Ihara zeta function. We demonstrated that the FKM
model on a regular graph has a small/large coupling duality associated with
the functional relation of the Ihara zeta function. We also found a dual
description of the Ihara zeta function on a general graph as a generalization
of the function relation. Using this description, we have shown that the FKM
model can be represented as the sum of all possible Wilson loops on the graph
even in the region $q>1$.

As discussed in [19], the FKM model in large $N_{c}$ undergoes the GWW phase
transition in general. When the fundamental cycles of minimal lengths are
uniquely determined, the phase structure of the model at large $\gamma$ can be
estimated from the results of the GWW model. However, in other cases, non-
trivial coupling terms cannot be ignored. We have checked the validity of this
insight by using the numerical simulations and have shown that the system
behaves as expected in the theoretical analysis. In particular, it is
interesting that the behavior of the model in large $q$ region coincides with
the theoretical prediction based on the dual description of the Ihara zeta
function. The origin of this duality lies in the fact that the Ihara zeta
function is represented through the sum of all possible Wilson loops. This
phenomenon is reminiscent of the duality of the string theory, and the FKM
model is expected to provide a new perspective on the relationship between the
string theory and the gauge theory.

While this paper primarily focuses on the Ihara zeta function, by
incorporating more general mass parameters, it is possible to explore a
broader class of graph zeta functions. If we set the mass parameter $m_{v}$ as
(2.16), the partition function of the FKM model is expressed by the Bartholdi
zeta function $\zeta_{G}(q,u)$ with the bump parameter $u$. It is
straightforward to extend the discussion in this paper to this case, where the
phase diagram should be drawn in $(q,u)$-plane. It will be discussed
separately in future.

As we have seen in Fig. 11, the behavior of the specific heat of the model on
$K_{4}$ is different from that of the other graphs, that is, the pointed
“elbow” appears at the phase transition point. As discussed already, it is
reasonable to think that the effect of the interaction term is not making the
order of the phase transition from third-order to second-order, but changing
the behavior of the specific heat with keeping the order of the phase
transition to be third-order. This observation also explains the small elbows
that appeared in the model on graphs other than $K_{4}$ for finite $\gamma$
(see Figs. 9,14,15 and 17). Although the effective actions of these models
reduce to the GWW model in large $\gamma$ where the elbow disappears and the
phase transition is third-order, the interaction terms derived from the
composite cycles remain finite for finite $\gamma$, which are expected to have
a similar effect to the term $U_{1}U_{2}U_{3}$ in the $K_{4}$ case. In order
to make this point clear, we should consider the unitary matrix model with the
action,

$S=\frac{N}{\lambda}\operatorname{Tr}\Bigl{[}U_{1}+U_{2}+U_{3}+g\,U_{1}U_{2}U_{3}+h.c.\Bigr{]}\,,$
(7.1)

in close. It will also be discussed separately in another paper in
preparation.

## Acknowledgments

The authors would like to thank S. Aoki, T. Misumi and T. Yoda for helpful
discussions. This work is supported in part by Grant-in-Aid for Scientific
Research (KAKENHI) (C), Grant Number 20K03934 (S. M.) and Grant Number
23K03423 (K. O.).

##  Appendix A Poles of the Ihara zeta function of a regular graph

In this appendix, we show that all poles of the Ihara zeta function
$\zeta_{G}(q)$ exist in the region $R_{G}\leq q\leq 1$ in general.

Let us recall that the Ihara zeta function is expressed as
$\zeta_{G}(q)=\det(1-qW)^{-1}$ where $W$ is the edge adjacency matrix (2.6).
Therefore, the poles of $\zeta_{G}(q)$ are given by the inverse of the
eigenvalues of $W$. In particular, the maximal convergence radius $R_{G}$ is
the inverse of the maximal eigenvalue of $W$.

From this fact, we can immediately conclude the inequality (3.13). The point
is that the matrix $W$ is a non-negative matrix. As a general property of a
non-negative matrix, the maximal eigenvalue $\lambda_{\rm max}$ of $W$
satisfies

$\min_{1\leq{\boldsymbol{e}}\leq
2n_{E}}\left(\sum_{{\boldsymbol{e}}^{\prime}=1}^{2n_{E}}W_{{\boldsymbol{e}}{\boldsymbol{e}}^{\prime}}\right)\leq|\lambda_{\rm
max}|\leq\max_{1\leq{\boldsymbol{e}}\leq
2n_{E}}\left(\sum_{{\boldsymbol{e}}^{\prime}=1}^{2n_{E}}W_{{\boldsymbol{e}}{\boldsymbol{e}}^{\prime}}\right)\,.$
(A.1)

See e.g. [30] for a proof. From the definition of the edge adjacency matrix
(2.6), $W$ satisfies

$\sum_{{\boldsymbol{e}}^{\prime}=1}^{2n_{E}}W_{{\boldsymbol{e}}{\boldsymbol{e}}^{\prime}}=\deg
s({\boldsymbol{e}})-1=t_{s({\boldsymbol{e}})}\,.$ (A.2)

Combining (A.1) and (A.2), we can conclude (3.13) since $R_{G}=|\lambda_{\rm
max}|^{-1}$.

To estimate the spectrum of $W$ in more detail, it is convenient to estimate
$WW^{T}$;

$\displaystyle\sum_{{\boldsymbol{f}}\in
E_{D}}W_{{\boldsymbol{e}}{\boldsymbol{f}}}W_{{\boldsymbol{e}}^{\prime}{\boldsymbol{f}}}$
$\displaystyle=\sum_{{\boldsymbol{f}}\in
E_{D}}\left(\delta_{t({\boldsymbol{e}})s({\boldsymbol{f}})-\delta_{{\boldsymbol{e}}{\boldsymbol{f}}^{-1}}}\right)\left(\delta_{t({\boldsymbol{e}}^{\prime})s({\boldsymbol{f}})-\delta_{{\boldsymbol{e}}^{\prime}{\boldsymbol{f}}^{-1}}}\right)$
$\displaystyle=\left(\deg
t(e)-2\right)\delta_{t({\boldsymbol{e}})t({\boldsymbol{e}}^{\prime})}+\delta_{{\boldsymbol{e}}{\boldsymbol{e}}^{\prime}}\,.$
(A.3)

If we define

$E_{v}\equiv\left\\{{\boldsymbol{e}}\in
E_{D}|t({\boldsymbol{e}})=v\right\\}\,,$ (A.4)

and arrange $E_{D}$ as the disjoint union of $\\{E_{v}|v\in V\\}$, $WW^{T}$
can be expressed as a block diagonal matrix,

$WW^{T}=\begin{pmatrix}F_{\deg v_{1}}&&\\\ &\ddots&\\\ &&F_{\deg
v_{n_{V}}}\end{pmatrix}\,,$ (A.5)

where $F_{d}$ is a matrix of size $d$ defined as

$F_{d}\equiv(d-2)\begin{pmatrix}1&\cdots&1\\\ \vdots&\ddots&\vdots\\\
1&\cdots&1\end{pmatrix}+\boldsymbol{1}_{d}\,.$ (A.6)

Since it is easy to show

${\rm Spec}\,F_{d}=\left\\{(d-1)^{2},1,\cdots,1\right\\}\,,$ (A.7)

we see that

${\rm Spec}\,WW^{T}=\left\\{(\deg v_{1}-1)^{2},\cdots,(\deg
v_{n_{V}}-1)^{2},1,\cdots,1\right\\}\,.$ (A.8)

In other words, the Rayleigh quotient of $WW^{T}$ with respect to a vector
$\vec{v}$,

$R(\vec{v})\equiv\frac{\vec{v}\cdot WW^{T}\vec{v}}{|\vec{v}|^{2}}\,,$ (A.9)

satisfies

$1\leq R(\vec{v})\leq\left(\max_{v\in V}\deg v-1\right)^{2}=t_{\rm
max}^{2}\,.$ (A.10)

If we choose $\vec{v}$ as an eigenvector of $W$ with the eigenvalue $\lambda$,
the corresponding Rayleigh quotient is $R(\vec{v})=\lambda^{2}$. Therefore,
$\lambda$ satisfies

$1\leq|\lambda|\leq t_{\rm max}\,,$ (A.11)

and then the poles of the Ihara zeta function are in the region $R_{G}\leq
q\leq 1$ since the maximal convergence radius is the norm of the minimal pole
of $\zeta_{G}(q)$ by definition.

## References

  * [1] K.G. Wilson, _Confinement of quarks_ , _Phys. Rev. D_ 10 (1974) 2445.
  * [2] S. Matsuura and K. Ohta, _Kazakov-Migdal model on the graph and Ihara zeta function_ , _JHEP_ 09 (2022) 178 [2204.06424].
  * [3] S. Matsuura and K. Ohta, _Graph zeta functions and Wilson loops in a Kazakov–Migdal model_ , _Progress of Theoretical and Experimental Physics_ 2022 (2022) 123B03.
  * [4] V. Kazakov and A. Migdal, _Induced gauge theory at large n_ , _Nuclear Physics B_ 397 (1993) 214.
  * [5] Y. Ihara, _On discrete subgroups of the two by two projective linear group over $p$-adic fields_, _J. Math. Soc. Japan_ 18 (1966) 219.
  * [6] J.-P. Serre, _Trees_ , Springer-Verlag, Berlin (1980).
  * [7] T. Sunada, _L-functions in geometry and some applications_ , in _Curvature and topology of Riemannian manifolds_ , pp. 266–284, Springer (1986).
  * [8] K. Hashimoto, _ON ZETA AND L-FUNCTIONS OF FINITE GRAPHS_ , _International Journal of Mathematics_ 01 (1990) 381.
  * [9] H. Bass, _The ihara-selberg zeta function of a tree lattice_ , _International Journal of Mathematics_ 3 (1992) 717.
  * [10] L. Bartholdi, _Counting paths in graphs_ , _arXiv preprint math/0012161_ (2000) .
  * [11] H. Mizuno and I. Sato, _Bartholdi zeta functions of digraphs_ , _European Journal of Combinatorics_ 24 (2003) 947.
  * [12] A. Terras, _Zeta Functions of Graphs: A Stroll through the Garden_ , Cambridge Studies in Advanced Mathematics, Cambridge University Press (2010), 10.1017/CBO9780511760426.
  * [13] I.I. Kogan, G.W. Semenoff and N. Weiss, _Induced qcd and hidden local z n symmetry_ , _Physical review letters_ 69 (1992) 3435.
  * [14] I.I. Kogan, A. Morozov, G. Semenoff and N. Weiss, _Area law and continuum limit in “induced qcd”_ , _Nuclear Physics B_ 395 (1993) 547.
  * [15] I.I. Kogan, A. Morozov, G. Semenoff and N. Weiss, _Continuum limits of “induced qcd”: Lessons of the gaussian model at d= 1 and beyond_ , _International Journal of Modern Physics A_ 8 (1993) 1411.
  * [16] A. Migdal, _Bose condensation and zn symmetry breaking in the mixed model of induced qcd_ , _Modern Physics Letters A_ 8 (1993) 259.
  * [17] J. Cline and S. Paban, _Is induced qcd really qcd? the preservation of asymptotic freedom by matter interactions_ , _Physics Letters B_ 316 (1993) 158.
  * [18] B. Balakrishna, _Difficulties in inducing a gauge theory at large n_ , _Physical Review D_ 49 (1994) R24.
  * [19] S. Matsuura and K. Ohta, _Gross-Witten-Wadia phase transition in induced QCD on the graph_ , _Phys. Rev. D_ 108 (2023) 054504.
  * [20] I.Y. Arefeva, _A Simple model inducing QCD_ , _Phys. Lett. B_ 308 (1993) 347 [hep-th/9301065].
  * [21] D.J. Gross and E. Witten, _Possible Third Order Phase Transition in the Large N Lattice Gauge Theory_ , _Phys. Rev. D_ 21 (1980) 446.
  * [22] S.R. Wadia, _$N$ = Infinity Phase Transition in a Class of Exactly Soluble Model Lattice Gauge Theories_, _Phys. Lett. B_ 93 (1980) 403.
  * [23] M. Creutz, _Gauge fixing, the transfer matrix, and confinement on a lattice_ , _Phys. Rev. D_ 15 (1977) 1128.
  * [24] H.M. Stark, _Multipath zeta functions of graphs_ , in _Emerging Applications of Number Theory_ , pp. 601–615, Springer (1999).
  * [25] M. Kotani and T. Sunada, _Zeta functions of finite graphs_ , _Journal of Mathematical Sciences-the University of Tokyo_ 7 (2000) 7.
  * [26] K. Ohta and S. Matsuura, _Supersymmetric gauge theory on the graph_ , _PTEP_ 2022 (2022) 043B01 [2111.00676].
  * [27] J. Yumoto and T. Misumi, _Lattice fermions as spectral graphs_ , _Journal of High Energy Physics_ 2022 (2022) 1.
  * [28] J. Yumoto and T. Misumi, _New conjecture on exact dirac zero-modes of lattice fermions_ , _Progress of Theoretical and Experimental Physics_ 2023 (2023) 093B01.
  * [29] J. Yumoto and T. Misumi, _Equivalence of lattice operators and graph matrices_ , _Progress of Theoretical and Experimental Physics_ (2024) ptae009.
  * [30] R.A. Horn and C.R. Johnson, _Matrix analysis_ , Cambridge university press (2012).