# Note on Disjoint Cycles in Multipartite Tournaments

Gregory Gutin Department of Computer Science. Royal Holloway University of
London<EMAIL_ADDRESS> Wei Li School of Mathematics and Statistics,
Northwestern Polytechnical University<EMAIL_ADDRESS> Shujing Wang
Corresponding author. School of Mathematics and Statistics, Central China
Normal University<EMAIL_ADDRESS> Anders Yeo Department of Mathematics
and Computer Science, University of Southern Denmark and Department of
Mathematics and Applied Mathematics, University of Johannesburg.
<EMAIL_ADDRESS> Yacong Zhou Department of Computer Science. Royal Holloway
University of London<EMAIL_ADDRESS>
###### Abstract

In 1981, Bermond and Thomassen conjectured that for any positive integer $k$,
every digraph with minimum out-degree at least $2k-1$ admits $k$ vertex-
disjoint directed cycles. In this short paper, we verify the Bermond-Thomassen
conjecture for triangle-free multipartite tournaments and 3-partite
tournaments. Furthermore, we characterize 3-partite tournaments with minimum
out-degree at least $2k-1$ ($k\geq 2$) such that in each set of $k$ vertex-
disjoint directed cycles, every cycle has the same length.

##  1 Introduction

In this paper, digraphs do not have loops or parallel arcs.111Undefined
terminology can be found later in this section or in the monograph [4]. For a
digraph $D$, $V(D)$ ($A(D)$, resp.) are the vertex set (arc set, resp.) of
$D$. A digraph $D$ is complete if $xy$ and $yx$ are arcs for every pair $x,y$
of distinct vertices, i.e., $|A(D)|=|V(D)|(|V(D)|-1).$ The minimum out-degree
of a vertex in $D$ is denoted by $\delta^{+}(D).$ An orientation of an
undirected graph $G$ is a digraph obtained by replacing each edge $\\{x,y\\}$
of $G$ with arc $xy$ or $yx$. A multipartite tournament is an orientation of a
complete multipartite graph. A multipartite tournament with $t\geq 2$ partite
sets is a $t$-partite tournament. Note that a tournament is a multipartite
tournament in which every partite set is a singleton.

Bermond and Thomassen [5] observed that a complete digraph $K^{*}_{n}$ on
$2k-1$ vertices has $k-1$ vertex-disjoint cycles, but no $k$ such cycles. Note
that $\delta^{+}(K^{*}_{n})=2k-2$. Let $f(k)$ be the minimum integer such that
if $\delta^{+}(D)\geq f(k)$ for a digraph $D$ then $D$ has $k$ vertex-disjoint
cycles (if such an $f(k)$ does not exist then we set $f(k)=\infty$). Bermond
and Thomassen [5] posed the following:

######  Conjecture 1.1.

For every $k\in\mathbb{Z}_{>0}$, $f(k)=2k-1.$

This conjecture is trivially true for $k=1$. Thomassen [10] showed it for
$k=2.$ The conjecture was proved for $k=3$ by Lichiardopol et al. [9]; a
shorter proof of this case is given by Bai et al. [2]. For every
$k\in\mathbb{Z}_{>0}$, Thomassen [10] proved that $f(k)\leq(k+1)!$. The upper
bound on $f(k)$ was drastically improved to $64k$ by Alon [1] and improved
further to $18k$ by Bucić [6] who combined the probabilistic approach of Alon
[1] and the result of Lichiardopol et al. [9].

In this paper, all cycles are directed. A cycle $C$ in a digraph $D$ is a
$p$-cycle if $C$ has $p$ arcs. A 3-cycle is called a triangle. A digraph is
triangle-free if it has no triangles. Conjecture 1.1 was proved for
tournaments by Bang-Jensen et al. [3] as follows.

######  Theorem 1.2.

Every tournament with $\delta^{+}(T)\geq 2k-1$ has $k$ vertex-disjoint
triangles.

A natural question one would ask is whether we can prove this conjecture for
some generalization of tournaments. The first generalization one may try is
extended tournaments, where we blow up every vertex of a tournament to an
arbitrary independent set. However, it is easy to see that the assertion that
the conjecture holds for extended tournaments is a simple corollary of Theorem
1.2. In fact, we can replace every partite set of an extended tournament with
an arbitrary transitive tournament on the same vertices and apply Theorem 1.2
to the obtained tournament. Then, we are done by the fact that every triangle
in the tournament is also a triangle in the extended tournament.

Extended tournaments are a special class of multipartite tournaments. Thus,
let us consider multipartite tournaments. Bai et al. [2] proved Conjecture 1.1
for bipartite tournaments.

######  Theorem 1.3.

Let $D$ be a bipartite tournament. If $\delta^{+}(D)\geq 2k-1$ then $D$ has at
least $k$ vertex-disjoint $4$-cycles.

In this paper, we prove that (i) a triangle-free multipartite tournament $D$
with $\delta^{+}(D)\geq 2k-1$ contains $k$ vertex-disjoint cycles, and (ii) a
$3$-partite tournament $D$ with $\delta^{+}(D)\geq 2k-1$ contains $k$ vertex-
disjoint cycles.

Chen and Chang [7] extended Theorems 1.2 and 1.3 where diversity of the
lengths of the cycles is taken into consideration. To state their results, we
will use the following notation. Let $\mathcal{C}$ be a set of $k$ vertex-
disjoint cycles of $D$. Let $\kappa^{k}(\cal C)$ be the maximum number of
cycles in $\cal C$ of distinct lengths. Let

$\kappa^{k}(D)=\max\\{\kappa^{k}(\mathcal{C}):\mathcal{C}\textit{~{}is a set
of k vertex-disjoint directed cycles in}~{}D\\}.$

For tournaments, Chen and Chang [7] proved the following:

######  Theorem 1.4.

Let $D$ be a tournament with $\delta^{+}(D)\geq 2k-1$ and let $k\geq 3$. Then
$\kappa^{k}(D)\geq 2$.

For two subsets $X,Y$ of $V(D)$, let $A_{D}(X,Y)$ be the set of arcs from $X$
to $Y$ (we will omit the subscript if it is clear from the context). For an
integer $k\geq 2$ and positive integers $n_{1},n_{2},\dots,n_{2k}$, a digraph
$D$ belongs to the digraph set $BT(n_{1},n_{2},\dots,n_{2k})$ if $D$ is a
bipartite tournament with partite sets $X,Y$ such that (i) $X$ is the disjoint
union of $X_{1},X_{2},\dots,X_{2k}$, where $|X_{i}|=n_{i}$ for every
$i\in[2k]$, and $Y$ is the disjoint union of $Y_{1},Y_{2},\dots,Y_{2k}$, where
$|Y_{i}|=1$ for every $i\in[2k]$, and (ii) the arc set of $D$ is

$(\cup_{1\leq i\neq j\leq 2k}A(X_{i},Y_{j}))\cup(\cup_{1\leq i\leq
2k}A(Y_{i},X_{i})).$

Chen and Chang [7] obtained the following characterization of strong bipartite
tournaments $D$ with $\kappa^{k}(D)=1$.

######  Theorem 1.5.

Let $D$ be a strong bipartite tournament with $\delta^{+}(D)\geq 2k-1,$ where
$k\geq 2$. If $\kappa^{k}(D)=1$, then $D$ is a member of
$BT(n_{1},n_{2},\ldots,n_{2k})$, where $n_{i}\geq 2k-1$ for $i=1,\ldots,k$.

In this paper, we obtain a characterization of 3-partite tournaments $D$ with
$\delta^{+}(D)\geq 2k-1$ and $\kappa^{k}(D)=1$.

Additional Terminology and Notation Sometimes, if $xy\in A(D)$, we will write
$x\rightarrow y$ when $D$ is known from the context. Let $C$ be a cycle and
let $x,y$ be two vertices in the cycle. We write $C[x,y]$ to denote the
$(x,y)$-path on $C$. Let $a_{D}(X,Y)=|A_{D}(X,Y)|$.

For $x\in V(D)$, $N^{+}(x)$ ($N^{-}(x)$, resp.) denote the out-neighbourhood
(in-neighbourhood, resp.) of $x$. We write $d^{+}(x)$ ($d^{-}(x)$, resp.) to
denote the out-degree (in-degree, resp.) of $x$, i.e., $d^{+}(x)=|N^{+}(x)|$
($d^{-}(x)=|N^{-}(x)|$, resp.). A vertex $v$ is a sink if $d^{+}(x)=0$. For a
subset $V^{\prime}\subseteq V(D)$, we use $d^{+}_{V^{\prime}}(x)$ to denote
the number of out-neighbours of $x$ in $V^{\prime}$, i.e.,
$d^{+}_{V^{\prime}}(x)=|N^{+}(x)\cap V^{\prime}|$.

##  2 Proofs

We need to use the following vertex-pancyclic-like result for multipartite
tournaments obtained by Guo and Volkmann [8]. Initially, we obtained a long
technical proof of Lemma 2.2 (and therefore Theorem 2.3) without using this
result. However, with the help of this result, we can now give a short proof
for this theorem.

######  Theorem 2.1.

[8] If $D$ is a strong $t$-partite $(t\geq 3)$ tournament, then every partite
set of $D$ has at least one vertex, which lies on an $m$-cycle for all
$m\in\\{3,4,\dots,t\\}$.

######  Lemma 2.2.

If $D$ is a triangle-free multipartite tournament with $\delta^{+}(D)\geq
2k-1$, then $D$ has $k$ vertex-disjoint cycles.

###### Proof.

Let $H$ be the terminal strong component of $D$. Clearly, $\delta^{+}(H)\geq
2k-1$. We claim that $H$ is a bipartite tournament. Otherwise, $H$ is a strong
$t$-partite tournament with $t\geq 3$, and therefore, by Theorem 2.1, $H$ has
a triangle, which contradicts the fact that $D$ is triangle-free. Thus, we are
done by applying Theorem 1.3 to $H$. ∎

######  Theorem 2.3.

Let $D$ be a $3$-partite tournament and $\delta^{+}(D)\geq 2k-1$. Then $D$
contains $k$ vertex-disjoint cycles.

###### Proof.

We prove this by induction on $k$. When $k=1$, this result clearly holds.
Assume that $k\geq 2$ and this theorem is true for $k-1$. If $D$ has a
triangle $C$, then one can observe that $\delta^{+}(D-V(C))\geq
2k-1-2=2(k-1)-1$. By the induction hypothesis, $D-V(C)$ has $k-1$ vertex-
disjoint cycles. Therefore, together with $C$, $D$ has at least $k$ vertex-
disjoint cycles. If $D$ is triangle-free, then we are done by Lemma 2.2. ∎

By a similar induction, we can prove the following result for all multipartite
tournaments.

######  Corollary 2.4.

Let $D$ be a multipartite tournament and $\delta^{+}(D)\geq 3k-2$. Then, $D$
contains $k$ vertex-disjoint cycles.

###### Proof.

We prove this by induction on $k$. When $k=1$, this result clearly holds.
Assume that $k\geq 2$ and it is true for $k-1$. We may assume that $D$ has a
triangle $C$ (for otherwise we are done by Lemma 2.2). Then,
$\delta^{+}(D-V(C))\geq 3k-2-3=3(k-1)-2$. And we are done by using the
induction hypothesis. ∎

######  Observation 2.5.

Let $D$ be a multipartite tournament and $S$ be its sink set, then all sinks
lie in the same partite set and every vertex in $V(D)\setminus S$ can reach
all vertices in $S$ by a path with length at most two.

######  Lemma 2.6.

Let $D$ be a multipartite tournament and let $C$, $C^{\prime}$ be two vertex-
disjoint $3$-cycles in $D$. If $a(V(C),V(C^{\prime}))>0$ and
$a(V(C^{\prime}),V(C))>0$, then $D[V(C)\cup V(C^{\prime})]$ contains a cycle
with length at least $4$.

###### Proof.

Let $V(C)=\\{x,y,z\\}$ and
$V(C^{\prime})=\\{x^{\prime},y^{\prime},z^{\prime}\\}$. As $D$ is a
multipartite tournament, $a(V(C),V(C^{\prime}))+a(V(C^{\prime}),V(C))\geq 6$.
Thus, we may assume without loss of generality that $a(V(C),V(C^{\prime}))\geq
3$ and $y^{\prime}\to y$. Since $a(V(C),V(C^{\prime}))\geq 3$, there is an arc
from $\\{x,z\\}$ to $\\{x^{\prime},z^{\prime}\\}$. Assume w.l.o.g. that $x\to
x^{\prime}$. Then,
$xx^{\prime}C^{\prime}[x^{\prime},y^{\prime}]y^{\prime}yC[y,x]x$ is a cycle
with length at least 4. ∎

######  Lemma 2.7.

Let $D$ be a $3$-partite tournament and $C_{1}$, $C_{2}$ and $C_{3}$ be three
vertex-disjoint $3$-cycles in $D$. If for all $i\in[3]$,
$a(V(C_{i+1}),V(C_{i}))=0$ where $C_{4}=C_{1}$, then $D[V(C_{1})\cup
V(C_{2})\cup V(C_{3})]$ has two vertex-disjoint cycles with different lengths.

###### Proof.

Let $V(C_{i})=\\{x_{i},y_{i},z_{i}\\}$. Assume w.l.o.g. that $x_{i}$’s
($y_{i}$’s, $z_{i}$’s, respectively) are in the same partite set and $x_{1}\to
y_{1}\to z_{1}\to x_{1}$ (for otherwise we can exchange $y_{i}$’s and
$z_{i}$’s). Thus, we have a 4-cycle $x_{1}y_{1}z_{2}y_{3}x_{1}$ and a 3-cycle
$z_{1}y_{2}x_{3}z_{1}$. ∎

######  Theorem 2.8.

Every $3$-partite tournament $D$ with $\delta^{+}(D)\geq 2k-1$ satisfies that
$\kappa^{k}(D)\geq 2$ unless $D$ is triangle-free.

###### Proof.

Denote the three partite sets of $D$ by $V_{1}$, $V_{2}$, and $V_{3}$. Assume
to the contrary that $D$ contains a triangle and $\kappa^{k}(D)=1$. By the
proof of Theorem 2.3, it suffice to consider the case when we have collection
of $k$ vertex-disjoint 3-cycles $\mathcal{C}=\\{C_{1},\ldots,C_{k}\\}$ where
$C_{i}=x_{i}y_{i}z_{i}x_{i}$ for all $i\in[k]$. Let
$V(\mathcal{C})=\cup_{i=1}^{k}V(C_{i})$, $U=D-V(\mathcal{C})$ and $S$ the set
of sinks in $U$ (i.e., $S=\\{x\in V(U):d^{+}_{V(U)}(x)=0\\}$). We will
consider the following two cases.

Case 1: $S\neq\emptyset$.

By Observation 2.5, we may assume without loss of generality that $S\subseteq
V_{3}$ and $z_{i}\in V_{3}$ for all $i\in[k]$. Note that for every $s\in S$,
we have that $d^{+}_{V(\cal C)}(s)=d^{+}(s)\geq 2k-1$. Since $(N^{+}(s)\cup
N^{-}(s))\cap V(\mathcal{C})\subseteq V(\mathcal{C})\cap(V_{1}\cup V_{2})$,
$d_{V(\mathcal{C})}^{+}(s)+d_{V(\mathcal{C})}^{-}(s)=2k$. Hence, for any $s\in
S$, we have that

$d^{-}_{V(\mathcal{C})}(s)\leq 1.$ (1)

Claim 1.A. For any $i\in[k]$, $d^{+}_{V(U)}(z_{i})=0$ and
$N^{+}(z_{i})=\\{x_{1},y_{1},\ldots,x_{k},y_{k}\\}\setminus\\{y_{i}\\}$.

###### Proof.

Assume that there exists $u\in V(U)$ such that $z_{i}u\in A(D)$. As $z_{i}\in
V_{3}$, we have that $u\notin V_{3}$ and therefore $us\in A(D)$ for all $s\in
S$. Note that by (1), we have that either $sx_{i}\in A(D)$ or $sy_{i}\in
A(D)$. Thus, either $z_{i}usx_{i}y_{i}z_{i}$ or $z_{i}usy_{i}z_{i}$ is a cycle
of $D$. This cycle together with cycles in $\mathcal{C}\setminus\\{C_{i}\\}$
forms a collection of $k$ vertex-disjoint cycles that have $k-1$ $3$-cycles
and one cycle of length at least 4, a contradiction. Thus,
$N^{+}(z_{i})\subseteq V({\mathcal{C}})$. Since $d^{+}(z_{i})\geq 2k-1$, we
have $N^{+}(z_{i})=\\{x_{1},y_{1},\ldots,x_{k},y_{k}\\}\setminus\\{y_{i}\\}$.
∎

Claim 1.B. For $i\in[k],S\subseteq N^{+}(y_{i})$.

###### Proof.

For any $i\in[k]$, by Claim 1.A, we have that
$d^{+}_{V(\mathcal{C})}(x_{i})\leq k$ and therefore $d^{+}_{V(U)}(x_{i})\geq
k-1>0$, which implies that $x_{i}$ has an out-neighbour in $V(U)$. Let $u\in
V(U)$ be the out-neighbour of $x_{i}$. If $u\in S$, then, by (1), we have that
$uy_{i}\in A(D)$. Thus we have that $C_{i}^{\prime}:=x_{i}uy_{i}z_{i}x_{i}$ is
a cycle of length 4 and
$\mathcal{C}^{\prime}:=(\mathcal{C}\setminus\\{C_{i}\\})\cup\\{C_{i}^{\prime}\\}$
is a collection of $k$ vertex-disjoint cycles with
$\kappa^{k}(\mathcal{C}^{\prime})=2$, a contradiction. Thus, $u\in
V(U)\setminus S$.

Now, we assume to the contrary that there is a vertex $s\in S$ such that
$sy_{i}\in A(D)$. By observation 2.5, there is an $(u,s)$-path with length at
most 2 in $U$. The existence of this path together with $s\to y_{i}\to
z_{i}\to x_{i}\to u$ implies a cycle with length at least 4 in $D[V(U)\cup
V(C_{i})]$ which together with cycles in $\mathcal{C}\setminus\\{C_{i}\\}$
implies that $\kappa^{k}(D)\geq 2$, a contradiction. ∎

By Claim 1.B, we have that $d^{-}_{V(\mathcal{C})}(s)\geq k\geq 2$ for all
$s\in S$, a contradiction to (1). This completes the proof of Case 1.

Case 2. $S=\emptyset$.

Let $D^{\prime}$ be the terminal strong component of $U$. Note that
$\delta^{+}(D^{\prime})=1$ for otherwise $\delta^{+}(D^{\prime})\geq 2$ and
therefore $D^{\prime}$ contains a cycle of length at least 4. Let
$C_{k+1}=x_{k+1}y_{k+1}z_{k+1}x_{k+1}$ be a cycle in $D^{\prime}$ that
includes a vertex of the minimum out-degree. Assume without loss of generality
that $d^{+}_{V(U)}(x_{k+1})=1$.
Claim 2.A. For $i,j\in[k+1]$, all arcs between $V(C_{i})$ and $V(C_{j})$ go in
the same direction.

###### Proof.

Assume to the contrary that $a(V(C_{i}),V(C_{j}))>0$ and
$a(V(C_{j}),V(C_{i}))>0$. By Lemma 2.6, there is a cycle of length at least 4
in $D[V(C_{i})\cup V(C_{j})]$ and therefore we are done by the fact that this
new cycle together with cycles
$(\mathcal{C}\cup\\{C_{k+1}\\})\setminus\\{C_{i},C_{j}\\}$ forms a set of $k$
cycles that have $k-1$ 3-cycles and one cycle with length at least 4, a
contradiction. ∎

Construct a digraph $T$ with $k+1$ vertices $c_{1}$, $\dots$, $c_{k+1}$ where
$c_{i}\to c_{j}$ if and only if $a(V(C_{i}),V(C_{j}))>0$. By Claim 2.A, $T$ is
a tournament.
Claim 2.B. $T$ is a transitive tournament.

###### Proof.

Suppose it is not. Then by Theorem 2.1, $T$ has a 3-cycle
$c_{i}c_{j}c_{h}c_{i}$, which, by Claim 2.A, implies that
$a(V(C_{j}),V(C_{i}))=a(V(C_{h}),V(C_{j}))=a(V(C_{i}),V(C_{h}))=0$. By Lemma
2.7, there are two vertex-disjoint cycles with different lengths in
$D[V(C_{i})\cup V(C_{j})\cup V(C_{h})]$. These two cycles, together with
cycles in $(\mathcal{C}\cup\\{C_{k+1}\\})\setminus\\{C_{i},C_{j},C_{h}\\}$,
forms $k$ cycles with at least two of them having different lengths, a
contradiction. ∎

By Claim 2.B, we assume w.l.o.g. that $c_{1}$, $c_{2}$, $\dots$, $c_{k}$ is an
acyclic ordering in $T-c_{k+1}$. Since $d^{+}_{V(U)}(x_{k+1})=1$,
$d^{+}_{V(\mathcal{C})}(x_{k+1})\geq 2(k-1)$. Thus, $d^{+}_{T}(c_{k+1})\geq
k-1$ and therefore either $c_{k+1}$, $c_{1}$, $c_{2}$, $\dots$, $c_{k}$ or
$c_{1}$, $c_{k+1}$, $c_{2}$, $\dots$, $c_{k}$ is an acyclic ordering of $T$.
In both case we have $c_{k+1}\to c_{k}$. Recall that
$C_{k}=x_{k}y_{k}z_{k}x_{k}$. Since $c_{k}$ is the sink in $T$,
$d^{+}_{V(\mathcal{C})}(x_{k})=d^{+}_{V(\mathcal{C})}(y_{k})=d^{+}_{V(\mathcal{C})}(z_{k})=1$.
In particular, $d^{+}_{V(U)}(x_{k})\geq 2(k-1)>0$. Therefore, as $c_{k+1}\to
c_{k}$, $x_{k}$ has an out-neighbour $u$ in $V(U)\setminus V(C_{k+1})$. On the
other hand, since $D$ is a multipartite tournament, $D^{\prime}$ is the unique
terminal component of $U$. Thus, there is an $(u,x_{k+1})$-path $P$ in $U$.
Since $c_{k+1}\to c_{k}$, $x_{k+1}$ dominates $y_{k}$ or $z_{k}$ (or both).
Thus, $x_{k}uPx_{k+1}y_{k}C_{k}[y_{k},x_{k}]x_{k}$ or
$x_{k}uPx_{k+1}z_{k}C_{k}[z_{k},x_{k}]x_{k}$ is a cycle with length at least 4
in $D[V(U)\cup V(C_{k})\cup V(C_{k+1})]$ which together with $k-1$ 3-cycles in
$\mathcal{C}\setminus\\{C_{k}\\}$ implies that $\kappa^{k}(D)\geq 2$, a
contradiction. This completes the proof. ∎

######  Theorem 2.9.

Let $D$ be a $3$-partite tournament with $\delta^{+}(D)\geq 2k-1$. Then
$\kappa^{k}(D)=1$ if and only if the only nontrivial strong component is the
terminal strong component $D^{\prime}$ and $D^{\prime}$ is a member of
$BT(n_{1},\ldots,n_{2k})$, where $n_{i}\geq 2k-1$ for $i=1,\ldots,2k$.

###### Proof.

Let $D$ be a 3-partite tournament with $\delta^{+}(D)\geq 2k-1$ and
$\kappa^{k}(D)=1$. By Theorem 2.8, $D$ is triangle-free. Let
$D_{1},\ldots,D_{s}$ be an acyclic ordering of the strong components of $D$.
Note that $\delta^{+}(D_{s})\geq\delta^{+}(D)\geq 2k-1$. Then by Theorem 1.5
and Theorem 2.1, $D_{s}$ is a member of $BT(n_{1},\ldots,n_{2k})$, where
$n_{i}\geq 2k-1$ for every $i\in[2k]$. One can observe that there exists a
collection of $k-1$ vertex-disjoint cycles $\mathcal{C}$ in $D_{s}$, where
$k-2$ of them are 4-cycles and one of them is 6-cycle. For $i\not=s$, if
$D_{i}$ is strong, then there exists a cycle $C^{\prime}$ in $D_{i}$, and
$\kappa^{k}(\mathcal{C}\cup\\{C^{\prime}\\})\geq 2$, a contradiction. Thus,
for every $i\neq s$, $D_{i}$ is a single vertex. This completes the proof. ∎

##  3 Discussion

We proved the Bermond-Thomassen conjecture for triangle-free multipartite
tournaments and 3-partite tournaments. We also characterized 3-partite
tournaments with minimum out-degree at least $2k-1$ ($k\geq 2$) in which every
$k$ vertex-disjoint directed cycles have the same length. This conjecture
seems to be very difficult even for multipartite tournaments. In fact, we were
unable to verify the conjecture even for 4-partite tournaments, in general. It
would be interesting to see a proof of this conjecture for multipartite
tournaments.

######  Problem 3.1.

Verify that every multipartite tournament $D$ with $\delta^{+}(D)\geq 2k-1$
has $k$ vertex-disjoint cycles.

Perhaps, it’d be useful to consider the following two subproblems as a
starting point. An undirected graph $G=(V,E)$ is a split graph if $V$ can be
partitioned into two non-empty subsets $X$ and $Y$ such that $G[X]$ is a
complete graph and $G[Y]$ has no edges. A split graph $G=(V,E)$ with partition
$X,Y$ of $V$ is a complete split graph if every pair of vertices $x\in X$ and
$y\in Y$ is adjacent. A complete split oriented graph is an orientation of a
complete split graph. Note that a complete split digraph is a multipartite
tournament.

######  Problem 3.2.

Verify that every complete split oriented graph $D$ with $\delta^{+}(D)\geq
2k-1$ has $k$ vertex-disjoint cycles.

######  Problem 3.3.

Verify that every $4$-partite tournament $D$ with $\delta^{+}(D)\geq 2k-1$ has
$k$ vertex-disjoint cycles.

### Acknowledgement

This research of Li was supported by China Scholarship Council (CSC) and the
Natural Science Basic Research Planin Shaanxi Province of China (Grant No.
2020JM-133). This research of Wang and Zhou was supported by China Scholarship
Council (CSC). This research of Yeo was supported by the danish research
council under grant number DFF-7014-00037B.

## References

  * [1] Noga Alon.  Disjoint directed cycles.  J. Comb. Theory, Ser. B, 68(2):167–178, 1996.
  * [2] Yandong Bai, Binlong Li, and Hao Li.  Vertex-disjoint cycles in bipartite tournaments.  Discret. Math., 338(8):1307–1309, 2015.
  * [3] Jørgen Bang-Jensen, Stéphane Bessy, and Stéphan Thomassé.  Disjoint 3-cycles in tournaments: A proof of the Bermond-Thomassen conjecture for tournaments.  J. Graph Theory, 75(3):284–302, 2014.
  * [4] Jørgen Bang-Jensen and Gregory Z. Gutin.  Basic terminology, notation and results.  In Jørgen Bang-Jensen and Gregory Z. Gutin, editors, Classes of Directed Graphs, Springer Monographs in Mathematics, pages 1–34. Springer, 2018.
  * [5] Jean-Claude Bermond and Carsten Thomassen.  Cycles in digraphs- a survey.  J. Graph Theory, 5(1):1–43, 1981.
  * [6] Matija Bucic.  An improved bound for disjoint directed cycles.  Discret. Math., 341(8):2231–2236, 2018.
  * [7] Bin Chen and An Chang.  Disjoint cycles in tournaments and bipartite tournaments.  J. Graph Theory, in press, 2023.
  * [8] Yubao Guo and Lutz Volkmann.  Cycles in multipartite tournaments.  J. Comb. Theory, Ser. B, 62(2):363–366, 1994.
  * [9] Nicolas Lichiardopol, Attila Pór, and Jean-Sébastien Sereni.  A Step toward the Bermond–Thomassen Conjecture about Disjoint Cycles in Digraphs.  SIAM J. Discret. Math., 23(2):979–992, 2009.
  * [10] Carsten Thomassen.  Disjoint cycles in digraphs.  Combinatorica, 3(3):393–396, 1983.


# Radiative Carrier Lifetime in Ge1-xSnx Mid-Infrared Emitters

Gérard Daligou  Department of Engineering Physics, École Polytechnique de
Montréal, C.P. 6079, Succ. Centre-Ville, Montréal, Québec, Canada H3C 3A7
Anis Attiaoui  Department of Engineering Physics, École Polytechnique de
Montréal, C.P. 6079, Succ. Centre-Ville, Montréal, Québec, Canada H3C 3A7
Simone Assali  Department of Engineering Physics, École Polytechnique de
Montréal, C.P. 6079, Succ. Centre-Ville, Montréal, Québec, Canada H3C 3A7
Patrick Del Vecchio  Department of Engineering Physics, École Polytechnique de
Montréal, C.P. 6079, Succ. Centre-Ville, Montréal, Québec, Canada H3C 3A7
Oussama Moutanabbir<EMAIL_ADDRESS>Department of Engineering
Physics, École Polytechnique de Montréal, C.P. 6079, Succ. Centre-Ville,
Montréal, Québec, Canada H3C 3A7

###### Abstract

Ge1-xSnx semiconductors hold the premise for large-scale, monolithic mid-
infrared photonics and optoelectronics. However, despite the successful
demonstration of several Ge1-xSnx-based photodetectors and emitters, key
fundamental properties of this material system are yet to be fully explored
and understood. In particular, little is known about the role of the material
properties in controlling the recombination mechanisms and their consequences
on the carrier lifetime. Evaluating the latter is in fact fraught with large
uncertainties that are exacerbated by the difficulty to investigate narrow
bandgap semiconductors. To alleviate these limitations, herein we demonstrate
that the radiative carrier lifetime can be obtained from straightforward
excitation power- and temperature-dependent photoluminescence measurements. To
this end, a theoretical framework is introduced to simulate the measured
spectra by combining the band structure calculations from the k.p theory and
the envelope function approximation (EFA) to estimate the absorption and
spontaneous emission. The model computes explicitly the momentum matrix
element to estimate the strength of the optical transitions in single bulk
materials, unlike the joint density of states (JDOS) model which assumes a
constant matrix element. Based on this model, the temperature-dependent
emission from Ge0.83Sn0.17 samples at a biaxial compressive strain of $-1.3\%$
was investigated. The simulated spectra reproduce accurately the measured data
thereby enabling the evaluation of the steady-state radiative carrier
lifetimes, which are found in the $3$-$22\text{\,}\mathrm{ns}\text{/}$ range
for temperatures between $10$ and $300\text{\,}\mathrm{K}\text{/}$ at an
excitation power of $0.9\text{\,}\mathrm{kW}\text{/}{\mathrm{cm}}^{2}$. For a
lower power of $0.07\text{\,}\mathrm{kW}\text{/}{\mathrm{cm}}^{2}$, the
obtained lifetime has a value of $1.9\text{\,}\mathrm{ns}\text{/}$ at
$4\text{\,}\mathrm{K}\text{/}$. The demonstrated approach yielding the
radiative lifetime from simple emission spectra will provide valuable inputs
to improve the design and modeling of Ge1-xSnx-based devices.

##  I INTRODUCTION

Ge1-xSnx alloys constitute an emerging class of group IV semiconductors
providing a tunable narrow bandgap, which has been highly attractive to
implement scalable, silicon-compatible mid-infrared photonic and
optoelectronic devices [1]. This potential becomes increasingly significant
with the recent progress in nonequilibrium growth processes enabling high Sn
content Ge1-xSnx layers and heterostructures leading to the demonstration of a
variety of monolithic mid-infrared emitters and detectors [4, 5, 6, 7, 8, 9,
10, 11, 2, 12, 13, 14, 15, 16, 17, 18, 3]. Notwithstanding the recent
developments in device engineering, the impact of structural characteristics
on the basic behavior of charge carriers is yet to be fully understood. This
includes the role of Sn content, lattice strain, and growth defects in shaping
the nature and magnitude of the recombination mechanisms and their
consequences on the carrier lifetime. Particularly, investigating the latter
remains a daunting task due to the lack of methods and tools that can be
applied to probe charge carriers in narrow bandgap materials. For instance,
time-resolved photoluminescence (PL) can hardly be applied to investigate
materials at emission wavelengths in the mid-infrared range as high-speed
detectors covering this range are not broadly available. Thus, the very few
reported time-resolved studies concern Ge1-xSnx emitting below
$2.3\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$ corresponding to a
relatively low Sn content and/or highly compressively strained materials [19,
20, 21].

In an attempt to circumvent the aforementioned limitations, a recent study
employed time-resolved PL with a nonlinear crystal allowing the up-conversion
of photons emitted to a shorter wavelength that can be detected by a
conventional silicon-based avalanche photodiode [19]. An effective carrier
lifetime of $217\text{\,}\mathrm{ps}\text{/}$ at
$20\text{\,}\mathrm{K}\text{/}$ was estimated for Ge0.875Sn0.125 with
$-0.55\%$ strain using this method [19]. Additionally, by investigating spin-
dependent optical transitions leveraging the Hanle effect under steady-state
excitation, systematic studies combining modeling and magneto-PL analysis of
pseudomorphic layers at a Sn content below $10\%$ reported a radiative
lifetime in the $0.5$-$2.5\text{\,}\mathrm{ns}\text{/}$ range at
$10\text{\,}\mathrm{K}\text{/}$ [20]. However, significantly higher carrier
lifetimes reaching $450\text{\,}\mathrm{ns}\text{/}$ were recently reported
for Ge1-xSnx ($x<0.06$) grown on InAlAs buffer layers as measured by
contactless microwave photoconductive decay [21]. This scarcity of studies on
carrier dynamics in narrow bandgap Ge1-xSnx semiconductors limits the
understanding of their fundamental behavior and burdens the development of
accurate and predictive models for Ge1-xSnx-based mid-infrared optoelectronic
devices.

In this work, we demonstrate that straightforward PL analyses along with the
proper theoretical framework are sufficient to alleviate these challenges and
extract the radiative carrier lifetime in Ge1-xSnx mid-infrared emitters and
evaluate its evolution as a function of temperature. The approach relies on
the simulation of the experimental PL spectra by combining the band structure
calculations using the k.p formalism together with the EFA to estimate the
absorption and spontaneous emission spectra. Unlike the JDOS model, in which
the momentum matrix element is considered constant, the oscillator strengths
are explicitly computed in this model.

In the following sections, the model is described followed by the experimental
demonstration using as-grown Ge0.83Sn0.17 layers, emitting at wavelengths
above $3\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$.

##  II Theoretical framework

The PL spectrum intensity is usually determined using the direct interband
emission theory and the spontaneous emission spectrum [22]. Indeed, by
considering a slab of homogeneously excited material, Lasher and Stern [23]
and Würfel [24] expressed the external flux of spontaneous radiative emission
in terms of the spectral absorptivity under non-equilibrium conditions in
terms of the quasi-Fermi level splitting, $\Delta\mu=\mu_{e}-\mu_{h}$. The
resulting “Lasher-Stern-Würfel” (LSW) equation is

$I_{\text{PL}}(E)=\frac{2\pi}{\hbar^{3}c^{2}}\frac{E^{2}a(E)}{\exp(\frac{E-\Delta\mu}{k_{\text{B}}T})-1}=\frac{a(E)}{\alpha(E)}\cdot\frac{r^{\text{sp}}(E)}{4n_{r}^{2}}$
(1)

where $r^{\text{sp}}$ is the internal spontaneous emission spectrum,
$\alpha(E)$ is the absorption spectrum and $n_{r}$ the refractive index of the
medium. $a(E)$ is the spectral absorptivity defined as expression (2) with $R$
the reflection from the outside onto the sample surface, and $d$ the thickness
of the conceptual slab [25].

$a(E)=(1-R(E))\Big{[}1-\exp\quantity(-\alpha(E)d)\Big{]}$ (2)

Note that $d$ is also considered as a characteristic length scale over which
carriers are generated, travel and recombine radiatively [25]. Based on this
definition, this parameter should be inversely proportional to the absorption
coefficient at the excitation wavelength i.e. $d\approx
1/\alpha(\lambda_{\text{laser}})$. However, the PL spectrum will mostly be
centered around the bandgap energy where the $\alpha(E)$ is at least one or
two orders of magnitude smaller than $\alpha(\lambda_{\text{laser}})$, in the
case of non-resonant excitation. Therefore $\alpha d<<1$ and the absorbance
$a(E)$ can be simplified by expanding the exponent with a Taylor series such
that $a(E)\approx A\alpha(E)$. In that case, the PL intensity from equation
(1) becomes

$I_{\text{PL}}(E)\approx\frac{2\pi
A}{\hbar^{3}c^{2}}\frac{E^{2}\alpha(E)}{\exp(\frac{E-\Delta\mu}{k_{\text{B}}T})-1}\approx\frac{A}{4n_{r}^{2}}\cdot
r^{\text{sp}}(E)$ (3)

With this approximation, The PL spectrum intensity is therefore entirely
defined by the internal spontaneous emission spectrum or the absorption
spectrum depending on the formula used.

The spontaneous emission spectrum $r^{\text{sp}}$ is calculated using the
Fermi’s golden rule [26, 27] and the perturbation theory as described in
equation (4):

$\displaystyle r^{\text{sp}}(\hbar\omega)$
$\displaystyle=\frac{\hbar\omega\Lambda}{V}\sum_{c,v,\overrightarrow{{\bf\it
k}}}\absolutevalue{\matrixelement{\Phi_{c}(\overrightarrow{{\bf\it
k}})}{\frac{\hbar}{m_{0}}\hat{{\bf\it\varepsilon}}.\overrightarrow{{\bf\it
p}}}{\Phi_{v}(\overrightarrow{{\bf\it
k}})}}^{2}f\left(\epsilon_{c},\mu_{e}\right)$
$\displaystyle\quad\times\delta\left(\epsilon_{c}(\overrightarrow{{\bf\it
k}})-\epsilon_{v}(\overrightarrow{{\bf\it
k}})-\hbar\omega\right)\left[1-f\quantity(\epsilon_{v},\mu_{h})\right],$ (4)

where, $\Lambda=n_{r}e^{2}/\pi c^{3}\varepsilon_{0}\hbar^{4}$ is a material-
related constant with $e$ the elementary charge, $n_{r}$ the refractive index
of the material, and $c$ the speed of light in vacuum. $V$ is the volume of
the states in the $\overrightarrow{{\bf\it k}}$-space. The summations are done
over the different values of $\overrightarrow{{\bf\it k}}$ in the Brillouin
zone (BZ) to account for the possible transitions between the conduction and
the valence bands. The Dirac delta distribution is used to limit the
transitions to those with an energy difference of $\hbar\omega$, the photon
energy. Moreover, $M^{2}_{i,f}(\overrightarrow{{\bf\it
k}})=\absolutevalue{\matrixelement{\Phi_{f}(\overrightarrow{{\bf\it
k}})}{\frac{\hbar}{m_{0}}\hat{{\bf\it\varepsilon}}.\overrightarrow{{\bf\it
p}}}{\Phi_{i}(\overrightarrow{{\bf\it k}})}}^{2}$ represents the strength of
the transition from the state $\ket{\Phi_{i}}$ to $\ket{\Phi_{f}}$ with
$\hat{{\bf\it\varepsilon}}$ the polarization unit vector and
$\overrightarrow{{\bf\it p}}$ the momentum matrix operator. Finally, the
Fermi-Dirac statistic is used to account for the probability occupation of the
different states, with $f$ given by
$f(\epsilon,\mu)=\left[1+\exp\left(\frac{\epsilon-\mu}{k_{\text{B}}T}\right)\right]^{-1}$
in which $\epsilon$ is the energy and $\mu$ the Fermi-level of the charge
carrier described by the function.

The computation of the spontaneous emission spectrum requires prior knowledge
of the band structure of the semiconductors, the momentum matrix elements, and
the quasi-Fermi levels, as seen in equation (4). In the current literature,
for a single bulk direct bandgap semiconductor, $r^{\text{sp}}$ is commonly
computed using the JDOS model [28, 29, 30]. This model relies on the parabolic
band approximation (PBA) which leads to a set of relatively easy analytical
formulas. It is mostly accurate for a non-degenerately doped semiconductor in
weak-injection regime with the quasi-Fermi levels lying within the bandgap and
away from the different band edges by several $k_{\text{B}}T$, where $k_{B}$
is the Boltzmann’s constant and $T$ the temperature ($\Delta\mu\approx 0$).
This model was extended in [31] to account for different excitation regimes by
explicitly evaluating the quasi-Fermi level splitting $\Delta\mu$, and the
non-equilibrium absorption spectrum in equation (3). However, it still relies
on the PBA which restricts the analysis. Indeed, for higher excitation power
and/or doping concentration, $\mu_{e}$ and $\mu_{h}$ would shift towards, and
even beyond, the band edges where the PBA should be less accurate. Besides,
for a biaxially strained semiconductor, the $\overrightarrow{{\bf\it k}}$
direction degeneracies in the BZ are expected to be broken. In this situation,
the band dispersion would be increasingly anisotropic, thus challenging one of
the core principles of the PBA. A more accurate theoretical framework is
therefore required for the computation of the spontaneous emission spectrum
$r^{\text{sp}}$, and the description of the measured PL spectra.

Figure 1: Different steps followed in the computation of the spontaneous
emission spectra

###  II.1 Spontaneous emission spectrum and eight-band k.p formalism

The spectrum $r^{\text{sp}}$ is computed using the eight-band k.p formalism
together with the EFA [32], following the simulation workflow summarized in
Fig. 1. The eight-band k.p Ge1-xSnx material parametrization is based on early
reports [33, 34, 35], while strain implementation is based on the Bir-Pikus
formalism [36]. To account for the inaccuracy of the Vegard’s law to estimate
the bandgaps of Ge1-xSnx alloys, bandgap bowing parameters are introduced for
L and $\Gamma$ high-symmetry points. Unlike the JDOS model (and all the
different models relying on the PBA), the evolution of the strength of the
optical transitions with the wave vector $\overrightarrow{{\bf\it k}}$ is
explicitly computed using the formalism developed by Szmulowicz [37]. If
$\ket{\Phi_{i}}$ and $\ket{\Phi_{f}}$ are the initial and final states within
the EFA, the strength of the transition is given by [37]:

$\displaystyle\matrixelement{\Phi_{i}}{\frac{\hbar}{m_{0}}\hat{\mathbf{e}}.\overrightarrow{{\bf\it
p}}}{\Phi_{f}}$
$\displaystyle=\sum_{\mu,\,\nu}\Phi_{i,\mu}^{*}(\overrightarrow{{\bf\it
k}})\left(\hat{\mathbf{\varepsilon}}\dotproduct\partialderivative{\mathcal{H}_{\mu\nu}(\overrightarrow{{\bf\it
k}})}{\overrightarrow{{\bf\it k}}}\right)\Phi_{f,\nu}(\overrightarrow{{\bf\it
k}})$
$\displaystyle=\sum_{\mu,\,\nu,\,l}\Phi_{i,\mu}^{*}(\overrightarrow{{\bf\it
k}})\left[\varepsilon_{l}\cdot\left(\partialderivative{\mathcal{H}_{\mu\nu}}{k_{l}}\right)\right]\Phi_{f,\nu}(\overrightarrow{{\bf\it
k}}),$ (5)

where $\Phi_{i,\mu}$ and $\Phi_{f,\nu}$ are the coefficients of the envelope
function vector related to the states $\ket{\Phi_{i}}$ and $\ket{\Phi_{f}}$,
respectively. The unit vector $\hat{\mathbf{\varepsilon}}$ gives the
polarization of the incident light, while
$\partial\mathcal{H}_{\mu\nu}(\overrightarrow{{\bf\it
k}})/\partial\overrightarrow{{\bf\it k}}$ is the derivative of the k.p
Hamiltonian with respect to the wave vector $\overrightarrow{{\bf\it k}}$. The
expressions of the different momentum matrices
$\partial\mathcal{H}_{\mu\nu}/\partial k_{l}$ can be found in the Supporting
information S1 111See Supplemental Material for more details about the
momentum matrix elements and the special points approximation .

For a given value of the optically injected carrier density $\Delta n$, if
$n_{0}$ and $p_{0}$ denote the total electrons and holes’ charge densities at
thermal equilibrium, the quasi-Fermi levels $\mu_{e}$ (for electrons) and
$\mu_{h}$ (for holes) are determined by solving the set of equations (6):

$\displaystyle\begin{aligned} n_{0}+\Delta
n&=\frac{1}{(2\pi)^{3}}\sum_{i\in\mathrm{CB}}\int_{\text{BZ}}\frac{d^{3}\overrightarrow{{\bf\it
k}}}{1+\exp\left(\frac{\epsilon_{i}(\overrightarrow{{\bf\it
k}})-\mu_{e}}{kT}\right)}\\\ p_{0}+\Delta
n&=\frac{1}{(2\pi)^{3}}\sum_{i\in\mathrm{VB}}\int_{\text{BZ}}\frac{d^{3}\overrightarrow{{\bf\it
k}}}{1+\exp\left(\frac{\mu_{h}-\epsilon_{j}(\overrightarrow{{\bf\it
k}})}{kT}\right)}\end{aligned}$ (6)

Herein, the conduction band electrons are assumed to be shared between the
$\Gamma$ and $L$ valleys. This assumption is only relevant when the energy
band offset between these valleys is relatively close to the thermal energy
$k_{\text{B}}T$ to enable the electrons to transition between them. The
carrier concentration $n_{0}$ and $p_{0}$ are evaluated after solving the
electroneutrality equation to estimate the thermal equilibrium Fermi level
$E_{f}$. Besides, the computation of the integrals over the BZ, required for
estimating the quasi-Fermi levels and $r^{\text{sp}}$, relies on the special-
lines approximation (SLA) [39]. Within this approximation, the three-
dimensional BZ integrals are replaced by a sum of one-dimensional integrals
over some characteristic directions (denoted as “special”) of the crystal
lattice. These directions could, for example, be the symmetry directions used
in the eight-band k.p formalism. If we denote by $\mathcal{L}$ the set of the
special directions, the electrons density from equation (6) becomes

$\displaystyle n$
$\displaystyle=\frac{1}{2\pi^{2}}\sum_{\begin{subarray}{c}D\in\mathcal{L}\\\
i\in\mathrm{CB}\end{subarray}}w_{D}\left(\int_{0}^{k_{\text{BZ}}}\frac{k^{2}_{D}\differential
k_{D}}{1+\exp\left(\frac{\epsilon_{i}(k_{D})-\mu_{e}}{kT}\right)}\right)$ (7)

with $w_{D}$ the weight of the direction $D$, $\epsilon_{i}(k_{\text{D}})$ the
energy of the conduction band $i$ at $k_{\text{D}}$ and $k_{\mathrm{BZ}}\sim
0.5$ (units of $\pi/a_{0}$, $a_{0}$ being the lattice constant of the
material) for the eight-band k.p model to still be accurate. Depending on the
computation, the exact value of the upper limit $k_{\text{BZ}}$ could be
neglected since the integrands are expected to vanish rapidly while increasing
the value of $k_{\text{D}}$. More information about all the different
directions considered in our framework can be found in the Supporting
information S2 [38]. Unlike the PBA, which leads to parabolic and isotropic-
like band structure, this method accounts for the anisotropy and the non-
parabolicity of the bands obtained with the k.p theory and approximates the
warping of real bands. It is, therefore, expected to be more accurate.

The theoretical spontaneous emission spectrum developed previously is often
insufficient to accurately describe the PL spectrum. Indeed, the sub-bandgap
emission resulting from carrier disorders and broadening mechanisms in the
materials is not accounted for in equation (4) [40, 41, 42]. To include these
contributions, the theoretical spectrum $r^{\text{sp}}_{\text{ideal}}$ from
equation (4) is convoluted with a broadening function $\mathcal{B}$, as
outlined in equation (8).

$r^{\text{sp}}(\hbar\omega)=\int_{-\infty}^{+\infty}r^{\text{sp}}_{\text{ideal}}(\epsilon)\cdot\mathcal{B}(\hbar\omega-\epsilon)\mathrm{d}\epsilon$
(8)

The broadening function $\mathcal{B}$ is usually chosen as a Gaussian or a
Lorentzian to account for the inhomogeneous and homogeneous broadening
mechanisms, respectively. However, the Lorentzian function was reported to
sometimes overestimate the effects of the homogeneous broadening due to its
slowly decaying tails. For that reason, it is usually replaced by a hyperbolic
secant function [43].

Figure 2: (a) Power-dependent photoluminescence at
$4\text{\,}\mathrm{K}\text{/}$ for the -1.27% biaxially strained Ge0.83Sn0.17.
The scatter points are from the measurements while the black lines are the
results from the simulations. (b) Evolution of the extracted quasi-Fermi
levels ($\mu_{e},\,\mu_{h}$) with the power density.

###  II.2 Spontaneous emission intensity and steady-state radiative carrier
lifetime

We shall use the symbol $R_{\text{sp}}^{\hat{\mathbf{\varepsilon}}}$ to denote
the rate of polarization-dependent spontaneous emission per unit volume, where
$\hat{\mathbf{\varepsilon}}$ gives the polarization of the incident light. The
quantity $R_{\text{sp}}^{\hat{\mathbf{e}}}$ is generally defined as the
integral of the polarization-dependent spontaneous emission spectrum over the
range of photon energy bigger than the band-gap of the material, and given by
equation (9) [26].

$R_{\text{sp}}^{\hat{\mathbf{\varepsilon}}}=\int_{0}^{+\infty}r^{\text{sp}}_{\hat{\mathbf{\varepsilon}}}(\hbar\omega)\differential{\hbar\omega}$
(9)

For an unpolarized incident light, the total spontaneous emission rate per
unit volume $R_{sp}$ is defined as the average of the contributions from the
three polarizations defined by the unit vectors
$\hat{\mathbf{\varepsilon}}_{x}=(1,0,0)$, $\hat{\mathbf{\varepsilon}}_{y}$ and
$\hat{\mathbf{\varepsilon}}_{z}$.

The steady-state radiative carrier lifetime $\tau_{\text{rad}}$ is determined
by the net rate of spontaneous emission $R_{\text{sp}}^{\text{net}}$ and the
density of photo-excited carriers $\Delta n$ (equation (10)).
$R_{\text{sp}}^{\text{net}}$ is defined as the amount by which the non-
equilibrium spontaneous recombination rate $R_{\text{sp}}$ exceeds the thermal
equilibrium generation rate $G_{0}$, which is the same as the thermal
equilibrium spontaneous emission rate.

$\displaystyle\tau_{\text{rad}}$ $\displaystyle=\frac{\Delta
n}{R_{\text{sp}}^{\text{net}}}$ (10)

$R_{\text{sp}}^{\text{net}}$ is usually estimated using equation (11), in
which $B$ is a material-dependent parameter known as the bimolecular
recombination coefficient. In that case, $\tau_{\text{rad}}$ becomes
relatively easy to compute, as presented in equation (12).

$R_{\text{sp}}^{\text{net}}=R_{\text{sp}}^{\text{neq}}-R_{\text{sp}}^{\text{eq}}=B(np-
n_{0}p_{0})$ (11) $\displaystyle\tau_{\text{rad}}$
$\displaystyle=\frac{1}{B\quantity(\Delta n+n_{0}+p_{0})}$ (12)

The bimolecular recombination coefficient $B$ is typically assumed to be
independent of $\Delta n$ (and, therefore, the quasi-Fermi levels). However,
this approximation is not always accurate. For example, $B$ was previously
shown to vary linearly with the excess carrier density $\Delta n$ in III-V
semiconductors [44, 45, 46]. For that reason, it is reasonable to rely only on
equations (9) and (10), which state the general case without any specific
approximations.

##  III Results and discussion

The accuracy of the established theoretical framework has been evaluated
through the analysis of the PL spectra recorded as a function of the
excitation power and temperature from Ge0.83Sn0.17 layers [47]. The epitaxial
growth of these layers was achieved using low-pressure chemical vapor
deposition (LP-CVD) starting from a $600$ \- $700\text{\,}\mathrm{nm}\text{/}$
Ge virtual substrate on a $4\text{\,}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{h}$
Si wafer. To ensure the growth of a Ge0.83Sn0.17 layer with a uniform Sn
composition, a multilayer heterostructure consisting of top layer (TL)/middle
layer (ML)/bottom layer (BL) was grown while the incorporation of Sn in each
layer is controlled by adjusting the growth temperature. More details on the
growth and characterization of Ge0.83Sn0.17 material can be found in [47].

In as-grown Ge0.83Sn0.17 layers, the band alignment favors the electrons and
holes diffusion to the TL, where they should recombine. Indeed, the PL spectra
are confirmed to originate from carrier recombination in this specific layer
[47]. Therefore, from a theoretical standpoint, it would be judicious to
analyze the PL results as if they were emitted by a bulk GeSn material with a
17 at. % Sn composition. On this basis, the different power-dependent PL
spectra, recorded at 4 K, were simulated by iteratively evaluating
$r^{\text{sp}}$ as well as the excess carrier concentration $\Delta n$ and
$\gamma$ the full width at half maximum (FWHM) of the broadening function. To
solve equations (6), and extract the quasi-Fermi levels couple
($\mu_{e},\,\mu_{h}$), the p-type background doping was considered to be
around ${10}^{15}\text{\,}{\mathrm{cm}}^{-3}$ at
$4\text{\,}\mathrm{K}\text{/}$. This value was chosen with reference to the
p-type background doping estimated between
$1\text{\times}{10}^{17}\text{\,}\text{/}$ and
$5\text{\times}{10}^{17}\text{\,}{\mathrm{cm}}^{-3}$ at 300 K [48].

Figure 3: (a) Evolution of the radiative carrier lifetime as a function of the
excitation power density $P_{\text{exc}}$, (b) Evolution of the $B$
coefficient as a function of the excess carrier density $\Delta n$. The solid
spheres represent the values of $B$ computed from the extracted $\Delta n$,
$R_{\text{sp}}$ using equation (11), while the black line is the result of a
fit using a rational function.

Fig. 2(a) displays the measured and simulated spectra for the as-grown
Ge0.83Sn0.17 material. For each power density, a coefficient of determination
($R^{2}$) of around $99.5\%$ is obtained, thus highlighting the accuracy of
the simulated spectra. Moreover, the evolution of the extracted quasi-Fermi
levels with the excitation power density $P_{\text{exc}}$ is outlined in Fig.
2(b).

Figure 4: Evolution of the radiative carrier lifetime as function of
temperature for the as-grown Ge0.83Sn0.17.

For a p-type background doping of ${10}^{15}\text{\,}{\mathrm{cm}}^{-3}$, the
thermal equilibrium Fermi level $E_{\text{F}}$ is about
$42.73\text{\,}\mathrm{meV}\text{/}$. As shown in Fig. 2(b), the non-
degenerate semiconductor approximation is not appropriate here since
$E_{\text{F}}$ is less than the top valence band edge located around
$42.93\text{\,}\mathrm{meV}\text{/}$. Starting from a power density of
$67.95\text{\,}\mathrm{W}\text{/}{\mathrm{cm}}^{2}$, both the electrons and
holes quasi-Fermi levels start to deviate from $E_{\text{F}}$. In fact, a
progressive increase from $450.02$ to $462.57\text{\,}\mathrm{meV}\text{/}$ is
observed for the quasi-Fermi level $\mu_{e}$, causing the electron
concentration to increase. Simultaneously, the holes quasi-Fermi level
$\mu_{h}$ decreases while remaining very close to the thermal equilibrium
level with a maximum offset of $6.88\text{\,}\mathrm{meV}\text{/}$ at
$5.4\text{\,}\mathrm{kW}\text{/}{\mathrm{cm}}^{2}$. While these variations may
be perceived as small, they are not insignificant. Indeed, with the thermal
energy of about $0.34\text{\,}\mathrm{meV}\text{/}$ at 4 K, one should expect
a noticeable increase in the spontaneous emission intensity $R_{\text{sp}}$.
Moreover, using equation (10), the steady-state radiative carrier lifetime
$\tau_{\text{rad}}$ was extracted (Fig. 3(a)) and shown to decrease from
$3.52\text{\,}\text{/}$ to $1.89\text{\,}\mathrm{ns}\text{/}$ in the range of
power density used in this study.

Besides, the accuracy of equation (11) was also evaluated using the different
parameters obtained from the analysis above. Indeed, the bimolecular
recombination coefficient $B$ was computed from the extracted values of
$R_{sp}$ and $\Delta n$ and its behavior is outlined in Fig. 3(b). Rather than
being constant, it decreases with $\Delta n$, as suggested earlier for III-V
semiconductors [44, 45]. However, its evolution for the as-grown Ge0.83Sn0.17
is not as linear as presented by Olshansky et al. for InGaAsP and AlGaAs light
sources [46]. In fact, after performing a fit of the obtained data, $B$ was
shown to evolve with $\Delta n$ following a rational function (Fig. 3(b)).
Additionally, for $\Delta n$ above ${10}^{15}\text{\,}{\mathrm{cm}}^{-3}$, the
values extracted were shown to be lower than the value of
$1.04\text{\times}{10}^{-7}\text{\,}{\mathrm{cm}}^{3}\text{/}\mathrm{s}$
computed assuming parabolic band dispersion and the non-degenerate
semiconductor approximation.

The impact of temperature on the steady-state radiative carrier lifetime has
also been investigated. Herein, assuming a $p$-type background doping of
$1\text{\times}{10}^{15}\text{\,}{\mathrm{cm}}^{-3}$ at
$4\text{\,}\mathrm{K}\text{/}$ and
$1\text{\times}{10}^{17}\text{\,}{\mathrm{cm}}^{-3}$ at
$300\text{\,}\mathrm{K}\text{/}$, which is in line with recent measurements
[48], the evolution of the doping with temperature was estimated. Using these
values, the temperature-dependent PL spectra were simulated with the
theoretical estimation of the spontaneous emission spectrum from the framework
described above, and the evolution of $\tau_{\text{rad}}$ was extracted for
the as-grown Ge0.83Sn0.17, as displayed in Fig. 4. Note that from this
analysis, a minimum $R^{2}$ factor of about $98\%$ was observed throughout the
$4$-$300\text{\,}\mathrm{K}\text{/}$ range. The estimated steady-state
radiative carrier lifetime $\tau_{\text{rad}}$ increases with the temperature
from $\sim$3.2\text{\,}\mathrm{ns}\text{/}$$ at
$10\text{\,}\mathrm{K}\text{/}$ to $\sim$22.2\text{\,}\mathrm{ns}\text{/}$$ at
$300\text{\,}\mathrm{K}\text{/}$. These values are very comparable to the
reported recombination lifetimes in literature for III-V compound
semiconductors, which are generally in the nanoseconds range [49, 50, 51].
They are also of the same order of magnitude as the values for Ge calculated
from first principles [52]. Indeed, the radiative lifetime for Ge in the
diamond structure was shown to be around $10\text{\,}\mathrm{ns}\text{/}$ for
$T$ below $300\text{\,}\mathrm{K}\text{/}$.

Finally, to appreciate the radiative emission strength of Ge1-xSnx with
respect to other direct bandgap semiconductors, we compare the radiative
emission rate or more precisely the bimolecular recombination coefficient $B$.
Using the same process as for the power-dependent PL, $B$ is extracted as a
function of the temperature. From this analysis, $B$ is found to evolve
following the allometric power law $aT^{b}$ with $b\approx-1.5143$, and
reaching
$3.81\text{\times}{10}^{-10}\text{\,}{\mathrm{cm}}^{3}\text{/}\mathrm{s}$ at
$240\text{\,}\mathrm{K}\text{/}$. This value is comparable to those extracted
at $300\text{\,}\mathrm{K}\text{/}$ for GaAs
($3.5\text{\times}{10}^{-10}\text{\,}{\mathrm{cm}}^{3}\text{/}\mathrm{s}$),
InP
($1.2\text{\times}{10}^{-10}\text{\,}{\mathrm{cm}}^{3}\text{/}\mathrm{s}$),
and hexagonal Si0.20Ge0.80
($$0.7\text{\times}{10}^{-10}\text{\,}\text{/}$-$11\text{\times}{10}^{-10}\text{\,}{\mathrm{cm}}^{3}\text{/}\mathrm{s}$$).

##  IV CONCLUSION

To circumvent the limitations in the experimental studies of carrier dynamics
in narrow bandgap Ge1-xSnx materials, this work demonstrates a straightforward
method to obtain the carrier radiative lifetime from simple PL spectra. The
approach relies on a theoretical framework combining the band structure
calculations using the k.p formalism together with the envelope function
approximation to estimate the absorption and spontaneous emission spectra.
This framework simulates accurately the experimental measurements thereby
allowing the evaluation of the steady-state radiative carrier lifetime from
the net rate of spontaneous emission and the density of photo-excited
carriers. For a Ge0.83Sn0.17 material under an in-plane biaxial compressive
strain $\varepsilon_{\parallel}=-1.3\%$, the analysis revealed a lifetime
$\tau_{\text{rad}}$ in the nanoseconds range increasing from 3 to
$22\text{\,}\mathrm{ns}\text{/}$ for temperatures between 10 and
$300\text{\,}\mathrm{K}\text{/}$. Additionally, the introduced model also
solves the restrictions that are inherent to the joint density of states
(JDOS) model resulting from the parabolic band approximation (PBA) and the
weak-injection approximation.

ACKNOWLEDGEMENTS
O.M. acknowledges support from NSERC Canada (Discovery, SPG, and CRD Grants),
Canada Research Chairs, Canada Foundation for Innovation, Mitacs, PRIMA
Québec, Defence Canada (Innovation for Defence Excellence and Security,
IDEaS), the European Union’s Horizon Europe research and innovation programme
under grant agreement No 101070700 (MIRAQLS), and the US Army Research Office
Grant No. W911NF-22-1-0277.

AUTHORS INFORMATION
Corresponding Author:
†<EMAIL_ADDRESS>
Notes:
The authors declare no competing financial interest.

## References

  * Moutanabbir  _et al._ [2021] O. Moutanabbir, S. Assali, X. Gong, E. O’Reilly, C. A. Broderick, B. Marzban, J. Witzens, W. Du, S.-Q. Yu, A. Chelnokov, D. Buca, and D. Nam, Applied Physics Letters 118, 110502 (2021).
  * Atalla  _et al._ [2021a] M. R. M. Atalla, S. Assali, A. Attiaoui, C. Lemieux-Leduc, A. Kumar, S. Abdi, and O. Moutanabbir, Advanced Functional Materials 31, 2006329 (2021a).
  * Atalla  _et al._ [2022] M. R. M. Atalla, S. Assali, S. Koelling, A. Attiaoui, and O. Moutanabbir, ACS Photonics 9, 1425 (2022).
  * Buca  _et al._ [2022] D. Buca, A. Bjelajac, D. Spirito, O. Concepción, M. Gromovyi, E. Sakat, X. Lafosse, L. Ferlazzo, N. von den Driesch, Z. Ikonic, D. Grützmacher, G. Capellini, and M. El Kurdi, Advanced Optical Materials 10, 2201024 (2022).
  * Chang  _et al._ [2022] C.-Y. Chang, P.-L. Yeh, Y.-T. Jheng, L.-Y. Hsu, K.-C. Lee, H. Li, H. H. Cheng, and G.-E. Chang, Photonics Research 10, 2278 (2022).
  * Chrétien  _et al._ [2019] J. Chrétien, N. Pauc, F. Armand Pilon, M. Bertrand, Q.-M. Thai, L. Casiez, N. Bernier, H. Dansas, P. Gergaud, E. Delamadeleine, R. Khazaka, H. Sigg, J. Faist, A. Chelnokov, V. Reboud, J.-M. Hartmann, and V. Calvo, ACS Photonics 6, 2462 (2019).
  * Chrétien  _et al._ [2022] J. Chrétien, Q. M. Thai, M. Frauenrath, L. Casiez, A. Chelnokov, V. Reboud, J. M. Hartmann, M. El Kurdi, N. Pauc, and V. Calvo, Applied Physics Letters 120, 051107 (2022).
  * Elbaz  _et al._ [2020] A. Elbaz, D. Buca, N. von den Driesch, K. Pantzas, G. Patriarche, N. Zerounian, E. Herth, X. Checoury, S. Sauvage, I. Sagnes, A. Foti, R. Ossikovski, J.-M. Hartmann, F. Boeuf, Z. Ikonic, P. Boucaud, D. Grützmacher, and M. El Kurdi, Nature Photonics 14, 375 (2020).
  * Joo  _et al._ [2021] H.-J. Joo, Y. Kim, D. Burt, Y. Jung, L. Zhang, M. Chen, S. J. Parluhutan, D.-H. Kang, C. Lee, S. Assali, Z. Ikonic, O. Moutanabbir, Y.-H. Cho, C. S. Tan, and D. Nam, Applied Physics Letters 119, 201101 (2021).
  * Jung  _et al._ [2022] Y. Jung, D. Burt, L. Zhang, Y. Kim, H.-J. Joo, M. Chen, S. Assali, O. Moutanabbir, C. S. Tan, and D. Nam, Photonics Research 10, 1332 (2022).
  * Li  _et al._ [2021] X. Li, L. Peng, Z. Liu, Z. Zhou, J. Zheng, C. Xue, Y. Zuo, B. Chen, and B. Cheng, Photonics Research 9, 494 (2021).
  * Liu  _et al._ [2022] X. Liu, J. Zheng, C. Niu, T. Liu, Q. Huang, M. Li, D. Zhang, Y. Pang, Z. Liu, Y. Zuo, and B. Cheng, Photonics Research 10, 1567 (2022).
  * Luo  _et al._ [2022] L. Luo, S. Assali, M. R. M. Atalla, S. Koelling, A. Attiaoui, G. Daligou, S. Martí, J. Arbiol, and O. Moutanabbir, ACS Photonics 9, 914 (2022).
  * Marzban  _et al._ [2022] B. Marzban, L. Seidel, T. Liu, K. Wu, V. Kiyek, M. H. Zoellner, Z. Ikonic, J. Schulze, D. Grützmacher, G. Capellini, M. Oehme, J. Witzens, and D. Buca, ACS Photonics , acsphotonics.2c01508 (2022).
  * Talamas Simola  _et al._ [2021] E. Talamas Simola, V. Kiyek, A. Ballabio, V. Schlykow, J. Frigerio, C. Zucchetti, A. De Iacovo, L. Colace, Y. Yamamoto, G. Capellini, D. Grützmacher, D. Buca, and G. Isella, ACS Photonics 8, 2166 (2021).
  * Tran  _et al._ [2019] H. Tran, T. Pham, J. Margetis, Y. Zhou, W. Dou, P. C. Grant, J. M. Grant, S. Al-Kabi, G. Sun, R. A. Soref, J. Tolle, Y.-H. Zhang, W. Du, B. Li, M. Mortazavi, and S.-Q. Yu, ACS Photonics 6, 2807 (2019).
  * Xu  _et al._ [2019] S. Xu, W. Wang, Y.-C. Huang, Y. Dong, S. Masudy-Panah, H. Wang, X. Gong, and Y.-C. Yeo, Optics Express 27, 5798 (2019).
  * Zhou  _et al._ [2020] Y. Zhou, Y. Miao, S. Ojo, H. Tran, G. Abernathy, J. M. Grant, S. Amoah, G. Salamo, W. Du, J. Liu, J. Margetis, J. Tolle, Y.-h. Zhang, G. Sun, R. A. Soref, B. Li, and S.-Q. Yu, Optica 7, 924 (2020).
  * Julsgaard  _et al._ [2020] B. Julsgaard, N. von den Driesch, P. Tidemand-Lichtenberg, C. Pedersen, Z. Ikonic, and D. Buca, Photon. Res. 8, 788 (2020).
  * Vitiello  _et al._ [2020] E. Vitiello, S. Rossi, C. A. Broderick, G. Gravina, A. Balocchi, X. Marie, E. P. O’Reilly, M. Myronov, and F. Pezzoli, Phys. Rev. Applied 14, 064068 (2020).
  * Hudait  _et al._ [2022] M. K. Hudait, S. W. Johnston, M. B. Clavel, S. Bhattacharya, S. Karthikeyan, and R. Joshi, Journal of Materials Chemistry C 10, 10530 (2022).
  * Saleh [2019] B. E. A. Saleh,  _Fundamentals of Photonics_ , Wiley Series in Pure and Applied Optics (Wiley, 2019).
  * Lasher and Stern [1964] G. Lasher and F. Stern, Phys. Rev. 133, A553 (1964).
  * Wurfel [1982] P. Wurfel, Journal of Physics C: Solid State Physics 15, 3967 (1982).
  * Katahara and Hillhouse [2014] J. K. Katahara and H. W. Hillhouse, Journal of Applied Physics 116, 173504 (2014).
  * Chuang [2009] S. L. Chuang,  _Physics of Photonic Devices_ , 2nd ed., Wiley Series in Pure and Applied Optics (Wiley, 2009).
  * Landsberg [1992] P. T. Landsberg,  _Recombination in Semiconductors_, 1st ed. (Cambridge University Press, 1992).
  * Schubert [2006] E. F. Schubert,  _Light-Emitting Diodes_, 2nd ed. (Cambridge University Press, 2006).
  * Stange [2019] D. Stange,  _Group IV (Si)GeSn light emission and lasing studies_ , Dissertation, RWTH Aachen University, Jülich (2019), druckausgabe: 2019. - Onlineausgabe: 2019. - Auch veröffentlicht auf dem Publikationsserver der RWTH Aachen University; Dissertation, RWTH Aachen University, 2019\.
  * Wirths [2016] S. Wirths,  _Group IV epitaxy for advanced nano- and optoelectronic applications_ , Dissertation, RWTH Aachen, Jülich (2016), druckausgabe: 2016. - Onlineausgabe: 2016. - Auch veröffentlicht auf dem Publikationsserver der RWTH Aachen University; Dissertation, RWTH Aachen, 2015.
  * Dijkstra [2021] A. Dijkstra,  _Optical properties of direct band gap group IV semiconductors_ , Ph.D. thesis, Applied Physics and Science Education (2021), proefschrift. - Embargo 1 year, pdf open access 28-5-2022.
  * Bahder [1990] T. B. Bahder, Physical Review B 41, 11992 (1990).
  * Chang  _et al._ [2010] G.-E. Chang, S.-W. Chang, and S. L. Chuang, IEEE Journal of Quantum Electronics 46, 1813 (2010).
  * Lu Low  _et al._ [2012] K. Lu Low, Y. Yang, G. Han, W. Fan, and Y.-C. Yeo, Journal of Applied Physics 112, 103715 (2012).
  * Polak  _et al._ [2017] M. P. Polak, P. Scharoch, and R. Kudrawiec, Journal of Physics D: Applied Physics 50, 195103 (2017).
  * Bir and Pikus [1974] G. L. Bir and G. E. Pikus,  _Symmetry and strain-induced effects in semiconductors_ (Wiley, 1974).
  * Szmulowicz [1995] F. Szmulowicz, Physical Review B 51, 1613 (1995).
  * Note [1] See Supplemental Material for more details about the momentum matrix elements and the special points approximation.
  * Enders [1996] P. Enders, Semiconductor Science and Technology 11, 187 (1996).
  * Schubert  _et al._ [1984] E. F. Schubert, E. O. Göbel, Y. Horikoshi, K. Ploog, and H. J. Queisser, Physical Review B 30, 813 (1984).
  * Ouadjaout and Marfaing [1990] D. Ouadjaout and Y. Marfaing, Physical Review B 41, 12096 (1990).
  * Asada [1989] M. Asada, IEEE Journal of Quantum Electronics 25, 2019 (1989).
  * Chow  _et al._ [1997] W. W. Chow, S. W. Koch, and M. Sargent, Semiconductor-laser physics (Springer-Verlag, Berlin ; New York, 1997) pp. 94–95.
  * Stern [1976] F. Stern, Journal of Applied Physics 47, 5382 (1976).
  * Su  _et al._ [1984] C. B. Su, R. Olshansky, J. Manning, and W. Powazinik, Applied Physics Letters 44, 732 (1984).
  * Olshansky  _et al._ [1984] R. Olshansky, C. Su, J. Manning, and W. Powazinik, IEEE Journal of Quantum Electronics 20, 838 (1984).
  * Assali  _et al._ [2021] S. Assali, A. Dijkstra, A. Attiaoui, E. Bouthillier, J. Haverkort, and O. Moutanabbir, Physical Review Applied 15, 024031 (2021).
  * Atalla  _et al._ [2021b] M. R. M. Atalla, S. Assali, A. Attiaoui, C. Lemieux-Leduc, A. Kumar, S. Abdi, and O. Moutanabbir, Advanced Functional Materials 31, 2006329 (2021b).
  * Feldmann  _et al._ [1987] J. Feldmann, G. Peter, E. O. Göbel, P. Dawson, K. Moore, C. Foxon, and R. J. Elliott, Phys. Rev. Lett. 59, 2337 (1987).
  * Bellessa  _et al._ [1998] J. Bellessa, V. Voliotis, R. Grousson, X. L. Wang, M. Ogura, and H. Matsuhata, Phys. Rev. B 58, 9933 (1998).
  * ’t Hooft  _et al._ [1985] G. ’t Hooft, M. Leys, and H. Talen-v.d. Mheen, Superlattices and Microstructures 1, 307 (1985).
  * Rödl  _et al._ [2019] C. Rödl, J. Furthmüller, J. R. Suckert, V. Armuzza, F. Bechstedt, and S. Botti, Phys. Rev. Mater. 3, 034602 (2019).


# Identifying latent classes with ordered categorical indicators

R. Noah Padgett
Department of Educational Psychology
Baylor University
Waco, TX 76706
<EMAIL_ADDRESS>Rebecca J. Tipton
Department of Educational Psychology
Baylor University
Waco, TX 76706

###### Abstract

A Monte Carlo simulation was used to determine which assumptions for ordered
categorical data, continuity vs. discrete categories, most frequently
identifies the underlying factor structure when a response variable has five
ordered categories. The impact of infrequently endorsed response categories
was also examined, a condition that has not been fully explored. The typical
method for overcoming infrequently endorsed categories in applied research is
to collapse response options with adjacent categories resulting in less
response categories that are endorsed more frequently, but this approach may
not necessarily provide useful information. Response category endorsement
issues have been studied in Item Response Theory, but this issue has not been
addressed in classification analyses nor has fit measure performance been
examined under these conditions. We found that the performance of commonly
used fit statistics to identify the true number of latent class depends on the
whether continuity is assumed, sample size, and convergence. Fit statistics
performed best when the five response options are assumed to be categorical.
However, in situations with lower sample sizes and when convergence is an
issue, assuming continuity and using the adjusted Lo-Mendell-Rubin likelihood
ratio test may be useful.

Classification is an integral part of modern society. The presence of
classification systems is ubiquitous; ranging from how groceries are organized
in a store to grouping children by ability levels in a classroom.
Classification, or clustering, of individuals in a population is oftentimes
required so treatments or programs can be implemented for identified
subgroups. Classification analyses aim to bring a more objective stance into
such groupings. The utility of interventions is dependent on accurate group
identification through appropriate use of classification analysis, thus
choosing the appropriate models and discarding those that are not a good fit
or representation of the sample is important [16, 5]. Many procedures and
frameworks exist to assess group belonging and how well a generated model, or
theoretical mathematical representation of relationships in a population, fits
the observed sample. One such framework is mixture modeling. Mixture modeling
works under the premise that observed data are from a heterogeneous population
and the modeling framework aims to identify groups of cases that are similar
based on observed characteristics into more homogeneous subgroups. The
resulting homogeneous populations are the mixtures that comprise the observed
data [13].

Two such analyses that aim to find homogeneous subgroups from data that are
assumed to be from a heterogeneous population are latent class analysis (LCA)
and latent profile analysis (LPA). These two analyses have the same goal; to
classify cases (individuals) into a small number of groups. Although, each
analysis makes different assumptions about the observed data, LCA assumes the
indicators are categorical while LPA assumes the indicators are continuous.
Prior methodological investigations have examined the effect of nonnormal
indicators have on model selection and particularly on the deciding the number
of latent classes [13, 15].

A type of nonnormally distributed indicator is when a response category is
infrequently selected; a condition that has not been fully explored. Consider
an item where all the respondents answer a question such that no missing data
exist for the item, but the responses could still likely be skewed and one of
the categories could be unendorsed or infrequently endorsed by chance or by
nature of the sample. For example, the resulting frequency of item category
responses would contain a category without sufficient endorsement in the
sample (i.e. sample size = 500; SD = 25, D = 0, N = 50, D = 200, SD = 225).
The typical method for overcoming this occurrence in applied research is to
collapse the response options with an adjacent category, resulting in a more
parsimonious solution but not necessarily one that is as representative or
useful. Collapsing categories is a simple solution to help with estimation of
parameters to make the number of categories one less. Infrequency of response
category endorsement is discussed in Item Response Theory, where the latent
variable of interest is continuous. “If a response category has zero or few
responses, a program to estimate parameters will not be able to estimate
[graded response model parameters]” [8]. Although, the issue has not been
fully explored in classification analyses nor has fit measure performance been
examined under these conditions.

In this study, we investigate what happens when ordered response categories
are assumed continuous and mixture modeling is employed. The assumption of
continuity is expected to aid in model selection when a response category is
unendorsed or infrequently endorsed. Furthermore, we investigate the
performance of fit statistics used in model selection under these conditions.
The remainder of our discussion structures as follows. We discuss the LCA/LPA
model considered in more detail, followed by an brief introduction to the
statistical fit measures for model selection included. We then describe the
Monte Carlo simulation study conditions. The results are in the third section
followed by a discussion that includes recommendations for use of statistical
fit indices for model selection.

### LCA and LPA

Classification analyses such as Latent Class Analysis (LCA) and Latent Profile
Analysis (LPA) are two types of mixture models used for classification. These
models are receiving growing attention in methodological research due to the
expanding application of such models [26, 14, 18, 22]. The goal of LCA and LPA
alike is to classify similar objects into one of K groups or classes of
unknown form and frequency, with the form of the group referring to cluster-
specific centroids, variances, and covariance, and the frequency referring to
the number of underlying groups present.

The model is briefly introduction below; however a more in-depth discussion of
the LCA/LPA models can be found in [13]. Let $\mathbf{y}_{i}$ represent a
vector of responses of individual $i$ on the $j^{th}$ item, where
$i=1,2,...,N$ for the number of individuals sampled and $j=1,...,K$ for $K$
being the number of latent classes specified. LCA is a probabilistic model for
each unique response pattern observed in a sample and is defined as

$P(\mathbf{y}_{i})=\sum_{k=1}^{K}{P(X=k)}\ \prod_{j=1}^{J}{P(y_{ij}\ |X=k)}$
(1)

where $P(\mathbf{y}_{i})$ is the probability of response pattern for
individual $i$, $P(X=k)$ is the probability of membership in class $k$ also
known as class prevalence which is the proportional size of class $k$, and
$P(y_{ij}|X=k)$ is the probability of the response of individual $i$ on the
$j^{th}$ item conditional on class membership. The difference between LCA and
LPA is in the form of the item response probability. In LCA, each class has
specific item response/endorsement probabilities that differential the
classes, whereas in LPA each class has specific response means and variances.

### Model Selection Assumptions

Models, by definition, are only approximations to unknown reality or truth.
George Box made the famous statement, “All models are wrong but some are
useful.” Models that perfectly reflect reality do not exist, but selecting
models that best represent relationships in the population can still be useful
in decision-making. Three general guiding principles that should be considered
during model selection are (a) parsimony, (b) multiplicity or alternative
model hypotheses, and (c) strength of evidence. Inference under models with
too few parameters can be biased, but models having too many parameters
results in poor precision and spurious conclusions. A balance must exist
between under- and over-specification during the model selection process. As
such, several working hypotheses or alternate models must be considered as
potential representations of the true population based on theory and previous
research, and assessed accordingly, but the number should be kept small and
should consider scientific context. Judgment on the part of the researcher is
critical during the hypothesis testing stage, and considering multiple sources
and strength of evidence collected in concert with previous research findings
is important in model selection.

### Fit Indices for Model Selection

Once models have been estimated, the researcher is tasked with selecting the
optimal fitting model. The model with the correct number of latent classes is
selected using a mix of statistical evidence and guiding theory. There are no
established guidelines to determine adequate fit of a theoretical model to
data, but ensuring that parameter estimates are within reasonable range and
standard errors for the estimates are not too large is an important component
of considering model fit [12]. The fit indices included in a study are
suggested to be chosen based upon the study being conducted, but cannot be
evaluated independently of one another [24]. In this article, we examine the
methods for determining fit using two major procedures of model selection
statistical using information criteria and likelihood ratio based tests.

#### Information Criteria Procedures

Measures that use the information criteria are of particular importance
because they provide a framework for comparing models with differing numbers
of parameters and different class enumerations. In general, the model that has
the lowest information criteria in terms of absolute value, is the model that
best approximates the relationships observed in the population, and models
that do not fit better than the baseline model can be dismissed. For this
study, we have chosen to assess Akaike’s Information Criterion (AIC),
corrected AIC (AICc), Bayesian Information Criterion (BIC), and sample size-
adjusted Bayesian Information Criterion (ssBIC). These criteria are sensitive
to parameters like sample size and parsimony, and model selection should be
informed by biases innate in the criteria selected. To date, there is not
common acceptance of the best criteria for determining the number of classes
in mixture modeling, despite various suggestions.

To see how well these models fit collected data, the log-likelihood value is
typically used. Due to the likelihood function taking on small values and the
resulting log likelihood being a negative value, a value close to 0 indicates
optimal fit. A log-likelihood closer to 0 is the same as the likelihood
function approaching 1, which would indicated that the model predicts these
data well. Prior examination of correct class enumeration based on the log-
likelihood values, however, has shown poor results [18, 14].

The Akaike Information Criterion (AIC) is a goodness-of-fit measure that
reflects the extent to which the observed covariance matrix varies from the
model or predicted covariance matrix [1, 2], with a lower AIC value for two
competing models indicating better model fit (Kaplan, 2000; Takane & Bozdogan,
1987). The AIC is defined as

$AIC=-2\log L+2p$ (2)

Where $L$ is value from the likelihood function and $p$ is the number of free
model parameters. The AIC value is penalized for complexity of the proposed
model, meaning use of the AIC is an attempt to minimize the overall error
caused by added parameters. The corrected AIC (AICc) is biased corrected
version of the AIC under small sample sizes [23, 6]. The AICc is considered to
be more stringent than the AIC, having a greater penalty for models having
larger numbers of parameters.

$AICc=-2\log L+2p+\frac{2p(p+1)}{n-p-1}$ (3)

The Bayesian Information Criteria (BIC) is similar to AIC in that it is also a
measure of goodness-of-fit as well as the extent to which the observed
covariance matrix varies from the predicted covariance matrix [21]. The BIC is
defined as

$BIC=-2\log L+p\log(n)$ (4)

A lower BIC value for two competing models indicates a better model fit; the
penalty of BIC increases with sample size. To accommodate for this penalty,
the sample size-adjusted Bayesian Information Criterion (ssBIC) replaces n in
the above equation with $n^{*}=(n+2)/24$

$ssBIC=-2\log L+p\log[(n^{*}+2)/24]$ (5)

For each of the information criteria, the optimal fitting model is chosen
based on which of the compared model have the lowest value.

#### Likelihood Ratio Tests

Other measures of fit include the Lo-Mendell-Rubin likelihood ratio tests. We
focus our attention on the adjusted Lo-Mendell-Rubin Likelihood Ratio Test
(aLMR) [11, 25]. The likelihood ratio test is able to compare models that
differ in the number of classes by indicating that the model with K-1 classes
should be rejected in favor of the model with K classes (Lo, Mendell, & Rubin,
2001; Muthen, 2004). The aLMR is sample-size dependent, meaning that a larger
sample size inflates the test statistics [11].

Bootstrapping procedures can also be used to supplement the likelihood ratio
test family by generating and using empirical distributions of likelihoods
[13]. Comparison of models using a likelihood-based technique can be done by a
parametric bootstrapping method where bootstrap samples are used to estimate
differences in the distribution and a p-value allows for comparison of the K-1
and K class models. The bootstrap likelihood ratio test (BLRT) can then be
used to determine if two competing models are significantly different, as
described in detail by [18, 4]. It has been suggested that the BLRT may be
more appropriate than LMR statistics when comparing LCA models due to BLRT’s
accuracy and tendency to produce nonsignificant findings for models with
increased K, whereas LMR findings tend to show significance, then
insignificance, and then significance again as k increases.

Both the aLMR and BLRT procedures will be generated in this study to provide a
value for assessment of whether the proposed K-1 class model should be
rejected in favor of the K class model. The K-1 class is selected as optimal
if the aLMR p-value for the test of K-1 to K is found to be greater than .05,
meaning that the K class does not provide a better fit than the K-1 class.

## Methods

### Population Model Structure

Conditions were selected that mirror empirical research situations as closely
as possible. The population structures for this study were replicated from
[14], where Morgan reviewed educational databases for studies using LCA or
similar mixture modeling techniques. The three factors varied were sample
size, class prevalence, and category response distribution. The included
sample sizes are 500, 1000, and 1500.

#### Class Prevalence

The class prevalence or class size is representative of various sizes of the
latent populations. The underlying populations may be approximately equal in
size or a class may be more prevalent than the others. The existence of a rare
class might also be speculated based on guiding theory; for example, in
psychiatric disorders the class exhibiting the disorder is likely a small
portion of the general population. Varying class prevalence changes the
probability of accurately defining a class, and classes that are rare are
harder to define or statistically justify. Three underlying classes are
assumed in our study and the class prevalences considered are 0.45-0.40-0.15,
0.59-0.26-0.15, and 0.89-0.08-0.03.

#### Category Response Distribution

The main scope of our investigation is to assess the category response
distribution when obtained data are ordered. We assumed the collection of data
from ten items each have five ordered categories. Response probability is the
primary method for classifying people based on like responses, so if a
category is unendorsed or infrequently endorsed, classification becomes more
difficult because less variability exists in response patterns. In our study,
three different conditions are compared. The first condition has a
distribution that contains an unendorsed response in the second the lowest
category. The second condition contains a category that is infrequently
endorsed with only 1% of subjects in each class endorsing that category. The
last condition is a control condition where each pf the categories are
endorsed in a normal distribution. The distribution of class category response
probabilities for the indicators are given in Table 1.

Table 1: Distribution of class response probabilities per condition

Response |  | Response Categories
---|---|---
Frequency1 | Class | SD | D | N | A | SA
Unendorsed | 1 | .05 | .00 | .05 | .10 | .80
| 2 | .05 | .00 | .75 | .15 | .05
| 3 | .90 | .00 | .05 | .025 | .025
Infrequent | 1 | .049 | .01 | .05 | .10 | .80
| 2 | .05 | .01 | .75 | .149 | .05
| 3 | .90 | .01 | .049 | .025 | .025
Control | 1 | .025 | .025 | .05 | .20 | .70
| 2 | .05 | .10 | .70 | .10 | .05
| 3 | .70 | .20 | .05 | .025 | .025

  * Note. 1condition designating how endorsed each item response is across conditions; 1) a response category is not endorsed, 2) infrequently endorsed response category, and 3) control with all response options selected.

### Summarizing Results

Mirroring real-life applied research problems when trying to identify the
unknown number of classes, we fit one through five class solutions for each
type of analysis. Model estimation was conducted in Mplus 8.1 using the
MIXTURE option [17]. The optimal number of classes for each condition was
based on the decision rule for each fit measure. For each fit measure, if a
three class solution is found to be optimal then the solution is flagged as
correct with a one while if another other class enumeration is found to be
optimal then the solution is flagged as incorrect with zero. We summarize the
proportion of replications that each fit indice correctly identified the
number of latent classes. Using logistic regression, we analyzed which
condition(s) and/or interaction of conditions had the largest effect on model
selection. The logistic regression model tested uses the three factors
described above and the analysis type (LCA vs. LPA) as indicators. [9] suggest
the use of inferential statistics in simulation to aid in result
interpretation. All analyses after model estimation was completed using R
[20].

## Results

The entirety of the results of this simulation study are available online
[19]. A total of 27,000 data sets were successfully generated, 27(3x3x3)
conditions by 1,000 replications. For each data set one- through five- class
solutions were tested using categorical and continuous indicator assumptions,
i.e. LCA and LPA. Testing these models for each of 1000 replicates in each of
27 conditions yields 270,000 models were estimated. Convergence was found to
be problematic when assuming categorical indicators, especially in the
infrequently endorsed category conditions. Models that did not converge were
flagged and the replication was deemed unusable, an approach taken in similar
simulation studies, and nonconvergence is typically due to an ill-defined
likelihood function and/or an insufficient number of random starts [18, 14].
The convergence rates for each condition are provided in Table 2.

Table 2: Convergence rate across replications Response | Class | N=500 | N=1000 | N=1500
---|---|---|---|---
Frequency | Prevalence | LCA | LPA | LCA | LPA | LCA | LPA
Unendorsed | .45-.40-.15 | 1.000 | 1.000 | .990 | 1.000 | .967 | 1.000
| .59-.26-.15 | .998 | 1.000 | .994 | 1.000 | .970 | 1.000
| .89-.08-.03 | .999 | 1.000 | .979 | 1.000 | .898 | 1.000
Infrequent | .45-.40-.15 | .996 | 1.000 | .935 | 1.000 | .826 | 1.000
| .59-.26-.15 | .998 | 1.000 | .967 | 1.000 | .875 | 1.000
| .89-.08-.03 | .996 | 1.000 | .938 | 1.000 | .736 | 1.000
Control | .45-.40-.15 | .999 | 1.000 | .981 | 1.000 | .931 | 1.000
| .59-.26-.15 | .999 | 1.000 | .985 | 1.000 | .949 | 1.000
| .89-.08-.03 | .999 | 1.000 | .958 | 1.000 | .801 | 1.000

### Information Criteria

In this study, four different information criteria were examined for model
selection, the Akaike information criterion (AIC), corrected Akaike
information criterion (AICc), Bayesian information criterion (BIC), and sample
size-adjusted (ssBIC). The proportion of valid replications that each index
correctly identified the number of latent classes are reported in Table 3. The
results from the log-likelihood statistic are not reported because under all
conditions the correct number of latent classes was not found.

The difference in correct enumeration selection is obvious between the
analyses for all information criteria considered. By assuming the indicators
are continuous, the four information criteria are typically unable to find the
correct number of latent classes under these conditions. The only exception is
the BIC under the control condition when the response distribution is fairly
normally distributed; yet, the identification rate is still low (see Table 3).
When the assumption of categorical indicators is used, the rate of correct
specification increases remarkably. For example, the ssBIC perfectly
identified the number of latent classes in all conditions under the assumption
of categorical indicators, but under the assumption of continuity, the ssBIC
performed poorly (see Table 3). These results lend evidence to how sensitive
the information criteria are to the assumption of continuity.

Table 3: Proportion of data sets each fit measure correctly identified a three
class solution as optimal by analysis type

|  |  | AIC | AICc | BIC | ssBIC | aLMR | BLRT
---|---|---|---|---|---|---|---|---
RF1 |  CP2 | N | LCA | LPA | LCA | LPA | LCA | LPA | LCA | LPA | LCA | LPA | LCA | LPA
U | 1 | 500 | .75 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .58 | .69 | .69 | .00
|  | 1000 | .74 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .43 | .27 | .61 | .00
|  | 1500 | .72 | .00 | .98 | .00 | 1.00 | .00 | 1.00 | .00 | .35 | .12 | .61 | .00
| 2 | 500 | .72 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .70 | .64 | .69 | .00
|  | 1000 | .70 | .00 | .99 | .00 | 1.00 | .00 | 1.00 | .00 | .47 | .15 | .62 | .00
|  | 1500 | .67 | .00 | .96 | .00 | 1.00 | .00 | 1.00 | .00 | .37 | .03 | .58 | .00
| 3 | 500 | .75 | .00 | 1.00 | .00 | .98 | .00 | 1.00 | .00 | .66 | .45 | .83 | .00
|  | 1000 | .70 | .00 | .99 | .00 | 1.00 | .00 | 1.00 | .00 | .49 | .07 | .76 | .00
|  | 1500 | .72 | .00 | .97 | .00 | 1.00 | .00 | 1.00 | .00 | .34 | .01 | .73 | .00
I | 1 | 500 | .86 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .01 | .69 | .67 | .00
|  | 1000 | .83 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .25 | .31 | .69 | .00
|  | 1500 | .81 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .56 | .12 | .70 | .00
| 2 | 500 | .90 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .07 | .66 | .67 | .00
|  | 1000 | .82 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .48 | .20 | .66 | .00
|  | 1500 | .76 | .00 | .99 | .00 | 1.00 | .00 | 1.00 | .00 | .77 | .04 | .60 | .00
| 3 | 500 | .89 | .00 | 1.00 | .00 | .96 | .00 | 1.00 | .00 | .09 | .46 | .78 | .00
|  | 1000 | .83 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .25 | .13 | .79 | .00
|  | 1500 | .80 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .38 | .02 | .76 | .00
C | 1 | 500 | .63 | .00 | 1.00 | .00 | 1.00 | .45 | 1.00 | .00 | .58 | .71 | .87 | .00
|  | 1000 | .65 | .00 | 1.00 | .00 | 1.00 | .04 | 1.00 | .00 | .78 | .43 | .83 | .00
|  | 1500 | .66 | .00 | .99 | .00 | 1.00 | .00 | 1.00 | .00 | .92 | .25 | .80 | .00
| 2 | 500 | .64 | .00 | 1.00 | .00 | 1.00 | .14 | 1.00 | .00 | .72 | .67 | .84 | .00
|  | 1000 | .62 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .89 | .21 | .81 | .00
|  | 1500 | .62 | .00 | .99 | .00 | 1.00 | .00 | 1.00 | .00 | .96 | .07 | .78 | .00
| 3 | 500 | .65 | .00 | 1.00 | .00 | .88 | .11 | 1.00 | .00 | .43 | .57 | .91 | .00
|  | 1000 | .64 | .00 | 1.00 | .00 | 1.00 | .00 | 1.00 | .00 | .48 | .16 | .90 | .00
|  | 1500 | .66 | .00 | .99 | .00 | 1.00 | .00 | 1.00 | .00 | .57 | .03 | .86 | .00

  * Note.Cells in boldface indicate that LPA correctly specified the correct number of class more often than LCA. LCA = latent class analysis; LPA = latent profile analysis; AIC = Akaike information criterion; AICc = corrected Akaike information criterion; BIC = Bayesian information criterion; ssBIC = sample size-adjusted Bayesian information criterion; aLMR = adjusted Lo - Mendell - Rubin likelihood ratio test; BLRT = bootstrap likelihood ratio test.

  * 1

RF is the response frequency condition described in Table 1, where U =
unendorsed category, I = infrequently endorsed category, and C = control
condition;

  * 2

CP is class prevalence condition, 1) .45, .40, .15; 2) .59, .26, .15; 3) .89,
.08, .03.

### Likelihood Ratio Tests

Two likelihood ratio tests for model selection were examined, the adjusted Lo-
Mendell-Rubin likelihood ratio test (aLMR) and the bootstrap likelihood ratio
test (BLRT). Both likelihood ratio tests performed worse than the information
criteria included. The assumption of continuity had the most effect on the
performance of the BLRT, where the BLRT did not identify the correct number of
latent classes under all conditions (see Table 3). Although, the BLRT
performed consistently when the indicators are assumed to be categorical. The
stark difference of BLRT performance between LCA and LPA is not found with the
aLMR.

The performance of aLMR is more nuanced in the studied conditions compared to
the other fit measures. We examined the performance of the aLMR with logistic
regression to help identify which factors and levels included in our study are
influential in determining correct specification. Notable results are evidence
of an interaction between analysis type-LPA and response distribution-
infrequently endorsed category (OR = 124.64); and an interaction among the
analysis type (LPA), response distribution-infrequently endorsed category,
class prevalence (.89 - .08 - .03) and sample size-1000 (OR = 8.17).

The complex interactions lend evidence that the aLMR performs better under
some conditions and when different assumptions of the continuity of the
indicators are made. These interaction are difficult to make sense of from
looking at the values reported in Table 3, so the proportion of replications
where the aLMR identified a three class solution as optimal is plotted by each
condition and analysis type in Figure 1 to aid interpretation. The assumption
of continuity does appear to help when a response category is infrequently
endorsed and sample size is small, but as sample size increases the assumption
of continuity does not appear to help with model selection.

Figure 1: Complex interactions of the aLMR results for estimating the correct
number of latent classes

Note. The solid line represents LCA while the dashed line represents LPA
results. The vertical panels are for each condition of class prevalence and
the horizontal panels are for the response frequency conditions.

## Discussion

The use of classification analyses is multidisciplinary so determining proper
uses of available techniques and how to evaluate obtained results is crucial.
The aim of our investigation was to find evidence of what occurs in model
selection when ordered response categories are assumed to be continuous. We
studied what statistical measures of model fit are useful when data are
ordered response categories and when a response category is unendorsed or
infrequently endorsed under the different assumptions of continuity of
indicators.

The use of ordered response categories can lead to possible extraction issues.
When ordered categories are assumed to be continuous, extraction of an
incorrect number of latent classes is more likely. These ordered categorical
data are highly prevalent in applied research in education and psychology, as
found in a review by [10]. An accepted guideline in methodological research is
that these data can be treated as continuous if there are at least five
categories. However, in conditions studied here for mixture modeling we found
evidence that five ordered response categories should be treated as ordered
categorical and not continuous. The assumption of continuity only seems
beneficial when convergence is problematic, such as when a response category
is infrequently endorsed. When convergence is not an issue, our results are
similar to previous simulation studies of these fit indices. The BIC and ssBIC
were found to consistently select the correct number of latent classes, a
result similar to other studies [14, 18, 6, 22]. The BLRT was also found to
perform poorly, which is counter to what [18] found but is similar to the
results of [7].

These statistical measures of model selection appears to be sensitive to the
assumption of continuity and the frequency of response category endorsement.
This implies that there are no strict rules or guidelines for use of these fit
indices, and researchers should be aware of possible issues of relying on
these indices too heavily. Despite the fact that these measures can be useful,
performance of these fit indices changes based on the distribution pf the
indicators and the assumptions of the researcher.

### Limitations and Delimitations

As with any simulation study, these results generalize only to the conditions
included in this study. We used default settings in Mplus 8.1 to try to
generalize as widely as possible. However, additional settings in Mplus could
help with issues that arise with estimation. For example, Bayesian estimation
procedures are a promising avenue [3]. The included statistical fit indices
are only a subset of possible choices available to practitioners; therefore
the generalizations of our investigation for model selection is limited to
these included measures.

### Recommendations for Use of Fit Indices in LCA & LPA

We recommend that when data collected are five ordered categories, such as
Likert-type data, and are being used for identification of homogeneous
subpopulations, these data should be considered categorical and not
continuous. Despite the advantage of convergence, the commonly available fit
statistics for model selection performed poorly under the assumption of
continuity with these data. One exception is the adjusted Lo-Mendell-Rubin
likelihood ratio statistic. The aLMR is the only fit measure examined in this
study to identify the correct number of latent classes when a response
category is infrequently or not endorsed and the indicators are assumed to be
continuous. The use of aLMR could help practitioners when model convergence is
an issue and sample size is small. Although, under most conditions we
recommend that the corrected AIC (AICc), Bayesian information criterion (BIC),
or the sample size-adjusted BIC (ssBIC) be used to aid model selection when
the indicators are five ordered categories along with substantial weight to
substantive theory and interpretation of results.

### Author Contributions

R.N. Padgett and R. J. Tipton jointly generated the idea for the study. R.N.
Padgett wrote the code for generating the data, estimating the models,
extracting the relevant results and analyzing the model summaries. R. J.
Tipton verified the accuracy of those analyses. R.N. Padgett and R. J. Tipton
collaboratively wrote the first draft of the manuscript, and both authors
critically edited it. Both authors approved the final submitted version of the
manuscript.

## References

  * [1] Hirotugu Akaike  “A New Look at the Statistical Model Identification”  In _IEEE Transactions on Automatic Control_ 19.6, 1974, pp. 716–723  DOI: 10.1109/TAC.1974.1100705
  * [2] Hirotugu Akaike  “Factor analysis and AIC”  In _Psychometrika_ 52.3, 1987, pp. 317–332  DOI: 10.1007/BF02294359
  * [3] Tihomir Asparouhov, Ellen L. Hamaker and Bengt Muthén  “Dynamic Latent Class Analysis”  In _Structural Equation Modeling_ 24.2, 2017, pp. 257–269  DOI: 10.1080/10705511.2016.1253479
  * [4] Tihomir Asparouhov and Bengt O Muthén  “Using Mplus TECH11 and TECH14 to test the number of latent classes”  In _Mplus Web Notes_ , 2012, pp. 1–17
  * [5] Leo Breinman  “Random forests”  In _Machine Learning_ 45.1, 2001, pp. 5–32  DOI: 10.1186/1478-7954-9-29
  * [6] Kenneth P Burnham and David R Anderson  “Multimodel inference: understanding AIC and BIC in model selection”  In _Sociological methods & research_ 33.2, 2004, pp. 261–304
  * [7] John J Dziak, Stephanie T Lanza and Xianming Tan  “Effect Size, Statistical Power, and Sample Size Requirements for the Bootstrap Likelihood Ratio Test in Latent Class Analysis”  In _Structural Equation Modeling: A Multidisciplinary Journal_ 21.4, 2014, pp. 534–552  DOI: 10.1080/10705511.2014.919819
  * [8] Susan E Embretson and Steven P Reise  “Item response theory for psychologists”  Psychology Press, 2000, pp. 83–85
  * [9] Gregory R. Hancock and Ralph O. Mueller  “Structural Equation Modeling: A Second Course”, Quantitative methods in education and the behavioral sciences  Charlotte: Information Age Publishing, 2013, pp. 625–663
  * [10] Michael R Harwell and Guido G Gatti  “Rescaling Ordinal Data to Interval Data in Educational Research”  In _Review of Educational Research Spring_ 71.1, 2001, pp. 105–131  DOI: 10.3102/00346543071001105
  * [11] Y Lo, Nancy R. Mendell and Donald B. Rubin  “Testing the number of components in a normal mixture”  In _Biometrika_ 88.3, 2001, pp. 767–778  DOI: 10.1093/biomet/88.3.767
  * [12] Herbert W Marsh and David Grayson  “Latent variable models of multitrait-multimethod data”  In _Structural equation modeling: Concepts, issues and applications_ Thousand Oaks, CA: Sage Publications, 1995, pp. 177–198
  * [13] Geoffrey McLachlan and David Peel  “Finite Mixture Models”  John Wiley & Sons, 2000
  * [14] Grant B. Morgan  “Mixed Mode Latent Class Analysis: An Examination of Fit Index Performance for Classification”  In _Structural Equation Modeling: A Multidisciplinary Journal_ 22.1, 2015, pp. 76–86  DOI: 10.1080/10705511.2014.935751
  * [15] Grant B. Morgan, Kari J. Hodge and Aaron R. Baggett  “Latent profile analysis with nonnormal mixtures: A Monte Carlo examination of model selection using fit indices”  In _Computational Statistics and Data Analysis_ 93, 2016, pp. 146–161  DOI: 10.1016/j.csda.2015.02.019
  * [16] Frederick Mosteller and W. Tukey, J.  “Data Analysis and Regression: A Second Course in Statistics”  In _Addison-Wesley Series in Behavioral Science: Quantitative Methods_ , 1977
  * [17] L.K. Muthén and B.O. Muthén  “Mplus User’s Guide”  Los Angeles, CA: Muthén & Muthén, 2017
  * [18] Karen L. Nylund, Tihomir Asparouhov and Bengt O. Muthén  “Deciding on the Number of Classes in Latent Class Analysis and Growth Mixture Modeling: A Monte Carlo Simulation Study”  In _Structural Equation Modeling: A Multidisciplinary Journal_ 14.4, 2007, pp. 535–569  DOI: 10.1080/10705510701575396
  * [19] R. Padgett and Rebecca J Tipton  “Monte Carlo Simulation Results Categorical Latent Class Analysis”  Texas Data Repository Dataverse, 2019  DOI: 10.18738/T8/CHW28Z
  * [20] R Core Team  “R: A Language and Environment for Statistical Computing”, 2017  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [21] Gideon Schwarz  “Estimating the Dimension of a Model”  In _The Annals of Statistics_ 6.2, 1978, pp. 461–464  DOI: 10.1214/aos/1176344136
  * [22] G Soromenho  “Comparing approaches for testing the number of components in a finite mixture model”  In _Computational Statistics_ 9, 1994, pp. 65–78
  * [23] Nariaki Sugiura  “Further analysts of the data by akaike’s information criterion and the finite corrections: Further analysts of the data by akaike’s”  In _Communications in Statistics-Theory and Methods_ 7.1, 1978, pp. 13–26
  * [24] J Tanaka  “Multifaceted conceptions of fit in structural equation models”  In _Testing structural equation models_ Newbury Park, CA: Sage, 1993, pp. 10–40
  * [25] Quang H Vuong  “Likelihood ratio tests for model selection and non-nested hypotheses”  In _Econometrica_ 57.2, 1989, pp. 307
  * [26] Glenn A Williams and Fraenze Kibowski  “Latent Class Analysis and Latent Profile Analysis in Handbook of methodological approaches to community-based research”  New York: Oxford University Press, 2016, pp. 143–151