# Limit Theorems for Patterns in Ranked Tree-Child Networks

Michael Fuchs, Hexuan Liu, Tsan-Cheng Yu
Department of Mathematical Sciences
National Chengchi University
Taipei 116
Taiwan Supported by MOST under the research grant MOST-109-2115-M-004-003-MY2.

###### Abstract

We prove limit laws for the number of occurrences of a pattern on the fringe
of a ranked tree-child network which is picked uniformly at random. Our
results extend the limit law for cherries proved by Bienvenu et al. (2022).
For patterns of height $1$ and $2$, we show that they either occur frequently
(mean is asymptotically linear and limit law is normal) or sporadically (mean
is asymptotically constant and limit law is Poisson) or not all (mean tends to
$0$ and limit law is degenerate). We expect that these are the only possible
limit laws for any fringe pattern.

Keywords: Phylogenetic network, tree-child network, ranked tree-child network,
pattern, limit law, method of moments.

Mathematics subject classification (2020): 05C20, 60C05, 60F05, 92D15.

##  1 Introduction

Studying properties of shape statistics for random models that are used to
describe the evolutionary relationship between species is an important topic
in biology. For phylogenetic trees, which are used to model non-reticulate
evolution, many such studies have been performed and the stochastic behavior
of, e.g., pattern counts are known in great detail; see [4, 6, 7, 8, 14, 15,
17, 19, 20]. On the other hand, for phylogenetic networks, which are used to
model reticulate evolution, very little is known about the number of
occurrences of patterns when the networks from a given class are randomly
sampled. This is due to the fact that even counting questions for phylogenetic
networks from a given class were still open until recently; see [3, 5, 16, 10,
11, 12, 13, 18] for progress on counting questions for some of the major
classes of phylogenetic networks.

One idea which made the above mentioned counting questions easier and also
allowed the investigation of stochastic properties of shape parameters was the
idea of ranking phylogenetic networks; see [1]. (In fact, as argued in [1],
ranked networks might be more important from a practical point of view because
these networks are obtained from an evolution process; see the definition
below.) More precisely, the author of [1] defined ranked tree-child networks
and proved several results for them, e.g., they obtained a Poisson limit law
for the number of cherries when a network is sampled uniformly at random. This
is, as far as we know, the first limit result for the number of occurrences of
a pattern in phylogenetic networks. (Note that cherries were also the first
pattern for which a limit law was proved for phylogenetic trees; see [17].)

The main goal of this work is to study limit laws for the number of
occurrences of other patterns on the fringe of ranked tree-child networks. Our
results might give an indication of what to expect when studying patterns for
other classes of phylogenetic networks.

Before explaining our results, we will give definitions and recall previous
results. We start with the definition of a (rooted, binary) phylogenetic
network which is a directed acyclic graph (DAG) without double edges such that
every node falls into one of the following four categories:

  * •

A (unique) root which has in-degree $0$ and out-degree $1$;

  * •

Leaves which have in-degree $1$ and out-degree $0$ and which are bijectively
labeled by $\\{1,\ldots,n\\}$ where $n$ is their number;

  * •

Tree nodes which are nodes of in-degree $1$ and out-degree $2$;

  * •

Reticulation nodes which are nodes of in-degree $2$ and out-degree $1$.

Figure 1: (a) A phylogenetic network which satisfies the tree-child property
and is rankable; (b) The rankable tree-child network from (a) with a ranking.

See Figure 1-(a) for an example. In fact, the network from this figure is even
a tree-child network.

######  Definition 1.

A phylogenetic network is called tree-child network if every non-leaf node has
at least one child which is not a reticulation node.

For a tree-child network, we call a tree-node a branching event and a
reticulation node with its two parents a reticulation event; see Figure 2 for
the graphical depiction of these two events that we are going to use in the
sequel. (Vertical edges in this depiction will subsequently be called
lineages.)

Figure 2: The branching and reticulation event used in the construction of
ranked tree-child networks.

The main object in this paper are ranked tree-child networks which will be
defined next.

######  Definition 2.

A tree-child network is called rankable if it has recursively evolved starting
from a branching event by attaching in each step either a branching event or a
reticulation event. A rankable tree-child network together with a ranking of
its events is called a ranked tree-child network.

See Figure 1-(b) for a ranked tree-child arising from the rankable tree-child
network from Figure 1-(a).

Ranked tree-child networks have been introduced in [1] where the authors
proved a wealth of combinatorial and stochastic properties about them. Two of
these properties concern the number of occurrences of patterns in a ranked
tree-child network which is picked uniformly at random from all ranked tree-
child networks with $n$ leaves. (We call such a network a random ranked tree-
child network in the sequel.) To recall these results, we need two
definitions. First, a cherry is a tree node with both children leaves (or
equivalently, a branching event with both outgoing lineages external); a
trident is a branching event with all three outgoing lineages external. For
instance, in the ranked tree-child network from Figure 1-(b), there are $2$
cherries and no trident.

Denote by $C_{n}$ resp. $T_{n}$ the number of cherries resp. tridents in a
random ranked tree-child network with $n$ leaves. The following limit law
result was proved in [1] for $C_{n}$: as $n\rightarrow\infty$, $C_{n}$ weakly
converges to the Poisson distribution with parameter $1/4$, i.e.,

$C_{n}\stackrel{{\scriptstyle d}}{{\longrightarrow}}{\rm
Poisson}(1/4),\qquad(n\rightarrow\infty).$ (1)

On the other hand, for the number of tridents, the authors in [1] just proved
a weak law of large numbers:

$\frac{T_{n}}{n}\stackrel{{\scriptstyle{\mathbb{P}}}}{{\longrightarrow}}\frac{1}{7},\qquad(n\rightarrow\infty).$
(2)

Our first result improves this to a central limit theorem.

######  Theorem 1.

For the number of tridents $T_{n}$ in a random ranked tree-child network with
$n$ leaves, we have

$\frac{T_{n}-n/7}{\sqrt{24n/637}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N(0,1),\qquad(n\rightarrow\infty),$ (3)

where $N(0,1)$ denotes the standard normal distribution.

Note that (1) and (3) are all the limit laws of patterns of height $1$, where
the height is defined as the number of steps in the evolution process from the
definition of ranked tree-child networks. We next turn to all patterns of
height $2$ which are listed in Figure 3. For the number of occurrences of
these patterns in random ranked tree-child networks with $n$ leaves, we have
the following result.

Figure 3: All patterns of height 2. The number of occurrences of these
patterns in a ranked tree-child networks with a large number of leaves is as
follows: (a) These two do not occur; (b) These five occur only sporadically;
(c) These two occur frequently. (For details see Theorem 2.)

######  Theorem 2.

Denote by $X_{n}$ the number of occurrences of a (fixed) pattern of height $2$
in a random ranked tree-child networks with $n$ leaves. Then, we have the
following limit law results.

(A)

For the patterns in Figure 3-(a), we have that the limit law of $X_{n}$ is
degenerate. More precisely,

$X_{n}\stackrel{{\scriptstyle
L_{1}}}{{\longrightarrow}}0,\qquad(n\rightarrow\infty).$

(B)

For the patterns in Figure 3-(b), we have

$X_{n}\stackrel{{\scriptstyle d}}{{\longrightarrow}}{\rm
Poisson}(\lambda),\qquad(n\rightarrow\infty),$

where

| (b-i) | (b-ii) | (b-iii) | (b-iv) | (b-v)
---|---|---|---|---|---
$\lambda$ | $1/8$ | $1/28$ | $1/56$ | $1/14$ | $1/28$

(C)

For the patterns in Figure 3-(c), we have

$\frac{X_{n}-\mu n}{\sigma\sqrt{n}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N(0,1),\qquad(n\rightarrow\infty),$

where $(\mu,\sigma^{2})=(4/77,4575916/137582445)$ and
$(\mu,\sigma^{2})=(2/77,2930764/137582445)$ for the patterns from Figure
3-(c-i) and Figure 3-(c-ii), respectively.

Thus, the only possible limit laws for patterns of height $1$ and height $2$
are either the normal law or the Poisson law or the degenerate law. In fact,
we believe that these are the only possible limit laws for the occurrence of
any (fixed) pattern; see the discussion in Section 4. Note that is a very
different behavior from the one observed for patterns in random phylogenetic
trees which are all known to be asymptotically normal distributed; see [4].

We conclude the introduction by (briefly) explaining the proof of the above
two theorems and at the same time presenting an outline of the paper. The
proofs use the same strategy that was used for (1) in [1], namely, coupling
with a Markov chain and applying the method of moments. In addition, an
important step in the proof of Theorem 1 (as well as in the central limit
theorems of Theorem 2) will be shifting-the-mean which will considerably
reduce the complexity of computing higher central moments. The proof of
Theorem 1 will be presented in the next section. The proof of Theorem 2, which
will be given in Section 3, will be more demanding since for all patterns from
Figure 3 except for the cases (b-i), (c-i) and (c-ii), we will need to
consider the pattern together with another pattern to be able to compute
moments; for the remaining cases, namely, Figure 3-(b-i), (c-i) and (c-ii), we
will need to consider each of them with two more patterns (one of which will
be a pattern of height $3$) and the description of the Markov chain will
require many cases; this will make the computation of the moments complicated
(and thus the computations will be done with the help of Maple). The paper
will be concluded with some remarks in Section 4.

##  2 Patterns of Height 1

In this section, we are going to prove Theorem 1. Before doing so, we will
need some preliminaries.

First, we recall a random process from [1] (which was called forward
construction in that paper). Start with a branching event and recursively
attach in each consecutive step either a branching or a reticulation event.
This is done as follows: in the $\ell-1$-st step a tuple $(\ell_{1},\ell_{2})$
of the $\ell$ external lineages is picked uniformly at random and a branching
event is attached to the lineage $\ell_{1}$ if $\ell_{1}=\ell_{2}$ or a
reticulation event with incoming lineages $\ell_{1}$ and $\ell_{2}$ is created
if $\ell_{1}\neq\ell_{2}$. The process stops if $n$ external lineages are
created. Note that the resulting network is not a ranked tree-child network
since (a) leaves are not labeled and (b) an order of the incoming lineages of
a reticulation event has been fixed. However, it was proved in [1] that the
random variable counting the number of tridents (or more generally any
pattern) in the network resulting from the the above process has the same
distribution as the corresponding random variable in random ranked tree-child
networks.

By analyzing what will happen with the number of tridents by performing one
more step in the above construction, the authors in [1] obtained the following
result.

######  Lemma 1 (Bienvenu et al.; [1]).

Define a Markov process by $T_{2}=0$ and

$(T_{n+1}|T_{n}=j)=\begin{cases}j-1,&\text{with probability}\
{\displaystyle\frac{3j(3j-2)}{n^{2}}};\vspace*{0.15cm}\\\ j,&\text{with
probability}\
{1-\displaystyle\frac{3j(3j-2)+(n-3j)(n-3j-1)}{n^{2}}};\vspace*{0.15cm}\\\
j+1,&\text{with probability}\
{\displaystyle\frac{(n-3j)(n-3j-1)}{n^{2}}}.\end{cases}$ (4)

Then, $T_{n}$ has the same distribution as the number of tridents in a random
ranked tree-child network with $n$ leaves.

This was then used to compute the mean and variance of $T_{n}$ (which then in
turn implied (2)). Since we will need the result for the mean, we recall it
below.

######  Lemma 2 (Bienvenu et al.; [1]).

Let $\mu_{n}:={\mathbb{E}}(T_{n})$. Then, $m_{n}$ satisfies the recurrence

$\mu_{n+1}=\left(1-\frac{3}{n}\right)^{2}\mu_{n}+1-\frac{1}{n}$

whose solution is given by

$\mu_{n}=\frac{(15n^{3}-85n^{2}+144n-71)n}{105(n-1)(n-2)(n-3)},\qquad(n\geq
4).$ (5)

Note that this result gives that $\mu_{n}\sim n/7$ as $n\rightarrow\infty$. An
asymptotic expansion for the variance was also obtained in [1] by first
computing the second moment (which also satisfies a similar first-order
recurrence as the mean) and then using the definition of the variance. It is
actually easier to work directly with the variance which also satisfies a
similar recurrence as the mean. (The reason why this is easier is that it
already incorporates the cancellations which arise from
${\mathbb{E}}(T_{n}^{2})\sim({\mathbb{E}}(T_{n}))^{2}$; this method is the
shifting-the-mean method mentioned in the last paragraph of the previous
section.) Moreover, this also extends to higher moments.

More precisely, we define

$\phi_{n,m}:={\mathbb{E}}(T_{n}-\mu_{n})^{m}.$

Then, this sequence satisfies the following recurrence.

######  Lemma 3.

For $m\geq 2$, we have

$\phi_{n+1,m}=\left(1-\frac{3m}{n}\right)^{2}\phi_{n,m}+\psi_{n,m}$

with

$\psi_{n,m}=\sum_{j=0}^{m-1}\phi_{n,j}\Lambda_{j}(n),$

where $\Lambda_{j}(n)$ admits the complete asymptotic expansion

$\Lambda_{j}(n)\sim\sum_{\ell=0}^{\infty}\frac{\lambda_{j,\ell}}{n^{\ell}},\qquad(n\rightarrow\infty)$

with $\lambda_{j,\ell}\in{\mathbb{R}}$; in particular, $\lambda_{m-1,0}=0$ and

$\lambda_{m-2,0}=\binom{m}{2}\times\frac{24}{49}.$

_Proof._ Set $\bar{T}_{n}:=T_{n}-\mu_{n}$. From (4), we have

$\displaystyle\phi_{n+1,m}=$ $\displaystyle\
{\mathbb{E}}(\bar{T}_{n}+\mu_{n}-\mu_{n+1}-1)^{m}\frac{3(\bar{T}_{n}+\mu_{n})(3(\bar{T}_{n}+\mu_{n})-2)}{n^{2}}$
$\displaystyle+{\mathbb{E}}(\bar{T}_{n}+\mu_{n}-\mu_{n+1})^{m}$
$\displaystyle\qquad\times\left(1-\frac{3(\bar{T}_{n}+\mu_{n})(3(\bar{T}_{n}+\mu_{n})-2)+(n-3(\bar{T}_{n}+\mu_{n}))(n-3(\bar{T}_{n}+\mu_{n})-1)}{n^{2}}\right)$
$\displaystyle+{\mathbb{E}}(\bar{T}_{n}+\mu_{n}-\mu_{n+1}+1)^{m}\frac{(n-3(\bar{T}_{n}+\mu_{n}))(n-3(\bar{T}_{n}+\mu_{n})-1)}{n^{2}}.$

From this, by expanding what is inside the means by the binomial theorem, we
see that

$\psi_{n,m}=\sum_{j=0}^{m+2}\phi_{n,j}\Lambda_{j}(n).$

However, $\Lambda_{m+2}(n)=0$ (because the probabilities in (4) sum up to $1$)
and straightforward computation (best done with a computer algebra system such
as Maple) shows that $\Lambda_{m+1}(n)=\Lambda_{m}(n)=0$. Next, since Lemma 2
implies that $\mu_{n}-\mu_{n+1}={\mathcal{O}}(1)$, we have

$\Lambda_{j}(n)\sim\sum_{\ell=0}^{\infty}\frac{\lambda_{j,\ell}}{n^{\ell}},\qquad(n\rightarrow\infty)$

and by some more computations (again best done with Maple), we obtain that
$\lambda_{m-1,0}=0$ and that $\lambda_{m-2,0}$ is as claimed.

Note that the above recurrence has the (general) form

$\phi_{n+1}=\left(1-\frac{\kappa}{n}\right)^{2}\phi_{n}+\psi_{n},\qquad(n\geq\kappa+1)$
(6)

with a suitable initial value $\phi_{\kappa+1}$, where $\kappa\in{\mathbb{N}}$
and $\\{\psi_{n}\\}_{n\geq\kappa+1}$ is a given sequence. We need a general
result for such a sequence. (This result was also implicitly contained in
[1].)

######  Lemma 4.

Assume that $\phi_{n}$ satisfies (6). If $\psi_{n}\sim cn^{\alpha}$ with
$\alpha>-2\kappa-1$ a real number, then

$\phi_{n}\sim\frac{c}{2\kappa+\alpha+1}n^{\alpha+1},\qquad(n\rightarrow\infty).$

######  Remark 1.

Throughout the paper, the notation $f(n)\sim cg(n)$ (with $g(n)>0$) means that
$f(n)=cg(n)+o(g(n))$. (Note that this is the usual meaning if $c\neq 0$;
however, if $c=0$, then $f(n)=o(g(n)$.)

_Proof._ Iterating (6) gives the solution

$\phi_{n}=\binom{n-1}{\kappa}^{-2}\left(\psi_{\kappa+1}+\sum_{\ell=\kappa+1}^{n-1}\binom{\ell}{\kappa}^{2}\psi_{\ell}\right)$
(7)

Note that

$\binom{\ell}{\kappa}^{2}\sim\frac{\ell^{2\kappa}}{\kappa!^{2}},\qquad(\ell\rightarrow\infty).$
(8)

Thus, if $\psi_{n}\sim cn^{\alpha}$, then

$\psi_{\kappa+1}+\sum_{\ell=\kappa+1}^{n-1}\binom{\ell}{\kappa}^{2}\psi_{\ell}\sim\frac{c}{\kappa!^{2}}\sum_{\ell=\kappa+1}^{n-1}\ell^{2\kappa+\alpha}\sim\frac{c}{\kappa!^{2}(2\kappa+\alpha+1)}n^{2\kappa+\alpha+1}.$

Plugging this into (7) and using once more (8) gives the claimed result.

Applying the last result to the recurrence for the central moments (Lemma 3)
and using induction gives the following asymptotic result for all central
moments of $T_{n}$. (This method is sometimes refered to as moment-pumping in
the literature; see Section VII.10.1 in [9].)

######  Proposition 1.

As $n\rightarrow\infty$, the central moments of $T_{n}$ satisfy

${\mathbb{E}}\left((T_{n}-\mu_{n})^{m}\right)\sim
g_{m}\left(\frac{24}{637}\right)^{m/2}n^{m/2}.$

Here, $g_{m}$ denotes the $m$-th moment of the standard normal distribution,
i.e.,

$g_{m}={\mathbb{E}}\left(N(0,1)^{m}\right)=\begin{cases}{\displaystyle\frac{m!}{2^{m/2}(m/2)!}},&\text{if}\
m\ \text{is even};\\\ 0,&\text{if}\ m\ \text{is odd}.\end{cases}$ (9)

_Proof._ The claim (trivially) holds for $m=0$ and $m=1$. Next, we assume that
the claim holds for all $m^{\prime}<m$. In order to show it for $m$, note that
from the induction hypothesis, we have for $\psi_{n,m}$ from Lemma 3:

$\psi_{n,m}\sim\binom{m}{2}\times\frac{24}{49}\times
g_{m-2}\left(\frac{24}{637}\right)^{m/2-1}n^{m/2-1}.$

Applying Lemma 4, we obtain that

$\phi_{n,m}\sim\binom{m}{2}\times\frac{24}{49}\times\frac{g_{m-2}}{6m+m/2}\left(\frac{24}{637}\right)^{m/2-1}n^{m/2}\sim(m-1)g_{m-2}\left(\frac{24}{637}\right)^{m/2}n^{m/2}$

from which the claimed result follows since $(m-1)g_{m-2}=g_{m}$.

Theorem 1 follows now from the last proposition by the Fréchet-Shohat Theorem;
see, e.g., Section 5.30 in [2].

##  3 Patterns of Height 2

In this section, we will prove Theorem 2. This will be done in the three
paragraphs below, one for each of the three cases in Theorem 2. The proof will
proceed along similar lines as the proof from the previous section. In
addition, it will make use of the expansions from Proposition 1 which imply
that

${\mathbb{E}}(T_{n}^{m})\sim\mu_{n}^{m}\sim\frac{n^{m}}{7^{m}},\qquad(n\rightarrow\infty)$
(10)

and

${\mathbb{E}}(C_{n}^{\underline{m}}):={\mathbb{E}}(C_{n}(C_{n}-1)\cdots(C_{n}-m+1))\sim\frac{1}{4^{m}},\qquad(n\rightarrow\infty)$
(11)

which was used in [1] to prove (1).

### Degenerate Limit Laws.

We will only consider the pattern in Figure 3-(a-i); the other pattern whose
limit law is degenerate, namely the one in Figure 3-(a-ii), is treated
similarly. Our method below will be slightly more general than needed; the
reason for this is that the other cases from Theorem 2 will be proved with
similar arguments (and in these cases, the generality below is needed).

We will consider the pattern in Figure 3-(a-i) (called a pattern of type $A$
in the sequel) together with a cherry which is not contained in a pattern of
type $A$ (called a pattern of type $B$ in the sequel). Assume that a random
ranked tree-child network with $n$ leaves contains $a$ pattern of type $A$ and
$b$ pattern of type $B$, respectively. Note that each external lineage belongs
either exactly to a pattern of type $A$ or a pattern of type $B$ or to neither
of these patterns (such a lineage will be called a pattern of type $C$ in the
sequel). We now carefully list what will happen with the number of patterns of
type $A$ and type $B$ if either a branching or a reticulation event is added
in the next step of the forward construction that was described at the
beginning of the last section.

First, if a branching event is added, then we have the following cases.

| type $A$ | type $B$ | probability
---|---|---|---
type $A$ | $-1$ | $+2$ | $a/n^{2}$
$0$ | $0$ | $2a/n^{2}$
type $B$ | $+1$ | $-1$ | $2b/n^{2}$
type $C$ | $0$ | $+1$ | $(n-3a-2b)/n^{2}$
Table 1: The change of the number of patterns of type $A$ and type $B$ (first
and second column) if a branching event is attached to an external lineage
belonging to a pattern of type $A$, type $B$ or type $C$ (rows). Here, $a$ and
$b$ denote the number of patterns of type $A$ and $B$, respectively. The
probability for each cases is given in the third column if one starts with a
ranked tree-child network with $n$ leaves (and thus the probability of picking
a fixed external lineage is $1/n^{2}$).

For instance, if the branching event is attached to an external lineage from a
pattern of type $A$ and that lineage does not belong to the cherry in that
pattern, then two patterns of type $B$ are created whereas one pattern of type
$A$ was destroyed; this is the first sub-row of the first row in the above
table and the probability that this happens is given by the number of possible
choices of such a lineage ($a$) divided by the number of choices of a pair of
two external lineages ($n^{2}$). Likewise, the second sub-row of the first row
is the case where the branching event is attached to an external lineage from
a pattern of type $A$ where this lineage now belongs to the cherry of that
pattern. Similarly, the remaining rows are explained.

Next, if a reticulation event is added, then we have the cases listed in Table
2.

| type $A$ | type $B$ | probability
---|---|---|---
$A$ | $-1$ | $0$ | $6a/n^{2}$
$A\ \&\ A$ | $-2$ | $0$ | $4a(a-1)/n^{2}$
$-2$ | $+1$ | $4a(a-1)/n^{2}$
$-2$ | $+2$ | $a(a-1)/n^{2}$
$B$ | $0$ | $-1$ | $2b/n^{2}$
$B\ \&\ B$ | $0$ | $-2$ | $4b(b-1)/n^{2}$
$C\ \&\ C$ | $0$ | $0$ | $(n-3a-2b)(n-3a-2b-1)/n^{2}$
$A\ \&\ B$ | $-1$ | $-1$ | $8ab/n^{2}$
$-1$ | $0$ | $4ab/n^{2}$
$A\ \&\ C$ | $-1$ | $0$ | $4a(n-3a-2b)/n^{2}$
$-1$ | $+1$ | $2a(n-3a-2b)/n^{2}$
$B\ \&\ C$ | $0$ | $-1$ | $4b(n-3a-2b)/n^{2}$
Table 2: Columns: The change of the number of patterns of type $A$ (first
column) and type $B$ (second column); the probabilities of these changes are
contained in the last column. Rows: a single letter means that both external
lineages are picked from a pattern of that type; two letters mean that the
external lineages are chosen from (different) patterns of these two types.

For instance, the first sub-row of the sixth row of this column is explained
as follows: if the reticulation event is attached to an external lineage
belonging to the cherry of a pattern of type $A$ and an external lineage from
a pattern of type $B$, then a pattern of type $A$ and type $B$ is destroyed
(and no new pattern of type $B$ is created as will happen if the third
external lineage from the pattern of type $A$ is chosen). The probability that
this will happen is given by the number of choices of the external lineages
divided by $n^{2}$; the number of choices equals $2\cdot 2\cdot 2\cdot ab$
since once the patterns of type $A$ and $B$ are chosen ($ab$ choices) there
are $2$ choices for the external lineages in the pattern of type $A$ and $2$
in the pattern of type $B$; moreover, the factor $2$ comes from symmetry.
Similarly, the other rows are explained.

Note that Table 1 and Table 2 give the transition probabilities of the Markov
chain $(X_{n},\tilde{C}_{n})$ where $X_{n}$ and $\tilde{C}_{n}$ are the number
of patterns of type $A$ and type $B$, respectively, in a random ranked tree-
child networks of $n$ leaves. With these probabilities, we obtain the
following result.

######  Lemma 5.

Let $X_{n}$ be the number of occurrences of the pattern from Figure 3-(a-i) in
a random ranked tree-child network with $n$ leaves. Then,

${\mathbb{E}}(X_{n+1})=\left(1-\frac{3}{n}\right)^{2}{\mathbb{E}}(X_{n})+\frac{2{\mathbb{E}}(C_{n})}{n^{2}}.$
(12)

_Proof._ For $X_{n}$ and $\tilde{C}_{n}$ (see the paragraph proceeding the
lemma), we have

${\mathbb{E}}(X_{n+1})=\left(1-\frac{6}{n}+\frac{11}{n^{2}}\right){\mathbb{E}}(X_{n})+\frac{2{\mathbb{E}}(\tilde{C}_{n})}{n^{2}}.$
(13)

Note that $C_{n}=X_{n}+\tilde{C}_{n}$. Consequently, we obtain the claimed
recurrence from the above one by replacing $\tilde{C}_{n}=C_{n}-X_{n}$.

######  Remark 2.

Note that for (13), we in fact only need the first columns of Table 1 and
Table 2. However, all columns will be needed for the patterns below.

The last lemma now implies the limit law result from Theorem 2-(a).

Proof of Theorem 2-(a) for the pattern in Figure 3-(a-i). Recurrence (12) has
the form (6) with

$\psi_{n}=\frac{2{\mathbb{E}}(C_{n})}{n^{2}}\sim\frac{1}{2n^{2}},$

where the last step follows from (11). Applying Lemma 4 gives

${\mathbb{E}}(X_{n})\sim\frac{1}{10n}$

which implies the claimed result.

### Poisson Limit Laws.

Here, we will consider the patterns from Figure 3-(b). Since the proof for all
patterns except for the pattern from Figure 3-(b-i) is the same, we will only
give details for one of these four patterns, namely, the pattern in Figure
3-(b-iv).

Similarly to the treatment of the pattern from the previous paragraph, we will
use the types of patterns depicted in Figure 4. Note that each external
lineage from a ranked tree-child network belongs to exactly one of these
types.

Assume now that a ranked tree-child network contains $a$ patterns of type $A$
and $b$ patterns of type $B$. Then, we again consider the changes and
corresponding probabilities of these numbers when a branching event or
reticulation event is added in the forward construction; see Table 3 and Table
4. This again gives the transition probabilities of a Markov chain which can
be used to compute mixed moments.

Figure 4: (a) The pattern from Figure 3-(b-iv); (b) A trident which is not contained in a pattern from (a); (c) The remaining external lineages. | type $A$ | type $B$ | probability
---|---|---|---
$A$ | $-1$ | $0$ | $3a/n^{2}$
$-1$ | $+1$ | $a/n^{2}$
$B$ | $0$ | $-1$ | $3b/n^{2}$
$C$ | $0$ | $0$ | $(n-4a-3b)/n^{2}$
Table 3: The change of the number of patterns of type $A$ and type $B$ (see Figure 4) when the next event is a branching event. | type $A$ | type $B$ | probability
---|---|---|---
$A$ | $-1$ | $+1$ | $8a/n^{2}$
$0$ | $0$ | $4a/n^{2}$
$A\ \&\ A$ | $-2$ | $+1$ | $9a(a-1)/n^{2}$
$-2$ | $+2$ | $6a(a-1)/n^{2}$
$-2$ | $+3$ | $a(a-1)/n^{2}$
$B$ | $0$ | $0$ | $2b/n^{2}$
$+1$ | $-1$ | $4b/n^{2}$
$B\ \&\ B$ | $0$ | $-1$ | $9b(b-1)/n^{2}$
$C\ \&\ C$ | $0$ | $+1$ | $(n-4a-3b)(n-4a-3b-1)/n^{2}$
$A\ \&\ B$ | $-1$ | $0$ | $18ab/n^{2}$
$-1$ | $+1$ | $6ab/n^{2}$
$A\ \&\ C$ | $-1$ | $+1$ | $6a(n-4a-3b)/n^{2}$
$-1$ | $+2$ | $2a(n-4a-3b)/n^{2}$
$B\ \&\ C$ | $0$ | $0$ | $6b(n-4a-3b)/n^{2}$
Table 4: The change of the number of patterns of type $A$ and type $B$ (see
Figure 4) when the next event is a reticulation event.

######  Lemma 6.

Let $X_{n}$ be the number of occurrences of the pattern from Figure 3-(b-iv)
in a random ranked tree-child network with $n$ leaves. Then, for all $r,s\geq
0$, we have

${\mathbb{E}}(X_{n+1}^{\underline{r}}T_{n+1}^{s})=\left(1-\frac{4r+3s}{n}\right)^{2}{\mathbb{E}}(X_{n}^{\underline{r}}T_{n}^{s})+\frac{4r}{n^{2}}{\mathbb{E}}(X_{n}^{\underline{r-1}}T_{n}^{s+1})-\frac{4r(r-1)}{n^{2}}{\mathbb{E}}(X_{n}^{\underline{r-1}}T_{n}^{s})+R_{n},$

where

$R_{n}=\sum_{j=0}^{s-1}\left(\lambda_{j,0}+\frac{\lambda_{j,1}}{n}+\frac{\lambda_{j,2}}{n^{2}}\right){\mathbb{E}}(X_{n}^{\underline{r}}T_{n}^{j}).$

Here, $\lambda_{j,i}$ are constants which only depend on $r$ and $s$ with
$\lambda_{s-1,0}=s$.

_Proof._ Assume that $X_{n}=c$ and $T_{n}=d$. Then, we have $c$ patterns of
type $A$ and $d-c$ patterns of type $B$. Using the cases from Table 3 and
Table 4, we can write $X_{n+1}^{\underline{r}}T_{n+1}^{s}$ given $X_{n}=c$ and
$T_{n}=d$ as a sum of terms $(X_{n}+k)^{\underline{r}}(T_{n}+\ell)^{s}$ which
are multiplied with the probabilities of every case and where
$k\in\\{-2,-1,0,1\\}$ is the value from the first column and
$\ell\in\\{-1,0,1\\}$ is the sum of the values of the first and second column
for every case. Next, we replace $(X_{n}+k)^{\underline{r}}$ for $k=-2$ and
$k=-1$ by

$(X_{n}-2)^{\underline{r}}=\frac{(X_{n}-r)(X_{n}-r-1)}{X_{n}(X_{n}-1)}X_{n}^{\underline{r}},\qquad(X_{n}-1)^{\underline{r}}=\frac{X_{n}-r}{X_{n}}X_{n}^{\underline{r}}$

and for $k=1$ by

$(X_{n}+1)^{\underline{r}}=\frac{X_{n}+1}{X_{n}-r+1}X_{n}^{\underline{r}}.$

Finally, we expand $(T_{n}+\ell)^{s}$ by the binomial theorem and simplify the
sum of all terms of the same order in this expansion with Maple. (Note that
since $\ell\in\\{-1,0,1\\}$ this only has to be done for the first three terms
in the expansion; the other terms are the same only multiplied with different
weights). This proves the claimed result.

We are now ready to prove the Poisson limit result for the pattern from Figure
3-(b-iv).

Proof of Theorem 2-(b) for the pattern in Figure 3-(b-iv). As above, let
$X_{n}$ denote the number of occurrences of the pattern from Figure 3-(b-iv)
in a random ranked tree-child network with $n$ leaves. We will use induction
to show that for all $r,s\geq 0$:

${\mathbb{E}}(X_{n}^{\underline{r}}T_{n}^{s})\sim\frac{n^{s}}{14^{r}7^{s}},\qquad(n\rightarrow\infty),$
(14)

where the induction is with respect to the lexicographic order of $(r,s)$.

Note that the base case, namely $(0,s)$ with $s\geq 0$, is implied by (10).

Next, by the above lemma ${\mathbb{E}}(X_{n}^{\underline{r}}T_{n}^{s})$
satisfies a recurrence of the form (6) with all terms in $\psi_{n}$ being of a
smaller lexicographic order. Thus, by using the induction hypothesis,

$\psi_{n}\sim\frac{4r}{n^{2}}{\mathbb{E}}(X_{n}^{\underline{r-1}}T_{n}^{s+1})+\lambda_{s-1,0}{\mathbb{E}}(X_{n}^{\underline{r}}T_{n}^{s-1})\sim\frac{8r+7s}{14^{r}7^{s}}n^{s-1}.$

The induction claim follows now from this by applying Lemma 4.

Finally, the claimed Poisson limit law for $X_{n}$ follows from (14).

Next, we turn to the remaining pattern from Figure 3-(b), namely the pattern
in Figure 3-(b-i). This pattern is different from the other four that satisfy
a Poisson limit law because two occurrences of this pattern in a ranked tree-
child network might overlap; see Figure 5-(a). (For the pattern in this
figure, we do not distinguish the ranks of the first two events.) Thus, we now
need to consider patterns of types $A$-$D$ in order to set up the Markov chain
for proving the Poisson limit law for this pattern; see Figure 5-(b).

Assume that the number of patterns of type $A$, type $B$ and type $C$ in a
ranked tree-child network with $n$ leaves is given by $a,b$ and $c$,
respectively. Then, the changes and probabilities when adding one more event
in the forward construction are listed in Table 5 and Table 6.

Figure 5: (a) The pattern of height $3$ which contains two overlapping patterns from Figure 3-(b-i); (b) The types of patterns considered in the proof of the Poisson limit law for the pattern in Figure 3-(b-i); in order that every external lineage belongs exactly to one type, a pattern of type $B$ is not allowed to be contained in a pattern of type $A$; also, the lineage in type $D$ is not an external lineage in $A$, $B$ or $C$. | type $A$ | type $B$ | type $C$ | probability
---|---|---|---|---
$A$ | $-1$ | $0$ | $+1$ | $3a/n^{2}$
$-1$ | $+1$ | $+1$ | $2a/n^{2}$
$B$ | $0$ | $-1$ | $+1$ | $4b/n^{2}$
$C$ | $0$ | $0$ | $0$ | $2c/n^{2}$
$D$ | $0$ | $0$ | $+1$ | $(n-5a-4b-2c)/n^{2}$
Table 5: The change of the number of patterns of type $A$, type $B$ and type
$C$ (see Figure 5-(b)) when the next event is a branching event.

Using the transition probabilities from Table 5 and Table 6, we obtain now the
following lemma.

######  Lemma 7.

Denote by $Y_{n}$ and $\tilde{X}_{n}$ the number of occurrences of patterns of
type $A$ and type $B$, respectively, in a random ranked tree-child network
with $n$ leaves. Then, for $r,s,t\geq 0$, we have

$\displaystyle{\mathbb{E}}(Y_{n+1}^{\underline{r}}\tilde{X}_{n+1}^{\underline{s}}C_{n+1}^{\underline{t}})=\left(1-\frac{5r+4s+2t}{n}\right)^{2}$
$\displaystyle{\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t}})+\frac{t}{n}{\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t-1}})+\frac{4s}{n}{\mathbb{E}}(Y_{n}^{\underline{r+1}}\tilde{X}_{n}^{\underline{s-1}}C_{n}^{\underline{t}})$
$\displaystyle+\frac{4s}{n}{\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s-1}}C_{n}^{\underline{t+1}})+\frac{R_{n}}{n^{2}},$
(15)

where $R_{n}$ is given by

$\displaystyle t$
$\displaystyle(2-5r-4s-2t){\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t-1}})+4r{\mathbb{E}}(Y_{n}^{\underline{r-1}}\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t+2}})-2s(1+10r+8s+4t){\mathbb{E}}(Y_{n}^{\underline{r+1}}\tilde{X}_{n}^{\underline{s-1}}C_{n}^{\underline{t}})$
$\displaystyle+2st{\mathbb{E}}(Y_{n}^{\underline{r+1}}\tilde{X}_{n}^{\underline{s-1}}C_{n}^{\underline{t-1}})-8s{\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s-1}}C_{n}^{\underline{t+2}})+4s(2-5r-4s-2t){\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s-1}}C_{n}^{\underline{t+1}})$
$\displaystyle+4s(s-1){\mathbb{E}}(Y_{n}^{\underline{r+2}}\tilde{X}_{n}^{\underline{s-2}}C_{n}^{\underline{t}})+8s(s-1){\mathbb{E}}(Y_{n}^{\underline{r+1}}\tilde{X}_{n}^{\underline{s-2}}C_{n}^{\underline{t+1}}).$

| type $A$ | type $B$ | type $C$ | probability
---|---|---|---|---
$A$ | $-1$ | $0$ | $0$ | $20an^{2}$
$A\ \&\ A$ | $-2$ | $0$ | $0$ | $9a(a-1)/n^{2}$
$-2$ | $+1$ | $0$ | $12a(a-1)/n^{2}$
$-2$ | $+2$ | $0$ | $4a(a-1)/n^{2}$
$B$ | $0$ | $-1$ | $0$ | $12b/n^{2}$
$B\&\ B$ | $0$ | $-2$ | $0$ | $16b(b-1)/n^{2}$
$C$ | $0$ | $0$ | $-1$ | $2c/n^{2}$
$C\ \&\ C$ | $+1$ | $0$ | $-2$ | $4c(c-1)/n^{2}$
$D\ \&\ D$ | $0$ | $0$ | $0$ | $(n-5a-4b-2c)(n-5a-4b-2c-1)/n^{2}$
$A\ \&\ B$ | $-1$ | $-1$ | $0$ | $24ab/n^{2}$
$-1$ | $0$ | $0$ | $16ab/n^{2}$
$A\ \&\ C$ | $-1$ | $+1$ | $-1$ | $12ac/n^{2}$
$-1$ | $+2$ | $-1$ | $8ac/n^{2}$
$A\ \&\ D$ | $-1$ | $0$ | $0$ | $6a(n-5a-4b-2c)/n^{2}$
$-1$ | $+1$ | $0$ | $4a(n-5a-4b-2c)/n^{2}$
$B\ \&\ C$ | $0$ | $0$ | $-1$ | $16bc/n^{2}$
$B\ \&\ D$ | $0$ | $-1$ | $0$ | $8b(n-5a-4b-2c)/n^{2}$
$C\ \&\ D$ | $0$ | $+1$ | $-1$ | $4c(n-5a-4b-2c)/n^{2}$
Table 6: The change of the number of patterns of type $A$ and type $B$ (see
Figure 5-(b)) when the next event is a reticulation event.

_Proof._ Using the results from Table 5 and Table 6,
${\mathbb{E}}(Y_{n+1}^{\underline{r}}\tilde{X}_{n+1}^{\underline{s}}C_{n+1}^{\underline{t}})$
given $Y_{n},\tilde{X}_{n}$ and $C_{n}$ can be written as a sum of terms of
the form

$(Y_{n}+k)^{\underline{r}}(\tilde{X}_{n}+\ell)^{\underline{s}}(C_{n}+m)^{\underline{t}}$

which are multiplied with the probabilities. Here, $k,\ell$ and $m$ are
suitable integers, e.g., for the contribution from the second sub-row of
column $A\ \&\ C$ in Table 6, we have $k=-1,\ell=2,m=-1$. Then, we rewrite,
e.g., this term as

$(Y_{n}-1)^{\underline{r}}(\tilde{X}_{n}+2)^{\underline{s}}(C_{n}-1)^{\underline{t}}=\frac{(Y_{n}-r)(\tilde{X}_{n}+2)(\tilde{X}_{n}+1)(C_{n}-t)}{Y_{n}(\tilde{X}_{n}-s)(\tilde{X}_{n}-s-1)C_{n}}Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t}}$

and similar for the other terms. The rest of the proof is just a long
computation (which is best done with the help of Maple).

The last lemma implies the following result which contains the Poisson limit
law for the pattern in Figure 3-(b-i).

######  Proposition 2.

  * (a)

Let $Y_{n}$ be the number of occurrences of the pattern from Figure 5-(a) in a
random ranked tree-child network with $n$ leaves. Then,

$Y_{n}\stackrel{{\scriptstyle
L_{1}}}{{\longrightarrow}}0,\qquad(n\rightarrow\infty).$

  * (b)

Let $X_{n}$ be the number of occurrences of the pattern from Figure 3-(b-i) in
a random ranked tree-child network with $n$ leaves. Then,

$(X_{n},C_{n})\stackrel{{\scriptstyle
d}}{{\longrightarrow}}(X,C),\qquad(n\rightarrow\infty),$

where $X$ and $C$ are independent Poisson random variables with parameters
$1/8$ and $1/4$, respectively.

_Proof._ We use induction with respect to the lexicographic order on $(s,r,t)$
to show that for all $r,s,t\geq 0$:

${\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t}})\sim
0^{r}\cdot\left(\frac{1}{8}\right)^{s}\cdot\left(\frac{1}{4}\right)^{t},$ (16)

where $\tilde{X}_{n}$ is the random variable from Lemma 7 and we use the
convention that $0^{0}:=1$. Note that the second, third and fourth term on the
right hand side of (15) and all terms in $R_{n}$ are of a smaller
lexicographic order than $(s,r,t)$. Also note that (15) has the form (6).

Now, first the induction base holds because of (11). Next, assume that the
claim holds for all sequences $(s^{\prime},r^{\prime},t^{\prime})$ which are
lexicographic smaller than $(s,r,t)$. In order to prove the claim for
$(s,r,t)$, we make a case distinction.

First, if $r>0$, then the induction hypothesis implies that (15) satisfies (6)
with $\psi_{n}=o(1/n)$. Thus, from Lemma 4, we obtain that
${\mathbb{E}}(Y_{n}^{\underline{r}}\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t}})=o(1)$
which proves the claim in this case.

Secondly, if $r=0$, then again by the induction claim, (15) is of the form (6)
with

$\psi_{n}\sim\frac{t}{n}{\mathbb{E}}(\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t-1}})+\frac{4s}{n}{\mathbb{E}}(\tilde{X}_{n}^{\underline{s-1}}C_{n}^{\underline{t+1}})\sim\frac{8s+4t}{8^{s}4^{t}n}.$

Thus, by Lemma 4

${\mathbb{E}}(\tilde{X}_{n}^{\underline{s}}C_{n}^{\underline{t}})\sim\frac{1}{8^{s}4^{t}}.$

which also proves the claim in this case.

Finally, observe that (16) implies part (a). Moreover, it also implies that

$(\tilde{X}_{n},C_{n})\stackrel{{\scriptstyle
d}}{{\longrightarrow}}(\tilde{X},C),$

where $\tilde{X}$ and $C$ are independent Poisson random variables with
parameters $1/8$ and $1/4$, respectively. Part (b) follows from this by part
(a), the relation $X_{n}=\tilde{X}_{n}+2Y_{n}$ and Slutsky’s theorem.

### Normal Limit Laws.

In this paragraph, we consider the patterns in Figure 3-(c) which according to
Theorem 2 are both normal distributed. We will only give details for the
pattern in Figure 3-(c-i), the other pattern in Figure 3-(c) is treated
similarly.

First note that two of the patterns in Figure 3-(c-i) again can overlap in a
ranked tree-child network; see Figure 6-(a). (We do not distinguish the ranks
of the first two events in this figure.) Thus, similar to the pattern in
Figure 3-(b-i), we will consider four types of patterns in the sequel; see
Figure 6-(b). Assume that a ranked tree-child network with $n$ leaves contains
$a,b$ and $c$ patterns of type $A,B$ and $C$, respectively. We again have to
list the changes and probabilities if one more event is attached in the
forward construction, however, this time there are many cases. We list all of
them in Tables 7-9.

Using these tables, we can first derive the means of the number of occurrences
of the patterns from Figure 3-(c-i) and Figure 6-(a).

######  Lemma 8.

Denote by $X_{n}$ and $Y_{n}$ the number of occurrences of the patterns from
Figure 3-(c-i) and Figure 6-(a), respectively. Let
$\rho_{n}:={\mathbb{E}}(X_{n})$ and $\tau_{n}:={\mathbb{E}}(Y_{n})$. Then,

$\rho_{n+1}=\left(1-\frac{5}{n}\right)^{2}\rho_{n}+\frac{4}{n}{\mathbb{E}}(T_{n})-\frac{12}{n^{2}}{\mathbb{E}}(T_{n})$

and

$\tau_{n+1}=\left(1-\frac{7}{n}\right)^{2}\tau_{n}-\frac{4}{n^{2}}{\mathbb{E}}(T_{n})+\frac{4}{n^{2}}{\mathbb{E}}(T_{n}^{2}).$

The solutions of the above two recurrences are given by

$\rho_{n}=\frac{(1080n^{5}-16668n^{4}+96992n^{3}-261735n^{2}+319471n-135654)n}{20790(n-1)(n-2)(n-3)(n-4)(n-5)},\qquad(n\geq
6)$ (17)

and, for $n\geq 7$,

$\tau_{n}=\frac{2(4290n^{7}-125730n^{6}+1509970n^{5}-9550275n^{4}+33968326n^{3}-66128140n-24510098)n}{1576575(n-1)(n-2)(n-3)(n-4)(n-5)(n-6)(n-7)}$.

(18)

Figure 6: (a) The pattern of height $3$ which contains two overlapping
patterns from Figure 3-(c-i); (b) The types of pattern considered in the proof
of the normal limit law for the pattern in Figure 3-(c-i); in order that every
external lineage belongs exactly to one type, a pattern of type $B$ resp. $C$
is not allowed to be contained in a pattern of type $A$ resp. $B$ and $A$;
also, the lineage in type $D$ is not an external lineage in $A$, $B$ or $C$.

_Proof._ Assume that $Y_{n}=d,X_{n}=e$ and $T_{n}=f$. Then, we have $d$
patterns of type $A$, $e-2d$ patterns of type $B$ and $f+e-d$ patterns of type
$C$.

The recurrence for $\mu_{n}$ and $\rho_{n}$ are derived by considering the
transition probabilities from the Tables 7-9. Note that both recurrences are
of type (6) whose exact solution is given by (7). Using the result from Lemma
2, a corresponding result for ${\mathbb{E}}(T_{n}^{2})$ (best computed with
Maple) and straightforward computations (again best done with Maple) gives the
claimed results for $\rho_{n}$ and $\tau_{n}$.

| type $A$ | type $B$ | type $C$ | probability
---|---|---|---|---
$A$ | $-1$ | $0$ | $0$ | $3a/n^{2}$
$-1$ | $+1$ | $0$ | $4a/n^{2}$
$B$ | $0$ | $-1$ | $0$ | $3b/n^{2}$
$0$ | $-1$ | $+1$ | $2b/n^{2}$
$C$ | $0$ | $0$ | $-1$ | $3c/n^{2}$
$D$ | $0$ | $0$ | $0$ | $(n-7a-5b-3c)/n^{2}$
Table 7: The change of the number of patterns of type $A,B$ and $C$ (see Figure 6-(b)) when the next event is a branching event. | type $A$ | type $B$ | type $C$ | probability
---|---|---|---|---
$A$ | $-1$ | $0$ | $+1$ | $14a/n^{2}$
$-1$ | $0$ | $+2$ | $8a/n^{2}$
$-1$ | $+1$ | $0$ | $16a/n^{2}$
$-1$ | $+1$ | $+1$ | $4a/n^{2}$
$A\ \&\ A$ | $-1$ | $0$ | $0$ | $4a(a-1)/n^{2}$
$-2$ | $0$ | $+1$ | $a(a-1)/n^{2}$
$-2$ | $+1$ | $0$ | $4a(a-1)/n^{2}$
$-2$ | $+1$ | $+1$ | $8a(a-1)/n^{2}$
$-2$ | $+2$ | $0$ | $16a(a-1)/n^{2}$
$-2$ | $+2$ | $+1$ | $16a(a-1)/n^{2}$
$B$ | $0$ | $0$ | $0$ | $8b/n^{2}$
$0$ | $-1$ | $+1$ | $10b/n^{2}$
$0$ | $-1$ | $+2$ | $2b/n^{2}$
$B\ \&\ B$ | $0$ | $-1$ | $0$ | $4b(b-1)/n^{2}$
$0$ | $-1$ | $+1$ | $8b(b-1)/n^{2}$
$0$ | $-2$ | $+1$ | $b(b-1)/n^{2}$
$0$ | $-2$ | $+2$ | $4b(b-1)/n^{2}$
$0$ | $-2$ | $+3$ | $4b(b-1)/n^{2}$
$+1$ | $-2$ | $0$ | $4b(b-1)/n^{2}$
$C$ | $0$ | $0$ | $0$ | $6c/n^{2}$
$C\ \&\ C$ | $0$ | $0$ | $-1$ | $c(c-1)/n^{2}$
$0$ | $+1$ | $-2$ | $4c(c-1)/n^{2}$
$+1$ | $0$ | $-2$ | $4c(c-1)/n^{2}$
$D\ \&\ D$ | $0$ | $0$ | $+1$ | $(n-7a-5b-3c)(n-7a-5b-3c-1)/n^{2}$
Table 8: The change of the number of patterns of type $A,B$ and $C$ (see Figure 6-(b)) when the next event is a reticulation event which is attached to one or two patterns of type $X$ with $X\in\\{A,B,C,D\\}$. | type $A$ | type $B$ | type $C$ | probability
---|---|---|---|---
$A\ \&\ B$ | $0$ | $-1$ | $0$ | $8ab/n^{2}$
$-1$ | $-1$ | $+1$ | $2ab/n^{2}$
$-1$ | $-1$ | $+2$ | $4ab/n^{2}$
$-1$ | $0$ | $0$ | $8ab/n^{2}$
$-1$ | $0$ | $+1$ | $16ab/n^{2}$
$-1$ | $0$ | $+2$ | $16ab/n^{2}$
$-1$ | $+1$ | $0$ | $16ab/n^{2}$
$A\ \&\ C$ | $0$ | $0$ | $-1$ | $8ac/n^{2}$
$-1$ | $0$ | $0$ | $2ac/n^{2}$
$-1$ | $+1$ | $-1$ | $8ac/n^{2}$
$-1$ | $+1$ | $0$ | $8ac/n^{2}$
$-1$ | $+2$ | $-1$ | $16ac/n^{2}$
$A\ \&\ D$ | $-1$ | $0$ | $+1$ | $2a(n-7a-5b-3c)/n^{2}$
$-1$ | $+1$ | $0$ | $4a(n-7a-5b-3c)/n^{2}$
$-1$ | $+1$ | $+1$ | $8a(n-7a-5b-3c)/n^{2}$
$B\ \&\ C$ | $0$ | $-1$ | $0$ | $2bc/n^{2}$
$0$ | $-1$ | $+1$ | $4bc/n^{2}$
$0$ | $0$ | $-1$ | $8bc/n^{2}$
$0$ | $0$ | $0$ | $8bc/n^{2}$
$+1$ | $-1$ | $-1$ | $8bc/n^{2}$
$B\ \&\ D$ | $0$ | $-1$ | $+1$ | $2b(n-7a-5b-3c)/n^{2}$
$0$ | $-1$ | $+2$ | $4b(n-7a-5b-3c)/n^{2}$
$0$ | $0$ | $0$ | $4b(n-7a-5b-3c)/n^{2}$
$C\ \&\ D$ | $0$ | $0$ | $0$ | $2c(n-7a-5b-3c)/n^{2}$
$0$ | $+1$ | $-1$ | $4c(n-7a-5b-3c)/n^{2}$
Table 9: The change of the number of patterns of type $A,B$ and $C$ (see
Figure 6-(b)) when the next event is a reticulation event which is attached to
a pattern of type $X$ and a pattern of type $Y$ with
$(X,Y)\in\\{(A,B),(A,C),(A,D),(B,C),(B,D),(C,D)\\}$.

We next shift the means and derive the recurrence for the mixed moments.

######  Lemma 9.

Denote by $X_{n}$ and $Y_{n}$ the number of occurrences of the patterns from
Figure 3-(c-i) and Figure 6-(a), respectively, in a random ranked tree-child
network with $n$ leaves. Moreover, set
$\mu_{n}:={\mathbb{E}}(T_{n}),\rho_{n}:={\mathbb{E}}(X_{n}),\tau_{n}:={\mathbb{E}}(Y_{n})$
and
$\bar{T}_{n}:=T_{n}-\mu_{n},\bar{X}_{n}:=X_{n}-\rho_{n},\bar{Y}_{n}:=Y_{n}-\tau_{n}$.
Then, for all $r,s,t\geq 0$, we have

${\mathbb{E}}(\bar{Y}_{n+1}^{r}\bar{X}_{n+1}^{s}\bar{T}_{n+1}^{t})=\left(1-\frac{7r+5s+3t}{n}\right)^{2}{\mathbb{E}}(\bar{Y}_{n}^{r}\bar{X}_{n}^{s}\bar{T}_{n}^{t})+R_{n}$
(19)

with

$R_{n}=\sum_{(s^{\prime},r^{\prime},t^{\prime})}{\mathbb{E}}(\bar{Y}_{n}^{r^{\prime}}\bar{X}_{n}^{s^{\prime}}\bar{T}_{n}^{t^{\prime}})\Lambda_{r^{\prime},s^{\prime},t^{\prime}}(n),$
(20)

where the sum runs over $(s^{\prime},r^{\prime},t^{\prime})$ which are of a
smaller lexicographic order than $(s,r,t)$ and
$\Lambda_{r^{\prime},s^{\prime},t^{\prime}}(n)$ admits the complete asymptotic
expansion:

$\Lambda_{r^{\prime},s^{\prime},t^{\prime}}(n)\sim\sum_{\ell=0}^{\infty}\frac{\lambda_{r^{\prime},s^{\prime},t^{\prime},\ell}}{n^{\ell}},\qquad(n\rightarrow\infty).$
(21)

Moreover, all terms in (20) with
$(r^{\prime}+s^{\prime}+t^{\prime})/2-\ell\geq(r+s+t)/2-1$ are given by

$\displaystyle\frac{4s}{n}$
$\displaystyle{\mathbb{E}}(\bar{Y}_{n}^{r}\bar{X}_{n}^{s-1}\bar{T}_{n}^{t+1})+\frac{8r}{7n}{\mathbb{E}}(\bar{Y}_{n}^{r-1}\bar{X}_{n}^{s}\bar{T}_{n}^{t+1})$
$\displaystyle+\binom{r}{2}\frac{80092}{540225}{\mathbb{E}}(\bar{Y}_{n}^{r-2}\bar{X}_{n}^{s}\bar{T}_{n}^{t})+\binom{s}{2}\frac{21916}{29645}{\mathbb{E}}(\bar{Y}_{n}^{r}\bar{X}_{n}^{s-2}\bar{T}_{n}^{t})+\binom{t}{2}\frac{24}{49}{\mathbb{E}}(\bar{Y}_{n}^{r}\bar{X}_{n}^{s}\bar{T}_{n}^{t-2})$
$\displaystyle-\frac{128}{539}st{\mathbb{E}}(\bar{Y}_{n}^{r}\bar{X}_{n}^{s-1}\bar{T}_{n}^{t-1})-\frac{32}{343}rt{\mathbb{E}}(\bar{Y}_{n}^{r-1}\bar{X}_{n}^{s}\bar{T}_{n}^{t-1})+\frac{712}{3773}rs{\mathbb{E}}(\bar{Y}_{n}^{r-1}\bar{X}_{n}^{s-1}\bar{T}_{n}^{t}).$
(22)

_Proof._ Assume that $Y_{n},X_{n}$ and $T_{n}$ are given. Then, using the
cases listed in the Tables 7-9, we can write the conditional expectation of
$\bar{Y}_{n+1}^{r}\bar{X}_{n+1}^{s}\bar{T}_{n+1}^{t}$ given $Y_{n},X_{n}$ and
$T_{n}$ as a sum of terms of the form

$(\bar{Y}_{n}+\tau_{n}-\tau_{n+1}+k)^{r}(\bar{X}_{n}+\rho_{n}-\rho_{n+1}+\ell)^{s}(\bar{T}_{n}+\mu_{n}-\mu_{n+1}+m)^{t}$
(23)

multiplied with the probabilities from the Tables 7-9 (where
$a=\bar{Y}_{n}+\tau_{n},b=\bar{X}_{n}+\rho_{n}-2(\bar{Y}_{n}+\tau_{n}),c=\bar{T}_{n}+\mu_{n}+\bar{X}_{n}+\rho_{n}-\bar{Y}_{n}-\tau_{n}$).
Here, $k,\ell,m$ are integers depending on the case considered. For
convenience, we will call this sum $S$ in the sequel.

Now, use the expansion:

$(\bar{X}_{n}+\rho_{n}-\rho_{n+1}+\ell)^{s}=\sum_{j=0}^{s}\binom{s}{j}\bar{X}_{n}^{j}(\rho_{n}-\rho_{n+1}+\ell)^{s-j}.$
(24)

First, by replacing the middle term in (23) by $\bar{X}_{n}^{s}$ (the first
term in the expansion (24)) and using Maple to sum up all these terms
multiplied with their probabilities (which contain at most $\bar{X}_{n}^{2}$),
we find that $\bar{X}_{n}^{s+2}$ and $\bar{X}_{n}^{s+1}$ do not appear in this
sum. Likewise, by replacing the middle term in (23) by
$s\bar{X}_{n}^{s-1}(\rho_{n}-\rho_{n+1}+\ell)$ (the second term in the
expansion (24)), we see that $\bar{X}_{n}^{s+1}$ does not appear. Thus, the
highest power of $\bar{X}_{n}$ in the sum $S$ is $\bar{X}_{n}^{s}$. Next, by
collecting these highest terms and repeating the above argument with the first
term in (23), we see that terms with $\bar{Y}_{n}^{r+2}\bar{X}_{n}^{s}$ and
$\bar{Y}_{n}^{r+1}\bar{X}_{n}^{s}$ do not occur in $S$. Finally, a similar
line of reasoning shows that terms with
$\bar{Y}_{n}^{r}\bar{X}_{n}^{s}\bar{T}_{n}^{t+2}$ and
$\bar{Y}_{n}^{r}\bar{X}_{n}^{s}\bar{T}_{n}^{t+1}$ do not appear in $S$ as
well. Overall, this shows that sum in (20) is over the indicated range.

Next, the claimed expansion for
$\Lambda_{r^{\prime},s^{\prime},t^{\prime}}(n)$ follows by expanding
$\tau_{n}-\tau_{n+1},\rho_{n}-\rho_{n+1}$ and $\mu_{n}-\mu_{n+1}$ (see (5),
(17) and (18), respectively; note that these expansions are all of the form
(21)) and pointing out that the probabilities might contain factors of the
form $\tau_{n}^{2},\rho_{n}^{2},\mu_{n}^{2},\tau_{n}\rho_{n},\tau_{n}\mu_{n},$
or $\rho_{n}\mu_{n}$ which are however divided by $n^{2}$ and thus expanding
them gives also terms of the form (21).

Finally, in order to find all terms with
$(r^{\prime}+s^{\prime}+t^{\prime})/2-\ell\geq(r+s+t)/2-1$, we proceed as
follows: if we expand the first factor in (23) and retain only the terms which
contain
$\bar{Y}_{n}^{r},\bar{Y}_{n}^{r-1},\bar{Y}_{n}^{r-2},\bar{Y}_{n}^{r-3}$ and
${\bar{Y}}_{n}^{r-4}$, then we see that we only loose terms
$\bar{Y}_{n}^{r^{\prime}}\bar{X}_{n}^{s^{\prime}}\bar{T}_{n}^{t^{\prime}}$ in
$S$ with

$(r^{\prime}+s^{\prime}+t^{\prime})/2-\ell\leq(r^{\prime}+s^{\prime}+t^{\prime})/2\leq(r-5+s+t+2)/2<(r+s+t)/2-1.$

Similarly, we only loose smaller order terms when we just keep the terms with
$\bar{X}_{n}^{s},\ldots,\bar{X}_{n}^{s-4}$ and the terms with
$\bar{T}_{n}^{t},\ldots,\bar{T}_{n}^{t-4}$ in the expansion of the second and
third factor in (23), respectively. Thus, we only need to retain a fixed
number (which does not depend on $r,s$ and $t$) of terms in each of the terms
of (23) from $S$. Then, the rest of the computation can then be done with
Maple.

The recurrence from the above lemma can now be used to prove the normal limit
law of the pattern in Figure 3-(c-i) (which completes the proof of Theorem 2).
In fact, we have a more general result.

######  Proposition 3.

Let the notation be as in Lemma 9. Then, as $n\rightarrow\infty$,

$\frac{1}{n}(Y_{n}-{\mathbb{E}}(Y_{n}),X_{n}-{\mathbb{E}}(X_{n}),T_{n}-{\mathbb{E}}(T_{n}))\stackrel{{\scriptstyle
d}}{{\longrightarrow}}N({\bf 0},\Sigma),$

where $N({\bf 0},\Sigma)$ denotes a trivariate normal distribution with mean
vector ${\bf 0}$ and covariance matrix

$\Sigma=\begin{pmatrix}\frac{1002796}{203664825}&\frac{433528}{62537475}&-\frac{32}{13377}\\\\[6.0pt]
\frac{433528}{62537475}&\frac{4575916}{137582445}&-\frac{608}{119119}\\\\[6.0pt]
-\frac{32}{13377}&-\frac{608}{119119}&\frac{24}{637}\end{pmatrix}.$

Proof of Proposition 3. We use induction with respect to the lexicographic
order of $(s,r,t)$ to show that for all $r,s,t\geq 0$:

${\mathbb{E}}(\bar{Y}_{n}^{r}\bar{X}_{n}^{s}\bar{T}_{n}^{t})\sim
c_{r,s,t}n^{(r+s+t)/2},$

where $c_{r,s,t}:={\mathbb{E}}(N_{1}^{r}N_{2}^{s}N_{3}^{t})$ with
$N=(N_{1},N_{2},N_{3})$ the trivariate normal distribution with mean vector
${\bf 0}$ and covariance matrix $\Sigma$.

First, the base case is implied by Proposition 1. Thus, we may assume that the
claim holds for all sequences $(s^{\prime},r^{\prime},t^{\prime})$ which are
lexicographic smaller than $(s,r,t)$.

In order to prove the claim for $(s,r,t)$, we observe that (19) has the form
(6) with the sequence $\psi_{n}$ satisfying
$\psi_{n}\sim\tilde{c}_{r,s,t}n^{(r+s+t)/2-1}$ where

$\displaystyle\tilde{c}_{r,s,t}:=$ $\displaystyle
4sc_{r,s-1,t+1}+\frac{8r}{7}c_{r-1,s,t+1}+\binom{r}{2}\frac{80092}{540225}c_{r-2,s,t}+\binom{s}{2}\frac{21916}{29645}c_{r,s-2,t}+\binom{t}{2}\frac{24}{49}c_{r,s,t-2}$
$\displaystyle\quad-\frac{128}{539}stc_{r,s-1,t-1}-\frac{32}{343}rtc_{r-1,s,t-1}+\frac{712}{3773}rsc_{r-1,s-1,t}.$
(25)

Now, recall that by Isserlis’ theorem, we have

$c_{r,s-1,t+1}=\sum_{{\mathcal{P}}}\prod_{\\{P_{i},P_{j}\\}\in{\mathcal{P}}}{\rm
Cov}(P_{i},P_{j}),$

where the sum runs over all partitions of the multiset
$\\{N_{1},\ldots,N_{1},N_{2},\ldots,N_{2},N_{3},\ldots,N_{3}\\}$ into pairs
where $N_{1}$ is repeated $r$ times, $N_{2}$ is repeated $s-1$ times and
$N_{3}$ is repeated $t+1$ times. By fixing one $N_{3}$ and considering its
pairing with all other elements, we obtain that

$c_{r,s-1,t+1}=t\Sigma_{3,3}c_{r,s-1,t-1}+(s-1)\Sigma_{3,2}c_{r,s-2,t}+r\Sigma_{3,1}c_{r-1,s-1,t}.$

Moreover, by a similar argument for $c_{r-1,s,t+1}$, we have

$c_{r-1,s,t+1}=t\Sigma_{3,3}c_{r-1,s,t-1}+(s-1)\Sigma_{3,2}c_{r-1,s-1,t}+r\Sigma_{3,2}c_{r-2,s,t}.$

Plugging this into (25), we obtain that

$\displaystyle\tilde{c}_{r,s,t}=$
$\displaystyle\binom{r}{2}\frac{1002796}{7022925}c_{r-2,s,t}+\binom{s}{2}\frac{4575916}{6551545}c_{r,s-2,t}+\binom{t}{2}\frac{24}{49}c_{r,s,t-2}$
$\displaystyle\quad-\frac{608}{7007}stc_{r,s-1,t-1}-\frac{32}{637}rtc_{r-1,s,t-1}+\frac{433528}{2501499}rsc_{r-1,s-1,t}.$
(26)

Next, again by Isserlis’ theorem, we have the recurrences:

$\displaystyle c_{r,s,t}$
$\displaystyle=(r-1)\Sigma_{1,1}c_{r-2,s,t}+s\Sigma_{1,2}c_{r-1,s-1,t}+t\Sigma_{1,3}c_{r-1,s,t-1};$
$\displaystyle c_{r,s,t}$
$\displaystyle=r\Sigma_{2,1}c_{r-1,s-1,t}+(s-1)\Sigma_{2,2}c_{r,s-2,t}+t\Sigma_{2,3}c_{r,s-1,t-1};$
$\displaystyle c_{r,s,t}$
$\displaystyle=r\Sigma_{3,1}c_{r-1,s,t-1}+s\Sigma_{3,2}c_{r,s-1,t-1}+(t-1)\Sigma_{3,3}c_{r,s,t-2}.$

Multiplying the first by $29r/2$, the second by $21s/2$ and the third by
$13t/2$ gives as right-hand side exactly (26). Thus,

$\tilde{c}_{r,s,t}=\left(\frac{29}{2}r+\frac{21}{2}s+\frac{13}{2}t\right)c_{r,s,t}.$

Overall, we have shown that (19) has the form (6) with

$\psi_{n}\sim\left(\frac{29}{2}r+\frac{21}{2}s+\frac{13}{2}t\right)c_{r,s,t}n^{(r+s+t)/2-1}.$

The induction claim follows now from this by applying Lemma 4. This completes
the proof.

##  4 Conclusion

The main purpose of this paper was to study the number of occurrences of
patterns on the fringe of ranked tree-child networks. More precisely, we
strengthened the weak law of large numbers for the number of tridents from [1]
by proving a central limit theorem (which completed the classification of the
limit laws of patterns of height $1$) and studied the limit laws of all
patterns of height $2$. There are only three possible limit laws: normal,
Poisson and degenerate; accordingly, pattern either occur frequently (normal
pattern), sporadically (Poisson pattern) or not all (degenerate pattern). We
think that this behavior will persist for patterns of any height.

In fact, we can formulate a conjecture for the limit law of any fringe pattern
which is defined as a connected substructure of a ranked tree-child network
which has entirely evolved from a fixed set of lineages by consecutively
adding branching and reticulation events. The limit law of such a pattern can
be obtained recursively as follows.

######  Conjecture 1.

Let $F$ be a fringe pattern. Denote by $P$ resp. $P_{1}$ and $P_{2}$ the
patterns which are obtained from it by removing the last event. (Here, the
second case is only possible if the last event is a reticulation event and the
pattern gets disconnected when this event is removed.) Then, we have the
following cases.

  * (a)

If $P$ is a normal pattern, then $F$ is a Poisson pattern; in all other cases
for $P$, the pattern $F$ is a degenerate pattern.

  * (b)

If $P_{1}$ and $P_{2}$ are both normal patterns, then $F$ is also a normal
pattern; if $P_{1}$ is a normal pattern and $P_{2}$ is a Poisson pattern or
vice versa, then $F$ is a Poisson pattern; in all remaining cases for $P_{1}$
and $P_{2}$, $F$ is a degenerate pattern.

This conjecture is easily seen to be consistent with the limit laws for the
patterns from Figure 3; see Theorem 2. It is also consistent with the limit
law results for the two patterns of height $3$ above: the pattern from Figure
5-(a) was shown to be degenerate (see Proposition 2, (a)) and indeed the
pattern splits into two Poisson patterns when the last event is removed;
likewise, the pattern from Figure 6-(a) is a normal pattern (see Proposition
3) and it splits into two normal patterns when the last event is removed.

For the proof of the results in this paper (Theorem 1 and Theorem 2) we
extended the approach from [1] which was based on coupling with a Markov chain
and the method of moments. The method can be applied to prove further cases of
the above conjecture; however, technical details become more and more
demanding for patterns of increasing height. Thus, in order to proof the
conjecture in its full generality, we think that a new (less computation-
intensive) approach has to be devised.

We conclude by pointing out that one of the main contributions of the current
study is that it constitutes the first such study for a class of phylogenetic
networks. Indeed, it would be interesting to undertake similar studies for
other classes of phylogenetic networks, e.g., the class of tree-child
networks. However, for this class, almost nothing is known so far about the
occurrence of patterns. In fact, the only result we are aware of was proved in
[16] where the authors showed that the number of cherries is $o(n)$ for almost
all tree-child networks with $n$ leaves. Thus, a (wild) guess would be that
the limit law of this pattern is again the Poisson law; however, tools for
proving this (as well as deriving limit laws for other patterns) are
completely lacking. We leave this as an open problem.

### Acknowledgments.

We thank Yu-Sheng Chang and Guan-Ru Yu for joining the discussion on the
research presented in this paper.

## References

  * [1] F. Bienvenu, A. Lambert, M. Steel. Combinatorial and stochastic properties of ranked tree-child networks, Random Struct. Algor., in press. https://doi.org/10.1002/rsa.21048
  * [2] P. Billingsley. Probability and Measure, third edition, Wiley Series in Probability and Mathematical Statistics, A Wiley-Interscience Publication, John Wiley & Sons, Inc., New York, 1995.
  * [3] G. Cardona and L. Zhang (2020). Counting tree-child networks and their subclasses, J. Comput. Syst. Sci., 114, 84–104.
  * [4] H. Chang and M. Fuchs (2010). Limit theorems for patterns in phylogenetic trees, J. Math. Biol., 60:4, 481–512.
  * [5] Y.-S. Chang, M. Fuchs, H. Liu, M. Wallner, G.-R. Yu. Enumeration of $d$-combining tree-child networks, submitted.
  * [6] K. P. Choi, G. Kaur, T. Wu (2021). On asymptotic joint distributions of cherries and pitchforks for random phylogenetic trees, J. Math. Biol., 83:4, Paper No. 40, 34 pp.
  * [7] K. P. Choi, A. Thompson, T. Wu (2020). On cherry and pitchfork distributions of random rooted and unrooted phylogenetic trees, Theor. Popul. Biol., 132, 92–104.
  * [8] F. Disanto and T. Wiehe (2013). Exact enumeration of cherries and pitchforks in ranked trees under the coalescent model, Math. Biosci., 242:2, 195–200.
  * [9] P. Flajolet and R. Sedgewick. Analytic Combinatorics, 1st edition, Cambridge University Press, Cambridge, 2009.
  * [10] M. Fuchs, H. Liu, G.-Ru. Yu. A short note on the exact counting of tree-child networks, arXiv:2110.03842.
  * [11] M. Fuchs, G.-R. Yu, L. Zhang (2021). On the asymptotic growth of the number of tree-child networks, European J. Combin., 93, 103278, 20pp.
  * [12] M. Fuchs, G.-R. Yu, L. Zhang (2022). Asymptotic enumeration and distributional properties of galled networks, J. Comb. Theory Ser. A., 189, 105599, 28 pages.
  * [13] A. D. M. Gunawan, J. Rathin, L. Zhang (2020). Counting and enumerating galled networks, Discrete Appl. Math., 283, 644–654.
  * [14] C. Holmgren C and S. Janson (2015). Limit laws for functions of fringe trees for binary search trees and recursive trees, Electron. J. Probab., 20, 1–51.
  * [15] G. Kaur, K. P. Choi, T. Wu. Distributions of cherries and pitchforks for the Ford model, arXiv:2110.02850.
  * [16] C. McDiarmid, C. Semple, D. Welsh (2015). Counting phylogenetic networks, Ann. Comb., 19:1, 205–224.
  * [17] A. McKenzie and M. A. Steel (2000). Distributions of cherries for two models of trees, Math. Biosci., 164:1, 81–92.
  * [18] M. Pons and J. Batle (2021). Combinatorial characterization of a certain class of words and a conjectured connection with general subclasses of phylogenetic tree-child networks, Scientific Reports, 11, Article number: 21875.
  * [19] N. A. Rosenberg (2006). The mean and variance of the numbers of $r$-pronged nodes and $r$-caterpillars in Yule generated genealogical trees, Ann. Comb., 10:1, 129–146.
  * [20] T. Wu and K. P. Choi (2016). On joint subtree distributions under two evolutionary models, Theor. Popul. Biol., 108, 13–23.


Current affiliation: ]Quantum Matter Institute, University of British
Columbia, Vancouver, BC V6T 1Z4, Canada

# Electronic and vibrational excitations on the surface of the three-
dimensional topological insulator Bi2Te3-xSex ($x=0,2,3$)

A. Lee<EMAIL_ADDRESS>Department of Physics & Astronomy, Rutgers
University, Piscataway, NJ 08854, USA  H.-H. Kung  [  Department of Physics &
Astronomy, Rutgers University, Piscataway, NJ 08854, USA  Xueyun Wang
Department of Physics & Astronomy, Rutgers University, Piscataway, NJ 08854,
USA  Rutgers Center for Emergent Materials, Rutgers University, Piscataway, NJ
08854, USA  S.-W. Cheong  Department of Physics & Astronomy, Rutgers
University, Piscataway, NJ 08854, USA  Rutgers Center for Emergent Materials,
Rutgers University, Piscataway, NJ 08854, USA  G. Blumberg
<EMAIL_ADDRESS>Department of Physics & Astronomy, Rutgers
University, Piscataway, NJ 08854, USA  National Institute of Chemical Physics
and Biophysics, 12618 Tallinn, Estonia

###### Abstract

We study surface states in the three-dimensional topological insulators
Bi2Te3-xSex ($x=0,2,3$) by polarization resolved resonant Raman spectroscopy.
By tracking the spectral intensity of the surface phonon modes with respect to
the incident photon energy, we show that the surface phonons are qualitatively
similar to their bulk counterparts. Using the resonant Raman excitation
profile, we estimated the energy gaps between the surface conduction bands and
the bulk conduction bands. In addition, we selectively excite the surface-to-
bulk electronic continuum near the Fermi energy in Bi2Se3 to determine the
strength of Fano interaction between the most prominent surface phonon and the
surface-to-bulk continuum.

##  I Introduction

Three-dimensional topological insulators (3DTIs) are materials that contain
gapless, Dirac-like electronic surface states [1, 2, 3] that are topologically
protected against elastic scattering by non-magnetic defects [4] making them
attractive candidates to the application of devices where coherent spin-
polarized transport is desired [5]. Given that the electronic surface states
in 3DTIs are confined to the topmost unit cells [6], complete understanding of
the lattice dynamics interacting with the surface states is critical. This
knowledge is essential not only because the abrupt surface termination in
materials can radically alter the crystal structure and electronic properties
at the surface, but also because inelastic scattering from surface to bulk
electronic states via phonons is expected to be the dominant low-energy
scattering mechanism at room temperature [7]. For these reasons, we
investigated the surface phonons in the 3DTIs Bi2Te3-xSex ($x=0,2,3$) using
polarization-resolved resonant Raman spectroscopy to elucidate their
interaction with the electronic states.

The bismuth-based 3DTIs Bi2Te3-xSex are extensively studied topological
insulator systems due to their stoichiometry, large bulk band gap, and simple
surface spectrum. Bi2Te3-xSex are composed of quintuple layers (QL) of
covalently bonded Bi, Te, and Se atoms stacked along the c-axis [8] with
$D_{3d}$ (R$\bar{3}$m) point group symmetry [8]. The symmetry of the topmost
QL is reduced from the 3D R$\bar{3}$m group to the two-dimensional (2D)
$C_{6v}$ (P6mm) wallpaper group. The QL are weakly bound together by van der
Waals forces resulting in a quasi-2d electronic band structure [9, 10, 11].
The principal bulk band gap is found at the $\Gamma$-point. The lowest bulk
conduction band is denoted as CB0. The highest occupied Dirac-cone, $SS1$, is
composed of electronic surface states, and lies within the principal bulk band
gap, where the bands cross at the $\bar{\Gamma}$-point [12, 10, 11]. There is
also an unoccupied Dirac cone that lies between the second- and third-lowest
bulk conduction bands [13].

In this work, we apply the polarization-resolved resonant Raman excitation
profile (RREP) spectroscopy of the bulk and surface phonon modes in the
bismuth based topological insulators Bi2Te3-xSex to study excitations between
the bulk and surface bands. The surface phonons were identified and
characterized by using polarization-resolved measurements at the incident
photon energies that maximize the surface phonon intensity. We detect the
relative energy gaps between the surface bands by analyzing the intensity of
the phonon modes as a function of incident photon energy. We selectively
excite the lowest energy electron-hole excitation continuum due to transitions
from surface to bulk electronic states that interacts with the surface
phonons. For in Bi2Se3 , by analyzing the interference Fano lineshape of the
surface phonons with the continuum as a function of excitation energy, we
establish the electron-phonon coupling strength between the surface phonons
and the electronic surface states near the Fermi energy. This work also
extends the earlier studies of surface phonons in Bi2Se3 [14] and Bi2Te3 [15]
as it includes resonant Raman enhancement effects, as well as results for
Bi2Te2Se .

##  II Experimental Details

Single crystals of Bi2Te3-xSex were grown by the modified Bridgman technique
[16]. Although the chemical potential in ideal topological insulators rests at
the intersection point of the Dirac cones within the principal bulk band gap,
the presence of natural defects typically raises the chemical potential into
the upper Dirac cone and CB0 [16, 17], thus, the lowest conduction bands for
all samples are assumed to be occupied. The chemical potential of the Bi2Se3
samples grown for this study were determined by scanning tunneling
spectroscopy to cross approximately 150 meV above the Dirac point of the
principal surface band [16] and was later confirmed also by the optical
studies [18, 19]. All samples were cleaved in a N2-rich glove bag and then in-
situ transferred to a continuous He-flow optical cryostat.

Table 1: The Raman selection rules in the bulk and on the surface of Bi2Te3-xSex ($x=0,2,3$). Scattering | Bulk | Surface
---|---|---
Geometry | ($D_{3d}$) | (C6v)
XX |  $A_{1g}$ \+ $E_{g}$ |  $A_{1}$ \+ $E$
XY |  $A_{2g}$ \+ $E_{g}$ |  $A_{2}$ \+ $E$
RR |  $A_{1g}$ \+ $A_{2g}$ |  $A_{1}$ \+ $A_{2}$
RL | $2E_{g}$ | $2E$
Figure 1: Raman response spectra of Bi2Te3-xSex in the XX scattering
polarization geometry taken from (a) Bi2Se3 , (b) Bi2Te2Se , and (c) Bi2Te3 .
(inset) Data from Bi2Se3 plotted against linear scale to illustrate signal-to-
noise ratio. All measurements were performed at 15 K with $\omega_{i}=$ 1.92
eV excitation. The bulk phonons are overlaid with a fit to Lorentzian
lineshape over a local linear background for visual aid. The mode labeled
A${}_{1g}^{*}$ is an antisite defect mode. The secondary emission spectra were
normalized for wavelength dependent scattering volume, the attenuation in the
sample, and reflectivity/transmission of light at the surface interface, see
Eqs. (II) and (2), using optical data from Refs. [20, 21, 22]; the Raman
susceptibility was derived by subtracting the PL signal.

Polarization-resolved Raman scattering measurements were performed in the
quasi-backscattering geometry. For excitation, we use Kr+ ion laser lines with
photon energies ranging from 1.55 to 3.05 eV. All spectra were acquired with
1800 grooves$/$mm diffraction gratings (spectral resolution $\sim 0.1$ meV)
using a triple stage spectrometer setup.

The scattering polarization geometries used in this paper are labeled as
$e_{i}e_{s}$, where $e_{i(s)}$ denotes the polarization of the incident
(collected) light. The incident and collected photons propagate along the
c-axis, with the polarization directed in-plane. The energy of the incident
(scattered) photons are denoted as $\omega_{i(s)}$. The four scattering
polarization geometries used in this paper are XX, YX, RR, and RL. The linear
light polarization directions, X and Y, lie within the ab-plane and are
orthogonal to one another. R and L denote right and left circularly polarized
light, respectively, such that $R(L)=X\pm iY$.

Raman selection rules are listed in Table 1. The bulk Raman excitations have
either $A_{1g}$, $A_{2g}$, or $E_{g}$ symmetry. The irreducible
representations of the Raman-active phonons at the $\Gamma$ point are
$2A_{1g}+2E_{g}$ [24, 25, 26]. In addition, there are four infrared (IR)
active phonons, whose irreducible representations are $2A_{2u}+2E_{u}$ [24].

Table 2:  Summary of the bulk and surface mode energies in Bi2Te2Se and Bi2Te3 . All values are given in units of meV. | Bi2Te2Se |  | Bi2Te3
---|---|---|---
Symmetry | This work | Literature |  | This work | Literature
$A_{1g}^{(1)}$ | 8.1 | 8.4 [27] |  | 7.8 | 7.8 [24]
$A_{1g}^{(2)}$ | 19.1 | 19.6 [27] |  | 17.2 | 16.6 [24]
$E_{g}^{(1)}$ | 4.5 | 4.3 [27] |  | 4.7 | 4.4 [25]
$E_{g}^{(2)}$ | 13.5 | 13.7 [27] |  | 13.1 | 12.8 [24]
$A_{2u}^{(1)}$ | – | 15.9 [28] |  | – | 11.7 [24]
$A_{2u}^{(2)}$ | – | 17.1 [28] |  | – | 14.9 [24]
$E_{u}^{(1)}$ | – | 7.7 [29] |  | – | 6.2 [24] 6.0 [15, 25]
$E_{u}^{(2)}$ | – | 14.5 [27, 30] |  | – | 11.8 [24]
$A_{1}^{(1)}$ | 7.4 | – |  | 7.5 | 7.7 [31]
$A_{1}^{(2)}$ | – | – |  | 16.7 | 17.4 [31]
$A_{1}^{(3)}$ | 15.2 | 14.8 [32] |  | 10.3 | –
$A_{1}^{(4)}$ | 16.7 | – |  | 14.0 | 14.1 [15] 14.8 [31]
$E^{(1)}$ | 4.18 | – |  | – | 5.0 [31]
$E^{(2)}$ | 13.0 | – |  | 12.3 | 12.4 [31]
$E^{(3)}$ | 7.8 | – |  | 7.1 | –
$E^{(4)}$ | 14.9 | – |  | 12.3 | –
Figure 2: Temperature dependence of the energy and linewidth of the
A${}_{1g}^{(2)}$ phonon mode of Bi2Te3-xSex . The energy (red) and linewidth
(black) values of the A${}_{1g}^{(2)}$ mode of (a) Bi2Se3 , (b) Bi2Te2Se , and
(c) Bi2Te3 are displayed in units of meV. The temperature dependence of the
energy and linewidth of the A${}_{1g}^{(2)}$ mode is fit to a 2-phonon
anharmonic decay model [23].

In Fig. 1, we plot on semi-log scale spectra of the Raman susceptibility for
Bi2Te3-xSex . Measurements were performed in the XX scattering geometry at 15
K. The major spectral features are fitted to a Lorentzian function and
overlaid with the experimental data to aid the reader’s eye.

The secondary emission spectra contains both Raman and photoluminescence (PL)
contributions. The measured secondary emission spectra,
$S(\omega_{i},\omega_{s})$, differs from the Raman susceptibility,
$\chi^{{}^{\prime\prime}}(\omega)$, by two factors

$\displaystyle S(\omega_{i},\omega_{s})\sim
C(\omega_{i},\omega_{s})[\\{1+n_{B}(\omega)\\}\chi^{{}^{\prime\prime}}(\omega_{i},\omega)$
$\displaystyle+PL(\omega_{i},\omega_{s})]$ (1)

where $n_{B}(\omega)$ is the Bose-Einstein distribution,
$C(\omega_{i},\omega_{s})$ is a correction factor that accounts for the
scattering volume and power loss due to the optical properties of the
materials, and $PL(\omega_{i},\omega_{s})$ is the contribution from
photoluminescence. The form of $C(\omega_{i},\omega_{s})$ for the
backscattering geometry was derived in Refs. [33, 34]

$C(\omega_{i},\omega_{s})=\frac{\alpha(\omega_{i})+\alpha(\omega_{s})}{T(\omega_{i})T(\omega_{s})}n(\omega_{s})^{2}$
(2)

where $\alpha$ is the optical extinction coefficient, $T$ is the transmission
coefficient at the crystal-vacuum interface, and $n(\omega_{s})$ is the index
of refraction of the material. The values of the optical properties of the
Bi2Te3-xSex crystals are taken from previously published data [20, 21, 22].

The Raman susceptibility spectra was derived from the renormalized secondary
emission spectra by subtracting the PL contribution. The PL contribution to
the secondary emission spectra is approximately a lineary background in the
narrow spectral range of the Raman excitations [19]. We approximated the
background to be the normalized spectral intensity 2.5 meV away from the laser
excitation line, where the Raman continuum is still negligibly weak.

##  III Results and Discussion

The most prominent spectral features in the Bi2Te3-xSex secondary emission
spectra are the Raman-active bulk phonons. The energy of the modes are
consistent with previously reported results, see Table 2. In Bi2Te2Se , there
is an additional broad mode, labeled $A_{1g}^{*}$, that is generally accepted
to be the result of an anti-site defect induced local vibration [32]. It
should be noted that the bulk phonons consistently shifts to higher energies
with Se concentration, see Fig. 3(e).

Figure 3: High resolution $\chi^{{}^{\prime\prime}}(\omega)$ spectra of Bi2Te2Se and Bi2Te3 . $\chi^{{}^{\prime\prime}}(\omega)$ spectra of Bi2Te2Se are taken in the (a) RR and (b) RL scattering geometries. The modes labeled $E^{(*)}$ were previously reported [27], though their origin is unclear. $\chi^{{}^{\prime\prime}}(\omega)$ spectra of Bi2Te3 are acquired in the (c) RR and (d) RL scattering geometries. In all cases, the spectral were acquired at 15 K with 1.92 eV excitation energy, and are plotted in semi-log scale. (e) Phonon energy of Bi2Te3-xSex as a function of Se concentration ($x$).  Table 3: Residual linewidths $\Gamma_{0}$ for the $A_{1g}^{(2)}$ phonon modes. All values are given in units of meV. Material | This work | Literature
---|---|---
Bi2Se3 | 0.08 | 0.06 [35]
Bi2Te2Se | 0.20 | 0.39 [32]
Bi2Te3 | 0.07 | 0.09 [35]

In Fig. 2, we show the evolution of the energy and linewidth of the
$A_{1g}^{(2)}$ mode in Bi2Te3-xSex as a function of temperature. The
$A_{1g}^{(2)}$ mode is used as an example of the other phonon modes, which
display qualitatively similar behavior. The phonons are fit using the Voigt
profile function to deconvolute the effects of the spectrometer resolution
from the phonon linewidths. The energy and linewidth of the mode are fit to
2-phonon anharmonic decay model [23]

$\displaystyle\omega(T)=\omega_{0}+\omega_{R}(2n_{B}(\omega_{0}/2,T)+1)$ (3)
$\displaystyle\Gamma(T)=\Gamma_{0}+\Gamma_{R}(2n_{B}(\omega_{0}/2,T)+1)$ (4)

where $\omega(T)$ and $\Gamma(T)$ are the energy and linewidth of the phonon
at temperature $T$, $\omega_{0}$ is the bare energy, $\Gamma_{0}$ is the
residual linewidth, and $\omega_{R}$ and $\Gamma_{R}$ are fitting parameters.
The $A_{1g}^{(2)}$ mode is well described by the standard 2-phonon anharmonic
decay model, with residual linewidths consistent with previous temperature
dependence measurements indicating good bulk crystal quality, see Table 3.

The discontinuity of the crystal structure at the surface of Bi2Te3-xSex has
several important consequences. Since inversion symmetry is broken at the
surface, the distinction between Raman and IR-active phonons is lifted and
both types of modes are observable using the Raman probe. The optical phonons
at the surface are described by the irreducible representations of the
$\bar{\Gamma}$ point. By making the approximation that the crystal structure
at the surface is comparable with the bulk crystal structure, the surface
phonons are classified according to the corresponding bulk phonons they branch
off from. Under this assumption, the irreducible representations of the
optical surface phonons are $4A_{1}+4E$. In the notation used throughout this
paper, $W^{(1,2)}$ refers to the Raman-related surface modes and $W^{(3,4)}$
refers to the IR-related surface modes, where $W=A_{1}$ or $E$; for example,
the surface mode counterparts of the $A_{2u}^{(1)}$ and $A_{2u}^{(2)}$ modes
are $A_{1}^{(3)}$ and $A_{1}^{(4)}$. A second important consequence is that
the change in the crystal field potential at the surface will be reflected in
the differences between the bulk and surface Raman spectra. For example, the
energy and linewidth of the surface modes with shift relative to their bulk
counterparts, indicating weaker or stronger bond strengths. Thus, information
about the electronic/vibration properties at the surface may be infer from the
surface Raman spectra.

The most distinctive surface related features are the $A_{1}^{(3)}$ and
$A_{1}^{(4)}$ surface phonons. In Bi2Se3 and Bi2Te3 , the $A_{1}^{(4)}$ mode
is the strongest surface mode [14, 15] whereas in Bi2Te2Se the $A_{1}^{(3)}$
mode is the strongest. However, the majority of surface phonon modes are more
than an order of magnitude weaker than their bulk counterparts, making their
observation difficult. In order to aid their observation, we used
polarization-resolved resonant Raman data to selectively probe the different
irreducible representations. The previous Raman study used this method to
identify the surface phonon modes in Bi2Se3 [14], therefore, we restrict our
analysis to the surface phonons in Bi2Te3 and Bi2Te2Se .

Figure 4: Resonant Raman excitation profile (RREP) of the major low energy
spectral features in Bi2Te3-xSex . (a)-(c) Integrated intensity of the major
bulk phonons, and (d)-(f) the $A_{1}^{(3)}$ (Bi2Se3 and Bi2Te3 ) or
$A_{1}^{(4)}$ (Bi2Te2Se ) surface phonon plotted against excitation energy
$\omega_{i}$. The dotted vertical lines are to guide the reader’s eye. (g)-(i)
Raman spectra for Bi2Te3-xSex in the spectral region of the $A_{1}^{(3)}$ and
$A_{1}^{(4)}$ surface phonons with optical corrections, see Eqs. (II) and (2).
The PL signal was approximated to a linear background and subtracted
accordingly. (j) Model of the electronic band structure near the $\Gamma$
point for Bi2Te3-xSex . The lowest energy surface states (SS1), and unoccupied
topological surface states (SS2) are depicted by red lines [36, 37, 13, 38,
18, 19]. The bulk bands are shown in gray. (k) Evolution of the $\Delta_{1}$
and $\Delta_{2}$ energy gaps as a function of Se concentration ($x$).

In Fig. 3(a-b), we show the Raman response function of Bi2Te2Se in the RR and
RL scattering geometries, respectively. The $A_{1}^{(1)}$ mode (7.3 meV) is
red-shifted ($\sim 0.7$ meV) from the $A_{1g}^{(1)}$ mode. The $A_{1}^{(2)}$
mode is obscured by the $A_{1g}^{*}$ feature, therefore cannot be detected by
this method. The $E^{(1)}$ and $E^{(2)}$ modes are not sufficiently red-
shifted from their bulk counterparts to be detected. Although to the best of
our knowledge there is no spectroscopic data of the $A_{2u}$ phonons in
Bi2Te2Se , lattice dynamics calculations of Bi2Te2Se using density functional
theory (DFT) [28] suggest that the $A_{1}$ features at 15.2 meV and 16.7 meV
are the $A_{1}^{(3)}$ and $A_{1}^{(4)}$ modes, respectively. Similarly, the
$E$ features at 7.8 meV and 14.9 meV are the $E^{(3)}$ and $E^{(4)}$ modes,
respectively [29, 30, 21, 27]. In addition, we detected several vibrational
modes with spectral intensities comparable to the surface phonons, labeled
$E^{*}$, that were previously reported from IR spectroscopy measurements [27],
though their origin remains unclear.

In Fig. 3(c-d), we show the Raman spectra of Bi2Te3 in the RR and RL
scattering geometries, respectively. The $A_{1}^{(1)}$ (7.5 meV) and
$A_{1}^{(2)}$ (16.7 meV) modes are red-shifted from their
$A_{1g}$-counterparts by 0.3 meV and 0.5 meV, respectively. The $A_{1}^{(3)}$
(10.3 meV) and $A_{1}^{(4)}$ (14.0 meV) modes red-shifted from their IR-
counterparts by 1.4 meV and 0.9 meV, respectively [24]. Similar to the case of
Bi2Te2Se , the $E^{(1)}$ and $E^{(2)}$ phonons are obscured by their bulk
counterparts. The $E^{(3)}$ and $E^{(4)}$ surface modes appear at 7.1 meV and
12.3 meV, respectively.

The shifts in energy of the surface phonons relative to their bulk
counterparts can be understood in the following way. The small red-shifting of
the $A_{1}$ modes is due to the lack of interlayer van der Waals restorative
force on the vacuum side of the surface layer [39]. The lack of energy
shifting by the $E$-modes suggests that the surface termination does not
significantly alter the in-plane bond dynamics. We may infer, based on the
similar character of the surface phonons to their bulk counter parts, that the
surface lattice dynamics of Bi2Te3-xSex is nearly identical to the bulk. This
is in stark contrast to conventional 3D materials where the surface
termination radically alters the structural and chemical properties of the
top-most atomic layers.

To further investigate the spectral characteristics of the surface and bulk
phonons we collected the data for varying excitation frequency $\omega_{i}$.
The relevant spectral characteristic of interest is the integrated intensity
of a particular mode as a function of $\omega_{i}$, $I_{ph}(\omega_{i})$. The
resonant enhancement factor of Raman active phonons is correlated to
transitions between the initial and intermediate electronic states. Therefore,
we may deduce the evolution of the electronic energy gaps in these materials
by tracking $I_{ph}(\omega_{i})$. Moreover, the Fano anti-resonance effect on
the surface phonons in Bi2Se3 [14, 15] can be studied in greater detail as
both the Raman coupling to the surface phonons and to the electronic continua
is enhanced when incoming photon energy $\omega_{i}$ resonates with the
energies of interband transitions. The surface phonons were fit using the
local fits of the bulk phonons in closest in proximity.

In Fig. 4, we show the resonant Raman excitation profile (RREP) for the major
spectral features in Bi2Te3-xSex and the corresponding electronic band gaps
that can be deduced from the peaks in $I_{ph}(\omega_{i})$. In Figs. 4(a-c),
we show the $I_{ph}(\omega_{i})$ for the $A_{1g}^{(1)}$, $A_{1g}^{(2)}$, and
$E_{g}^{(2)}$ bulk phonon modes from $\omega_{i}=1.65-3.05$ eV; the
$E_{g}^{(1)}$ phonon mode is not included as its spectral intensity was
insufficient for proper analysis. In Fig. 4(d-f), we show the
$I_{ph}(\omega_{i})$ of either the $A_{1}^{(3)}$ or $A_{1}^{(4)}$ surface
modes as the excitation profile of the other mode could not be determined. The
dashed lines between the figures are meant to guide the reader’s eye. In Figs.
4(g-i), we show Raman spectra for several excitations used in the RREP. The
electronic band structure of Bi2Te3-xSex near the $\Gamma$-point is
illustrated in Fig. 4(j) [40], with the corresponding values of the electronic
band gaps $\Delta_{1}$ and $\Delta_{2}$ plotted in Fig. 4(k).

The RREPs of the phonon modes can be understood by considering the evolution
of $\Delta_{1}$ and $\Delta_{2}$ with Se-concentration. The two resonances in
the excitation profile in Bi2Se3 at 1.9 eV and 2.6 eV correspond to the band
gaps between CB0 and the third- and fourth-lowest conduction bands, CB2 and
CB3, respectively [41, 13]. In contrast, Bi2Te2Se and Bi2Te3 only possess a
single peak in their excitation profiles at $\sim 2.2$ eV and $\sim 1.9$ eV,
respectively, corresponding to the $\Delta_{2}$ band gap. This is because the
$\Delta_{1}$ band gap is smaller than 1.6 eV in these materials [13]. These
results are consistent with previous excitation dependent measurements of
Bi2Se3 and Bi2Te3 [15]. However, previous studies only employed 3 laser
energies, making it difficult to find the resonance peak and further correlate
the resonance factor to the band structure [42].

Figure 5: Fano resonance profile of Bi2Se3 (a) Fano resonance of the
A${}_{1}^{(4)}$ surface phonon in Bi2Se3 for different $\omega_{i}$. The data
is shifted by a constant for each spectra for readability. (b) $t_{e}^{2}$ and
$t_{ph}^{2}$ as a function of excitation energy $\omega_{i}$. Values are
derived from the fit to the Fano lineshape. The dotted vertical lines
represent $\Delta_{1}$ and $\Delta_{2}$ energies, see Fig. 4(j).

Finally, we determined the spectral properties of the electronic continuum
that induces the Fano resonance of the $A_{1}^{(4)}$ surface mode of Bi2Se3 .
When an electronic continuum interacts with a discrete phonon state, the Fano
effect induces an asymmetrical spectral lineshape of the phonon [43, 44]. The
degree of asymmetry depends on several factors including the electron-phonon
interaction strength, $v$, the electronic density of states, $\rho(\omega)$,
the phonon-light-interaction strength, $t_{ph}$, and the continuum-light-
interaction strength, $t_{e}$ [44, 45]. Therefore, one may extract information
about the continuum by tracking the phonon lineshape as a function of
excitation energy $\omega_{i}$.

In principle, a surface phonon can couple to any electronic continuum with the
same symmetry where the transitioning electrons interact with the surface
phonon. The “electronic continuum” in this case is the Raman continuum of of
electron-hole excitations that the surface phonon energy coincides with. The
strength of the interaction depends on several factors; in the case of Bi2Se3
only a single continuum can couple to the $A_{1}^{(4)}$ mode. There are three
low-energy electronic continua in Bi2Se3 : the intraband surface continuum of
SS1${}_{-}\rightarrow$ SS1+ transitions (where SS1± refers to the upper
(lower) branches of SS1), the interband surface-to-bulk continua of
SS1${}_{-}\rightarrow$ CB0, and SS1${}_{+}\rightarrow$ CB0 transitions
continuum. Since the chemical potential crosses 150 meV above the Dirac point,
the onset of the intraband surface continuum is roughly 300 meV, which is well
above $A_{1}^{(4)}$ mode energy. The onset of the SS1${}_{-}\rightarrow$ CB0
continuum will be approximately 150 meV, with transitions between the Dirac
point and the CB0 minimum [12]. Thus, the $A_{1}^{(4)}$ mode falls only within
the SS1${}_{+}\rightarrow$ CB0 continuum. Furthermore, if the
SS1${}_{+}\rightarrow$ CB0 continuum is selectively enhanced using resonant
Raman scattering, then the effects of the other two continua are suppressed.

The lineshape of the Raman response for a phonon mode interacting with an
electronic continuum, $\chi^{{}^{\prime\prime}}(\omega)$, is described in the
following way [44, 34, 14, 45]

$\displaystyle\chi^{{}^{\prime\prime}}(\omega)\sim\frac{t_{e}^{2}\pi\rho(\omega_{0}-\omega-
vt_{ph}/t_{e})^{2}}{(\omega-\omega_{0})^{2}+(\gamma-v^{2}\rho\pi)^{2}}$ (5)

where $t_{ph}$ and $t_{e}$ denotes the vertices for the light scatting
processes, $v$ is the electron-phonon interaction, $\omega_{0}$ is the bare
phonon frequency, $\gamma$ is the phonon half-width at half-maximum, and
$\rho$ is the density of the relevant continuum states. In Eq. 5, we have
assumed that $t_{e}$ and $v$ are frequency independent constants for the given
continuum. We have also neglected the real part of the electronic self-energy
and assumed that $v^{2}\rho>>\gamma$ [14].

To determine the spectral properties of the SS1${}_{+}\rightarrow$ CB0
transition, we fit the $A_{1}^{(4)}$ mode Raman spectra taken with different
$\omega_{i}$ using Eq. 5, where $\omega_{0}$, $\gamma$, $v$, and $\rho$ are
assumed to be $\omega_{i}$-independent constants. The factor $t_{e}$ strongly
depends on $\omega_{i}$ as it undergoes enhancement via the resonant Raman
process through SS1$\rightarrow$SS2 dipole transitions.

In Fig. 5(a), we display fits to the A${}_{1}^{(4)}$ surface mode data in
Bi2Se3 for $\omega_{i}=1.55-2.65$ eV excitations. The fits capture the strong
Fano lineshape of the mode in the proximity of $\omega_{i}=1.65$ eV
excitation, the reduction of mode’s asymmetry away from the resonance, and the
excitation dependence of the phonon-light interaction strength. We obtain
electron-phonon interaction strength $v=0.28$ meV, which is similar to the
value from earlier study [14] (0.32 meV).

In Fig. 5(b), we show the evolution of $t_{e}^{2}$ and $t_{ph}^{2}$ with
excitation energy $\omega_{i}$. There are two resonances in the $t_{e}^{2}$
profile: one at $\omega_{i}=$1.65 eV due to SS1-SS2 transitions, and another
at $\omega_{i}=$2.5 eV due to SS1-CB3 transitions. Similarly, in the
$t_{ph}^{2}$ profile there is a weak resonance at 1.83 eV and a stronger
resonance at 2.55 eV. The onset of the SS1-SS2 transitions is slightly below
the $\Delta_{1}$ transition, which is consistent with the model of the band
structure [Fig. 4(j)]. The offset between the lower energy phononic and
electronic resonances is likely because the surface phonon couples stronger to
SS1-CB2 transitions than to SS1-SS2 transitions. The pronounced Fano lineshape
of the A${}_{1}^{(4)}$ surface mode at $\omega_{i}=$1.65 eV is result of the
strong resonance in $t_{e}^{2}$ and the relatively weak resonance in
$t_{ph}^{2}$ at this energy. Similarly, the A${}_{1}^{(4)}$ surface mode
develops a more Lorentzian lineshape at $\omega_{i}=$2.5 eV due to the broader
and weaker resonance in $t_{e}^{2}$ and the strong resonance in $t_{ph}^{2}$
at this energy.

##  IV Conclusions

In conclusion, we have probed the surface phonons and the continuum of
electronic transitions between surface conduction bands in Bi2Te3-xSex using
resonant Raman spectroscopy. By using polarization-resolved resonant Raman
spectroscopy, we observed the surface counterparts of the bulk IR active
phonons and the $A_{1}$ Raman surface phonons in all three materials. The
surface counterparts of the bulk Raman active $E_{g}$ phonons were not
sufficiently red-shifted to be resolved for adequate analysis. We determined
that the energies of all phonon modes consistently harden with increased Se
concentration. In Bi2Se3 , the A${}_{1}^{(4)}$ surface phonon undergoes a
strong Fano anti-resonance when the principal surface bands resonate with the
next-lowest surface bands. This effect is not observed in Bi2Te2Se and Bi2Te3
as the band gap in these materials is below the lowest laser energy of our
setup. The quantitative similarities between the resonance profiles of the
bulk and surface phonon modes is good evidence that the surface modes in the
materials are related to the corresponding bulk phonons with renormalized
energies due to their surface discontinuity.

##  V Acknowledgements

###### Acknowledgements.

We thank R. Merlin and I. Boulares for discussions. The spectroscopic work
(A.C.L., H.-H.K. and G.B.) was supported by the NSF under Grant No.
DMR-2105001. The sample growth (X.W. and S.-W.C.) were supported by the DOE
under Grant No. DE-FG02-07ER46382. The work at NICPB was supported by the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme Grant Agreement No. 885413.

## References

  * Moore [2010] J. E. Moore, The birth of topological insulators, Nature 464, 194 (2010).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Fu  _et al._ [2007] L. Fu, C. L. Kane, and E. J. Mele, Topological insulators in three dimensions, Phys. Rev. Lett. 98, 106803 (2007).
  * Roushan  _et al._ [2009] P. Roushan, J. Seo, C. V. Parker, Y. S. Hor, D. Hsieh, D. Qian, A. Richardella, M. Z. Hasan, R. J. Cava, and A. Yazdani, Topological surface states protected from backscattering by chiral spin texture, Nature 460, 1106 (2009).
  * Žutić  _et al._ [2004] I. Žutić, J. Fabian, and S. Das Sarma, Spintronics: Fundamentals and applications, Rev. Mod. Phys. 76, 323 (2004).
  * Linder  _et al._ [2009] J. Linder, T. Yokoyama, and A. Sudbø, Anomalous finite size effects on surface states in the topological insulator ${\text{Bi}}_{2}{\text{Se}}_{3}$, Phys. Rev. B 80, 205401 (2009).
  * Park  _et al._ [2010] S. R. Park, W. S. Jung, C. Kim, D. J. Song, C. Kim, S. Kimura, K. D. Lee, and N. Hur, Quasiparticle scattering and the protected nature of the topological states in a parent topological insulator ${\text{Bi}}_{2}{\text{Se}}_{3}$, Phys. Rev. B 81, 041405 (2010).
  * Nakajima [1963] S. Nakajima, The crystal structure of Bi2Te3-xSex, Journal of Physics and Chemistry of Solids 24, 479 (1963).
  * Ribak  _et al._ [2016] A. Ribak, K. B. Chashka, E. Lahoud, M. Naamneh, S. Rinott, Y. Ein-Eli, N. C. Plumb, M. Shi, E. Rienks, and A. Kanigel, Internal pressure in superconducting Cu-intercalated ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. B 93, 064505 (2016).
  * Arakane  _et al._ [2012] T. Arakane, T. Sato, S. Souma, K. Kosaka, K. Nakayama, M. Komatsu, T. Takahashi, Z. Ren, K. Segawa, and Y. Ando, Tunable Dirac cone in the topological insulator Bi2-xSbxTe3-ySey, Nature communications 3, 636 (2012).
  * Chen  _et al._ [2009] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, Experimental Realization of a Three-Dimensional Topological Insulator, Bi2Te3, Science 325, 178 (2009).
  * Xia  _et al._ [2009] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, _et al._ , Observation of a large-gap topological-insulator class with a single Dirac cone on the surface, Nature physics 5, 398 (2009).
  * Niesner  _et al._ [2012] D. Niesner, T. Fauster, S. V. Eremeev, T. V. Menshchikova, Y. M. Koroteev, A. P. Protogenov, E. V. Chulkov, O. E. Tereshchenko, K. A. Kokh, O. Alekperov, A. Nadjafov, and N. Mamedov, Unoccupied topological states on bismuth chalcogenides, Phys. Rev. B 86, 205403 (2012).
  * Kung  _et al._ [2017a] H.-H. Kung, M. Salehi, I. Boulares, A. F. Kemper, N. Koirala, M. Brahlek, P. Lošťák, C. Uher, R. Merlin, X. Wang, S.-W. Cheong, S. Oh, and G. Blumberg, Surface vibrational modes of the topological insulator ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$ observed by Raman spectroscopy, Phys. Rev. B 95, 245406 (2017a).
  * Boulares  _et al._ [2018] I. Boulares, G. Shi, E. Kioupakis, P. Lošt’ák, C. Uher, and R. Merlin, Surface phonons in the topological insulators Bi2Se3 and Bi2Te3, Solid State Communications 271, 1 (2018).
  * Dai  _et al._ [2016] J. Dai, D. West, X. Wang, Y. Wang, D. Kwok, S.-W. Cheong, S. B. Zhang, and W. Wu, Toward the Intrinsic Limit of the Topological Insulator ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. Lett. 117, 106401 (2016).
  * Wang  _et al._ [2013] L.-L. Wang, M. Huang, S. Thimmaiah, A. Alam, S. L. Bud’ko, A. Kaminski, T. A. Lograsso, P. Canfield, and D. D. Johnson, Native defects in tetradymite Bi2(TexSe3-x) topological insulators, Phys. Rev. B 87, 125303 (2013).
  * Kung  _et al._ [2017b] H.-H. Kung, S. Maiti, X. Wang, S.-W. Cheong, D. L. Maslov, and G. Blumberg, Chiral Spin Mode on the Surface of a Topological Insulator, Phys. Rev. Lett. 119, 136802 (2017b).
  * Kung  _et al._ [2019] H.-H. Kung, A. P. Goyal, D. L. Maslov, X. Wang, A. Lee, A. F. Kemper, S.-W. Cheong, and G. Blumberg, Observation of chiral surface excitons in a topological insulator Bi2Se3, Proceedings of the National Academy of Sciences 116, 4006 (2019).
  * McIver  _et al._ [2012] J. W. McIver, D. Hsieh, S. G. Drapcho, D. H. Torchinsky, D. R. Gardner, Y. S. Lee, and N. Gedik, Theoretical and experimental study of second harmonic generation from the surface of the topological insulator Bi2Se3, Phys. Rev. B 86, 035327 (2012).
  * Aleshchenko  _et al._ [2014] Y. A. Aleshchenko, A. V. Muratov, V. V. Pavlova, Y. G. Selivanov, and E. G. Chizhevskii, Infrared spectroscopy of Bi2Te2Se, JETP Letters 99, 187 (2014).
  * Cui  _et al._ [1999] H. Cui, I. Bhat, and R. Venkatasubramanian, Optical constants of Bi2Te3 and Sb2Te3 measured using spectroscopic ellipsometry, Journal of Electronic Materials 28, 10.1007/s11664-999-0247-z (1999).
  * Klemens [1966] P. G. Klemens, Anharmonic decay of optical phonons, Phys. Rev. 148, 845 (1966).
  * Richter and Becker [1977] W. Richter and C. Becker, A Raman and far-infrared investigation of phonons in the rhombohedral V2–VI3 compounds Bi2Te3, Bi2Se3, Sb2Te3 and Bi2(Te1-xSex)3 (0$<$ x$<$ 1),(Bi1-ySby)2Te3 (0$<$ y$<$ 1), physica status solidi (b) 84, 619 (1977).
  * Wang and Zhang [2012] B.-T. Wang and P. Zhang, Phonon spectrum and bonding properties of Bi2Se3: Role of strong spin-orbit interaction, Applied Physics Letters 100, 082109 (2012).
  * Cheng and Ren [2011] W. Cheng and S.-F. Ren, Phonons of single quintuple Bi2Te3 and Bi2Se3 films and bulk materials, Phys. Rev. B 83, 094301 (2011).
  * Akrap  _et al._ [2012] A. Akrap, M. Tran, A. Ubaldini, J. Teyssier, E. Giannini, D. van der Marel, P. Lerch, and C. C. Homes, Optical properties of Bi2Te2Se at ambient and high pressures, Phys. Rev. B 86, 235207 (2012).
  * Shi  _et al._ [2015] H. Shi, D. Parker, M.-H. Du, and D. J. Singh, Connecting Thermoelectric Performance and Topological-Insulator Behavior: ${\mathrm{Bi}}_{2}{\mathrm{Te}}_{3}$ and ${\mathrm{Bi}}_{2}{\mathrm{Te}}_{2}\mathrm{Se}$ from First Principles, Phys. Rev. Applied 3, 014004 (2015).
  * Reijnders  _et al._ [2014] A. A. Reijnders, Y. Tian, L. J. Sandilands, G. Pohl, I. D. Kivlichan, S. Y. F. Zhao, S. Jia, M. E. Charles, R. J. Cava, N. Alidoust, S. Xu, M. Neupane, M. Z. Hasan, X. Wang, S. W. Cheong, and K. S. Burch, Optical evidence of surface state suppression in Bi-based topological insulators, Phys. Rev. B 89, 075138 (2014).
  * Di Pietro  _et al._ [2012] P. Di Pietro, F. M. Vitucci, D. Nicoletti, L. Baldassarre, P. Calvani, R. Cava, Y. S. Hor, U. Schade, and S. Lupi, Optical conductivity of bismuth-based topological insulators, Phys. Rev. B 86, 045439 (2012).
  * Ren  _et al._ [2012] L. Ren, X. Qi, Y. Liu, G. Hao, Z. Huang, X. Zou, L. Yang, J. Li, and J. Zhong, Large-scale production of ultrathin topological insulator bismuth telluride nanosheets by a hydrothermal intercalation and exfoliation route, J. Mater. Chem. 22, 4921 (2012).
  * Tian  _et al._ [2016] Y. Tian, G. B. Osterhoudt, S. Jia, R. J. Cava, and K. S. Burch, Local phonon mode in thermoelectric Bi2Te2Se from charge neutral antisites, Applied Physics Letters 108, 041911 (2016).
  * Loudon [1965] R. Loudon, Theory of the resonance raman effect in crystals, Journal de Physique 26, 677 (1965).
  * Blumberg  _et al._ [1994] G. Blumberg, M. Klein, L. Börjesson, R. Liang, and W. Hardy, Investigation of the temperature dependence of electron and phonon Raman scattering in single crystal YBa2Cu3O6.952, Journal of superconductivity 7, 445 (1994).
  * Tian  _et al._ [2017] Y. Tian, S. Jia, R. J. Cava, R. Zhong, J. Schneeloch, G. Gu, and K. S. Burch, Understanding the evolution of anomalous anharmonicity in ${\mathrm{Bi}}_{2}{\mathrm{Te}}_{3-x}{\mathrm{Se}}_{x}$, Phys. Rev. B 95, 094104 (2017).
  * Soifer  _et al._ [2019] H. Soifer, A. Gauthier, A. F. Kemper, C. R. Rotundu, S.-L. Yang, H. Xiong, D. Lu, M. Hashimoto, P. S. Kirchmann, J. A. Sobota, and Z.-X. Shen, Band-resolved imaging of photocurrent in a topological insulator, Phys. Rev. Lett. 122, 167401 (2019).
  * Eremeev  _et al._ [2013] S. Eremeev, I. Silkin, T. Menshchikova, A. P. Protogenov, and E. V. Chulkov, New topological surface state in layered topological insulators: unoccupied Dirac cone, JETP letters 96, 780 (2013).
  * Nurmamat  _et al._ [2013] M. Nurmamat, E. E. Krasovskii, K. Kuroda, M. Ye, K. Miyamoto, M. Nakatake, T. Okuda, H. Namatame, M. Taniguchi, E. V. Chulkov, K. A. Kokh, O. E. Tereshchenko, and A. Kimura, Unoccupied topological surface state in Bi2Te2Se, Phys. Rev. B 88, 081301 (2013).
  * Zhang  _et al._ [2011] J. Zhang, Z. Peng, A. Soni, Y. Zhao, Y. Xiong, B. Peng, J. Wang, M. S. Dresselhaus, and Q. Xiong, Raman spectroscopy of few-quintuple layer topological insulator Bi2Se3 nanoplatelets, Nano letters 11, 2407 (2011).
  * Zhang  _et al._ [2009] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Topological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surface, Nature physics 5, 438 (2009).
  * Sobota  _et al._ [2013] J. A. Sobota, S.-L. Yang, A. F. Kemper, J. J. Lee, F. T. Schmitt, W. Li, R. G. Moore, J. G. Analytis, I. R. Fisher, P. S. Kirchmann, T. P. Devereaux, and Z.-X. Shen, Direct Optical Coupling to an Unoccupied Dirac Surface State in the Topological Insulator ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. Lett. 111, 136802 (2013).
  * Weber and Merlin [2000] W. H. Weber and R. Merlin,  _Raman scattering in materials science_, Vol. 42 (Springer Science & Business Media, 2000) Chap. 1.2.
  * Fano [1961] U. Fano, Effects of configuration interaction on intensities and phase shifts, Phys. Rev. 124, 1866 (1961).
  * [44] M. edited by Cardona,  _Electronic Raman scattering, in Light Scattering in Solids I_ (Springer-Verlag, Berlin, 1983) pp. 147–204.
  * Wu  _et al._ [2020] S.-F. Wu, W.-L. Zhang, L. Li, H.-B. Cao, H.-H. Kung, A. S. Sefat, H. Ding, P. Richard, and G. Blumberg, Coupling of fully symmetric As phonon to magnetism in $\mathrm{Ba}{({\mathrm{Fe}}_{1-x}{\mathrm{Au}}_{x})}_{2}{\mathrm{As}}_{2}$, Phys. Rev. B 102, 014501 (2020).