#  $OnthecomplexityofCayleygraphonadihedralgroup.$ ††thanks: Supported by

###### Abstract

In this paper, we investigate the complexity of an infinite family of Cayley
graphs
$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}})$
on the dihedral group $\mathbb{D}_{n}=\langle
a,b|a^{2}=1,b^{n}=1,(a\,b)^{2}=1\rangle$ of order $2n.$

We obtain a closed formula for the number $\tau(n)$ of spanning trees in
$\mathcal{D}_{n}$ in terms of Chebyshev polynomials, investigate some
arithmetical properties of this function, and find its asymptotics as
$n\to\infty.$ Moreover, we show that the generating function
$F(x)=\sum\limits_{n=1}^{\infty}\tau(n)x^{n}$ is a rational function with
integer coefficients.

AMS classification: 05C30, 39A12
Keywords: spanning tree, Cayley graph, dihedral group, Chebyshev polynomial

On the complexity of Cayley graphs on a dihedral group

Bobo Hua,111Fudan University<EMAIL_ADDRESS>A. D. Mednykh,222Sobolev
Institute of Mathematics, Novosibirsk State University<EMAIL_ADDRESS>I. A.
Mednykh,333Sobolev Institute of Mathematics, Novosibirsk State University,
<EMAIL_ADDRESS>Lili Wang,444School of Mathematics and Statistics & Key
Laboratory of Analytical Mathematics and Applications (Ministry of Education)
& Fujian Provincial Key Laboratory of Statistics and Artificial Intelligence &
Fujian Key Laboratory of Analytical Mathematics and Applications (FJKLAMA) &
Center for Applied Mathematics of Fujian Province (FJNU), Fujian Normal
University<EMAIL_ADDRESS>
## Introduction

Let $G$ be a finite connected graph. The notion of the complexity of a graph
can be defined in several different ways. One can consider the number of edges
or vertices, the number of spanning trees or rooted spanning forests. The
important characteristic of a graph is the Kirchhoff index defined as the sum
of resistance distances between vertices. All the above-mentioned values can
be expressed in terms of the Laplacian spectrum of a graph. In particular, by
the famous Kirchhoff Matrix-Tree Theorem the number of spanning trees in a
connected graph is equal to the product of all non-zero eigenvalues of its
Laplacian matrix divided by the number of vertices.

The study of such invariants usually leads to the following question: how to
find the product of eigenvalues of the Laplacian matrix? If the size (number
of vertices) of a graph is small, it is an easy task. However, the most
interesting cases involve the family of graphs with increasing number of
vertices. The direct calculation of this product becomes tedious and
unmanageable when the number of vertices $n$ of the graph tends to infinity.
To solve this problem, we use the techniques developed in previous papers by
the authors [8], [10] and [14].

As a result, one can find a closed formula which is the product of a bounded
number of factors, each given by the $n$-th Chebyshev polynomial of the first
kind evaluated at the roots of some polynomial of prescribed degree. This
paves the way to investigate arithmetical properties and asymptotics.

The complexity of a graph plays an important role in statistic physics, where
the graphs with arbitrarily large number of vertices are considered ([9],
[20], [21]). With increasing number of vertices the structure of the Laplacian
characteristic polynomial becomes quite complicated. In this case, the most
interesting invariants are given by their asymptotics. See surveys for
counting spanning trees [4], [18] and the references therein.

The aim of the present paper is to produce explicit analytic formulas for the
number of spanning trees in a Cayley graph $\mathcal{D}_{n},$ see (1) for the
definition, on a dihedral group. Such formulas help in the investigation of
number–theoretical properties and asymptotics of several spectral invariants
of the graph. This research expands series of publications by various authors
on the complexity of circulant graphs ([1], [7], [8], [12], [15], [16], [22],
[23], [24]). Note that the circulant graph is a Cayley graph on a cyclic
group.

The methods laid out in the present paper can be equally used to find explicit
formulas for the number of rooted spanning forests and Kirchhoff index.

The paper is organized as follows. Some preliminary results and basic
definitions are given in Section 1. Section 2 is devoted to the investigation
of spectral properties of the Laplacian matrix of the Cayley graph
$\mathcal{D}_{n}$ on the dihedral group $\mathbb{D}_{n}.$ In Section 3, we
present explicit formulas for the number $\tau(n)$ of spanning trees in graph
$\mathcal{D}_{n}.$ The formulas will be given in terms of Chebyshev
polynomials evaluated in the roots of Laurent polynomial associated with
$\mathcal{D}_{n}.$ In Section 4, we provide some arithmetical properties of
function $\tau(n)$ for the family $\mathcal{D}_{n}.$ More precisely, we show
that the number of spanning trees can be represented in the form
$\tau(n)=p\,n\,a(n)^{2},$ where $a(n)$ is an integer sequence and $p$ is a
prescribed number depending only on parameters of $\mathcal{D}_{n}$ and the
parity of $n.$ In Section 5, we use explicit formulas for the complexity to
produce its asymptotic. Section 6 describes the structure of the generation
function $F(x)=\sum_{n=1}^{\infty}\tau(n)x^{n}.$ We show that $F(x)$ is a
rational function with integer coefficients satisfying some symmetry property.
In the last section, we illustrate the obtained results by a series of
examples.

##  1 Basic definitions and preliminary facts

Consider a finite connected graph $G$ with possibly multiple edges, but
without loops. We denote by $V(G)$ and $E(G)$ the set of vertices and the set
of edges of $G$ respectively.

A tree is a connected undirected graph without cycles. A spanning tree in a
graph $G$ is a subgraph that is a tree and contains all the vertices of $G.$

Given $u,v\in V(G),$ we write $a_{uv}$ as the number of edges between vertices
$u$ and $v.$ The matrix $A=A(G)=\\{a_{uv}\\}_{u,v\in V(G)}$ is called the
adjacency matrix of the graph $G.$ The degree $d_{v}$ of a vertex $v\in V(G)$
is defined by $d_{v}=\sum_{u\in V(G)}a_{uv}.$ Let $D=D(G)$ be the diagonal
matrix indexed by the elements of $V(G)$ with $d_{vv}=d_{v}.$ The matrix
$L=L(G)=D(G)-A(G)$ is called the Laplacian matrix, or simply Laplacian, of the
graph $G.$

In what follows, we denote by $I_{n}$ the identity matrix of order $n.$

We say that an $n\times n$ matrix is circulant, denoted by
$circ(a_{0},a_{1},\ldots,a_{n-1}),$ if it is of the form

$circ(a_{0},a_{1},\ldots,a_{n-1})=\left(\begin{array}[]{ccccc}a_{0}&a_{1}&a_{2}&\ldots&a_{n-1}\\\
a_{n-1}&a_{0}&a_{1}&\ldots&a_{n-2}\\\ &\vdots&&\ddots&\vdots\\\
a_{1}&a_{2}&a_{3}&\ldots&a_{0}\\\ \end{array}\right).$

Recall [5] that the eigenvalues of the matrix
$C=circ(a_{0},a_{1},\ldots,a_{n-1})$ are given by the following simple
formulas $\lambda_{j}=p(\varepsilon^{j}_{n}),\,j=0,1,\ldots,n-1$ where
$p(x)=a_{0}+a_{1}x+\ldots+a_{n-1}x^{n-1}$ and $\varepsilon_{n}$ is an order
$n$ primitive root of the unity. For any $i=0,\ldots,n-1$, let
$w_{i}=(1,\varepsilon_{n}^{i},\varepsilon_{n}^{2i},\ldots,\varepsilon_{n}^{(n-1)i})^{t}$
be a column vector of length $n.$ Then all $n\times n$ circulant matrices
share the same set of linearly independent eigenvectors
$w_{0},w_{1},\ldots,w_{n-1}.$ Hence, any set of $n\times n$ circulant matrices
can be simultaneously diagonalizable. Moreover, let
$T_{n}=circ(0,1,0,\ldots,0)$ be the matrix representation of the shift
operator
$T_{n}:(x_{0},x_{1},\ldots,x_{n-2},x_{n-1})\rightarrow(x_{1},x_{2},\ldots,x_{n-1},x_{0}).$
Then
$T_{n}w_{i}=\varepsilon_{n}w_{i},\,T_{n}^{-1}w_{i}=\varepsilon_{n}^{-1}w_{i},$
and for any Laurent polynomial $P(z)$ one has
$P(T_{n})w_{i}=P(\varepsilon_{n})w_{i}.$

Let $D$ be a group, and let $S$ be a subset of $D,$ which doesn’t contain the
identity element $1.$ The Cayley digraph associated with $(D,S)$ is then
defined as the directed graph with the set of vertices $D$ and the set of
edges

$\\{(g,h):g,h\in D,gh^{-1}\in S\\}.$

The Cayley graph depends on the choice of a generating set $S,$ and is
connected if and only if $S$ generates $D$ (i.e., the set $S$ are group
generators of $D$). We deal with undirected graphs, and always assume that
$S=S^{-1},$ where $S^{-1}=\\{s:s^{-1}\in S\\}.$

Let $\mathbb{D}_{n}=\langle a,b|a^{2}=1,b^{n}=1,(a\,b)^{2}=1\rangle$ be
dihedral group of order $2n.$ We arrange the elements of the group
$\mathbb{D}_{n}$ as $V=\\{1,b,\ldots,b^{n-1},a,ba,\ldots,b^{n-1}a\\}$ and
consider the Cayley graph

$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}})$
(1)

with the generating set
$S=\\{b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}}\\}$
for some integers $\beta_{1},\beta_{2},\cdots,$ and
$\gamma_{1},\gamma_{2},\cdots.$ We suppose that $\mathbb{D}_{n}$ acts on $V$
by the rule: $g\in\mathcal{D}_{n}$ sends a vertex $v\in V$ to the vertex $vg.$
Then the set of oriented edges of $\mathcal{D}_{n}$ can be describe as
follows. Given $j\in\\{\pm\beta_{1},\pm\beta_{2},\ldots,\pm\beta_{s}\\}$ there
is an edge $b^{k}\xrightarrow{b^{j}}b^{k+j}$ and an edge
$b^{k}a\xrightarrow{b^{j}}b^{k-j}a$ for any $k=0,1,\ldots,n-1;$ given
$j\in\\{\gamma_{1},\gamma_{2},\ldots,\gamma_{t}\\}$ there is an edge
$b^{k}\xrightarrow{ab^{j}}b^{k-j}a$ and an edge
$b^{k}a\xrightarrow{ab^{j}}b^{k+j}$ for any $k=0,1,\ldots,n-1.$ Noting that
$(ab^{j})^{-1}=b^{-j}a=ab^{j},$ we have $S=S^{-1}.$ Hence $\mathcal{D}_{n}$ is
an undirected graph.

We always restrict to the case
$0<\beta_{1}<\beta_{2}<\cdots<\beta_{s}<\frac{n}{2}$ and
$0\leq\gamma_{1}<\gamma_{2}<\cdots<\gamma_{t}\leq{n-1}.$ We suppose that
$s\geq 0$ and $t\geq 1.$ Then the graph $\mathcal{D}_{n}$ has no loops and
multiple edges.

Now we introduce a necessary and sufficient condition for the connectedness of
the graph $\mathcal{D}_{n}.$

###### Lemma 1.1

The graph
$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}})$
is connected if and only if $\gcd(n,\,\beta_{j},\,1\leq j\leq
s,\,\gamma_{j}-\gamma_{k},\,1\leq j<k\leq t)$ is equal to $1.$

Proof: Since our graph $\mathcal{D}_{n}$ is a Cayley graph associated with
$(\mathbb{D}_{n},\,S),$ it is connected if and only if the set $S$ generates
group $\mathbb{D}_{n}.$

Suppose that $\gcd(n,\,\beta_{j},\,1\leq j\leq
s,\,\gamma_{j}-\gamma_{k},\,1\leq j<k\leq t)=1.$ Then, there exist integer
numbers $\\{c_{j}\\}_{1\leq j\leq s},\,\\{d_{j,k}\\}_{1\leq j<k\leq t},\,e$
such that

$\sum_{j=1}^{s}c_{j}\beta_{j}+\sum_{1\leq j<k\leq
t}d_{j,k}(\gamma_{j}-\gamma_{k})+e\,n=1.$

(See, for example, [2], p. 21). One easily sees that

$b=b^{\sum_{j=1}^{s}c_{j}\beta_{j}+\sum_{1\leq j<k\leq
t}d_{j,k}(\gamma_{j}-\gamma_{k})+e\,n}=b^{\sum_{j=1}^{s}c_{j}\beta_{j}}b^{\sum_{1\leq
j<k\leq t}d_{j,k}(\gamma_{j}-\gamma_{k})}$

We rewrite $b^{\sum_{1\leq j<k\leq t}d_{j,k}(\gamma_{j}-\gamma_{k})}$ in the
form

$b^{\sum_{1\leq j<k\leq t}d_{j,k}(\gamma_{j}-\gamma_{k})}=\prod_{1\leq j<k\leq
t}\underbrace{b^{-\gamma_{k}}b^{\gamma_{j}}\cdot\ldots\cdot
b^{-\gamma_{k}}b^{\gamma_{j}}}_{d_{j,k}\text{times}}=$ $\prod_{1\leq j<k\leq
t}\underbrace{b^{-\gamma_{k}}\,a\,a\,b^{\gamma_{j}}\cdot\ldots\cdot
b^{-\gamma_{k}}\,a\,a\,b^{\gamma_{j}}}_{d_{j,k}\text{times}}=\prod_{1\leq
j<k\leq t}(a\,b^{\gamma_{k}}\cdot a\,b^{\gamma_{j}})^{d_{j,k}}.$

In other words, we can express the generator $b$ in terms of elements of $S$

$b=\prod_{j=1}^{s}(b^{\beta_{j}})^{c_{j}}\prod_{1\leq j<k\leq
t}(a\,b^{\gamma_{k}}\cdot a\,b^{\gamma_{j}})^{d_{j,k}}.$

The set $S$ contains a term $a\,b^{\gamma_{1}}.$ So that
$a=a\,b^{\gamma_{1}}\cdot b^{-\gamma_{1}}$ also is the product of elements of
$S.$

On the other hand, let $\mathcal{D}_{n}$ be connected. Then the set $S$
generates full group $\mathbb{D}_{n}.$ Hence, there exist words in $S$
representing $b.$ We can always choose the representation of $b$ in the form

$b=\prod_{j=1}^{s}(b^{\beta_{j}})^{c_{j}}\prod_{k=1}^{N}a\,b^{\gamma_{i_{k}}},$

where $N\geq 0,\ c_{j}$ are integers for all $j=1,\ldots,s$ and each $i_{k}$
lies in the set $\\{1,2,\ldots,t\\}$ for all $k=1,\ldots,N.$ In fact, using
the fact that any element $a\,b^{\gamma_{j}}$ is an involution and the
following equalities in $\mathbb{D}_{n}$

$a\,b^{\gamma_{i}}\cdot b^{\beta_{j}}=a\,b^{\beta_{j}}\cdot
b^{\gamma_{i}}=b^{-\beta_{j}}\cdot a\,b^{\gamma_{i}},$

one can reduce any word in terms of $S$ to a given form.

We note that $N$ is even. Otherwise, if $N$ is odd then using above arguments
we deduce that $a$ is a proper power of $b.$ This is impossible in the
dihedral group. So

$b=\prod_{j=1}^{s}(b^{\beta_{j}})^{c_{j}}\prod_{k=1}^{N}a\,b^{\gamma_{i_{k}}}=\prod_{j=1}^{s}(b^{\beta_{j}})^{c_{j}}\prod_{k=1}^{N/2}a\,b^{\gamma_{i_{2k-1}}}\cdot
a\,b^{\gamma_{i_{2k}}}=$
$\prod_{j=1}^{s}(b^{\beta_{j}})^{c_{j}}\prod_{k=1}^{N/2}b^{-\gamma_{i_{2k-1}}}\,a\cdot
a\,b^{\gamma_{i_{2k}}}=\prod_{j=1}^{s}(b^{\beta_{j}})^{c_{j}}\prod_{k=1}^{N/2}b^{-\gamma_{i_{2k-1}}+\gamma_{i_{2k}}}.$

So that

$b=b^{\sum_{j=1}^{s}{\beta_{j}}{c_{j}}+\sum_{k=1}^{N/2}(\gamma_{i_{2k}}-\gamma_{i_{2k-1}})}.$

Or, (since $b$ is an element of order $n$)

$\sum_{j=1}^{s}c_{j}\beta_{j}+\sum_{k=1}^{N/2}(\gamma_{i_{2k}}-\gamma_{i_{2k-1}})\equiv
1(\hbox{\rm mod}\;n).$

Hence, there exists an integer $e$ such that

$\sum_{j=1}^{s}c_{j}\beta_{j}+\sum_{k=1}^{N/2}(\gamma_{i_{2k}}-\gamma_{i_{2k-1}})+e\,n=1.$

Therefore, $\gcd(n,\,\beta_{j},\,1\leq j\leq s,\,\gamma_{j}-\gamma_{k},\,1\leq
j<k\leq t)=1$ and the lemma is proved.

The adjacency matrix of the graph $\mathcal{D}_{n}$ is given by the $2n\times
2n$ block matrix

$A=\left(\begin{array}[]{cc}\sum\limits_{j=1}^{s}(T_{n}^{\beta_{j}}+T_{n}^{-\beta_{j}})&\sum\limits_{j=1}^{t}T_{n}^{-\gamma_{j}}\\\
\sum\limits_{j=1}^{t}T_{n}^{\gamma_{j}}&\sum\limits_{j=1}^{s}(T_{n}^{\beta_{j}}+T_{n}^{-\beta_{j}})\end{array}\right),$

where $T_{n}=circ(0,1,0,\ldots,0).$ The corresponding degree matrix is
$D=\left(\begin{array}[]{cc}(2s+t)I_{n}&0\\\ 0&(2s+t)I_{n}\end{array}\right),$
where $I_{n}$ is the $n\times n$ identity matrix. Since the Laplacian of
$\mathcal{D}_{n}$ is $L=D-A,$ we have

$L=\left(\begin{array}[]{cc}(2s+t)I_{n}-\sum\limits_{j=1}^{s}(T_{n}^{\beta_{j}}+T_{n}^{-\beta_{j}})&-\sum\limits_{j=1}^{t}T_{n}^{-\gamma_{j}}\\\
-\sum\limits_{j=1}^{t}T_{n}^{\gamma_{j}}&(2s+t)I_{n}-\sum\limits_{j=1}^{s}(T_{n}^{\beta_{j}}+T_{n}^{-\beta_{j}})\end{array}\right).$

##  2 Spectrum of a Cayley graph on the dihedral group

Let
$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}})$
be the Cayley graph on a dihedral group $\mathbb{D}_{n}.$ We introduce the
following Laurent polynomials

${\cal A}(z)=2s+t-\sum\limits_{i=1}^{s}(z^{\beta_{i}}+z^{-\beta_{i}}),\,{\cal
B}(z)=-\sum\limits_{i=1}^{t}z^{\gamma_{i}}\text{ and }P(z)={\cal A}(z){\cal
A}(z^{-1})-{\cal B}(z){\cal B}(z^{-1}).$

We will refer to $P(z)$ as the Laurent polynomial associated with the graph
$\mathcal{D}_{n}.$ Then the Laplacian of the graph $\mathcal{D}_{n}$ is given
by the following $2n\times 2n$ block matrix $L=\left(\begin{array}[]{cc}{\cal
A}(T)&{\cal B}(T^{-1})\\\ {\cal B}(T)&{\cal A}(T^{-1})\end{array}\right),$
where $T=circ(0,1,0,\ldots,0).$

Recall that $T$ is conjugated to
$\mathbb{T}=diag(1,\varepsilon,\ldots,\varepsilon^{n-1}),$ where
$\varepsilon=\varepsilon_{n}=\exp(2\pi\mathtt{i}/n).$ So that $L$ is
conjugated to the matrix $\mathbb{L}=\left(\begin{array}[]{cc}{\cal
A}(\mathbb{T})&{\cal B}(\mathbb{T}^{-1})\\\ {\cal B}(\mathbb{T})&{\cal
A}(\mathbb{T}^{-1})\end{array}\right).$ The spectra of $L$ and $\mathbb{L}$
are the same. It can be found by solving the system of linear equations

$\begin{cases}{\cal A}(\mathbb{T})\textbf{x}+{\cal
B}(\mathbb{T}^{-1})\textbf{y}=\lambda\textbf{x}\\\ {\cal
B}(\mathbb{T})\textbf{x}+{\cal
A}(\mathbb{T}^{-1})\textbf{y}=\lambda\textbf{y},\end{cases}$

where $\textbf{x},\textbf{y}\in\mathbb{C}^{n}$ and
$(\textbf{x},\textbf{y})\neq(\textbf{0},\textbf{0}).$ Since the matrix
$\mathbb{T}$ is diagonal, the system of equations splits into $n$ scalar
linear systems

$\begin{cases}{\cal A}(\varepsilon^{j})x+{\cal B}(\varepsilon^{-j})y=\lambda
x\\\ {\cal B}(\varepsilon^{j})x+{\cal A}(\varepsilon^{-j})y=\lambda
y,\end{cases}$

where $j=0,1\ldots,n-1.$ Hence, $\lambda$ is a root of the quadratic equation

$\lambda^{2}-({\cal A}(\varepsilon^{j})+{\cal
A}(\varepsilon^{-j}))\lambda+{\cal A}(\varepsilon^{j}){\cal
A}(\varepsilon^{-j})-{\cal B}(\varepsilon^{j}){\cal B}(\varepsilon^{-j})=0.$

The solutions of this equation are $\lambda_{j,1}=\Re({\cal
A}(\varepsilon^{j}))+\sqrt{-\Im({\cal A}(\varepsilon^{j}))^{2}+|{\cal
B}(\varepsilon^{j})|^{2}}$ and $\lambda_{j,2}=\Re({\cal
A}(\varepsilon^{j}))-\sqrt{-\Im({\cal A}(\varepsilon^{j}))^{2}+|{\cal
B}(\varepsilon^{j})|^{2}}.$ We note that L. Babai [3] found similar formulas
for eigenvalues for the adjacency matrices of $\mathcal{D}_{n}$ by making use
of the representation theory for finite groups. The corresponding eigenvectors
of the operator $\mathbb{L}$ are $u_{j,1}={\cal
B}(\varepsilon^{j})\textbf{e}_{j+1}+(\lambda_{j,1}-{\cal
A}(\varepsilon^{j}))\textbf{e}_{j+1}$ and $u_{j,2}={\cal
B}(\varepsilon^{j})\textbf{e}_{j+1}+(\lambda_{j,2}-{\cal
A}(\varepsilon^{j}))\textbf{e}_{j+1},$ where $\textbf{e}_{j}$ is the $j$-th
basic vector in the $\mathbb{C}^{2n}.$

For the special case $j=0,$ we have $\lambda^{2}-2{\cal A}(1)\lambda+{\cal
A}(1)^{2}-{\cal B}(1)^{2}=0.$ Hence

$\lambda_{0,1}={\cal A}(1)+{\cal B}(1)=0,\,\lambda_{0,2}={\cal A}(1)-{\cal
B}(1)=2{\cal A}(1).$

We note also that $\lambda_{j,1}\lambda_{j,2}={\cal A}(\varepsilon^{j}){\cal
A}(\varepsilon^{-j})-{\cal B}(\varepsilon^{j}){\cal
B}(\varepsilon^{-j})=P(\varepsilon^{j}).$ We suppose that the graph
$\mathcal{D}_{n}$ is connected. Then all the Laplacian eigenvalues, except of
$\lambda_{0,1}=0,$ are non-zero. They are

$\lambda_{0,2},\,\lambda_{j,1},\,\lambda_{j,2},\,j=1,2,\ldots,n-1.$

By the Kirchhoff’s theorem we have the following formula for the number of
spanning trees in $\mathcal{D}_{n}.$

$\tau(n)=\frac{\lambda_{0,2}}{2n}\prod\limits_{j=1}^{n-1}\lambda_{j,1}\lambda_{j,2}=\frac{{\cal
A}(1)}{n}\prod\limits_{j=1}^{n-1}P(\varepsilon^{j}).$

Since ${\cal A}(1)=t,$ we get the following result.

###### Theorem 2.1

Let
$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}})$
be a Cayley graph on the dihedral group $\mathbb{D}_{n}.$ Then the number of
spanning trees of $\mathcal{D}_{n}$ is given by the formula

$\tau(n)=\frac{t}{n}\prod\limits_{j=1}^{n-1}P(\varepsilon_{n}^{j}),$

where $\varepsilon_{n}=\exp(2\pi\mathtt{i}/n),\,P(z)={\cal A}(z){\cal
A}(z^{-1})-{\cal B}(z){\cal B}(z^{-1}),\,{\cal
A}(z)=2s+t-\sum\limits_{i=1}^{s}(z^{\beta_{i}}+z^{-\beta_{i}})$ and ${\cal
B}(z)=-\sum\limits_{i=1}^{t}z^{\gamma_{i}}.$

The following two lemmas give more properties of the associated polynomial
$P(z).$

###### Lemma 2.2

Let the graph $\mathcal{D}_{n}$ be connected. Then $P(1)=0,P^{\prime}(1)=0$
and $P^{\prime\prime}(1)<0.$

Proof: By direct calculation we have $P(1)=P^{\prime}(1)=0$ and

$P^{\prime\prime}(1)=-4t\sum_{j=1}^{s}\beta_{j}^{2}+2(\sum_{j=1}^{t}\gamma_{j})^{2}-2t\sum_{j=1}^{t}\gamma_{j}^{2}.$

We use Lagrange’s identity

$(\sum_{j=1}^{t}\xi_{j}^{2})(\sum_{j=1}^{t}\zeta_{j}^{2})-(\sum_{j=1}^{t}\xi_{j}\zeta_{j})^{2}=\sum\limits_{1\leq
j<k\leq t}(\xi_{j}\zeta_{k}-\xi_{k}\zeta_{j})^{2}.$

Setting $\xi_{j}=1$ and $\zeta_{j}=\gamma_{j},$ we get

$P^{\prime\prime}(1)=-2(2t\sum_{j=1}^{s}\beta_{j}^{2}+\sum_{1\leq j<k\leq
t}(\gamma_{j}-\gamma_{k})^{2})<0.$

In fact, it is impossible that all elements $\\{\beta_{i}\\}_{i=1,\ldots,s}$
and $\\{\gamma_{j}-\gamma_{k}\\}_{j,k=1,\ldots,t}$ are equal to $0,$ since the
graph $\mathcal{D}_{n}$ is connected. See Lemma 1.1.

###### Lemma 2.3

Suppose that the numbers $\\{\beta_{j},\,1\leq j\leq
s,\,\gamma_{j}-\gamma_{k},\,1\leq j<k\leq t\\}$ are relatively prime. Then for
any $\varphi\in\mathbb{R}$, we have $P(e^{\mathtt{i}\varphi})\geq 0.$
Furthermore, $P(e^{\mathtt{i}\varphi})=0$ if and only if
$e^{\mathtt{i}\varphi}=1.$

Proof: We have $P(e^{\mathtt{i}\varphi})={\cal A}(e^{\mathtt{i}\varphi}){\cal
A}(e^{-\mathtt{i}\varphi})-{\cal B}(e^{\mathtt{i}\varphi}){\cal
B}(e^{-\mathtt{i}\varphi})=|{\cal A}(e^{\mathtt{i}\varphi})|^{2}-|{\cal
B}(e^{\mathtt{i}\varphi})|^{2}.$ Also,

$|A(e^{\mathtt{i}\varphi})|^{2}=(s+t-\sum\limits_{j=1}^{s}\cos\beta_{j}\varphi)^{2}+(\sum\limits_{j=1}^{s}\sin\beta_{j}\varphi)^{2}=(t+\sum\limits_{j=1}^{s}(1-\cos\beta_{j}\varphi))^{2}+(\sum\limits_{j=1}^{s}\sin\beta_{j}\varphi)^{2}\geq
t^{2}$

and

$|{\cal
B}(e^{\mathtt{i}\varphi})|^{2}=|\sum_{j=1}^{t}e^{\mathtt{i}\gamma_{j}}|^{2}=\sum\limits_{1\leq
j,\,k\leq t}\cos(\gamma_{j}-\gamma_{k})\varphi=t+2\sum\limits_{1\leq j<k\leq
t}\cos(\gamma_{j}-\gamma_{k})\varphi\leq t^{2}.$

Hence, $P(e^{\mathtt{i}\varphi})=|{\cal A}(e^{\mathtt{i}\varphi})|^{2}-|{\cal
B}(e^{\mathtt{i}\varphi})|^{2}\geq 0.$ We have the equality if and only if

$(t+\sum\limits_{j=1}^{s}(1-\cos\beta_{j}\varphi))^{2}+(\sum\limits_{j=1}^{s}\sin\beta_{j}\varphi)^{2}=t^{2}\text{
and }t^{2}=t+2\sum\limits_{1\leq j<k\leq
t}\cos(\gamma_{j}-\gamma_{k})\varphi.$

The first equation is equivalent to

$\cos\beta_{j}\varphi=1,1\leq j\leq s,$

and the second one is equivalent to

$\cos(\gamma_{j}-\gamma_{k})\varphi=1,\,1\leq j<k\leq t.$

Since the numbers $\\{\beta_{j},\,1\leq j\leq
s,\,\gamma_{j}-\gamma_{k},\,1\leq j<k\leq t\\}$ are relatively prime, we
obtain $\cos\varphi=1,$ or $e^{\mathtt{i}\varphi}=1.$

##  3 Counting spanning trees

The main results of this section paper are the following theorems.

###### Theorem 3.1

Let
$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}})$
be a Cayley graph on the group $\mathbb{D}_{n}$ and $P(z)$ be the associated
Laurent polynomial for $\mathbb{D}_{n}.$ Then the number of spanning trees
$\tau(n)$ in the graph $\mathcal{D}_{n}$ is given by the formula

$\tau(n)=\frac{n\,t\,|\eta|^{n}}{q}\prod\limits_{\begin{subarray}{c}P(z)=0\\\
z\neq 1\end{subarray}}|z^{n}-1|,$

where the product is taken over all the roots, different from $1,$ of the
associated Laurent polynomial $P(z),$ $\eta$ is the leading coefficient of
$P(z),\text{ and }q=2t\sum_{j=1}^{s}\beta_{j}^{2}+\sum_{1\leq j<k\leq
t}(\gamma_{j}-\gamma_{k})^{2}.$

Proof: By Theorem 2.1 we already have

$\tau(n)=\frac{t}{n}\prod\limits_{j=1}^{n-1}P(\varepsilon_{n}^{j}).$ (2)

Denote by $\eta$ the leading coefficient of $P(z).$ Since $P(z)=P(1/z)$ we can
present the polynomial in the form

$P(z)=\eta\,z^{-r}+a_{1}z^{-r+1}+\ldots+a_{r}+\ldots+a_{1}z^{r-1}+\eta\,z^{r},$

for some $r\geq 0.$ To continue the proof, we replace the Laurent polynomial
$P(z)$ by $\widetilde{P}(z)=\frac{z^{r}}{\eta}P(z).$ Then $\widetilde{P}(z)$
is a monic polynomial of the degree $2r$ with the same roots as $P(z).$ We
note that

$\prod\limits_{j=1}^{n-1}\widetilde{P}(\varepsilon_{n}^{j})=\frac{\varepsilon_{n}^{\frac{(n-1)n}{2}r}}{\eta^{n-1}}\prod\limits_{j=1}^{n-1}P(\varepsilon_{n}^{j})=\frac{(-1)^{r(n-1)}}{\eta^{n-1}}\prod\limits_{j=1}^{n-1}P(\varepsilon_{n}^{j}).$
(3)

By Lemma 2.2 the polynomial $\widetilde{P}(z)$ has two roots equal to $1$ and
all the other roots different from $1.$ Also, we recognize the complex numbers
$\varepsilon_{n}^{j},\,j=1,\ldots,n-1$ as the roots of polynomial
$\frac{z^{n}-1}{z-1}.$ By the basic properties of resultant we have ([19], Ch.
1.3)

$\displaystyle\prod\limits_{j=1}^{n-1}\widetilde{P}(\varepsilon_{n}^{j})=\textrm{Res}(\widetilde{P}(z),\frac{z^{n}-1}{z-1})=\textrm{Res}(\frac{z^{n}-1}{z-1},\widetilde{P}(z))=\prod\limits_{z:\,\widetilde{P}(z)=0}\frac{z^{n}-1}{z-1}$
$\displaystyle=\prod\limits_{z:\,P(z)=0}\frac{z^{n}-1}{z-1}=(\lim_{z\to
1}\frac{z^{n}-1}{z-1})^{2}\prod\limits_{\begin{subarray}{c}P(z)=0\\\ z\neq
1\end{subarray}}\frac{z^{n}-1}{z-1}=n^{2}\prod\limits_{\begin{subarray}{c}P(z)=0\\\
z\neq 1\end{subarray}}\frac{z^{n}-1}{z-1}.$ (4)

Combining (2), (3) with (3), we have the following formula for the number of
spanning trees

$\displaystyle\tau(n)=(-1)^{r(n-1)}\eta^{n-1}n\,t\prod\limits_{\begin{subarray}{c}P(z)=0\\\
z\neq 1\end{subarray}}\frac{z^{n}-1}{z-1}.$ (5)

To finish the proof we need to evaluate the product
$\prod\limits_{\begin{subarray}{c}P(z)=0\\\ z\neq
1\end{subarray}}(z-1)=\prod\limits_{\begin{subarray}{c}\widetilde{P}(z)=0\\\
z\neq 1\end{subarray}}(z-1).$ To do this, we use the following property of
resultants ([19], Ch. 1.3, property 3). If $B$ is a monic polynomial, then for
any polynomials $A$ and $C$ we have

$\textrm{Res}(B,A+CB)=\textrm{Res}(B,A).$

We use the Taylor expansion of polynomial $\widetilde{P}(z)$ at the point
$z=1,$ which has the form

$\widetilde{P}(z)=\widetilde{P}(1)+\widetilde{P}^{\prime}(1)(z-1)+\frac{\widetilde{P}^{\prime\prime}(1)}{2}(z-1)^{2}+\frac{\widetilde{P}^{\prime\prime\prime}(1)}{6}(z-1)^{3}+\ldots.$

Since $\widetilde{P}(1)=\widetilde{P}^{\prime}(1)=0,$ we have

$\prod\limits_{\begin{subarray}{c}\widetilde{P}(z)=0\\\ z\neq
1\end{subarray}}(z-1)=\textrm{Res}(z-1,\frac{\widetilde{P}(z)}{(z-1)^{2}})=\textrm{Res}(z-1,\frac{\widetilde{P}^{\prime\prime}(1)}{2}+\frac{\widetilde{P}^{\prime\prime\prime}(1)}{6}(z-1)+\ldots)=\textrm{Res}(z-1,\frac{\widetilde{P}^{\prime\prime}(1)}{2}).$

The last term is equal to
$\frac{\widetilde{P}^{\prime\prime}(1)}{2}=\frac{1}{2}(\frac{z^{r}}{\eta}P(z))^{\prime\prime}_{z=1}=\frac{P^{\prime\prime}(1)}{2\eta}=-\frac{q}{\eta}$
(see Lemma 2.2). As a result, we get

$\displaystyle\tau(n)=\frac{(-1)^{r(n-1)+1}n\,t\,\eta^{n}}{q}\prod\limits_{\begin{subarray}{c}P(z)=0\\\
z\neq 1\end{subarray}}(z^{n}-1).$ (6)

Since $\tau(n)$ is a positive integer, the statement of theorem follows.

Remark. According to Lemma 2.2, for the calculation of $q$ one can get the
formula $q=-\frac{1}{2}P^{\prime\prime}(1).$

The associated polynomial $P(z)$ is a palindromic Laurent polynomial,
satisfying the property $P(z)=P(1/z).$ In other words, $P(z)$ has the form
$P(z)=c_{0}+\sum_{k=1}^{r}c_{k}(z^{k}+z^{-k}).$ Denote by
$T(n,x)=\cos(n\arccos(x))$ the $n$-th Chebyshev polynomial of the first kind.
The following equality is known
$T(n,\frac{z+z^{-1}}{2})=\frac{z^{n}+z^{-n}}{2}.$ Other basic properties of
Chebyshev polynomials can be found in [13].

From here, we conclude that

$P(z)=Q(\frac{z+z^{-1}}{2}),\text{ where
}Q(w)=c_{0}+\sum_{k=1}^{r}2c_{k}T(k,w).$

We refer to the polynomial $Q(w)$ as the Chebyshev transform of $P(z).$ These
two polynomials are also related as $Q(w)=P(w+\sqrt{w^{2}-1}).$ It is easier
to deal with $Q(w)$ since it is an ordinary polynomial with degree twice less
than $P(z).$

The following theorem is a direct consequence of Theorem 3.1.

###### Theorem 3.2

Let $P(z)$ be the associated Laurent polynomial of the Cayley graph
$\mathcal{D}_{n}.$ Denote by $Q(w)$ the Chebyshev transform of $P(z)$ and let
$r$ be the degree of polynomial $Q(w).$ The number of spanning trees $\tau(n)$
in the graph $\mathcal{D}_{n}$ is given by the formula

$\tau(n)=\frac{n\,t|\eta|^{n}}{q}\prod_{p=1}^{r-1}|2T(n,w_{p})-2|,$

where $w_{p},\,p=1,2,\ldots,r-1$ are roots, different from $1,$ of the
algebraic equation $Q(w)=0,$ and $T(n,w)$ is the Chebyshev polynomial of the
first kind, $\eta$ is the leading coefficient of $P(z),$ and
$q=2t\sum_{j=1}^{s}\beta_{j}^{2}+\sum_{1\leq j<k\leq
t}(\gamma_{j}-\gamma_{k})^{2}.$

Proof: Since $P(z)=Q((z+1/z)/2),$ by Lemma 2.2, the roots of polynomials
${P}(z)$ and $Q(w)$ are
$1,1,z_{1},1/z_{1},\ldots,z_{r-1},1/z_{r-1},\,z_{j}\neq 1\textrm{ and
}1,\,\frac{1}{2}(z_{j}+z_{j}^{-1})=w_{j}\neq 1,\,j=1,\ldots,r-1,$
respectively. Also $T(n,w_{j})=\frac{z_{j}^{n}+z_{j}^{-n}}{2}.$ Hence,

$\prod\limits_{\begin{subarray}{c}P(z)=0\\\ z\neq
1\end{subarray}}(z^{n}-1)=\prod\limits_{j=1}^{r-1}({z_{j}^{n}-1})({z_{j}^{-n}-1})=\prod\limits_{j=1}^{r-1}(2-{z_{j}^{n}-z_{j}^{-n}})=-\prod\limits_{j=1}^{r-1}({2T(n,w_{j})-2}).$

Substituting the latter in Theorem 3.1 we get the result.

##  4 Arithmetical properties of complexity for the graph $\mathcal{D}_{n}$

The main result of this section is the following theorem.

###### Theorem 4.1

Let $\tau(n)$ be the number of spanning trees in the graph

$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}}).$

Denote by $\beta_{odd}$ ($\gamma_{odd}$ resp.) the quantity of odd numbers in
the sequences $\beta_{1},\beta_{2},\ldots,\beta_{s}$
($\gamma_{1},\gamma_{2},\ldots,\gamma_{t}$ resp.). Also, denote by
$\gamma_{even}$ the quantity of even numbers in the sequence
$\gamma_{1},\gamma_{2},\ldots,\gamma_{t}.$ Let $\delta$ be the square free
part of the integer
$\xi=(2\beta_{odd}+\gamma_{odd})(2\beta_{odd}+\gamma_{even}).$ Then there
exists an integer sequence $a(n)$ such that

  1. $1^{\circ}$

$\tau(n)=n\,t\,a(n)^{2},$ if $n$ is odd;

  2. $2^{\circ}$

$\tau(n)=n\,t\,\delta\,a(n)^{2},$ if $n$ is even.

Proof: We consider the associated polynomial $P(z)$ for the graph
$\mathcal{D}_{n}.$ Now we will express the value $P(-1)$ through basic
parameters of $\mathcal{D}_{n}.$ Denote by $\beta_{even}$ ($\beta_{odd}$
resp.) the number of even numbers (odd numbers resp.) in the sequence
$\beta_{1},\beta_{2},\ldots,\beta_{s}.$ Denote by $\gamma_{even}$
($\gamma_{odd}$ resp.) the number of even numbers (odd numbers resp.) in the
sequence $\gamma_{1},\gamma_{2},\ldots,\gamma_{t}.$ One easily sees that
$s=\beta_{even}+\beta_{odd}$ and $t=\gamma_{even}+\gamma_{odd}.$ Since
$P(z)=\mathcal{A}(z)\mathcal{A}(z^{-1})-\mathcal{B}(z)\mathcal{B}(z^{-1}),$ we
get

$\displaystyle P(-1)$ $\displaystyle=$
$\displaystyle\mathcal{A}(-1)^{2}-\mathcal{B}(-1)^{2}=(2s+t-\sum_{j=1}^{s}((-1)^{\beta_{j}}+(-1)^{-\beta_{j}}))^{2}-(\sum_{j=1}^{t}(-1)^{-\gamma_{j}})^{2}$
$\displaystyle=$
$\displaystyle(2s+t-2\sum_{j=1}^{s}(-1)^{\beta_{j}})^{2}-(\sum_{j=1}^{t}(-1)^{-\gamma_{j}})^{2}$
$\displaystyle=$
$\displaystyle(2(\beta_{even}+\beta_{odd})+\gamma_{even}+\gamma_{odd}-2(\beta_{even}-\beta_{odd}))^{2}-(\gamma_{even}-\gamma_{odd})^{2}$
$\displaystyle=$ $\displaystyle
4(2\beta_{odd}+\gamma_{even})(2\beta_{odd}+\gamma_{odd}).$

As a consequence, we have $P(-1)=4\xi=\delta(2\omega)^{2},$ where $\delta$ is
the square free part of $\xi$ and $\omega$ is some integer.

By formula (2) we have
$n\,\tau(n)=t\prod_{j=1}^{n-1}\lambda_{j,1}\lambda_{j,2}.$ Note that
$\lambda_{j,1}\lambda_{j,2}=P(\varepsilon_{n}^{j})=P(\varepsilon_{n}^{n-j})=\lambda_{n-j,1}\lambda_{n-j,2}.$
Define $c(n)=\prod\limits_{j=1}^{\frac{n-1}{2}}\lambda_{j,1}\lambda_{j,2},$ if
$n$ is odd, and
$d(n)=\prod\limits_{j=1}^{\frac{n}{2}-1}\lambda_{j,1}\lambda_{j,2},$ if $n$ is
even. Following [11] we note that each algebraic number $\lambda_{i,j}$ comes
into both products $\prod_{j=1}^{(n-1)/2}\lambda_{j,1}\lambda_{j,2}$ and
$\prod_{j=1}^{n/2-1}\lambda_{j,1}\lambda_{j,2}$ with all of its Galois
conjugate elements. Therefore, both products $c(n)$ and $d(n)$ are integer
numbers. Moreover, if $n$ is even, then we get
$\lambda_{\frac{n}{2},1}\lambda_{\frac{n}{2},2}=P(-1).$ Now we have

  * $1^{\circ}$

$n\,\tau(n)=t\,c(n)^{2}$ if $n$ is odd,

  * $2^{\circ}$

$n\,\tau(n)=t\,P(-1)\,d(n)^{2}=4t\,\delta\,\omega^{2}\,d(n)^{2}$ if $n$ is
even.

Note that using formula (5) from the proof of Theorem 3.1 we conclude that
$\frac{\tau(n)}{n\,t}$ is an integer. Indeed, since $\widetilde{P}(z)$ and
$P(z)$ share the roots we have

$\tau(n)=(-1)^{r(n-1)}\eta^{n-1}n\,t\prod\limits_{\begin{subarray}{c}\widetilde{P}(z)=0\\\
z\neq 1\end{subarray}}\frac{z^{n}-1}{z-1}.$

The last product is equal to resultant of two integer polynomials
$\frac{\widetilde{P}(z)}{(z-1)^{2}}$ and $\frac{z^{n}-1}{z-1}$ and, hence, it
is an integer number. So that $\frac{\tau(n)}{n\,t}$ is also an integer. We
get

  * $1^{\circ}$

$\frac{\tau(n)}{n\,t}=(\frac{c(n)}{n})^{2}$ if $n$ is odd,

  * $2^{\circ}$

$\frac{\tau(n)}{n\,t}=\delta\,(\frac{2\omega\,d(n)}{n})^{2}$ if $n$ is even.

As $\frac{\tau(n)}{n\,t}$ is an integer and $\delta$ is square free, all
squared rational numbers in $1^{\circ}$ and $2^{\circ}$ are integers. We set
$a(n)=\frac{c(n)}{n}$ if $n$ is odd and $a(n)=\frac{2\omega\,d(n)}{n}$ if $n$
is even. This proves the theorem.

##  5 Asymptotic formulas for the number of spanning trees

In this section we find the asymptotics for the number of spanning trees in
the graph

$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}}).$

To do this we suppose that parameters
$\beta_{1},\ldots,\beta_{s},\gamma_{1},\ldots,\gamma_{t}$ are fixed, and the
inequalities $0<\beta_{1}<\beta_{2}<\cdots<\beta_{s}<\frac{n}{2}$ and
$0\leq\gamma_{1}<\gamma_{2}<\cdots<\gamma_{t}\leq{n-1}$ hold for all
sufficiently large values $n.$ We suppose also that graphs $\mathcal{D}_{n}$
are connected.

According to Lemma 1.1 we have $(n,d)=1,$ where $d=\gcd(\beta_{j},\,1\leq
j\leq s,\,\gamma_{j}-\gamma_{k},\,1\leq j<k\leq t).$ Since $n$ is an
arbitrarily large number and $d$ is fixed, we can choose $n$ to be a multiple
of $d.$ Then $d=(n,d)=1.$ So that the condition of Lemma 2.3 is satisfied.
This gives the following result.

###### Theorem 5.1

Let
$\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm\beta_{1}},b^{\pm\beta_{2}},\ldots,b^{\pm\beta_{s}},ab^{\gamma_{1}},ab^{\gamma_{2}},\ldots,ab^{\gamma_{t}})$
be an infinite family of connected graphs. Then the asymptotic behaviour for
the number of spanning trees $\tau(n)$ for the graph $\mathcal{D}_{n}$ is
given by the formula

$\tau(n)\sim\frac{n\,t}{q}A^{n},\,n\to\infty,$

where $A=\exp\left({\int\limits_{0}^{1}\log
P(e^{2\pi\texttt{i}t})\textrm{d}t}\right)$ and
$q=2t\sum_{j=1}^{s}\beta_{j}^{2}+\sum_{1\leq j<k\leq
t}(\gamma_{j}-\gamma_{k})^{2}.$

Proof: By Theorem 3.1 we have
$\tau(n)=\frac{n\,t\,|\eta|^{n}}{q}\prod\limits_{j=1}^{r-1}|{2T(n,w_{j})-2}|,$
where $w_{j},\,j=1,2,\ldots,r-1$ are roots, different from $1,$ of the
Chebyshev transform of $P(z).$

By Lemma 2.3, $T(n,w_{j})=\frac{z_{j}^{n}+z_{j}^{-n}}{2},$ where the $z_{j}$
and $1/z_{j}$ are roots of the polynomial $P(z)$ with the property
$|z_{j}|\neq 1,\,j=1,2,\ldots,r-1.$ Replacing $z_{j}$ by $1/z_{j},$ if it is
necessary, we can assume that $|z_{j}|>1$ for all $j=1,2,\ldots,r-1.$ Then
$T(n,w_{j})\sim\frac{1}{2}z_{j}^{n}$ and $|2T(n,w_{j})-2|\sim|z_{j}|^{n}$ as
$n\to\infty.$ Hence

$\frac{n\,t|\eta|^{n}}{q}\prod_{j=1}^{r-1}|2T_{n}(w_{j})-2|\sim\frac{n\,t|\eta|^{n}}{q}\prod_{j=1}^{r-1}|z_{j}|^{n}=\frac{n\,t}{q}|\eta|^{n}\prod\limits_{\begin{subarray}{c}P(z)=0,\\\
|z|>1\end{subarray}}|z|^{n}=\frac{n\,t\,A^{n}}{q},$

where $A=|\eta|\prod\limits_{P(z)=0,\,|z|>1}|z|$ is the Mahler measure of the
polynomial $P(z).$ By ([6], p. 67), we have
$A=\exp\left(\int_{0}^{1}\log|P(e^{2\pi it})|\textrm{d}t\right).$ The theorem
is proved.

##  6 Generating function for the number of spanning trees

In this section, our aim is to prove the following result.

###### Theorem 6.1

Let $\tau(n)$ be the number of spanning trees in the graph $\mathcal{D}_{n}.$
Then $F(x)=\sum\limits_{n=1}^{\infty}\tau(n)x^{n}$ is a rational function with
integer coefficients. Moreover, $F(\eta\,x)=F(\frac{1}{\eta\,x}),$ where
$\eta$ is the leading coefficient of the associated polynomial $P(z).$ The
latter allows to represent $F(x)$ as a rational function of
$u=\frac{1}{2}(\eta\,x+\frac{1}{\eta\,x}).$

The proof of Theorem 6.1 is based on the following proposition proved in [17].

###### Proposition 6.2

Let $R(z)$ be a degree $2s$ polynomial with integer coefficients. Suppose that
all the roots of the polynomial $R(z)$ are
$\xi_{1},\xi_{2},\ldots,\xi_{2s-1},\xi_{2s}.$ Then

$F(x)=\sum\limits_{n=1}^{\infty}\Big{(}n\prod\limits_{j=1}^{2s}(\xi_{j}^{n}-1)\Big{)}x^{n}$

is a rational function with integer coefficients.

Moreover, if $\xi_{j+s}=\xi_{j}^{-1},\,j=1,2,\ldots,s,$ then $F(x)=F(1/x).$

Proof of Theorem 6.1. By formula (6) we have

$F(x)=\sum\limits_{n=1}^{\infty}\tau(n)x^{n}=\sum\limits_{n=1}^{\infty}\Big{(}\frac{(-1)^{r(n-1)+1}n\,t\,\eta^{n}}{q}\prod\limits_{\begin{subarray}{c}P(z)=0\\\
z\neq 1\end{subarray}}(z^{n}-1)\Big{)}x^{n}.$

Since all the roots of $P(z),$ different from $1,$ are
$z_{1},1/z_{1},\ldots,z_{r-1},1/z_{r-1},$ we can rewrite the latter as

$F(x)=\frac{(-1)^{-r+1}t}{q}\sum\limits_{n=1}^{\infty}\Big{(}n\prod\limits_{j=1}^{r-1}(z_{j}^{n}-1)(z_{j}^{-n}-1)\Big{)}((-1)^{r}\eta\,x)^{n}.$

Since $t$ and $q$ are rational numbers, by Proposition 6.2, $F(x)$ is a
rational function with integer coefficients satisfying
$F((-1)^{r}\eta\,x)=F(\frac{1}{(-1)^{r}\eta\,x}).$ Hence,
$F(\eta\,x)=F(\frac{1}{\eta\,x}).$

##  7 Examples

###  7.1 Prism graph $\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm 1},a).$

$\textbf{1}^{\circ}.$ The number of spanning trees. The associated Laurent
polynomial and its Chebyshev transform are

$P(z)=z^{-2}-6z^{-1}+10-6z+z^{2}\text{ and }Q(w)=4(w-2)(w-1).$

Here $t=1,\,\eta=1,\,q=2.$ Hence, by Theorem 3.2 we have

$\tau(n)=\frac{n\,t\,\eta^{n}}{q}(2T(n,2)-2)=n(T(n,2)-1).$

This coincides with the well-known result in [4].

$\textbf{2}^{\circ}.$ The asymptotics of $\tau(n).$ By Theorem 5.1,
$\tau(n)\cong\frac{n}{2}A^{n},$ where $A=2+\sqrt{3}.$

$\textbf{3}^{\circ}.$ The generating function of $\tau(n).$ Theorem 6.1 gives

$F(x)=\sum\limits_{n=1}^{\infty}\tau(n)x^{n}=\frac{-3+u+u^{2}}{2(-2+u)^{2}(-1+u)},$

where $u=\frac{1}{2}(x+\frac{1}{x}).$

$\textbf{4}^{\circ}.$ Divisibility by squares. To see the divisibility by
squares, consider a few terms of generating function

$F(x)=x+12x^{2}+75x^{3}+384x^{4}+1805x^{5}+8100x^{6}+35287x^{7}+150528x^{8}+\ldots.$

By Theorem 4.1, we have
$\xi=(2\beta_{odd}+\gamma_{odd})(2\beta_{odd}+\gamma_{even})=(2\cdot 1)(2\cdot
1+1)=6.$ Hence, $\delta=6.$ So that there exists an integer sequence $a(n)$
such that $\tau(u)=n\,a(n)^{2}$ if $n$ is odd and $\tau(u)=6n\,a(n)^{2}$ if
$n$ is even.

###  7.2 Dihedral graph $\mathcal{D}_{n}=Cay(\mathbb{D}_{n},b^{\pm 1},b^{\pm
2},ab,ab^{3},ab^{5}).$

$\textbf{1}^{\circ}.$ The number of spanning trees. By Theorem 3.1

$\tau(n)=\frac{n\,t\,\eta^{n}}{q}|(-2)^{n}-1|\cdot|(-1/2)^{n}-1|\cdot|(4-\sqrt{15})^{n}-1|\cdot|(4+\sqrt{15})^{n}-1|,$

where $t=3,\,\eta=2,\,q=54.$ Equivalently, by Theorem 3.2 we get

$\tau(n)=\frac{n\,2^{n}}{18}|2T(n,-\frac{5}{4})-2|\cdot|2T(n,4)-2|.$

$\textbf{2}^{\circ}.$ The asymptotics of $\tau(n).$ By Theorem 5.1, we have
$\tau(n)\cong\frac{n}{18}A^{n},$ where $A=4(4+\sqrt{15}).$

$\textbf{3}^{\circ}.$ The generating function of $\tau(n).$ From Theorem 6.1,
we get

$F(x)=\sum\limits_{n=1}^{\infty}\tau(n)x^{n}$

$=\frac{6(-1745300+4540750u-3003815u^{2}+346990u^{3}+171265u^{4}-47660u^{5}+4840u^{6}-272u^{7}+16u^{8})}{(2+u)(8+u)^{2}(-5+2u)^{2}(265-80u+4u^{2})^{2}},$

where $u=\frac{1}{2}(2x+\frac{1}{2x}).$

$\textbf{4}^{\circ}.$ Divisibility by squares. To see the divisibility by
squares we consider a few terms of generating function

$F(x)=3x+60x^{2}+6561x^{3}+192000x^{4}+9149415x^{5}+315059220x^{6}+\cdots.$

By Theorem 4.1, we have
$\xi=(2\beta_{odd}+\gamma_{odd})(2\beta_{odd}+\gamma_{even})=(2\cdot
1+3)(2\cdot 1+0)=10.$ Hence, $\delta=10.$ So that there exists an integer
sequence $a(n)$ such that $\tau(u)=3n\,a(n)^{2}$ if $n$ is odd and
$\tau(u)=30n\,a(n)^{2}$ if $n$ is even.

## ACKNOWLEDGMENTS

The first author was supported by NSFC, no. 11831004 and Shanghai Science and
Technology Program [Project No. 22JC1400100]. The second and the third authors
were supported by Mathematical Center in Akademgorodok under agreement No.
075-15-2019-1613 with the Ministry of Science and Higher Education of the
Russian Federation. The fourth author is supported by NSFC, no.12101125,
Natural Science Foundation of Fujian Province [Project No.2021J05035].

## References

  * [1] Abrosimov N. V., Baigonakova G. A., Mednykh I. A., Counting spanning trees in cobordism of circulant graphs, Sib. Electron. Mat. Rep., 15 (2018), 1145–1157.
  * [2] Apostol T.M., Introduction to Analytic Number Theory, Springer-Verlag, New York, 1976.
  * [3] Babai L., Spectra of Cayley Graphs, J. Combin. Theory Series B.. 27 (1979), 180–189.
  * [4] Boesch F. T., Prodinger H., Spanning tree formulas and Chebyshev polynomials, Graphs and Combinatorics,. 2(1) (1986), 191–200.
  * [5] Davis P.J., Circulant Matrices, AMS Chelsea Publishing, 1994.
  * [6] Everest G., Ward T., Heights of polynomials and entropy in algebraic dynamics. Springer Science & Business Media, 2013.
  * [7] Grunwald L.A., Mednykh I.A., The number of rooted forests in circulant graphs, Ars Math. Contemp., 22(4) (2022), #10, 12 pp.
  * [8] Grunwald, L. A., Kwon, Y. S., Mednykh, I. A., Counting rooted spanning forests for circulant foliation over a graph, Tohoku Math. J. (2)., 74(4) (2022), 535–548.
  * [9] Jacobsen J.L., Salas J. and Sokal A.D., Spanning forests and the $q$-State Potts model in the limit $q\to 0,$ J. Stat. Phys., 119 (2005), 1153–1281.
  * [10] Kwon Y.S., Mednykh A.D., Mednykh I.A., Complexity of the circulant foliation over a graph, J. Algebraic Combin., 53(1) (2021), 115–129.
  * [11] Lorenzini D., Smith normal form and Laplacians, J. Combin. Theory Ser. B., 98(6) (2008), 1271–1300.
  * [12] Louise J., A formula for the number of spanning trees in circulant graphs with nonfixed generators and discrete tori, Bull. Australian Math. Soc., 92(3) (2015), 365–373.
  * [13] Mason J.C., Handscomb D.C., Chebyshev Polynomials, CRC Press, Boca Raton, 2003.
  * [14] Mednykh I. A., On Jacobian group and complexity of $I$-graph $I(n,k,l)$ through Chebyshev polynomials, Arc Math. Contemp., 15(2) (2018), 467–485.
  * [15] Mednykh, A. D., Mednykh, I. A., The number of spanning trees in circulant graphs, its arithmetic properties and asymptotic, Discrete Math., 342(6) (2019), 1772–1781.
  * [16] Mednykh, A. D., Mednykh, I. A., Kirchhoff index for circulant graphs and its asymptoticsm, Dokl. Math., 102(2) (2020), 392–395.
  * [17] Mednykh, A. D., Mednykh, I. A., On rationality of generating function for the number of spanning trees in circulant graphs, Algebra Colloq., 27(1) (2020), 87–94.
  * [18] Mednykh A.D., Mednykh I.A., Cyclic coverings of graphs. Counting rooted spanning forests and trees, Kirchhoff index, and Jacobians, Russian Math. Surveys 78:3 (2023), 501–548.
  * [19] Prasolov V. V., Polynomials, Series: Algorithms and Computation in Mathematics, 11, Springer, Berlin Heidelberg, 2004.
  * [20] Shrock R., Wu F.Y., Spanning trees on graphs and lattices in $d$-dimensions, J. Phys. A, 33 (2000), 3881–3902.
  * [21] Wu F.Y., Number of spanning trees on a lattice, J. Phys. A: Math. Gen., 10 (1977), L113–L115.
  * [22] Chen Xiebin, Qiuying Lin, Fuji Zhang, The number of spanning trees in odd valent circulant graphs, Discrete Math., 282(1) (2004), 69–79.
  * [23] Zhang Yuanping, Yong Xuerong, Golin M. J., The number of spanning trees in circulant graphs, Discrete. Math., 223(1) (2000), 337–350.
  * [24] Zhang Yuanping, Xuerong Yong, Golin M. J., Chebyshev polynomials and spanning tree formulas for circulant and related graphs, Discrete Math., 298(1) (2005), 334–364.


# Constraining cosmological phase transitions with the Parkes Pulsar Timing
Array

Xiao Xue  CAS Key Laboratory of Theoretical Physics, Institute of Theoretical
Physics, Chinese Academy of Sciences, Beijing 100190, China  School of
Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049,
China  II. Institute of Theoretical Physics, Universität Hamburg, 22761
Hamburg, Germany  Ligong Bian<EMAIL_ADDRESS>Department of Physics,
Chongqing University, Chongqing 401331, China  Chongqing Key Laboratory for
Strongly Coupled Physics, Chongqing 401331, China  Jing Shu
<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics, Insitute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  CAS Center for Excellence in Particle Physics, Beijing
100049, China  School of Fundamental Physics and Mathematical Sciences,
Hangzhou Institute for Advanced Study, University of Chinese Academy of
Sciences, Hangzhou 310024, China  Center for High Energy Physics, Peking
University, Beijing 100871, China  International Center for Theoretical
Physics Asia-Pacific, Beijing/Hangzhou, China  Qiang Yuan<EMAIL_ADDRESS>Key
Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory,
Chinese Academy of Sciences, Nanjing 210023, China  School of Astronomy and
Space Science, University of Science and Technology of China, Hefei 230026,
China  Center for High Energy Physics, Peking University, Beijing 100871,
China  Xingjiang Zhu<EMAIL_ADDRESS>School of Physics and Astronomy, Monash
University, Clayton, VIC 3800, Australia  OzGrav: The ARC Centre of Excellence
for Gravitational Wave Discovery, Hawthorn, VIC 3122, Australia  Advanced
Institute of Natural Sciences, Beijing Normal University at Zhuhai 519087,
China  N. D. Ramesh Bhat  International Centre for Radio Astronomy Research,
Curtin University, Bentley, WA 6102, Australia  Shi Dai  Western Sydney
University, Locked Bag 1797, Penrith South DC, NSW 1797, Australia  Yi Feng
University of Chinese Academy of Sciences, Beijing 100049, China  Boris
Goncharov  School of Physics and Astronomy, Monash University, Clayton, VIC
3800, Australia  OzGrav: The ARC Centre of Excellence for Gravitational Wave
Discovery, Hawthorn, VIC 3122, Australia  George Hobbs  CSIRO Astronomy and
Space Science, P.O. Box 76, Epping, NSW 1710, Australia  Eric Howard  CSIRO
Astronomy and Space Science, P.O. Box 76, Epping, NSW 1710, Australia
Macquarie University, Department of Physics and Astronomy, Sydney, NSW, 2109,
Australia  Richard N. Manchester  CSIRO Astronomy and Space Science, P.O. Box
76, Epping, NSW 1710, Australia  Christopher J. Russell  CSIRO Scientific
Computing, Australian Technology Park, Locked Bag 9013, Alexandria, NSW 1435,
Australia  Daniel J. Reardon  OzGrav: The ARC Centre of Excellence for
Gravitational Wave Discovery, Hawthorn, VIC 3122, Australia  Centre for
Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box
218, Hawthorn, VIC 3122, Australia  R. M. Shannon  OzGrav: The ARC Centre of
Excellence for Gravitational Wave Discovery, Hawthorn, VIC 3122, Australia
Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia  Renée Spiewak
Jodrell Bank Centre for Astrophysics, University of Manchester, Manchester M13
9PL, UK  Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia  Nithyanandan
Thyagarajan  CSIRO Astronomy and Space Science (CASS), P. O. Box 1130,
Bentley, WA 6102, Australia  Jingbo Wang  Xinjiang Astronomical Observatory,
Chinese Academy of Sciences, 150 Science 1-Street, Urumqi, Xinjiang 830011,
China

###### Abstract

A cosmological first-order phase transition is expected to produce a
stochastic gravitational wave background. If the phase transition temperature
is on the MeV scale, the power spectrum of the induced stochastic
gravitational waves peaks around nanohertz frequencies, and can thus be probed
with high-precision pulsar timing observations. We search for such a
stochastic gravitational wave background with the latest data set of the
Parkes Pulsar Timing Array. We find no evidence for a Hellings-Downs spatial
correlation as expected for a stochastic gravitational wave background.
Therefore, we present constraints on first-order phase transition model
parameters. Our analysis shows that pulsar timing is particularly sensitive to
the low-temperature ($T\sim 1-100$ MeV) phase transition with a duration
$(\beta/H_{*})^{-1}\sim 10^{-2}-10^{-1}$ and therefore can be used to
constrain the dark and QCD phase transitions.

##  I introduction

The stochastic gravitational wave background (SGWB) is an important target of
gravitational-wave astronomy after the successful detection of gravitational
waves from 50 compact binary mergers by the LIGO and Virgo detectors
abbott2016observation ; abbott2019gwtc ; GWTC-2 . The SGWB may come from the
superposition of numerous unresolvable binary coalescences, or from the
cosmological inflation, cosmological phase transition, or topological defects
expected to form during cosmological phase transitions, such as cosmic
strings. Since the discovery of gravitational waves in 2015, the study of SGWB
has gained increasingly broad interest because of its close connection with
particle physics and cosmology. Different gravitational wave experiments, such
as LIGO and the upcoming Laser Interferometer Space Antenna (LISA), aim at the
detection of SGWB at different frequency bands Aasi:2013wya ; Audley:2017drz ;
Caprini:2019egz .

The pulsar timing array (PTA) experiments measure the arrival times of radio
pulses from a number of highly stable millisecond pulsars in the Milky Way,
which enables the detection of gravitational waves of nanohertz frequencies.
The three major PTA experiments, the Parkes PTA (PPTA; manchester2013parkes ),
the European PTA EPTA16 , and NANOGrav NANOGrav12.5data  have been monitoring
the times of arrival (ToAs) of dozens of pulsars on weekly to monthly cadence
for over 10 years, and provide unique probes of the SGWB in the nanohertz band
Lentati:2015qwp ; Caprini:2018mtu ; Arzoumanian:2018saf ; Garcia-
Bellido:2017aan ; Barack:2018yly ; Caprini:2010xv ; Breitbach:2018ddu .

Whereas the observation of the 125 GeV Higgs boson at the Large Hadron
Collider indicates that the phase transition in the Standard Model of particle
physics is cross-over DOnofrio:2014rug , a first-order phase transition (FOPT)
is predicted in many new physics models beyond the Standard Model, such as the
real and complex extended standard models Espinosa:2011ax ; Profumo:2014opa ;
Jiang:2015cwa 111See Ref. Profumo:2007wc  for the early study on basic
features of the Electroweak phase transition in the singlet model and its
collider phenomenology. , two-Higgs doublet models Andersen:2017ika ;
Bernon:2017jgv ; Dorsch:2013wja ; Basler:2016obg , next-to-minimal
supersymmetric model Kozaczuk:2014kva ; Huber:2015znp ; Bian:2017wfv ,
composite Higgs models Espinosa:2011eu ; Bian:2019kmg , or the hidden sector
Schwaller:2015tja . If the FOPT occurs at a relatively low temperature, e.g.
around $T_{*}\sim$ MeV, a SGWB with a peak frequency $f\sim 10^{-9}-10^{-8}$
Hz is expected Caprini:2010xv ; Bian:2020urb , which falls into the sensitive
region of PTA experiments foster1990constructing ; hellings1983upper . In this
work, we perform a search for the SGWB from the first-order phase transition
(FOPT) of the early Universe using the latest PPTA data set. We show that the
PPTA data can probe a considerable proportion of the parameter space of the
FOPT. We ignore contribution from other sources to the SGWB; thus, constraints
on FOPT derived in this work can be considered as being conservative.

##  II Gravitational waves from FOPT

There are three main sources of gravitational waves produced during the FOPT:
bubble collisions Caprini:2007xq ; Huber:2008hg , sound waves
Hindmarsh:2017gnf ; Hindmarsh:2013xza ; Hindmarsh:2016lnk ; Cutting:2019zws ;
hindmarsh2015numerical , and MHD turbulence Caprini:2009yp . It is believed
that sound waves are dominate in many phase transition models. In this work we
only consider the sound wave contribution to the energy density of the SGWB
from the FOPT, for the sake of conservativeness. Sound waves contribution to
the SGWB over a characteristic frequency band $f$ can be expressed as
Caprini:2015zlo ; Hindmarsh:2017gnf ; Ellis:2018mja ; Ellis:2019oqb

$\displaystyle\Omega_{sw}(f)h^{2}$ $\displaystyle=2.65\times
10^{-6}v_{w}\left(\frac{H_{*}}{\beta}\right)\left(\frac{\kappa\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}$
$\displaystyle\times\left[(f/f_{sw})^{3}\left(\frac{7}{4+3(f/f_{sw})^{2}}\right)^{7/2}\right]$
$\displaystyle\times\min\left[1,~{}(8\pi)^{1/3}v_{w}\left(\frac{H_{*}}{\beta}\right)\left(\frac{3}{4}\frac{\kappa\alpha}{1+\alpha}\right)^{-1/2}\right],$
(1)

where $g_{*}$ is the effective number of relativistic degrees of freedom;
$\kappa$ is the fraction of the latent heat of phase transition transferred
into the kinetic energy of plasma Espinosa:2010hh ; $v_{w}$ is the bubble wall
velocity(we take $v_{w}=1$ in this work); $H_{*}$ is the Hubble parameter at
the phase transition temperature $T_{*}$; $\alpha$ is phase transition
strength, which largely determine the amplitude of the GW spectrum; $\beta$ is
the inverse duration of the phase transition (or the phase transition rate
parameter), which together with $T_{*}$ set the peak frequency of the GW
spectrum. $f_{sw}$ is the peak frequency of the sound wave defined as

$\displaystyle f_{sw}$ $\displaystyle\simeq\frac{1.9\times
10^{-5}~{}\mathrm{Hz}}{v_{w}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100\textrm{
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{1/6}\;.$ (2)

We fix the value of $g_{*}$ at different $T_{*}$: $g_{*}=100$ when $T_{*}>0.2$
GeV, $g_{*}=10$ when $0.1$ MeV $<T_{*}<0.2$ GeV, and $g_{*}=2$ when
$T_{*}<0.1$ MeV husdal2016effective . The set of parameters $\alpha$ and
$\beta$ are defined at the nucleation temperature and entirely determine the
gravitational wave spectrum. We note that, despite $g_{*}$ changing
drastically at the critical temperatures, the accuracy of our model is not
significantly affected due to the small power in Eqs. (1) and (2).

##  III Data analysis

The PTA experiments record pulse ToAs from highly stable millisecond pulsars.
The SGWB can be detected by searching for specific statistical correlations
among timing residuals of different pulsars allen1999detecting ;
maggiore2000gravitational . The noise model we use is based on a detailed
single pulsar analysis goncharov2020identifying , and it is publicly available
from Github repositories
ppta_dr2_noise_analysis222https://github.com/bvgoncharov/ppta_dr2_noise_analysis
and enterprise_warp333https://github.com/bvgoncharov/enterprise_warp. We use
the TEMPO2 software hobbs2006tempo2 ; edwards2006tempo2  to calculate the
timing residuals, and enterprise444https://github.com/nanograv/enterprise,
PTMCMC555https://github.com/jellis18/PTMCMCSampler,
enterprise_extensions666https://github.com/nanograv/enterprise_extensions for
the Bayesian analysis. We also use the Jet Propulsion Laboratory solar-system
ephemeris DE436, and the TT (BIPM18) reference timescale published by the
International Bureau of Weights and Measures (BIPM).

####  III.0.1 Noise Model

We briefly introduce our noise models here; we refer interested readers to
Ref. goncharov2020identifying  for more details. The noise components can be
categorized as white noise and red noise. We use white noise parameters EFAC
and EQUAD to account for additional white noise that is not included in ToA
uncertainties for each signal processor/receiver system van2011placing ;
liu2011prospects ; liu2012profile . For some pulsars with multi-channel ToAs,
white noise parameter ECORR is also included to account for the pulse phase
jitter noise arzoumanian2015nanograv . For red noise, there are power-law spin
noise, dispersion measure (DM) noise, chromatic red noise, and band noise
van2013understanding ; van2015low ; caballero2016noise ; lentati2016spin ;
you2007dispersion ; keith2013measurement . The exact composition of red noise
terms differs from one pulsar to another. The formalism of the pulsar
intrinsic red noises can be found in the Supplemental Material. The
“exponential dip” events are included for 4 pulsars. For each pulsar, we apply
a single pulsar noise analysis first. Then in the global Bayesian analysis, we
fix the white noise parameters at their best-fit values from the single pulsar
analyses, and vary the red noise and Bayes ephemeris parameters. The priors of
the noise parameters can be found in Ref. goncharov2020identifying .

Recently, the NANOGrav group found strong statistical support for a common
power-law (CPL) spectrum process in their 12.5-year data set
Arzoumanian:2020vkk . Following the same analysis procedure/assumptions, such
a CPL process was also identified in the PPTA data Boris21PPTAcrn . However,
it was pointed out that similar noise processes that are intrinsic to
individual pulsars might lead to a false detection of a common noise. In both
studies, evidence for the Hellings-Downs (HD) spatial correlation as expected
from a SGWB is insignificant, and the nature of this CPL process is yet to be
determined. We refer interested readers to Ref. Boris21PPTAcrn  for discussion
on the possible common noise process. Following these two analyses, we
introduce a common-spectrum process that is described by a power-law spectrum
with the same amplitude and power-law index across different pulsars. Note
that this work searches for SGWB from the FOPT process, which differs from
Ref. Boris21PPTAcrn  both physically and technically.

####  III.0.2 Correlation search for gravitational waves

The energy density of a SGWB, $\Omega_{gw}h^{2}$, is related with the one-
sided power spectral density, $S_{gw}(f)$, as

$S_{gw}(f)=\frac{1}{12\pi^{2}}\frac{1}{f^{5}}\frac{3H_{100}^{2}}{2\pi^{2}}\Omega_{gw}(f)h^{2},$
(3)

where $H_{100}=100{~{}\rm km~{}s^{-1}~{}Mpc^{-1}}$ and $h\approx 0.7$. The
statistical correlation between timing residuals $\delta t$ of pulsar $I$ and
pulsar $J$ is

$\langle\delta t_{i}^{I}~{}\delta
t_{j}^{J}\rangle=\int\mathrm{d}f~{}S(f)~{}\Gamma^{IJ}(\theta_{IJ})~{}\cos\left[2\pi
f\left(t_{i}-t_{j}\right)\right],$ (4)

where $t_{i}$ and $t_{j}$ are ToAs, $\Gamma^{IJ}(\theta_{IJ})$ is the HD
correlation, and $\theta_{IJ}$ is the angle between two pulsars (when $I=J$,
$\theta_{IJ}=0$, $\Gamma^{IJ}=1$). Eq. (4) is valid under the assumption that
the gravitational wave is isotropic and unpolarized, and only the Earth terms
are considered. The time-frequency method is adopted, where a finite number of
Fourier frequencies are chosen to approximate the power spectrum for an
acceptable computation efficiency. For the gravitational wave spectrum, we
also consider the case when the pulsar auto-correlation is subtracted. Namely,
we use an off-diagonal HD correlation $\tilde{\Gamma}^{IJ}(\theta_{IJ})$, and
when $I=J$, we have $\theta_{IJ}=0$, $\tilde{\Gamma}^{IJ}_{0}=0$. The
correlation between a pulsar’s own ToAs is referred to as pulsar auto-
correlation. We use terms “full HD” and “no-auto HD” to differentiate these
two cases.

####  III.0.3 SGWB spectrum

We first consider that the power spectrum of the SGWB is a free function at 20
logarithmically spaced frequency modes between $10^{-9}$ and $10^{-8}$ Hz.
Therefore the integral in Eq. (4) is calculated as a discrete summation of 20
elements van2015low . We focus on the low-frequency range because PTAs are
more sensitive in this range. We treat $\sqrt{\Omega_{gw}(f_{i})h^{2}}$ at
each $f_{i}$ as an independent parameter, for which we adopt a uniform prior
between $10^{-10}$ and $1$. The red noise parameters and the Bayes ephemeris
parameters are also treated as free parameters in the Bayesian analysis. The
Markov Chain Monte Carlo (MCMC) method is employed to compute the posterior
distribution of the parameters. The prior ranges of model parameters are
summarized in Table S1 of the Supplemental Material.

The FOPT model is roughly a broken power law with a peak frequency given by
Eq. (2). Similar to the free spectrum case, we also use 20 logarithmically
spaced frequency modes between $10^{-9}$ and $10^{-8}$ Hz in Eq. (4), but with
$\Omega_{gw}(f_{i})h^{2}$ as predicted by parameters $T_{*}$, $\alpha$, and
$(\beta/H_{*})^{-1}$.

The Savage-Dickey formula dickey1971weighted  is used to calculate the Bayes
factor, $\rm{BF}$, for nested models; in our case we wish to calculate the BF
between the signal hypothesis and the noise hypothesis:

$\mathrm{BF}\equiv\frac{\mathcal{Z}_{1}}{\mathcal{Z}_{0}}=\frac{\mathrm{P}(\varphi=\varphi_{0})}{\mathrm{P}(\varphi=\varphi_{0}|\boldsymbol{\rm
D})},$ (5)

where $\mathcal{Z}_{1}$ ($\mathcal{Z}_{0}$) is the model evidence for the
signal (noise) hypothesis, $\boldsymbol{\rm D}$ and $\varphi$ denote the data
and the signal parameters, respectively, and $\varphi_{0}$ represents the
noise-only model that is nested in the signal model (i.e., the signal has zero
amplitude). $\rm{P}(\varphi=\varphi_{0}|\boldsymbol{\rm D})$ and
$\rm{P}(\varphi=\varphi_{0})$ are the posterior and prior probability of
noise-only hypothesis, respectively. A nested model requires signal parameters
$\varphi$ to be independent of noise parameters $\vartheta$, which is true for
our FOPT model.

##  IV Results

Figure 1: The 95% credibility Bayesian upper limits on the energy density of a free-spectrum SGWB based on the PPTA data. The solid, dotted and dashed black curves indicate results under different assumptions (hypotheses H2, H3 and H4 in Table 1, respectively, where the FOPT spectrum is substituted with the free spectrum). The region excluded by BBN cyburt2005new  is shaded grey. Red curves with filled markers show examples of the spectrum predicted in the FOPT model, with $(\beta/H_{*})^{-1}=0.1$, and varying parameters $(\alpha,\,T_{*})=\\{(0.2,\,3~{}\textrm{MeV})\ \textrm{(circle)},\,(0.5,3~{}\textrm{MeV})\ \textrm{(triangle)},(0.5,\,30~{}\textrm{MeV})\ \textrm{(square)}\\}$, respectively. Table 1: Description of hypotheses tested in this work and the Bayes factors between them. Hypothesis |  Pulsar noise |  CPL process |  HD process FOPT spectrum | Bayes Factors | Parameter Estimation (median and 1-$\sigma$ interval)
---|---|---|---|---|---
$T_{*}$/MeV, $\alpha\times 10^{3}$, $\beta/H_{*}$ |  $A_{\rm CPL}$, $\gamma_{\rm CPL}$
H0:Pulsar Noise | ✓ |  |  |  |  |
H1:CPL | ✓ | ✓ |  |  $10^{2.9}$ (against H0) |  |  $-14.45_{-0.64}^{+0.62}$,$3.31_{-1.53}^{+1.36}$
H2:FOPT | ✓ |  | ✓(full HD) |  $10^{1.8}$ (against H0) |  $7.4_{-4.7}^{+11.9}$, $271_{-92}^{+165}$, $9.9^{+11.4}_{-5.4}$ |
H3:FOPT1 | ✓ | ✓ | ✓(full HD) | 1.04 (against H1) |  $9.6_{-9.2}^{+232.2}$, $3.8_{-3.4}^{+27.9}$, $854^{+9622}_{-782}$ |  $-14.51_{-0.68}^{+0.64}$,$3.36_{-1.54}^{+1.39}$
H4:FOPT2 | ✓ | ✓ | ✓(no-auto HD) | 0.96 (against H1) |  $10.9_{-10.6}^{+290.5}$, $3.2_{-2.8}^{+19.9}$, $1053^{+11256}_{-962}$ |  $-14.45_{-0.64}^{+0.62}$,$3.27_{-1.54}^{+1.37}$

We show in Fig. 1 the 95% credibility upper limits on the energy density of a
SGWB in the $10^{-9}-10^{-8}$ Hz frequency range. The SGWB is characterised
with a free spectrum. We note that the upper limits are almost unaffected by
the inclusion of a CPL process. This is due to the choice of priors for the
free-spectrum SGWB amplitude (uniform) and the CPL amplitude (log-uniform);
see Table S1 in the Supplemental Material for details. We choose a uniform
amplitude prior for the SGWB because it results in more conservative upper
limits than the log-uniform prior. The upper limits are slightly more
constraining for the case with only no-auto HD correlation included. Our
results improve significantly upon the Big Bang nucleosynthesis (BBN
cyburt2005new ) limit by two orders of magnitude in the frequency range
considered in this work.

In the search for a FOPT signal, we consider three cases, listed as H2, H3 and
H4 in Table 1, respectively. In H2, we assume that there is no CPL process. In
both H3 and H4, a CPL process is included. In H3, we include the full HD
correlation induced by a SGWB, whereas in H4 only the off-diagonal correlation
(denoted as no-auto HD) is used. For comparison, we also consider H1 where we
assume there is only a CPL process.

We summarize the hypotheses tested in this work and the Bayes factors between
them in Table 1. First, being consistent with a separate analysis
Boris21PPTAcrn , we find a strong evidence ($\rm{BF}\sim 10^{2.9}$) for a CPL
process, if we follow the same analysis procedure/assumptions in Ref
Arzoumanian:2020vkk . We then search for evidence of phase-transition
gravitational waves. Without adding a CPL term (H2), we obtain a BF of
$10^{1.8}$ in favour of the FOPT signal, which is less significant than that
of CPL process (H1). Nevertheless, we show the posterior distributions of the
FOPT model parameters in Fig. S1 of the Supplemental Material.

Since the origin of the apparent CPL process is still an open question, we
take a conservative approach by assuming that there are unknown systematic
errors that can be approximated by a CPL process777In the case that the
apparent CPL process is due to a supermassive binary black hole (SMBBH)
background, we find that the constraints on FOPT parameters are essentially
unchanged; see the Supplemental Material.. Under hypotheses H3 and H4, the
evidence of the FOPT signal is negligible ($\rm{BF}\sim 1$).

In Fig. 2, we show the exclusion contours (at the 95% credible level) on FOPT
model parameters for $\alpha=$ 0.2, 0.5, and 1. There are a few noteworthy
features. First, the constraints become tighter for a larger value of $\alpha$
(i.e., larger phase transition strength). Second, the excluded range of phase-
transition temperature $T_{*}$ is greater for longer phase-transition duration
(i.e., larger $(\beta/H_{*})^{-1}$). Third, the excluded parameter space is
larger for the case of no-auto HD correlation. This is because the FOPT models
is slightly disfavoured by the data when auto-correlation of the SGWB is
subtracted (see Table 1). Compared with the BBN bound on the phase transition
temperature ($T_{*}>1$ MeV), our results are more sensitive when $1\ {\rm
MeV}<T_{*}<100\ {\rm MeV}$.

Figure 2: The PPTA exclusion contours (at 95% C.L.) on FOPT parameters
$\log_{10}(T_{*}/{\rm GeV})$ and $\log_{10}[(\beta/H_{*})^{-1}]$, for
$\alpha=0.2$, 0.5, and 1.0.

##  V Conclusion

Pulsar timing arrays provide a unique way to detect gravitational waves in the
nanohertz frequency range ($10^{-9}$ Hz to $10^{-7}$ Hz), enabling the search
for an SGWB potentially produced from cosmological phase transitions at low
phase-transition temperatures ($\sim 1$ MeV). In the PPTA second Data Release
kerr2020parkes , we search for the SGWB induced by first-order phase
transitions of the early Universe. Assuming there is an unknown systematic
error that can be approximated by a CPL process, we find no evidence for a
Hellings-Downs spatial correlation induced by the phase-transition
gravitational wave signal. We are able to exclude a significant parameter
space of cosmological phase transitions that is not probed by previous
observations. We find that phase transitions can be effectively constrained by
the PPTA data for phase-transition temperatures of $\sim 1-100$ MeV and
durations $(\beta/H_{*})^{-1}\sim 10^{-2}-10^{-1}$, which are within
predictions from low-scale dark phase-transition models Ratzinger:2020koh ;
Bian:2020urb ; Breitbach:2018ddu  and cosmological first-order QCD phase
transitions Caprini:2010xv . The results therefore place constraints on the
new physics admitting a first-order phase transition at low-energy scales in
the early Universe, and provide an alternative way to probe particle physics
at low-energy scales.

Note added.—Recently the NANOGrav collaboration also performed a search for
the gravitational-wave background from the FOPT in their 12.5-year data set
NANOGrav21fopt . In that work, the gravitational-wave signal from phase
transitions was proposed as a possible interpretation of the detected CPL
process. As pointed out in Ref. Boris21PPTAcrn , the detection of a CPL
process might suffer from model misspecification. Therefore, we consider the
feature found in our data to be a noise component of unknown origin and derive
constraints on the FOPT parameters. Better understanding of the nature of the
CPL process requires longer data sets and more pulsars which might be
available through the International Pulsar Timing Array IPTAdr2  in the
future.

##  VI acknowledgements

We thank the referees for very useful comments. This work has been carried out
by the Parkes Pulsar Timing Array, which is part of the International Pulsar
Timing Array. The Parkes radio telescope (“Murriyang”) is part of the
Australia Telescope, which is funded by the Commonwealth Government for
operation as a National Facility managed by CSIRO. This work is supported in
part by the National Key Research and Development Program of China under Grant
No. 2020YFC2201501. XX is supported by Deutsche Forschungsgemeinschaft under
Germany’s Excellence Strategy EXC2121 “Quantum Universe” - 390833306. LB is
supported by the National Natural Science Foundation of China (NSFC) under the
grants No. 12075041 and No. 12047564, the Fundamental Research Funds for the
Central Universities of China (No. 2021CDJQY-011 and No. 2020CDJQY-Z003), and
Chongqing Natural Science Foundation (Grant No. cstc2020jcyj-msxmX0814). JS is
supported by the NSFC under Grants No. 12025507, No. 11690022, No. 11947302,
by the Strategic Priority Research Program and Key Research Program of
Frontier Science of the Chinese Academy of Sciences (CAS) under Grants No.
XDB21010200, No. XDB23010000, No. XDPB15, No. ZDBS-LY-7003, and by the CAS
Project for Young Scientists in Basic Research under Grant No. YSBR-006. QY is
supported by CAS and the Program for Innovative Talents and Entrepreneur in
Jiangsu. XZ was supported by ARC CE170100004. SD is the recipient of an
Australian Research Council Discovery Early Career Award, project DE210101738
funded by the Australian Government. RMS is the recipient of Australian
Research Council Future Fellowship FT190100155.

## References

  * (1) B. P. Abbott et al., “Observation of gravitational waves from a binary black hole merger,” Physical review letters, vol. 116, no. 6, p. 061102, 2016.
  * (2) B. Abbott et al., “Gwtc-1: a gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs,” Physical Review X, vol. 9, no. 3, p. 031040, 2019.
  * (3) R. Abbott et al., “GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run,” Physical Review X, vol. 11, p. 021053, Apr. 2021.
  * (4) B. Abbott et al., “Prospects for Observing and Localizing Gravitational-Wave Transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Rev. Rel., vol. 21, no. 1, p. 3, 2018.
  * (5) P. Amaro-Seoane et al., “Laser Interferometer Space Antenna,” arxiv:1702.00786, 2017.
  * (6) C. Caprini et al., “Detecting gravitational waves from cosmological phase transitions with LISA: an update,” JCAP, vol. 03, p. 024, 2020.
  * (7) R. N. Manchester et al., “The parkes pulsar timing array project,” Publications of the Astronomical Society of Australia, vol. 30, 2013.
  * (8) G. Desvignes et al., “High-precision timing of 42 millisecond pulsars with the European Pulsar Timing Array,” Mon. Not. Roy. Astron. Soc. , vol. 458, pp. 3341–3380, May 2016.
  * (9) M. F. Alam et al., “The NANOGrav 12.5 yr Data Set: Observations and Narrowband Timing of 47 Millisecond Pulsars,” Astrophys. J. Supp. , vol. 252, p. 4, Jan. 2021.
  * (10) L. Lentati et al., “European Pulsar Timing Array Limits On An Isotropic Stochastic Gravitational-Wave Background,” Mon. Not. Roy. Astron. Soc., vol. 453, no. 3, pp. 2576–2598, 2015.
  * (11) C. Caprini and D. G. Figueroa, “Cosmological Backgrounds of Gravitational Waves,” Class. Quant. Grav., vol. 35, no. 16, p. 163001, 2018.
  * (12) Z. Arzoumanian et al., “The NANOGrav 11-year Data Set: Pulsar-timing Constraints On The Stochastic Gravitational-wave Background,” Astrophys. J., vol. 859, no. 1, p. 47, 2018.
  * (13) J. Garcia-Bellido, M. Peloso, and C. Unal, “Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter,” JCAP, vol. 09, p. 013, 2017.
  * (14) L. Barack et al., “Black holes, gravitational waves and fundamental physics: a roadmap,” Class. Quant. Grav., vol. 36, no. 14, p. 143001, 2019\.
  * (15) C. Caprini, R. Durrer, and X. Siemens, “Detection of gravitational waves from the QCD phase transition with pulsar timing arrays,” Phys. Rev. D, vol. 82, p. 063511, 2010.
  * (16) M. Breitbach, J. Kopp, E. Madge, T. Opferkuch, and P. Schwaller, “Dark, Cold, and Noisy: Constraining Secluded Hidden Sectors with Gravitational Waves,” JCAP, vol. 07, p. 007, 2019.
  * (17) M. D’Onofrio, K. Rummukainen, and A. Tranberg, “Sphaleron Rate in the Minimal Standard Model,” Phys. Rev. Lett., vol. 113, no. 14, p. 141602, 2014.
  * (18) J. R. Espinosa, T. Konstandin, and F. Riva, “Strong Electroweak Phase Transitions in the Standard Model with a Singlet,” Nucl. Phys. B, vol. 854, pp. 592–630, 2012.
  * (19) S. Profumo, M. J. Ramsey-Musolf, C. L. Wainwright, and P. Winslow, “Singlet-catalyzed electroweak phase transitions and precision Higgs boson studies,” Phys. Rev. D, vol. 91, no. 3, p. 035018, 2015.
  * (20) M. Jiang, L. Bian, W. Huang, and J. Shu, “Impact of a complex singlet: Electroweak baryogenesis and dark matter,” Phys. Rev. D, vol. 93, no. 6, p. 065032, 2016.
  * (21) S. Profumo, M. J. Ramsey-Musolf, and G. Shaughnessy, “Singlet Higgs phenomenology and the electroweak phase transition,” JHEP, vol. 08, p. 010, 2007.
  * (22) J. O. Andersen, T. Gorda, A. Helset, L. Niemi, T. V. I. Tenkanen, A. Tranberg, A. Vuorinen, and D. J. Weir, “Nonperturbative Analysis of the Electroweak Phase Transition in the Two Higgs Doublet Model,” Phys. Rev. Lett., vol. 121, no. 19, p. 191802, 2018.
  * (23) J. Bernon, L. Bian, and Y. Jiang, “A new insight into the phase transition in the early Universe with two Higgs doublets,” JHEP, vol. 05, p. 151, 2018\.
  * (24) G. Dorsch, S. Huber, and J. No, “A strong electroweak phase transition in the 2HDM after LHC8,” JHEP, vol. 10, p. 029, 2013.
  * (25) P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt, and A. Wlotzka, “Strong First Order Electroweak Phase Transition in the CP-Conserving 2HDM Revisited,” JHEP, vol. 02, p. 121, 2017.
  * (26) J. Kozaczuk, S. Profumo, L. S. Haskins, and C. L. Wainwright, “Cosmological Phase Transitions and their Properties in the NMSSM,” JHEP, vol. 01, p. 144, 2015.
  * (27) S. J. Huber, T. Konstandin, G. Nardini, and I. Rues, “Detectable Gravitational Waves from Very Strong Phase Transitions in the General NMSSM,” JCAP, vol. 03, p. 036, 2016.
  * (28) L. Bian, H.-K. Guo, and J. Shu, “Gravitational Waves, baryon asymmetry of the universe and electric dipole moment in the CP-violating NMSSM,” Chin. Phys. C, vol. 42, no. 9, p. 093106, 2018.  [Erratum: Chin.Phys.C 43, 129101 (2019)].
  * (29) J. R. Espinosa, B. Gripaios, T. Konstandin, and F. Riva, “Electroweak Baryogenesis in Non-minimal Composite Higgs Models,” JCAP, vol. 01, p. 012, 2012.
  * (30) L. Bian, Y. Wu, and K.-P. Xie, “Electroweak phase transition with composite Higgs models: calculability, gravitational waves and collider searches,” JHEP, vol. 12, p. 028, 2019.
  * (31) P. Schwaller, “Gravitational Waves from a Dark Phase Transition,” Phys. Rev. Lett., vol. 115, no. 18, p. 181101, 2015.
  * (32) L. Bian, R.-G. Cai, J. Liu, X.-Y. Yang, and R. Zhou, “Evidence for different gravitational-wave sources in the NANOGrav dataset,” Phys. Rev. D, vol. 103, no. 8, p. L081301, 2021.
  * (33) R. S. Foster and D. Backer, “Constructing a pulsar timing array,” The Astrophysical Journal, vol. 361, pp. 300–308, 1990.
  * (34) R. Hellings and G. Downs, “Upper limits on the isotropic gravitational radiation background from pulsar timing analysis,” The Astrophysical Journal, vol. 265, pp. L39–L42, 1983.
  * (35) C. Caprini, R. Durrer, and G. Servant, “Gravitational wave generation from bubble collisions in first-order phase transitions: An analytic approach,” Phys. Rev. D, vol. 77, p. 124015, 2008.
  * (36) S. J. Huber and T. Konstandin, “Gravitational Wave Production by Collisions: More Bubbles,” JCAP, vol. 09, p. 022, 2008.
  * (37) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Shape of the acoustic gravitational wave power spectrum from a first order phase transition,” Phys. Rev. D, vol. 96, no. 10, p. 103520, 2017.  [Erratum: Phys.Rev.D 101, 089902 (2020)].
  * (38) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Gravitational waves from the sound of a first order phase transition,” Phys. Rev. Lett., vol. 112, p. 041301, 2014.
  * (39) M. Hindmarsh, “Sound shell model for acoustic gravitational wave production at a first-order phase transition in the early Universe,” Phys. Rev. Lett., vol. 120, no. 7, p. 071301, 2018.
  * (40) D. Cutting, M. Hindmarsh, and D. J. Weir, “Vorticity, kinetic energy, and suppressed gravitational wave production in strong first order phase transitions,” Phys. Rev. Lett., vol. 125, no. 2, p. 021302, 2020.
  * (41) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Numerical simulations of acoustically generated gravitational waves at a first order phase transition,” Physical Review D, vol. 92, no. 12, p. 123009, 2015\.
  * (42) C. Caprini, R. Durrer, and G. Servant, “The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition,” JCAP, vol. 12, p. 024, 2009.
  * (43) C. Caprini et al., “Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions,” JCAP, vol. 04, p. 001, 2016.
  * (44) J. Ellis, M. Lewicki, and J. M. No, “On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal,” JCAP, vol. 04, p. 003, 2019.
  * (45) J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, “Gravitational wave energy budget in strongly supercooled phase transitions,” JCAP, vol. 06, p. 024, 2019.
  * (46) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, “Energy Budget of Cosmological First-order Phase Transitions,” JCAP, vol. 06, p. 028, 2010\.
  * (47) L. Husdal, “On effective degrees of freedom in the early universe,” Galaxies, vol. 4, no. 4, p. 78, 2016.
  * (48) B. Allen and J. D. Romano, “Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities,” Physical Review D, vol. 59, no. 10, p. 102001, 1999.
  * (49) M. Maggiore, “Gravitational wave experiments and early universe cosmology,” Physics Reports, vol. 331, no. 6, pp. 283–367, 2000.
  * (50) B. Goncharov et al., “Identifying and mitigating noise sources in precision pulsar timing data sets,” Mon. Not. Roy. Astron. Soc. , vol. 502, pp. 478–493, Mar. 2021.
  * (51) G. Hobbs, R. Edwards, and R. Manchester, “Tempo2, a new pulsar-timing package–i. an overview,” Monthly Notices of the Royal Astronomical Society, vol. 369, no. 2, pp. 655–672, 2006.
  * (52) R. T. Edwards, G. Hobbs, and R. Manchester, “Tempo2, a new pulsar timing package–ii. the timing model and precision estimates,” Monthly Notices of the Royal Astronomical Society, vol. 372, no. 4, pp. 1549–1574, 2006.
  * (53) R. van Haasteren, Y. Levin, G. Janssen, K. Lazaridis, M. Kramer, B. Stappers, G. Desvignes, M. Purver, A. Lyne, R. Ferdman, et al., “Placing limits on the stochastic gravitational-wave background using european pulsar timing array data,” Monthly Notices of the Royal Astronomical Society, vol. 414, no. 4, pp. 3117–3128, 2011.
  * (54) K. Liu, J. Verbiest, M. Kramer, B. Stappers, W. van Straten, and J. Cordes, “Prospects for high-precision pulsar timing,” Monthly Notices of the Royal Astronomical Society, vol. 417, no. 4, pp. 2916–2926, 2011.
  * (55) K. Liu, E. Keane, K. Lee, M. Kramer, J. Cordes, and M. Purver, “Profile-shape stability and phase-jitter analyses of millisecond pulsars,” Monthly Notices of the Royal Astronomical Society, vol. 420, no. 1, pp. 361–368, 2012\.
  * (56) Z. Arzoumanian, A. Brazier, S. Burke-Spolaor, S. Chamberlin, S. Chatterjee, B. Christy, J. M. Cordes, N. Cornish, K. Crowter, P. B. Demorest, et al., “The nanograv nine-year data set: observations, arrival time measurements, and analysis of 37 millisecond pulsars,” The Astrophysical Journal, vol. 813, no. 1, p. 65, 2015.
  * (57) R. van Haasteren and Y. Levin, “Understanding and analysing time-correlated stochastic signals in pulsar timing,” Monthly Notices of the Royal Astronomical Society, vol. 428, no. 2, pp. 1147–1159, 2013.
  * (58) R. van Haasteren and M. Vallisneri, “Low-rank approximations for large stationary covariance matrices, as used in the Bayesian and generalized-least-squares analysis of pulsar-timing data,” Mon. Not. Roy. Astron. Soc. , vol. 446, pp. 1170–1174, Jan. 2015.
  * (59) R. Caballero, K. Lee, L. Lentati, G. Desvignes, D. Champion, J. Verbiest, G. Janssen, B. Stappers, M. Kramer, P. Lazarus, et al., “The noise properties of 42 millisecond pulsars from the european pulsar timing array and their impact on gravitational-wave searches,” Monthly Notices of the Royal Astronomical Society, vol. 457, no. 4, pp. 4421–4440, 2016.
  * (60) L. Lentati, R. M. Shannon, W. A. Coles, J. P. Verbiest, R. van Haasteren, J. Ellis, R. Caballero, R. N. Manchester, Z. Arzoumanian, S. Babak, et al., “From spin noise to systematics: stochastic processes in the first international pulsar timing array data release,” Monthly Notices of the Royal Astronomical Society, vol. 458, no. 2, pp. 2161–2187, 2016.
  * (61) X. You, G. Hobbs, W. A. Coles, R. N. Manchester, R. Edwards, M. Bailes, J. Sarkissian, J. P. Verbiest, W. Van Straten, A. Hotan, et al., “Dispersion measure variations and their effect on precision pulsar timing,” Monthly Notices of the Royal Astronomical Society, vol. 378, no. 2, pp. 493–506, 2007.
  * (62) M. Keith, W. Coles, R. Shannon, G. Hobbs, R. Manchester, M. Bailes, N. Bhat, S. Burke-Spolaor, D. Champion, A. Chaudhary, et al., “Measurement and correction of variations in interstellar dispersion in high-precision pulsar timing,” Monthly Notices of the Royal Astronomical Society, vol. 429, no. 3, pp. 2161–2174, 2013.
  * (63) Z. Arzoumanian et al., “The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background,” Astrophys. J. Lett. , vol. 905, p. L34, Dec. 2020.
  * (64) B. Goncharov et al., “On the evidence for a common-spectrum process in the search for the nanohertz gravitational-wave background with the Parkes Pulsar Timing Array,” arXiv e-prints, p. arXiv:2107.12112, July 2021.
  * (65) J. M. Dickey, “The weighted likelihood ratio, linear hypotheses on normal location parameters,” The Annals of Mathematical Statistics, pp. 204–223, 1971.
  * (66) R. H. Cyburt, B. D. Fields, K. A. Olive, and E. Skillman, “New bbn limits on physics beyond the standard model from 4he,” Astroparticle Physics, vol. 23, no. 3, pp. 313–323, 2005.
  * (67) M. Kerr et al., “The Parkes Pulsar Timing Array project: second data release,” Publications of the Astronomical Society of Australia, vol. 37, p. e020, June 2020.
  * (68) W. Ratzinger and P. Schwaller, “Whispers from the dark side: Confronting light new physics with NANOGrav data,” SciPost Phys., vol. 10, p. 047, 2021.
  * (69) Z. Arzoumanian et al., “Searching For Gravitational Waves From Cosmological Phase Transitions With The NANOGrav 12.5-year dataset,” arXiv e-prints, p. arXiv:2104.13930, Apr. 2021.
  * (70) B. B. P. Perera et al., “The International Pulsar Timing Array: second data release,” Mon. Not. Roy. Astron. Soc. , vol. 490, pp. 4666–4687, Dec. 2019.

## Supplemental Material

###  .1 Formalism of Pulsar Intrinsic Red Noise

The spin noise (for all pulsar ToAs) is parameterized as a power-law of
frequency

$\displaystyle S_{\rm spin}(f)=\frac{A_{s}^{2}}{12\pi^{2}}\left(\frac{f}{{\rm
yr}^{-1}}\right)^{-\gamma_{s}}{\rm yr}^{3}.$ (6)

The band noise is the same as spin noise, but only for ToAs of a certain
photon frequency band,

$\displaystyle S_{\rm band}(f)=\frac{A_{b}^{2}}{12\pi^{2}}\left(\frac{f}{{\rm
yr}^{-1}}\right)^{-\gamma_{b}}{\rm yr}^{3}.$ (7)

The dispersion measure (DM) noise applies to all pulsar ToAs, but correlates
between different observing photon frequencies, as

$\displaystyle S^{ij}_{\rm DM}(f)=\frac{A_{\rm
DM}^{2}}{12\pi^{2}}\left(\frac{f}{{\rm yr}^{-1}}\right)^{-\gamma_{\rm DM}}{\rm
yr}^{3}\left(\frac{\nu_{i}}{\nu_{\rm
ref}}\right)^{-2}\left(\frac{\nu_{j}}{\nu_{\rm ref}}\right)^{-2},$ (8)

where $\nu_{i}$ and $\nu_{j}$ are photon frequencies of the $i$th and $j$th
ToA, $\nu_{\rm ref}=1400$ MHz is a reference frequency. The chromatic red
noise is

$\displaystyle S^{ij}_{\rm
chrom}(f)=\frac{A_{c}^{2}}{12\pi^{2}}\left(\frac{f}{{\rm
yr}^{-1}}\right)^{-\gamma_{c}}{\rm yr}^{3}\left(\frac{\nu_{i}}{\nu_{\rm
ref}}\right)^{-n}\left(\frac{\nu_{j}}{\nu_{\rm ref}}\right)^{-n},$ (9)

where $n$ variates from pulsar to pulsar.

The CPL process for all pulsars is

$\displaystyle S_{\rm CPL}(f)=\frac{A_{\rm
CPL}^{2}}{12\pi^{2}}\left(\frac{f}{{\rm yr}^{-1}}\right)^{-\gamma_{\rm
CPL}}{\rm yr}^{3}.$ (10)

The correlation $C_{ij}$ between ToA residuals is written as

$\displaystyle C^{ij}_{\rm red}=\int_{f_{L}}{\tt d}f\left(\sum
S^{ij}(f)\right)\cos\left(2\pi f(t_{i}-t_{j})\right).$ (11)

where $\sum$ represents the summation over all red contributions, including
noises, CPL process and gravitational wave background signals. For pulsar
noises Eq. (6-9), $f_{L}=1/T$, $T$ is the duration of observation for the
pulsar, The time-frequency method is adopted to approximate the integrals, and
30 linearly spaced frequency modes are chosen, i.e., $f_{i}=i/T$, $i=1...30$.
For the CPL process Eq. (10) and all gravitational wave background signals, we
use 20 logarithmically spaced frequency modes between $10^{-9}$ Hz and
$10^{-8}$ Hz.

Table 2: Description and priors of the model parameters. The model of pulsar intrinsic noise we use is the same as ref. goncharov2020identifying , therefore no redundant description will be given here. |  Parameter |  Prior |  Description
---|---|---|---
CPL Process Parameters |  $A_{\rm CPL}$ |  log-U[$10^{-18}$,$10^{-11}$] |  one per PTA
|  $\gamma_{\rm CPL}$ |  U[$0$,$7$] |  one per PTA
FOPT Parameters |  $\alpha$ |  log-U[$10^{-5}$,$1$] |  one per PTA
$(\beta/H_{*})^{-1}$ |  log-U[$10^{-6}$,$1$] |  one per PTA
|  $T_{*}$ (GeV) |  log-U[$10^{-6}$,$10^{2}$] |  one per PTA
Free Spectrum Parameters |  $\sqrt{\Omega_{gw}(f_{i})h^{2}}$ |  U[$10^{-10}$,$10^{0}$] |  20 per PTA
BayesEphem [12] Parameters |  $z_{\rm drift}$ |  U[$-10^{-9}$,$10^{-9}$] |  one per PTA
$\Delta M_{\rm Jupiter}$ |  N$(0,1.5\times 10^{-11})$ |  one per PTA
$\Delta M_{\rm Saturn}$ |  N$(0,8.2\times 10^{-12})$ |  one per PTA
$\Delta M_{\rm Uranus}$ |  N$(0,5.7\times 10^{-11})$ |  one per PTA
$\Delta M_{\rm Neptune}$ |  N$(0,7.9\times 10^{-11})$ |  one per PTA
${\rm PCA}_{i}$ |  U[$-0.05$,$0.05$] |  six per PTA
Figure 3: Posterior probability distributions of the FOPT parameters assuming
no CPL process (under the H2 hypothesis in Table I). For each panel the other
parameter is marginalized within its prior range. The excluded region from BBN
cyburt2005new  is shown as a grey shaded band.

###  .2 On the SMBBH Interpretation for the Common Power-Law Excess

The constraints on the FOPT model parameters in the SMBBH + FOPT hypothesis
are found to be similar with those of H3 and H4, as can be seen in Fig. 4. The
power spectral density of the SMBBH signal is written as

$\displaystyle S_{\rm SMBBH}(f)=\frac{A_{\rm
SMBBH}^{2}}{12\pi^{2}}\left(\frac{f}{{\rm yr}^{-1}}\right)^{-\gamma_{\rm
SMBBH}}{\rm yr}^{3},$ (12)

where the prior for $\log_{10}(A_{\rm SMBBH})$ is uniformly distributed from
$-18$ and $-14$ and $\gamma_{\rm SMBBH}$ is fixed at $13/3$.

Figure 4: The PPTA exclusion contours (at 95% C.L.) on FOPT parameters
$\log_{10}(T_{*}/{\rm GeV})$ and $\log_{10}[(\beta/H_{*})^{-1}]$ in the
presence of a gravitational-wave background from supermassive binary black
holes (SMBBH; in yellow) or a CPL process (green), for $\alpha=0.2$, 0.5, and
1.0. In both cases, we make use of the full-HD correlation in the analysis.