# Masked Matrix Multiplication for Emergent Sparsity

Brian Wheatman<EMAIL_ADDRESS>Johns Hopkins University ,  Meghana
Madhyastha<EMAIL_ADDRESS>Johns Hopkins University and  Randal Burns
<EMAIL_ADDRESS>Johns Hopkins University

###### Abstract.

Artificial intelligence workloads, especially transformer models, exhibit
emergent sparsity in which computations perform selective sparse access to
dense data. The workloads are inefficient on hardware designed for dense
computations and do not map well onto sparse data representations. We build a
vectorized and parallel matrix-multiplication system $(A\times B=C)$ that
eliminates unnecessary computations and avoids branches based on a runtime
evaluation of sparsity. We use a combination of dynamic code lookup to adapt
to the specific sparsity encoded in the B matrix and preprocessing of sparsity
maps of the A and B matrices to compute conditional branches once for the
whole computation. For a wide range of sparsity, from 60% to 95% zeros, our
implementation performs fewer instructions and increases performance when
compared with Intel MKL’s dense or sparse matrix multiply routines. Benefits
can be as large as 2 times speedup and 4 times fewer instructions.

Matrix Multiplication, Sparse Computation

††copyright: none††ccs: Mathematics of computing Computations on
matrices††ccs: Computing methodologies Shared memory algorithms††ccs:
Computing methodologies Machine learning algorithms††ccs: Theory of
computation Algorithm design techniques

##  1\. Introduction

Recent results have observed that sparsity is an emergent phenomena in AI and
leveraging this sparsity is a powerful tool in reducing computational cost.
The Lazy Neuron Phenomena (Li et al., 2023) shows that the activations of ReLU
in multi-layer perceptrons in transformer models have greater than 90%
sparsity. They also observed that all neurons are active for some activations.
Essentially, the data are dense, all neurons are needed, and their usage is
sparse. Contextual sparsity (Liu et al., 2023) shows that using activations to
dynamically select a sparse subset of a model can approximate the output of a
dense model when using only 5% of MLP parameters and 20% of attention heads.
The potential savings are huge and neither approach realizes speedup
proportional to sparsity. The problem lies in mapping sparse computation to
the dense matrix multiplication primitives supported by hardware. The Lazy
Neuron authors “emphasize that this result is far from fully realizing the
benefit of sparse activation, due to a lack of hardware support for sparse
computation.” (Li et al., 2023)

We present a matrix multiplication system that better maps sparse computations
to hardware designed for dense, vectorized computation. We call this Masked
Matrix Multiplication (MMM), because it operates on data that are stored
densely but are selected sparsely at runtime. Specifically, $M(X)\times
N(Y)=C$ where $X$ and $Y$ are dense matrices and $M$ and $N$ are masking
functions that zero out large fractions of the input and are evaluated at
runtime. For ease of presentation, we define $A=M(X)$ and $B=N(Y)$. To be
effective, MMM requires sparsity in both matrices and block-sparsity in $B$.
Block sparsity means that there are small, dense sequential blocks of non-
zeros on which we perform vector operations. If only one matrix is sparse,
simpler techniques, such as Intel MKL’s sparse-dense approach are more
efficient.

The system works best for the types of sparsity observed in transformer models
and at the fraction of non-zeros 60-95% seen in the data. The specific
combination of bit-level and block-sparsity occurs in the Deja Vu system (Liu
et al., 2023). For the most part, these dynamic sparse-sparse multiplications
have not been explored, because they are no more efficient than sparse-dense
or dense-dense multiplications without our techniques. We show that these
computations exist and that they can be leveraged to improve performance. We
expect that as machine-learning workloads continue to explore sparsity, MMM
will become a tool to reduce cost, power, and time. Sparsity exists or can be
leveraged in many settings, including switch transformers (Fedus et al.,
2021), sparse transformers (Jaszczur et al., 2021), and reformers with sparse
attention (Kitaev et al., 2020).

We use dynamic code execution to adapt to the sparsity encoded in the $B$
matrix and single-pass preprocessing of $A$ and $B$ to minimize branching
instructions. The challenge lies in that there is sparsity in both matrices,
also called dual-sided sparsity (Wang et al., 2021). With a single sparse
matrix, one can just skip elements of that sparse matrix in the outer loop and
perform dense computations in the inner loops. Our system deals with sparsity
at a finer-grain in inner loops by selecting code from a lookup table that
matches the block-sparse patterns seen during preprocessing. Preprocessing is
parallel, cache efficient, and vectorizable.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$0$$20$$40$$60$$80$2x
speedupSparsitytime (ms)MMMmkl densemkl sparse Figure 1. Runtime for a
2048x2048 matrix multiplication with varying sparsity. MMM outperforms Intel
MKLs best algorithm for intermediate levels of sparsity, providing 2 times
speedup at 80% zeros.
$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$1$$2$$3$Sparsityinstructions(B)MMMmkl densemkl
sparse Figure 2. Number of instructions to perform a 2048x2048 matrix
multiplication with varying sparsity. MMM reduces instructions be a factor of
4 at 70% zeros.

Overall, the system provides speedup over Intel MKL’s sparse-sparse, sparse-
dense, and dense-dense matrix multiplication for a wide-range of sparsity
(Figure 1), exhibiting two times speedup at 80% sparsity in both matrices.
This is a summary result in that there are many parameters that affect
performance, including $A$ sparsity, $B$ sparsity, matrix size, and
architecture. The total number of instructions shows how MMM avoids
computation on 0s. MMM can perform up to 4 times fewer instructions than MKL
sparse or dense (Figure 2). The combination of results also indicates that
Intel MKL is realizing higher instructions per cycle (IPC) than MMM. This
paper will not answer the question: Can the gap between performance and
instructions be narrowed in MMM with respect to Intel MKL? We believe that the
sophisticated locality, blocking, and register management used in Intel MKL
(Wang et al., 2014) are compatible with MMM and would close this gap. However,
our implementation does not explore this and a much greater engineering effort
would be needed.

We present results for Intel CPUs only. GPUs are largely preferred for AI
training and a GPU implementation remains as future work. GPUs have the same
vector computing and sequential memory properties that MMM leverages. CPUs
remain important for inference: 70 percent of inference in the datacenter,
including the hyperscalers and cloud builders, runs on Intel Xeon CPUs
(Morgan, 2023; Corporation, 2022). Inference at the IoT edge largely uses CPUs
(Hao et al., 2023).

##  2\. Background

We describe the matrix multiplication on which MMM builds, starting with the
simple implementation of dense multiplication for row-major matrices. The
basic algorithm (Figure 3) uses three-nested loops that iterate over the
inputs $A$ and $B$ sequentially and writes partial sums to the output matrix
$C$ sequentially.

⬇

void MatrixMultiplySimple(Matrix &C,

const Matrix &A, const Matrix &B) {

for (int i = 0; i < C.rows; ++i)

for (int k = 0; k < B.rows; ++k)

for (int j = 0; j < C.cols; ++j)

C[i, j] += A[i, k] * B[k, j];

}

Figure 3. Simple dense matrix multiplication.

Recursive calls to the simple algorithm convert the computation into a cache-
oblivious algorithm that fits into all levels of cache simultaneously. This
breaks each matrix into 4 sub-matrices and performs 8 matrix multiplications
among the sub-matrices as seen in Figure 4. Recursive division continues until
a fixed size (threshold) that fits into processor caches.

⬇

void recursiveMM(Matrix &C, const Matrix &A,

const Matrix &B) {

if (size <= threshold) {

MatrixMultiplySimple(C, A, B);

return;

}

recursiveMM(C.UL(), A.UL(), B.UL());

recursiveMM(C.LL(), A.LL(), B.UL());

recursiveMM(C.UR(), A.UL(), B.UR());

recursiveMM(C.LR(), A.LL(), B.UR());

recursiveMM(C.UL(), A.UR(), B.LL());

recursiveMM(C.LL(), A.LR(), B.LL());

recursiveMM(C.UR(), A.UR(), B.LR());

recursiveMM(C.LR(), A.LR(), B.LR());

}

Figure 4. Simple recursive matrix multiplication implementation. UL, UR, LL,
and LR stand for the upper left, upper right, lower left and lower right sub-
matrices.

This formulation implements the $O(n^{3})$ algorithm for $n\times n$ matrices.
Larger matrices can benefit from Strassen’s $O(n^{2.807})$ (Strassen, 1969)
algorithm that avoids redundant computation at the expense of numerical
stability. Strassen’s uses a recursive reduction and can use MMM at the base
of recursion. Algorithms with lower asymptotic complexity ($O(n^{2.371552})$
(Williams et al., 2023)) are not used in practice. MMM does not implement
Strassen’s and is not optimized for larger matrices. Our effective performance
range lies between $1024^{2}$ and $8192^{2}$. However, Strassen’s algorithm
uses a similar recursive structure and requires these smaller block matrix
multiplications.

Sparse Generalized Matrix Multiplication (SpGEMM) (Gustavson, 1978) multiplies
two matrices that are stored in a compressed format. The goal is to have the
computational work scale with the number of non-zeros in the matrix, rather
than the matrix size. SpGEMM stores the indexes and values of the non zeros
compactly, this reduces computation, but loses some of the efficiency gained
from aligned vector operations. For this reason, SpGEMM does not outperform
dense dense matrix multiplication until high levels of sparsity. This is true
when SpGEMM matrices are already encoded sparsely, i.e. there is no conversion
cost. MMM will adapt some data structures from the SpGEMM literature.

Random SparsityBlock Random SparsityRandom Column SparsityBlock Random Column
Sparsity Figure 5. Types of sparsity. In random sparsity, any bit may be non-
zero. Block-random sparsity assumes aligned, contiguous regions of each row
are sparse or not sparse. The degree of block-random sparsity may be lower
than the fraction of sparsity when 0s are not fully aligned. We also show
random column and random block column sparsity. Block column patterns arise in
Deja Vu (Liu et al., 2023).

##  3\. Masked Matrix Multiplication

We provide improvements to matrix multiplication for emergent and contextual
sparsity by detecting zero-values at runtime to avoid unnecessary
computations. The challenge lies in detecting zero values efficiently while
preserving the multicore and vector parallelism of the underlying hardware.
The implementation sits in between dense matrix multiplication and SpGEMM.
Dense matrix multiplication maximizes memory throughput and hardware
parallelism, but computes all elements. SpGEMM only computes non-zero
elements, but uses data structures that require branching instructions and
perform more random memory accesses.

We call this Masked Matrix Multiplication (MMM) because the sparsity arises at
runtime. The underlying data are dense so that sparsity cannot be evaluated at
compile time and the data cannot be stored in sparse formats. We express this
as $M(X)\times N(Y)=C$ in which $M$ and $N$ are masking functions over dense
arrays $X$ and $Y$, we define $A=M(X)$ and $B=N(Y)$. Performance benefits
arise when both $A$ and $B$ are sparse; if only one is sparse, simpler
techniques are more efficient. More specifically, we require random bit
sparsity in the $A$ matrix evaluated in the outer loop and block random
sparsity in $B$ (Figure 5).

⬇

for (int i = 0; i < C.rows; ++i) {

for (int k = 0; k < B.rows; ++k) {

auto A_val = A[i, k];

if (A_val == 0) {

continue;

}

for (int j = 0; j < C.cols; ++j) {

C[i, j] += A_val * B[k, j];

}

}

}

Figure 6. A simple optimization to skip computation when A is sparse. The if
statement skips over rows of B when the corresponding element of A is 0.

Random Sparsity in $A$. We evaluate elements of $A$ one at a time in our
middle loop and avoid computation when elements are zero. This is a standard
technique. It is customary when choosing loop orders for tensors operations
(Kjolstad et al., 2017) or in hardware implementations of sparse accelerators
(Wu et al., 2021) to lift the sparsest matrix to the outer loop for exactly
this reason. A branch instruction, if (A_val == 0), avoids evaluation of the
inner loop (Figure 6).

This skips an entire row of $B$ and an entire row of partial writes to $C$. In
later steps, we will optimize the evaluation of the branch instruction.

Block-Random Sparsity in $B$. Eliminating computation in the inner loop is
more complicated, because efficient implementations preserve vector processing
of multiple elements. We cannot eliminate computation on a single 0 element.
The compiler turns the inner loop into a vector computation for which the
psuedocode can be seen in Figure 7. The inner operation is a vector fused
multiply add (FMA) and most architectures support this instruction natively,
including our targets Intel AVX-2 and AVX-512. There are masked variants of
FMA that skip elements in the vector. However, masked operations have the same
cost as FMA on the entire vector.

⬇

for (index_t j = 0; j < C.cols; j += vector_size) {

C[i, j:j+vector_size] +=

{A_val,…,A_val} * B[k, j:j+vector_size];

}

Figure 7. The base vectorized row code to help compute matrix multiplication.
Figure 8. Preprocessing patterns of length 4 in $B$. Grayed entries are non-
zero. A compare of each vector extracts a 0/1 bit indicating whether any words
in the vector are non-zero. A pattern worth of bits is used to lookup a
function pointer which is stored in an array or list of functions pointers.
The function pointer array is used during multiplication to directly reference
code that avoids computation on 0 vectors.

The goal is to identify entire vectors of 0s on which can skip computation. A
simple approach checks that elements $B[k,j:j+vector\\_size]$ are all zero
prior to performing the vector FMA. Unfortunately, checking that the vector is
zeros costs about as many cycles as the FMA. So, no savings are possible.

We will leverage that the rows of $B$ are used $C.cols$ times to preprocess
patterns of sparsity, paying a one-time cost to encode sparsity and then incur
only one memory reference of overhead per pattern to avoid computing on
vectors of all zeros. A pattern describes whether each vector in a fixed-width
set of vectors within a row are all zeros. At compile time, we generate
optimized code for all possible patterns that we store in a code table.

Figure 8 gives a schematic overview of the computation. Patterns consist of
multiple vectors. We preprocess the rows of the $B$ matrix to extract patterns
that we encode as bit strings. In this example, patterns of 4 in which $0$
indicates all zero values and $1$ indicates some non-zero values. Using the
bit string, we build an array of function pointers to the optimized code for
each pattern. During matrix multiplication, we will iterate over the $B$
matrix and the function pointers in lock step to select the optimized code for
each pattern.

Compile-Time Code Generation: We build a code table for all possible vector
patterns in which the instructions compute the non-zero vectors and skip over
the zero vectors. There are no branch instructions. Figure 9 gives the full
set of AVX-2 256-bit vector instructions for patterns of size 8 on 32-bit
floats. Each specific subversion takes the lines A-H corresponding to the non-
zero pattern bits, e.g. pattern 01010000 would take lines BDBDBD which would
conduct load, fused multiply add, and store on only the second and fourth
vector in the pattern. All patterns are listed in a code table and we build a
map of function pointers for each pattern. Code is not of fixed size. With
AVX2, functions are as small as 1 byte (the empty case 00000000) to 166 bytes
(the complete case 1111111) for a total size of 21376 bytes.

Choosing the size of the pattern has tradeoffs that affect compile time,
instruction cache performance, and the granularity at which we can detect
sparsity. The number of patterns grows as $2^{size}$. Our implementation has 8
bits which results in 256 patterns. Large patterns have larger code tables,
take more compile time, and perform fewer memory accesses. Smaller patterns
tend to fit into L1i cache so that the memory indirection is less costly. We
found that a pattern size of 8 works well on recent generations of Intel
processors. Generating the code table takes a few seconds and performance
nearly matches that when using 16 bit patterns111Generating the code table for
a pattern size of 16 can take more than an hour and starts running into
implementation limits for some compilers..

⬇

A: vmovups ymm1, ymmword ptr [rsi]

B: vmovups ymm2, ymmword ptr [rsi + 32]

C: vmovups ymm3, ymmword ptr [rsi + 64]

D: vmovups ymm4, ymmword ptr [rsi + 96]

E: vmovups ymm5, ymmword ptr [rsi + 128]

F: vmovups ymm6, ymmword ptr [rsi + 160]

G: vmovups ymm7, ymmword ptr [rsi + 192]

H: vmovups ymm8, ymmword ptr [rsi + 224]

A: vfmadd213ps ymm1, ymm0, ymmword ptr rdi]

B: vfmadd213ps ymm2, ymm0, ymmword ptr rdi + 32]

C: vfmadd213ps ymm3, ymm0, ymmword ptr rdi + 64]

D: vfmadd213ps ymm4, ymm0, ymmword ptr rdi + 96]

E: vfmadd213ps ymm5, ymm0, ymmword ptr rdi + 128]

F: vfmadd213ps ymm6, ymm0, ymmword ptr rdi + 160]

G: vfmadd213ps ymm7, ymm0, ymmword ptr rdi + 192]

H: vfmadd213ps ymm8, ymm0, ymmword ptr [rdi + 224]

A: vmovups ymmword ptr [rdi], ymm1

B: vmovups ymmword ptr [rdi + 32], ymm2

C: vmovups ymmword ptr [rdi + 64], ymm3

D: vmovups ymmword ptr [rdi + 96], ymm4

E: vmovups ymmword ptr [rdi + 128], ymm5

F: vmovups ymmword ptr [rdi + 160], ymm6

G: vmovups ymmword ptr [rdi + 192], ymm7

H: vmovups ymmword ptr [rdi + 224], ymm8

ret

Figure 9. Full pattern multiplication code. Specific patterns will use a
subset of lines A-H. At start ymm0 has the element from A in all positions,
rsi is the address of the row of B, rdi is the address to the row of C.

Preprocessing: We take a single pass over each sub-matrix of the matrix
evaluating sparsity for each pattern and build an array of function pointers
to the optimized code. A sub-matrix is the unit and the bottom of the
recursion that actually computes the matrix multiplication as seen in Figure
4.

This process is efficient, because it accesses memory sequentially,
comparisons are vector aligned, and it writes function pointers sequentially.
We preprocess sub-matrices of B in the same order that we multiply in the next
phase so that function pointers are stored in the order that they will be
accessed. The preprocessing pseudocode (Figure 10) looks at a pattern at a
time (outer for loop) and does a compare to determine if each vector is all 0s
(inner for loop) to extract a bit for each vector. The pattern_index bitstring
is used to lookup the function pointer and write it to the fps function
pointer array. This computation is fully aligned. This is in contrast to the
computation needed to convert a dense matrix to a sparse matrix format (Buluç
and Gilbert, 2012) that requires packing elements together. Conversion to
sparse formats is difficult to vectorize.

This preprocessing code is called a single time for each sub-matrix at a cost
of $O(n^{2})$ for $n\times n$ dense matrices. The cost of preprocessing
becomes negligible as the matrix size grows. Empirically we find that this
crossover point to be around $1024^{2}$. This code can also be trivially
parallelized over the different sub-matrices.

⬇

# function pointers for sub-matrix

fps = []

for r in range(rows):

for c in range(0, cols, vector_size*pattern_size):

pattern = 0

pattern_index = 1 << pattern_size

for p_index in range(c,

c+vector_size*pattern_size, vector_size):

# vector comparison

if submatrix(r, p_index:p_index+vector_size) == 0:

pattern |= pattern_index

pattern_index >>= 1

fps[submatrix_index].append(

function_pointer_table[pattern])

Figure 10. Preprocessing function pointer arrays that encode pattern sparsity.
All computations are vector aligned.

Matrix Multiplication (Figure 11) loops over the elements by regions of size
vector_size * pattern_size and calls the corresponding optimized function
fps[fp_index]. This pays a single indirect call for each region. The optimized
functions skip over all zero vectors using no branch instructions.

⬇

step_size = vector_size * pattern_size

for j in range(A.rows):

fp_index = 0

for k in range(B.rows):

if A[i, k] == 0:

# skip region of fps corrosponding to the zero element of A

fp_index += B.cols / step_size

continue

for j in range(0, B.cols, step_size):

fps[fp_index](&C[i,j], A[i,k], &B[k,j])

fp_index+=1

Figure 11. MMM Pseudocode. The MMM computation has the same iteration
structure as the dense matrix multiply code, but skips zeros with the
optimized row functions.

### Optimizations

Our first optimization helps avoid branch instructions as $A$ gets more
sparse. We preprocess $A$ on $256^{2}$ blocks to encode sparse indexes (Figure
12). We build a fully aligned data structure of the same size that encodes the
count of non-zero elements and their offsets. This allows us to replace the
branch if A[i, k] == 0: with a for loop over the-nonzero elements so that the
number of conditional branches scales with the number of non-zeros, rather
than with the number of columns. This is a memory aligned version of the
compressed sparse row format (Buluç and Gilbert, 2012).

Another optimization increases the size of the base vector on which we
operate. It is possible to have each pattern represent more than 1 vector,
e.g. 2 or 4 vectors. Doubling the base vector size halves the number of calls
to the optimized sub-matrix computation, but also doubles the size of the
smallest block pattern we capture. We see benefits here for some workloads.

Figure 12. Blocked sparse row preprocessing of sub-matrices of $A$. In the
input (left), grayed entries are $V$ indicate a value at a location. The
output (right) encodes the count of non-zero elements and their offsets. The
offsets in a full row (count of 8) do not need to be stored explicitly.

##  4\. Evaluation

We explore the potential benefits of MMM in comparison with dense-dense,
sparse-dense, and sparse-sparse multiplication in Intel’s MKL library. There
is an extremely complicated performance landscape among matrix size,
variations is sparsity in A and B matrices, parallelism, and architecture.
Experiments against synthetic data on a consumer grade machine vary these
parameters to characterize this landscape and show the regions in which MMM
outperforms all MKL techniques. These will be mid-ranges of sparsity, mid-size
matrices, and multi-core parallelism. We extend these results to server-class
machines and on data sets extracted from Deja Vu (Liu et al., 2023).
Experiments show improvements on select cases. We do not implement hardware
optimizations (register reuse, pipelining) beyond vectorization nor algorithms
for reducing complexity on large matrices, which limits the parameter values
on which we improve performance.

#### System Setup

MMM is implemented in C++ for shared memory parallel x86 CPUs. We evaluate on
a variety of x86 cpus from Intel (Table 1). We focus on Intel CPUs to compare
with Intel’s math kernel library, which has been highly optimized. The code
was hand vectorized using Intel AVX intrinsics. Two of the machines only
support AVX2 while the others support AVX512. The code is parallelized with
Cilk (Intel Corporation, 2010) and the Tapir (Schardl et al., 2019; Schardl
and Lee, 2023) branch of the LLVM (Lattner, 2002; Lattner and Adve, 2004)
compiler (version 16).

Machine | CPU | avx | cpus | threads
---|---|---|---|---
i9 | Intel Core i9-13900H | AVX2 | 1 | 20
c4 | Intel Xeon E5-2666 v3 | AVX2 | 2 | 36
c5 | Intel Xeon Platinum 8275CL | AVX512 | 2 | 96
c6 | Intel Xeon Platinum 8375C | AVX512 | 2 | 128
Table 1. Computing platforms.

#### Random Matrices

Randomly generated matrices vary the random sparsity in $A$, fraction of
zeros, and block random sparsity in $B$, fraction of 0 valued vectors. Varying
these quantities allows us to build maps that demonstrate speedup regions for
MMM.

Figure 1 shows the relative speedup to the Intel MKL on a consumer-grade
machine (i9-13900H). This experiment multiplies $2048^{2}$ matrices using all
threads. We report the median of 10 trials. We then run these same
multiplications using the valgrind simulator and find that MMM executes
significantly fewer instructions than either MKL approach in this range
(Figure 2). The reduction in instructions (4 times) exceeds the reduction in
runtime (2 times), indicating that dense MKL realizes more instructions per
cycle.

Figure 13 shows the same speedup experiment over all values of $A$ and $B$
sparsity. MMM provides speedup in configurations with $A$ sparsity higher than
0.5 and less than 0.95. The interesting region is when $B$ is also between 0.5
and 0.95 where we avoid computations in both matrices and achieve a speedup of
two times in many configurations. The region of high A sparsity and low B
sparsity is less significant as all the benefit comes from skipping rows.
Looking at the multiplication method chosen by Intel MKL reveals that MMM fits
in an intermediate space between Intel MKLs different methods (Figure 14). Our
best speedups lie in the sparsity region in which MKL’s optimizer switches
from dense-dense, to sparse-dense, to sparse-sparse.

$0$$0.2$$0.4$$0.6$$0.8$$1$1.8.6.4.20B SparsityA Sparsity Figure 13. Heatmap of
the performance of MMM versus Intel MKL’s best algorithm. Numbers are speedup
and green indicates better performance. $0$$0.2$$0.4$$0.6$$0.8$$1$1.8.6.4.20B
SparsityA Sparsity Figure 14. Intel’s fastest multiplication algorithm–dense-
dense (DD), sparse-dense(SD), and sparse-sparse(SS)–as a function of sparsity.

Our next experiment varies the size of the matrix to characterize
preprocessing overheads and scaling. Table 2 shows the performance on
different sized matrices at 80% sparsity in both $A$ and $B$. Once the
matrices are at least 2048 in size we achieve substantial speedup. For small
matrices, the $O(n^{2})$ preprocessing to build function pointers consumes a
larger fraction of overall $O(n^{3})$ computing time and reduces MMM’s
speedup. We note that the MMM approach only impacts the base multiplication of
a single block and does not affect the high level algorithm or the
parallelization strategy.

n | MMM | mkl dense |  | mkl dense
---
MMM
| mkl sparse |  | mkl sparse
---
MMM
256 | 737 | 116 | 0.16 |  | 207 | 0.28
512 | 1930 | 766 | 0.40 |  | 688 | 0.36
1024 | 5922 | 6188 | 1.04 |  | 8311 | 1.40
2048 | 26613 | 52787 | 1.98 |  | 47784 | 1.80
4096 | 214471 | 327315 | 1.53 |  | 320535 | 1.49
8192 | 1608555 | 2511919 | 1.56 |  | 2877960 | 1.79
Table 2. Runtime (s) and speedup that MMM achieves with different matrix sizes
for a sparsity of .8. Larger matrices amortize preprocessing costs in MMM.

We present the scalability of MMM from 1 to 20 threads (Figure 15). MMM scales
better. However, MKL starts slower at 1 thread. MMM is faster than MKL only at
8 or more threads. We hypothesize that MMM requires concurrency to overlap
memory accesses and computation. That is, that sparse regions lead to memory-
bound, low-compute periods in each thread.

$0$$5$$10$$15$$20$$2$$4$$6$ThreadsspeedupMMMmkl densemkl sparse Figure 15.
Multicore scalability of MMM on $8192^{2}$ matrices.

#### Server Architectures

We run experiments on three generations of Intel server processors to show the
potential of MMM on different architectural and memory configurations. These
experiments both (1) show that MMM reduces instruction counts and (2) verify
that there are parameters and configurations on which MMM provides speedup.
There are many limitations to using whole, multi-socket machines with MMM’s
current implementation. Most importantly, MMM doesn’t have optimizations for
large matrices. To create a workload that uses the entire machine and keeps
matrices mid-size, we run $n$ independent processes, each of which runs a few
matrix multiplies and sum up the total time. This roughly corresponds to a
transformer/MLP workload in which the dimensionality of any operation is
limited by the layer. MMM does not work well on single massive matrix
multiplies that would use the entire machine.

Table 3 shows that MMM achieves reasonable speedups on all three generations
of CPUs and does better when computing floats than doubles. We also measure
the number of instructions each computation takes and report the ratios. MMM
performs many fewer instructions that MKL dense in all cases. We even perform
fewer instructions than MKL sparse in all cases. This is just the
multiplication and does not include conversion costs, i.e. MKL starts with
sparse matrices. Instruction counts were measured using performance counters
on the c5 and c6 and valgrind on the C4 (did not finish in larger cases).
Additionally, we find that the instruction count reduction is substantially
higher on the c4 in AVX2 with 256-bit vectors. MMM does not yet fully exploit
the additional registers available to AVX-512.

Machine | Size |  | data
---
type
Sparsity | speedup |  | I dense
---
I MMM
| I sparse
---
I MMM
c4 | 2048 | float | 0.88 | 1.26 | 6.75 | 2.18
c4 | 2048 | double | 0.86 | 1.33 | 7.69 | 2.14
c4 | 4096 | float | 0.88 | 1.27 | DNF | DNF
c4 | 4096 | double | 0.85 | 1.25 | DNF | DNF
c5 | 2048 | float | 0.9 | 1.12 | 5.02 | 1.35
c5 | 2048 | double | 0.88 | 1.29 | 6.06 | 1.29
c5 | 4096 | float | 0.92 | 1.21 | 7.47 | 1.10
c5 | 4096 | double | 0.88 | 1.35 | 7.67 | 1.18
c6 | 2048 | float | 0.91 | 0.80 | 4.66 | 1.14
c6 | 2048 | double | 0.89 | 1.11 | 6.68 | 1.13
c6 | 4096 | float | 0.93 | 0.82 | 8.13 | 1.13
c6 | 4096 | double | 0.88 | 1.22 | 7.83 | 1.21
Table 3. Best speedup (varying sparsity) over MKL for each machine, size, data
type combination. $\frac{\text{I dense}}{\text{I MMM}}$ and $\frac{\text{I
sparse}}{\text{I MMM}}$ gives instruction count ratios for MKL dense MKL
sparse respectively. This was measured using valgrind for the c4 and
performance counters for c5 and c6.

#### Matrices from ML workloads

We measure MMM’s performance on matrices gathered from transformer inference
after sparsification of the MLP layers in Deja Vu (Liu et al., 2023). We used
the 1.3 billion parameter Open Pre-Trained Transformer Language (OPT 1.3b)
from Hugging Face (Zhang et al., 2022). The MLP layer computes:

(1) $MLP_{SM}(Y)=\sigma(yW^{1}_{SM})(W^{2}_{SM})^{T}$

in which $\sigma$ is the GeLU activation function. $W^{1},W^{2}\in R^{d\times
4d}$ and $d=2048$. $SM\subseteq[4d]$ is the set of neurons that are activated.
The remaining neurons are zero resulting in column sparse matrices. The first
matrix $\sigma(yW^{1}_{SM})$ is column sparse. The second matrix
($(W^{2}_{SM})^{T}$) is row sparse since it is the transpose of a column
sparse matrix.

This is one specific instance of sparse-sparse multiplication in Deja Vu.
There are many sparsified matrices that we do not use, including column-sparse
transformer heads and random sparse activations. However, the block sparse
attention heads are multiplied by a projection matrix that is dense. The
projection matrix cannot be sparsified without losing information because it
combines information from all the attention heads. We expect many more
operations to have two-sided sparsity, but that requires integrating sparsity
in the model design (Wang et al., 2021).

MMM achieves speedup on many of the matrices from Deja Vu’s MLP layers. The
degree of speedup depends on the sparsity parameter, which Deja Vu varies.
Table 4 shows three examples at a favorable sparsity (88%) within Deja Vu’s
target range.

A sparsity | B sparsity | MMM | MKL SD | speedup
---|---|---|---|---
0.88 | 0.88 | 13102 | 19678 | 1.5
0.88 | 0.88 | 14247 | 20104 | 1.4
0.88 | 0.88 | 13045 | 20065 | 1.5
Table 4. Speedup over MKL sparse on three matrices from ML workloads. Sparse-
dense was the fastest MKL approach.

##  5\. Related Work

Nvidia sparse tensor cores (Mishra et al., 2021) pursue a similar goal of
mapping sparse computations onto dense matrix hardware. They map 50% sparse
computations onto dense hardware using a 2:4 strategy, operating on two non-
zero elements of 4 adjacent elements. This is effective on a wide-variety of
models that can be sparsified so that a 2:4 encoding preserves the accuracy of
the dense model. The 2:4 principle has been generalized to M:N sparsity and
further work shows how to quantify performance loss, prune models, and use M:N
sparsity during training (Zhou et al., 2021; Sun et al., 2021). MMM makes no
assumptions about fine-grained sparsity.

The Tensor Algebra compiler (TACO) (Kjolstad et al., 2017) produces optimized
computation kernels for general tensor operations on many sparse data formats.
Techniques have been extended to GPUs (Senanayake et al., 2020), distributed
memory and heterogeneous architectures (Yadav et al., 2022), and dynamic data
structures that support updates (Chou and Amarasinghe, 2022). Compiled kernels
match or approach the performance of hand-optimized code, easing development
effort and making code robust to changing data and computing environments.
Tensor compilation could be extended to use just-in-time compilation or
dynamic code evaluation to adapt to sparsity at runtime. However, all
compilation techniques target a sparse data structure. MMM targets mid-ranges
of sparsity in which there are not enough zeros to justify sparse data
structures as demonstrated by our performance against Intel’s sparse-sparse
and sparse-dense methods. Also, compiled kernels are arbitrarily complex and
MMM supports only GeMM.

The principles of tensor compilation and sparse data structures have been
extended to the design of sparse accelerators. Dual-Side Sparse Tensor Core
(Wang et al., 2021) uses bitmap representations to perform outer-product
SpGEMM in hardware. It captures sparsity in both activations and weights.
SparseLoop (Wu et al., 2021) builds a cost-driven optimizer that generalizes
the principles in DSTC and NVidia 2:4. It generates designs for multiple
sparse representations and compression formats. An abstract dataflow machine
model for sparse tensor algebra leads to a compiler and design of
reconfigurable and fixed-function accelerators (Hsu et al., 2023).

##  6\. Conclusion and Discussion

MMM shows the potential to exploit emergent sparsity for matrix multiplication
on transformers. It demonstrates a reduction of the work needed to perform
matrix multiplication on sparse matrices.

Future work involves involves improvements to more classes of matrices such as
larger or denser, as well as extend it to work on different architectures. The
low-level code for multiplying small regions could be improved to increase the
register reuse and could outperform MKL dense for a larger range of
sparsities. Other algorithmic style optimizations could be used such as using
a JIT to compile the exact forms needed for the existing patterns, instead of
pre-compiling all possible forms. The code could also be extended to
seamlessly choose the optimal base function for each recursive block between
both MMM approaches as well as different standard approaches from a vendor
library.

We plan to extend this work to GPUs. This will involve some changes to the
required sparsity structure. For example, to maximize the performance and
align the different warps, all of the threads in a warp need to compute on the
same pattern of rows of B and elements from A. This would require vertical
blocks in A in addition to horizontal blocks in B. The $A$ matrices of MLP
computation were column matrices and implement this pattern. We believe this
approach will extend to GPUs.

## References

  * (1)
  * Buluç and Gilbert (2012) Aydin Buluç and John R Gilbert. 2012.  Parallel sparse matrix-matrix multiplication and indexing: Implementation and experiments.  _SIAM Journal on Scientific Computing_ 34, 4 (2012), C170–C191.
  * Chou and Amarasinghe (2022) Stephen Chou and Saman Amarasinghe. 2022.  Compilation of dynamic sparse tensor algebra.  _Proceedings of the ACM on Programming Languages_ 6, OOPSLA, Article 175 (2022).
  * Corporation (2022) Intel® Corporation. 2022\.  _Optimize Inference with Intel® CPU Technology_ (1st ed.).  Intel® Corporation, Santa Clara, CA.  Available at https://www.intel.com/content/dam/www/central-libraries/us/en/documents/2022-06/optimize-inference-cpu-technology-sb.pdf.
  * Fedus et al. (2021) William Fedus, Barret Zoph, and Noam Shazeer. 2021\.  Switch Transformers: Scaling to Trillion Parameter Models with Simple and Efficient Sparsity.  _Journal of Machine Learning Research_ 23 (2021).
  * Gustavson (1978) Fred G Gustavson. 1978\.  Two fast algorithms for sparse matrices: Multiplication and permuted transposition.  _ACM Transactions on Mathematical Software (TOMS)_ 4, 3 (1978), 250–269.
  * Hao et al. (2023) Jianwei Hao, Piyush Subedi, Lakshmish Ramaswamy, and In Kee Kim. 2023.  Reaching for the Sky: Maximizing Deep Learning Inference Throughput on Edge Devices with AI Multi-Tenancy.  _ACM Transactions on Internet Technology_ 23, 1 (2023), 1–33.
  * Hsu et al. (2023) Olivia Hsu, Maxwell Strange, Ritvik Sharma, Jaeyeon Won, Kunle Olukotun, Joel S. Emer, Mark A. Horowitz, and Fredrik Kjølstad. 2023.  The Sparse Abstract Machine. In _Proceedings of the ACM International Conference on Architectural Support for Programming Languages and Operating Systems (ASPLOS)_.
  * Intel Corporation (2010) Intel Corporation. 2010\.  _Intel Cilk Plus Language Specification_.  Document Number: 324396-001US. Available from http://software.intel.com/sites/products/cilk-plus/cilk_plus_language_specification.pdf.
  * Jaszczur et al. (2021) Sebastian Jaszczur, Aakanksha Chowdhery, Afroz Mohiuddin, Łukasz Kaiser, Wojciech Gajewski, Henryk Michalewski, and Jonni Kanerva. 2021\.  Sparse is Enough in Scaling Transformers.
  * Kitaev et al. (2020) Nikita Kitaev, Lukasz Kaiser, and Anselm Levskaya. 2020\.  Reformer: The Efficient Transformer. In _International Conference on Learning Representations (ICLR)_.
  * Kjolstad et al. (2017) Fredrik Kjolstad, Shoaib Kamil, Stephen Chou, David Lugato, and Saman Amarasinghe. 2017.  The tensor algebra compiler.  _Proceedings of the ACM on Programming Languages_ 1, OOPSLA (2017), 1–29.
  * Lattner (2002) Chris Lattner. 2002\.  _LLVM: An Infrastructure for Multi-Stage Optimization_.  Master’s thesis. Computer Science Dept., University of Illinois at Urbana-Champaign, Urbana, IL.  See http://llvm.cs.uiuc.edu..
  * Lattner and Adve (2004) Chris Lattner and Vikram Adve. 2004.  LLVM: A compilation framework for lifelong program analysis & transformation. In _International symposium on code generation and optimization, 2004. CGO 2004._ IEEE, 75–86.
  * Li et al. (2023) Zonglin Li, Chong You, Srinadh Bhojanapalli, Daliang Li, Ankit Singh Rawat, Sashank Reddi, Ke Ye, Felix Chern, Felix Yu, Ruiqi Guo, et al. 2023\.  The lazy neuron phenomenon: On emergence of activation sparsity in transformers. In _Conference on Parsimony and Learning (Recent Spotlight Track)_.
  * Liu et al. (2023) Zichang Liu, Jue Wang, Tri Dao, Tianyi Zhou, Binhang Yuan, Zhao Song, Anshumali Shrivastava, Ce Zhang, Yuandong Tian, Christopher Re, et al. 2023\.  Deja vu: Contextual sparsity for efficient llms at inference time. In _International Conference on Machine Learning_. PMLR, 22137–22176.
  * Mishra et al. (2021) Asit Mishra, Jorge Albericio Latorre, Jeff Pool, Darko Stosic, Dusan Stosic, Ganesh Venkatesh, Chong Yu, and Paulius Micikevicius. 2021.  Accelerating Sparse Deep Neural Networks.  arXiv:2104.08378 [cs.LG]
  * Morgan (2023) Timothy Prickett Morgan. 2023\.  Why AI Inference Will Remain Largely On The CPU.  Morgan, T. P. (2023, April 5). https://www.nextplatform.com/2023/04/05/why-ai-inference-will-remain-largely-on-the-cpu/.  [Accessed 08-02-2024].
  * Schardl and Lee (2023) Tao B. Schardl and I-Ting Angelina Lee. 2023.  OpenCilk: A Modular and Extensible Software Infrastructure for Fast Task-Parallel Code _(PPoPP ’23)_. 189––203.  https://doi.org/10.1145/3572848.3577509
  * Schardl et al. (2019) Tao B Schardl, William S Moses, and Charles E Leiserson. 2019\.  Tapir: Embedding recursive fork-join parallelism into LLVM’s intermediate representation.  _ACM Transactions on Parallel Computing (TOPC)_ 6, 4 (2019), 1–33.
  * Senanayake et al. (2020) Ryan Senanayake, Changwan Hong, Ziheng Wang, Amalee Wilson, Stephen Chou, Shoaib Kamil, Saman Amarasinghe, and Fredrik Kjolstad. 2020.  A sparse iteration space transformation framework for sparse tensor algebra.  _Proceedings of the ACM on Programming Languages_ 4, OOPSLA, Article 158 (2020).
  * Strassen (1969) V. Strassen. 1969\.  Gaussian elimination is not optimal.  _Numer. Math._ 13, 4 (1969).
  * Sun et al. (2021) Wei Sun, Aojun Zhou, Sander Stuijk, Rob G. J. Wijnhoven, Andrew Nelson, Hongsheng Li, and Henk Corporaal. 2021.  DominoSearch: Find layer-wise fine-grained N: M sparse schemes from dense neural networks. In _Neural Information Processing Systems_.
  * Wang et al. (2014) Endong Wang, Qing Zhang, Bo Shen, Guangyong Zhang, Xiaowei Lu, Qing Wu, and Yajuan Wang. 2014.  _Intel Math Kernel Library_.  Springer International Publishing, Cham, 167–188.  https://doi.org/10.1007/978-3-319-06486-4_7
  * Wang et al. (2021) Yang Wang, Chen Zhang, Zhiqiang Xie, Cong Guo, Yunxin Liu, and Jingwen Leng. 2021\.  Dual-side Sparse Tensor Core. In _International Symposium on Computer Architecture (ISCA)_.
  * Williams et al. (2023) V. V. Williams, Y. Xu, X. Xu, and R. Zhou. 2023\.  New Bounds for Matrix Multiplication: from Alpha to Omega.  _arXiv:2307.07970_ (2023).
  * Wu et al. (2021) Yannan Nellie Wu, Po-An Tsai, Angshuman Parashar, Vivienne Sze, and Joel S. Emer. 2021.  Sparseloop: An Analytical, Energy-Focused Design Space Exploration Methodology for Sparse Tensor Accelerators. In _2021 IEEE International Symposium on Performance Analysis of Systems and Software (ISPASS)_.
  * Yadav et al. (2022) Rohan Yadav, Alex Aiken, and Fredrik Kjolstad. 2022.  DISTAL: the distributed tensor algebra compiler. In _Proceedings of the 43rd ACM SIGPLAN International Conference on Programming Language Design and Implementation (PLDI)_.
  * Zhang et al. (2022) Susan Zhang, Stephen Roller, Naman Goyal, et al. 2022\.  OPT: Open Pre-trained Transformer Language Models.  arXiv:2205.01068 [cs.CL]
  * Zhou et al. (2021) Aojun Zhou, Yukun Ma, Junnan Zhu, Jianbo Liu, Zhijie Zhang, Kun Yuan, Wenxiu Sun, and Hongsheng Li. 2021\.  Learning N:M Fine-grained Structured Sparse Neural Networks From Scratch.


Presently at ]Research Division Neutrons and Muons, Paul Scherrer Institute,
CH-5232 Villigen, Switzerland

# Magnetic correlations in triangular antiferromagnet FeGa2S4

K. Guratinder  Laboratory for Neutron Scattering and Imaging, Paul Scherrer
Institute, CH-5232 Villigen, Switzerland  Department of Quantum Matter
Physics, University of Geneva, CH-1211 Geneva, Switzerland  M. Schmidt  Max-
Planck-Institut für Chemische Physik fester Stoffe, 01187 Dresden, Germany  H.
C. Walker, R. Bewley  Rutherford Appleton Laboratory, ISIS Facility, Chilton,
Didcot, Oxon OX11 0QX, United Kingdom  M. Wörle  Laboratorium für Anorganische
Chemie, CH-8093 Zürich, Switzerland  D. Cabra  Instituto de Física de Líquidos
y Sistemas Biológicos, CCT La Plata, CONICET and Departamento de Física,
Facultad de Ciencias Exactas, Universidad Nacional de La Plata, C.C. 67, 1900
La Plata, Argentina  S. A. Osorio  Instituto de Física de Líquidos y Sistemas
Biológicos, CCT La Plata, CONICET and Departamento de Física, Facultad de
Ciencias Exactas, Universidad Nacional de La Plata, C.C. 67, 1900 La Plata,
Argentina  M. Villalba  Instituto de Física de Líquidos y Sistemas Biológicos,
CCT La Plata, CONICET and Departamento de Física, Facultad de Ciencias
Exactas, Universidad Nacional de La Plata, C.C. 67, 1900 La Plata, Argentina
A. K. Madsen  Quantum Criticality and Dynamics Group, Paul Scherrer Institute,
CH-5232 Villigen, Switzerland  L. Keller  Laboratory for Neutron Scattering
and Imaging, Paul Scherrer Institute, CH-5232 Villigen , Switzerland  A.
Wildes  Institut Laue-Langevin, 156X, 38042 Grenoble Cédex, France  P. Puphal
Laboratory for Multiscale Materials Experiments, Paul Scherrer Institute,
CH-5232 Villigen, Switzerland  A. Cervellino  Swiss Light Source, Paul
Scherrer Institute, CH-5232 Villigen, Switzerland  Ch. Rüegg  Laboratory for
Neutron Scattering and Imaging, Paul Scherrer Institute, CH-5232 Villigen,
Switzerland  Department of Quantum Matter Physics, University of Geneva,
CH-1211 Geneva, Switzerland  [  O. Zaharko<EMAIL_ADDRESS>Laboratory
for Neutron Scattering and Imaging, Paul Scherrer Institute, CH-5232 Villigen,
Switzerland.

###### Abstract

The crystal structure and magnetic correlations in triangular antiferromagnet
FeGa2S4 are studied by x-ray diffraction, magnetic susceptibility, neutron
diffraction and neutron inelastic scattering. We report significant mixing at
the cation sites and disentangle magnetic properties dominated by major and
minor magnetic sites. The magnetic short-range correlations at 0.77 Å-1
correspond to the major sites and being static at base temperature they evolve
into dynamic correlations around 30 - 50 K. The minor sites contribute to the
magnetic peak at 0.6 Å-1, which vanishes at 5.5 K.
Our analytical studies of triangular lattice models with bilinear and
biquadratic terms provide the ratios between exchanges for the proposed
ordering vectors. The modelling of the inelastic neutron spectrum within
linear spin wave theory results in the set of exchange couplings $J_{1}=1.7$
meV, $J_{2}=0.9$ meV, $J_{3}=0.8$ meV for the bilinear Heisenberg Hamiltonian.
However, not all features of the excitation spectrum are explained with this
model.

neutron scattering, frustrated magnetism, triangular lattice

##  I Introduction

The interest in the MT2X4 family of compounds (M = Ni, Fe; T = Ga, Al; X = S,
Se) arises due to the quasi- 2-dimensional (2D) triangular geometry of the
magnetic M-sublattice. This family, especially NiGa2S4, are strong candidates
as an experimental realisation of the triangular lattice antiferromagnet
(TLAFM), the extensively studied theoretical model in frustrated magnetism.
For the spin S=1/2 TLAFM case a resonating valence bond spin liquid ground
state was proposed by AndersonAnderson1973 , which was later shown instead to
have a classical-like 120 $\deg$ planar ground stateJolicoeur1990  in spite of
the quantum fluctuations that only reduce the value of the order parameter.
For Ising spins, on the other hand, it has been shown that the first nearest-
neighbor ($J_{1}$) Ising TLAFM model stays disordered even down to zero
temperature. Additional exchange interactions spanning over the second
($J_{2}$) and third ($J_{3}$) nearest-neighbors could enhance frustration and
for certain combinations of the exchange parameters exotic ground states are
predicted. For example for $J_{1}/J_{3}$ = - $\frac{1}{3}$ or $J_{1}/J_{2}$ =
- $\frac{1}{4}$ a skyrmion crystal multi-$k$ state should be realisedOkubo2012
. Such swirling configurations with topological properties are of interest for
the fast emerging field of spintronics.
MGa2S4 compounds with M = Ni, Fe crystallise in the trigonal space group
$P\bar{3}m1$Dogguy1980 . All atoms are located on triangular layers, which
generally could be stacked in the A-, B- or C-fashion (Fig. 1a). In MGa2S4 the
sulphur layers stack in the BCBC sequence providing octahedral and tetrahedral
voids. The octahedral voids are occupied by the smaller M-cations, and the
tetrahedral voids by Ga (Fig. 1b). The unit cell of MGa2S4 is terminated by an
empty void layer. Our x-ray diffraction study detected significant mixed
occupancies (inversion) of the cation sites which we associate with highly
nontrivial magnetic behaviour of the MGa2S4 series. This feature has not been
observed before and it is one of the important outcomes of our work.

Figure 1: a) Projection of the A-, B-, C-layers on the $ab$-plane. b) View on
the MGa2S4 unit cell with labeled (from left to right) sequence of S, Ga and M
atoms, coordination of the voids between the S-layers and order of the A, B, C
layers.

The published bulk propertiesNakatsuji2007 ; Myoung2008 ; Zhao2012  suggest
stronger exchange couplings for FeGa2S4 (Fe2+, S=2) compared to NiGa2S4 (Ni2+,
S=1), but the main observations remain very similar. The Curie-Weiss
temperature is $\theta_{CW}$=160 K for FeGa2S4 ($\theta_{CW}$=80 K for the Ni-
analogue) and there is a bifurcation between field-cooled (FC) and zero-field-
cooled (ZFC) dc magnetic susceptibilities at Tf(Fe)=16 K (Tf(Ni)=8-9 K). No
long-range order has been detected so far for both systems, but an inelastic
neutron scattering (INS) studyStock2010  for NiGa2S4 reported incommensurate
(ICM) quasielastic short-range magnetic correlations in the vicinity of
Q=(0.155 0.155 0) and spin waves propagating from this Q-point up to 3.5 meV.
The dynamic muon relaxation rates Reotier2012 ; Zhao2012  show an anomaly at
T∗(Fe)=33 K (T∗(Ni)=10 K) approaching correlation times ${\tau_{c}}$(Fe)=10-6
s (${\tau_{c}}$(Ni)=10-7 s). These processes are much slower than the time
scale of exchange interaction (10-12 s) and the corresponding states are
interpreted as a slowly fluctuating ’spin gel’. Surprisingly, no anomaly at Tf
was observed by 57Fe Mössbauer, $\mu$SR and specific heat CM, while the T∗
anomaly is not present in magnetic bulk propertiesNakatsuji2007 ; Myoung2008 .
The magnetic component of the specific heat CM(T) has a double-peak structure
for both compounds with a first maximum at 10 K for both and a second maximum
at 60 K (100 K) for FeGa2S4Nakatsuji2007  (NiGa2S4Nakatsuji2005 ). CM(T) has a
T2 dependence at low temperatures, which indicates gapless and linearly
dispersive modes in 2D, and it is insensitive to magnetic fields up to 7 T.
Our samples show the same magnetic behaviour and we complement experimental
information on the FeGa2S4 system by neutron experiments, both diffraction and
inelastic scattering.
From the theoretical side a number of extensions to the TLAFM model were
proposed to explain the observed magnetic behaviour of the Ni- and Fe-
analogues. A strong, dominant third-nearest-neighbor antiferromagnetic
exchange interaction is supported by DFT calculations Mazin2007  and could
explain the ICM magnetic propagation vector of NiGa2S4Stock2010 . Biquadratic
interactions, which would yield spin-nematic quadrupolar correlations, account
for the double-peak structure in specific heatTsunetsugu2006 ; Lauchli2006 ;
Bhattacharjee2006 . The origin of the spin nematic could be the coupling of
magnetic correlations to phonons observed by Raman scatteringValentine2020 .
The $Z_{2}$ vortex model with binding-unbinding transitionTsunetsugu2007
could also explain the high temperature anomalies in magnetic properties,
while linearly dispersing Halperin-Saslow modes in the absence of long-range
magnetic order could elucidate the anomalous low-temperature
propertiesPodolsky2009 . We scrutinize the TLAFM model with $J_{1}$, $J_{2}$,
$J_{3}$ billinear exchanges and show that it is not sufficient to model our
INS data. We associate the unconventional magnetic properties with the non-
ideal structure of the MGa2S4 family.

##  II Crystal structure of FeGa2S4 and NiGa2S4

###  II.1 Study of structure imperfection by X-ray diffraction

First we check thoroughly the crystal structure of our polycrystalline FeGa2S4
samplenote1  analysing synchrotron x-ray diffraction powder pattern collected
at 300 K. These data were well refined within the crystal structure published
by Dogguy SimiriDogguy1980 : the Ga and two S atoms located at the 2(d) sites
and Fe atoms at the 1(b) site of the $P\bar{3}m1$ space group (RBragg=3.64 %,
Rf=1.99 %, ${\chi^{2}}$=47.3). The fits became even better (RBragg=2.98 %,
Rf=1.89 %, ${\chi^{2}}$=46.2) when the inversion between the Fe and Ga atoms
was introduced with the total occupancy constrained to the stoichiometric
composition FeGa2S4. According to this refinement the 1(b) site is occupied
20% by gallium and 80% by iron, whereas the 2(d) site - by 10% iron and 90%
gallium, respectively. As the model with partial inversion is only nominally
better and a small amount of Ga2O3 impurity phase was detected in the powder
pattern, we verified this result by x-ray diffraction on single crystals.

Figure 2: Top: The hk0 (left) and h0l (right) reciprocal layers reconstruction
based on the X-ray single crystal diffraction of FeGa2S4. Bottom: (left)
scanning electron microscopy image of a single crystal of FeGa2S4 and (right)
zoom into its left marked corner.

For several studied crystals (even very small ones) the shape of reflections
was not exceptionally good - they had ’tails’ typical for small arbitrary
rotations of macroscopic domains around the [001] axis, as can be seen in the
reciprocal hk0-layer reconstruction presented in Fig. 2 (top). Nevertheless,
the consequent data reduction lead to good data with Rint=3.69 % for 520
unique reflections for the Laue group $\bar{3}m1$. A refinement assuming the
non-disordered FeGa2S4 structure with full occupancies converged at acceptable
R-values: R1 = 2.83 %, wR2 = 8.21 %, GooF = 1.184 for 514 $F_{o}>$ 4 $\sigma
F_{o}$. The difference Fourier maps, as in the case of powder pattern,
contained relatively strong residual electron density peaks: positive at the
1(b) Fe site (corresponding to around 3 e/Å3) and negative (corresponding to
around -1.7 e/Å3) in the vicinity of the 2(d) Ga sites. Along with this no
significant electron density occurred between the atomic positions of the
ideal structure. The residual electron density could be accounted for by the
partial inversion at the 1(b) and 2(d) sites. Refinement converged to 21.5(3)%
Ga (and 78.4(3)% Fe) on the 1(b) Wyckoff position and 10.8(1)% Fe (and
89.2(1)% Ga) on the Wyckoff position 2(d) with improved figure of merit: R1 =
2.04 %, Rw2 = 6.08 %, GooF = 1.367.
We took the indication of partial inversion from x-ray diffraction with great
caution. Bulk properties (magnetic susceptibility, specific heat) measured by
us and by other groups are very similar suggesting that partial inversion
might improve the understanding of magnetic properties of the whole MGa2S4
family. To clarify this, we synthesised and characterised isostructural
NiGa2S4. Refinement of NiGa2S4 single crystal x-ray diffraction data revealed
the same degree of inversion as for FeGa2S4. This complies with nuclear
quadrupole resonance (NQR) studies of polycrystalline NiGa2S4 samples
Takeya2008 ; Nambu2009 , which reported two Ga signals. Our results imply that
these NQR signals could arise from gallium atoms located on two different
crystallographic sites. We also considered whether other imperfections
discussed in the literature could explain our diffraction data. For example,
high-resolution electron microscopy detected the presence of the stacking
faults along the [001]-directionNambu2009 . Any significant fraction of
stacking faults generated by a shift of the layers, as realised for example in
the case of the rhombohedral structure of MgAl2S4 or by introducing layers
found in Fe2Ga2S5 or in GaS with the hexagonal $\beta$-InSe structure, would
lead to residual electron density off the atomic positions of the ideal
structure. The same holds for the stacking faults along the slip planes
(11-20), (10-11) and (11-23), which sometimes occur in the hexagonal close
packing structures. Diffuse scattering should also be observed, but this is
not the case in our diffraction patterns presented in Fig. 2, top. We also
checked for twinning, finding no indication of twofold rotations around [110],
[1-10] and [001], the usual twin laws occurring in trigonal systems. Our
scanning electron microscopy images clearly show that crystals consist of
layers with a thickness on the micrometer scale, which are rotated against
each other by small angles (Fig. 2 bottom). One more sample imperfection was
reported by Nambu et al. Nambu2009  based on energy-dispersive X-ray
spectroscopy (EDXS) - nonstoichiometry in the sulfur content. Our EDXS
analysis show no sulfur deficiency within the detection limit.
Thus we conclude that polycrystalline and single crystal MGa2S4 samples have
significant inversion on the cation sites. In the following text we correlate
this partial inversion with magnetic bulk properties and neutron scattering
results, labeling the 1(b) site as the ’major’-site and the 2(d) site as the
’minor’-site, since they are occupied by 80% and 10% with the magnetic ion,
respectively.

###  II.2 Temperature evolution of crystal structure

To identify possible changes to the crystal structure of FeGa2S4 we measured
and refined synchrotron x-ray powder diffraction patterns of a polycrystalline
sample in the 5 K- 300 K temperature range. The obtained temperature evolution
of the lattice constants, atomic positions and thermal parameters is presented
in (Fig. 3). No significant structural anomalies could be detected in these
X-ray patterns. The lattice constants show a typical Debye thermal expansion,
the slope of the ratio c/a changes below 50 K and near 10 K. The interatomic
distances dFeS2, dGaS1 decrease, while dGaS2 increases with lower temperature.
The anisotropic thermal parameters of all atoms decrease very much alike, only
the thermal parameters $\beta_{11}$ and $\beta_{3}$ of the 2d site decrease
faster and slower, respectively.
These parameters of FeGa2S4 do not point to anomalies in vibrations of the
sulfur atom S2 detected in the sister compound NiGa2S4 by Raman
scatteringValentine2020 .

Figure 3: Temperature evolution of the lattice constants (a), selected bond
distances (b) and thermal parameters (c, d) of FeGa2S4 from synchrotron x-ray
powder diffraction.

##  III Magnetic properties

###  III.1 Magnetic bulk properties

Figure 4: a) Magnetic susceptibility $\chi$ measured at 0.1 T on a single
crystal of FeGa2S4 (mass m=2.87 mg) along the three principal directions
$[100]$, $[-120]$ and $[001]$. Zero field cooled (ZFC) and field cooled (FC)
$\chi$ are plotted by empty and filled circles, respectively. Inset: zoom of
the low temperature region, b) derivative of ZFC magnetic susceptibility
highlighting deflection points. c) Inverse susceptibility ${1/\chi}$ for
H$//[100]$ (blue) and H$//[001]$ (red) and fit to the Curie-Weiss law in
yellow. d) - the ratio of ZFC $\chi_{ab}/\chi_{c}$.

We briefly report magnetic bulk properties of our FeGa2S4 samples. Fig. 4
(left) presents the magnetic susceptibility of a FeGa2S4 single crystal
measured at 0.1 T in the range 1.8 K - 350 K for the three principal
orthogonal crystallographic directions: $[001]$, $[100]$ and $[-120]$. At high
temperatures (200 - 350 K) the curves almost overlap and fits to the Curie-
Weiss law yield antiferromagnetic Curie-Weiss temperatures $\Theta_{\rm CW}$=
-165(11) K for the $[001]$ direction and -151(9) K for $[100]$ and $[-120]$
and effective magnetic moments $m_{eff}$=5.53 $\mu_{B}$ and 5.08 $\mu_{B}$,
respectively. The susceptibility in the $ab$-plane and along the $c$-axis
deviate below 200 K and this anisotropy reaches 20% near 16 K. At Tf=16 K the
susceptibility has a kink and a bifurcation between the field-cooled (FC) and
zero-field- cooled (ZFC) sample for all measured directions. The difference
between the FC and ZFC curves is tiny, 0.004 $m_{B}$/Fe. Careful inspection of
the derivative of ${\chi}$ (Fig. 4, right inset) reveals an additional change
of the slope at 5.6 K. We associate it with long-range ordering of the minor
Fe-spins by powder neutron diffraction (see Section IV).
For our polycrystalline FeGa2S4 sample the Curie-Weiss temperature is
-158.0(9) K. This results in the effective exchange constant

$J/{k_{B}}=-\frac{3\Theta_{\rm CW}}{zS(S+1)}\sim 13~{}K\sim 1.1~{}meV,$

here $k_{B}$ is the Boltzmann constant, $z=6$ is the coordination number and
$S=2$ is the spin of the Fe2+-ions.

###  III.2 Magnetic static correlations from neutron diffraction

Figure 5: Neutron diffraction DMC patterns of FeGa2S4 at 1.3 K and 50 K (a).
Temperature dependence of the magnetic peak at Q1=0.6 Å-1 (b).

Figure 6: Magnetic diffuse scattering measured at 1.5 K (red) and 60 K (cyan)
on D7 (Ei=3.8 meV) and fit to the Warren function a) and calculated by the RMC
method (b) with resulting radial spin-correlation function $<S(0)\cdot S(r)>$
shown in insert.

It is crucial to understand, if the magnetic properties of FeGa2S4 samples
could be disentangled into major- and minor- sites properties. First we
address this question with powder neutron diffraction (PND).
In PND patterns collected on DMC two different features of magnetic origin
appear. Prominent in Fig. 5a is the peak at Q1=0.6 Å-1 below TN=5.5 K. It
could be indexed with the incommensurate propagation vector
$\bf{k}$=(0.1737(1) 0.1737(1) 0) in analogy to $\bf{k}$=(0.153 0.153 0)
reported for NiGa2S4. However, no long-range magnetic arrangement akin to only
the major site reproduces the diffraction pattern. Models with Fe only on the
1(b) site require the second peak at (0.1737 0.1737 1) with Q=0.79 Å-1 to have
significant intensity. The single (0.1737 0.1737 0) peak can be reproduced by
a model with the Fe moments on the 2(d) site. The moments are sinusoidally
modulated along the $z$-axis. The moment amplitude 7.0(2) $\mu_{B}$/Fe
(maximal value of the moment) would be however unreasonably large for the 10 %
occupied 2(d) site. So the Mz component on the minor 2(d) site is most
probably significant but not the only ’participant’ in the magnetic ordering
at TN=5.5 K and Q1=0.6 Å-1.
The second, less evident feature in PND is diffuse magnetic scattering
centered at Q2=0.77 Å-1 (Fig. 5a). The XYZ-polarization analysis of powder
neutron diffraction patterns collected on D7 verifies that this is purely
magnetic scattering (Fig. 6). In NiGa2S4 a similar diffuse bump was reported
beneath the $\bf{k}$=(0.153 0.153 0) magnetic peak (Fig. 4 in Ref.
$[$Nakatsuji2005, $]$) and it was interpreted in terms of 2D correlations
extending within 25(3) Å. In the case of FeGa2S4 the diffuse feature is
broader. To estimate the 2D correlation length we first used the Warren
functionWarren1941  initially developed to describe scattered intensity from
parallel, equidistant layers randomly rotated along the normal direction.
Modelling of the magnetic PND pattern presented in Fig. 6 (a) suggests an in-
plane correlation length of 10 Å, which does not exceed the 4th nearest
neighbor Fe-Fe distance. As an alternative route, we used the reverse Monte-
Carlo (RMC) algorithm implemented in the SPINVERT programPaddison2013 . We fit
the magnetic PND pattern to a large array (30x30x10 unit cells) of spin
vectors located on the major 1(b) sites (Fig. 6 b). The spin-pair correlation
function in real space $<\bf{S}$(0)$\cdot$ $\bf{S}$(r)$>$ is calculated from
the fitted spin configurations by the program SPINCORRELPaddison2013 . The
negative sign of correlations for the in-plane distances d${}_{1}^{1b}$=3.66
Å, d${}_{2}^{1b}$=6.35 Å, d${}_{3}^{1b}$=7.33 Å (three first points in Fig. 6b
insert) suggests antiparallel alignment of spins within a triangular layer.
RMC runs with additional spins on the minor 2(d) sites did not improve the
fit, so the minor spins are not decisive for the Q2 feature.
We summarise our PND results as follows: long-range correlations resulting in
the Q1=0.6 Å-1 peak arise mostly from the minor Fe-spins, while the short-
range correlations at Q2=0.77 Å-1 are magnetic in-plane correlations of the
major Fe-spins. To distinguish whether these are these static or dynamic we
turned to inelastic neutron scattering, presented below.

###  III.3 Magnetic dynamic correlations from inelastic neutron scattering

The excitation spectrum of polycrystalline FeGa2S4 was measured on the LET
spectrometerBewleyLET . Fig. 7 displays the wide-ranging S(Q, $\omega$)
accessed with the $E_{i}$= 24.3 meV rep-mode at 2 K and 180 K. Two prominent
INS features are evident.

Figure 7: Top: Excitation spectrum of FeGa2S4 measured on LET with $E_{i}$ =
24.3 meV at the lowest T= 2 K (a) and the highest 180 K (b) temperatures.
Bottom: Temperature evolution of the low-Q feature with signal integrated over
Q=[0.3$\ldots$1.5] Å-1 with the $E_{i}$= 24.3 meV rep-mode ($\delta E_{0meV}$=
1.15 meV) (c) and over Q=[0.4$\ldots$0.8] Å-1, $E_{i}$ = 7.7 meV ($\delta
E_{0meV}$ = 0.22 meV) (d). The Al can was subtracted and data were corrected
for the Bose factor.

The first one is peaked at Q2 = 0.77 Å-1 and E = 5 meV. It decreases towards
the elastic line, suggesting that these might be gapped spin waves, but it is
hard to identify the value of the anticipated gap and temperature of its
softening even from the lower energy $E_{i}$= 7.7 meV rep-mode spectrum (Fig.
7 d). The second, high-energy feature at 13 meV spans over the whole measured
interval of 1 $<$ Q $<$ 6 Å-1. The low Q-part (Q $<$ 4.5 Å-1) reduces with
temperature, thus it has magnetic origin. The high Q part (Q $>$ Q=4.5 Å-1)
persists, so it has significant phonon contribution. These two features reside
on a very broad background which remains even at 180 K, the highest
temperature reached in the experiment.
The anticipated gap region for the first feature is magnified in the top panel
of Fig. 8, where S(Q, $\omega$) maps measured with the $E_{i}$ = 3.8 meV rep-
mode at three representative temperatures 2 K, 38 K and 60 K are presented. At
2 K the correlations are predominantly static. In Fig. 8d the counts summed
within 0 $<$ E${}_{1}<$ 0.2 meV are divided by a factor 50 to match the scale
of the intensity within the gap 0.2 $<$ E${}_{2}<$ 0.5 meV and of the
excitations above the gap 0.5 $<$ E${}_{2}<$ 3 meV. Interestingly, at 2 K a
weak dispersive feature starting from Q1 could also be identified. At higher
temperatures it is completely obscured by a continuum centered at Q2=0.77 Å-1.
The continuum gains intensity in the 30-40 K region and is reduced with
further increasing temperature.
It would be of interest to compare the temperature evolution of various static
and dynamic correlations observed on LET and in bulk magnetic measurements.
Unfortunately the insufficient resolution does not allow to separate the Q1
feature from the nearby nuclear peak Fig. 8d. Therefore we present the
temperature dependence of the elastic intensity of the predominantly Q2
feature summed within 0.6 $<$ Q $<$ 1 Å-1 in Fig. 8e. It decreases sharply
around 50 K. The two prominent INS signals centered at 5 and 13 meV, in
contrast are maximal at 50 K. The 5 meV feature is integrated within the
$E_{i}$= 3.8 meV setup and this allows the documentation of the deflection
point at 10 K (close to Tf).
The contribution of spin waves from Q1=0.6 Å-1 is minute in the collected INS
signal and cannot be disentangled from the broad major feature centered at Q2
= 0.77 Å-1. The origin of the Q2 feature evolves from static into dynamic with
temperature increase and reaches maximum intensity at 50 K, which coincides
with the broad maximum in the magnetic part of the specific heat.Nakatsuji2007
We presume it originates from the major spins.

Figure 8: Top: Part of the excitation spectrum of FeGa2S4 measured on LET with
$E_{i}$= 3.8 meV ($\delta E_{0meV}$= 0.1 meV) at T= 2 K (a), 38 K (b), 60 K
(c). Bottom: (d) 2 K elastic signal summed over E0= [0$\ldots$0.2] meV (green
symbols, signal divided by 50) and inelastic signals E1= [0.2$\ldots$0.5] meV
(blue symbols), E2= [0.5$\ldots$3] meV (red symbols). (e) Temperature
evolution of the elastic (blue squares) and inelastic (red and black circles)
features. The elastic signal is integrated over Q=[0.6$\ldots$1] Å-1 and
E0=[-0.1$\ldots$0.1] meV, the inelastic signals are corrected for the Bose
factor and integrated over Q=[0.44$\ldots$2] Å-1, E1= [0.5$\ldots$2] meV and
Q=[0.3$\ldots$1.5] Å-1, E2= [11$\ldots$13] meV. The last curve is extracted
from the $E_{i}$= 24.3 meV rep-mode.

In brief: The low energy part of INS data is dominated by the (5 meV, 0.77
Å-1) feature, the high energy part has a feature at 13 meV which is magnetic
for Q $<$ 4.5 Å-1 and has a phonon contribution for Q $>$ 4.5 Å-1. The
magnetic INS signal has maximum intensity at 50 K, where the elastic short-
range correlations diminish. The broad scattering background beneath these
features persist even at 180 K.

##  IV Modelling of magnetic Hamiltonian

We used several theoretical approaches to explain these variegated
experimental observations. We start with the ideal triangular lattice and
present analytical results for models with increasing number of bilinear
exchanges and with the biquadratic exchange (Section. IV.1). Then we test the
model with three bilinear exchanges against our INS data (Section. IV.2).

###  IV.1 Analytical results for models with bilinear and biquadratic
exchanges

First, we try to establish the range of exchange parameters which could give
rise to magnetic correlations observed at Q1=0.6 Å-1 for FeGa2S4 and 0.57 Å-1
for NiGa2S4Nakatsuji2005  within the triangular lattice Heisenberg bilinear
model. In this Section we use the theoretical notation (TN) for the reciprocal
space (see Appendix VII for the definition of TN and its relation to the
experimental notation (EN)).
We inspect analytically, with rising complexity, the single $J_{i=1,3}$, the
bilinear $J_{i=1,3}-J_{j=1,3}$ and the trilinear $J_{1}-J_{2}-J_{3}$ models on
the two-dimensional triangular lattice described by the Hamiltonian:

$H=J_{1}\sum_{<i,j>}({\bf{S}}_{i}\cdot{\bf{S}}_{j})+J_{2}\sum_{<<i,j>>}({\bf{S}}_{i}\cdot{\bf{S}}_{j})+J_{3}\sum_{<<<i,j>>>}({\bf{S}}_{i}\cdot{\bf{S}}_{j})$
(1)

Since we are interested in the region of the $J_{1}$, $J_{2}$, $J_{3}$ phase
diagram with the incommensurately modulated ground states, we express the spin
$\mathbf{S}_{i}$ on site $\mathbf{r}_{i}$ asTamura2008 :

$\mathbf{S}_{i}=\mathbf{n}_{1}\cos(\mathbf{k}\cdot\mathbf{r}_{i})-\mathbf{n}_{2}\sin(\mathbf{k}\cdot\mathbf{r}_{i}),$
(2)

where $\mathbf{n}_{1}$ and $\mathbf{n}_{2}$ are orthonormal vectors. For the
Hamiltonian in Eq. 1 we obtain the expression of the energy as a function of
$\mathbf{k}$:

$\displaystyle E(\mathbf{k})$ $\displaystyle=$ $\displaystyle
2J_{1}N\left[\cos(k_{x})+2\cos(\frac{k_{x}}{2})\cos(\frac{\sqrt{3}k_{y}}{2})\right]+$
$\displaystyle+$ $\displaystyle
2J_{2}N\left[2\cos(\frac{3k_{x}}{2})\cos(\frac{\sqrt{3}k_{y}}{2})+\cos(\sqrt{3}k_{y})\right]+$
$\displaystyle+$ $\displaystyle
2J_{3}N\left[\cos(2k_{x})+2\cos(k_{x})\cos(\sqrt{3}k_{y})\right].$

Figure 9: The relation between (a) $J_{1}/J_{3}$ and $J_{2}/J_{3}$, and (b)
$J_{1}/J_{3}$ and K${}_{B}/J_{3}$ to obtain the ordering wave-vectors for
NiGa2S4 (blue) and FeGa2S4 (red and green) within the $J_{1}-J_{2}-J_{3}$ and
$J_{1}-J_{3}-K_{B}$ models. The critical points corresponding to a local
maximum, a local minimum and a global minimum of the energy at the given
values of the ordering wave-vectors are distinguished by the dotted, dashed
and solid lines, respectively.

It is important to note that there is no explicit dependence of the energy on
the vectors $\mathbf{n}_{1}$ and $\mathbf{n}_{2}$, reflecting the SO(3) global
symmetry of the model. In order to find the possible ordering vectors we
minimize $E(\mathbf{k})$ with respect of $\mathbf{k}$. Depending on the values
of the parameters $J_{1}$, $J_{2}$ and $J_{3}$, the minima of this model lie
along the following two sets of directions: (i): $\pm(1,0)$,
$\pm(\frac{1}{2},\frac{\sqrt{3}}{2})$, $\pm(-\frac{1}{2},\frac{\sqrt{3}}{2})$,
and (ii): $\pm(\frac{\sqrt{3}}{2},\frac{1}{2})$, $\pm(0,1)$ and
$\pm(\frac{\sqrt{3}}{2},\frac{1}{2})$. Since the ordering vector
$\bf{k}\approx$($\frac{1}{6}$ $\frac{1}{6}$ 0) in EN observed in MGa2S4 (M=
Ni, Fe) corresponds to the set (i), we focus on this direction and search for
a minimum along the $\mathbf{\hat{x}}$ direction (or equivalently along any
direction of this set). With this condition the value of $k$ should satisfy:

$\frac{J_{1}}{J_{3}}=-\frac{2(\sin(k)+\sin(2k))+3\frac{J_{2}}{J_{3}}\sin(\frac{3k}{2})}{\sin(\frac{k}{2})+\sin(k)}.$
(4)

From this equation we obtain linear relations between $J_{1}$, $J_{2}$ and
$J_{3}$ presented in Fig. 9 a. This linear relation expresses the values of
the parameters for which the energy develops a critical point. The critical
points gain stability when we move along each line, running from a local
maximum (represented by dots) to a local minimum (dashed segment) and finally
to a global minimum (solid line). For FeGa2S4 the Q2 vector (green line)
becomes a global minimum for $J_{1}/J_{3}>$ 1.72. Further restrictions on
$J_{1}/J_{3}$ should be included for these critical points to belong in the
set (i): for NiGa2S4 $J_{1}/J_{3}>$ -1 and to obtain either the Q1 or Q2
vectors for FeGa2S4 requires $J_{1}/J_{3}<$ 0.6 and $J_{1}/J_{3}<$ 4.3,
respectively.

The sign of the slope for a modulated state with $k\approx 0.155\times 4\pi$
for NiGa2S4 (blue) and with $k\approx 0.1737\times 4\pi$ for FeGa2S4 (red) is
different. This change of the sign takes place at the critical wave vector
$\bf{k}$=($\frac{1}{6}$ $\frac{1}{6}$ 0).
It is important to note that the single $J_{i=1,3}$ \- exchange models and the
bilinear $J_{1}-J_{2}$ and $J_{1}-J_{3}$ exchange models are not sufficient to
explain the Q1=0.6 Å-1 feature observed for FeGa2S4. For example, the ordering
vector $k\approx 0.155\times 4\pi$ for NiGa2S4 can be obtained from the
$J_{1}-J_{3}$ modelTamura2011 , this is not the case for $k\approx
0.1737\times 4\pi$ for FeGa2S4. In those models the ordering vector implies
either a uniform (ferromagnetic) state, a commensurate 120 $\deg$ structure or
it does not match the direction (the set indicated in (i)) signalled by the
experimental results. For the $J_{1}-J_{2}-J_{3}$ exchange model such
solutions, however, are available.
Now we consider the ordering vector with the length $k\approx 0.222\times
4\pi$, which could correspond to the center of the diffuse feature at Q2=0.77
Å-1 observed for FeGa2S4. The parameters of the ideal triangular lattice model
are related through a linear equation depicted in Fig. 9 by the green line.
These solutions strongly differ from the Q1 lines and require significant
$J_{1}$ coupling.
We proceed with the ’model engineering’ and raise the complexity by adding the
biquadratic interaction $K_{B}$ between first nearest neighbors
$K_{B}\sum_{<i,j>}({\bf{S}}_{i}\cdot{\bf{S}}_{j})^{2}$. We repeat the previous
procedures for two models: $J_{1}-J_{3}-K_{B}$ and $J_{1}-J_{2}-J_{3}-K_{B}$.
The energy of a modulated state for the $J_{1}-J_{3}-K_{B}$ model reads:

$\displaystyle E(\mathbf{k})$ $\displaystyle=$ $\displaystyle
2J_{1}N\left[\cos(k_{x})+2\cos(\frac{k_{x}}{2})\cos(\frac{\sqrt{3}k_{y}}{2})\right]+$
$\displaystyle+$ $\displaystyle
2J_{3}N\left[\cos(2k_{x})+2\cos(k_{x})\cos(\sqrt{3}k_{y})\right]+$
$\displaystyle+$ $\displaystyle
2K_{B}N\left[3+\cos(2k_{x})+2\cos(k_{x})\cos(\sqrt{3}k_{y})\right],$

As we did before we look for the minimum of the energy along the
$\mathbf{\hat{x}}$ direction. With this condition the value of $k$ should
satisfy:

$\frac{J_{1}}{J_{3}}=-\frac{2(\frac{K_{B}}{J_{3}}+1)(\sin(k)+\sin(2k))}{\sin(\frac{k}{2})+\sin(k)}.$
(6)

From this equation we obtain the relations between the values of $J_{1}$,
$J_{3}$ and $K_{B}$ (Fig. 9) right for which modulated states with $k\approx
0.155\times 4\pi$ for NiGa2S4, $k\approx 0.1737\times 4\pi$ and $k\approx
0.222\times 4\pi$ for FeGa2S4 are the ground states. The role of KB is
significant, but in the real system we expect KB to be a perturbative and not
the leading parameter.

Finally, if we add the $J_{2}$ exchange interaction to the previous model, the
energy as a function of $\mathbf{k}$ is:

$\displaystyle E(\mathbf{k})$ $\displaystyle=$ $\displaystyle
2J_{1}N\left[\cos(k_{x})+2\cos(\frac{k_{x}}{2})\cos(\frac{\sqrt{3}k_{y}}{2})\right]+$
$\displaystyle+$ $\displaystyle
2J_{2}N\left[2\cos(\frac{3k_{x}}{2})\cos(\frac{\sqrt{3}k_{y}}{2})+\cos(\sqrt{3}k_{y})\right]+$
$\displaystyle+$ $\displaystyle
2J_{3}N\left[\cos(2k_{x})+2\cos(k_{x})\cos(\sqrt{3}k_{y})\right]+$
$\displaystyle+$ $\displaystyle
2K_{B}N\left[3+\cos(2k_{x})+2\cos(k_{x})\cos(\sqrt{3}k_{y})\right].$

Since we are interested in the states with $\mathbf{k}$ along the
$\mathbf{\hat{x}}$ direction, we look for the minimum of the energy along this
direction. With this condition the value of $k$ should satisfy:

$\frac{J_{1}}{J_{3}}=-\frac{2(\frac{K_{B}}{J_{3}}+1)(\sin(k)+\sin(2k))+3\frac{J_{2}}{J_{3}}\sin(\frac{3k}{2})}{\sin(\frac{k}{2})+\sin(k)}.$
(8)

From this linear equation the values of $J_{1}$, $J_{2}$ and $K_{B}$ for which
the ICM ground states with $k\approx 0.155\times 4\pi$, $k\approx 0.1737\times
4\pi$ and $k\approx 0.222\times 4\pi$ are obtained.
Our analytical approach shows that the parameter space of the ideal triangular
model which fits the experimental ordering wave vectors is quite extended,
however constrained by Eq.8, and additional input is necessary to select the
unique set of exchange parameters. Despite the broad range of parameters, we
emphasize on the relevance of the constraints imposed by Eqs.4, 6 and 8 for
searching the correct set of couplings.

###  IV.2 Extraction of exchange parameters from the INS spectrum using linear
spin wave theory

Figure 10: Maximum log-likelihood after optimising the magnitude of the spin
wave contributions.

Figure 11:  Predicted $E_{i}=24.3$ meV (left) and $E_{i}=7.7$ meV (right) INS
spectra at 2 K, constructed from an empirical background model measured at 180
K and spin wave simulations for the set of the exchange values $J_{1}=1.7$
meV, $J_{2}=0.9$ meV, $J_{3}=0.8$ meV. The margin plots show comparisons with
data measured at 2 K.

We narrow the exchange parameter space calculating the spin excitation
spectrum and fitting it to the experimental S(Q, $\omega$). We consider the
triangular lattice Heisenberg bilinear model of Eq.1 and focus on the Q2=0.77
Å-1 feature observed in the INS spectra assuming that it emerges from the spin
waves. We take the 180 K data after the Bose factor correction as an empirical
background (even though it contains a featureless magnetic contribution) and
attribute the rest of the 2 K data to spin wave excitations. Using the ground
state of Eq.IV.1, we calculate spin wave excitations with the SpinW
programSpinW  and compare the resulting spectra with the INS data measured at
2 K. To avoid regions in data that are saturated by elastic scattering, only
energy transfers above 2.4 meV and 0.75 meV are considered for the two setups
of $E_{i}=24.3$ meV and $E_{i}=7.7$ meV, respectively. In the simulation, the
experimental energy resolution of each setup is approximated by a constant
corresponding to its nominal (energy-dependent) value at the peak of the
magnetic excitation. No intrinsic broadening is considered. The nominal
uncertainties on the measured intensity are obtained from the data itself and
are assumed to be underestimated for bins with low neutron counts. Therefore,
these uncertainties are thresholded to a minimum value determined by the
actual variance across low-intensity bins in regions of data without large
scale features. The model is evaluated on a discretised grid of values ranging
from -2 meV to 2 meV in increments of 0.1 meV for each of exchange parameters
$J_{1}$, $J_{2}$, $J_{3}$. Figure 10 shows the maximum log-likelihood at each
point after optimising the magnitude of the spin wave contributions. The
maximum log-likelihood suggests that $J_{1}$, $J_{2}$ and $J_{3}$ are
antiferromagnetic. The best fit is obtained for $J_{1}=1.7$ meV, $J_{2}=0.9$
meV, $J_{3}=0.8$ meV, which is consistent with the relation imposed by Eq. 4.
We note that the general locations of the observed excitations can be
approximately reproduced by this model, however the detailed structure of the
simulated spin waves is not evident in the data, and this set of parameters is
out the range indicated in Figure 9 a. In particular the $E_{i}=7.7$ meV setup
reveals that the predicted drop in intensity at low energy transfer is less
rapid than what can be observed in data. These discrepancies are not removed
by inclusion of the single-ion anisotropy $D_{n}\sum_{i}(S_{i})^{2}$, which
should be 20% of exchange and planar according to the magnetic susceptibility
measurements.
We conclude that the best fit to our powder INS data provides the set of
exchange parameters within the $J_{1}-J_{2}-J_{3}$ bilinear model but this
model is not sufficient to explain the whole observed excitation spectrum.

##  V Summary and Discussion

Our extensive experimental study of the triangular antiferromagnet FeGa2S4
implies that this material, and most probably also NiGa2S4, has significant
mixing of the M (M = Ni, Fe) and Ga cations between the octahedrally
coordinated 1(b) and tetrahedrally coordinated 2(d) positions. Nearly 20% of
the 1(b) sites forming magnetic triangular layers are occupied by nonmagnetic
ions and nearly 10% of the 2(d) sites forming non-magnetic triangular layers
host magnetic M-ions.

This mixing is an important ingredient of highly nontrivial magnetic
properties. It introduces a long-rage magnetic order at TN=5.5 K revealed by
the Q1=0.6 Å-1 peak in neutron powder diffraction, which we associated with
the $\bf{k}$=(0.1737(1) 0.1737(1) 0) propagation vector. This magnetic order
is dominated by the amplitude modulated $z$-component of the minor Fe-spins at
the 2(d) sites. The amount of the minor Fe-ions are well below the percolation
limit of the stacked triangular antiferromagnetic lattice
(p${}_{c}^{side}$=0.26240(5)Schrenk2014 ), so the layers of the major Fe ions
possibly mediate this exchange.

The broad feature at Q2=0.77 Å-1 in neutron powder diffraction corresponds to
short-range magnetic correlations within the triangular layers of the major Fe
sites. These in-plane correlations are static at low temperatures, evolve into
dynamic correlations at elevated temperatures 30 - 50 K and fade away at
temperatures comparable to $|\Theta|_{\rm CW}$=160 K. These dynamic
correlations could be confined spin waves, the so-called Halperin-Saslow
modesHalperin1977 ; Podolsky2009  of the major Fe moments, which freeze in due
to numerous nonmagnetic impurities. The nonmagnetic defects possibly trigger
formation-dissociation of the Z2 vortices. According to the MC
simulationsAjiro1988  the characteristic temperature TKM related to the
pairing-dissociation of the Z2 vortices is correlated with the effective
exchange as TKM =0.66 $J$ S2. For our case $J$ = 13 K (estimated from
susceptibility) and S=2, which this results in TKM = 34 K. This temperature
corresponds to the anomalies found in muon spectroscopy (T*= 33 K) Reotier2012
; Zhao2012  and to the change between static and dynamic nature of magnetic
soft range correlations in our neutron scattering. Thus, even if the MGa2S4
family is not an ideal TLAFM system, it might realize interesting phenomena
predicted theoretically for this model.
Our extensive theoretical study of the ideal triangular lattice bilinear-
biquadratic Heisenberg model has value by itself. The analytical study
provided in Section IV.1 could be applied to a generic triangular system and
as such serves as guide to study other compounds. It rationalizes the ratios
$J_{1}/J_{3}$, $J_{2}/J_{3}$ and $K_{B}/J_{3}$ for a provided ground state
propagation vector. For the specific case of FeGa2S4 it notifies that the Q1
and Q2 features require possibly different sets of exchange parameters as they
reside on the lines with one intercept point. We reduce the range of possible
exchange parameters for FeGa2S4 by modelling the distinct features of the
experimental low-temperature INS spectra as conventional spin waves within the
$J_{1}-J_{2}-J_{3}$ bilinear Hamiltonian. The best fit is obtained for the set
of AF couplings $J_{1}$=1.7 meV, $J_{2}$=0.9 meV, $J_{3}$=0.8 meV, which is
close to the analytical solutions for Q2. However, only the gross features of
the INS spectra are explained. The origin of the broad scattering background
remains elusive and for its understanding it is important to disentangle
magnetic properties of inversion-free and inversion-containing triangular
antiferromagnetic materials.
Our results imply that the short-range correlations at Q2 and the anomalies
around 30 K - 50 K are the properties of the major magnetic sites of the
FeGa2S4 system. The main exchange couplings are antiferromagnetic and
$J_{1}>J_{3}\approx J_{2}$. How strongly the behaviour of major Fe triangular
sublattice is influenced by the minor Fe sites remains unclear. Theoretical
predictions are presently limited to only small amount of impurities and to
the first nearest-neighborsWollny2011 ; Maryasin2013 . Extensions to the
experimentally found values can bring answers to this question. Preliminary
Monte Carlo simulations with the experimentally observed amount of impurity
seem to reproduce the Q1 and Q2 features observed in inelastic neutron
scattering. However further exploration would be necessary. Novel synthetic
routes to obtain inversion free members of the MGa2S4 family could in the
future facilitate experimental validation of the ideal TLAFM model.

###### Acknowledgements.

This work was performed at SINQ, Paul Scherrer Institute, Villigen,
Switzerland with financial support of the Swiss National Science Foundation
(Grant No. 200020-182536) and Centro Lanitoamericano-Suizo (Seed money Grant
No. SMG1811) and the Swiss Data Science Center (Project BISTOM C17-12). We
acknowledge F. Krumeich for recording SEM images and EDXS analysis.

##  VI Appendix A: Materials and methods

Polycrystalline FeGa2S4 and NiGa2S4 samples were prepared by solid state
synthesis from FeS (NiS) and GaS starting materials. The mixture was annealed
for one week at 900 deg C in an evacuated quartz ampule. Then the sample was
homogenised by grinding, resealed, annealed again and cooled within 24 hours.
Single crystals were grown by the chemical vapor transport method using iodine
as a transport agent and a temperature gradient 1000/900 deg C. The crystals
were then quenched.

The purity of the samples and the details of the crystal structure for FeGa2S4
in the temperature range 5 K - 300 K were determined for polycrystalline
material on the x-ray powder diffractometer of the X04SA Material Science
beamline at the SLS synchrotron ($\lambda$=0.641090 Å) and for FeGa2S4 and
NiGa2S4 single crystals on a Bruker SMART diffractometer equipped with an Apex
I detector and an Apex II detector (Mo K$\alpha$=0.71073 Å) at ETH Zürich at
300 and 100 K, respectively. CSD 2078453 and CSD 2078454 contain the
supplementary crystallographic data for FeGa2S4 and NiGa2S4, respectively.
These data can be obtained free of charge from FIZ Karlsruhe via
www.ccdc.cam.ac.uk/structures. Powder patterns were refined using the Fullprof
suitefullprof1993 , for single crystal refinements Shelxshelx  was used.
Scanning electron microscopy (SEM) images of the as-obtained samples were
recorded with a secondary electron detector on a Gemini 1530 (Zeiss)
microscope. Energy-dispersive X-ray spectroscopy (EDXS) was performed with a
SSD system (Noran) attached to this microscope.

Magnetic susceptibility was measured in the temperature range 1.8 K - 350 K in
an applied magnetic field of 0.1 T using a Magnetic Properties Measurement
System superconducting quantum interference device (SQUID) magnetometer,
Quantum Design.

Neutron powder diffraction (NPD) patterns for FeGa2S4 were measured on the DMC
diffractometer at SINQ ($\lambda$=2.4575 Å). Data were collected at several
temperatures between 1.3 K and 180 K with the measuring time up to 12 h per
pattern. XYZ-polarized neutron powder diffuse scattering was measured on the
D7 spectrometer at ILL (Ei=3.8 meV) at 1.4 K (4 h) and 60 K (10 h). The
magnetic scattering was extracted from isotropic magnetic scattering the
magnetic ${({\frac{d\sigma}{d\Omega}})}_{mag}$, incoherent
${({\frac{d\sigma}{d\Omega}})}_{inc}$ and nuclear
${({\frac{d\sigma}{d\Omega}})}_{nuc}$ cross sections were evaluated by the
following equationsStewart2008 :

${({\frac{d\sigma}{d\Omega}})}_{mag}=2{({\frac{d\sigma}{d\Omega}})}^{x}_{sf}+2{({\frac{d\sigma}{d\Omega}})}^{y}_{sf}-4{({\frac{d\sigma}{d\Omega}})}^{z}_{sf}$
(9)

where $x,y,z$ refer to the direction of the incident polarization, $sf$ and
$nsf$ stands for spin-flip and non-spin-flip.

Time-of-flight inelastic neutron scattering (INS) was measured for FeGa2S4 on
the LET spectrometerBewleyLET  at ISIS. We used the Rep-Rate Multiplication
mode with the Fermi chopper rotating at 240 Hz in the high-flux setup, which
allowed in a single time frame to collect data from incident neutron energies
of 1.5, 2.24, 3.38, 7.7 and 24.3 meV. The Al can was subtracted and data were
corrected for the Bose factor.

##  VII Appendix B: Conventions and notations for the reciprocal space

Two different notations can be encountered in the literature depending on the
type of article (experimental or pure theoretical). We distinguish them by the
names “experimental notation” (EN) and “theoretical notation” (TN).

In the EN notation a vector $\mathbf{k}$ in the reciprocal space is expressed
in terms of the basis vectors $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ as follows
$\mathbf{k}=k_{1}\mathbf{b}_{1}+k_{2}\mathbf{b}_{2}$. In the TN notation a
vector $\mathbf{k}$ is expressed in terms of the basis vectors
$\hat{\mathbf{x}}$ and $\hat{\mathbf{y}}$ as follows
$\mathbf{k}=k_{x}\hat{\mathbf{x}}+k_{y}\hat{\mathbf{y}}$.

Figure 12:  The red arrows represent the $\hat{\mathbf{x}}$ and
$\hat{\mathbf{y}}$ unit vectors in all the cases. (a) Basis vectors
$\\{\mathbf{a}_{1},\mathbf{a}_{2}\\}$ (in green) for the real space of the
triangular lattice. (b) The first Brillouin zone for the triangular lattice.
The M point (green circles), K point (purple circle) and $\Gamma$ point (black
circle). The blue arrows represent the basis vectors $\mathbf{b}_{1}$ and
$\mathbf{b}_{2}$ in the reciprocal space. (c) Some examples of interest, from
the $\Gamma$ point (0, 0, 0), in light-blue, to the K point
($\frac{1}{3}$,$\frac{1}{3}$,0), in light-green, we have the points (expressed
in CN): ($\frac{1}{6}$,$\frac{1}{6}$,0), in purple, (0.1737, 0.1737,0), in
red, and (0.222, 0.222,0), in brown.

The connection between both notations is given by a change of basis (from
$\\{\mathbf{b}_{1},\mathbf{b}_{2}\\}$ to
$\\{\hat{\mathbf{x}},\hat{\mathbf{y}}\\}$). The basis vectors $\mathbf{b}_{1}$
and $\mathbf{b}_{2}$ (in blue) are expressed in terms of the vectors
$\hat{\mathbf{x}}$ and $\hat{\mathbf{y}}$ (in red) as follows:

$\displaystyle\mathbf{b}_{1}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{a}\hat{\mathbf{x}}+\frac{2\pi}{\sqrt{3}a}\hat{\mathbf{y}},$
(10) $\displaystyle\mathbf{b}_{2}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{a}\hat{\mathbf{x}}-\frac{2\pi}{\sqrt{3}a}\hat{\mathbf{y}}.$
(11)

So, for instance, in EN the state $({{\frac{1}{3}}},{{\frac{1}{3}}},0)$ (the
K-point) corresponds to the triplet $(\frac{4\pi}{3a},0,0)$ in TN,

$\frac{1}{3}\mathbf{b}_{1}+\frac{1}{3}\mathbf{b}_{2}=\frac{4\pi}{3}\hat{\mathbf{x}}.$
(12)

The M-point in EN is identified with ({$\frac{1}{2}$, 0, 0), and in the TN
with the triplet $(\frac{\pi}{a},\frac{\pi}{\sqrt{3}a},0)$,

$\frac{1}{2}\mathbf{b}_{1}=\frac{\pi}{a}\hat{\mathbf{x}}+\frac{\pi}{\sqrt{3}a}\hat{\mathbf{y}}.$
(13)

Thus, when in the EN notation we refer to the $(0.155,0.155)$ wave vector, it
corresponds to the wave vector $k\approx 0.155\times 4\pi$ in the TN notation
(taking $a=1$).

## References

  * (1) P. W. Anderson, Mat. Res. Bull 8, 153 (1973).
  * (2) Th. Jolicoeur, E. Dagotto, R. Gagliano, S. Bacci, Phys. Rev. B 42, 4800(R) (1990).
  * (3) T. Okubo, S. Chung, and H. Kawamura, Phys. Rev. Lett. 108, 017206 (2012).
  * (4) L. D. Smiri, N. H. Dung, and M. P. Pardo, Materials Research Bulletin 15 861 (1980).
  * (5) S. Nakatsuji, H. Tonomura, K. Onuma, Y. Nambu, O. Sakai, Y. Maeno, R. T. Macaluso, and J. Y. Chan, Phys. Rev. Lett. 99, 157203 (2007).
  * (6) B. R. Myoung, S. J. Kim, and C.S. Kim, J. Korean Phys. Soc. 53, 750 (2008).
  * (7) S. Zhao, P. Dalmas de Réotier, A. Yaouanc, D. E. MacLaughlin, J. M. Mackie, O. O. Bernal, Y. Nambu, T. Higo, and S. Nakatsuji, Phys. Rev. B 86, 064435 (2012).
  * (8) C. Stock, S. Jonas, C. Broholm, S. Nakatsuji, Y. Nambu, K. Onuma, Y. Maeno, and J. H. Chung, Phys. Rev. Lett. 105, 037402 (2010).
  * (9) P. Dalmas de Réotier, A. Yaouanc, D. E. MacLaughlin, S. Zhao, T. Higo, S. Nakatsuji, Y. Nambu, C. Marin, G. Lapertot, A. Amato, and C. Baines, Phys. Rev. B 85, 140407(R) (2012).
  * (10) S. Nakatsuji, Y. Nambu, H. Tonomura, O. Sakai, S. Jonas, C. Broholm, H. Tsunetsugu, Y. Qiu, and Y. Maeno, Science 309, 1697 (2005).
  * (11) I. I. Mazin, Phys. Rev. B 76, 140406(R) (2007).
  * (12) H. Tsunetsugu and M. Arikawa, J. Phys. Soc. Jpn 75, 083701 (2006).
  * (13) A. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. 97, 087205 (2006).
  * (14) S. Bhattacharjee, V. B. Shenoy, and T. Senthil, Phys. Rev. B 74, 092406 (2006).
  * (15) M. E. Valentine, T. Higo, Y. Nambu, D. Chaudhuri, J. Wen, C. Broholm, S. Nakatsuji, and N. Drichko, Phys. Rev. Lett. 125, 197201 (2020).
  * (16) H. Tsunetsugu and M. Arikawa, J. Phys.: Condens. Matter 19,145248 (2007).
  * (17) D. Podolsky and Y.B. Kim, Phys. Rev. B 79, 140402(R) (2009).
  * (18) Details of sample preparation and experimental setups are documented in Appendix VI.
  * (19) H. Takeya, K. Ishida, K. Kitagawa, Y. Ihara, K. Onuma, Y. Maeno, Y. Nambu, S. Nakatsuji, D. E. MacLaughlin, A. Koda, and R. Kadono, Phys. Rev. B 77, 054429 (2008).
  * (20) Y. Nambu, R.T. Macaluso, T. Higo, K. Ishida, and S. Nakatsuji, Phys. Rev. B 79, 214108 (2009).
  * (21) B. E. Warren, Phys. Rev. 59, 693 (1941).
  * (22) J. A. M. Paddison, J. R. Stewart, and A. L. Goodwin, J. Phys.: Condens. Matter 25,454220 (2013).
  * (23) R. Tamura and N. Kawashima, J. Phys. Soc. Jpn 77, 103002 (2008).
  * (24) R. Tamura and N. Kawashima, J. Phys. Soc. Jpn 80, 074008 (2011).
  * (25) https://spinw.org
  * (26) K. J. Schrenk, N. A. M. Araújo, and H. J. Herrmann, Phys. Rev. E 87, 032123 (2013).
  * (27) B. I. Halperin and W. M. Saslow, Phys. Rev. B 16, 2154 (1977).
  * (28) Y. Ajiro, H. Kikuchi, S. Sugiyama, T. Nakashima, S. Shamoto, N. Nakayama, M. Kiyama, N. Yamamoto, and Y. Oka, J. Phys. Soc. Jpn 57, 2268 (1988).
  * (29) A. Wollny, L. Fritz, and M. Vojta, Phys. Rev. Lett. 107, 137204 (2011).
  * (30) V. S. Maryasin and M. E. Zhitomirsky, Phys. Rev. Lett. 111, 247201 (2013).
  * (31) J. Rodríguez-Carvajal, Phys. B: Condens. Matt. 192, 55 (1993).
  * (32) G. M. Sheldrick, Acta Cryst. C71, 3(2015).
  * (33) J. R. Stewart, P. P. Deen, K. H. Andersen, H. Schober, J. F. Barthelemy, J. M. Hillier, A. P. Mjurani, T. Hayes, and B. Lindenau, J Appl. Cryst. 42, 69 (2009).
  * (34) R. I. Bewley, J. W. Taylor, S. M. Bennington, Nucl. Instrum. Methods Phys. Res. Sect. A 637, 128 (2011).
  * (35) H. Yamaguchi, S. Kimura, M. Hagiwara, Y. Nambu, S. Nakatsuji, Y. Maeno, A. Matsuo, and K. Kindo, J. Phys. Soc. Jpn 79, 054710 (2010).