©2017 International World Wide Web Conference Committee (IW3C2), published
under Creative Commons CC BY 4.0 License.

# Leveraging Large Amounts of Weakly Supervised Data for Multi-Language
Sentiment Classification

Jan Deriu
Aurelien Lucchi
Valeria De Luca
ZHAW  ETH Zurich  ETH Zurich  Aliaksei Severyn
Simon Müller
Mark Cieliebak
Google Research  SpinningBytes AG  SpinningBytes AG  Thomas Hofmann
Martin Jaggi
ETH Zurich  EPFL

(3 April 2017)

###### Abstract

This paper presents a novel approach for multi-lingual sentiment
classification in short texts. This is a challenging task as the amount of
training data in languages other than English is very limited. Previously
proposed multi-lingual approaches typically require to establish a
correspondence to English for which powerful classifiers are already
available. In contrast, our method does not require such supervision. We
leverage large amounts of weakly-supervised data in various languages to train
a multi-layer convolutional network and demonstrate the importance of using
pre-training of such networks. We thoroughly evaluate our approach on various
multi-lingual datasets, including the recent SemEval-2016 sentiment prediction
benchmark (Task 4), where we achieved state-of-the-art performance. We also
compare the performance of our model trained individually for each language to
a variant trained for all languages at once. We show that the latter model
reaches slightly worse - but still acceptable - performance when compared to
the single language model, while benefiting from better generalization
properties across languages.

###### keywords:

Sentiment classification; multi-language; weak supervision; neural networks

††conference: WWW 2017, April 3–7, 2017, Perth, Australia.

##  1 Introduction

Automatic sentiment analysis is a fundamental problem in natural language
processing (NLP). A huge volume of opinionated text is currently available on
social media. On Twitter alone, 500 million tweets are published every day.
Being able to manually process such a high volume of data is beyond our
abilities, thus clearly highlighting the need for automatically understanding
the polarity and meaning of these texts. Although there have been several
progresses towards this goal, automatic sentiment analysis is still a
challenging task due to the complexity of human language, where the use of
rhetorical constructions such as sarcasm and irony easily confuse sentiment
classifiers. Contextualization and informal language, which are often adopted
on social media, are additional complicating factors. The Internet is also
multi-lingual and each language has its own grammar and syntactic rules.

Given all these difficulties, it is not surprising that the performance of
existing commercial systems is still rather poor, as shown in several recent
studies [6, 27]. The benchmark work of Ribeiro et al. [27] showed that even
the performance of the best systems largely varies across datasets and overall
leaves much room for improvement. Hence it is important to design a method
that generalizes well to different domains and languages.

Contributions. The majority of current research efforts in sentiment analysis
focuses on the English language. This is partially due to the large number of
resources available in English, including sentiment dictionaries, annotated
corpora and even benchmark datasets. An example is the SemEval competition,
which is one of the largest competitions on semantic text evaluation and
covers several tasks for sentiment analysis [24].

However, only 26.3% of the total number of internet users in 2016 are English
speakers [15] and only 34% of all tweets are written in English [21]. Hence
there is a strong incentive to develop methods that work well with other
languages. In this work, we focus on the question of how sentiment analysis
can be done for multiple languages by leveraging existing technologies. Our
method is the state-of-the-art approach for sentiment analysis on Twitter data
which recently won the SemEval-2016 competition [9]. Here we additionally
explore how to best adapt this approach to other languages. The core component
of our system is a multi-layer convolutional neural network (CNN), trained in
three phases: i) unsupervised phase, where word embeddings are created on a
large corpus of unlabeled tweets; ii) distant supervised phase, where the
network is trained on a weakly-labeled dataset of tweets containing emoticons;
and iii) supervised phase, where the network is finally trained on manually
annotated tweets. For English, this system achieved an F1-score of 62.7% on
the test data of SemEval-2016 [9].

Although existing CNN approaches [31, 9] can a-priori be trained on any
language other than English, these nevertheless require a large amount of
training data. Yet resources in languages other than English are lacking, and
manually labeling tweets is a time-consuming and expensive process. Two
straightforward solutions that do not require manual work can be envisioned:
(1) automatically translate the data into English and run the existing English
classifier; or (2) train a CNN using only weakly-labeled tweets without using
any supervised data. It is expected that a fully-trained CNN would perform
better than the aforementioned cases. However, it is unclear if such
improvement is significant and justifies the need of manually labeling
thousands of tweets.

In this paper, we investigate how to effectively train and optimize a CNN for
multi-lingual sentiment analysis. We compare the performance of various
approaches for non-English texts. In details, our main contributions are:

  * •

An evaluation of the state-of-the-art CNN approach similar to the one proposed
in [9] on three new languages, namely French, German and Italian

  * •

A thorough analysis of the influence of network parameters (number of layers,
hyper-parameters) and other factors, e.g. the amount of distant-supervised and
supervised data, on end-to-end performance

  * •

For each language, a comparison of various approaches for sentiment analysis:
(i) full training of the CNN for the considered language; and (ii)
automatically translating the texts into a language (English) where a
sentiment classifier already exists. Other baseline methods, described in the
experimental section, are also compared

  * •

In addition, we show that a single CNN model can be successfully trained for
the joined task on all languages, as opposed to separate networks for each
individual language. This approach has the advantages of removing the reliance
on (possibly inaccurate) language identification systems and it can be easily
extended to new languages and multi-language texts. We provide detailed
comparison to similar per-language models, and show that the proposed joint
model still performs relatively well

  * •

Public release of the source code as well as pre-trained models for all
languages tested in this paper, on
http://github.com/spinningbytes/deep-mlsa

##  2 Related work

In the following, we provide an overview of the most relevant works, related
to the application of neural networks to sentiment classification, distant
supervision and training multi-lingual text classifiers.

Neural networks. Neural networks have shown great promise in NLP over the past
few years. Examples are in semantic analysis [33], machine translation [12]
and sentiment analysis [34]. In particular, shallow CNNs have recently
improved the state-of-the-art in text polarity classification demonstrating a
significant increase in terms of accuracy compared to previous state-of-the-
art techniques [18, 17, 10, 32, 16, 28, 9]. These successful CNN models are
characterized by a set of convolution filters acting as a sliding window over
the input word sequence, typically followed by a pooling operation (such as
max-pooling) to generate a fixed-vector representation of the input sentence.

CNNs vs RNNs. Recently, recurrent neural network architectures (RNNs), such as
long short-term memory networks (LSTMs), have received significant attention
for various NLP tasks. Yet these have so far not outperformed convolutional
architectures on polarity prediction [29, Table 4]. This has been evidenced by
the recent SemEval-2016 challenge [24], where systems relying on convolutional
networks rank at the top. In fact, long-term relationships captured well by
LSTMs are of minor importance to the sentiment analysis of short tweets. On
the contrary, learning powerful $n$-gram feature extractors (which
convolutional networks handle very well) contributes much more to the
discriminative power of the model, since these are able to effectively detect
sentiment cues. Additionally, LSTMs are much more computationally expensive
than CNNs, preventing their application to very large collections like the one
used in this paper (hundreds of millions tweets).

Distant-supervised learning. The use of semi-supervised or unsupervised
learning has been an active research direction in machine learning and
particularly for various NLP applications. There is empirical evidence that
unsupervised training can be beneficial for supervised machine learning tasks
[11]. In this paper, we consider a variant of unsupervised learning named
distant pre-training which consists in inferring weak labels from data without
manual labels. This approach has been used for text polarity classification
where significantly larger training sets were generated from texts containing
emoticons [13, 32]. Severyn and Moschitti [32] have shown that training a CNN
on these larger datasets, followed by additional supervised training on a
smaller set of manually annotated labels, yields improved performance on
tweets.

Multi-language sentiment classification. Sentiment classification has drawn a
lot of attention in the past few years both in industry and academia [24, 6].
Yet most of the research effort has been focusing on tweets written in one
language (mostly English). One exception is the work of Boiy and Moens [4]
that studied the portability of a learned sentiment classification model
across domains and languages. They focused on French, Dutch and English, and
showed that significant disparities between these languages can severely
hinder the performance of a classifier trained on hand-crafted features.

The major factor that limits the development of accurate models for multi-
lingual sentiment analysis is the lack of supervised corpora [2, 7]. Most of
the existing approaches addressing this problem [22, 1] try to transfer
knowledge from English – for which tools, labelled data and resources are
abundant – to other languages for which resources are rather limited. An
example is the approach introduced in [22], which transfers hand-crafted
subjectivity annotation resources - such as a per-word sentiment lexicon -
from English to Romanian. A similar approach introduced in [1] consists in
translating a target language to English and then to use an English sentiment
classifier rather than one specific to the target language. Several approaches
have also been proposed to build distributed representations of words in
multiple languages. The work of Wick et al. [35] used a Wikipedia corpus of
five languages to train word embeddings, and then used anchor terms (names,
cross-lingual words) to align the embeddings. Gouws et al. [14] proposed a
method to create bilingual word vectors by requiring words that are related
across the two languages.

All the aforementioned approaches rely on having access to a set of
correspondences between English and the target language. Some of these methods
also require translating the target language to English. Yet machine
translation is a very challenging task in NLP and represents an additional
source of error in the classification system, due to various problems such as
sparseness and noise in the data [7]. Furthermore, such methods crucially rely
on accurate language identification, which is a very difficult task,
especially on short texts. See e.g. [20, 19] for an overview of these methods
and their limitations in generalizing to different domains.

In this work, we also investigate the performance of a language-independent
classifier consisting of a CNN trained on all languages at once. This approach
is similar to the Naïve Bayes classifier proposed in [25], excepts that it
relies on simple hand-crafted word-level features instead of the CNN
architecture used in this work.

Figure 1: Architecture of the proposed CNN model with 2 convolutional layers

##  3 Model

Our model follows a multi-layer CNN architecture, which we firstly introduced
in [9]. Given an input sequence of words, the corresponding sequence of word
embeddings is fed as input to the first 1d convolutional layer. Each
convolutional filter here operates in a sliding window fashion along the input
dimension (details are described below). This layer is followed by a max-
pooling operation whose output is then fed into the next convolutional layer.
We extend a single-layer CNN, originally proposed in [32, 18, 17], to two
convolutional and pooling layers. The resulting network architecture is
illustrated in Figure 1 and in its basic variant consists of two consecutive
pairs of convolutional-pooling layers followed by a single hidden layer and a
soft-max output layer. In the following, we describe in detail each layer and
corresponding parameters.

###  3.1 Convolutional Neural Network

Embedding layer. Each word is associated with a $d$-dimensional vector
(embedding). An input sequence of $n$ words is represented by concatenating
their embeddings, yielding a sentence matrix
$\mathbf{X}\in\rm{\mathbb{R}}^{d\times n}$. $\mathbf{X}$ is used as input to
the network.

Convolutional layer. This layer applies a set of $m$ convolutional filters of
length $h$ over the matrix $\mathbf{X}$. Let $\mathbf{X}_{[i:i+h]}$ denote the
concatenation of word vectors $\mathbf{x}_{i}$ to $\mathbf{x}_{i+h}$. A
feature $c_{i}$ is generated for a given filter $\mathbf{F}$ by:

$c_{i}:=\sum_{k,j}(\mathbf{X}_{[i:i+h]})_{k,j}\cdot\mathbf{F}_{k,j}\vspace{-1pt}$
(1)

The concatenation of all vectors in a sentence defines a feature vector
$\mathbf{c}\in\rm{\mathbb{R}}^{n-h+1}$. The vectors $\mathbf{c}$ are then
aggregated from all $m$ filters into a feature map matrix
$\mathbf{C}\in\rm{\mathbb{R}}^{m\times(n-h+1)}$. The filters are learned
during the training phase of the neural network, as described in Section 3.2.
The output of the convolutional layer is passed through a non-linear
activation function, before entering a pooling layer.

Pooling layer. The pooling layer aggregates the input vectors by taking the
maximum over a set of non-overlapping intervals. The resulting pooled feature
map matrix has the form:
$\mathbf{C_{pooled}}\in{\mathbb{R}}^{m\times\frac{n-h+1}{s}}$, where $s$ is
the length of each interval. In the case of overlapping intervals with a
stride value $s_{t}$, the pooled feature map matrix has the form
$\mathbf{C_{pooled}}\in{\mathbb{R}}^{m\times\frac{n-h+1-s}{s_{t}}}$. Depending
on whether the boundaries are included or not, the result of the fraction is
rounded up or down respectively.

Hidden layer. A fully connected hidden layer computes the transformation
$\alpha(\mathbf{W}*\mathbf{x}+\mathbf{b})$, where
$\mathbf{W}\in\rm{\mathbb{R}}^{m\times m}$ is the weight matrix,
$\mathbf{b}\in\rm I\\!R^{m}$ the bias, and $\alpha$ the rectified linear
(`relu`) activation function [23]. The output
$\mathbf{x}\in\rm{\mathbb{R}}^{m}$ of this layer can be seen as an embedding
of the input sentence.

Softmax. Finally, the outputs of the previous layer
$\mathbf{x}\in\rm{\mathbb{R}}^{m}$ are fully connected to a soft-max
regression layer, which returns the class $\hat{y}\in[1,K]$ with largest
probability, i.e.,

$\displaystyle\begin{split}\hat{y}&=\arg\max_{j}\,P(y=j\ |\
\mathbf{x},\mathbf{w},\mathbf{a})\\\
&=\arg\max_{j}\,\frac{e^{\mathbf{x}^{\intercal}\mathbf{w}_{j}+a_{j}}}{\sum_{k=1}^{K}e^{\mathbf{x}^{\intercal}\mathbf{w}_{k}+a_{j}}},\end{split}$
(2)

where $\mathbf{w}_{j}$ denotes the weights vector of class $j$, from which the
dot product with the input is formed, and $a_{j}$ the bias of class $j$.

Network Parameters. The following parameters of the neural network are learned
during training: $\mathbf{\theta}=$
$\\{\mathbf{X},\mathbf{F}_{1},\mathbf{b}_{1},$
$\mathbf{F}_{2},\mathbf{b}_{2},\mathbf{W},\mathbf{a}\\}$, with $\mathbf{X}$
the word embedding matrix, where each row contains the $d$-dimensional
embedding vector for a specific word; $\mathbf{F}_{i},\mathbf{b}_{i}$ the
filter weights and biases of convolutional layers; and $\mathbf{W}$ and
$\mathbf{a}$ the weight-matrix for output classes in the soft-max layer.

###  3.2 Learning the Model Parameters

The model parameters are learned using the following three-phase procedure:
(i) creation of word embeddings; (ii) distant-supervised phase, where the
network parameters are tuned by training on weakly labelled examples; and
(iii) final supervised phase, where the network is trained on the supervised
training data.

Preprocessing and Word Embeddings. The word embeddings are learned on an
unsupervised corpus containing $300M$ tweets. We apply a skip-gram model of
window-size 5 and filter words that occur less than 15 times [32]. The
dimensionality of the vector representation is set to $d=52$. Our experiments
showed that using a larger dimension did not yield any significant
improvement.

Training. During the first distant-supervised phase, we use emoticons to infer
noisy labels on tweets in the training set [26, 13]. Details of the training
set obtained from this procedure are listed in Table 1, "Pre-training" rows.
Note that the data used for this pre-training phase do not necessarily overlap
with the data used to create the word embeddings. The neural network was
trained on these datasets for one epoch, before finally training on the
supervised data for about 20 epochs. The word embeddings
$\mathbf{X}\in{\mathbb{R}}^{d\times n}$ are updated during both the distant-
and supervised training phases by applying back-propagation through the entire
network. The resulting effect on the embeddings of words are illustrated in
Figure 7 and discussed in Section 5.

Optimization. During both training phases, the network parameters are learned
using _AdaDelta_ [36]. We compute the score on the validation set at fixed
intervals and select the parameters achieving the highest score.

Figure 2 shows a complete overview of the three phases of the learning
procedure.

Figure 2: Training Phases Overview.

##  4 Experimental Setting

Table 1: Data used for training the CNN model. Language | Dataset | Total | Neutral | Neg. | Pos.
---|---|---|---|---|---
English | Word embeddings | 300M | - | - | -
Pre-training | 60M | - | 30M | 30M
Training | 18044 | 7544 | 2748 | 7752
Validation | 2390 | 987 | 365 | 1038
Test | 20632 | 10342 | 3231 | 7059
French | Word embeddings | 300M | - | - | -
Pre-training | 60M | - | 20M | 40M
Training | 9107 | 4095 | 2120 | 2892
Test | 3238 | 1452 | 768 | 1018
German | Word embeddings | 300M | - | - | -
Pre-training | 40M | - | 8M | 32M
Training | 8955 | 5319 | 1443 | 2193
Test | 994 | 567 | 177 | 250
Italian | Word embeddings | 300M | - | - | -
Pre-training | 40M | - | 10M | 30M
Training | 6669 | 2942 | 2293 | 1434
Test | 741 | 314 | 250 | 177

###  4.1 Data

We used a large set of $300M$ tweets to create the word embeddings for each
language, as well as a distant supervision corpus of $40$-$60M$ tweets for
each language, where each tweet contained at least one emoticon (positive or
negative smiley). Smileys were used to automatically infer weak labels and
subsequently removed from the tweets. This idea of distant-supervised learning
was described in [13, 32, 31]. For the final supervised phase, we used
publicly available labeled datasets for English [24], Italian [30] and French
[8]. The German corpus was newly created by the authors and is available at
http://spinningbytes.com/resources. An overview of the datasets used in our
experiment, including the number of labelled tweets per dataset, is given in
Table 1.

Data Preparation. Each tweet was preprocessed in three steps: (i) URLs and
usernames were substituted by a replacement token, (ii) the text was
lowercased and (iii) finally tokenized using the NLTK tokenizer.

###  4.2 Sentiment Analysis Systems

In our experiments, we compare the performance of the following sentiment
analysis systems:

  * •

Random forest (RF) as a common baseline classifier. The RF was trained on
$n$-gram features, as described in [25]

  * •

Single-language CNN (SL-CNN). The CNN with three-phase training, as described
in Section 3, is trained for each single language. In a set of experiments,
the amount of training in the three phases is gradually reduced. The system
using all available training data for one language is also referred to as
’fully-trained CNN’

  * •

Multi-language CNN (ML-CNN), where the distant-supervised phase is performed
jointly for all languages at once, and the final supervised phase
independently for each language. For the pre-training phase, we used a
balanced set of $300M$ that included all four languages, see Table 1, ’Pre-
training’

  * •

Fully multi-language CNN (FML-CNN), where all training phases were performed
without differentiation between languages. The pre-training data is the same
as in ML-CNN

  * •

SemEval benchmark. In addition, results on the English dataset were compared
to the best known ones from the SemEval
benchmark111http://alt.qcri.org/semeval2016/. For the data sets in the other
three languages, no public benchmark results could be found in literature

  * •

Translate: this approach uses Google Translate222https://translate.google.com/
(as of Oct 2016) to translate each input text from a source language to a
target language. It then uses the SL-CNN classifier trained for the target
language to classify the tweets

###  4.3 Performance Measure

We evaluate the performance of the proposed models using the metric of
SemEval-2016 challenge, which consists in averaging the macro F1-score of the
positive and negative classes333Note that this still takes into account the
prediction preformance on the neutral class.. Each approach was trained for a
fixed number of epochs. We then selected the results that yielded the best
results on a separate validation set.

For French, German and Italian, we created a validation set by randomly
sampling $10\%$ of the data. For English we used the test2015 set as
validation set and test2016 for testing from the SemEval-2016 challenge, see
Validation set in Table 1.

###  4.4 Implementation Details

The core routines of our system are written in `Theano` [3] exploiting GPU
acceleration with the `CuDNN` library [5]. The whole learning procedure takes
approximately 24-48 hours to create the word embeddings, 20 hours for the
distant-supervised phase with $160M$ tweets and only 30 minutes for the
supervised phase with $35K$ tweets.

Experiments were conducted on ’g2.2xlarge’ instances of _Amazon Web Services_
(AWS) with _GRID K520_ GPU having 3072 CUDA cores and 8 GB of RAM.

##  5 Results

In this section, we summarize the main results of our experiments.

The F1-scores of the proposed approach and competing baselines are summarized
in Table 2. The fully-trained SL-CNNs significantly outperforms the other
methods in all four languages. The best F1-score was achieved for Italian
(67.79%), followed by German (65.09%) and French (64.79%), while the system
for English reached only 62.26%. The proposed SL-CNNs outperform the
corresponding baselines from literature and RF.

Table 2: F1-scores of compared methods on the test sets. The highest scores among the three proposed models are highlighted in bold face. ML-CNN and FML-CNN are two variants of the method presented in Section 5.1. Method | Language
---|---
| English | French | German | Italian
SL-CNN | 63.49 | 64.79 | 65.09 | 67.79
ML-CNN | 61.61 | - | 63.62 | 64.73
FML-CNN | 61.03 | - | 63.19 | 64.80
RF | 48.60 | 53.86 | 52.40 | 52.71
SENSEI-LIF [24] | 62.96 | - | - | -
UNIMELB [24] | 61.67 | - | - | -

Figure 3: Results obtained by varying the amount of data during the distant
supervised phase. Each CNN was trained for one epoch.

Leveraging Distant Training Data. We increased the amount of data for the
distant-supervised phase for SL-CNN. Figure 3 compares F1-scores for each
language when changing the amount of tweets from 0 to 40M. The scores without
distant supervision are the lowest for all languages. We observe a general
increase of F1-score when increasing the amount of training data. The
performance gain for English, Italian and German is around 3%, while it is
more moderate for French.

Figure 4: Results obtained by initializing the CNNs with different word
embeddings. The fully trained variant typically performs better than the other
three variants, thus demonstrating the importance of initializing the words
vectors as well as performing distant-supervised training.

Supervised data. In Figure 5 we report the F1-scores of each model for
increasing amount of supervised data. We observe a score increase of 2-4% when
using 100% of the available data instead of 10%.

Figure 5: Results obtained by varying the amount of supervised data. The
maximum number on the y-axis corresponds to the total number of training
tweets given in Table 1.

Word Embeddings. We investigate the importance of initialization of word
embeddings and the interaction of the latter with the distant supervised phase
in four scenarios: (i) using randomly initialized word embedding weights, not
updated during the distant-supervised phase (named Full Random CNN), (ii)
using randomly initialized word embeddings, updated during the distant-
supervised phase (Random Word Embeddings), (iii) using word2vec embeddings
without distant supervision (No Distant Supervision) and (iv) using word2vec
embeddings, updated during the distant-supervised phase using $160M$ tweets
(Fully trained CNN). Results in Figure 4 demonstrate that the Fully trained
CNN approach performed the best in almost all cases. These results prove that
the quality of initialization as well as updating the large number of word
vector parameters during training of the network yield significant
improvements.

Table 3: Summary of the network parameters used for experimental results. L2 is the default architecture used by our approach. The alternative L1 and L3 architectures are discussed in Figure 8.  | Number of layers | Number of filters |  Filter window size $h$ |  Size of max-pooling window $w$ and striding $st$
---|---|---|---|---
L1 | $1$ | $300$ | $h_{1}=5$ | None
L2 | $2$ | $200$ |  $h_{1}=4$, $h_{2}=3$ |  $w_{1}=4$, $st_{1}=2$
L3 | $3$ | $200$ |  $h_{1}=4$, $h_{2}=3$, $h_{3}=2$ |  $w_{1}=4$, $st_{1}=2$, $w_{2}=3$, $st_{2}=1$
|
---|---
Figure 6: Results obtained by varying the amount of data during the distant
supervised phase. Each CNN was trained for one epoch. We rescaled the curve
for SL-CNN to match the amount of data used per language by the multi-language
approaches. For example, while the multi-language approaches were trained with
$40M$ tweets ($10M$ tweets per language), each SL-CNN model was trained with
$10M$ tweets from the same language. Each experiment set, up to $160M$ tweets,
contains the same number of tweets per language.
---
(a)
(b)
Figure 7: Word embeddings projected onto a 2-dimensional space using PCA
before (a) and after (b) the distant-supervised training phase.

Figure 8: F1-score on the English corpus as a function of the number of
network layers. These results were obtained by training the SL-CNN models with
different number of layers, as specified in Table 3, on different amounts of
data during the distant-supervised phase. Each CNN was trained for one distant
epoch.

Figure 7 illustrates the effect of the distant-supervised phase on the word
embeddings. For visualization purposes, principal component analysis (PCA) was
used to project the word embeddings onto two dimensions. We see that the
geometry of the word embeddings reflects the distance in terms of sentiment
between pairs of words. Figure 7(a) shows the initial word embeddings created
by word2vec, before the distant-supervised phase. Taking as an example the
pair of words "good" and "bad", it is clear that these two words often appear
in the same context and are thus close to each other in the embedded space.
The similarity score of these two vectors is $0.785$. After the distant-
supervised phase, the semantic of the space is changed and the distance
between words come to reflect the difference in terms of sentiment. As shown
in Figure 7(b), negative and positive words are neatly separated into two
clusters. In this case, the similarity score between the word embeddings of
"good" and "bad" becomes $-0.055$. Finer grained clusters are also revealed in
the second embedding. For example, words that convey sadness are close
together.

Comparing Network Architectures. One common question asked by practitioners
relates to the influence of the number of layers on the performance of a
neural network. We thus evaluated the performance of various architectures
with multiple layers. In order to reduce the number of experiments, we
evaluated the performance of the single-language model SL-CNN on the English
set and varied the number of convolutional/pooling layer pairs from 1 to 3. We
evaluated a total of 12 networks. Here, we only report the best set of
parameters for each number of layers in Table 3 and corresponding F1-scores in
Figure 8. The network performance generally improves with the number of
layers, if a sufficient amount of training data is used in the distant-
supervised phase. For the task of sentiment classification, current recurrent
architectures such at LSTMs still do not perform as well as CNNs, see e.g.
[29, Table 4] and the discussion in the related work section.

Table 4: Translation experiment. We translated each source language to each target language and used the model trained on the target language to predict tweets polarity. Source | Target
---|---
| English | Italian | German | French
English | 63.49 | 59.87 | 57.53 | 58.47
Italian | 64.37 | 67.79 | 61.57 | 60.19
German | 61.86 | 61.66 | 65.09 | 61.22
French | 65.68 | 63.23 | 61.5 | 64.79

Translation Approach. In Table 4 we report results of the translation
experiment described in Section 4.2. The F1-score is higher when not
translating tweets to another language for English, Italian and German. As an
exception, we obtained better results when translating French to English and
using the English model to predict sentiments.

###  5.1 Comparison to multi-language classifiers

Figure 6 summarizes F1-scores of the three CNN variants described in Section
4.2, namely SL-, ML- and FML-CNN, when varying the amount of distant-
supervised phase. When comparing the three CNN variants, we see that SL-CNN
gets slightly better scores than ML-CNN and FML-CNN. The difference in
performance between the single and multi-language models is around 2% on
average. However, one benefit of the multi-language models over the single-
language ones is their ability to deal with text in mixed languages. To check
this hypothesis, we used the langpi tool [20] to extract a set of 300 tweets
from the German corpus containing English words. Although these tweets were
classified by Twitter as German, they contain a significant number of English
words (some of them entirely written in English). We also manually inspected
this set and discarded tweets that did not contain English. We then retrained
the two models on the training set from which we first removed the set of 300
tweets. When evaluating on this subset, ML-CNN obtained an F1-score of 68.26
while SL-CNN obtained 64.07. When manually inspecting the results, we clearly
observed that ML-CNN was better at classifying tweets that were entirely in
English or contained several English words. The effect of using different word
embedding initializations in the multilingual networks is summarized in Figure
9.

Figure 9: Results obtained by initializing the CNNs with different word
embeddings. As for the results obtained with the SL-CNN model, the fully
trained variant typically performs better than the other three variants, both
for the ML-CNN and FML-CNN models.

##  6 Conclusion

We described a deep learning framework to predict the sentiment polarity of
short texts written in multiple languages. In contrast to most existing
methods, our approach does not rely on establishing a correspondence to
English but instead exploits large amounts of weakly supervised data to train
a multi-layer CNN directly in the target language. Through a thorough
experimental evaluation, we addressed some fundamental questions around the
performance of such model. First, we demonstrated that the strategy used to
train these models plays an important role in the obtained performance. Two
important factors are a good initialization for the word vectors as well as
pre-training using large amounts of weakly supervised data. Second, we
compared the performance of a single-language and a multi-language approach.
The single-language model reaches the best performance and it even outperforms
existing state-of-the-art methods on all the datasets of the SemEval-2016
competition. The multi-language approach performs comparably well or slightly
worse than its single-language counterpart, while exhibiting several
advantages: it does not need to know a priori the language(s) used in each
tweet; the model can be easily extended to more languages; and it can cope
with texts written in multiple languages.

#### Acknowledgments

This research has been funded by Commission for Technology and Innovation
(CTI) project no. 18832.1 PFES-ES, and by SpinningBytes AG, Switzerland.

## References

  * [1] M. Araujo, J. Reis, A. Pereira, and F. Benevenuto.  An evaluation of machine translation for multilingual sentence-level sentiment analysis.  In Proceedings of the 31st Annual ACM Symposium on Applied Computing, pages 1140–1145, 2016.
  * [2] A. Balahur and M. Turchi.  Multilingual sentiment analysis using machine translation?  In Proceedings of the 3rd Workshop in Computational Approaches to Subjectivity and Sentiment Analysis, 2012.
  * [3] J. Bergstra, O. Breuleux, F. F. Bastien, P. Lamblin, R. Pascanu, G. Desjardins, J. Turian, D. Warde-Farley, and Y. Bengio.  Theano: a CPU and GPU math compiler in Python.  Proceedings of the Python for Scientific Computing Conference (SciPy), pages 1–7, 2010.
  * [4] E. Boiy and M.-F. Moens.  A machine learning approach to sentiment analysis in multilingual web texts.  Information retrieval, 12(5):526–558, 2009.
  * [5] S. Chetlur and C. Woolley.  cuDNN: Efficient Primitives for Deep Learning.  arXiv preprint, pages 1–9, 2014.
  * [6] M. Cieliebak, O. Dürr, and F. Uzdilli.  Potential and limitations of commercial sentiment detection tools.  In ESSEM@ AI* IA, pages 47–58, 2013.
  * [7] K. Dashtipour, S. Poria, A. Hussain, E. Cambria, A. Y. A. Hawalah, A. Gelbukh, and Q. Zhou.  Multilingual sentiment analysis: State of the art and independent comparison of techniques.  Cognitive Computation, 2016.
  * [8] DEFT.  https://deft.limsi.fr/2015/, 2015.
  * [9] J. Deriu, M. Gonzenbach, F. Uzdilli, A. Lucchi, V. De Luca, and M. Jaggi.  Swisscheese at semeval-2016 task 4: Sentiment classification using an ensemble of convolutional neural networks with distant supervision.  Proceedings of SemEval, pages 1124–1128, 2016.
  * [10] C. N. dos Santos and M. Gatti.  Deep Convolutional Neural Networks for Sentiment Analysis of Short Texts.  In COLING, pages 69–78, 2014.
  * [11] D. Erhan, Y. Bengio, A. Courville, P.-A. Manzagol, P. Vincent, and S. Bengio.  Why does unsupervised pre-training help deep learning?  The Journal of Machine Learning Research, 11:625–660, 2010.
  * [12] J. Gao, X. He, W.-t. Yih, and L. Deng.  Learning continuous phrase representations for translation modeling.  In ACL, 2014.
  * [13] A. Go, R. Bhayani, and L. Huang.  Twitter Sentiment Classification using Distant Supervision.  Technical report, The Stanford Natural Language Processing Group, 2009\.
  * [14] S. Gouws, Y. Bengio, and G. Corrado.  Bilbowa: Fast bilingual distributed representations without word alignments.  In Proceedings of the 32nd International Conference on Machine Learning, 2015.
  * [15] IWS.  http://www.internetworldstats.com/stats7.htm, 2016.
  * [16] R. Johnson and T. Zhang.  Semi-supervised Convolutional Neural Networks for Text Categorization via Region Embedding.  In Advances in Neural Information Processing Systems 28, pages 919–927, 2015.
  * [17] N. Kalchbrenner, E. Grefenstette, and P. Blunsom.  A Convolutional Neural Network for Modelling Sentences.  In ACL - Proceedings of the 52nd Annual Meeting of the Association for Computational Linguistics, pages 655–665, 2014.
  * [18] Y. Kim.  Convolutional Neural Networks for Sentence Classification.  In EMNLP 2014 - Empirical Methods in Natural Language Processing, pages 1746–1751, 2014.
  * [19] M. Lui and T. Baldwin.  Cross-domain feature selection for language identification.  In In Proceedings of 5th International Joint Conference on Natural Language Processing, 2011.
  * [20] M. Lui and T. Baldwin.  langid. py: An off-the-shelf language identification tool.  In Proceedings of the ACL 2012 system demonstrations, pages 25–30. Association for Computational Linguistics, 2012.
  * [21] Mashable.  http://mashable.com/2013/12/17/twitter-popular-languages/#gaFEjWnHPkql, 2013\.
  * [22] R. Mihalcea, C. Banea, and J. M. Wiebe.  Learning multilingual subjective language via cross-lingual projections.  In ACL, 2007.
  * [23] V. Nair and G. E. Hinton.  Rectified linear units improve restricted boltzmann machines.  In Proceedings of the 27th International Conference on Machine Learning, pages 807–814, 2010.
  * [24] P. Nakov, A. Ritter, S. Rosenthal, V. Stoyanov, and F. Sebastiani.  SemEval-2016 task 4: Sentiment analysis in Twitter.  In Proceedings of the 10th International Workshop on Semantic Evaluation, SemEval ’16, June 2016.
  * [25] S. Narr, M. Hulfenhaus, and S. Albayrak.  Language-independent twitter sentiment analysis.  Knowledge Discovery and Machine Learning (KDML), LWA, pages 12–14, 2012.
  * [26] J. Read.  Using emoticons to reduce dependency in machine learning techniques for sentiment classification.  In Proceedings of the ACL student research workshop, pages 43–48. Association for Computational Linguistics, 2005.
  * [27] F. N. Ribeiro, M. Araújo, P. Gonçalves, F. Benevenuto, and M. A. Gonçalves.  A benchmark comparison of state-of-the-practice sentiment analysis methods.  arXiv preprint arXiv:1512.01818, 2015.
  * [28] S. Rothe, S. Ebert, and H. Schutze.  Ultradense Word Embeddings by Orthogonal Transformation.  arXiv, 2016.
  * [29] S. Semeniuta, A. Severyn, and E. Barth.  Recurrent Dropout without Memory Loss.  arXiv, 2016.
  * [30] sentipolc.  http://www.di.unito.it/~tutreeb/sentipolc-evalita16/data.html, 2016\.
  * [31] A. Severyn and A. Moschitti.  Twitter Sentiment Analysis with Deep Convolutional Neural Networks.  In 38th International ACM SIGIR Conference, pages 959–962, 2015\.
  * [32] A. Severyn and A. Moschitti.  UNITN: Training Deep Convolutional Neural Network for Twitter Sentiment Classification.  In SemEval 2015 - Proceedings of the 9th International Workshop on Semantic Evaluation, 2015.
  * [33] Y. Shen, X. He, J. Gao, L. Deng, and G. Mesnil.  A latent semantic model with convolutional-pooling structure for information retrieval.  In Proceedings of the 23rd ACM International Conference on Conference on Information and Knowledge Management, pages 101–110, 2014.
  * [34] R. Socher, A. Perelygin, J. Y. Wu, J. Chuang, C. D. Manning, A. Y. Ng, and C. Potts.  Recursive deep models for semantic compositionality over a sentiment treebank.  In Proceedings of the conference on empirical methods in natural language processing, volume 1631, page 1642, 2013.
  * [35] M. Wick, P. Kanani, and A. Pocock.  Minimally-constrained multilingual embeddings via artificial code-switching.  2015\.
  * [36] M. D. Zeiler.  ADADELTA: An Adaptive Learning Rate Method.  arXiv, 2012.


# Interplay between Kitaev interaction and single ion anisotropy in
ferromagnetic CrI3 and CrGeTe3 monolayers

Changsong Xu1 Junsheng Feng2,3 Hongjun Xiang2,4,∗ & L. Bellaiche1,∗

###### Abstract

Magnetic anisotropy is crucially important for the stabilization of two-
dimensional (2D) magnetism, which is rare in nature but highly desirable in
spintronics and for advancing fundamental knowledge. Recent works on CrI3 and
CrGeTe3 monolayers not only led to observations of the long-time-sought 2D
ferromagnetism, but also revealed distinct magnetic anisotropy in the two
systems, namely Ising behavior for CrI3 versus Heisenberg behavior for
CrGeTe3. Such magnetic difference strongly contrasts with structural and
electronic similarities of these two materials, and understanding it at a
microscopic scale should be of large benefits. Here, first-principles
calculations are performed and analyzed to develop a simple Hamiltonian, to
investigate magnetic anisotropy of CrI3 and CrGeTe3 monolayers. The
anisotropic exchange coupling in both systems is surprisingly determined to be
of Kitaev-type. Moreover, the interplay between this Kitaev interaction and
single ion anisotropy (SIA) is found to naturally explain the different
magnetic behaviors of CrI3 and CrGeTe3. Finally, both the Kitaev interaction
and SIA are further found to be induced by spin-orbit coupling of the heavy
ligands (I of CrI3 or Te of CrGeTe3) rather than the commonly believed 3d
magnetic Cr ions.

Physics Department and Institute for Nanoscience and Engineering, University
of Arkansas, Fayetteville, Arkansas 72701, USA

Key Laboratory of Computational Physical Sciences (Ministry of Education),
State Key Laboratory of Surface Physics, and Department of Physics, Fudan
University, Shanghai, 200433, China

School of Electronic and Information Engineering, Hefei Normal University,
Hefei 230601, P. R. China

Collaborative Innovation Center of Advanced Microstructures, Nanjing 210093,
P. R. China

INTRODUCTION
Two-dimensional (2D) magnetic materials are receiving a lot of attention, due,
e.g., to the search for long range ferromagnetism (FM)[1, 2], which can
facilitate various applications from sensing to data storage[3, 4]. According
to Mermin and Wagner s theorem[2], however strong the short-range isotropic
couplings are, the realization of 2D magnetism relies on magnetic anisotropy,
as a result of spin-orbit coupling (SOC). The requirement of strong magnetic
anisotropy in low dimensional systems therefore explains the rareness of 2D FM
materials.

The recent observation of ferromagnetism in monolayers made of CrI3 and
CrGeTe3 [5, 6, 7] therefore opens a new chapter in the field of 2D materials.
The chromium in both compounds share the same valence state of Cr3+, with the
3$d^{3}$ configuration and S = $\frac{3}{2}$[8, 9, 7, 10]. Ferromagnetism
arises there from the super exchange between nearest neighbor Cr ions, that
are linked by I or Te ligands through nearly 90∘ angles[10, 11]. CrI3 has been
demonstrated to be well described by the Ising behavior[5, 12, 1], for which
the spins can point up and down along the out-of-plane z-direction. In
contrast, the magnetic anisotropy of CrGeTe3 was determined to be consistent
with the Heisenberg behavior[12, 7, 2], for which the spins can freely rotate
and adopt any direction in the three-dimensional space. Interestingly,
structural and electronic similarities between these two compounds strongly
contrast with their difference in magnetic behaviors, which implies subtle
origins for their magnetic anisotropy. A recent theoretical work adopted the
XXZ model, for which the exchange coupling is identical between the in-plane
x- and y-directions but different along the out-of-plane z-direction, to
explain the out-of-plane magnetization of CrI3[8]. However, there is no
definite proof that the XXZ model is accurate enough to describe the magnetic
anisotropy of CrI3, and there is a current paucity of knowledge for the
mechanism responsible for the magnetic anisotropy of CrGeTe3. Hence, a
thorough microscopic understanding of the difference between the Ising
behavior of CrI3 and the Heisenberg behavior of CrGeTe3 is highly desired.

In particular, it is tempting to investigate if the Kitaev interaction[13],
which is a specific anisotropic exchange coupling, can also be significant in
CrI3 and CrGeTe3. This temptation is mainly based on the fact that these two
materials adopt a honeycomb lattice and edge-sharing octahedra, exactly as the
layered Na2IrO3 and $\alpha$-RuCl3 compounds which exhibit magnetic behaviors
that are close to spin liquids[14] –as a result of significant Kitaev
interactions. Interestingly, finding finite Kitaev interaction in Cr-3$d$
based CrI3 and CrGeTe3 compounds would enlarge the types of systems possessing
such interaction, not only from $4d$ or $5d$ to $3d$ transition-metal-based
insulators, but also from S = 1/2 to S = 3/2 systems. Such broadening is in-
line with recent theoretical predictions of Kitaev interaction in $d^{7}$ or
$3d$ systems[15, 16].

The main goal of this manuscript is to report results of first-principles
calculations, along with the concomitant development of a simple but
predictive Hamiltonian, to demonstrate that: (i) significant Kitaev
interaction does exist in CrI3 and CrGeTe3 (which also invalidates the XXZ
model in these two compounds); and (ii) the different interplay between this
Kitaev interaction and the single ion anisotropy (SIA) naturally explains the
observed magnetic anisotropy in these 2D ferromagnetic materials. Another
surprising result is that the Kitaev and SIA anisotropies are both dominantly
induced by the SOC of the heavy ligand elements rather than the $3d$ element
Cr.

RESULTS
To precisely describe the magnetic anisotropy and explore differences between
CrI3 and CrGeTe3, we consider a Hamiltonian containing both exchange coupling,
$\mathcal{H}_{ex}$, and SIA, $\mathcal{H}_{si}$, terms:

$\mathcal{H}=\mathcal{H}_{ex}+\mathcal{H}_{si}=\frac{1}{2}\sum_{i,j}\bm{{\rm
S}}_{i}{\cdot}\mathcal{J}_{ij}{\cdot}\bm{{\rm S}}_{j}+\sum_{i}\bm{{\rm
S}}_{i}{\cdot}\mathcal{A}_{ii}{\cdot}\bm{{\rm S}}_{i}$ (1)

where $\mathcal{J}_{ij}$ and $\mathcal{A}_{ii}$ are 3$\times$3 matrices
gathering exchange and SIA parameters, respectively. The sum over $i$ in Eq.
(1) runs over all Cr sites, while the sums over $i,j$ run over all nearest
neighbor Cr pairs(note that the anisotropy in exchange coupling between more
distant Cr neighbors is at least an order smaller and is thus negligible, see
details in Supplementary Materials (SM)[17]). DFT calculations are performed
on CrI3 and CrGeTe3 monolayers to extract the components of $\mathcal{J}$ and
$\mathcal{A}$ using a precise four-states method[18] (see SM for details[17]).
Note that (i) all the results shown below are based on the use of an effective
Hubbard $U$ = 0.5 eV, unless stated (see the effects of the choice of other
$U$’s and details of method in SM[17]); and (ii) the Dzyaloshinski-Moriya (DM)
interaction is absent in our studied systems because of the existence of an
inversion center between nearest neighbor Cr ions[19].

Figure 1: Schematization of the CrI3 and CrGeTe3 structures, as well as the
different coordinate systems indicated in the text. Note that Ge of CrGeTe3 is
not shown for simplicity.

Let us first focus on the exchange coupling $\mathcal{J}$ matrix for CrI3 and
CrGeTe3. The $\mathcal{J}$ matrix is expressed in the {xyz} basis, for which
the x-y plane is the film plane while the z-axis is the out-of-plane direction
of the film. We choose the Cr0-Cr1 pair (see Fig. 1) to calculate the exchange
coupling parameters, from which the parameters for Cr0-Cr2 and Cr0-Cr3 pairs
can be deduced via three-fold rotational symmetry. It is numerically found
that this matrix is symmetric, i.e., $J_{xy}$=$J_{yx}$, $J_{yz}$=$J_{zy}$ and
$J_{xz}$=$J_{zx}$, which is consistent with the fact that there is no DM
interaction in our investigated compounds. As shown in Table 1, $J_{xx}$,
$J_{yy}$ and $J_{zz}$ of CrI3 possess quite different values of -2.29, -1.93
and -2.23 meV, respectively, while the off-diagonal elements of $\mathcal{J}$
in the {xyz} basis are smaller but non-negligible. Such results contrast with
the XXZ model adopted in Ref.[8] and that assumes that (i) $J_{xx}=J_{yy}\neq
J_{zz}$; and (ii) $J_{xy}$, $J_{yz}$ and $J_{xz}$ can all be neglected.
Different schemes and strategies, such as changing the value of $U$, using
experimental structures rather than the computationally optimized ones and
even replacing other Cr ions by nonmagnetic Al, are used to check their
influence on $J_{xx}$, $J_{yy}$ and $J_{zz}$ of CrI3. It is numerically found
that they all qualitatively give the same results (as detailed in SM[17]) in
the sense that the aforementioned assumption (i) of the XXZ model providing
equality between $J_{xx}$ and $J_{yy}$ is not satisfied, which automatically
implies that such latter model is not accurate enough to precisely describe
magnetic anisotropy in CrI3 and CrGeTe3 systems.

Table 1: Matrix components of $\mathcal{J}$, as well as the $J$ and $K$ parameters, of the Cr0-Cr1 pair, and the SIA coefficient $A_{zz}$ for the Cr0 ion, as given by DFT calculations with $S$ = $\frac{3}{2}$. Corresponding coefficients for the Cr0-Cr2 and Cr0-Cr3 pairs can be deduced via three-fold rotational symmetry. Note that the global coordinate system {xyz} and the local coordinate system {$\alpha\beta\gamma$} of each Cr-Cr pair are shown in Fig. 1. The units of the parameters indicated here is meV. |  $J_{xx}$ |  $J_{yy}$ |  $J_{zz}$ |  $J_{xy}$ |  $J_{yz}$ |  $J_{xz}$
---|---|---|---|---|---|---
CrI3 |  -2.29 |  -1.93 |  -2.23 |  0.30 |  0.29 |  0.17
CrGeTe3 |  -6.54 |  -6.37 |  -6.64 |  0.15 |  -0.04 |  -0.02
|  $J_{\alpha}$ |  $J_{\beta}$ |  $J_{\gamma}$ |  $J$ |  $K$ |  $A_{zz}$
CrI3 |  -2.46 |  -2.41 |  -1.59 |  -2.44 |  0.85 |  -0.26
CrGeTe3 |  -6.65 |  -6.63 |  -6.28 |  -6.64 |  0.36 |  0.25
Table 2: Parameters of Eq. (5) and $\Delta{\mathcal{E}}$ from both Eq. (5) and DFT for CrI3 and CrGeTe3. $\Delta{\mathcal{E}}$ defines the energy difference between the energy for an out-of-plane magnetization and the averaged energy of ferromagnetic states having in-plane magnetization.  | Parameters (meV) |  |  $\Delta{\mathcal{E}}$ (meV/f.u.)
---|---|---|---
$bK$ |  $\frac{2}{3}A_{zz}$ |  $bK$+$\frac{2}{3}A_{zz}$ |  |  Eq.(5) |  DFT
CrI3 |  -0.16 |  -0.17 |  -0.33 |  |  -1.11 |  -0.82
CrGeTe3 |  -0.17 |  0.17 |  -0.003 |  |  -0.01 |  0.02

The symmetric $\mathcal{J}$ matrix is then diagonalized to obtain its
eigenvalues (to be denoted as $J_{\alpha}$, $J_{\beta}$ and $J_{\gamma}$) and
corresponding eigenvectors (to be coined $\alpha$, $\beta$ and $\gamma$) for
the Cr0-Cr1 pair. As seen in Table 2, $J_{\alpha}$ and $J_{\beta}$ are the
strongest eigenvalues in magnitude and are close to each other in CrI3 (-2.46
meV and -2.41 meV, respectively), while $J_{\gamma}$ is smaller in magnitude
by about 1 meV. The same hierarchy exists between $J_{\alpha}$, $J_{\beta}$
and $J_{\gamma}$ in CrGeTe3, but with $J_{\alpha}$ and $J_{\beta}$ being now
stronger in strength (about -6.65 meV), while $J_{\gamma}$ = -6.28 meV is less
than 0.4 meV smaller in magnitude than the other two exchange coefficients. As
shown in Fig. 1, the $\alpha$-axis points from Cr1 to Cr0 in both systems, and
therefore belongs to the x-y plane. On the other hand, the $\beta$-axis is
roughly along the direction joining the two ligands bridging Cr0 and Cr1, and
thus does not belong to the x-y plane. Similarly, the $\gamma$-axis, which is
perpendicular to both the $\alpha$\- and $\beta$-axes, does not lie in the x-y
plane. It is important to realize that the {$\alpha$$\beta$$\gamma$} basis
diagonalizing the $\mathcal{J}$ matrix is specific to each considered Cr
pairs, unlike the “global” {xyz} basis. In other words, the
{$\alpha$$\beta$$\gamma$} basis differs for the Cr0-Cr1, Cr0-Cr2 and Cr0-Cr3
pairs, as shown in Fig. 1 by means of red, green and blue arrows.

The $\mathcal{H}_{ex}$ exchange coupling Hamiltonian can now be rewritten,
gathering these three local orthogonal {$\alpha\beta\gamma$} coordinate bases
(one for each Cr pair) and assuming that $J_{\alpha}$ = $J_{\beta}$ (as
consistent with the aforementioned calculations[20]), as

$\displaystyle\mathcal{H}_{ex}=$
$\displaystyle\frac{1}{2}\sum_{i,j}(J_{\alpha}S_{i}^{\alpha}S_{j}^{\alpha}+J_{\beta}S_{i}^{\beta}S_{j}^{\beta}+J_{\gamma}S_{i}^{\gamma}S_{j}^{\gamma})$
(2) $\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{i,j}(J\bm{{\rm
S}}_{i}{\cdot}\bm{{\rm S}}_{j}+KS_{i}^{\gamma}S_{j}^{\gamma})$

where $J$ = $J_{\alpha}$ = $J_{\beta}$ is the isotropic exchange coupling and
$K$ = $J_{\gamma}-J_{\alpha}>0$ is the so-called Kitaev interaction that
characterizes the anisotropic contribution. Table 1 provides the values of
both $J$ and $K$ and, in particular, indicates that the Kitaev interaction can
not be neglected in CrI3 and CrGeTe3.

Let us now investigate the other energy of Eq. (1), that is the SIA, which
involves the $\mathcal{A}$ matrix. For that, one needs to go back to the
global {xyz} basis, since only the $A_{zz}$ term can be finite by symmetry. It
is numerically found that $A_{zz}$ = -0.26 meV in CrI3, while it adopts a
similar magnitude but with a change of sign in CrGeTe3 (since $A_{zz}$ = 0.25
meV there). Such significant values of $A_{zz}$ (which is of the same order of
magnitude than the $K$ Kitaev parameter) implies that SIA is not negligible,
which contrasts with the results in Ref.[8]. $\mathcal{H}_{si}$ can thus be
simplified as:

$\mathcal{H}_{si}=\sum_{i}A_{zz}{\rm S}^{z}_{i}{\rm S}^{z}_{i}$ (3)

The total Hamiltonian of Eq. (1) can then be rewritten by combining Eqs. (2)
and (3) as:

$\mathcal{H}=\frac{1}{2}\sum_{i,j}(J\bm{{\rm S}}_{i}{\cdot}\bm{{\rm
S}}_{j}+KS_{i}^{\gamma}S_{j}^{\gamma})+\sum_{i}A_{zz}{\rm S}^{z}_{i}{\rm
S}^{z}_{i}$ (4)

This simplified Hamiltonian gathers (i) isotropic exchange coupling from $J$;
(ii) anisotropic Kitaev interaction from $K$ in the different local
{$\alpha\beta\gamma$} bases; and (iii) SIA in the global {xyz} basis. Let us
now try to express the total energy associated with magnetism in an unified
coordinate system. Equation (4) shows that the anisotropic part of the
exchange energy (arising from $K$) is only related to the projections of spins
on the three different $\gamma$-axes (one for each Cr-Cr pair). Due to the
fact that these three $\gamma$-axes (to be denoted as $\gamma_{1}$,
$\gamma_{2}$ and $\gamma_{3}$, respectively) are normally not perpendicular to
each other, we now orthogonalize them using the Löwdin’s symmetric
orthogonalization scheme[21]. The resulting orthogonal axes form the {XYZ}
coordinate system that is shown in Fig. 1. In this global {XYZ} basis, the
out-of-plane z-axis of the film is along the [111] direction, and
$\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ can be expressed as (1,a,a),
(a,1,a) and (a,a,1), where a$\in$[-$\frac{1}{2}$, 1]. As consistent with the
relatively large-in-magnitude and negative value of the isotropic exchange
coupling $J$ (see Table 1), ferromagnetic states are considered here. When
expressing their spin in the {XYZ} basis, i.e., $\bm{{\rm S}}$($S_{X}$,
$S_{Y}$, $S_{Z}$), and considering a magnitude $|\bm{{\rm S}}|=\frac{3}{2}$,
it is straightforward to prove that the energy per Cr ion associated with Eq.
(4) can then be rewritten as:

$\displaystyle\mathcal{E}=(bK+\frac{2}{3}A_{zz})(S_{X}S_{Y}+S_{Y}S_{Z}+S_{Z}S_{X})+{\rm
C}$ (5)

where $b=\frac{a^{2}+2a}{2a^{2}+1}$ and C =
$\frac{9}{8}(3J+K)+\frac{3}{4}A_{zz}$ are independent of the spin direction
(see SM[17] for details). One can also easily demonstrate that the symmetric
form of $(S_{X}S_{Y}+S_{Y}S_{Z}+S_{Z}S_{X})$ implies that the magnetization
within the x-y plane of the film (for which $S_{X}+S_{Y}+S_{Z}=0$) is fully
isotropic. In other words, any direction of the spin within this x-y plane
generates the same energy. Note that we further conducted DFT calculations
(not shown here) that indeed numerically confirm that such in-plane isotropy
is mainly obeyed in CrI3 and CrGeTe3 (the maximal energetic difference we
found between in-plane directions of spins is 0.006 meV/f.u. and 0.004
meV/f.u. in CrI3 and CrGeTe3, respectively), which attests of the relevance
and accuracy of the simple Hamiltonian of Eq. (4) and the resulting energy of
Eq. (5). Note also that, although such isotropy in the x-y plane is in line
with the results of Ref.[8], its origin is totally different: here it arises
from the Kitaev interaction and its subsequent frustration, while in Ref.[8],
the isotropy in the x-y plane lies in the assumption of the XXZ model (see a
detailed comparison between the two models in SM[17])

It is also worthwhile to emphasize that Ref.[8] assumed that SIA is negligible
small, while according to Eq. (5) and as we will show below, both the Kitaev
interaction ($K$) and SIA ($A_{zz}$) play an important role on the overall
magnetic anisotropy of CrI3 and CrGeTe3. To demonstrate such fact, one can
realize that Eq. (5) involves $(S_{X}S_{Y}+S_{Y}S_{Z}+S_{Z}S_{X})$, which
adopts (i) its maximum when $S_{X}=S_{Y}=S_{Z}=\frac{\sqrt{3}}{2}$, which
corresponds to spins being aligned along the out-of-plane z-direction; versus
(ii) a minimum when $S_{X}+S_{Y}+S_{Z}=0$, that is when spins are lying within
the x-y plane. The sign and value of the $bK+\frac{2}{3}A_{zz}$ coefficient
appearing in front of $(S_{X}S_{Y}+S_{Y}S_{Z}+S_{Z}S_{X})$ in Eq. (5) should
therefore determine the magnetic anisotropy: a negative $bK+\frac{2}{3}A_{zz}$
favors an easy axis along the out-of-plane direction while a positive
$bK+\frac{2}{3}A_{zz}$ will encourage spins to lie within the x-y plane. To
characterize the strength of such anisotropy between the out-of-plane
direction and the x-y plane, we also computed the energy difference,
$\Delta{\mathcal{E}}$, between the energy of the state having a fully out-of-
plane magnetization and the averaged energy of states having in-plane
magnetization.

In the case of CrI3, the $b$ parameter is numerically found to be negative.
Together with the positive $K$ and negative $A_{zz}$ from Table 1, both $bK$
and $\frac{2}{3}A_{zz}$ are thus negative. They are determined to be -0.16 meV
and -0.17 meV, respectively, as shown in Table 2. Such negative values
indicate that both Kitaev interaction and SIA lead to an out-of-plane easy
axis, which is further confirmed by the negative value of -1.11 meV/f.u for
$\Delta{\mathcal{E}}$, as calculated from Eq. (5). Such value is not only
consistent with the result of -0.82 meV/f.u. obtained from DFT calculations
(confirming once again the validity and accuracy of our rather simple Eqs. (4)
and (5)) but also explains the previously determined Ising behavior of CrI3
favoring the out-of-plane direction for the magnetization[5, 12, 1].

In the case of CrGeTe3, the $b$ parameter is also found to be negative and
leads to $bK$ adopting a negative $-0.17$ meV value that is similar to the one
of CrI3. On the other hand, $\frac{2}{3}A_{zz}$ is positive and yields
$\simeq$ 0.17 meV, which therefore results in a nearly vanishing
$bK$+$\frac{2}{3}A_{zz}$ and thus to a $\Delta{\mathcal{E}}$ being nearly zero
– that is, -0.01 meV/f.u. according to Eq. (5), which also compares well with
the result of 0.02 meV/f.u. directly obtained from DFT (the difference in sign
between the energies from Eq. (5) and DFT can be overlooked since both
calculations provide vanishing $\Delta{\mathcal{E}}$). Such nearly zero value
for $bK$+$\frac{2}{3}A_{zz}$ therefore implies that Eq. (5) predicts that
CrGeTe3 is basically isotropic in the whole three-dimensional space, that is
any spatial direction of the magnetization (in-plane, out-of-plane or even
combination between in-plane and out-of-plane components) should provide
similar magnetic energy. Such finding is fully consistent with the observed
Heisenberg behavior of CrGeTe3[12, 7, 2]. Note that, although the isotropic
Heisenberg behavior of CrGeTe3 is confirmed by both previous experiments and
the present computational work, the anisotropy (Kitaev interaction and SIA)
still play a crucial important role in stabilization of the long-range
magnetic ground state, as indicated by Mermin and Wagner s theorem[2].

Figure 2: Atomic dependencies of the Kitaev parameter $K$ and SIA coefficient
$A_{zz}$. Panels a and b show the dependency of $K$ of CrI3 and CrGeTe3,
respectively, as a function of the SOC strength. Panels c and d display the
dependency of $A_{zz}$ for CrI3 and CrGeTe3, respectively, as a function of
this SOC strength. Note that the value of 1.0 (0, respectively) for this SOC
strength corresponds to the actual strength (no SOC, respectively) of the
considered element.

It is also worthwhile to realize that both the Kitaev parameter $K$ and the
SIA coefficient $A_{zz}$ originate from spin-orbit coupling. One may wonder
what specific ions contribute to such coefficients via SOC. To address such
issue, Figs. 2a and b display the atomically resolved contribution of the $K$
parameter as a function of the SOC strength in CrI3 and CrGeTe3, respectively.
Such parameter mainly arises from the SOC of the heavy ligands (namely, I of
CrI3 or Te of CrGeTe3) in a quadratic way, while the SOC of Cr has almost no
effect on $K$ (the effect from SOC of Ge on $K$ of CrGeTe3 is also
negligible). Such predictions for CrI3 and CrGeTe3 are consistent with the
results of Ref.[8] that anisotropy of CrI3 mainly arises from the SOC of I
ligand. This can be understood by the fact that the super exchange between
nearest neighbor Cr sites is mostly mediated by these ligands. To further
confirm such results, we developed a tight-binding model (see details in
SM[17]) that contains only two Cr ions and two bridging ligands, which form
the $x^{\prime}y^{\prime}$ plane (see Fig. 3a for details). Such model
confirms that $K>$ 0 and $K\varpropto\lambda^{2}$, where $\lambda$ is the SOC
strength of ligands.

Figure 3: Panel a, the {x’y’z’} coordinate system and the configuration of the
considered Cr2L2 cluster in the tight binding model, where L is a ligand ion.
Panel b, schematization of the extra hopping path related to spins being along
the $z^{\prime}$ ($\uparrow$) direction. Panel c, the forbidden hopping paths
related to spins lying in the $x^{\prime}y^{\prime}$ (–) plane, respectively.

We now work on understanding the analytical results that $K>0$ and
$K\varpropto\lambda^{2}$. If SOC is not considered, the magnetic coupling is
isotropic with the strength of $J$. When the SOC of ligands is included, an
extra hopping path emerges for spins that are along the $z^{\prime}$
direction, as shown in Figs. 3a and b, since $\langle
p_{x^{\prime}\uparrow}|\mathbf{L}\cdot\mathbf{S}|p_{y^{\prime}\uparrow}\rangle\neq
0$. Such an extra path provides an additional energy term $K$ to the magnetic
coupling, as
$J_{z^{\prime}z^{\prime}}=J+K=\dfrac{1}{2S^{2}}(E_{FM,z^{\prime}}-E_{AFM,z^{\prime}})$.

In the antiferromagnetic case, the two Cr have opposite spins. The electrons
can hop from the occupied $t_{2g}$ orbitals of one Cr to the unoccupied
$t_{2g}$ orbitals of another Cr, which can lower the $E_{AFM,z^{\prime}}$ by
the amount of $K\varpropto\lambda^{2}$. The form of $\lambda^{2}$ can be
understood as that the whole hopping procedure includes two times of ligands’
SOC effects, as shown in Fig. 3b. On the other hand, in the ferromagnetic
case, $t_{2g}$ orbitals of both Cr are occupied with electrons having the same
spin direction. In such case, although the aforementioned extra hopping path
still exists, it can not lead to the energy lowering of $E_{FM,z^{\prime}}$.
In contrast, there is no such extra hopping path when spins lie in the Cr2L2
($x^{\prime}y^{\prime}$) plane (denoted with the “–” mark), since $\langle
p_{x^{\prime}-}|\mathbf{L}\cdot\mathbf{S}|p_{y^{\prime}-}\rangle=0$. Such
effects are further illustrated in Fig. 3c. As a result, the total effective
$J_{z^{\prime}z^{\prime}}=J+K$ is larger than that
$J_{x^{\prime}x^{\prime}}=J_{y^{\prime}y^{\prime}}=J$, i.e., $K>$ 0\. The
findings reported here thus demonstrate, for the first time, that (i) the
Kitaev interaction not only exists in 4$d$ or 5$d$ transition metal
insulators, but also can occur in 3$d$ systems; and (ii) the SOC of the
ligands can play a crucial role on that interaction. Such findings can
facilitate the ongoing efforts to realize Kitaev-type interactions in 3$d$
systems[15, 16].

Regarding SIA, Fig. 2c shows that the SOC of I is basically responsible for
the negative $A_{zz}$ of CrI3, and that Cr does not significantly contribute
to such parameter. Such finding of a ligand-induced SIA is also novel, since
SIA is typically believed to arise from the transition metal ion[22, 23].
Similarly and as indicated by Fig. 2d, the Te ligand in CrGeTe3 produces a
rather strong (and positive, in that case) $A_{zz}$, which is even twice as
large as the one resulting from the SOC of all considered ions. It is in fact
the combination of SOC from both Te and Cr that provides a value of $A_{zz}$
that is close to the total one in CrGeTe3 (note that Ge does not
significantely contribute to this total SIA).

DISCUSSION
Explanations on magnetic behaviors of related 2D systems. Another 2D
ferromagnetic system, CrBr3, was also studied. It is numerically found (not
shown here) that, as similar to CrI3, both the Kitaev interaction and SIA
favor an out-of-plane easy axis in this material. Such finding is consistent
with measurements determining that the net anisotropy of CrBr3 is out-of-
plane[24]. Interestingly, Eq. (5) of the manuscript can also be useful to shed
some light into controversial issues, such as which magnetic model is more
pertinent to CrSiTe3. As a matter of fact, early neutron work suggested an
Ising-like model[25] for this system, while recent measurements argue between
an Heisenberg-like model[26] and an Ising model coupled with long-range
interaction[27]. CrSiTe3 shares similarities with CrGeTe3 in the sense that it
has a negative $bK$ associated with Kitaev interaction, as well as a positive
$\frac{2}{3}A_{zz}$ induced by SIA. However, the former is equal to $-$0.21
meV and is larger in magnitude than the latter (that is equal to +0.11 meV) in
CrSiTe3, therefore leading to a less negative $bK$+$\frac{2}{3}A_{zz}$ of
$-$0.07 meV and thus slightly tipping the balance towards out-of-plane
magnetization (as confirmed by our DFT results providing
$\Delta{\mathcal{E}}$=-0.10 meV for the difference between the energy for an
out-of-plane magnetization and the averaged energy of ferromagnetic states
having in-plane magnetization). One can thus propose that the correct magnetic
model for CrSiTe3 should be related to a slight perturbation of the Ising
model, in order to account for possible (weaker) in-plane components of the
magnetization in addition to (stronger) out-of-plane ones.

Applications of the present general model and the XXZ model. The present
general model (Eqs. 1, 4 and 5) adopts the most generalized form of the
$\mathcal{J}$ and $\mathcal{A}$ matrices, which can capture the microscopic
details of different anisotropy. Such model is powerful, as it explains the
origin of Kitaev interaction, as well as the competition and collaboration
between Kitaev interaction and SIA. It also allows for antisymmetric exchange
coupling, i.e. the so called Dzyaloshinskii-Moriya interaction, when the
inversion centers between Cr-Cr pair are somehow removed. In contrast, the XXZ
model is more macroscopic in nature, since it starts from the overall effects
of the frustration among Cr-Cr pairs. Nevertheless, this XXZ model can still
somehow describe the competition and collaboration between Kitaev interaction
(since $\mathcal{J}_{x}$ and $\mathcal{J}_{z}$ are different from each other)
and SIA. As a result, the XXZ model can be technically applied to both CrI3
and CrGeTe3, as well as to the aforementioned related systems, but, as
documented in section 4 of the Supplemental Materials, one really has to
include SIA there to be more accurate.

Implication of potential Kitaev-type quantum spin liquid. The predicted
presence of Kitaev interaction in CrI3 and CrGeTe3 systems hints towards the
possibility of realizing quantum spin liquid state in 3$d$ systems. As a
matter of fact, additional calculations we performed (not shown here) indicate
that varying in-plane strain can make the isotropic exchange coupling
vanishing while the Kitaev interaction remains finite, which is promising to
realize quantum spin liquids state in CrI3 and CrGeTe3 systems. Moreover, our
predictions also provide another way to enhance Kitaev interaction in the
“traditional” 4$d$ and 5$d$ systems, that is, for example, to substitute the
light Cl ligand with the heavier I ion in RuCl3 compound. The hybrid source of
SOC to produce strong Kitaev interaction should result in interesting physics
and phenomena in related systems.

We hope that our first-principles calculations and concomitant development of
a simple insightful Hamiltonian (see Eqs. (4) and (5)), along with a tight-
binding model, deepens the understanding of magnetic anisotropy in low-
dimensional systems. The decomposition of the total magnetic anisotropy into
Kitaev and SIA effects further sheds light into the behaviors of other related
systems, such as CrBr3 and CrSiTe3 (see SM[17]), therefore further
demonstrating its relevance and importance.

See supplementary materials[17] for additional details about the methods used
here, as well as further information about our predictions.

All data generated or analysed during this study are included in this
published article (and its supplementary information files).

We thank Zhenglu Li and Jianfeng Wang for useful discussion. This work is
supported by the Office of Basic Energy Sciences under contract ER-46612. H.X.
is also supported by NSFC (11374056), the Special Funds for Major State Basic
Research (2015CB921700), Program for Professor of Special Appointment (Eastern
Scholar), Qing Nian Ba Jian Program, and Fok Ying Tung Education Foundation.
The Arkansas High Performance Computing Center (AHPCC) is also acknowledged.

C.X. and J.F contributed equally to this work. H.X and C.X. conceived the
idea. H.X. and L.B. supervised this work. C.X. performed the DFT calculations
and related analysis, J.F. worked on the tight-binding model. All authors
contribute to the discussion of the results and the writing of the manuscript.

The Authors declare no Competing Financial or Non-Financial Interest.

Correspondence should be addressed to H.X. (email<EMAIL_ADDRESS>and
L.B. (email: laurent@uark.edu).

## References

  * [1] Griffiths, R. B.  Peierls proof of spontaneous magnetization in a two-dimensional Ising ferromagnet.  _Physical Review_ 136, A437 (1964).
  * [2] Mermin, N. D. & Wagner, H.  Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models.  _Physical Review Letters_ 17, 1133 (1966).
  * [3] Soumyanarayanan, A., Reyren, N., Fert, A. & Panagopoulos, C.  Emergent phenomena induced by spin-orbit coupling at surfaces and interfaces.  _Nature_ 539, 509–517 (2016).
  * [4] Miao, N., Xu, B., Zhu, L., Zhou, J. & Sun, Z.  2D intrinsic ferromagnets from Van der Waals antiferromagnets.  _Journal of the American Chemical Society_ 140, 2417–2420 (2018).
  * [5] Huang, B. _et al._ Layer-dependent ferromagnetism in a Van der Waals crystal down to the monolayer limit.  _Nature_ 546, 270–273 (2017).
  * [6] Gong, C. _et al._ Discovery of intrinsic ferromagnetism in two-dimensional Van der Waals crystals.  _Nature_ 546, 265–269 (2017).
  * [7] Xing, W. _et al._ Electric field effect in multilayer Cr2Ge3Te6: a ferromagnetic 2D material.  _2D Materials_ 4, 024009 (2017).
  * [8] Lado, J. L. & Fernandez-Rossier, J.  On the origin of magnetic anisotropy in two dimensional CrI3.  _2D Materials_ 4, 035002 (2017).
  * [9] Li, X. & Yang, J.  CrXTe3 (X= Si, Ge) nanosheets: two dimensional intrinsic ferromagnetic semiconductors.  _Journal of Materials Chemistry C_ 2, 7071–7076 (2014).
  * [10] Carteaux, V., Brunet, D., Ouvrard, G. & Andre, G.  Crystallographic, magnetic and electronic structures of a new layered ferromagnetic compound Cr2Ge2Te6.  _Journal of Physics: Condensed Matter_ 7, 69 (1995).
  * [11] McGuire, M. A., Dixit, H., Cooper, V. R. & Sales, B. C.  Coupling of crystal structure and magnetism in the layered, ferromagnetic insulator CrI3.  _Chemistry of Materials_ 27, 612–620 (2015).
  * [12] Samarth, N.  Condensed-matter physics: Magnetism in flatland.  _Nature_ 546, 216–218 (2017).
  * [13] Kitaev, A.  Anyons in an exactly solved model and beyond.  _Annals of Physics_ 321, 2–111 (2006).
  * [14] Jackeli, G. & Khaliullin, G.  Mott insulators in the strong spin-orbit coupling limit: from Heisenberg to a quantum compass and Kitaev models.  _Physical Review Letters_ 102, 017205 (2009).
  * [15] Liu, H. & Khaliullin, G.  Pseudospin exchange interactions in d 7 cobalt compounds: Possible realization of the kitaev model.  _Physical Review B_ 97, 014407 (2018).
  * [16] Sano, R., Kato, Y. & Motome, Y.  Kitaev-heisenberg hamiltonian for high-spin d 7 mott insulators.  _Physical Review B_ 97, 014408 (2018).
  * [17] See Supplemental Material at for more details, which includes ref.[28, 29, 30, 31, 18, 32, 8, 33, 34, 35, 36, 37, 24, 25, 26, 27, 9, 38, 39, 10, 7, 40] .
  * [18] Xiang, H., Lee, C., Koo, H.-J., Gong, X. & Whangbo, M.-H.  Magnetic properties and energy-mapping analysis.  _Dalton Transactions_ 42, 823–853 (2013).
  * [19] Moriya, T.  Anisotropic superexchange interaction and weak ferromagnetism.  _Physical Review_ 120, 91 (1960).
  * [20] Note that the fact that $J_{\alpha}$ = $J_{\beta}\neq J_{\gamma}$ may appear that the XXZ model is accurate when working in the {$\alpha\beta\gamma$} basis. However, one has to realize that, unlike the XXZ model adopted in Ref.[8], (i) such basis is local in the sense that it depends on the considered Cr pair; and (ii) only the $\alpha$ axis lies in the x-y plane of the film (i.e., the $\beta$ and $\gamma$ axes have out-of-plane components) .
  * [21] Löwdin, P.-O.  On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals.  _The Journal of Chemical Physics_ 18, 365–375 (1950).
  * [22] Xiang, H., Wei, S.-H. & Whangbo, M.-H.  Origin of the structural and magnetic anomalies of the layered compound SrFeO2: a density functional investigation.  _Physical Review Letters_ 100, 167207 (2008).
  * [23] Hu, J. & Wu, R.  Giant magnetic anisotropy of transition-metal dimers on defected graphene.  _Nano Letters_ 14, 1853–1858 (2014).
  * [24] McGuire, M. A. _et al._ Magnetic behavior and spin-lattice coupling in cleavable van der waals layered crcl 3 crystals.  _Physical Review Materials_ 1, 014001 (2017).
  * [25] Carteaux, V., Moussa, F. & Spiesser, M.  2D Ising-like ferromagnetic behaviour for the lamellar Cr2Si2Te6 compound: a neutron scattering investigation.  _Euro. Physics Letters_ 29, 251 (1995).
  * [26] Williams, T. J. _et al._ Magnetic correlations in the quasi-two-dimensional semiconducting ferromagnet CrSiTe3.  _Physical Review B_ 92, 144404 (2015).
  * [27] Liu, B. _et al._ Critical behavior of the quasi-two-dimensional semiconducting ferromagnet CrSiTe3.  _Scientific Reports_ 6 (2016).
  * [28] Kresse, G. & Joubert, D.  From ultrasoft pseudopotentials to the projector augmented-wave method.  _Phys. Rev. B_ 59, 1758 (1999).
  * [29] Blöchl, P. E.  Projector augmented-wave method.  _Phys. Rev. B_ 50, 17953 (1994).
  * [30] Kohn, W. & Sham, L. J.  Self-consistent equations including exchange and correlation effects.  _Physical Review_ 140, A1133 (1965).
  * [31] Momma, K. & Izumi, F.  VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data.  _J. Appl. Crystallogr._ 44, 1272–1276 (2011).
  * [32] Xu, C., Xiang, H. & Bellaiche, L.  Novel multiferroic phases and phenomena in epitaxial (111) BiFeO3 films.  _Advanced Electronic Materials_ 1700332–n/a (2017).
  * [33] Ji, H. _et al._ A ferromagnetic insulating substrate for the epitaxial growth of topological insulators.  _Journal of Applied Physics_ 114, 114907 (2013).
  * [34] Carteaux, V., Ouvrard, G., Grenier, J. & Laligant, Y.  Magnetic structure of the new layered ferromagnetic chromium hexatellurosilicate Cr2Si2Te6.  _Journal of Magnetism and Magnetic Materials_ 94, 127–133 (1991).
  * [35] Feng, J. & Xiang, H.  Anisotropic symmetric exchange as a new mechanism for multiferroicity.  _Physical Review B_ 93, 174416 (2016).
  * [36] Katsura, H., Nagaosa, N. & Balatsky, A. V.  Spin current and magnetoelectric effect in noncollinear magnets.  _Physical Review Letters_ 95, 057205 (2005).
  * [37] Xiang, H., Kan, E., Zhang, Y., Whangbo, M.-H. & Gong, X.  General theory for the ferroelectric polarization induced by spin-spiral order.  _Physical Review Letters_ 107, 157202 (2011).
  * [38] Liu, J., Sun, Q., Kawazoe, Y. & Jena, P.  Exfoliating biocompatible ferromagnetic Cr-trihalide monolayers.  _Physical Chemistry Chemical Physics_ 18, 8777–8784 (2016).
  * [39] Löwdin, P.-O.  A note on the quantum-mechanical perturbation theory.  _The Journal of Chemical Physics_ 19, 1396–1401 (1951).
  * [40] Marzari, N. & Vanderbilt, D.  Maximally localized generalized Wannier functions for composite energy bands.  _Physical review B_ 56, 12847 (1997).