# Words, Subwords, and Morphemes: What Really Matters in the Surprisal-Reading
Time Relationship?

Sathvik Nair  Philip Resnik
Linguistics & Institute for Advanced Computer Studies
University of Maryland
<EMAIL_ADDRESS>
###### Abstract

An important assumption that comes with using LLMs on psycholinguistic data
has gone unverified. LLM-based predictions are based on subword tokenization,
not decomposition of words into morphemes. Does that matter? We carefully test
this by comparing surprisal estimates using orthographic, morphological, and
BPE tokenization against reading time data. Our results replicate previous
findings and provide evidence that _in the aggregate_ , predictions using BPE
tokenization do not suffer relative to morphological and orthographic
segmentation. However, a finer-grained analysis points to potential issues
with relying on BPE-based tokenization, as well as providing promising results
involving morphologically-aware surprisal estimates and suggesting a new
method for evaluating morphological prediction.

##  1 Introduction

There is widespread consensus that human sentence processing includes word-
level prediction (Ehrlich and Rayner, 1981); see Staub (2015) for a review. A
growing body of research is making use of language models as computational
proxies for human prediction at the word level, including traditional $n$-gram
models (Smith and Levy, 2013), syntax-based models Hale (2001); Demberg and
Keller (2008), and more recent work that makes use of neural language models
Goodkind and Bicknell (2018); Wilcox et al. (2020); Shain et al. (2022).

As an overall paradigm, research in this area generally correlates _surprisal_
, computed using corpus-based probability estimates
($-log\Pr(\mbox{word}|\mbox{context})$), against measurable indices of human
processing effort. These include measurement of processing activity using fMRI
Henderson et al. (2016); Bhattasali and Resnik (2021); Shain et al. (2020),
MEG Donhauser and Baillet (2020); Brodbeck et al. (2022), EEG Rabovsky et al.
(2018); Li and Ettinger (2023); Michaelov et al. (2023), and reading times
Demberg and Keller (2008); Smith and Levy (2013); Goodkind and Bicknell
(2018); Wilcox et al. (2020); Shain et al. (2022). This correlation paradigm
has produced useful insights about the role of prediction in language
comprehension. In addition, correlations between language model surprisal and
human processing activity have been taken to be an indication that language
models are capturing relevant aspects of how human language works Ryu and
Lewis (2021); Goldstein et al. (2022).111These two kinds of results are
important to distinguish. In the first case, human beings and their cognitive
systems are the object of study. The second is an instance of what Simon
(1969) referred to as “science of the artificial”, where in this case the
constructed computational system is itself the object of study.

Lost in the shuffle of this research progress, however, is the question of
what, exactly, “prediction at the word level” is supposed to mean. Within
current linguistic theory, the very construct of “word” is receiving critical
attention: albeit controversially, convincing arguments exist that the term
_word_ lacks a consistent, scientifically valid definition that holds
consistently across the full range of human languages Halle and Marantz
(1993); Haspelmath (2017); see Krauska and Lau (2023) for an overview
connecting these theoretical claims to psycholinguistics and neurolinguistics.
Even setting aside that theoretical debate, however, the move to
psycholinguistics and neurolinguistics using large language models has
generated a mismatch between the units analyzed in human studies — typically
orthographic words — and the subword tokens over which LLMs operate via
approaches like WordPiece and Byte-Pair Encoding Wu et al. (2016); Sennrich et
al. (2016).

To take an example, a typical LLM tokenization of the word _decomposition_
using a widely used BPE tokenization Sennrich et al. (2016) yields subword
units dec, om, position.222We use the GPT2 implementation of BPE throughout.
In a typical human subjects experiment, measurements such as first-pass
reading time would typically involve a region from the beginning of the whole
word to its end. More generally, the measurement of activity for a word $w$ in
the human experimentation is related to language model predictions using $w$’s
subword units $s_{1}s_{2}\ldots s_{n}$. How is a correlation computed between
measurements over different units?

The solution adopted by most researchers (Hao et al., 2020; Wilcox et al.,
2020; Stanojević et al., 2022, and others) is to compute surprisal separately
for the $s_{i}$, and then approximate the model’s surprisal for $w$ as the
_sum_ of those individual subword surprisals. The logic behind this choice is
that, if the linking hypothesis behind the work connects surprisal with
cognitive effort, the effort for the entire word should be the sum of the
effort on each of the parts (Smith and Levy, 2013).333In principle one could
be more sophisticated by calculating the uniqueness point within the word in a
given lexicon (Luce, 1986) and only summing up assigned complexity values for
wordpieces that precede that point.

This, however, leads to another question: is that a reasonable approximation,
given that the subwords produced by LLM tokenization bear no clear
correspondence to the subword decompositions in human processing? Consider
again the word _decomposition_. A large body of theoretical and empirical work
would suggest that to the extent subword effort takes place, it would involve
_morphological_ units, in this case de, compos(e), and (i)tion Gwilliams
(2020).444We discuss another illustrative example in Appendix A and Section
2.1 provides details on the morphological segmenter we use in our
experimentation.

The question we set out to answer in this paper, then, is whether the
divergence between LLM subword tokenization and human morphological
decomposition is something to worry about in computational psycholinguistics
and computational neurolinguistics research. Operationalizing the question,
would the sum-of-surprisals approach with morphologically valid units yield a
better correspondence with human measurements than the standard approach using
LLM subword tokens? We would argue that the result is important regardless of
which way the experimentation goes. If morphological units turn out to be a
significantly better fit for human measurements, then cognitive researchers
using LLMs should be using them — which could potentially raise real
challenges given the widespread use of off-the-shelf pretrained LLMs. If
statistically-driven subword units work just as well, then we have checked an
important, previously unchecked box in terms of validating their use.

This work is very much in the spirit of Wilcox et al. (2020), who evaluated a
“fleet” of language models across architectures, plus orthographic $n$-gram
models, against eyetracking and self-paced reading data. However to our
knowledge, this study is the first to consider the assumption that LLM-
subwords can be used in lieu of morphological units.

##  2 Methods

We trained three $n$-gram models under different word segmentation methods and
evaluated them against reading time data from psycholinguistic experiments
conducted in English. Our choice of evaluation corpora and metrics are
consistent with previous literature, such that we were only examining the
effect of word segmentation without model architecture as a confound. Our
implementation, along with instructions for accessing the associated data, is
available at https://github.com/sathvikn/dl-psych-tokenization/.

###  2.1 Language Models

Each $n$-gram model was a 5-gram model trained on the publicly available
section of the Corpus of Contemporary American English (COCA, Davies (2010)).
The models were trained under KenLM Heafield (2011) using modified Kneser-Ney
smoothing Ney et al. (1994). As a control, we used a model trained on COCA
which was trained to predict the next orthographic word without any subword
tokenization. To test BPE-based tokenization, we used the Huggingface
implementation of the GPT2 tokenizer Wolf et al. (2020); Radford et al. (2018)
for each sentence in the corpus, and trained the $n$-gram model on individual
tokens rather than words. Most current psycholinguistic work uses off-the-
shelf GPT2 implementations and GPT2 (and variants) have been shown to be
better fits to reading time data than larger models Shain et al. (2022); Oh
and Schuler (2023). Finally, for a more linguistically informed approach to
word segmentation, we trained an $n$-gram model based on the output of a
morphological transducer Wehrli et al. (2022) that is far more accurate than
WordPiece tokenization at word and sentence-level morpheme segmentation tasks
in a variety of languages, including English Batsuren et al. (2022).555This
analyzer came in second in the SIGMORPHON competition. We chose it because the
first place system’s implementation was not publicly available. The
morphological transducer was based on an encoder-decoder architecture, which
used a two-layer stacked LSTM as the encoder and performed greedy decoding. We
judged the corpus too small for retraining a GPT-style model, and we did not
train an LSTM because Wilcox et al. (2020) conclusively showed 5-gram models
are stronger predictors of results in broad-coverage psycholinguistic
experiments.666In this study we used a publicly available subset of the COCA
corpus for replicability, and it is less than 2/3 the size of the training
dataset for the smallest orthographic GPT2 model in Wilcox et al. (2020). A
near-term aim for future work is to replicate this experimentation with the
full COCA corpus, which requires a license.

###  2.2 Psycholinguistic Evaluation

Once we trained the models, we computed their surprisal estimates for words in
eyetracking and self-paced reading corpora and fit regression models
evaluating surprisal as a predictor of the reading times from these corpora.

####  2.2.1 Corpora

We used averaged eye movement data from the Dundee corpus Kennedy et al.
(2003) and self-paced reading times from the Natural Stories corpus Futrell et
al. (2018) made available by Wilcox et al. (2020). Both corpora are
representative of the material in COCA. In the Dundee corpus, each word’s
fixation time in milliseconds was averaged across 10 English-speaking
participants reading newspaper articles. The Natural Stories corpus consists
of sentences from narrative texts that were edited to include syntactic
constructions that are rare in spoken and written English for psycholinguistic
analysis. The self-paced reading times were recorded from 181 English speakers
who read the texts as they were presented word-by-word. We used the per-word
presentation times that were averaged across participants.

####  2.2.2 Measuring Predictive Power of Surprisal

To compare the reaction times to the models, we computed the surprisal for
each word under each $n$-gram model. The model trained on orthographic words
generated a surprisal for each word, but since the BPE tokenizer and the
morphological transducer used subword information, we tokenized the texts from
the behavioral experiments and computed each token’s surprisal. If a word was
split into multiple tokens, its surprisal under the other two models was the
sum of the tokens’ individual surprisals.777We excluded data for certain words
from our analysis following Goodkind and Bicknell (2018). These were words
preceding and following punctuation, words that contained non-alphabetical
characters, and words that were out of vocabulary for the language models. If
any token under a language model was not in its vocabulary, the entire word
was excluded. We then fit regression models predicting reading time based on
surprisal.

For each word segmentation method, we compared the per-token log likelihood
($\Delta LogLik$) under two multiple linear regression models, following
previous literature to quantify how much surprisal contributes to reading time
prediction, independent of other predictors. One model used the control
features of word length and log unigram frequency to predict reading times as
a baseline model, and the other used these factors in conjunction with
surprisal. If the regression model with surprisal-based features generated
more accurate predictions, the difference between the log likelihoods would be
above zero.

Although $\Delta LogLik$ as the measure of predictive power is standard for
this literature, we highlight two specific methodological details from Wilcox
et al. (2020). First, the predictive power normalizes the regression models’
aggregate log likelihood since we are comparing this metric on corpora with
different sizes. Second, we used 10-fold cross validation to report $\Delta
LogLik$ on a held-out test set. The training and testing data were consistent
across all models for each fold. Reporting the value of a cross-validated
regression model is important to ensure that the predictive power measures
computed on the complete dataset are not the result of overfitting (Yarkoni
and Westfall, 2017). Due to spillover effects from previous words Smith and
Levy (2013), we also included the surprisal, length, and log frequency of the
previous word as predictors in the the regression models for the Dundee
corpus, and similarly for the previous three words for the Natural Stories
corpus.888This difference arises because the corpora were used with different
psycholinguistic tasks Wilcox et al. (2020); we report results for each corpus
separately for this reason, and because the corpora have different sizes and
material.

##  3 Results

All our results show surprisal improves predictions of reading time relative
to the control features, comparably to the 5-gram models’ results in Goodkind
and Bicknell (2018) and Wilcox et al. (2020). Table 1 reports the difference
in per-token log likelihood between the linear regression models predicting
words’ reading times using surprisal as a feature and the models which simply
used length and log frequency as features. We also report a more conventional
measure of effect size using Cohen’s $f^{2}$ in Table 5. The surprisals of the
previous and current word were statistically significant predictors of reading
time for all models on all corpora $(p<0.001)$.999For Natural Stories, the
tokens at $w_{i-3}$ and $w_{i-2}$ also had some predictive power ($p<0.01$ and
$p<0.1$, respectively).

Tokenization Method |  $\Delta LogLik$ for Dundee |  $\Delta LogLik$ for Natural Stories
---|---|---
Orthographic |  0.01 |  0.009
BPE |  0.01 |  0.01
Morphological |  0.01 |  0.008
Table 1: Per-token $\Delta LogLik$ showing surprisal-based regression models
improve over controls.

We also report $\Delta LogLik$ on a held-out test set using 10-fold cross
validation (Figure 1). The training and testing data were consistent across
all models for each fold. Reporting the value of a cross-validated regression
model is important to ensure that the predictive power measures computed on
the complete dataset are not the result of overfitting. We compared the
distribution of these values under a Wilcoxon rank-sum test (Table 6). The
predictive power of surprisal under the models using BPE and the morphological
transducer’s output did not show a statistically significant difference from
the model using orthographic words. For the Natural Stories corpus, the
predictive power of surprisal was lower than the Dundee corpus, which is
expected since the Natural Stories corpus contained rare syntactic
constructions.

Figure 1: Distribution of predictive power of surprisal under models trained
under each tokenization scheme, evaluated using 10-fold cross-validation for
each corpus. There is no major difference in predictive power associated with
tokenization.

On the face of it, these results seem to show that LLM-style tokenization may
not be an issue in psycholinguistic modeling. However, finer-grained analyses
suggest otherwise. First, few words in the psycholinguistic corpora were split
by the BPE tokenizer in the first place. As it turns out, the BPE tokenizer
only split 5% of the tokens in the psycholinguistic corpora (11% when ignoring
stopwords), compared to 25% and 44% respectively for the morphological
analyzer (complete results in Table 3). Moreover, the standard linking theory
for surprisal suggests that effort for the entire word should be the sum of
the subword efforts (Smith and Levy, 2013), and therefore that processing
effort should increase incrementally with the number of units a word is
segmented into. But this appears to be true only for the morphological
tokenization: for BPE tokenization there is a sharp jump from low surprisal
with unsplit words to essentially equal surprisal for words split into 2, 3,
and 4 tokens (Figure 2). The data therefore suggest that surprisal based on
BPE tokenization is less cognitively realistic than surprisal over
morphological units. In addition, replicating the entire analysis separately
for non-segmented and segmented words, we find that the predictive power of
the BPE-based model is significantly worse for words that do get split by the
tokenizer, and this is not true for the morpheme-based model (Figure 3). We
conclude that, despite the aggregate results in Figure 1, LLM-style
tokenization should be viewed with caution in cognitive investigations.

Figure 2: Distributions of surprisal of words with different numbers of
subword tokens, split by corpus and segmentation method.

Figure 3: Replicating the cross-validated analysis of predictive power
separately for whole and split words under the different word segmentation
methods.

##  4 Conclusions

This study was a small, focused contribution that tackled an essential
question previously unaddressed in psycholinguistics research that uses LLMs
and their subword tokenizations. Previous work has demonstrated a linear
relationship between LLM surprisal and human processing as measured by reading
time, but there is good evidence that aspects of human processing are mediated
by _morphological_ units. For cognitive research, mightn’t it be important to
model surprisal using morphological units, rather than distributionally
derived subword tokens that often deviate drastically from a word’s
morphological decomposition?

On the one hand, our work offers the first comparison demonstrating that using
distributionally derived rather than morphological subwords does not affect
_aggregate_ correlations between surprisal and reading times, a widely used
behavioral measurement of human processing. On the other hand, given that
psycholinguistics work is increasingly using proprietary pre-trained models
with non-morphological subword tokenization, far beyond the scale available
for academic model training, our finer-grained analysis indicates that a
degree of caution is warranted for more fine-grained studies. Perhaps our most
important take-away is that cognitive investigations require a careful look at
the cognitive plausibility of the models’ units of analysis.

Our results also suggest new directions for more cognitively realistic models
of prediction in language comprehension. We view the results from our finer-
grained analysis as a step in this direction, and they also suggest going more
deeply into the role of morphological segmentation on an item-by-item basis;
for example, training an LLM on morphological units and evaluating it on
diagnostics for morphosyntactic generalization. Our work here also introduces
_morphological surprisal_ computed automatically using a morphological
segmenter, and validates its predictive power. This would be a natural fit for
further work on morphological prediction at the neural level (Gwilliams,
2020), including looking at the role of morphemes in phoneme prediction
Ettinger et al. (2014), and as a new representational level within integrative
processing models that take phonological, word-level, and sentence-level
contexts into account (Brodbeck et al., 2022). This implementation could be
refined through inferring hierarchical structure over morphological units in
the style of Oseki and Marantz (2020) to conduct larger-scale analyses.

Finally, regarding broader theoretical discussion, we note that surprisal (as
operationalized using an LLM or any other probability estimates) generally
contributes to explanations at Marr’s “computational” level Marr
(1982).101010As an interesting and promising exception, Futrell et al. (2020)
take a step in the direction of the algorithmic/representational level by
bringing memory considerations into the model. Moreover, LLM predictions
represent a black-box combination of categories of information that both
theoretical and experimental considerations suggest are processed in distinct
ways Oseki and Marantz (2020); Oh et al. (2021); Krauska and Lau (2023). We
would therefore argue that, despite their undeniable convenience and power,
the widespread use of LLMs as probabilistic predictors deserves drastically
more careful consideration than it has received if the field is to move in the
direction of deeper insights into human representations and mechanisms in
language processing.

##  5 Limitations

The study was limited to English, and it is possible different results might
obtain in languages with other morphological structures. However, the
morphological transducer can be trained on any language for which a suitable
morphologically segmented corpus is available, and it has already been
evaluated on a multilingual test suite Batsuren et al. (2022), so this is a
promising topic for future work. This is an active area of research,
especially because not all languages have the same notion of what counts as a
“word.” Existing work de Varda and Marelli (2022); Wilcox et al. (2023) has
evaluated predictions on a multilingual eyetracking corpus with some
typological diversity, but still trains transformer language models on subword
tokenization. More work is also needed to see if results vary across mono- and
multi-morphemic words; see Appendix B for an indication that LLM subword
tokenizations can still be problematic at the level of individual predictions,
even for words that do not include much morphological complexity. We also note
that the publicly available version of COCA we used was preprocessed by Yang
et al. (2022), and this may have led to some small discrepancies with the
results from previous studies trained on other academically feasible datasets.

Perhaps our most significant limitation was in using $n$-gram versus state-of-
the-art LLM architectures for our comparison, which in principle may not
generalize to the best models. We would strongly encourage those who are able
to train LLMs at scale to consider offering models with morphologically valid
segmentation, both to facilitate more extensive language model comparisons,
and to support scientific studies involving morphological representations as
articulated in Section 4.

##  6 Ethical Considerations

All data we used are publicly available. Human experimentation was approved by
the institutions who conducted the research, including our own. The software
we used was publicly available, and the trained morphological segmenter was
distributed by the authors of the paper, so our implementation and data
analyses do not require specialized computing hardware.

##  7 Acknowledgements

This material is based upon work supported by ONR MURI Award N00014–18–1–2670.
We would like to thank Silvan Wehrli for providing the English morphological
segmentation model. John Hale, Tal Linzen, Brian Dillon, and Allyson Ettinger
provided invaluable discussion as we formulated our research question, and
Marine Carpuat, Utku Turk, and other members of the Computational Linguistics
& Information Processing and Psycholinguistics groups at UMD provided
insightful feedback on this work during various stages. Last, we thank the
reviewers for suggesting improvements and clarifications to the paper,
particulary comments and questions that motivated our finer-grained analysis.

## References

  * Batsuren et al. (2022) Khuyagbaatar Batsuren, Gábor Bella, Aryaman Arora, Viktor Martinovic, Kyle Gorman, Zdeněk Žabokrtský, Amarsanaa Ganbold, Šárka Dohnalová, Magda Ševčíková, Kateřina Pelegrinová, Fausto Giunchiglia, Ryan Cotterell, and Ekaterina Vylomova. 2022\.  The SIGMORPHON 2022 shared task on morpheme segmentation.  In _Proceedings of the 19th SIGMORPHON Workshop on Computational Research in Phonetics, Phonology, and Morphology_ , pages 103–116, Seattle, Washington. Association for Computational Linguistics.
  * Bhattasali and Resnik (2021) Shohini Bhattasali and Philip Resnik. 2021.  Using surprisal and fMRI to map the neural bases of broad and local contextual prediction during natural language comprehension.  In _Findings of the Association for Computational Linguistics: ACL-IJCNLP 2021_ , pages 3786–3798, Online. Association for Computational Linguistics.
  * Brodbeck et al. (2022) Christian Brodbeck, Shohini Bhattasali, Aura AL Cruz Heredia, Philip Resnik, Jonathan Z Simon, and Ellen Lau. 2022.  Parallel processing in speech perception with local and global representations of linguistic context.  _ELife_ , 11:e72056.
  * Davies (2010) Mark Davies. 2010.  The corpus of contemporary american english as the first reliable monitor corpus of english.  _Literary and linguistic computing_ , 25(4):447–464.
  * de Varda and Marelli (2022) Andrea de Varda and Marco Marelli. 2022.  The effects of surprisal across languages: Results from native and non-native reading.  In _Findings of the Association for Computational Linguistics: AACL-IJCNLP 2022_ , pages 138–144, Online only. Association for Computational Linguistics.
  * Demberg and Keller (2008) Vera Demberg and Frank Keller. 2008.  Data from eye-tracking corpora as evidence for theories of syntactic processing complexity.  _Cognition_ , 109(2):193–210.
  * Donhauser and Baillet (2020) Peter W Donhauser and Sylvain Baillet. 2020.  Two distinct neural timescales for predictive speech processing.  _Neuron_ , 105(2):385–393.
  * Ehrlich and Rayner (1981) Susan F. Ehrlich and Keith Rayner. 1981.  Contextual effects on word perception and eye movements during reading.  _Journal of Verbal Learning and Verbal Behavior_ , 20(6):641–655.
  * Ettinger et al. (2014) Allyson Ettinger, Tal Linzen, and Alec Marantz. 2014.  The role of morphology in phoneme prediction: Evidence from meg.  _Brain and language_ , 129:14–23.
  * Futrell et al. (2020) Richard Futrell, Edward Gibson, and Roger P Levy. 2020.  Lossy-context surprisal: An information-theoretic model of memory effects in sentence processing.  _Cognitive science_ , 44(3):e12814.
  * Futrell et al. (2018) Richard Futrell, Edward Gibson, Harry J. Tily, Idan Blank, Anastasia Vishnevetsky, Steven Piantadosi, and Evelina Fedorenko. 2018.  The natural stories corpus.  In _Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC 2018)_ , Miyazaki, Japan. European Language Resources Association (ELRA).
  * Goldstein et al. (2022) Ariel Goldstein, Eric Ham, Samuel A Nastase, Zaid Zada, Avigail Grinstein-Dabus, Bobbi Aubrey, Mariano Schain, Harshvardhan Gazula, Amir Feder, Werner Doyle, et al. 2022.  Correspondence between the layered structure of deep language models and temporal structure of natural language processing in the human brain.  _bioRxiv_ , pages 2022–07.
  * Goodkind and Bicknell (2018) Adam Goodkind and Klinton Bicknell. 2018.  Predictive power of word surprisal for reading times is a linear function of language model quality.  In _Proceedings of the 8th Workshop on Cognitive Modeling and Computational Linguistics (CMCL 2018)_ , pages 10–18, Salt Lake City, Utah. Association for Computational Linguistics.
  * Gwilliams (2020) Laura Gwilliams. 2020.  How the brain composes morphemes into meaning.  _Philosophical Transactions of the Royal Society B_ , 375(1791):20190311.
  * Hale (2001) John Hale. 2001.  A probabilistic Earley parser as a psycholinguistic model.  In _Second Meeting of the North American Chapter of the Association for Computational Linguistics_.
  * Halle and Marantz (1993) M. Halle and A. Marantz. 1993.  Distributed morphology and the pieces of inflection.  In K. Hale and S.J. Keyser, editors, _The view from Building 20: Essays in honour of Sylvain Bromberger_ , volume 20 of _Current studies in linguistics series_ , pages 111–176. Cambridge, MA: MIT Press, Cambridge, Massachusetts.
  * Hao et al. (2020) Yiding Hao, Simon Mendelsohn, Rachel Sterneck, Randi Martinez, and Robert Frank. 2020.  Probabilistic predictions of people perusing: Evaluating metrics of language model performance for psycholinguistic modeling.  In _Proceedings of the Workshop on Cognitive Modeling and Computational Linguistics_ , pages 75–86, Online. Association for Computational Linguistics.
  * Haspelmath (2017) Martin Haspelmath. 2017.  The indeterminacy of word segmentation and the nature of morphology and syntax.  _Folia linguistica_ , 51(s1000):31–80.
  * Heafield (2011) Kenneth Heafield. 2011.  KenLM: Faster and smaller language model queries.  In _Proceedings of the Sixth Workshop on Statistical Machine Translation_ , pages 187–197, Edinburgh, Scotland. Association for Computational Linguistics.
  * Henderson et al. (2016) John M Henderson, Wonil Choi, Matthew W Lowder, and Fernanda Ferreira. 2016.  Language structure in the brain: A fixation-related fmri study of syntactic surprisal in reading.  _Neuroimage_ , 132:293–300.
  * Kennedy et al. (2003) Alan Kennedy, Robin Hill, and Joël Pynte. 2003.  The dundee corpus.  In _Proceedings of the 12th European conference on eye movement_.
  * Krauska and Lau (2023) Alexandra Krauska and Ellen Lau. 2023.  Moving away from lexicalism in psycho- and neuro-linguistics.  _Frontiers in Language Sciences_ , 2.
  * Li and Ettinger (2023) Jiaxuan Li and Allyson Ettinger. 2023.  Heuristic interpretation as rational inference: A computational model of the n400 and p600 in language processing.  _Cognition_ , 233:105359.
  * Luce (1986) Paul A Luce. 1986.  A computational analysis of uniqueness points in auditory word recognition.  _Perception & Psychophysics_, 39(3):155–158.
  * Marr (1982) David Marr. 1982.  _Vision: A Computational Investigation into the Human Representation and Processing of Visual Information_.  Henry Holt and Co., Inc., USA.
  * Michaelov et al. (2023) James A Michaelov, Megan D Bardolph, Cyma K Van Petten, Benjamin K Bergen, and Seana Coulson. 2023.  Strong prediction: Language model surprisal explains multiple n400 effects.  _Neurobiology of Language_ , pages 1–71.
  * Ney et al. (1994) Hermann Ney, Ute Essen, and Reinhard Kneser. 1994.  On structuring probabilistic dependences in stochastic language modelling.  _Computer Speech & Language_, 8(1):1–38.
  * Oh et al. (2021) Byung-Doh Oh, Christian Clark, and William Schuler. 2021.  Surprisal estimators for human reading times need character models.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 3746–3757, Online. Association for Computational Linguistics.
  * Oh and Schuler (2023) Byung-Doh Oh and William Schuler. 2023.  Why does surprisal from larger transformer-based language models provide a poorer fit to human reading times?  _Transactions of the Association for Computational Linguistics_ , 11:336–350.
  * Oseki and Marantz (2020) Yohei Oseki and Alec Marantz. 2020.  Modeling morphological processing in human magnetoencephalography.  _Proceedings of the Society for Computation in Linguistics_ , 3(1):209–219.
  * Rabovsky et al. (2018) Milena Rabovsky, Steven S Hansen, and James L McClelland. 2018.  Modelling the n400 brain potential as change in a probabilistic representation of meaning.  _Nature Human Behaviour_ , 2(9):693–705.
  * Radford et al. (2018) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2018.  Language models are unsupervised multitask learners.
  * Ryu and Lewis (2021) Soo Hyun Ryu and Richard Lewis. 2021.  Accounting for agreement phenomena in sentence comprehension with transformer language models: Effects of similarity-based interference on surprisal and attention.  In _Proceedings of the Workshop on Cognitive Modeling and Computational Linguistics_ , pages 61–71, Online. Association for Computational Linguistics.
  * Sennrich et al. (2016) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016.  Neural machine translation of rare words with subword units.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1715–1725, Berlin, Germany. Association for Computational Linguistics.
  * Shain et al. (2020) Cory Shain, Idan Asher Blank, Marten van Schijndel, William Schuler, and Evelina Fedorenko. 2020.  fmri reveals language-specific predictive coding during naturalistic sentence comprehension.  _Neuropsychologia_ , 138:107307.
  * Shain et al. (2022) Cory Shain, Clara Meister, Tiago Pimentel, Ryan Cotterell, and Roger Philip Levy. 2022.  Large-scale evidence for logarithmic effects of word predictability on reading time.
  * Simon (1969) Herbert A. Simon. 1969.  _The sciences of the artificial._ The sciences of the artificial. The MIT Press, Cambridge, MA, US.  Pages: xii, 123.
  * Smith and Levy (2013) Nathaniel J Smith and Roger Levy. 2013.  The effect of word predictability on reading time is logarithmic.  _Cognition_ , 128(3):302–319.
  * Stanojević et al. (2022) Miloš Stanojević, Jonathan R Brennan, Donald Dunagan, Mark Steedman, and John T Hale. 2022.  Modeling structure-building in the brain with ccg parsing and large language models.  _arXiv preprint arXiv:2210.16147_.
  * Staub (2015) Adrian Staub. 2015.  The Effect of Lexical Predictability on Eye Movements in Reading: Critical Review and Theoretical Interpretation: Predictability and Eye Movements.  _Language and Linguistics Compass_ , 9(8):311–327.
  * Wehrli et al. (2022) Silvan Wehrli, Simon Clematide, and Peter Makarov. 2022.  CLUZH at SIGMORPHON 2022 shared tasks on morpheme segmentation and inflection generation.  In _Proceedings of the 19th SIGMORPHON Workshop on Computational Research in Phonetics, Phonology, and Morphology_ , pages 212–219, Seattle, Washington. Association for Computational Linguistics.
  * Wilcox et al. (2020) Ethan Gotlieb Wilcox, Jon Gauthier, Jennifer Hu, Peng Qian, and Roger Levy. 2020\.  On the predictive power of neural language models for human real-time comprehension behavior.  _arXiv preprint arXiv:2006.01912_.
  * Wilcox et al. (2023) Ethan Gotlieb Wilcox, Tiago Pimentel, Clara Meister, Ryan Cotterell, and Roger P Levy. 2023.  Testing the predictions of surprisal theory in 11 languages.  _arXiv preprint arXiv:2307.03667_.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander Rush. 2020.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online. Association for Computational Linguistics.
  * Wu et al. (2016) Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V Le, Mohammad Norouzi, Wolfgang Macherey, Maxim Krikun, Yuan Cao, Qin Gao, Klaus Macherey, et al. 2016.  Google’s neural machine translation system: Bridging the gap between human and machine translation.  _arXiv preprint arXiv:1609.08144_.
  * Yang et al. (2022) Jinbiao Yang, Antal Van den Bosch, and Stefan L Frank. 2022.  Unsupervised text segmentation predicts eye fixations during reading.  _Frontiers in Artificial Intelligence_ , 5:11.
  * Yarkoni and Westfall (2017) Tal Yarkoni and Jacob Westfall. 2017.  Choosing prediction over explanation in psychology: Lessons from machine learning.  _Perspectives on Psychological Science_ , 12(6):1100–1122.

##  Appendix A Comparing Word Segmentation Methods

In this illustrative example, the BPE-based tokenizer fails to split up most
multimorphemic words. When it succeeds, the words are not segmented by
morpheme (relegated, fringes). The morphological transducer is able to make
cognitively plausible choices involving tense (relegates), possessives (its),
and plurals (fringes). It also splits up more complex words into roots,
prefixes, and suffixes (coverage, reporters, journalistic, community). Both
tokenizers marked word boundaries in their output, although they are not shown
in this example.

Segmentation Method |  Sentence
---|---
Orthographic Words |  the sporadic nature of press coverage of the court often relegates its reporters to the fringes of the journalistic community
BPE Tokenizer Output |  the sporadic nature of press coverage of the court often releg ates its reporters to the fr inges of the journalistic community
Morphological Transducer Output |  the sporadic nature of press cover age of the court often relegate s it s re port er s to the fringe s of the journal istic commune ity
Table 2: Illustrative example of the same sentence from COCA tokenized
orthographically, morphologically, and using BPE.

In Tables 3 and 4 we report the number of words that were split by both the
BPE tokenizer and the morphological segmenter in the psycholinguistic corpora.

Corpus |  Tokens Per Word |  Percent of Corpus |  Percent of Corpus Excluding Stopwords
---|---|---|---
Dundee |  1 |  94.4 |  88.5
2 |  4.19 |  8.68
3 |  1.22 |  2.53
4 |  0.104 |  0.217
5 |  0.005 |  0.011
Natural Stories |  1 |  95.2 |  89.9
2 |  3.85 |  8.01
3 |  0.971 |  2.02
4 |  0.016 |  0.03
Table 3: Words in the psycholinguistic corpora split into different numbers of tokens by the BPE tokenizer. Corpus |  Tokens Per Word |  Percent of Corpus |  Percent of Corpus Excluding Stopwords
---|---|---|---
Dundee |  1 |  75.7 |  55
2 |  21 |  38.3
3 |  3 |  6.18
4 |  0.218 |  0.451
5 |  0.011 |  0.022
Natural Stories |  1 |  76.9 |  58.3
2 |  20.9 |  37.1
3 |  2.05 |  4.27
4 |  0.125 |  0.26
Table 4: Words in the psycholinguistic corpora split into different numbers of
tokens by the morphological segmenter.

##  Appendix B Examples of Surprisal Differences for Morphemes and BPE tokens

Figures 4 and 5 provide two illustrations of surprisal differences between
subword segmentations. Note the major difference in the surprisal of bulb when
summed over BPE tokens when compared to morphological units bulb s. In the
sentence in Figure 4, the GPT2 tokenizer split _tulips_ into tul ips and did
not split _bulbs_. It is reasonable for a human comprehender to expect the
word bulb immediately after tulip since they would cooccur frequently in text,
but it is less predictable after lip. This is reflected in the higher
surprisal under the BPE-based model.

Figure 4: Differences in surprisal of a plural _(bulbs)_ under $n$-gram models
using morphological and BPE-based tokenization, from a sentence in the Natural
Stories corpus.

To take a more morphologically complex example, this difference also occurs
for _carefully_ in Figure 5, which is not split up by the BPE tokenizer. The
morphological transducer, on the other hand, split _carefully_ into care ful
ly, producing a much lower surprisal. This suggests that further item-wise
comparisons involving words with more morphologically relevant units may be
worth investigating.

Figure 5: Differences in surprisal of an adverb _(carefully)_ under $n$-gram
models using morphological and BPE-based tokenization, from a sentence in the
Dundee corpus.

##  Appendix C Statistical Testing for Predictive Power of Surprisal

###  C.1 Effect Sizes

In Table 5 we report effect sizes for the regression models trained on the
full psycholinguistic corpora via the widely used Cohen’s $f^{2}$. However,
since our work was replicating a subliterature that almost exclusively uses
$\Delta LogLik$, we used that measurement in the main body of the paper. We
refer readers to Goodkind and Bicknell (2018) for further statistical
justification of that choice.

Tokenization Method |  $f^{2}$ for Dundee |  $f^{2}$ for Natural Stories
---|---|---
Orthographic | 0.021 |  0.018
BPE | 0.021 |  0.02
Morphological | 0.021 |  0.017
Table 5: Effect size comparing surprisal as a feature of regression models
relative to controls.

###  C.2 Statistical Significance Testing

For the aggregate analysis (Figure 1, we used a Wilcoxon rank-sum test to
compute significance. We find no statistically significant difference between
the $\Delta LogLik$ estimated for the folds for the two subword tokenization
methods relative to predictions over orthographic words.

Tokenization | Corpus | $W$ | $p$
---|---|---|---
BPE | Dundee | 43 | 0.63
Morphological | Dundee | 50 | 1
BPE | Natural Stories | 41 | 0.53
Morphological | Natural Stories | 53 | 0.85
Table 6: Wilcoxon rank-sum test comparing distributions of $\Delta LogLik$
from cross-validation results of morphological and BPE versus orthographic
surprisal.


# A Study of Highly Frustrated Spin Systems with mixed PEPS in Infinite
Honeycomb Lattice

Huan He  Key Laboratory of Quantum Information, CAS, University of Science and
Technology of China  Zhen Wang  Key Laboratory of Quantum Information, CAS,
University of Science and Technology of China  Chuanfeng Li  Key Laboratory of
Quantum Information, CAS, University of Science and Technology of China
YongJian Han 111Email address<EMAIL_ADDRESS>Key Laboratory of Quantum
Information, CAS, University of Science and Technology of China  Guangcan Guo
Key Laboratory of Quantum Information, CAS, University of Science and
Technology of China

###### Abstract

Highly frustrated spin systems represent a central and challenging problem in
condensed mater physics. To this problem, we introduce an algorithm based on
mixed projected entangled pair states (_m_ -PEPS), which is a novel type of
tensor network. We use the famous Kitaev model on an infinite honeycomb
lattice, which can be solved exactly, as a benchmark. With very limited
parameters and finite scaling, our calculation results are in good agreement
with the exact results, indicating the efficiency of our algorithm. After
presenting the benchmark, we investigate the Kitaev-Heisenberg model, which
was proposed to describe the effective magnetic momentum interaction in
iridate Na2IrO3 which may be used to realize the spin liquid phase. However,
our calculations suggest that the gapless spin liquid phase is not robust at
the thermodynamic limit, and thus this phenomenon is very difficult to
observe.

###### pacs:

02.70.-c, 75.10.Kt, 03.67.Mn, 05.50.+q, 75.10.Jm

Frustrated spin systems play an essential role in two-dimensional many-body
physics and have attracted the attention of condensed matter physicists for
decadesbalents (2). In particular, some exotic ground states, such as the
quantum spin liquid state (which has recently been strongly supported
experimentally in the Kagome lattice by Y. Leeexperi (3) _et al_), may play a
key role in explaining high-$T_{c}$ superconductivity, as proposed by P. W.
Andersonand (1). The Kitaev model on a honeycomb latticekitaev (4), which can
be solved exactly and supports a spin liquid phase, is an important frustrated
model in two-dimensional physics. In addition, this proposed model was
proposed that it may be realized in a cold atom systemluming (5), or in the
iridate Na2IrO3 KHmodel (19). The exact solution shows that the spin liquid
phase is gapless, and a gap opens when a magnetic field, which breaks the
time-reversal symmetry, is applied to this model as a perturbation.
Unfortunately, general two-dimensional frustrated spin models - even slightly
modified Kitaev models such as the Kitaev model perturbed with a magnetic
field - cannot be solved exactly.

Due to the lack of analytical methods for general many-body systems, numerical
analysis remains the main method for understanding such physics. However,
traditional algorithms, such as the quantum Monte Carlo (QMC) method and
density matrix renormalization group (DMRG) method, have failed to simulate
frustrated systems in two dimensions: QMC suffers from the notorious ‘sign
problem’ sign (15) while DMRG is limited to one dimension and generally lose
its power in higher dimensional systemsdmrg1 (6, 7). Fortunately, the recently
developed tensor network (TN) algorithms, such as, the algorithm based on
projected entangled pair states (PEPS)PEPS1 (16, 17, 18) which is a natural
generalization of the DMRG algorithm to higher dimensions, has shown great
potential for addressing this problem. The algorithm based on infinite PEPS
has been used to obtain a reasonable phase diagram of the frustrated
antiferromagnetic $J_{1}-J_{2}$ Heisenberg model on a checkerboard latticechan
(14). More recently, a spin liquid phase was claimed near the maximally
frustrated region ($J_{2}\sim 0.5J_{1}$) for the spin $1/2$ $J_{1}-J_{2}$
antiferromagnetic Heisenberg model on a square lattice for the TN
algorithmwangling (12, 13). These results show that the TN algorithms are
powerful tools for exploring the frustrated systems and for evaluating some
new physics beyond the traditional methods.

Despite the success of the TN algorithm, the study of highly frustrated large
systems, such as those represented by the Kitaev model and related models,
remains a challenge. The general TN algorithm, such as the algorithm based on
PEPS, cannot give the appropriate results for the Kitaev model, because the
limited ability of the calculation constrains the bond dimension $D$ to a
relatively small number (typically $<10$). Thus new methods must be developed
to highly frustrated systems with the current calculation abilities. In this
letter, we introduce an algorithm based on infinite mPEPS, which is a mixture
of the bosonic PEPS (bPEPS)PEPS2 (17) and the fermionic PEPS (fPEPS)fPEPS1 (8,
9). We will demostrate that this novel method can be used to efficiently study
highly frustrated systems with a relatively small parameter $D$. We use the
Kitaev model as a benchmark and obtain satisfactory results, for both the
energy, and the correlations of this model. After the benchmark calculation,
we apply this method to the Kitaev-Heisenberg model, which is proposed to
describe the effective magnetic momentum interaction model of the iridate
Na2IrO3 KHmodel (19, 10), to obtain a phase diagram. Our calculations indicate
that the gapless spin liquid region is very narrow in the thermodynamic limit
and that would thus be very difficult to observe in the iridate Na2IrO3, even
with the Kitaev -Heisenberg model.

Mixed Projected Entangled-Pair States The TN algorithm is a variational method
based on TN states with special structures. A general TN state is defined in
two steps: first, we describe a configuration of entangled states between
virtual particles (we call each entangled state a bond, connected by a line,
as shown in Fig.1), which determines the structure and the upper bound of the
entanglement of a TN state; second, we define projectors to project the
virtual-particle space into real physical space. The different TN algorithms
are distinguished by the structure of the TN state, such as PEPS, or the
string bond state. The variational parameter space of TN algorithms is
determined by the parameters of the projectors, the number of which is
dependent on a polynomial of the dimension of the bond $D$, the dimension of
the physical space $d$ and the number of projectors (as determined by the
structure of the lattice and the symmetry of the state). For the special TN
state, PEPS, the configuration of the bonds is completely determined by the
lattice where each edge of the lattice corresponds to a bond, and each site
has an on-site projector. Without a loss of generality, we use the honeycomb
lattice as an example. In general, the entangled pair can be written as:
$\mbox{EPR}_{e}=\sum_{i}{|i\rangle_{e_{1}}|i\rangle_{e_{2}}}$, where
$|i\rangle$ is the state of virtual particle with dimension $D$ and $e$
denotes an edge connecting sites $e_{1}$ and $e_{2}$ in the lattice.
Algebraically, the state $|i\rangle$ can be represented by local creation
operators from the vacuum through the second quantization method. In the
original PEPS case, the entangled pair state can be represented as bosonic
operator $a^{{\dagger}}$ as:
$\mbox{EPR}_{e}=\sum_{i}{a^{{\dagger}i}_{e_{1}}a^{{\dagger}i}_{e_{2}}|Vac\rangle}$.

Similarly, the projector of each site on a honeycomb lattice can also be
represented as bosonic operators as:
$P^{[l_{1}l_{2}]}=\sum_{p,i,j,k}T^{[l_{1}l_{2}]}_{pijk}b^{{\dagger}p}a_{x}^{i}a_{y}^{j}a_{z}^{k}|vac\rangle\langle
Vac|$, where $[l_{1}l_{2}]$ is the coordinate of the site, $b^{\dagger}$ is
the bosonic creation operator of the physical particle, and
$a_{x}$,$a_{y}$,$a_{z}$ are the bosonic annihilation operators of the virtual
particles connected to edges $x$,$y$ and $z$, respectively (see Fig.1). In
this case $|vac\rangle$ is the vacuum of the physical space, and $|Vac\rangle$
is the vacuum of the virtual space. Using the representation of the entangled
pairs and projectors, the PEPS state can be represented as:
$|\psi\rangle=\prod_{l_{1},l_{2}}P^{[l_{1}l_{2}]}\prod_{e}\mbox{EPR}_{e}$,
where all of the virtual particles are projected as physical particles.

In general, the virtual particles in an EPR pair are not necessarily bosons.
In solving the many-body physics of fermionic systems, it is natural to extend
the bPEPS formulism to the fPEPS formulism. This can be completed by replacing
the bosonic operators in the entangled pairs and projectors, as shown in
Eq.(1) with fermionic operators. In the fPEPS formulism, the parity of the
projectors at different sites of the lattice, which determine the exchange
character of the projectors, should be the same even. This requirement
guarantees that the fPEPS is well defined, that is, the order of the
projectors is not important for the state beyond the global phase.

Here we focus on a frustrated spin system, which exhibits bosonic statistics
in physical space. Thus, we propose an algorithm based on the mixed PEPS
(mPEPS), which introduces fermionic statistics to the virtual-particle space
(same as fPEPS) and bosonic statistics in the physical-particle space
(different from fPEPS). Therefore, the projectors of the mPEPS should be
expressed in the second quantized form as:
$\sum_{p,i,j,k}T_{pijk}b^{{\dagger}p}a_{x}^{i}a_{y}^{j}a_{z}^{k}|vac\rangle\langle
Vac|$, where $b^{\dagger}$ satisfies: $[b_{i},b^{\dagger}_{j}]=\delta_{ij}$
and $a_{x}$,$a_{y}$,$a_{z}$ satisfy:
$\\{a_{i},a^{\dagger}_{j}\\}=\delta_{ij}$. To validate this definition, the
parity of the fermionic operators in the projector should also be same for
every site (for example, it could be even). Although there is some overlap,
the state spaces of the PEPS, fPEPS and the mPEPS are expected to be
different, even for the same parameters.

After the defining mPEPS, the ground state of a system in variational space
can be found by the imaginary evolution method which is performed in the same
manner as in the PEPS method; in particularly, we use the simple update method
to approximate the environment in the same manner as PEPS in two
dimensionsxiangtao (24). When the ground state is stable under the imaginary
evolution, the physical quantities of the system can be calculated through a
complex contraction process (for more details, see the Supplementary
Information). In the contraction, the bosonic and fermionic operators must be
carefully handled(see in the Supplementary Information).

Figure 1: (Color on line) Representation of the mPEPS on a honeycomb lattice:
each line denotes an entangled pair state; the blue balls denote virtual
fermions; and the red balls denote bosons; the circle around each site
represents an on-site projector $P$. Different link directions are denoted by
$x$,$y$ and $z$, which are used in the Kitaev model. In the following
calculations, we use an infinite lattice with translational symmetry and six
different projectors, $P_{1},P_{2},\cdots,P_{6}$ Figure 2: (Color on line)
Nearest Neighbor(NN) correlation $\langle\sigma^{z}\sigma^{z}\rangle_{z}$ of
the Kitaev model along the line $J_{x}=J_{y}$, which is indicated in the left
upper inset by a dashed line. The red dots represent the results obtained by
the mPEPS method with finite scaling of $D$; the black line shows the exact
results of $\langle\sigma^{z}\sigma^{z}\rangle$ which is given as the formula
$\langle\sigma^{z}_{i}\sigma^{z}_{j}\rangle=\frac{\sqrt{3}}{16\pi^{2}}\int_{BZ}cos\theta(k_{1},k_{2})dk_{1}dk_{2}$baskaran
(11); the blue triangles display the results obtained by the PEPS method using
the same parameters of mPEPS. The right lower inset shows the
$\langle\frac{\partial^{2}H}{\partial J^{2}}\rangle$, which displays the
transition point near 0.46.

Kitaev Model as the Benchmark Here, we use the Kitaev model, which is a highly
frustrated spin model in two dimensions, as a benchmark of the mPEPS
algorithm. The Kitaev model is defined on an infinite honeycomb latticekitaev
(4), and its interaction along different directions (see Fig.1) vary strongly
as: $H=\sum_{\gamma-link}J_{\gamma}\sigma^{\gamma}_{i}\sigma^{\gamma}_{j}$.
This model can be exactly solved with two resultant phases: the gapped phase
and the gapless spin liquid. There are two types of quasi-particles in the
gapped phase: fermions and vortices. The vortices associated with a $Z_{2}$
gauge field and a factor $-1$ will appear when the fermion moves around the
vortex. A gap opens when a magnetic field is applied to the gapless spin
liquid phase. The non-abelian anyon, which may play an important role in
topological quantum computation, will also appear in this phase. In our
calculation, we normally set the parameters of the model to satisfy
$J_{x}+J_{y}+J_{z}=1$ and $J_{x}=J_{y}$, which gives a line connecting a
vertex to the middle point of the corresponding edge (see the red dashed line
in the inset of Fig.2). After conducting the stability tests of the algorithm
with respect to the parameters (particularly, the truncation parameter
$D_{cut}$ in the contraction process, see details in the Supplementary
Information), we compare the energy of the system along the chosen line with
the exact solution. Result show good agreement (see the Supplementary
Information). In addition, the nearest neighbor(NN) correlation
$\langle\sigma^{z}\sigma^{z}\rangle$ which corresponds to
$\langle\frac{\partial H}{\partial Jz}\rangle$ along the line is also shown in
the Fig.2. The black line shows the exact value of this correlation as derived
inbaskaran (11), which is in good agreement with our calculation (red dots).
The inset shows the $\langle\frac{\partial^{2}H}{\partial Jz^{2}}\rangle$,
which exhibits a sharp transition point near 0.46 and indicates a phase
transition in this model, although it slightly deviates from the exact result
of 0.5.

We must emphasize that the virtual dimension of the system $D$ is rather
limited by our calculation ability. We can only calculate results for the
dimension up to 8. And we conduct the finite scaling for this limited
dimension in this model, and the points shown in Fig.2 were obtained with
those finite scaling. To compare our findings with the traditional PEPS case,
we also present the results from the PEPS with the same parameters and finite
scaling employed in Fig.2 (blue triangles). It is clear that the present
method improves the results obtain with limited calculation ability. More
importantly, the PEPS method does not give an indication of the phase
transition, while our mPEPS technique can. However, in simpler cases, such as
two-dimensional Ising model with a transverse magnetic field on the honeycomb
lattice, the mPEPS and PEPS methods give nearly the equivalent results.

Analytical resultkitaev (4) suggest that the gapless spin liquid will exhibit
a gap upon the addition of a magnetic field in the direction
$\textbf{e}_{x}+\textbf{e}_{y}+\textbf{e}_{z}$ ($(1,1,1)$ direction) as a
perturbation, which breaks the time-reversal symmetry. Unfortunately, the
Kitaev model with magnetic field cannot be solved exactly. We can only study
this transition by numerical methods. Our mPEPS calculations show this
transition directly by $-\frac{\partial{H}}{\partial{\vec{h}}}$ where $h$ is
the magnitude of the magnetic field (see Fig.3), and these calculations
strongly suggest a first order phase transition occurring at very low magnetic
fields, $h<0.005$.

Application to the Kitaev-Heisenberg Model Having verified the mPEPS method in
the Kitaev model, we extend the application of this method to the Kitaev-
Heisenberg model. The Hamiltonian of this model can be written as:
$H_{HK}=(1-\alpha)H_{H}+2\alpha H_{K}$, where $H_{H}$ is the standard
Heisenberg interaction: $\sum_{\gamma=x,y,z}\sigma^{\gamma}\sigma^{\gamma}$
and the $H_{K}=\sum_{\gamma-link}\sigma^{\gamma}_{i}\sigma^{\gamma}_{j}$ is
Kitaev interaction. With the Khaliullin transformationKHmodel (19), which
performs rotations on different sites, this model can be exactly solved at the
point $\alpha=0.5$ and $\alpha=1$. When $\alpha=0.5$, the ground state of the
transformed system can be exactly solved, and is found to be a ferromagnetic
state. With the inverse Khaliullin transformation, the ferromagnetic state
will be transformed to the stripy anti-ferromagnetic state. When $\alpha=1$,
this model is equivalent to the Kitaev model, with $J_{x}=J_{y}=J_{z}$, whose
ground state can be exactly solved as a gapless spin liquid state. For the
parameters beyond $\alpha=0.5$ and $\alpha=1$, the model can not be exactly
solved and can only be examined by numerical calculations. Former studies have
employed several numerical methods, such as exact diagonalizationKHmodel (19,
21) and DMRGJiang (20), and have found three different phases: the Neel order
phase, the stripy anti-ferromagnetic phase and the spin liquid phase. These
previous calculations, which were based on the systems of relatively small
finite sizeKHmodel (19) or quasi one-dimensionalJiang (20), showed that phase
transition points occur near $\alpha=0.4$ and $\alpha=0.8$ at zero
temperature.

We applied the mPEPS method to this model on an infinite two-dimensional
lattice which naturally satisfies the thermodynamic limit. We directly
calculate the spin configuration of the ground state with $\alpha=0.5$ which
is indeed a ferromagnetic state, as predicted by the exact results. As with
the benchmark calculation of the Kitaev model with $J_{x}=J_{y}=J_{z}=1/3$
(corresponding to $\alpha=1$), here we also give $\frac{dH}{d\alpha}$, which
can be used to determine the phase transition point, as shown in Fig.4.

Figure 3: (Color on line) The transition according to the modulation of
$-\frac{\partial{H}}{\partial{\vec{h}}}$ with a magnetic field $h$. The Kitaev
model is initially in the spin liquid phase with parameters $J_{x}=Jy=Jz=1/3$,
as denoted by the green point. The red dots indicate results for the parameter
$D=4$, and the black squares display the results for $D=6$. Although there is
numerical difference between the two parameters (without finite scaling for
$D$), the transition is sufficiently sharp enough and at the same position.
This finding provides strong evidence that there is a first phase transition
occurs near $h=0.005$ under finite scaling.

In our infinite system with six-site periodicity, the second transition point
is dramatically different from the former result($\alpha=0.8$), while the
first transition point, at $\alpha=0.46$, is slightly different from
$\alpha=0.4$. To confirm our calculation, we also determine the NN correlation
$\langle\vec{\sigma}\vec{\sigma}\rangle$, which is used as an indication of
the phase transition in referenceKHmodel (19). The results also support the
occurrence of a phase transition at the same position. In our calculations,
the second phase transition point occurs at $\alpha\approx 0.98$ for $D=4$,
and at $\alpha\approx 0.99$ for $D=6$. This findings suggest that the
transition point will tend to $\alpha=1$ with the finite scaling for $D$. Our
results indicate that the gapless spin liquid phase will be destroyed by the
Heisenberg perturbation very quickly, and thus this phase is not robust for
this type (Heisenberg) of perturbation. The direct implication of our results
indicates that it would be very difficult to observe a Kitaev-type gapless
spin liquid in iridate Na2IrO3. The findings of a decreasing spin liquid phase
in the infinite system (at the thermodynamics limit), compared to the finite
exact diagonalization and DMRG in quasi one-dimension with cylinder boundary
condition, is not so surprising. The size of the system plays a subtle role in
the spin liquid phase. Z. Meng _et al_ have reportedZYMeng (22) a spin liquid
phase between the semimetal phase and an antiferromagnetic Mott insulator in
the Fermi-Hubbard model on the honeycomb lattice, with the inclusion of 648
sites, by using projective QMC method. However, when S. Sorella _et al_
japanese (23) reconsidered this problem by including up to 2592 sites, the
evidence of spin liquid was lost. In our case, the infinite lattice may play a
similar role.

Figure 4: (Color online) The values of $\frac{dH}{d\alpha}$ (a) and the NN
correlation $\langle\vec{\sigma}\vec{\sigma}\rangle$ (b) versus $\alpha$ for
$D=4$ (black and green squares) and $D=6$ (red and blue dots). The left lower
inset shows the magnification near $\alpha=1$ and the right upper inset
displays the configuration of the spin under the Khaliullin transformation for
$\alpha=0.5$. It is clear that two phase transitions occur. The first
transition point is near $\alpha=0.46$ which is only slightly different with
$0.4$. However, the second phase transition occurs at $\alpha\approx 0.98$
which is qualitatively different from $0.8$. Our results show that the spin
liquid region is very narrow and that the gapless spin liquid will be rapidly
destroyed by the Heisenberg perturbation.

To summarize, we have introduced a novel type of TN (mPEPS) algorithm to
attack the highly frustrated spin systems. Our method is successful in
simulating the Kitaev model, with respect to both the energy and the
correlations of the ground state. Having established this method, we applied
it to the Kitaev-Heisenberg model. Our calculations show that the spin liquid
region is very narrow, and thus this state would be very difficult to observe
experimentally. The success of our method opens the possibility of considering
virtual particles in a tensor network with more complex statistical
characteristics, such as the fractional statistics of anyons. These new
statistics may be convenient for studying the system which excitations have
the similar statistics.

We acknowledge support from the Chinese National Fundamental Research Program
2011CB921200 and the Fundamental Research Funds for Central Universities Nos.
WK2470000006, WK2470000004, WJ2470000007 and NSFC11105135.

## References

  * (1) P.W. Anderson, Science 235, 1196 (1987).
  * (2) Leon Balents, Nature 464,(2010).
  * (3) Tian-Heng Han _et al_ Nature (2012).
  * (4) A. Kitaev, Ann. Phys. (N.Y.) 321, 2 (2006).
  * (5) L.-M. Duan, E. Demler, M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
  * (6) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * (7) U. Schollwock, Rev. Mod. Phys. 77, 259 (2005).
  * (8) C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. A 81, 052338 (2010).
  * (9) P. Corboz, R. Orus, B. Bauer, G. Vidal, Phys. Rev. B 81, 165104 (2010).
  * (10) Actually, the Kitaev-Heisenberg model is not the only model proposed for describing the effective magnetic momentum interaction model of iridate Na2IrO3. Other models have proposed. Such as that described in: Y. Singh _et al_ , Phys. Rev. Lett. 108, 127203 (2012). In this work we only consider the Kitaev-Heisenberg model.
  * (11) G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev. Lett. 98, 247201 (2007)
  * (12) L. Wang, D. Poilblanc, Z.-C. Gu, X.-G. Wen, and F. Verstraete, arXiv:1301.4492.
  * (13) H. C. Jiang, Hong Yao and L. Balents, Phys. Rev. B 86,024424(2012).
  * (14) Y. H. Chan, Y.J. Han, L. M. Duan, Phys. Rev. B 84, 224407(2012).
  * (15) E. Y. Loh, _et al_ Phys. Rev. B 41, 9301(1990)
  * (16) J. Jordan, R. Orus, G. Vidal, F. Verstraete and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
  * (17) F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066.
  * (18) F. Verstraete, J.I. Cirac, V. Murg, Adv. Phys. 57,143 (2008).
  * (19) J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. Lett. 105, 027204 (2010)
  * (20) H.-C. Jiang, Z.-C. Gu, X.-L. Qi, and S. Trebst, Phys. Rev. B 83, 245104 (2011)
  * (21) R. Schaffer, S. Bhattacharjee and Y. B. Kim, arXiv:1206.5814.
  * (22) Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad and A. Muramatsu, Nature 464, 847(2010)
  * (23) S. Sorella, Y. Otsuka and S. Yunoki, Sci. Reports 2, 992 (2012).
  * (24) H. C. Jiang, Z. Y. Weng, and T. Xiang, PRL 101, 090603 (2008).

##  .1 Supplementary

###  .1.1 Imaginary Time Evolution

After chosen the variational space of the state, i.e., mPEPS with given
parameter $D$ and $d$, we use the imaginary time evolution methodPEPS1 (16,
17, 18) to find the ground state. The imaginary time evolution could be
summarized as the following formula:

$|\psi_{g}\rangle=\lim_{t\rightarrow\infty}\frac{e^{-Ht}|\psi_{0}\rangle}{\langle\psi_{0}|e^{-Ht}e^{-Ht}|\psi_{0}\rangle}$

The $|\psi_{0}\rangle$ is the initial state, which can be chosen randomly,
after long enough imaginary time evolution, the final normalized state will be
close to the ground state with well controlled precision, as long as the
initial state is not orthogonal with the ground state (with zero measurement
to chosen these states). Generally, we can only deal with the local evolution
operators (the time evolution of a system is a global operator), and the
Trotter expansion can replace the global operators by local operators with
controllable precisionPEPS1 (16, 17, 18). For example, in Kitaev model,
Hamiltonian $H$ could be split as $H=H_{x}+H_{y}+H_{z}$, where
$H_{i},(i=x,y,z)$ is the part of Hamiltonian acting on $i$-bond. With the help
of Trotter expansion, we could split the evolution operator in second order
as:

$\displaystyle e^{-Hdt}=e^{-(H_{x}+H_{y}+H_{z})dt}$
$\displaystyle=e^{-H_{x}dt}e^{-H_{y}dt}e^{-2H_{z}dt}e^{-H_{y}dt}e^{-H_{x}dt}+O(dt^{3})$

With the Trotter expansion, it seems that we just need to act the local
evolution operators on the initial state with the given order and obtain the
ground state finally. However, the dimension of the bond $D$ will increase
exponentially with the evolution time $t$, if local operators are applied
without approximation. Therefore, as in the general TN algorithm,
approximations should be made after each step of evolution. One of the
simplest way to cut the tensors back to the initial parameter space is the
singular value decomposition (SVD), which maintain the important part of the
tensor (see Fig.A1). It is straightforward to use SVD in the case of PEPS,
while it needs extra operations in the case of mPEPS (also in fPEPS) with
fermionic statistics. Since the tensors in mPEPS are set with even parity, the
matrix in the SVD process can be re-indexed as block diagonal matrix according
to parity of the tensor’s indices (one block with odd parity and the other
with even parity), then we can apply the SVD and cut off in each of the block.

Figure A1: SVD in imaginary time evolution: (a) Act the evolution operator
$e^{-dtH}$ (the parallelogram blue tensor) on two contracted tensors. Then we
obtain a bigger tensor as indicated by (b). (c) Use SVD to cut the big tensor
and (d) reserve only D biggest singular values to make the tensors back to
their original dimension. For convenience, before any contract processes of
the tensors, we contract the diagonized matrix, which is the approximation of
the environment, into the nearest tensors.

###  .1.2 Contract

After the imaginary time evolution, we suppose that we have found the ground
state. Then the key task is to find the values of some physical quantities
which can be compared with the experiment. This task is realized by the
contract processes, which including the contract of tensors and entangled
pairs. The physical quantities in a given ground state $|\psi_{g}\rangle$ can
be expressed as

$\langle
O\rangle=\frac{\langle\psi_{g}|\hat{O}|\psi_{g}\rangle}{\langle\psi_{g}|\psi_{g}\rangle}$

For simplicity, we use $\langle\psi_{g}|\psi_{g}\rangle$ as an example to show
how to contract to calculate the physical quantities. The calculation of
physical quantities is similar. Suppose we have a 2D honeycomb lattice and
each site is labeled with two indices i,j. We have a state in the form of
projected entangled-pair State:

$|\psi\rangle=\prod_{ij}P^{[ij]}\prod_{i_{1}j_{1}i_{2}j_{2}}\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}\,,$

where
$\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}=\sum_{v}|v\rangle_{[i_{1}j_{1}]}|v\rangle_{[i_{2}j_{2}]}$
and $[i_{1}j_{1}][i_{2}j_{2}]$ denotes the edge e, linking site $[i_{1}j_{1}]$
and site $[i_{2}j_{2}]$. Normally, translational invariance of tensors will be
imposed in infinite lattice. Actually, we have the corresponding bra vector as

$\langle\psi|=\prod_{i_{1}j_{1}i_{2}j_{2}}{\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}^{{\dagger}}}\prod_{ij}{P^{[ij]{{\dagger}}}}\,,$

where $\mbox{EPR}^{{\dagger}}_{[i_{1}j_{1}][i_{2}j_{2}]}=\sum_{v}\langle
v|_{i_{2}j_{2}}\langle v|_{i_{1}j_{1}}$. The product order is reversed which
will make a difference in the contract procedure below. The scalar product of
bra vector and ket vector is:

$\displaystyle\langle\psi|\psi\rangle=\prod_{i_{1}j_{1}i_{2}j_{2}}{\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}^{{\dagger}}}\prod_{ij}{P^{[ij]{{\dagger}}}}$
$\displaystyle\prod_{ij}{P^{[ij]}}\prod_{i_{1}j_{1}i_{2}j_{2}}{\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}}$

Because we are considering mPEPS where tensors have even parity, it is free to
change the orders. We would like to rearrange the order as:

$\displaystyle\langle\psi|\psi\rangle=\prod_{i_{1}j_{1}i_{2}j_{2}}{\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}^{{\dagger}}}$
$\displaystyle\prod_{ij}{P^{[ij]{{\dagger}}}P^{[ij]}}\prod_{i_{1}j_{1}i_{2}j_{2}}{\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}}$

The product of $P^{[ij]{\dagger}}P^{[ij]}$ will contract over physical indices
and we symbolize it as $K_{ij}$.

$\displaystyle K_{ij}=P^{[ij]{\dagger}}P^{[ij]}$ $\displaystyle=\sum
T^{[ij]\ast}_{\overline{p}\overline{v_{x}}\overline{v_{y}}\overline{v_{z}}}|\overline{v_{z}}\overline{v_{y}}\overline{v_{x}}\rangle\langle\overline{p}|^{[ij]}_{pv_{x}v_{y}v_{z}}|p\rangle\langle
v_{x}v_{y}v_{z}|$ $\displaystyle=\sum
T^{[ij]\ast}_{p\overline{v_{x}}\overline{v_{y}}\overline{v_{z}}}T^{[ij]}_{pv_{x}v_{y}v_{z}}|\overline{v_{z}}\overline{v_{y}}\overline{v_{x}}\rangle\langle
v_{x}v_{y}v_{z}|$ $\displaystyle\rightarrow\sum
K^{[ij]}_{\overline{v_{x}}v_{x}\overline{v_{y}}v_{y}\overline{v_{z}}v_{z}}|\overline{v_{x}}\rangle\langle
v_{x}|\overline{v_{y}}\rangle\langle v_{y}|\overline{v_{z}}\rangle\langle
v_{z}|\,,$

where repeated indices mean contraction. Notice that
$K^{[ij]}_{\overline{v_{x}}v_{x}\overline{v_{y}}v_{y}\overline{v_{z}}v_{z}}$
is not just the contraction of
$T^{[ij]\ast}_{p\overline{v_{x}}\overline{v_{y}}\overline{v_{z}}}$ and
$T^{[ij]}_{pv_{x}v_{y}v_{z}}$, but with extra minus signs resulting from
commutation of operators.

Therefore, the scalar product is

$\langle\psi|\psi\rangle=\prod_{i_{1}j_{1}i_{2}j_{2}}{\mbox{EPR}^{{\dagger}}_{[i_{1}j_{1}][i_{2}j_{2}]}}\prod_{ij}{K_{ij}}\prod_{i_{1}j_{1}i_{2}j_{2}}{\mbox{EPR}_{[i_{1}j_{1}][i_{2}j_{2}]}}$

What we should do next is to contract over all virtual indices, ie. to product
all the $K$ tensors by EPRs. Even though we consider here is a infinite
lattice, we suppose there are boundaries in very remote areas. We contract row
by row many times from boundaries until the boundary tensors are stable. After
we contract the boundary with a row of tensors, the boundary tensors will get
thicker, like the case of imaginary time evolution, and we need to do
approximations to cut them back to their initial shape. Normally we use a
specific form of boundary tensors to make the approximation, for the sake of
simplicity (see Fig.A2 for details).

Figure A2: contract procedures: (1) Contract boundary tensors with one row of
K tensors; (2) trace the new row of tensors with their conjugate; (3) Seek
left and right environment tensors by iterative product; (4) Obtain L and R
tensor by cutting left and right environment using eigenvalue decomposition,
Contract L, the fixed tensor, and R together; (5) Use SVD to make the
approximation; (6) Contract L’s inverse and R’s inverse tensors respectively.

The main tool used to make this approximation is still SVD: We fix two
tensors, and trace over left and right parts of boundary tensors to gain the
left and right environments of the two fixed tensors. Now we contract over
left environment tensor, fixed tensors, right environment tensors together to
become a bigger tensor. We use SVD to cut the bigger tensor, then cancel the
environment effect by contracting inverse tensors of the left and right
environments. It is extremely obvious if you write down the approximation
judgement explicitly:

$\displaystyle\min_{B}(\mbox{Tr}(B-B_{0})^{{\dagger}}(B-B_{0}))$
$\displaystyle=\min_{U^{\prime}_{1}U^{\prime}_{2}}\mbox{Tr}((G_{L}U^{\prime}_{1}U^{\prime}_{2}G_{R}-G_{L}GG_{R})^{{\dagger}}(G_{L}U^{\prime}_{1}U^{\prime}_{2}G_{R}-G_{L}GG_{R}))$
$\displaystyle=\min_{U^{\prime}_{1}U^{\prime}_{2}}(\mbox{Tr}(G_{R}^{{\dagger}}U^{\prime{\dagger}}_{2}U^{\prime{\dagger}}_{1}G_{L}^{{\dagger}}G_{L}U^{\prime}_{1}U^{\prime}_{2}G_{R})$
$\displaystyle-\mbox{Tr}(G_{R}^{{\dagger}}U^{\prime{\dagger}}_{2}U^{\prime{\dagger}}_{1}G_{L}^{{\dagger}}G_{L}GG_{R})\mbox{Tr}(G_{R}^{{\dagger}}G^{{\dagger}}G_{L}^{{\dagger}}G_{L}U^{\prime}_{1}U^{\prime}_{2}G_{R})+$
$\displaystyle\mbox{Tr}(G_{R}^{{\dagger}}G^{{\dagger}}G_{L}^{{\dagger}}G_{L}GG_{R}))$
$\displaystyle=\min_{U^{\prime}_{1}U^{\prime}_{2}}(\mbox{Tr}(R^{{\dagger}}U^{\prime{\dagger}}_{2}U^{\prime{\dagger}}_{1}LU^{\prime}_{1}U^{\prime}_{2})\mbox{Tr}(R^{{\dagger}}U^{\prime{\dagger}}_{2}U^{\prime{\dagger}}_{1}LG)$
$\displaystyle-\mbox{Tr}(R^{{\dagger}}G^{{\dagger}}LU^{\prime}_{1}U^{\prime}_{2})+\mbox{Tr}(R^{{\dagger}}G^{{\dagger}}LG))\,,$

where $B_{0}$ is the old boundary contracting with a row of tensors, B is the
new boundary with only two tensors different from the old boundary. G is the
unit cell of $B_{0}$. $G_{L}$ $G_{R}$ are perspectively the left and right
part of G in the row of $B_{0}$. L is just $G_{L}^{{\dagger}}G_{L}$ and R is
just $G_{R}^{{\dagger}}G_{R}$. The above equation is just binary form of
$U_{1}^{\prime}$ and $U_{2}^{\prime}$. Obviously,
$\mbox{Tr}(B-B_{0})^{{\dagger}}(B-B_{0})$ reaches its minimum when
$\sqrt{R^{{\dagger}}}U_{1}^{\prime}U_{2}^{\prime}\sqrt{L}$ best approximates
$\sqrt{R^{{\dagger}}}G\sqrt{L}$. See also Fig.A3.

This example of $\langle\psi|\psi\rangle$ showed our method for contract. It
is similar to calculate average value of local operators
$\langle\psi|\hat{O}|\psi\rangle$, except that the sites where operators lie
should be kept until final contract.

There is one more point worth mentioning: When we use mPEPS, we choose the
virtual annihilation and creation operators to be anti-commutative. So we must
be extremely careful about the minus sign when we contract over virtual
indices. Similarly and more easily to overlook, in SVD, minus signs also
appear when we change the order of virtual indices in order to obtain the
original shape of tensors. It could be explained as equations below. We take a
simplified version as an example. T is a big tensor which needs to be split by
SVD.

$\displaystyle T_{\bar{u}u,\bar{v}v}$
$\displaystyle=\sum_{\bar{s}s}U^{\prime}_{\bar{u}u,\bar{s}s}S_{\bar{s}s}V^{\prime}_{\bar{s}s,\bar{v}v}$
$\displaystyle=\sum_{\bar{s}s}U_{\bar{u}u,\bar{s}s}V_{\bar{s}s,\bar{v}v}$
$\displaystyle=\sum_{\overline{s1}s1\overline{s2}s2}U_{\bar{u}u,\overline{s1}s1}\delta_{\overline{s1}s1,\overline{s2}s2}V_{\overline{s2}s2,\bar{v}v}$
$\displaystyle=\overline{\mbox{EPR}}^{{\dagger}}(\sum_{\overline{s1}s1\overline{s2}s2}U_{\bar{u}u,\overline{s1}s1}|\overline{s2}\overline{s1}\rangle\langle
s1s2|_{\overline{s2}s2,\bar{v}v})\mbox{EPR}$
$\displaystyle=\overline{\mbox{EPR}}^{{\dagger}}(\sum_{\overline{s1}s1\overline{s2}s2}U_{\bar{u}u,\overline{s1}s1}(-)^{(p(\overline{s1})+p(s1))p(\overline{s2})}$
$\displaystyle|\overline{s1}\rangle\langle s1||\overline{s2}\rangle\langle
s2|V_{\overline{s2}s2,\bar{v}v})\mbox{EPR}$
$\displaystyle=\overline{\mbox{EPR}}^{{\dagger}}(\sum_{\overline{s1}s1\overline{s2}s2}U_{\bar{u}u,\overline{s1}s1}(-)^{(p(\overline{s1})+p(s1))p(\overline{s1})}$
$\displaystyle|\overline{s1}\rangle\langle s1||\overline{s2}\rangle\langle
s2|V_{\overline{s2}s2,\bar{v}v})\mbox{EPR}$
$\displaystyle=\overline{\mbox{EPR}}^{{\dagger}}(\sum_{\overline{s1}s1\overline{s2}s2}U^{\prime\prime}_{\bar{u}u,\overline{s1}s1}$
$\displaystyle|\overline{s1}\rangle\langle s1||\overline{s2}\rangle\langle
s2|V_{\overline{s2}s2,\bar{v}v})\mbox{EPR}$ $\displaystyle=contract\
U^{\prime\prime}\ tensor\ with\ V\ tensor$

The first equation is normal SVD; the second is to absorb singular values; and
rest equations are to extract EPR from tensors. EPR is
$\sum_{s2=s1}|s2\rangle|s1\rangle$, correspondingly
$\mbox{EPR}^{{\dagger}}=\sum_{s2=s1}\langle s1|\langle s2|$, and so as their
Hermitian conjugate. In the equations, you can see that a minus sign
$(-)^{(p(\bar{s1})+p(s1))p(\bar{s1})}$ show up. Note that $p(\bar{s2})$ in the
fifth equation is replaced by $p(\bar{s1})$ since it is not zero only when
$s1=s2$ and $\bar{s1}=\bar{s2}$. p(*) is the function that maps quantum
numbers to parity. U and V is the tensors that $U^{\prime}$ and $V^{\prime}$
absorb square of the singular values respectively. $U^{\prime\prime}$ is the
tensor after U absorbs $(-)^{(p(\bar{s1})+p(s1))p(\bar{s1})}$.

Figure A3: The scheme of approximation in the language of variational method.
Use the second row of boundary tensors to approximate the first row.
Substitute the first row with the third row. And keep on this procedures until
stable.

###  .1.3 Additional Calculation of the Kitaev Model

In order to shows the stability of our results, we increase the dimension of
the truncation n the contract from $D_{cut}=3^{2}$ to $D_{cut}=7^{2}$. We can
see in the following figure that when $D_{cut}>4^{2}$ with $D=8$, the NN
correlation of the system is relatively stable.

Figure A4: The diagram shows that when $D=8$, the NN correlations
$\langle\sigma^{z}\sigma^{z}\rangle$ are stable as $D_{cut}$ increases from
$D_{cut}=4^{2}$ to $D_{cut}=7^{2}$.

With the stable algorithm, we calculate the energy of the Kitaev model which
is well agree with the exact energy

Figure A5: The diagram shows the ground state energy of Kitaev model by
fPEPS.(ST stands for exact results. FS stands for the results after finite
scaling of D.)

Since the limit of our calculation, we should do some finite scaling. However,
the finite scaling of D is not clear and well defined. Here we use linear of
polynomial fit with $1/D$, as showed below.

Figure A6: Finite D Scaling of $\langle\sigma^{z}\sigma^{z}\rangle_{z}$ at
Jz=0.4. We choose it to linear-fit with $1/D$. Figure A7: Khaliullin
Transformation: The spins on square sites are fixed. The spins on circle sites
are rotated by $\pi$ along x bond. The spins on triangle sites are rotated by
$\pi$ along y bond. The spins on pentagon sites are rotated by $\pi$ along z
bond.

###  .1.4  Khaliullin Transformation for Kitaev-Heisenberg Model

This transformation first proposed by Khaliullin and Chaloupka not only makes
the anti-ferromagnetic state to ferromagnetic state, but also makes Kitaev-
Heisenberg model exactly solvable at $\alpha=1/2$. This transformation
requires different spins in the lattice to rotate about different axis. To be
specific, we first choose a set of spins which are positioned on third nearest
neighbor sites at opposite corners of the hexagons throughout the lattice, and
hold these spins fixed. We next rotate the three nearest spins by $\pi$ about
the spin axis corresponding to the bond which connects it to the fixed spin.
The effect of this operation therefore transforms the Heisenberg Hamiltonian
to a mixture of Heisenberg Hamiltonian and Kitaev Model, and Kitaev
Hamiltonian invariant.

$\displaystyle H_{H}\rightarrow-H_{H}+2H_{K}$ $\displaystyle H_{K}\rightarrow
H_{K}$

Thus the Kitaev-Heisenberg model is transformed as:

$H\rightarrow-(1-\alpha)H_{H}-4(\alpha-1/2)H_{K}$ (A1)

It is obvious to find out in the transformed Hamiltonian, when $\alpha=1/2$,
it is just $H=-1/2H_{H}$ Heisenberg Hamiltonian with a negative coupling
coefficient. So the state at $\alpha=1/2$ after the transformation is a
ferromagnetic state. We can operate an inverse transformation to send the
ferromagnetic state back to an anti-ferromagnetic state. As we have explained
before, the anti-ferromagnetic state is harder to calculate, since stable
results are not easily required. So in the letter, we just use the Hamiltonian
after the transformation to discuss the quantum phase diagram of Kitaev-
Heisenberg model.

In order to weaken the influence by Heisenberg interaction and make the spin
liquid region broader and easier to see, we re-parameterize the Kitaev-
Heisenberg model to:

$H\rightarrow-\frac{1-\alpha^{\prime}}{3}H_{H}-4(\alpha^{\prime}-1/2)H_{K}$
(A2)

The new parameter $\alpha^{\prime}$ and the old one $\alpha$ can be related by
a simple formula, if $H_{K}$ is proportional to $H_{H}$ to the same extent.
Besides, the new Hamiltonian Eq.A2 preserves the same behaviors as Eq.A1 at
$\alpha^{\prime}=1$ and $\alpha^{\prime}=1/2$. In the paper, we discuss Eq.A2,
and note that the spin liquid region of Eq.A1 is actually much much narrower.
The $\alpha^{\prime}$ region of spin liquids is $0.98<\alpha^{\prime}<1$, as
the part below suggests. Correspondingly, the real parameter region of spin
liquids before re-parameterization are around $0.997\langle\alpha\langle 1$.
In the paper, we just use the new Hamiltonian Eq.A2 to further our discussion.

The configuration of the ground state of the Kitaev-Heisenberg model with
$\alpha=0.5$ can be calculated by our PEPS method as Fig.A8

Figure A8: The above figure shows the spin orientation of $\alpha=1/2$ Kitaev-
Heisenberg model after Kahliullin transformation. And with the inverse
transformation, the spin pattern returns to stripy anti-ferromagnetic order.