# Are Cluster Validity Measures (In)valid?

Marek Gagolewski  Maciej Bartoszuk  Anna Cena  Deakin University, School of
Information Technology, Geelong, VIC 3220, Australia  Warsaw University of
Technology, Faculty of Mathematics and Information Science,
ul. Koszykowa 75, 00-662 Warsaw, Poland  Systems Research Institute, Polish
Academy of Sciences
ul. Newelska 6, 01-447 Warsaw, Poland

###### Abstract

Internal cluster validity measures (such as the Caliński–Harabasz, Dunn, or
Davies–Bouldin indices) are frequently used for selecting the appropriate
number of partitions a dataset should be split into. In this paper we consider
what happens if we treat such indices as objective functions in unsupervised
learning activities. Is the optimal grouping with regards to, say, the
Silhouette index really meaningful? It turns out that many cluster
(in)validity indices promote clusterings that match expert knowledge quite
poorly. We also introduce a new, well-performing variant of the Dunn index
that is built upon OWA operators and the near-neighbour graph so that
subspaces of higher density, regardless of their shapes, can be separated from
each other better.

Keywords: clustering methodology, cluster validity index, Dunn index, nearest
neighbours (NNs), ordered weighted averaging (OWA) operator, no free lunch

††journal: Information Sciences

##  1 Introduction

An internal cluster validity index (CVI for short; see, e.g., [43, 42, 29, 2,
58]) is – in theory – a measure of how well a given partitioning of a dataset
reflects the underlying structure of the modelled domain. CVIs are frequently
employed as tools for selecting the appropriate number of clusters a dataset
should be segmented into [43]. By re-applying some algorithm (e.g., $k$-means
or spectral methods), one can determine the splits into 2-, 3-, 4-, …,
disjoint and nonempty subsets, compute the corresponding CVIs, and select the
partition that maximises a chosen utility measure.

Here we shall focus on the other popular use case thereof. Some practitioners
utilise CVIs to compare the outputs of different algorithms on the same
dataset. Is the partition that the average linkage method returned better than
that yielded by the DBSCAN algorithm (provided that they are of equal
cardinality)? Similarly, researchers use CVIs for evaluating new clustering
algorithms the same way: a new method $X^{++}$ produces partitions that have a
higher average Caliński–Harabasz indices (on some benchmark datasets) than
procedures $X$, $Y$, and $Z$, thus “proving” its superiority. However, we
would like to call this methodology into question, especially because CVIs in
general constitute an extremely diverse set of measures.

Thus, we shall be interested in determining which of the popular CVIs are
particularly suitable or unsuitable for judging the quality of different
partitions of the same cardinality. Does a high value of a CVI make sense at
all? Can it really be treated as an indicator of a useful clustering result?

To address these questions, we shall find the partitions that yield the
highest possible index values, for a large number of datasets and CVIs. In
other words, we will treat each CVI as an objective function to be maximised
over the whole space of all possible clusterings.

Our assumption here is that a cluster validity measure can only be considered
meaningful whenever it is maximised at the partitions closely resembling the
reference ones. Otherwise stated, good CVIs should promote results that agree
with expert knowledge.

In the course of our study, which we of course detail in the sequel, we have
discovered that this is often very much not the case – see Figures 1 and 2 for
two quite representative graphical examples. It turns out that some CVIs
promote highly overlapping groupings while other ones work better as outlier
detectors. One should thus not uncritically believe that a high value of,
e.g., a generalised Dunn index GDunn_d2_D1 (see Section 2) is better than a
lower one; at the bottom-right subfigures we see that this index promotes some
rather random partitions as “best”.

Figure 1:  fcps/engytime dataset: Optimal clusters as seen by 6 different
cluster validity indices (see Section 2 for more details). The reference
partitions consist of two (clearly visible if viewed in colour) Gaussian
blobs; the Caliński–Harabasz index (top-left subfigure) identifies them quite
correctly. Some indices promote very peculiar, overlapping groupings (e.g.,
GDunn_d2_D1), other ones should rather be employed as outlier detectors (like
GDunn_d1_D1).

Figure 2:  fcps/wingnut dataset: Optimal clusters as seen by 6 different
cluster validity indices (discussed in Section 2). The two point clouds are
well-separable; GDunn_d1_D1 (top-left subfigure) identified them correctly.
However, certain CVIs favour some rather unusual cluster shapes instead.

Let us note that papers introducing new CVIs are plentiful and new ones are
being published on a regular basis, see, e.g., [40, 39] for some recent
examples. More often than not, such indices are designed to be “the best” for
a given particular situation and/or they aim to “eliminate” certain
deficiencies with the previous measures. Because of this, the number of CVIs
to choose from can be overwhelming, see Section 2 and, e.g., [58, 56], for up-
to-date overviews.

This is why comprehensive, in-depth comparative studies are so important.
There are quite a few interesting surveys of various properties of cluster
validity measures in different contexts and settings, but we have only found
[59] somewhat methodologically relevant to the task at hand. There, the
performance of 8 CVIs was studied but the authors’ focus seems slightly
shifted towards the analysis of the stability of the solutions generated by
what they call differential-evolution–particle-swarm-optimisation algorithms,
which of course do not necessarily guarantee finding patterns that are optimal
in the eye of a specific cluster validity measure.

Further, in [2] and [42], the authors have shown that most of the CVIs work
well with spherical clusters but fail in other types of data. A CVI may fail
to assign the highest evaluation to the partition that fits the data best,
e.g., the reference one. The conclusion is that there is no single internal
cluster validation index that outperforms the other indices everywhere.
Similar conclusions were reached in [7], where a model-based study of the
correlations of CVIs with carefully chosen error rates across the outputs of
different clustering algorithms was conveyed. For more studies in similar
spirit, see [14, 33, 43, 38]. Note again that our task is to find the optimal
partition under the guidance of a given CVI (amongst the set of all
partitions), and not to assess a fixed set of particular clusterings.

We are well aware of the fact that the sheer act of optimising of CVIs is,
overall, not a new idea. For instance, the within-cluster sum of squares
(WCSS), which is the basis for the Caliński–Harabasz index is used as the
objective function in the $k$-means [41] algorithm.

Other algorithms employing some goodness-of-split measures include, e.g.,
fuzzy $c$-means [5] that features a smoothened variant of WCSS, ITM
(Information-Theoretic-MST; [46]) which applies a divisive scheme over an
Euclidean minimum spanning tree to optimise an information-theoretic
criterion, finding Gaussian mixtures via expectation-maximisation, and the
generalisation of the Ward linkage in the form of the Lance–Williams formulae
[35].

Nevertheless, as $n$ data items enjoy $\Omega(k^{n})$ possible
$k$-partitionings (Stirling number of the second kind), the problem of
optimising a general CVI over the whole search space is very difficult
computationally (provably hard in the case of the said WCSS, see [25, 1],
amongst many others).

Because of this, many algorithms which are defined as minimisers/maximisers of
some CVI, use a variety of simplifications, approximations, or heuristics,
e.g., they optimise the objective over reduced search spaces or apply some
greedy strategies. For instance, to find a partition optimising the
aforementioned WCSS, Ward in [57] suggested to build a hierarchy of clusters
in an agglomerative way, and Edwards and Sforza in [19] as well as Caliński
and Harabasz in [8] proposed to employ some divisive schemes. However, such
heuristics are usually limited to specific CVIs; in this paper we would like
to go far beyond that.

Also notice that in the literature we may of course find numerous techniques
constructed as a combination of different metaheuristics-based optimisation
procedures and chosen cluster validity measures playing the role of objective
functions, see, e.g., [30, 34, 52, 61, 48, 13]. Their respective authors often
claim that this way they create “new” clustering algorithms (e.g., optimise
WCSS using particle swarms vs by means of differential evolution), but it is
semantically not quite appropriate. Rather, they should be thought of as ways
to test the performance of the optimisation algorithms themselves (with CVIs
computed over particular datasets serving as benchmark objectives as in [31]).

Due to the intrinsic difficulty of optimising CVIs, it is no wonder that such
a comprehensive study as the current one has not been performed yet.
Interestingly, it will turn out that the methodology we have employed, despite
its still being based on some approximations, is sufficient for achieving our
goal.

Also it is worth mentioning that, up to date, the variety of studies of the
various aspects of CVIs (which we shall review in the next section), was quite
limited because of the lack of a larger, standardised benchmark set batteries.
Most studies considered few datasets, either synthetic, or inherently
difficult to partition (such as the UCI [17] datasets). Luckily, thanks to the
recent notable efforts by the authors of [55, 28, 21] and our new battery that
aggregates and extends them [24], studies such as this one finally become
possible. Furthermore, we propose a unique approach where we account for the
fact that a dataset can exhibit multiple equally valid clusterings (as
discussed also in [11]).

Let us cast a glance at the structure of this paper. In Section 2 we review
some of the most notable cluster validity indices. Furthermore, we propose a
new, wide class of CVIs generalising the Dunn index which is based on the
notion of ordered weighted averaging (OWA) operators and near-neighbour (NN)
graphs. In Section 3 we describe the methodology we have applied in order to
answer our main research question, including the description of the benchmark
datasets used, the approach to identify the partitions that are optimal from
the perspective of a given cluster validity measure, and ways to determine the
extent to which they agree with expert knowledge. In Section 4 we present the
results of our empirical study, e.g., explore the relationship between cluster
compactness or separability and what is considered a good clustering by
experts. Also, we perform a cluster analysis of clustering algorithms to
determine which methods are most similar to each other. We conclude the paper
in Section 5.

##  2 Cluster validity indices

Let $\mathbf{X}\in\mathbb{R}^{n\times d}$ denote the input dataset comprised
of $n$ points in a $d$-dimensional Euclidean space, with
$\mathbf{x}_{i}=(x_{i,1},\dots,x_{i,d})$ denoting the $i$-th point,
$i\in[1:n]=\\{1,2,\dots,n\\}$.

We shall be looking for a partition of $\mathbf{X}$ into $k\geq 2$ nonempty,
mutually disjoint clusters, with $k$ fixed in advance. Note that a
$k$-partition $\\{X_{1},\dots,X_{k}\\}$ of a set
$\\{\mathbf{x}_{1},\dots,\mathbf{x}_{n}\\}$ can be encoded by means of a
surjection $C:[1:n]\stackrel{{\scriptstyle\text{onto}}}{{\to}}[1:k]$, where
$C(i)\in[1:k]$ gives the cluster number of the $i$-th point. Let us denote the
set of all such possible mappings with $\mathcal{C}_{k}$.

For the sake of clarity and simplicity, we will only be focused on cluster
validity indices based on Euclidean distances between all pairs of points,
$\|\mathbf{x}_{i}-\mathbf{x}_{j}\|$, or the input points and some other
pivots, such as their corresponding cluster centroids,
$\|\mathbf{x}_{i}-\boldsymbol{\mu}_{j}\|$, where
$\mu_{j,l}=\frac{1}{|X_{j}|}\sum_{\mathbf{x}_{i}\in X_{j}}x_{i,l}$. The
fixation of the distance metric is not at all restrictive, as various
transformations can be applied onto $\mathbf{X}$ at the data pre-processing
stage, including variable selection, standardisation, outlier removal, feature
engineering (by means of spectral/kernel-based methods), etc., see [15, 44,
16], amongst others.

We shall consider 52 different internal cluster validity indices like
$I:\mathcal{C}_{k}\to\mathbb{R}$. Apart from the most popular, classical CVIs,
we also bring forth our own proposal.

###  2.1 CVIs based on cluster centroids

##### 1,2) BallHall, CalińskiHarabasz

Let $\boldsymbol{\mu}$ denote the centroid of the whole $\mathbf{X}$. The two
following indices are based on within-cluster sum of squares (WCSS), which
itself can be rewritten in terms of the squared Euclidean distances between
the points and their respective centroids.

The Ball–Hall index [3] is the WCSS weighted by the cluster cardinality:

$\mathrm{BallHall}(C)=-\sum_{i=1}^{n}\frac{1}{|X_{C(i)}|}\|\mathbf{x}_{i}-\boldsymbol{\mu}_{C(i)}\|^{2}.$
(1)

Note the minus that accounts for the fact that in Section 3.3 we want all the
indexes be maximised.

Then the Caliński–Harabasz index [8, Eq. (3)] (“variance ratio criterion”) is
given by:

$\mathrm{Cali\'{n}skiHarabasz}(C)=\frac{n-k}{k-1}\frac{\sum_{i=1}^{n}\|\boldsymbol{\mu}-\boldsymbol{\mu}_{C(i)}\|^{2}}{\sum_{i=1}^{n}\|\mathbf{x}_{i}-\boldsymbol{\mu}_{C(i)}\|^{2}}.$
(2)

It may be shown that the task of minimising the (unweighted) WCSS is
equivalent to maximising the Caliński–Harabasz index. Hence, this index is
precisely the objective function in $k$-means [41] and the algorithms by Ward,
Edwards and Cavalli-Sforza, etc., see [57, 19, 8].

##### 3) DaviesBouldin

The Davies–Bouldin [12, Def. 5] index also refers to the notion of cluster
centroids. It is given as the average similarity between each cluster and its
most similar counterpart (note the minus sign again):

$\mathrm{DaviesBouldin}(C)=-\frac{1}{k}\sum_{i=1}^{k}\left(\max_{j\neq
i}\frac{s_{i}+s_{j}}{m_{i,j}}\right),$ (3)

where $s_{i}$ is the dispersion of the $i$-th cluster: if $|X_{i}|>1$, it is
given by $s_{i}=\frac{1}{|X_{i}|}\sum_{\mathbf{x}_{u}\in
X_{i}}\|\mathbf{x}_{u}-\boldsymbol{\mu}_{i}\|$ and otherwise we set
$s_{i}=\infty$. Furthermore, $m_{i,j}$ is the intra-cluster distance,
$m_{i,j}=\|\boldsymbol{\mu}_{i}-\boldsymbol{\mu}_{j}\|$. In [12], other
choices of $s_{i}$ and $m_{i,j}$ are also suggested; here, we choose the most
popular setting (used, e.g., in [2]).

###  2.2 Silhouettes

##### 4, 5) Silhouette, SilhouetteW

In [54, Sec. 2], Rousseeuw proposes the notion of a silhouette as a graphical
aid in cluster analysis.

Denote the average dissimilarity between the $i$-th point and all other points
in its own cluster with:

$a_{i}=\frac{1}{|X_{C(i)}|-1}\sum_{\mathbf{x}_{u}\in
X_{C(i)}}\|\mathbf{x}_{i}-\mathbf{x}_{u}\|$ (4)

and the average dissimilarity between the $i$-th point and all other entities
in the “closest” cluster with:

$b_{i}=\min_{j\neq C(i)}\left(\frac{1}{|X_{j}|}\sum_{\mathbf{x}_{v}\in
X_{j}}\|\mathbf{x}_{i}-\mathbf{x}_{v}\|\right).$ (5)

Then the Silhouette index is defined as the average silhouette score:

$\mathrm{Silhouette}(C)=\frac{1}{n}\sum_{i=1}^{n}\frac{b_{i}-a_{i}}{\max\\{a_{i},b_{i}\\}},$
(6)

with convention $\pm\infty/\infty=0$.

The same paper also defines what we call here the SilhouetteW index, being the
mean of the cluster average silhouette widths:

$\mathrm{SilhouetteW}(C)=\frac{1}{k-s}\sum_{i=1}^{n}\frac{1}{|X_{C(i)}|}\frac{b_{i}-a_{i}}{\max\\{a_{i},b_{i}\\}},$
(7)

where $s$ is the number of singletons. Note that SilhouetteW, just like
BallHall, employs weighting by cluster cardinalities.

###  2.3 Generalised Dunn indices

##### 6–20) GDunn_dX_DY

In [18, Eq. (3)], Dunn proposed an index defined as the ratio between the
smallest between-cluster distance and the largest cluster diameter. It has
been generalised by Bezdek and Pal in [6] as:

$\mathrm{GDunn}(C)=\frac{\min_{i\neq
j}d\left(X_{i},X_{j}\right)}{\max_{i}D\left(X_{i}\right)}.$ (8)

The numerator measures the between-cluster separation whilst the denominator
quantifies the cluster compactness.

Function $d$ was assumed in [6] one of:

  * 1.

$d_{1}(X_{i},X_{j})=\mathrm{Min}\left(\left\\{\|\mathbf{x}_{u}-\mathbf{x}_{v}\|:\mathbf{x}_{u}\in
X_{i},\mathbf{x}_{v}\in X_{j}\right\\}\right)$,

  * 2.

$d_{2}(X_{i},X_{j})=\mathrm{Max}\left(\left\\{\|\mathbf{x}_{u}-\mathbf{x}_{v}\|:\mathbf{x}_{u}\in
X_{i},\mathbf{x}_{v}\in X_{j}\right\\}\right)$,

  * 3.

$d_{3}(X_{i},X_{j})=\mathrm{Mean}\left(\left\\{\|\mathbf{x}_{u}-\mathbf{x}_{v}\|:\mathbf{x}_{u}\in
X_{i},\mathbf{x}_{v}\in X_{j}\right\\}\right)$,

  * 4.

$d_{4}(X_{i},X_{j})=\|\boldsymbol{\mu}_{i}-\boldsymbol{\mu}_{j}\|$,

  * 5.

$d_{5}(X_{i},X_{j})=\frac{|X_{i}|\,\mathrm{Mean}\left(\left\\{\|\mathbf{x}_{u}-\boldsymbol{\mu}_{i}\|:\mathbf{x}_{u}\in
X_{i}\right\\}\right)+|X_{j}|\,\mathrm{Mean}\left(\left\\{\|\mathbf{x}_{v}-\boldsymbol{\mu}_{j}\|:\mathbf{x}_{v}\in
X_{j}\right\\}\right)}{|X_{i}|+|X_{j}|}$.

Bezdek and Pal in [18] considered also $d_{6}$ based on the Hausdorff metric
but this will be omitted here as it turned out too slow to compute.

On the other hand, $D$ was chosen amongst:

  * 1.

$D_{1}(X_{i})=\mathrm{Max}\left(\left\\{\|\mathbf{x}_{u}-\mathbf{x}_{v}\|:\mathbf{x}_{u},\mathbf{x}_{v}\in
X_{i}\right\\}\right)$,

  * 2.

$D_{2}(X_{i})=\mathrm{Mean}\left(\left\\{\|\mathbf{x}_{u}-\mathbf{x}_{v}\|:\mathbf{x}_{u},\mathbf{x}_{v}\in
X_{i}\right\\}\right)$,

  * 3.

$D_{3}(X_{i})=\mathrm{Mean}\left(\left\\{\|\mathbf{x}_{u}-\boldsymbol{\mu}_{i}\|:\mathbf{x}_{u}\in
X_{i}\right\\}\right)$.

There are 15 different combinations of the possible numerators and
denominators in our study, hence 15 different CVIs, which we will denote as
GDunn_dX_DY. In particular, GDunn_d1_D1 gives the original Dunn [18] index.

###  2.4 CVIs based on near-neighbour graphs

Let $\mathrm{NN}_{M}(i)=\\{{j_{1}},\dots,{j_{M}}\\}$ denote the set of the
$i$-th point’s $M$ nearest neighbours,
$0<\|\mathbf{x}_{i}-\mathbf{x}_{j_{1}}\|<\dots<\|\mathbf{x}_{i}-\mathbf{x}_{j_{M}}\|$
(assuming there are no tied distances, otherwise, some small random noise can
be added).

##### 21–50) DuNN_OWAs_OWAc

Note that the original Dunn index (denoted GDunn_d1_D1 above) can be viewed
as:

$\mathrm{Dunn}(C)=\frac{\mathrm{Min}\left(\left\\{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|:C(i)\neq
C(j)\right\\}\right)}{\mathrm{Max}\left(\left\\{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|:C(i)=C(j)\right\\}\right)}.$
(9)

Here we propose the following generalisation of the above – a generalised
Dunn-type index based on the notion of the $M$-near-neighbour graph and
ordered weighted averaging [60] operators – convex combinations (weighted
sums) of ordered inputs. Namely:

$\mathrm{DuNN}(C)=\frac{\mathrm{OWA}_{s}\left(\left\\{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|:C(i)\neq
C(j),i\in\mathrm{NN}_{M}(j)\text{ or
}j\in\mathrm{NN}_{M}(i)\right\\}\right)}{\mathrm{OWA}_{c}\left(\left\\{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|:C(i)=C(j),i\in\mathrm{NN}_{M}(j)\text{
or }j\in\mathrm{NN}_{M}(i)\right\\}\right)}.$ (10)

As a measure of cluster separation we aggregate the ordered between-point
distances but only provided that they are part of the near-neighbour graph.
This will enable us to take the local point density into account and detect
well-separable clusters of even quite sophisticated shapes. In a similar
manner, cluster compactness will be based on the nearest neighbours as well.

Below we shall study pairs of $\mathrm{OWA}_{s}$ and $\mathrm{OWA}_{c}$ chosen
amongst:

  * 1.

$\mathrm{Min}$,

  * 2.

$\mathrm{Max}$,

  * 3.

$\mathrm{Mean}$,

  * 4.

$\mathrm{SMin}_{\delta}(q_{1},q_{2},\dots,q_{z})=\sum_{i=1}^{z}w_{i,z}q_{(i)}$,
with
$w_{i,z}=\frac{\psi(i;z,\delta)}{\sum_{j=z-3\delta+1}^{z}\psi(j;z,\delta)}$
for $i>z-3\delta$ and 0 otherwise (“smooth minimum”),

  * 5.

$\mathrm{SMax}_{\delta}(q_{1},q_{2},\dots,q_{z})=\sum_{i=1}^{z}w_{i,z}q_{(i)}$,
with $w_{i,z}=\frac{\psi(i;z,\delta)}{\sum_{j=1}^{3\delta}\psi(j;1,\delta)}$
for $i\leq 3\delta$ and 0 otherwise (“smooth maximum”),

where $q_{(1)}\geq q_{(2)}\geq\dots\geq q_{(z)}$ and $\psi(\cdot;\mu,\sigma)$
denotes the probability density function of the normal distribution with
expectation $\mu$ and standard deviation $\sigma$, see also [10].

For instance, DuNN_25_SMin:5_Max denotes a generalised Dunn index based on
each point’s 25 nearest neighbours. It uses $\mathrm{SMin}_{5}$ as a
separation measure (computed over a subset of $25n$ distances restricted to
the pair of points belonging to different clusters) and $\mathrm{Max}$ as a
measure of compactness (the remainder of the $25n$ distances comprised of
point pairs belonging to the same clusters). Moreover, we will study indices
like DuNN_25_Mean_Const, where the denominator is fixed at 1.

In the sequel we will consider $M=5$ and $M=25$. In order to keep the number
of cases within reasonable limits, we will restrict ourselves to 30 different
CVIs of this type (see Table 3 for a complete listing).

##### 51,52) WCNN_M

The within-cluster near-neighbours (WCNN) index is parametrised by $M\geq 1$.
It aims to reflect how many nearest neighbours of every point actually belong
to the very same cluster:

$\mathrm{WCNN}(C)=\frac{\left|C(i)=C(j):j\in\mathrm{NN}_{M}(i)\right|}{nM}.$
(11)

Ideally, $\mathrm{WCNN}(C)=1$. Hence, this is a measure of how well the
clusters are separated from each other.

Additionally, to prevent the formation of small clusters, we will assume
$\mathrm{WCNN}(C)=-\infty$ whenever there is a cluster of cardinality $\leq
M$. Similarly as above, we shall consider $M\in\\{5,25\\}$.

Table 1:  Benchmark datasets studied, see [24] and https://github.com/gagolews/clustering_benchmarks_v1 for their visual depictions; $l$ gives the number of reference partitions and $k$s denote their possible cardinalities. | dataset |  $k$s |  $l$ |  | dataset |  $k$s | $l$
---|---|---|---|---|---|---|---
1. | fcps/atom | 2 |  1 | 32. | sipu/spiral | 3 | 1
2. | fcps/chainlink | 2 |  1 | 33. | sipu/unbalance | 8 | 1
3. | fcps/engytime | 2 |  2 | 34. | uci/ecoli | 8 | 1
4. | fcps/hepta | 7 |  1 | 35. | uci/ionosphere | 2 | 1
5. | fcps/lsun | 3 |  1 | 36. | uci/sonar | 2 | 1
6. | fcps/target | 2, 6 |  2 | 37. | uci/statlog | 7 | 1
7. | fcps/tetra | 4 |  1 | 38. | uci/wdbc | 2 | 1
8. | fcps/twodiamonds | 2 |  1 | 39. | uci/wine | 3 | 1
9. | fcps/wingnut | 2 |  1 | 40. | uci/yeast | 10 | 1
10. | graves/dense | 2 |  1 | 41. | wut/circles | 4 | 1
11. | graves/fuzzyx | 2, 4, 5 |  6 | 42. | wut/cross | 4 | 1
12. | graves/line | 2 |  1 | 43. | wut/graph | 10 | 1
13. | graves/parabolic | 2, 4 |  2 | 44. | wut/isolation | 3 | 1
14. | graves/ring | 2 |  1 | 45. | wut/labirynth | 6 | 1
15. | graves/ring_noisy | 2 |  1 | 46. | wut/mk1 | 3 | 1
16. | graves/ring_outliers | 2, 5 |  2 | 47. | wut/mk2 | 2 | 1
17. | graves/zigzag | 3, 5 |  2 | 48. | wut/mk3 | 3 | 1
18. | graves/zigzag_noisy | 3, 5 |  2 | 49. | wut/mk4 | 3 | 1
19. | graves/zigzag_outliers | 3, 5 |  2 | 50. | wut/olympic | 5 | 1
20. | other/chameleon_t4_8k | 6 |  1 | 51. | wut/smile | 4, 6 | 2
21. | other/chameleon_t5_8k | 6 |  1 | 52. | wut/stripes | 2 | 1
22. | other/hdbscan | 6 |  1 | 53. | wut/trajectories | 4 | 1
23. | other/iris | 3 |  1 | 54. | wut/trapped_lovers | 3 | 1
24. | other/iris5 | 3 |  1 | 55. | wut/twosplashes | 2 | 1
25. | other/square | 2 |  1 | 56. | wut/windows | 5 | 1
26. | sipu/aggregation | 7 |  1 | 57. | wut/x1 | 3 | 1
27. | sipu/compound | 4, 5, 6 |  5 | 58. | wut/x2 | 3 | 1
28. | sipu/flame | 2 |  2 | 59. | wut/x3 | 4 | 1
29. | sipu/jain | 2 |  1 | 60. | wut/z1 | 3 | 1
30. | sipu/pathbased | 3, 4 |  2 | 61. | wut/z2 | 5 | 1
31. | sipu/r15 | 8, 9, 15 |  3 | 62. | wut/z3 | 4 | 1

##  3 Method

###  3.1 What is a valid cluster validity index?

As we have proclaimed in the introduction, our key assumption in this paper is
that a meaningful cluster validity measure $I$ should be high whenever it is
asked to assess the quality of one of the reference partitions, and lower if
it is applied on other clusterings. In other words, useful CVIs should
encourage the results that agree with expert knowledge.

In order to be able to answer our main research question, i.e., which cluster
validity measures are valid, we need the following components:

  * 1.

benchmark data sets for evaluating the methods (Section 3.2),

  * 2.

a procedure for finding the partition that maximises a given CVI on each
dataset (Section 3.3),

  * 3.

a measure for quantifying the degree of agreement between what a CVI thinks is
a good partition vs what experts have to say on this matter (Section 3.4).

###  3.2 Benchmark Datasets

We shall use an extensive battery of clustering benchmarks [24]111 Available
at https://github.com/gagolews/clustering_benchmarks_v1. , which not only
combines data that have already been used in a number of studies [55, 17, 28,
21, 20, 32], but also features new test sets.

Most importantly, each benchmark dataset comes with a set of $l\geq 1$
reference labels that were assigned by experts. The case $l>1$ reflects the
situation where there might be multiple valid/plausible/useful partitions
(compare, e.g., [11]); we are dealing with an unsupervised learning problem
after all.

The original benchmark battery consists of 79 data instances, however 16
datasets are accompanied by labels that yield $n(k-1)>50{,}000$; they were
omitted for their computation would be too lengthy (namely: mnist/digits,
mnist/fashion, other/chameleon_t7_10k, other/chameleon_t8_8k, sipu/a1,
sipu/a2, sipu/a3, sipu/birch1, sipu/birch2, sipu/d31, sipu/s1, sipu/s2,
sipu/s3, sipu/s4, sipu/worms_2, sipu/worms_64). Also uci/glass has been
removed as one of its 25-near-neighbour graph’s connected components was too
small for the NN-based methods to succeed. This leaves us with 62 datasets in
total, see Table 1.

Further, all columns of 0 variance were removed and a tiny amount of noise
(Gaussian with $\mu=0$ and $\sigma$ equal to $10^{-6}$ of each column’s sample
standard deviation) was added so as to assure the uniqueness of the clustering
results.

###  3.3 Finding optimal partitions (w.r.t. a given CVI)

From now on we assume that the reader is familiar with the basics of the
language of mathematical programming, see, e.g., [49, 37] for a comprehensive
overview.

For a predefined $\mathbf{X}$ and $k$, let us fix a cluster validity measure
$I:\mathcal{C}_{k}\to\mathbb{R}$. Without loss in generality, we assume that
the higher the $I$, the more useful the partition. This is because we can
always take $I:=-I$, as we have done with the Ball–Hall and Davies–Bouldin
indices above.

For a given $k$-partition $C$, let $\mathrm{NEIGHBOURS}(C)$ denote the set of
all surjections like
$C^{\prime}:[1:n]\stackrel{{\scriptstyle\text{onto}}}{{\to}}[1:k]$ with
$C^{\prime}(i)\neq C(i)$ for some $i$ and $C^{\prime}(j)=C(j)$ for all $j\neq
i$. In other words, it is the set of all $k$-partitions that can be obtained
from $C$ by relocating a single point to some other cluster.

We are interested in finding a partition which is a solution to the
optimisation problem:

$\operatorname*{maximise}_{C\in\mathcal{C}_{k}}I(C),$ (12)

i.e., $C^{*}\in\mathcal{C}_{k}$ such that $I(C^{*})\geq I(C)$ for all
$C\in\mathcal{C}_{k}$.

######  Remark 1.

The solution to (12) is not unique; clusterings are defined up to a
permutation of the cluster numbers (IDs). For example, a 2-partition of a
4-ary set encoded like $(C(1),C(2),C(3),C(4))=(1,1,2,1)$ is semantically
equivalent to $(2,2,1,2)$. Moreover, it might happen that a dataset exhibits a
number of equally good splits. This is exactly the case when we apply WCNN_$M$
on datasets whose $M$-near-neighbour graphs are disconnected and the number of
connected components is greater than $k$. For instance, assuming
$Y_{1},Y_{2},Y_{3}$ are disconnected, in this setting the 2-partition
$\\{Y_{1}\cup Y_{2},Y_{3}\\}$ is as good as $\\{Y_{1},Y_{2}\cup Y_{3}\\}$.

In general, the combinatorial optimisation problem (12) is extremely difficult
to solve in practice. Enumerating all the possible solutions is virtually
impossible as the number of possible partitions is equal to the Stirling
number of the second kind, $S(n,k)=\frac{1}{k!}\sum_{i=0}^{k}(-1)^{i}{k\choose
i}(k-i)^{n}$ which is $O(k^{n})$, and note that in our case $2\leq k\ll n$.

In this paper, however, we shall make reasonable efforts towards finding the
maximum of the objective (12). In essence, for each dataset we will generate
dozens of “interesting” partitions (using existing state-of-the art clustering
algorithms and evolutionary-based heuristics, see Section 3.5) each of which
we shall then try to improve with an expansive variant of the steepest ascent
hill climbing (with tabu [27] search-like memoisation) that itself guarantees
to land in a local maximum of the objective (12).

The maximum of the objective (12) will be sought by means of the following
variant of the hill climbing scheme.

######  Algorithm 1.

With $\\{C_{1},\dots,C_{m}\\}$ let us denote the set of initial candidate
solutions (see Section 3.5 for more details) ordered in such a way that
$I(C_{1})\geq\dots\geq I(C_{m})$. For brevity of notation, we assume that
$I(\emptyset)=-\infty$.

  1. In:

$I:$, $C_{1},C_{2},\dots,C_{m}$, $P\in\mathbb{N}$;

  2. 1.

$\mathcal{T}=\emptyset$; (a “tabu” list)

  3. 2.

$C^{*}=C_{1}$; (best solution so far)

  4. 3.

for $C=C_{1},C_{2},\dots,C_{m}$ do: ($I(C_{1})\geq\dots\geq I(C_{m})$)

    1. 3.1.

$p=1$;

    2. 3.2.

$C^{+}=\emptyset$;

    3. 3.3.

for each $C^{\prime}\in\mathrm{NEIGHBOURS}(C)$ do:

      1. 3.3.1.

if $C^{\prime}\not\in\mathcal{T}$ and $I(C^{\prime})>I(C^{+})$, then
$C^{+}=C^{\prime}$;

    4. 3.4.

if $C^{+}=\emptyset$ then continue to step 3; (cannot improve further)

    5. 3.5.

$\mathcal{T}=\mathcal{T}\cup\\{C^{+}\\}$; (never visit $C^{+}$ again)

    6. 3.6.

$C=C^{+}$;

    7. 3.7.

if $I(C)>I(C^{*})$, then $C^{*}=C$, else $p=p+1$;

    8. 3.8.

if $p\leq P$, then go to step 3.2; (try to improve current $C^{+}$ next)

  5. 4.

return $C^{*}$;

It is easily seen that the algorithm guarantees that the solution returned
cannot be further improved by relocating an individual point to a different
cluster. Hence, the return value is definitely a local maximum, however there
is of course no guarantee that the identified optimum is global. As we argue
below, though, it will turn out sufficient for our purposes.

We shall set the upper bound for the number of iterations without improvement,
$P$, to $250$ (we have rarely seen any improvements beyond $100$, though).
This allows for the procedure to explore the area around the candidate
solutions quite broadly. Note that the $\mathcal{T}$ set, which guarantees
that no partition is considered twice, is shared across all the iterations so
that the visited subspace is even broader. Hence, the search is more
comprehensive than if we had restarted the whole procedure independently for
each $C_{1},\dots,C_{m}$ and then chose the best amongst the identified local
maxima.

######  Remark 2.

Note that the number of points in $\mathrm{NEIGHBOURS}(C)$ is $O(n(k-1))$.
Overall, the procedure for certain CVIs can be sped up by computing $I(C^{+})$
incrementally based on $I(C^{\prime})$ and the knowledge of which point is
being relocated to which cluster. For instance, the Silhouette index only
requires an $O(nk)$ update instead of a full recompute worth of $O(n^{2})$
time. The CVIs we have considered gave the time complexity of steps 3.1–3.8
most often lying between $O(n^{2}k^{2})$ and $O(n^{3}k^{2})$ (with constants
$d$ and $M$ having some obvious influence as well). The typical size of
$\mathcal{T}$ at the end of the algorithm’s run (i.e., the number of
executions of step 3.5) when started from $m=5$ random points was 1000–2000.

###  3.4 Measuring similarity to reference partitions

For a given benchmark dataset $\mathbf{X}$, let
$C_{1}^{\$},C_{2}^{\$},\dots,C_{l}^{\$}$ be the reference partitions and
$k_{1},\dots,k_{l}$ be their respective cardinalities.

Note that any clustering method $\mathfrak{c}$ (for example, the maximiser of
the Caliński–Harabasz index or the Ward linkage) can be thought of as a
function that takes $\mathbf{X}$ and $k_{i}$ on input and yields a
$k_{i}$-partition of $\mathbf{X}$ on output, i.e.,
$\mathfrak{c}(\mathbf{X},k_{i})\in\mathcal{C}_{k_{i}}$.

We will use the adjusted Rand index (ARI) [36, 53] to measure the similarity
between $\mathfrak{c}(\mathbf{X},k_{i})$ and $C_{i}^{\$}$. Recall that two
equivalent partitions yield ARI equal to 1. Moreover, two “independent” (see
[26] for discussion) clusterings have the expected ARI of 0. Negative ARIs
will be replaced with 0 for better interpretability of the results.

In order to quantify the quality of the method $\mathfrak{c}$ on $\mathbf{X}$,
we will evaluate its outputs against all the available reference labellings
and choose the highest ARI in result:

$Q_{\mathbf{X}}(\mathfrak{c})=\max\left\\{\mathrm{ARI}\left(\mathfrak{c}(\mathbf{X},k_{1}),C_{1}^{\$}\right),\dots,\mathrm{ARI}\left(\mathfrak{c}(\mathbf{X},k_{l}),C_{l}^{\$}\right)\right\\}.$
(13)

This is to account for the fact that there might be many equally valid
partitions and the (unsupervised) method $\mathfrak{c}$ should be rewarded if
it identifies one of them (does not matter which one).

It is worth noting that only 13 datasets have $l>1$, 11 of which come with
reference labellings that do not have identical cardinalities. Also, some
reference partitions include noise points – these were excluded during the
computations of the ARIs (after the output of $\mathfrak{c}$ was determined,
as none of the clustering methods studied features a noise point detector).
Overall, this validation methodology conforms with [24].

###  3.5 Candidate solutions

To generate the list of candidate (initial) solutions used in Algorithm 1, we
will apply many different clustering algorithms on each dataset, including:

  * 1.

the most popular hierarchical clustering methods (single, average, Ward,
centroid, complete linkage),

  * 2.

Genie [23] (with different thresholds),

  * 3.

information-theoretic algorithms (ITM [46] as well as IcA and GIc [9]),

  * 4.

other methods in the well-established sklearn [50] package for Python:
$k$-means, Gaussian mixtures, spectral clustering with different kernels,
Birch (with a range of parameter values).

This gives 87 different combinations of algorithms and their setups. In
Section 4 we provide some technical details about their implementations. Note
that 12 of them will constitute the baseline in our empirical study below.

Also, we shall utilise the following heuristic solvers:

  * 1.

particle swarm optimisation (via R package pso [4]),

  * 2.

“global” optimisation by differential evolution [51] (DEoptim [45] in R).

They pinpoint local maxima based on 3–5 restarts from different initial
candidate solutions. Both of them search over the continuous space
$\mathbb{R}^{(Vk)\times d}$ in such a way that the clusters are represented by
means of $Vk$ vantage points. $V$ vantage points represent one cluster
(empirically, we have determined $V=5$ be a good compromise between quality
and speed). In every iteration, each point is assigned to its closest vantage
point, and, as a consequence, to a cluster which is represented by this
vantage point. This approach allows to determine clusters of more
sophisticated shapes than when simply $V=1$ is utilised (as with $V>1$ we
consider different unions of cells in a Voronoi diagram).

Additionally, to broaden the search space even further, we will pick 5
partitions completely at random, i.e., each $C\in\mathcal{C}_{k}$ being such
that $C(1),\dots,\allowbreak C(n)$ being independent random variables from the
discrete uniform distribution on $[1:k]$. Nevertheless, starting from a random
partition never turned out better than an assisted initialisation based on one
of the aforementioned candidate solutions.

Most importantly, as each dataset comes with a set of reference labels given
by experts, these shall be considered as well.

Overall, for each dataset, we have obtained $m\simeq 100$ different
clusterings222All results are available at
https://github.com/gagolews/clustering_results_v1., however, quite often there
were duplicated entries, hence the effective $m$ was in the range 30–50.

An ideal index, if it existed, would be 100% concordant with expert labels.
That is, it would be impossible for the hill climbing method to improve them
any further. Note that when we maximise $I$, the reference partitions are
always amongst the initial candidate solutions which are fed to Algorithm 1.
Therefore, our procedure guarantees that all the good combinations of indices
and datasets must be identified. If the hill climbing method converges to a
different solution, it means that $I$ promotes some points that are less
compatible with experts’ opinion.

Let us stress that neither $I$ itself, nor the procedure for maximising it, is
“aware” of the existence of any external labels: only $\mathbf{X}$ and $k$ are
input to $\mathfrak{c}$, not $C_{i}^{\$}$; this is still an unsupervised
learning method. Algorithm 1 converges where it converges; the starting points
are plentiful and there is a great variety of them. Reference partitions are
only used at the final evaluation stage.

##  4 Experiments

###  4.1 Implementation

Experiments, data analysis, and visualisation tasks were performed using
Python 3.8.6 (PyPI packages: numpy 1.19.0, scipy 1.5.1, pandas 1.0.3,
matplotlib 3.3.3, seaborn 0.11.1) and R 4.0.3 (CRAN packages: DEoptim 2.2-5
[45], pso 1.0.3 [4]).

The correctness of our C++ implementations333Available at
https://github.com/gagolews/optim_cvi. of the cluster validity indices was
verified against R packages clusterCrit 1.2.8 and clusterSim 0.49-2 (wherever
applicable). Our library turned out significantly faster than the two
reference ones. Moreover, we allowed for the computing of the indices
incrementally (as mentioned above), which was particularly beneficial in terms
of the run-time of Algorithm 1.

Overall, the computations took ca. 3 months of computing time with the use of
2 computer clusters (within the allocated resource limits we have been granted
by ICM UW/PL-Grid and the School of IT at Deakin University).

###  4.2 Which index best agrees with expert knowledge?

Let us proceed with the evaluation of the agreement between the cluster
validity indices and expert knowledge.

####  4.2.1 GDunn_dX_DY

We first focus on the 15 generalised Dunn indices [6]. To recall, GDunn
indices are defined as the ratio of cluster separation ($d$) and compactness
($D$).

Figure 3 gives the box-and-whisker plots for the Adjusted Rand indices across
the benchmark datasets studied. We clearly see that $d_{1}$, i.e., the
pairwise minimum distance (used in the original Dunn index) outperforms the
other measures. In this scenario, the Wilcoxon signed-rank test does not find
the choice of $D$ significant ($\alpha=0.05$).

Also, the measures based on $d_{5}$ are significantly worse than all other
ones. Perhaps it would be better if $d_{5}$ was defined as the average squared
point-centroid distance, not just average raw distance; recall that a centroid
is the point that minimises exactly the square of the Euclidean metric.

Figure 3:  ARI: Generalised Dunn indices. We see that the choice of the
cluster compactness measure $D$ is rather negligible. On the other hand,
separation measure $d_{1}=\mathrm{Min}$ performs best whilst $d_{5}$ (averaged
distance to cluster centres) is subpar.

####  4.2.2 DuNN

Figure 4 shows the empirical distribution of ARIs in the case of the near-
neighbour versions of the Dunn index.

For a fixed separation measure $\mathrm{OWA}_{s}$, $\mathrm{OWA}_{c}$ equal to
$\mathrm{Min}$ and $\mathrm{Const}$ is never significantly worse (one-sided
Wilcoxon test, $\alpha=0.05$) than $\mathrm{Max}$ and $\mathrm{Mean}$.
Moreover, there is no significant difference between $\mathrm{Min}$ and
$\mathrm{Const}$, therefore, applying Ockham’s razor, we conclude that the
cluster compactness could be omitted whatsoever (at least as far as our
selection of aggregation functions is concerned).

Setting $\mathrm{OWA}_{c}$ at $\mathrm{Const}$, interestingly,
DuNN_25_SMin:5_Const significantly outperforms all the variants except
DuNN_5_Mean_Const.

Moreover, DuNN_25_Min_Const is better than DuNN_5_Min_Const. Also note that
the behaviour of Max or its smoothened version is particularly poor.

Figure 4:  ARI: Near-neighbour-based DuNN indices. Disabling the use of a
compactness measure whatsoever (*_Const) might be preferred. Also, it is
better to have the closest pairs of points from different clusters as far away
from each other as possible.

####  4.2.3 Other methods

As a base line, the above and remaining CVIs will be compared against the
outputs of 12 clustering algorithms:

  1. 1–5)

Average, Centroid, Complete, Ward, Single – classical agglomerative
hierarchical clustering algorithms;

  2. 6–9)

Genie_G0.1, Genie_G0.3, Genie_G0.5, Genie_G0.7 – the robust hierarchical
clustering algorithm Genie that we have proposed in [23], with different
thresholds for the Gini index of the inequity in cluster sizes;

  3. 10)

ITM – greedy divisive minimiser of an information theoretic criterion over
minimum spanning trees [46];

  4. 11)

GaussMix – expectation-maximisation (EM) for Gaussian mixtures with 100
restarts and each cluster having its own covariance matrix;

  5. 12)

KMeans – Lloyd-like $k$-means algorithm with 10 restarts (note that this is a
heuristic to optimise the Caliński–Harabasz index/within-cluster sum of
squares).

Their implementations are included in Python packages (available via PyPI; see
their respective API documentation for more details on algorithms and default
values of their parameters in place) fastcluster 1.1.26 (Average, Centroid,
Complete, Ward; [47]), genieclust 0.9.4 (Genie_G0.x, Single; [22]), sklearn
0.23.1 (GaussMix, KMeans; [50]). Moreover, the implementation of ITM [46] is
available from GitHub444 See https://github.com/amueller/information-
theoretic-mst; git commit 178fd43..

Tables 2 and 3 give the basic summary statistics on the empirical distribution
of the ARIs across the 62 benchmark datasets and all the methods studied.
Moreover, Figure 5 displays the boxplots.

We observe what follows:

  * 1.

Our Genie algorithm [22] outperforms other methods. Note that it is
significantly faster than most other algorithms as it is based on the minimum
spanning tree of the pairwise distance graph.

  * 2.

The lesser-known ITM method [46] performs relatively well.

  * 3.

The near-neighbour-based DuNN indices that we have proposed in this paper are
much better than GDunn.

  * 4.

The difference between DuNN_25_SMin:5_Const and WCNN_25 is insignificant.

  * 5.

Overall, the algorithms based on the near-neighbourhood (minimum spanning
trees can be considered a variant thereof, in some sense), seem to be much
more valid for clustering tasks.

  * 6.

Next in line are Gaussian mixtures that detect clusters of specific
(spherical) shapes.

  * 7.

Of course, $k$-means give similar results to the Caliński–Harabasz optimiser
as the former is a heuristic to optimise the latter as the objective. Also
Ward is a greedy agglomerative maximiser of the same objective.

  * 8.

Most other cluster validity measures seem to promote the clusterings that are
not concordant to expert knowledge which calls their relevance into question.

  * 9.

SilhouetteW and BallHall – the weighted-by-cluster-cardinality versions of
Silhouette and CalińskiHarabasz, respectively, perform worse than their
unadjusted counterparts.

  * 10.

The poor performance of some methods may be partially explained by the
inequality of the cluster sizes they output – some of them are prone to
generating few very large clusters and a number of very small ones (perhaps
even being singleton objects). This includes Single (median Gini index of the
cluster sizes=0.85), DaviesBouldin (0.95), and SilhouetteW (0.98). On the
other hand, Dunn_25_Max_Mean (median Gini index of the cluster sizes=0.09),
Dunn_5_Max_Mean (0.09), CalińskiHarabasz (0.13), KMeans (0.15), ITM (0.17),
and Genie_G0.1 (0.17) produced the least skewed partition sizes. However, let
us note that this is not necessarily an accurate predictor of the clustering
quality (see also [22] for discussion).

  * 11.

Note that some datasets are inherently hard to cluster (the outputs of no
algorithm matches the reference partition well); these include uci/sonar (max
ARI=0.036), uci/yeast (max ARI=0.181), and uci/ionosphere (max ARI=0.401).

###  4.3 Clustering of clustering algorithms

Let us perform an interesting exercise where we determine a grouping of the
clustering methods by means of the overall similarity of the results they
generate on all the benchmark datasets. This way, we will know which methods
(and CVIs) are “semantically” similar to each other.

We have computed the AR indices between all pairs of label vectors generated
by all the methods (this time, reference/expert labels were not used). We used
the mean, the median, or the 3rd quartile of $(1.0-\textrm{ARI})$ to obtain a
single number that summarises the “distance” between the algorithms.

We have applied the agglomerative hierarchical clustering algorithm with
complete linkage so that the resulting dendrograms, which are depicted in
Figure 6, are more interpretable (this will give us the maximal aggregated
dissimilarities between all methods in a cluster). Note that we will be only
interested in groups of algorithms that have small pairwise distances.

The majority vote of the results obtained by means of all the three
dissimilarity measures gives the following “consensus” clusters, where the
algorithms have quite high overall degree of similarity:

  * 1.

CalińskiHarabasz, KMeans,

  * 2.

DuNN_25_Mean_Const, DuNN_25_Mean_Min,

  * 3.

DuNN_5_Min_Const, DuNN_5_Min_Min,

  * 4.

DuNN_25_Min_Const, DuNN_25_Min_Min, DuNN_25_SMin:5_Const,

  * 5.

DuNN_5_Mean_Const, DuNN_5_Mean_Min,

  * 6.

GDunn_d1_D1, GDunn_d1_D2, GDunn_d1_D3, Single.

Let us note that:

  * 1.

The above indicates again that in our task, the denominator (compactness
measure) in all the generalisations of the Dunn index has no significant
impact on the results. However, when CVIs are applied for the purpose of
selecting the optimal number of clusters, the conclusions could be much
different.

  * 2.

GDunn and DuNN are not so similar, despite they both generalise the same
index, Dunn. On the other hand, the latter (to recall, it is based on
$d_{1}=\mathrm{Min}$) is similar to Single linkage (which is determined
through a greedy (agglomerative) consumption of the nearest pairs of points;
it can be computed based on a minimum spanning tree).

  * 3.

Ward is quite similar to CalińskiHarabasz, and KMeans, which was to be
expected as they tend to optimise the same objective function.

##  5 Conclusion and Future Work

We have studied whether cluster validity indices really promote partitions
that reflect what experts judge as meaningful. While some measures could still
be considered relevant in the task of selecting the right number of clusters,
it is better not to treat them as objective functions for identifying good
partitions. This is particularly the case with the Davies–Bouldin, Silhouette,
Ball–Hall, and the no-near-neighbour-based versions of the generalised Dunn
index.

In the future, we will verify the usability of our new near-neighbour-based
generalisations of the Dunn index in the problem of choosing the meaningful
number of clusters. Their advantage is that they take into account the
locality of the input points as well as the relative density of the points’
distribution. Certainly, more combinations of OWA operators should be studied.

## Acknowledgements

This research was supported by the Australian Research Council Discovery
Project ARC DP210100227 as well as the PL-Grid Infrastructure.

## Conflict of interest

The authors declare that they have no conflict of interest to disclose.

Figure 5:  ARI: All methods (ordered by the average ARI). The Genie algorithm outperforms other clustering approaches. Near-neighbour-based cluster validity measures reflect expert knowledge quite well. Table 2:  ARI: Basic summary statistics; part I.  Method | Mean | St.Dev. | Q1 | Median | Q3
---|---|---|---|---|---
Average | 0.47 | 0.38 | 0.10 | 0.44 | 0.88
Centroid | 0.44 | 0.39 | 0.05 | 0.41 | 0.86
Complete | 0.43 | 0.33 | 0.17 | 0.36 | 0.72
GaussMix | 0.63 | 0.38 | 0.36 | 0.76 | 0.99
Genie_G0.1 | 0.73 | 0.32 | 0.50 | 0.86 | 1.00
Genie_G0.3 | 0.78 | 0.29 | 0.57 | 0.98 | 1.00
Genie_G0.5 | 0.77 | 0.32 | 0.59 | 0.97 | 1.00
Genie_G0.7 | 0.68 | 0.38 | 0.29 | 0.96 | 1.00
ITM | 0.68 | 0.28 | 0.53 | 0.67 | 1.00
KMeans | 0.50 | 0.35 | 0.18 | 0.48 | 0.84
Single | 0.49 | 0.47 | 0.00 | 0.52 | 1.00
Ward | 0.50 | 0.35 | 0.19 | 0.50 | 0.83
BallHall | 0.40 | 0.40 | 0.01 | 0.29 | 0.85
CalińskiHarabasz | 0.49 | 0.35 | 0.18 | 0.48 | 0.84
DaviesBouldin | 0.31 | 0.37 | 0.00 | 0.12 | 0.51
Silhouette | 0.44 | 0.37 | 0.11 | 0.43 | 0.72
SilhouetteW | 0.15 | 0.29 | 0.00 | 0.00 | 0.10
WCNN_25 | 0.71 | 0.38 | 0.43 | 1.00 | 1.00
WCNN_5 | 0.67 | 0.40 | 0.32 | 0.99 | 1.00
GDunn_d1_D1 | 0.58 | 0.44 | 0.01 | 0.74 | 1.00
GDunn_d1_D2 | 0.59 | 0.45 | 0.01 | 0.77 | 1.00
GDunn_d1_D3 | 0.59 | 0.45 | 0.01 | 0.77 | 1.00
GDunn_d2_D1 | 0.40 | 0.34 | 0.07 | 0.37 | 0.63
GDunn_d2_D2 | 0.46 | 0.32 | 0.14 | 0.48 | 0.70
GDunn_d2_D3 | 0.46 | 0.33 | 0.13 | 0.46 | 0.72
GDunn_d3_D1 | 0.40 | 0.35 | 0.02 | 0.38 | 0.63
GDunn_d3_D2 | 0.43 | 0.36 | 0.08 | 0.44 | 0.79
GDunn_d3_D3 | 0.42 | 0.36 | 0.07 | 0.41 | 0.77
GDunn_d4_D1 | 0.40 | 0.36 | 0.01 | 0.38 | 0.63
GDunn_d4_D2 | 0.42 | 0.36 | 0.07 | 0.43 | 0.73
GDunn_d4_D3 | 0.41 | 0.36 | 0.03 | 0.39 | 0.72
GDunn_d5_D1 | 0.13 | 0.20 | 0.02 | 0.07 | 0.14
GDunn_d5_D2 | 0.16 | 0.19 | 0.02 | 0.09 | 0.23
GDunn_d5_D3 | 0.12 | 0.28 | 0.00 | 0.00 | 0.01
Table 3:  ARI: Basic summary statistics; part II.  Method | Mean | St.Dev. | Q1 | Median | Q3
---|---|---|---|---|---
DuNN_5_Max_Const | 0.47 | 0.44 | 0.09 | 0.26 | 1.00
DuNN_5_Mean_Const | 0.70 | 0.38 | 0.44 | 0.99 | 1.00
DuNN_5_Min_Const | 0.69 | 0.39 | 0.34 | 0.97 | 1.00
DuNN_25_Max_Const | 0.37 | 0.39 | 0.05 | 0.20 | 0.60
DuNN_25_Mean_Const | 0.67 | 0.38 | 0.31 | 0.92 | 1.00
DuNN_25_Min_Const | 0.72 | 0.36 | 0.42 | 0.97 | 1.00
DuNN_25_SMax:5_Const | 0.35 | 0.39 | 0.06 | 0.20 | 0.55
DuNN_25_SMin:5_Const | 0.72 | 0.37 | 0.42 | 0.99 | 1.00
DuNN_5_Max_Min | 0.48 | 0.43 | 0.09 | 0.32 | 1.00
DuNN_5_Mean_Min | 0.69 | 0.38 | 0.39 | 0.99 | 1.00
DuNN_5_Min_Min | 0.69 | 0.39 | 0.34 | 0.97 | 1.00
DuNN_25_Max_Min | 0.37 | 0.39 | 0.05 | 0.20 | 0.65
DuNN_25_Mean_Min | 0.67 | 0.38 | 0.31 | 0.92 | 1.00
DuNN_25_Min_Min | 0.72 | 0.35 | 0.42 | 0.97 | 1.00
DuNN_25_SMax:5_SMin:5 | 0.36 | 0.38 | 0.05 | 0.24 | 0.39
DuNN_25_SMax:5_Min | 0.34 | 0.39 | 0.01 | 0.20 | 0.48
DuNN_5_Max_Max | 0.40 | 0.47 | 0.00 | 0.06 | 1.00
DuNN_5_Mean_Max | 0.57 | 0.44 | 0.06 | 0.61 | 1.00
DuNN_5_Min_Max | 0.68 | 0.40 | 0.31 | 0.97 | 1.00
DuNN_25_Max_Max | 0.38 | 0.41 | 0.00 | 0.19 | 0.89
DuNN_25_Mean_Max | 0.64 | 0.39 | 0.26 | 0.82 | 1.00
DuNN_25_Min_Max | 0.71 | 0.36 | 0.47 | 0.97 | 1.00
DuNN_25_SMin:5_Max | 0.70 | 0.37 | 0.41 | 0.99 | 1.00
DuNN_25_SMin:5_SMax:5 | 0.71 | 0.37 | 0.46 | 0.98 | 1.00
DuNN_5_Max_Mean | 0.40 | 0.47 | 0.00 | 0.04 | 1.00
DuNN_5_Mean_Mean | 0.60 | 0.44 | 0.04 | 0.93 | 1.00
DuNN_5_Min_Mean | 0.66 | 0.41 | 0.25 | 0.95 | 1.00
DuNN_25_Max_Mean | 0.27 | 0.42 | 0.00 | 0.02 | 0.30
DuNN_25_Mean_Mean | 0.57 | 0.43 | 0.08 | 0.60 | 1.00
DuNN_25_Min_Mean | 0.69 | 0.38 | 0.37 | 0.92 | 1.00

Figure 6:  Complete linkage with Median(1.0-ARI), Q3(1.0-ARI), Mean(1.0-ARI),
respectively. Only clusters that are formed at low dissimilarity levels
(branches with connectors in the right parts of each figure) can be considered
meaningful.

## References

  * Aloise et al. [2009] D. Aloise, A. Deshpande, P. Hansen, P. Popat, NP-hardness of Euclidean sum-of-squares clustering, Machine Learning 75 (2009) 245–248.
  * Arbelaitz et al. [2013] O. Arbelaitz, I. Gurrutxaga, J. Muguerza, J.M. Pérez, I. Perona, An extensive comparative study of cluster validity indices, Pattern Recognition 46 (2013) 243–256.
  * Ball and Hall [1965] G. Ball, D. Hall, ISODATA: A novel method of data analysis and pattern classification, Technical Report AD699616, 1965\.
  * Bendtsen [2012] C. Bendtsen, pso: Particle Swarm Optimization, 2012. R package version 1.0.3; https://CRAN.R-project.org/package=pso.
  * Bezdek et al. [1984] J. Bezdek, R. Ehrlich, W. Full, FCM: The fuzzy c-means clustering algorithm, Computers & Geosciences 10 (1984) 191–203.
  * Bezdek and Pal [1998] J. Bezdek, N. Pal, Some new indexes of cluster validity, IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 28 (1998) 301–315.
  * Brun et al. [2007] M. Brun, C. Sima, J. Hua, J. Lowey, B. Carroll, E. Suh, E.R. Dougherty, Model-based evaluation of clustering validation measures, Pattern Recognition 40 (2007) 807–824.
  * Caliński and Harabasz [1974] T. Caliński, J. Harabasz, A dendrite method for cluster analysis, Communications in Statistics 3 (1974) 1–27.
  * Cena [2018] A. Cena, Adaptive hierarchical clustering algorithms based on data aggregation methods, Ph.D. thesis, Systems Research Institute, Polish Academy of Sciences, 2018. In Polish.
  * Cena and Gagolewski [2020] A. Cena, M. Gagolewski, Genie+OWA: Robustifying hierarchical clustering with OWA-based linkages, Information Sciences 520 (2020) 324–336.
  * Dasgupta and Ng [2009] S. Dasgupta, V. Ng, Single data, multiple clusterings, in: Proc. NIPS Workshop Clustering: Science or Art? Towards Principled Approaches, 2009\. http://clusteringtheory.org.
  * Davies and Bouldin [1979] D.L. Davies, D.W. Bouldin, A cluster separation measure, IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI–1 (1979) 224–227.
  * Dhal et al. [2021] K.G. Dhal, A. Das, S. Ray, J. Gálvez, Randomly attracted rough firefly algorithm for histogram based fuzzy image clustering, Knowledge-Based Systems 216 (2021) 106814.
  * Dimitriadou et al. [1999] E. Dimitriadou, S. Dolnicar, F. Leisch, A. Weingessel, More insight into clustering: Comparison of cluster algorithms and evaluation of indexes for determining the correct number of clusters, Methods of Psychological Research 4 (1999) 65–66.
  * Dinh et al. [2021] D.T. Dinh, V.N. Huynh, S. Sriboonchitta, Clustering mixed numerical and categorical data with missing values, Information Sciences 571 (2021) 418–442.
  * Du et al. [2021] M. Du, R. Wang, R. Ji, X. Wang, Y. Dong, ROBP a robust border-peeling clustering using Cauchy kernel, Information Sciences 571 (2021) 375–400.
  * Dua and Graff [2021] D. Dua, C. Graff, UCI Machine Learning Repository, 2021. http://archive.ics.uci.edu/ml.
  * Dunn [1974] J. Dunn, A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters, Journal of Cybernetics 3 (1974) 32–57.
  * Edwards and Cavalli-Sforza [1965] A.W.F. Edwards, L.L. Cavalli-Sforza, A method for cluster analysis, Biometrics 21 (1965) 362–375.
  * Fränti et al. [2016] P. Fränti, R. Mariescu-Istodor, C. Zhong, XNN graph, Lecture Notes in Computer Science 10029 (2016) 207–217.
  * Fränti and Sieranoja [2018] P. Fränti, S. Sieranoja, K-means properties on six clustering benchmark datasets, Applied Intelligence 48 (2018) 4743–4759.
  * Gagolewski [2021] M. Gagolewski, genieclust: Fast and robust hierarchical clustering, SoftwareX 15 (2021) 100722.
  * Gagolewski et al. [2016] M. Gagolewski, M. Bartoszuk, A. Cena, Genie: A new, fast, and outlier-resistant hierarchical clustering algorithm, Information Sciences 363 (2016) 8–23.
  * Gagolewski et al. [2020] M. Gagolewski, et al., Benchmark suite for clustering algorithms – version 1, 2020. https://github.com/gagolews/clustering_benchmarks_v1, doi:10.5281/zenodo.3815066.
  * Garey et al. [1982] M. Garey, D. Johnson, H. Witsenhausen, The complexity of the generalized Lloyd–Max problem, IEEE Transactions on Information Theory 28 (1982) 255–256.
  * Gates and Ahn [2017] A.J. Gates, Y.Y. Ahn, The impact of random models on clustering similarity, Journal of Machine Learning Research 18 (2017) 1–28.
  * Glover [1986] F. Glover, Future paths for integer programming and links to artificial intelligence, Computers & Operations Research 13 (1986) 533–549.
  * Graves and Pedrycz [2010] D. Graves, W. Pedrycz, Kernel-based fuzzy clustering: A comparative experimental study, Fuzzy Sets and Systems 161 (2010) 522–543.
  * Halkidi et al. [2001] M. Halkidi, Y. Batistakis, M. Vazirgiannis, On clustering validation techniques, Journal of Intelligent Information Systems (2001) 107–145.
  * Isimeto et al. [2017] R. Isimeto, C. Yinka-Banjo, C.O. Uwadia, D.C. Alienyi, An enhanced clustering analysis based on glowworm swarm optimization, in: 2017 IEEE 4th International Conference on Soft Computing Machine Intelligence (ISCMI), pp. 42–49.
  * Jamil and Yang [2013] M. Jamil, X.S. Yang, A literature survey of benchmark functions for global optimization problems, International Journal of Mathematical Modelling and Numerical Optimisation 4 (2013).
  * Karypis et al. [1999] G. Karypis, E. Han, V. Kumar, CHAMELEON: Hierarchical clustering using dynamic modeling, Computer 32 (1999) 68–75.
  * Kim and Ramakrishna [2005] M. Kim, R. Ramakrishna, New indices for cluster validity assessment, Pattern Recognition Letters 26 (2005) 2535–2363.
  * Kuo et al. [2021] R. Kuo, Y. Zheng, T.P.Q. Nguyen, Metaheuristic-based possibilistic fuzzy k-modes algorithms for categorical data clustering, Information Sciences 557 (2021) 1–15.
  * Lance and Williams [1967] G. Lance, W. Williams, A general theory of classification sorting strategies: 1. Hierarchical systems, Computer Journal (1967) 373–380.
  * Lawrence and Phipps [1985] H. Lawrence, A. Phipps, Comparing partitions, Journal of Classification 2 (1985) 193–218.
  * Lee [2011] J. Lee, A First Course in Combinatorial Optimisation, Cambridge University Press, 2011\.
  * Li et al. [2016] H. Li, S. Zhang, X. Ding, C. Zhang, P. Dale, Performance evaluation of cluster validity indices (cvis) on multi/hyperspectral remote sensing datasets, Remote Sensing 8 (2016).
  * Liang et al. [2020] S. Liang, D. Han, Y. Yang, Cluster validity index for irregular clustering results, Applied Soft Computing 95 (2020) 106583.
  * Liu et al. [2021] Y. Liu, Y. Jiang, T. Hou, F. Liu, A new robust fuzzy clustering validity index for imbalanced data sets, Information Sciences 547 (2021) 579–591.
  * Lloyd [1982] S. Lloyd, Least squares quantization in PCM, IEEE Transactions on Information Theory 28 (1957 (1982)) 128–137. Originally a 1957 Bell Telephone Laboratories Research Report; republished in 1982.
  * Maulik and Bandyopadhyay [2002] U. Maulik, S. Bandyopadhyay, Performance evaluation of some clustering algorithms and validity indices, IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (2002) 1650–1654.
  * Milligan and Cooper [1985] G.W. Milligan, M.C. Cooper, An examination of procedures for determining the number of clusters in a data set, Psychometrika 50 (1985) 159–179.
  * Mishra et al. [2021] G. Mishra, A.K. Kar, A.C. Mishra, S.K. Mohanty, M. Panda, SEND: A novel dissimilarity metric using ensemble properties of feature space for clustering numerical data, Information Sciences 574 (2021) 279–296.
  * Mullen et al. [2011] K. Mullen, D. Ardia, D. Gil, D. Windover, J. Cline, DEoptim: An R package for global optimization by differential evolution, Journal of Statistical Software 40 (2011) 1–26.
  * Müller et al. [2012] A. Müller, S. Nowozin, C. Lampert, Information theoretic clustering using minimum spanning trees, in: Proc. German Conference on Pattern Recognition, 2012. https://github.com/amueller/information-theoretic-mst.
  * Müllner [2013] D. Müllner, fastcluster: Fast hierarchical, agglomerative clustering routines for R and Python, Journal of Statistical Software 53 (2013) 1–18.
  * Nanda and Panda [2014] S.J. Nanda, G. Panda, A survey on nature inspired metaheuristic algorithms for partitional clustering, Swarm and Evolutionary Computation 16 (2014) 1–18.
  * Nocedal and Wright [2006] J. Nocedal, S.J. Wright, Numerical Optimization, Springer, 2006\.
  * Pedregosa et al. [2011] F. Pedregosa, et al., Scikit-learn: Machine learning in Python, Journal of Machine Learning Research 12 (2011) 2825–2830.
  * Price et al. [2006] K.V. Price, R.M. Storn, J.A. Lampinen, Differential Evolution – A Practical Approach to Global Optimization, Springer-Verlag, 2006.
  * Qaddoura et al. [2020] R. Qaddoura, H. Faris, I. Aljarah, An efficient evolutionary algorithm with a nearest neighbor search technique for clustering analysis, Ambient Intell Human Comput (2020).
  * Rezaei and Fränti [2016] M. Rezaei, P. Fränti, Set matching measures for external cluster validity, IEEE Transactions on Knowledge and Data Engineering 28 (2016) 2173–2186.
  * Rousseeuw [1987] P.J. Rousseeuw, Silhouettes: A graphical aid to the interpretation and validation of cluster analysis, Journal of Computational and Applied Mathematics 20 (1987) 53–65.
  * Ultsch [2005] A. Ultsch, Clustering with SOM: U*C, in: Workshop on Self-Organizing Maps, WSOM 2005, 2005, pp. 75–82.
  * Vij and Khandnor [2016] A. Vij, P. Khandnor, Validity of internal cluster indices, in: International Conference on Computational Systems for Sustainable Solutions, pp. 388–395.
  * Ward Jr. [1963] J.H. Ward Jr., Hierarchical grouping to optimize an objective function, Journal of the American Statistical Association 58 (1963) 236–244.
  * Xu et al. [2020] Q. Xu, Q. Zhang, J. Liu, B. Luo, Efficient synthetical clustering validity indexes for hierarchical clustering, Expert Systems with Applications 151 (2020) 113367\.
  * Xu et al. [2012] R. Xu, J. Xu, D.C. Wunsch, A comparison study of validity indices on swarm-intelligence-based clustering, IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 42 (2012) 1243–1256.
  * Yager [1988] R.R. Yager, On ordered weighted averaging aggregation operators in multicriteria decision making, IEEE Transactions on Systems, Man, and Cybernetics 18 (1988) 183–190.
  * Zhu et al. [2020] S. Zhu, L. Xu, E.D. Goodman, Evolutionary multi-objective automatic clustering enhanced with quality metrics and ensemble strategy, Knowledge-Based Systems 188 (2020) 105018.


# Energy-gap driven low-temperature magnetic and transport properties in
Cr${}_{1/3}M$S2 ($M$ = Nb or Ta)

T. J. Hicken  Current address: Department of Physics, Royal Holloway,
University of London, Egham, TW20 0EX, United Kingdom  Centre for Materials
Physics, Department of Physics, Durham University, Durham, DH1 3LE, United
Kingdom  Z. Hawkhead  Centre for Materials Physics, Department of Physics,
Durham University, Durham, DH1 3LE, United Kingdom  M. N. Wilson  Centre for
Materials Physics, Department of Physics, Durham University, Durham, DH1 3LE,
United Kingdom  B. M. Huddart  Centre for Materials Physics, Department of
Physics, Durham University, Durham, DH1 3LE, United Kingdom  A. E. Hall
Department of Physics, University of Warwick, Coventry, CV4 7AL, United
Kingdom  G. Balakrishnan  Department of Physics, University of Warwick,
Coventry, CV4 7AL, United Kingdom  C. Wang  Laboratory for Muon-Spin
Spectroscopy, Paul Scherrer Institut, Forschungsstrasse 111, 5232 Villigen
PSI, Switzerland  F. L. Pratt  ISIS Pulsed Neutron and Muon Facility, STFC
Rutherford Appleton Laboratory, Harwell Oxford, Didcot, OX11 OQX, United
Kingdom  S. J. Clark  Centre for Materials Physics, Department of Physics,
Durham University, Durham, DH1 3LE, United Kingdom  T. Lancaster  Centre for
Materials Physics, Department of Physics, Durham University, Durham, DH1 3LE,
United Kingdom

###### Abstract

The helimagnets Cr${}_{1/3}M$S2 ($M$ = Nb or Ta) have attracted renewed
attention due to the discovery of a chiral soliton lattice (CSL) stabilized in
an applied magnetic field, but reports of unusual low-temperature transport
and magnetic properties in this system lack a unifying explanation. Here we
present electronic structure calculations that demonstrate the materials are
half-metals. There is also a gap-like feature (width in range 40–100 meV) in
the density of states of one spin channel. This electronic structure explains
the low-temperature electronic and magnetic properties of Cr${}_{1/3}M$S2 ($M$
= Nb or Ta), with the gap-like feature particularly important for explaining
the magnetic behavior. Our magnetometry measurements confirm the existence of
this gap. Dynamic spin fluctuations driven by excitations across this gap are
seen over a wide range of frequencies (0.1 Hz to MHz) with AC susceptibility
and muon-spin relaxation ($\mu^{+}$SR) measurements. We show further how
effects due to the CSL in Cr1/3NbS2, as detected with $\mu^{+}$SR, dominate
over the gap-driven magnetism when the CSL is stabilized as the majority
phase.

Transition metal dichalcogenides ($MX_{2}$, $M=$ transition metal and $X=$
chalcogen) have long been the subject of much research effort [1], the reduced
dimensionality [2] of these semiconductors leads to a suite of exciting
electrical, optical and mechanical properties [3, 4, 5]. More recently,
Cr${}_{1/3}M$S2, $M$ = Nb or Ta (based on the same structure as
superconducting NbS2, but with Cr atoms intercalated between $M$ layers) have
attracted attention as rare examples of 2D materials with helical magnetic
ordering, thought to occur due to the competition between exchange, the
Dzyaloshinksii-Moriya interaction (DMI) and spin-orbit coupling (SOC) [6, 7].
Application of a magnetic field perpendicular to the $c$ axis causes the
helical ground state to transform to the chiral soliton lattice (CSL) [8, 9],
consisting of domain-wall-like 360$\degree$ magnetic kinks separated by
regions of ferromagnetic ordering, with a field along a different direction
tilting the CSL [10, 11]. These magnetic solitons, reduced-dimensional
analogues of skyrmions [12], are topological spin structures with many
potential applications [8]. The existence of the CSL in Cr1/3NbS2 is well
established [8], but has only recently been reported in Cr1/3TaS2 [9].

There have been multiple reports of unusual low-temperature transport and
magnetic properties in Cr1/3NbS2, including a sizable change in the thermal
conductivity and Seebeck coefficient around $T=40$ K [13], and a rapid
increase in the Hall coefficient below 50 K [14], all suggesting modifications
to the electronic behavior. Additionally, muon-spin rotation ($\mu^{+}$SR)
identifies an increase in oscillation frequency below 50 K [15], suggesting an
alteration to the local magnetic-field distribution. The above observations
have led to contrasting explanations, including a change in the dominant
transport mechanism due to spin-orbit-coupling-induced alterations of the
electronic structure [13], or an increased helical length at low $T$ due to a
decrease in the magnitude of the DMI [15]. Such explanations do not account
for all of the observations, and while an increase in the helical length below
90 K is observed in lamellae [16], this is not found in bulk [17]. In this
Letter we explain the low-temperature behavior of Cr${}_{1/3}M$S2 ($M$ = Nb or
Ta) using electronic structure calculations supported by AC susceptibility and
$\mu^{+}$SR, and demonstrate how the presence of an energy-gap-like feature in
the density of states can account for the observed phenomena over a wide range
of energy scales.

The electronic structure of Cr${}_{1/3}M$S2 ($M$ = Nb or Ta) was determined
using the castep implementation of spin-polarized density functional theory
(DFT) [18, 19, 20, 21]. Spectral calculations were performed to generate band
structures and densities of states (DOS) for Cr${}_{1/3}M$S2 (Fig. 1). In both
materials there is splitting of the spin-up and spin-down bands with
ferromagnetic moments on the Cr ions. The largest contribution to the DOS at
the Fermi level comes from the spin-up channel, with little contribution from
spin-down. Strongly peaked bands at the $\Gamma$-point in the spin-down
channel suggest delocalized electrons; more localized character in the spin-up
channel is indicated by flatter bands [21]. The lack of spin-down density
around the Fermi level, in addition to the integer spin per unit cell [22],
suggests both materials are half-metallic [23]. (Previously Cr1/3NbS2 has been
reported as a low carrier concentration metal or heavily doped semiconductor
[13].) Additionally, in both materials, there is a reduction to almost zero in
the spin-up channel of the DOS at the Fermi level with a full-width half-
maxima of approximately $\Delta E=80$–$90$ meV. This reduction has previously
been observed in Cr1/3NbS2 and termed a pseudogap [13], but the influence of
this feature was not investigated. Our calculations performed with the
inclusion of SOC have no resolvable impact on the pseudogap [21].

Figure 1: (a,b) DOS of $d$-orbital contributions from transition metals in
Cr1/3NbS2 and Cr1/3TaS2. (c,d) Regions around the Fermi level showing the
pseudogap feature in the spin-up channel. (Spin-up densities are positive and
spin-down densities are negative.)

To understand the low-temperature transport properties in Cr${}_{1/3}M$S2 we
calculate the Seebeck coefficient $S_{\mu\nu}$ from DFT [21] using the
boltztrap code [24]. Figure 2(a–b) shows
$S=\mathrm{Tr}\left(S_{\mu\nu}\right)/3$. This prediction of $S$ for Cr1/3NbS2
qualitatively matches the measurements in the low-temperature regime [13],
although around $T_{\text{c}}$ the computed behavior deviates from experiment
due to the onset of the magnetic transition, which is not captured in our
model. The good agreement between DFT and experiment verifies the calculated
electronic structure. The pseudogap found in DFT is slightly larger than the
measured gap (see below), suggesting there are small differences in the energy
scale around the Fermi energy. This results in the features in $S$ to occur at
higher $T$ than found experimentally. To better compare to the experimental
values, we can adjust the energy scale in the calculation to that of the
measured gap [$T\rightarrow T\Delta E_{\text{exp}}/\Delta E_{\text{DFT}}$,
Fig. 2(a)], which then provides a better agreement with the experimental peak
position. We find similar behavior of $S$ for Cr1/3TaS2 [Fig 2(b)], as
expected from the similar electronic structures. Our DFT calculations can also
reproduce the observed rapid increase in Hall coefficient on decreasing $T$ in
Cr1/3NbS2 below 50 K [14, 21].

Figure 2: Seebeck coefficient $S=\mathrm{Tr}\left(S_{\mu\nu}\right)/3$ for (a)
Cr1/3NbS2 and (b) Cr1/3TaS2. In (a) $S$ is compared to experiment [13], and is
also shown with temperatures scaled according to the gap size (see text). In
(b) the range of the possible values of the Seebeck coefficient due to
uncertainty in the Fermi energy is shown. (c–d) Temperature evolution of the
change in $m$ for Cr1/3NbS2 (c) and Cr1/3TaS2 (d). (Total change, and change
in Cr and $M$ $d$-orbitals is shown.) (e–f) Effect on the field at the muon-
site in (e) Cr1/3NbS2 and (f) Cr1/3TaS2. Field for Cr moment is shown with a
vertical offset. Inset: measured field at the muon-site [15] [same axes as
(e), dashed line is a guide to the eye].

The magnetic moment $m$ can be estimated from the difference in the DOS in the
spin-up and spin-down channels. We find, for both materials, $d$-orbitals
contribute most significantly to $m$, with around 95% of the total moment
coming from the Cr $d$-orbitals. We find good agreement between the saturation
moment of our calculations and that found experimentally. The $T$ dependence,
$m(T)$, caused by thermally-activated changes in band occupation (but
excluding critical behavior due to the phase transition at $T_{\text{c}}$ [13,
25]) can be estimated by integrating the DOS weighted by the Fermi-Dirac
distribution $f_{\mathrm{FD}}(E,T)$ at a given temperature,
$m(T)=\mu_{\text{B}}\int_{-\infty}^{\infty}\,\left[g^{\uparrow}(E)-g^{\downarrow}(E)\right]f_{\mathrm{FD}}(E,T)~{}\mathrm{d}E$,
where $g^{\uparrow(\downarrow)}(E)$ are the spin-up (spin-down) components of
the partial DOS. The change in $m(T)$ (compared to $T=0$) is shown in Fig.
2(c–d). As $T$ is increased, electrons start to populate excited states,
leading to a steady increase in $m(T)$, with the change being dominated by the
$M$ = Nb or Ta $d$-orbitals [Fig. 2(c–d)], despite their small overall
contribution to the total moment.

Although these changes in $m(T)$ are likely too small to be detected by a bulk
measurement, $\mu^{+}$SR, a local probe sensitive to the magnitude of the
internal field at the muon site, can detect such small changes. Despite the
lack of a magnetic phase transition at low $T$, on the basis of continuous
source zero-field (ZF) and transverse-field (TF) measurements, Ref. [15]
reports an increase in the local field at the muon site [inset of Fig. 2(e)].
We calculated the muon stopping site in these materials using DFT [21, 26, 27]
(finding different sites to those determined using nuclear dipole field
calculations [15]). The measured field can be compared to predictions using
our candidate stopping site to calculate the expected $T$-dependent dipole
field at the muon site [28] due to our predicted $m(T)$ [Fig. 2(c–d)] [21].
The field at the most probable muon site is highly sensitive to both the
amount of Cr intercalated in the material, and to the precise location of the
spin density with respect to the muon, neither of which are easy to accurately
estimate. Therefore, our calculation of the absolute magnitude of the field at
the muon site, as well as the size of the increase with temperature, does not
quantitatively match experiment. Some of the experimentally observed change
may be due to the hyperfine interaction, which typically follows the dipolar
field, or due to delocalization of spin; it is not possible to accurately
compute these contributions [21]. Despite these limitations, we predict the
same trend as found in Ref. [15] for Cr1/3NbS2 [Fig. 2(e)], with a sharp
increase in the field at the muon-site at low $T$. It is notable that the
moment of all elements must be considered to obtain this effect; considering
Cr alone does not reproduce the trend.

To confirm the influence of the pseudogap on the magnetism, we performed
magnetization and AC susceptibility measurements on polycrystalline samples of
Cr${}_{1/3}M$S2 ($M$ = Nb or Ta) [21]. The critical temperature $T_{\text{c}}$
of each material matches previous reports [7], although the field required to
induce the forced-ferromagnetic state is lower for both materials [29],
implying a smaller DMI (likely due to different levels of Cr intercalation
affecting the bonding and charge transfer between the $M$ ions, and hence the
magnetic properties). Despite these differences, single-crystal samples of
Cr1/3NbS2 prepared in the same way as our samples show a CSL in Lorentz
transmission electron microscopy measurements [30].

Both Cr1/3NbS2 and Cr1/3TaS2 exhibit three regimes of behavior: (i) $T\lesssim
T^{*}$, (ii) $T^{*}\lesssim T<T_{\text{c}}$, and (iii) $T>T_{\text{c}}$, where
the low-temperature crossover $T^{*}\simeq 20$ K for Cr1/3NbS2, and
$T^{*}\simeq 45$ K for Cr1/3TaS2. In region (iii) the materials are
paramagnetic, transitioning to long-range order when the samples are cooled
into region (ii). On cooling to region (i), changes are seen most clearly in
the phase lag $\phi$ between the applied time-varying magnetic field and the
AC susceptibility response. This suggests dissipative processes [31]
determined by fluctuations across the pseudogap. Cr1/3NbS2 shows a peak in
$\phi$ around $T^{*}$, whereas Cr1/3TaS2 exhibits an increase in this
quantity. For Cr1/3NbS2 the behavior is well described by a Gaussian peak with
a vertical offset, whereas a sigmoid function with a vertical offset is most
appropriate for Cr1/3TaS2. Both functions allow the extraction of a
characteristic temperature $T^{*}$ (the peak location of the Gaussian, or the
inflection point of the sigmoid). Fig. 3(b,d) shows that $T^{*}$ follows an
Arrhenius law, with the gradient giving an activation energy $E_{\text{a}}$
consistent with the full-width half-maxima of the pseudogap $\Delta E=80$–$90$
meV identified in the spin-up DOS. [In Cr1/3NbS2 the gap is slightly smaller
than identified via DFT, and this is the value used to rescale $S$ in Fig.
2(a)]. The observed dynamics are therefore consistent with moment fluctuations
caused by temperature driven transitions across the gap.

Figure 3: (a–b) Example measurements of the phase lag $\phi$ of AC
susceptibility measurements on Cr${}_{1/3}M$S2 ($M$ = Nb or Ta). Solid lines
are fits to the data as described in the text. The central position of these
fits is shown in (c–d).

To further probe dynamic fluctuations around the pseudogap, we performed
pulsed source longitudinal field (LF) $\mu^{+}$SR measurements on
Cr${}_{1/3}M$S2 ($M$ = Nb or Ta), which are sensitive predominantly to MHz to
GHz dynamics (i.e. much higher frequencies than the $<150$ Hz dynamics probed
with AC susceptibility). In a LF $\mu^{+}$SR experiment [32, 33, 21] spin-
polarized muons are implanted in a sample with the initial muon-spin direction
parallel to the externally-applied magnetic field. Time-varying internal
magnetic fields at the stopping site flip the muon spin, reducing the average
muon-spin polarization. The measured quantity of interest is the positron-
decay asymmetry, which is directly proportional to the muon-spin polarization.
Both materials show changes in the initial asymmetry $A_{0}$ near
$T_{\text{c}}$ for all applied fields, indicating the magnetic transition.

We first look at Cr1/3NbS2. A small applied field will decouple the nuclear
contribution to the spectra leaving us sensitive to electronic fluctuations.
An example spectrum for $T<T_{\mathrm{c}}$ is shown in Fig. 4, where
$\mu_{0}H_{\text{ext}}$=2 mT. A reduction in $A_{0}$ is observed at low $T$
[Fig. 4(a)]. Extrapolating the Arrhenius behavior seen with AC susceptibility
to frequencies appropriate for $\mu^{+}$SR
($\simeq\gamma_{\mu}B_{\text{ext}}$, see below), gap driven effects should be
observed around $T^{*}=35$ K (marked as a dashed line in Fig. 4), and it is
below this temperature where changes in $A_{0}$ are observed. A reduction in
$A_{0}$ indicates a larger proportion of the muons stopping in sites where
there is a component of the magnetic field perpendicular to the initial muon-
spin, causing precession too rapid to be observed in the experiment. This
suggests an increase in the internal field at low $T$, as seen with ZF
$\mu^{+}$SR [15] in Fig. 2(e). As the external field is increased, the
internal field has less influence on the field at the muon site, which is
reflected in the reduced resolution of this effect in $A_{0}$ [21].

The LF $\mu^{+}$SR spectra measured on Cr1/3NbS2 at all applied fields were
fitted [34] with
$A\left(t\right)=a_{1}\exp\left(-\sigma^{2}t^{2}\right)\cos\left(\gamma_{\mu}Bt\right)+a_{2}\exp\left(-\lambda
t\right)+a_{\text{b}}\exp\left(-\lambda_{\text{b}}t\right)$, where the first
term accounts for muons in sites where the average internal field is not
parallel to the initial muon-spin, leading to precession in the effective
local field $B$, the second term accounts for muons in a purely longitudinal
average magnetic field, or those where the static field largely cancels,
leading to relaxation from dynamic fluctuations, and the third term accounts
for muons that stop outside the sample. This model should only be applicable
below $T_{\text{c}}$, however we find it describes the data well over the
entire temperature range since it also approximates the expected LF Kubo-
Toyabe behavior above $T_{\text{c}}$ that occurs due to relaxation from
disordered electronic moments. The two regimes are separated by a peak in
$\lambda$ around $T_{\text{c}}$. We find at each field, $\lambda_{\text{b}}$
and $\sigma$ are $T$-independent, therefore we fix them to
$\lambda_{\text{b}}=0.005$ $\mu$s-1 and $\sigma=0.28$, $0.55$, $0.52$ $\mu$s-1
for 2, 4, 10 mT respectively. The $T$-independence of
$\sigma\propto\sqrt{\langle\left(B-\langle B\rangle\right)^{2}\rangle}$
suggests the distribution of magnetic fields at the muon site does not
significantly change. Fitted parameters extracted from the 2 mT measurements
are shown in Fig. 4(b–d).

Figure 4: (a–d) LF $\mu^{+}$SR results for Cr1/3NbS2. (a) Initial asymmetry.
Inset shows an example of the spectra and fit for 2 mT and 15 K. Fitted
parameters using the model presented in the text are shown in (b–d). The solid
blue line in (c) is a guide to the eye. (e–f) Parameters from fitting TF
$\mu^{+}$SR measurements of Cr1/3TaS2 with the model presented in the text.
The inset in (e) shows an example of the spectra for 2 mT at 20 K. High $T$
dashed lines indicate $T_{\text{c}}$, with low $T$ dashed lines indicating
$T^{*}$ (extrapolated from AC susceptibility). For TF $\mu^{+}$SR measurements
a gray region indicates $T^{*}$, accounting for the field dependence.

The main result from our $\mu^{+}$SR measurement of Cr1/3NbS2 is that the
observed MHz dynamics are consistent with the pseudogap-driven magnetism. As
shown in Fig. 4(c), there is a clear peak in $\lambda$ at $T^{*}$ despite no
reported change in magnetic structure [17]. This peak likely occurs due to the
same dynamic fluctuations of the moment size as are responsible for the
behavior observed in AC susceptibility measurements. As $T$ is increased the
fluctuation rate $\nu$ increases in frequency, passing through the muon time-
window. In general, the relaxation rate will be enhanced at $T^{*}$ when
$\gamma_{\mu}B\approx\nu\left(T^{*}\right)$ [21], where
$\gamma_{\mu}=2\pi\times 135.5$ MHzT-1 is the gyromagnetic ratio of the muon.
The observed peak coincides with the prediction from the Arrhenius behavior of
the AC susceptibility when
$\nu\left(T^{*}\right)\approx\gamma_{\mu}B_{\text{ext}}$. This is consistent
with the muon sites most sensitive to dynamics being those where the effective
field is close to the applied field (i.e. the average static field is small).
This interpretation is also supported by the fitted value of the field $B$
[Fig. 4(b)] [21].

We have performed additional LF $\mu^{+}$SR measurements in applied field of 4
mT (where we expect a significant volume fraction of the sample in the CSL
state), and 10 mT (where the solitons in the CSL have large separations).
Measurements at these fields were found to be well described by the same model
as discussed above. At 10 mT, the behavior of all fitted parameters was found
to be very similar to those at 2 mT (see Supplemental Material [21]), however
more significant changes are observed at 4 mT. Most notably, the dynamic
changes at $T^{*}$ observed through $\lambda$ are not seen at 4 mT [Fig.
4(d)], but are recovered at 10 mT. It is possible that since the
topologically-protected CSL has a large energy barrier preventing changes to
the static structure, the more robust magnetism in this phase might prevent
the effects of energy-gap driven magnetism being measurable. Other differences
are seen in the effective local field $B$ [Fig. 4(b)], which at 2 mT remains
fairly constant at slightly less than the applied external field (with a
deviation at $T_{\text{c}}$). However, at 4 mT there is a sizable reduction
compared to the applied field, with $B$ decreasing continuously with
decreasing temperature, perhaps suggesting a continuous evolution in the size
of the static field at the muon site in the CSL. Below $T_{\text{c}}$ at 10
mT, where the CSL is expected to occur with reduced volume fraction, $B$ is
approximately constant at a value significantly reduced compared to the
applied field [21].

We have also performed continuous source ZF $\mu^{+}$SR measurements on
Cr${}_{1/3}M$S2 ($M$ = Nb or Ta). In both cases a fast relaxation of
$A\left(t\right)$ is observed at early times, consistent with the oscillations
detected in single crystal samples of Cr1/3NbS2 [15], however full
oscillations are not resolved in our polycrystalline samples, likely due to
the different, complex field distribution arising from a combination of the
domain and magnetic structures. Continuous source TF $\mu^{+}$SR measurements,
performed on the same Cr1/3NbS2 sample, show very similar early-time behavior
to the LF measurements (as expected for an ordered magnet with a large
internal field compared to the applied field). Pulsed source LF $\mu^{+}$SR
measurements of Cr1/3TaS2 indicate the field at the muon site is considerably
larger than in Cr1/3NbS2, hence TF $\mu^{+}$SR measurements are most useful.
To analyze our continuous source TF $\mu^{+}$SR measurements, the number of
counts in each detector is simultaneously refined to a single
$P_{x}\left(t\right)$, with an appropriate phase offset applied for the
detector geometry. In these measurements we find
$P_{x}\left(t\right)=\sum_{i=1}^{2}a_{i}\cos\left(\gamma_{\mu}B_{i}t+\phi_{i}\right)\exp\left(-\Lambda_{i}t\right)$
well describes the data. The first component captures muons that stop in sites
predominantly sensitive to the large internal field and get rapidly relaxed
out the spectra. In this discussion we focus on the second component [Fig.
4(d–e)], capturing muons that are predominantly sensitive to the applied field
and a small internal field (it is these sites that are most similar to those
to which we are sensitive in the LF $\mu^{+}$SR measurements). We find
$\phi_{2}\simeq 0$, as expected for this class of muon sites. We see an
enhancement of the field at the muon site $B_{2}$ below $T^{*}$, consistent
with similar measurements of Cr1/3NbS2 in Ref. [15]. At high $T$,
$\Lambda_{2}$ appears to scale with the average moment as one would expect,
however the behavior of this parameter is more complicated below $T^{*}$. In
this configuration, $\Lambda_{2}$ is sensitive to both the width of the
distribution of fields at the muon site, and dynamic fluctuations, hence it is
difficult to predict the expected behavior as done for LF $\mu^{+}$SR. Despite
this complexity, the change in behavior below $T^{*}$ suggests the energy
scale set by the pseudogap is once again the important parameter, with the
significant changes below $T^{*}$ likely occurring due to the fluctuating
moment on the muon timescale. Hence, our TF $\mu^{+}$SR results on Cr1/3TaS2
are consistent with the pseudogap picture already discussed.

In conclusion, the results of electronic structure calculations provide an
explanation of the low-$T$ transport and magnetic properties of half-metallic
Cr${}_{1/3}M$S2 ($M$ = Nb or Ta); all published low-$T$ experimental
observations can be explained from this electronic structure, often by
considering electrons transitioning across the small pseudogap in the spin-up
channel. These effects do not require the inclusion of SOC in the calculation,
and do not imply a change of DMI, or a change in the helical length. We have
measured the pseudogap in both materials using AC susceptibility, finding good
agreement with DFT, and observe dynamic transitions across the gap occur over
a wide range of fluctuation frequencies (0.1 Hz to MHz) observed with multiple
techniques. In Cr1/3NbS2, $\mu^{+}$SR suggests that the CSL masks the low-$T$
effects. The similar electronic and magnetic behavior of Cr${}_{1/3}M$S2 ($M$
= Nb or Ta) is consistent with the realization of a CSL in Cr1/3TaS2, a recent
report [9] has provided the first direct evidence of this. Our work
demonstrates that, even in materials where topological magnetism leads to
desirable properties for applications, the magnetic properties cannot be
considered in isolation, and a thorough understanding of the electronic
structure is essential.

Part of this work was carried out at the STFC ISIS Facility, UK
(https://doi.org/10.5286/ISIS.E.RB2010327), and S$\mu$S, Paul Scherrer
Institut, Switzerland; we are grateful for the provision of beamtime. We are
grateful for computational support from both Durham Hamilton HPC and from the
UK national high performance computing service, ARCHER, for which access was
obtained via the UKCP consortium and funded by EPSRC (Grant No. EP/K013564/1).
The project was funded by EPSRC (UK) (Grant Nos.: EP/N032128/1 and
EP/N024028/1). M. N. Wilson acknowledges the support of the Natural Sciences
and Engineering Research Council of Canada (NSERC). Research data from this
paper will be made available via Durham Collections at XXX.

## References

  * Dickinson and Pauling [1923] R. G. Dickinson and L. Pauling, The crystal structure of molybdenite, Journal of the American Chemical Society 45, 1466 (1923).
  * Manzeli  _et al._ [2017] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, 2d transition metal dichalcogenides, Nature Reviews Materials 2, 1 (2017).
  * Wilson and Yoffe [1969] J. A. Wilson and A. D. Yoffe, The transition metal dichalcogenides discussion and interpretation of the observed optical, electrical and structural properties, Advances in Physics 18, 193 (1969).
  * Wang  _et al._ [2012] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Nature Nanotechnology 7, 699 (2012).
  * Choi  _et al._ [2017] W. Choi, N. Choudhary, G. H. Han, J. Park, D. Akinwande, and Y. H. Lee, Recent development of two-dimensional transition metal dichalcogenides and their applications, Materials Today 20, 116 (2017).
  * Moriya and Miyadai [1982] T. Moriya and T. Miyadai, Evidence for the helical spin structure due to antisymmetric exchange interaction in Cr1/3NbS2, Solid State Communications 42, 209 (1982).
  * Kousaka  _et al._ [2016] Y. Kousaka, T. Ogura, J. Zhang, P. Miao, S. Lee, S. Torii, T. Kamiyama, J. Campo, K. Inoue, and J. Akimitsu, Long periodic helimagnetic ordering in Cr$M_{3}$S6 ($M$=Nb and Ta), in  _Journal of Physics: Conference Series_, Vol. 746 (IOP Publishing, 2016) p. 012061.
  * Togawa  _et al._ [2012] Y. Togawa, T. Koyama, K. Takayanagi, S. Mori, Y. Kousaka, J. Akimitsu, S. Nishihara, K. Inoue, A. Ovchinnikov, and J.-i. Kishine, Chiral magnetic soliton lattice on a chiral helimagnet, Physical Review Letters 108, 107202 (2012).
  * Zhang  _et al._ [2021] C. Zhang, J. Zhang, C. Liu, S. Zhang, Y. Yuan, P. Li, Y. Wen, Z. Jiang, B. Zhou, Y. Lei, _et al._ , Chiral helimagnetism and one-dimensional magnetic solitons in a Cr-intercalated transition metal dichalcogenide, Advanced Materials , 2101131 (2021).
  * Yonemura  _et al._ [2017] J. I. Yonemura, Y. Shimamoto, T. Kida, D. Yoshizawa, Y. Kousaka, S. Nishihara, F. J. T. Goncalves, J. Akimitsu, K. Inoue, M. Hagiwara, _et al._ , Magnetic solitons and magnetic phase diagram of the hexagonal chiral crystal CrN3S6 in oblique magnetic fields, Physical Review B 96, 184423 (2017).
  * Mayoh  _et al._ [2021] D. A. Mayoh, J. Bouaziz, A. E. Hall, J. B. Staunton, M. R. Lees, and G. Balakrishnan, submitted (2021).
  * Lancaster [2019] T. Lancaster, Skyrmions in magnetic materials, Contemporary Physics 60, 246 (2019).
  * Ghimire  _et al._ [2013] N. J. Ghimire, M. A. McGuire, D. S. Parker, B. Sipos, S. Tang, J.-Q. Yan, B. C. Sales, and D. Mandrus, Magnetic phase transition in single crystals of the chiral helimagnet Cr1/3NbS2, Physical Review B 87, 104403 (2013).
  * Bornstein  _et al._ [2015] A. C. Bornstein, B. J. Chapman, N. J. Ghimire, D. G. Mandrus, D. S. Parker, and M. Lee, Out-of-plane spin-orientation dependent magnetotransport properties in the anisotropic helimagnet Cr1/3NbS2, Physical Review B 91, 184401 (2015).
  * Braam  _et al._ [2015] D. Braam, C. Gomez, S. Tezok, E. V. L. de Mello, L. Li, D. Mandrus, H.-Y. Kee, and J. E. Sonier, Magnetic properties of the helimagnet Cr1/3NbS2 observed by $\mu$SR, Physical Review B 91, 144407 (2015).
  * Togawa  _et al._ [2019] Y. Togawa, J. Kishine, P. A. Nosov, T. Koyama, G. W. Paterson, S. McVitie, Y. Kousaka, J. Akimitsu, M. Ogata, and A. S. Ovchinnikov, Anomalous temperature behavior of the chiral spin helix in CrN3S6 thin lamellae, Physical Review Letters 122, 017204 (2019).
  * Borisenko  _et al._ [2017] E. B. Borisenko, V. A. Berezin, N. N. Kolesnikov, V. K. Gartman, D. V. Matveev, and O. F. Shakhlevich, Structural and magnetic ordering of CrNb3S6 single crystals grown by gas transport method, Physics of the Solid State 59, 1310 (2017).
  * Clark  _et al._ [2005] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson, and M. Payne, First principles methods using CASTEP, Z. Kristall. 220, 567 (2005).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Physical Review B 13, 5188 (1976).
  * [21] See supplemental material for computational and experimental details, sample synthesis and characterization, and complementary $\mu^{+}$SR measurements.
  * Coey and Sanvito [2004] J. Coey and S. Sanvito, Magnetic semiconductors and half-metals, Journal of Physics D: Applied Physics 37, 988 (2004).
  * Coey and Venkatesan [2002] J. Coey and M. Venkatesan, Half-metallic ferromagnetism: Example of CrO2, Journal of Applied Physics 91, 8345 (2002).
  * Madsen and Singh [2006] G. K. Madsen and D. J. Singh, BoltzTraP. a code for calculating band-structure dependent quantities, Computer Physics Communications 175, 67 (2006).
  * Mito  _et al._ [2019] M. Mito, H. Ohsumi, T. Shishidou, F. Kuroda, M. Weinert, K. Tsuruta, Y. Kotani, T. Nakamura, Y. Togawa, J. Kishine, _et al._ , Observation of orbital angular momentum in the chiral magnet CrNb3S6 by soft x-ray magnetic circular dichroism, Physical Review B 99, 174439 (2019).
  * Franke  _et al._ [2018] K. J. A. Franke, B. M. Huddart, T. J. Hicken, F. Xiao, S. J. Blundell, F. L. Pratt, M. Crisanti, J. A. T. Barker, S. J. Clark, A. Štefančič, _et al._ , Magnetic phases of skyrmion-hosting GaV4S8-ySey ($y$=0, 2, 4, 8) probed with muon spectroscopy, Physical Review B 98, 054428 (2018).
  * Momma and Izumi [2011] K. Momma and F. Izumi, VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data, Journal of applied crystallography 44, 1272 (2011).
  * Bonfà  _et al._ [2018] P. Bonfà, I. J. Onuorah, and R. De Renzi, Introduction and a quick look at MUESR, the Magnetic structure and mUon Embedding Site Refinement suite, in  _Proceedings of the 14th International Conference on Muon Spin Rotation, Relaxation and Resonance ( $\mu$SR2017)_ (2018) p. 011052.
  * Han  _et al._ [2017] H. Han, L. Zhang, D. Sapkota, N. Hao, L. Ling, H. Du, L. Pi, C. Zhang, D. G. Mandrus, and Y. Zhang, Tricritical point and phase diagram based on critical scaling in the monoaxial chiral helimagnet Cr1/3NbS2, Physical Review B 96, 094439 (2017).
  * [30] A. E. Hall _et al._ , in preparation .
  * Topping and Blundell [2018] C. V. Topping and S. J. Blundell, AC susceptibility as a probe of low-frequency magnetic dynamics, Journal of Physics: Condensed Matter 31, 013001 (2018).
  * Blundell [1999] S. J. Blundell, Spin-polarized muons in condensed matter physics, Contemporary Physics 40, 175 (1999).
  * Blundell  _et al._ [2021] S. J. Blundell, R. De Renzi, T. Lancaster, and F. L. Pratt,  _Muon Spectroscopy: An Introduction_ (Oxford University Press, Oxford, 2021).
  * Pratt [2000] F. L. Pratt, WIMDA: a muon data analysis program for the Windows PC, Physica B: Condensed Matter 289, 710 (2000).