# Data Science Kitchen at GermEval 2021: A Fine Selection of Hand-Picked
Features, Delivered Fresh from the Oven

Niclas Hildebrandt, Benedikt Boenninghoff, Dennis Orth, Christopher Schymura
Data Science Kitchen
<EMAIL_ADDRESS>
###### Abstract

This paper presents the contribution of the Data Science Kitchen at GermEval
2021 shared task on the identification of toxic, engaging, and fact-claiming
comments. The task aims at extending the identification of offensive language,
by including additional subtasks that identify comments which should be
prioritized for fact-checking by moderators and community managers. Our
contribution focuses on a feature-engineering approach with a conventional
classification backend. We combine semantic and writing style embeddings
derived from pre-trained deep neural networks with additional numerical
features, specifically designed for this task. Ensembles of Logistic
Regression classifiers and Support Vector Machines are used to derive
predictions for each subtask via a majority voting scheme. Our best submission
achieved macro-averaged F1-scores of 66.8%, 69.9% and 72.5% for the
identification of toxic, engaging, and fact-claiming comments.

##  1 Introduction

In the early years after establishing social media platforms, setting up
online discussion forums and installing comment areas on newspapers’ websites,
a door into this new digital world has been opened, allowing people to
interconnect all over the world. Various communication platforms and social
networks enabled new ways of sharing information with followers, exchanging
opinions between politically interested people, and encouraging debates with
the readers. Unfortunately, recent trends revealed the ugly face and adverse
effects of these platforms when an increasing number of users make improper,
illegal, or abusive use of such digital services Mathew et al. (2019).

Nowadays, social media platforms are notorious for spreading toxic comments,
in which the writers justify violence and discrimination against a person or
groups of persons Munn (2020). Additionally, a second steadily growing trend
is producing and sharing fake news or misinformation, seeking to dominate
current discussions, and frame public debates Mahid et al. (2018).

Both hate speech, fake news and their impact have become very prominent in
recent years. However, the tremendous amount of shared and distributed toxic
messages on social media platforms make it utterly infeasible to identify and
tag or delete poisonous comments manually. The GermEval 2021 shared task tries
to encounter this negative trend and motivates participants to work on
automated solutions towards safer and more reliable digital rooms of
interaction Risch et al. (2021).

Therefore, the organizers of the task increased the difficulty of the
competition by expanding the focus not only on the identification of toxic
messages in online discussions but also on distinguishing between engaging and
fact-claiming comments. The first task is similar to the GermEval tasks in
2018 Wiegand et al. (2019) and 2019 Strus̈ et al. (2019) and deals with
identifying toxic comments, including offensive, hateful and vulgar language
or ruthless cynism. As novel subtasks, the participants are also invited to
identify two additional categories of comments: The second category defines
engaging comments, which are annotated as highly relevant contributions by the
moderators. The third category concentrates on finding fact-claiming comments
that should be considered for a manual fact-check with a higher priority.

##  2 Task and Data Description

Each subtask of GermEval 2021 is defined as a binary classification problem
and all tasks share the same training and test data. The set of training data
consists of 3,244 Facebook comments from a German news broadcast page. The
anonymized comments were posted in the time span from February to July 2019
and were labeled by trained annotators. Binary labels were provided for each
of the three categories. The test data is also extracted from Facebook
discussions and include 944 comments. However, these comments had a different
discussion topic than the training data. Precision, recall, and macro-averaged
F1-score were defined as the relevant evaluation metrics.

###  2.1 Subtask 1: Toxic Comment Classification

Toxic comments are characterized by their offensive and hateful language,
intended to blame other people or groups. For social media and content
providers, it is important to detect such comments in a highly automated and
scalable way. An example of a toxic comment from the training data of the
GermEval shared task is: “Na, welchem tech riesen hat er seine Eier
verkauft..?”. However, some of the comments which have been labeled as toxic
can be quite hard to detect. Examples of such cases are: “@USER eididei sieh
mal an” or “ein schöner VW Golf Diesel..”. Difficulties occur due to irony,
subtle overtones and missing contextual information.

###  2.2 Subtask 2: Engaging Comment Classification

Engaging comments encourage other users to join the discussion, express their
opinions and share ideas regarding the topic. They are characterised by being
rational, respectful, and reciprocal and hence can foster a constructive and
fruitful discussion. The comment “Wie wär’s mit einer Kostenteilung.
Schließlich haben beide Parteien (Verkäufer und Käufer) etwas von der
Tätigkeit des Maklers. Gilt gleichermassen für Vermietungen. Die Kosten werden
so oder soweit verrechnet, eine Kostenreduktion ist somit nicht zu erwarten.”
is an example of an engaging comment from the training data.

###  2.3 Subtask 3: Fact-Claiming Comment Classification

If a platform provider has to prevent the spread of fake news and
misinformation, there is the demand of automatically identifying fact-claiming
comments to assess their truthfulness. An example of a fact-claiming comment
from the training data is the comment “Dummerweise haben wir in der EU und in
der USA einen viel höheren CO2 Fußabdruck als z.B. die Afrikaner oder Inder.”.

##  3 System Overview

The general system architecture is shown in Fig. 1. As the number of samples
in the training dataset provided for GermEval 2021 is rather small, our
proposed framework focuses on suitable feature engineering with a conventional
classification backend. These features and further implementation details of
our system are described in the following.

Figure 1: General architecture of the proposed framework to detect toxic,
engaging and fact-claiming comments. Yellow boxes denote non-trainable,
computational functions and transformations, green boxes represent pretrained
models utilized for feature extraction. Trainable models whose parameters are
optimized using the challenge dataset are shown in blue.

###  3.1 Preprocessing

Raw input text is preprocessed in three different processing streams that are
handled in parallel. The first stream utilizes the tokenizer of a German BERT
model Chan et al. (2020) and crops the corresponding input text at a maximum
length of 512 tokens. The second stream uses the SoMaJo tokenizer Proisl and
Uhrig (2016) for German language and the third stream passes the raw text to
the subsequent feature extraction stage without any preprocessing.

###  3.2 Feature Extraction

The feature extraction stage focuses on embedding-based features, as well as
manually selected, numerical feature representations. Specific feature types
are computed using one of the three preprocessing streams described in Sec.
3.1. This specific feature extraction setup was chosen to efficiently combine
embedding representations that capture linguistic properties with “hand-
crafted” features specifically designed for the GermEval 2021 tasks.

####  3.2.1 Semantic Embeddings

The first kind of embeddings used in our framework are document embeddings
derived from a pre-trained German BERT model Chan et al. (2020). Specifically,
we used the bert-base-german-cased implementation from
Huggingface111https://huggingface.co/dbmdz/bert-base-german-cased. This model
was trained on a German Wikipedia dump, the OpenLegalData dump Ostendorff et
al. (2020) and news articles. Average pooling was used to compute
768-dimensional document-level embeddings from the BERT model output.

####  3.2.2 Writing Style Embeddings

Besides semantic document embeddings, we additionally experimented with neural
stylometric embeddings that have been automatically extracted from the
comments. More precisely, we used an extended framework of AdHominem
Boenninghoff et al. (2019) which outperformed all other systems that
participated in the PAN 2020 and 2021 authorship verification tasks
Boenninghoff et al. (2021).

The overall framework consists of three components: In a first step, we
perform neural feature extraction and deep metric learning (DML) to encode the
writing style characteristics of a pair of raw documents into a pair of fixed-
length representations, which is realized in the form of a Siamese network.
Inspired by Hochreiter and Schmidhuber (1997), the Siamese network consists of
a hierarchical LSTM-based topology. Next, the obtained representations are fed
into a Bayes factor scoring (BFS) layer to compute the posterior probability
for this trial. The idea of this second component is to take into account
both, the similarity between the questioned documents and the typicality
w.r.t. the relevant population represented by the training data. The third
component is given by an uncertainty adaptation layer (UAL) aiming to correct
possible misclassifications and to return corrected and calibrated posteriors.
More details can be found in Boenninghoff et al. (2021).

To train the model, we prepared a large dataset of Zeit-Online forum
comments222www.zeit.de. Altogether, we collected 9,812,924 comments written by
204,779 authors. Afterwards we split the dataset into training and validation
sets. We took 10% of the authors to build the validation set and removed all
comments with less than $60$ tokens. Due to the fact that the provided dataset
of the shared task also contains concise comments, we decided to leave all
short comments in the training set. As a result, the datasets are disjoint
w.r.t. the authors, i.e., all authors in the validation set have been removed
from the training set. During training, we perform data augmentation by
resampling new same-author and different-authors pairs in each epoch.
Contrary, the pairs of the validation set are sampled once and then kept
fixed. Since some authors contribute with hundreds of comments, we limited
their influence by sampling not more than 20 comments per author. In summary,
the training set contains approximately 234,500 same-author and 244,200
different-authors pairs in each epoch, where, on average, each comment
consists of $75.90\pm 68.07$ tokens. The validation set contains 15,125 same-
author and 18,740 different-authors pairs, where, on average, each comment
consists of $126.51\pm 65.31$ tokens. Hence, both datasets are nearly
balanced.

Table 1: Results for PAN 2021 evaluation metrics.

Model | PAN 2021 Evaluation Metrics
---|---
AUC | c@1/acc | f_05_u | F1 | Brier | Overall
DML | $87.4$ | $79.3$ | $81.7$ | $81.0$ | $85.1$ | $82.9$
BFS | $87.4$ | $79.5$ | $80.5$ | $82.0$ | $85.5$ | $83.0$
UAL | $87.6$ | $79.5$ | $81.6$ | $81.4$ | $85.6$ | $83.2$

Table 2: Overview of all features utilized in this work that are not based on
embeddings.

Feature name | Dim. | Description
---|---|---
NumCharacters | 1 | Total number of characters, including white spaces.
NumTokens | 1 | Total number of tokens, after splitting at white spaces.
AverageTokenLength | 1 | Average number of characters in all tokens.
TokenLengthStd | 1 | Standard deviation of the number of characters in all tokens.
StopwordRatio | 1 | Number of stop words divided by the number of tokens.
ExclamationMarkRatio | 1 | Number of exclamation marks divided by the number of characters.
NumReferences | 1 | Number of hyperlinks in the comment.
NumMediumAdressed | 1 | Number of @MEDIUM mentions in the comment.
NumUserAdressed | 1 | Number of @USER mentions in the comment.
AverageEmojiRepetition | 1 | Average repetition number of emojis used in the comment.
SpellingMistakes | 17 | Number of specific grammar and spelling mistakes, cf. Sec. 3.2.3.
SentimentBERT | 3 | Sentiment scores of a pre-trained BERT model Guhr et al. (2020).

We choose the PAN 2021 evaluation metrics to evaluate the performance as
described in Kestemont et al. (2021). Table. 1 summarizes the results, where
all three system components are evaluated separately. It can be seen that we
achieved overall scores between $82.9$ and $83.2$ for the components, which is
mainly supported by higher values for the AUC and Brier scores. Comparing the
c@1, F1 or f_05_u metrics, we generally obtained error rates of approximately
$20\%$ on this challenging dataset for a fixed threshold. After training, one
part of the neural feature extraction component within the Siamese network is
then used to extract the 100-dimensional writing style embeddings for the
shared task data.

####  3.2.3 Additional Numerical Features

In addition to the semantic and writing style embeddings, we integrated a set
of specifically designed numerical features into our framework. An overview of
these features, their dimensionality and corresponding descriptions is given
in Tab. 2. We applied the natural logarithm to all strictly-positive numerical
features.

The first group of features, NumCharacters, NumTokens, AverageTokenLength and
TokenLengthStd, were chosen to reflect general structural properties of the
comments in the dataset. In addition, we use the StopwordRatio and
ExclamationMarkRatio features to explicitly reflect task-related semantic
properties in the dataset. These task-specific features are accompanied by
additional count-based features NumMediumAdressed, NumUserAddressed,
NumReferences and AverageEmojiRepetition. We also included the scores
(corresponding to the classes “positive”, “neutral” and “negative”) of a BERT
model for sentiment classification trained on 1,834 million German-language
samples derived from various sources Guhr et al. (2020) as a dedicated
SentimentBERT feature.

Lastly, we included an 17-dimensional feature denoted as SpellingMistakes into
our set of additional features. This feature represents spelling and grammar
mistakes from 17 different categories. We used a Python wrapper from the open-
source grammar checker LanguageTool333https://languagetool.org/ to derive this
feature. In particular, the following classes of mistakes were considered:
Typography, punctuation, grammar, upper/lowercase, support in punctuation,
colloquialism, compounding, confused words, redundancy, typos, style, proper
nouns, idioms, recommended spelling, miscellaneous, double punctuation, double
exclamation mark. For every category, we counted the number of mistakes and
divided them by the number of tokens in the respective comment.

###  3.3 Classification Pipeline

The classification pipeline used in this work is depicted in Fig. 1. The
semantic and writing style embedding features described in Secs. 3.2.1 and
3.2.2 are concatenated, yielding a 868-dimensional joint embedding vector. A
truncated singular-value decomposition (SVD) Halko et al. (2011) is applied on
this vector to reduce its dimensionality for subsequent processing. The number
of dimensions kept is treated as a hyperparamter during training, cf. Sec. 4.
The reduced joint embedding vector is then concatenated with the
28-dimensional vector of additional numerical features. The resulting vector
is standardized to zero-mean and unit variance and serves as input to the
classification stage.

We use Logistic Regression Berkson (1944); Haggstrom (1983) and Support Vector
Machines (SVMs) Boser et al. (1992) with radial basis function (RBF) kernel as
base classifiers within individual ensembles. One ensemble of binary
classifiers is utilized for each subtask. Each classifier in the ensembles is
trained using a subset of the provided training data via a cross-validation
setup, cf. Sec. 4. A hard majority-voting scheme is used in each ensemble to
obtain the predicted labels.

(a) Subtask 1

(b) Subtask 2

(c) Subtask 3

Figure 2: Reliability diagrams of our first submission for all three subtasks
(see Section 4.1). The red bars are drawn darker for bins with a higher number
of samples.

##  4 Evaluation

Our framework is trained using a specific cross-validation and hyperparameter
tuning scheme, which is described in the following.

###  4.1 Evaluation Metrics

Precision, recall and macro average F1-score are used for model evaluation
since they represent the evaluation metrics of the GermEval 2021 shared task.
Additionally, we assess the calibration properties of our model by determining
the expected calibration error (ECE) as well as the maximum calibration error
(MCE), where the confidence interval is discretized into a fixed number of $M$
bins Naeini et al. (2015). The ECE is then computed as the weighted macro-
averaged absolute error between confidence and accuracy of all bins,

$\mathrm{ECE}=\sum_{m=1}^{M}\frac{|B_{m}|}{N}\left|\mathrm{acc}(B_{m})-\mathrm{conf}(B_{m})\right|,$
(1)

where $N$ is the total number of samples and
$\mathrm{acc}(B_{m})-\mathrm{conf}(B_{m})$ is the difference between the
actual accuracy and classifier confidence within a fixed-size bin $B_{m}$ in
the confidence interval. Note that all confidence values lie within the
interval $[0.5,1]$,since we are dealing with binary classification tasks.
Hence, to obtain confidence scores, the output predictions $p$ are transformed
w.r.t. to the estimated subtask label, showing $\mathrm{conf}=p$ if the
$\mathrm{acc}>=0.5$ and $\mathrm{conf}=1-p$ if $\mathrm{acc}<0.5$. The MCE
returns the maximum absolute error, given as

$\mathrm{MCE}=\max_{m\in{1,...,M}}\left|\mathrm{acc}(B_{m})-\mathrm{conf}(B_{m})\right|.$
(2)

We further display the reliability diagrams in Fig. 2 which will be discussed
in Section 5.

###  4.2 Experimental Setup

Table 3: Final submission results on the test set including the calibration
metrics for the first submission.

Run | Subtask 1 | Subtask 2 | Subtask 3
---|---|---|---
P | R | F1 | ECE | MCE | P | R | F1 | ECE | MCE | P | R | F1 | ECE | MCE
Submission 1 | $65.95$ | $63.67$ | $64.79$ | $5.5$ | $~{}~{}8.0$ | $69.70$ | $67.78$ | $68.72$ | $6.9$ | $10.4$ | $73.25$ | $71.44$ | $72.34$ | $3.5$ | $6.0$
Submission 2 | $64.89$ | $62.71$ | $63.78$ | $-$ | $-$ | $69.26$ | $67.43$ | $68.33$ | $-$ | $-$ | $73.39$ | $71.52$ | $72.44$ | $-$ | $-$
Submission 3 | $66.98$ | $66.73$ | $66.85$ | $-$ | $-$ | $71.71$ | $68.34$ | $69.98$ | $-$ | $-$ | $73.03$ | $72.08$ | $72.55$ | $-$ | $-$

Our experimental setup involves dedicated model selection and hyperparameter
tuning. The training set performance is evaluated in a stratified $K$-Fold
cross validation setup preserving the class label distribution among all
folds. One of the $K$ folds is used as validation set. We utilized a 7-fold
cross-validation scheme and computed the evaluation metrics described in Sec.
4.1 on the validation set of each fold.

For submission one and two there are 7 logistic regression models for each
subtask trained on different folds and stacked together in a voting ensemble
returning the prediction of the majority. On each fold the L2-regularisation
strength C and the number of features coming from the SVD dimension reduction
are tuned with respect to the macro averaged F1 score over all subtasks. This
means that hyperparameters may be slightly different from fold to fold but all
three models trained on the same fold get the same hyperparameters –
regardless the classification task.

Submission three uses a similar approach but the logistic regression models
are replaced by SVMs having the same fold-wise hyperparameter tuning as
mentioned above. In addition, task-wise tuned SVMs are added to the ensemble.
Doubling the number of models and including a higher level of customisation to
the task. The task-wise tuning includes optimisation of kernel, L2
regularisation strength C, class weight(whether or not to weight C with the
class label distribution) and the kernel coefficient gamma as defined in the
sklearn library Pedregosa et al. (2011).

Hyperparameter tuning is performed with Bayesian optimisation using the Optuna
library Akiba et al. (2019). The macro average F1-score is chosen as
optimisation target and the best hyperparameters among 100 trails are used in
the ensemble.

(a) Subtask 1

(b) Subtask 2

(c) Subtask 3

Figure 3: Confusion matrices for submission 3.

##  5 Results and Discussion

The final submission results are provided in Table 3. Unexpectedly, the
identification of toxic comments turns out to be the most challenging subtask
while the detection of fact-claiming comments achieved the highest F1-score.
This confirms our observations during hyperparameter tuning. For instance, the
F1-score for the third submission after cross validation are given by
$66.31\pm 1.76$, $75.12\pm 2.07$ and $74.68\pm 2.67$ for subtasks 1-3,
respectively. A comparison of our cross validation performance with the
results on the test set shows two interesting findings: On the one side, we
obtained very robust results of subtasks 1 and 3. On the other side, subtask 2
struggles with over-fitting effects.

In addition, Fig. 3 displays the confusion matrices of our third submission
(representative for all submissions). It can be seen for all subtasks that the
ratio of wrongly classified positively labeled samples is significantly larger
than for negatively labeled samples. This behavior is supported by the
reliabilty diagrams444https://github.com/hollance/reliability-diagrams in Fig.
2, where our submission delivers over-confident scores (i.e. conf $>$ acc) in
nearly all bins. As a results, the higher proportion of wrongly classified
comments for positively labeled comments leads to a lower performance in terms
of the F1-score.

Finally, we visualize an estimated probability density function of the first
submission using a non-parametric Gaussian kernel density
estimator555https://scikit-
learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html in
Fig. 4. Ideally, we would expect a bimodal probability density function.
However, the plot shows that the system clearly tends towards self-confident
predictions close to zero. But in regions closer to one, the systems behave
more hesitant. This effect can be explained by the imbalanced distribution of
the class labels.

Figure 4: Gaussian kernel density estimates for the distributions of the first
submission (bandwidth$=0.08$).

##  6 Conclusions

Within this contribution to the shared task of the GermEval 2021 we have
developed a modular feature extraction scheme which incorporates semantic and
writing style embeddings as well as task specific numerical features. Less
complex algorithms like logistic regression models and SVMs converge converge
faster than complex models like deep neural networks and therefore need less
training data. The combination with automated hyperparameter tuning and
dimension reduction as well as the final agglomeration of multiple models in
voting ensembles allow to achieve an macro-averaged F1-scores of 66.8%, 69.9%
and 72.5% for the identification of toxic, engaging, and fact-claiming
comments.

## References

  * Akiba et al. (2019) Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. 2019\.  Optuna: A next-generation hyperparameter optimization framework.  In _Proc. of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 2623–2631, New York, NY, USA. Association for Computing Machinery.
  * Berkson (1944) Joseph Berkson. 1944.  Application of the logistic function to bio-assay.  _Journal of the American Statistical Association_ , 39(227):357–365.
  * Boenninghoff et al. (2019) Benedikt Boenninghoff, Steffen Hessler, Dorothea Kolossa, and Robert M. Nickel. 2019\.  Explainable authorship verification in social media via attention-based similarity learning.  In _2019 IEEE International Conference on Big Data_ , pages 36–45.
  * Boenninghoff et al. (2021) Benedikt Boenninghoff, Dorothea Kolossa, and Robert M. Nickel. 2021.  Self-calibrating neural-probabilistic model for authorship verification under covariate shift.  In _12th International Conference of the CLEF Association (CLEF 2021)_. Springer.
  * Boser et al. (1992) Bernhard E. Boser, Isabelle M. Guyon, and Vladimir N. Vapnik. 1992.  A training algorithm for optimal margin classifiers.  In _Proceedings of the Fifth Annual Workshop on Computational Learning Theory_ , page 144–152. Association for Computing Machinery.
  * Chan et al. (2020) Branden Chan, Stefan Schweter, and Timo Möller. 2020.  German’s next language model.  In _Proceedings of the 28th International Conference on Computational Linguistics_ , pages 6788–6796. International Committee on Computational Linguistics.
  * Guhr et al. (2020) Oliver Guhr, Anne-Kathrin Schumann, Frank Bahrmann, and Hans Joachim Böhme. 2020\.  Training a broad-coverage German sentiment classification model for dialog systems.  In _Proceedings of the 12th Language Resources and Evaluation Conference_ , pages 1627–1632. European Language Resources Association.
  * Haggstrom (1983) Gus W. Haggstrom. 1983.  Logistic regression and discriminant analysis by ordinary least squares.  _Journal of Business & Economic Statistics_, 1(3):229–238.
  * Halko et al. (2011) N. Halko, P. G. Martinsson, and J. A. Tropp. 2011.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  _SIAM Review_ , 53(2):217–288.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long short-term memory.  _Neural Computation_ , 9(8):1735–1780.
  * Kestemont et al. (2021) Mike Kestemont, Efstathios Stamatatos, Enrique Manjavacas, Janek Bevendorff, Martin Potthast, and Benno Stein. 2021.  Overview of the authorship verification task at pan 2021.  In _CLEF 2021 Labs and Workshops, Notebook Papers_. CEUR-WS.org.
  * Mahid et al. (2018) Zaitul Iradah Mahid, Selvakumar Manickam, and Shankar Karuppayah. 2018.  Fake news on social media: Brief review on detection techniques.  In _2018 Fourth International Conference on Advances in Computing, Communication Automation (ICACCA)_ , pages 1–5.
  * Mathew et al. (2019) Binny Mathew, Ritam Dutt, Pawan Goyal, and Animesh Mukherjee. 2019.  Spread of hate speech in online social media.  In _Proceedings of the 10th ACM Conference on Web Science_ , page 173–182. Association for Computing Machinery.
  * Munn (2020) Luke Munn. 2020.  Angry by design: Toxic communication and technical architectures.  _Humanities and Social Sciences Communications_ , 7(53):229–238.
  * Naeini et al. (2015) Mahdi Pakdaman Naeini, Gregory F. Cooper, and Milos Hauskrecht. 2015.  Obtaining well calibrated probabilities using bayesian binning.  In _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ , page 2901–2907. AAAI Press.
  * Ostendorff et al. (2020) Malte Ostendorff, Till Blume, and Saskia Ostendorff. 2020.  Towards an open platform for legal information.  In _Proceedings of the ACM/IEEE Joint Conference on Digital Libraries in 2020_ , page 385–388. Association for Computing Machinery.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay. 2011.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12(85):2825–2830.
  * Proisl and Uhrig (2016) Thomas Proisl and Peter Uhrig. 2016.  SoMaJo: State-of-the-art tokenization for German web and social media texts.  In _Proceedings of the 10th Web as Corpus Workshop (WAC-X) and the EmpiriST Shared Task_ , pages 57–62. Association for Computational Linguistics (ACL).
  * Risch et al. (2021) Julian Risch, Anke Stoll, Lena Wilms, and Michael Wiegand. 2021.  Overview of the GermEval 2021 shared task on the identification of toxic, engaging, and fact-claiming comments.  In _Proceedings of the GermEval 2021 Shared Task on the Identification of Toxic, Engaging, and Fact-Claiming Comments co-located with KONVENS_ , pages 1–12.
  * Strus̈ et al. (2019) Julia Maria Strus̈, Melanie Siegel, Josef Ruppenhofer, Michael Wiegand, and Manfred Klenner. 2019.  Overview of GermEval task 2, 2019 shared task on the identification of offensive language.  pages 352–363. German Society for Computational Linguistics & Language Technology.
  * Wiegand et al. (2019) Michael Wiegand, Melanie Siegel, and Josef Ruppenhofer. 2019.  Overview of the germeval 2018 shared task on the identification of offensive language.  pages 1–10. Austrian Academy of Sciences.


* [399] J.-P. Lasota, P. Haensel, M. A. Abramowicz, Fast rotation of neutron stars, Astrophys. J. 456 (1996) 300.  arXiv:astro-ph/9508118, doi:10.1086/176650.
  * [400] C. Breu, L. Rezzolla, Maximum mass, moment of inertia and compactness of relativistic stars, Mon. Not. Roy. Astron. Soc. 459 (1) (2016) 646–656.  arXiv:1601.06083, doi:10.1093/mnras/stw575.
  * [401] T. W. Baumgarte, S. L. Shapiro, M. Shibata, On the maximum mass of differentially rotating neutron stars, Astrophys. J. Lett. 528 (2000) L29.  arXiv:astro-ph/9910565, doi:10.1086/312425.
  * [402] D. Gondek-Rosińska, I. Kowalska, L. Villain, M. Ansorg, M. Kucaba, A new view on the maximum mass of differentially rotating neutron stars, Astrophys. J. 837 (1) (2017) 58.  arXiv:1609.02336, doi:10.3847/1538-4357/aa56c1.
  * [403] L. R. Weih, E. R. Most, L. Rezzolla, On the stability and maximum mass of differentially rotating relativistic stars, Mon. Not. Roy. Astron. Soc. 473 (1) (2018) L126–L130.  arXiv:1709.06058, doi:10.1093/mnrasl/slx178.
  * [404] J. W. T. Hessels, S. M. Ransom, I. H. Stairs, P. C. C. Freire, V. M. Kaspi, F. Camilo, A radio pulsar spinning at 716-hz, Science 311 (2006) 1901–1904.  arXiv:astro-ph/0601337, doi:10.1126/science.1123430.
  * [405] J. E. Grindlay, A Neutron Star in F-sharp, Science 311 (2006) 1876.  arXiv:astro-ph/0605117, doi:10.1126/science.1125125.
  * [406] P. S. Koliogiannis, C. C. Moustakidis, Effects of the equation of state on the bulk properties of maximally-rotating neutron stars, Phys. Rev. C 101 (1) (2020) 015805.  arXiv:1907.13375, doi:10.1103/PhysRevC.101.015805.
  * [407] Y. Fujii, Dilaton and possible non-newtonian gravity, Nature 234 (1971) 5–7.
  * [408] Y. Fujii, On the Theoretical Background of the Fifth Force, Int. J. Mod. Phys. A 6 (1991) 3505–3557.  doi:10.1142/S0217751X91001714.
  * [409] M. I. Krivoruchenko, F. Šimkovic, A. Faessler, Constraints for weakly interacting light bosons from existence of massive neutron stars, Phys. Rev. D 79 (12) (2009) 125023.  arXiv:0902.1825, doi:10.1103/PhysRevD.79.125023.
  * [410] R. Machleidt, The Meson theory of nuclear forces and nuclear structure, Adv. Nucl. Phys. 19 (1989) 189–376.
  * [411] D.-H. Wen, B.-A. Li, L.-W. Chen, Super-soft symmetry energy encountering non-Newtonian gravity in neutron stars, Phys. Rev. Lett. 103 (2009) 211102.  arXiv:0908.1922, doi:10.1103/PhysRevLett.103.211102.
  * [412] D.-R. Zhang, P.-L. Yin, W. Wang, Q.-C. Wang, W.-Z. Jiang, Effects of a weakly interacting light U boson on the nuclear equation of state and properties of neutron stars in relativistic models, Phys. Rev. C 83 (2011) 035801.  arXiv:1103.2403, doi:10.1103/PhysRevC.83.035801.
  * [413] H. Zheng, L.-W. Chen, Effects of the U-boson on the inner edge of neutron star crusts, Phys. Rev. D 85 (2012) 043013.  arXiv:1111.0883, doi:10.1103/PhysRevD.85.043013.
  * [414] W. Lin, B.-A. Li, L.-W. Chen, D.-H. Wen, J. Xu, Breaking the EOS-Gravity Degeneracy with Masses and Pulsating Frequencies of Neutron Stars, J. Phys. G 41 (2014) 075203.  arXiv:1311.3216, doi:10.1088/0954-3899/41/7/075203.
  * [415] D.-R. Zhang, W.-Z. Jiang, S.-N. Wei, R.-Y. Yang, Q.-F. Xiang, Equation of state and hybrid star properties with the weakly interacting light U-boson in relativistic models, Eur. Phys. J. A 52 (5) (2016) 142.  arXiv:1607.01610, doi:10.1140/epja/i2016-16142-x.
  * [416] M. I. Haftel, Two-nucleon potentials and nuclear saturation, Ph.D. thesis, University of Pittsburgh (1969).
  * [417] M. I. Haftel, F. Tabakin, Nuclear Saturation and the Smoothness of Nucleon-Nucleon Potentials, Nucl. Phys. A 158 (1970) 1–42.  doi:10.1016/0375-9474(70)90047-3.
  * [418] M. Hjorth-Jensen, T. Kuo, E. Osnes, Realistic effective interactions for nuclear systems, Phys. Rept. 261 (1995) 125–270.  doi:10.1016/0370-1573(95)00012-6.
  * [419] M. Alford, M. Braby, M. W. Paris, S. Reddy, Hybrid stars that masquerade as neutron stars, Astrophys. J. 629 (2005) 969–978.  arXiv:nucl-th/0411016, doi:10.1086/430902.
  * [420] W. Wei, M. Salinas, T. Klähn, P. Jaikumar, M. Barry, Lifting the Veil on Quark Matter in Compact Stars with Core $g$-mode Oscillations, Astrophys. J. 904 (2) (2020) 187.  arXiv:1811.11377, doi:10.3847/1538-4357/abbe02.
  * [421] P. Jaikumar, A. Semposki, M. Prakash, C. Constantinou, $g$-mode oscillations in hybrid stars: A tale of two sounds, Phys. Rev. D 103 (12) (2021) 123009.  arXiv:2101.06349, doi:10.1103/PhysRevD.103.123009.
  * [422] C. Constantinou, S. Han, P. Jaikumar, M. Prakash, g modes of neutron stars with hadron-to-quark crossover transitions, Phys. Rev. D 104 (12) (2021) 123032\.  arXiv:2109.14091, doi:10.1103/PhysRevD.104.123032.
  * [423] A. Del Popolo, M. Deliyergiyev, M. Le Delliou, Solution to the hyperon puzzle using dark matter, Phys. Dark Univ. 30 (2020) 100622.  arXiv:2011.00984, doi:10.1016/j.dark.2020.100622.
  * [424] G. Raffelt, L. Stodolsky, Mixing of the Photon with Low Mass Particles, Phys. Rev. D 37 (1988) 1237.  doi:10.1103/PhysRevD.37.1237.
  * [425] G. G. Raffelt, Astrophysical methods to constrain axions and other novel particle phenomena, Phys. Rept. 198 (1990) 1–113.  doi:10.1016/0370-1573(90)90054-6.
  * [426] J.-F. Fortin, K. Sinha, Photon-dark photon conversions in extreme background electromagnetic fields, JCAP 11 (2019) 020.  arXiv:1904.08968, doi:10.1088/1475-7516/2019/11/020.
  * [427] B.-Q. Lu, C.-W. Chiang, Probing dark gauge boson with observations from neutron stars (7 2021).  arXiv:2107.07692.
  * [428] S. Knapen, T. Lin, K. M. Zurek, Light Dark Matter: Models and Constraints, Phys. Rev. D 96 (11) (2017) 115021.  arXiv:1709.07882, doi:10.1103/PhysRevD.96.115021.
  * [429] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, J. March-Russell, String Axiverse, Phys. Rev. D 81 (2010) 123530.  arXiv:0905.4720, doi:10.1103/PhysRevD.81.123530.
  * [430] R. Brito, V. Cardoso, P. Pani, Superradiance: New Frontiers in Black Hole Physics, Lect. Notes Phys. 906 (2015) pp.1–237.  arXiv:1501.06570, doi:10.1007/978-3-319-19000-6.
  * [431] M. Baryakhtar, R. Lasenby, M. Teo, Black Hole Superradiance Signatures of Ultralight Vectors, Phys. Rev. D 96 (3) (2017) 035019.  arXiv:1704.05081, doi:10.1103/PhysRevD.96.035019.
  * [432] V. Cardoso, O. J. C. Dias, G. S. Hartnett, M. Middleton, P. Pani, J. E. Santos, Constraining the mass of dark photons and axion-like particles through black-hole superradiance, JCAP 03 (2018) 043.  arXiv:1801.01420, doi:10.1088/1475-7516/2018/03/043.
  * [433] M. J. Stott, Ultralight Bosonic Field Mass Bounds from Astrophysical Black Hole Spin (9 2020).  arXiv:2009.07206.
  * [434] D. Ghosh, D. Sachdeva, Constraining Light Dark Photons from GW190517 and GW190426_152155, Phys. Rev. D 103 (9) (2021) 095028.  arXiv:2102.08857, doi:10.1103/PhysRevD.103.095028.
  * [435] N. Arkani-Hamed, L. Motl, A. Nicolis, C. Vafa, The String landscape, black holes and gravity as the weakest force, JHEP 06 (2007) 060.  arXiv:hep-th/0601001, doi:10.1088/1126-6708/2007/06/060.
  * [436] S. Alexander, E. McDonough, R. Sims, N. Yunes, Hidden-Sector Modifications to Gravitational Waves From Binary Inspirals, Class. Quant. Grav. 35 (23) (2018) 235012.  arXiv:1808.05286, doi:10.1088/1361-6382/aaeb5c.
  * [437] L. Sagunski, J. Zhang, M. C. Johnson, L. Lehner, M. Sakellariadou, S. L. Liebling, C. Palenzuela, D. Neilsen, Neutron star mergers as a probe of modifications of general relativity with finite-range scalar forces, Phys. Rev. D 97 (6) (2018) 064016.  arXiv:1709.06634, doi:10.1103/PhysRevD.97.064016.
  * [438] D. Croon, A. E. Nelson, C. Sun, D. G. E. Walker, Z.-Z. Xianyu, Hidden-Sector Spectroscopy with Gravitational Waves from Binary Neutron Stars, Astrophys. J. Lett. 858 (1) (2018) L2.  arXiv:1711.02096, doi:10.3847/2041-8213/aabe76.
  * [439] H. G. Choi, S. Jung, New probe of dark matter-induced fifth force with neutron star inspirals, Phys. Rev. D 99 (1) (2019) 015013.  arXiv:1810.01421, doi:10.1103/PhysRevD.99.015013.
  * [440] D. Krause, H. T. Kloor, E. Fischbach, Multipole radiation from massive fields: Application to binary pulsar systems, Phys. Rev. D 49 (1994) 6892–6906.  doi:10.1103/PhysRevD.49.6892.
  * [441] T. Kumar Poddar, S. Mohanty, S. Jana, Vector gauge boson radiation from compact binary systems in a gauged $L_{\mu}-L_{\tau}$ scenario, Phys. Rev. D 100 (12) (2019) 123023.  arXiv:1908.09732, doi:10.1103/PhysRevD.100.123023.
  * [442] J. A. Dror, R. Laha, T. Opferkuch, Probing muonic forces with neutron star binaries, Phys. Rev. D 102 (2) (2020) 023005.  arXiv:1909.12845, doi:10.1103/PhysRevD.102.023005.
  * [443] W. Hu, R. Barkana, A. Gruzinov, Cold and fuzzy dark matter, Phys. Rev. Lett. 85 (2000) 1158–1161.  arXiv:astro-ph/0003365, doi:10.1103/PhysRevLett.85.1158.
  * [444] A. E. Nelson, J. Scholtz, Dark Light, Dark Matter and the Misalignment Mechanism, Phys. Rev. D 84 (2011) 103501.  arXiv:1105.2812, doi:10.1103/PhysRevD.84.103501.
  * [445] P. W. Graham, J. Mardon, S. Rajendran, Vector Dark Matter from Inflationary Fluctuations, Phys. Rev. D 93 (10) (2016) 103520.  arXiv:1504.02102, doi:10.1103/PhysRevD.93.103520.
  * [446] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, F. Takahashi, Relic Abundance of Dark Photon Dark Matter, Phys. Lett. B 801 (2020) 135136.  arXiv:1810.07188, doi:10.1016/j.physletb.2019.135136.
  * [447] L. Hui, J. P. Ostriker, S. Tremaine, E. Witten, Ultralight scalars as cosmological dark matter, Phys. Rev. D 95 (4) (2017) 043541.  arXiv:1610.08297, doi:10.1103/PhysRevD.95.043541.
  * [448] K. Schutz, Subhalo mass function and ultralight bosonic dark matter, Phys. Rev. D 101 (12) (2020) 123026.  arXiv:2001.05503, doi:10.1103/PhysRevD.101.123026.
  * [449] M. Benito, J. C. Criado, G. Hütsi, M. Raidal, H. Veermäe, Implications of Milky Way substructures for the nature of dark matter, Phys. Rev. D 101 (10) (2020) 103023.  arXiv:2001.11013, doi:10.1103/PhysRevD.101.103023.
  * [450] S. Gardner, S. D. McDermott, B. Yanny, The Milky Way, coming into focus: Precision astrometry probes its evolution and its dark matter, Prog. Part. Nucl. Phys. 121 (2021) 103904.  arXiv:2106.13284, doi:10.1016/j.ppnp.2021.103904.
  * [451] L. Hui, Wave Dark Matter (1 2021).  arXiv:2101.11735, doi:10.1146/annurev-astro-120920-010024.
  * [452] A. Pierce, K. Riles, Y. Zhao, Searching for Dark Photon Dark Matter with Gravitational Wave Detectors, Phys. Rev. Lett. 121 (6) (2018) 061102.  arXiv:1801.10161, doi:10.1103/PhysRevLett.121.061102.
  * [453] H.-K. Guo, K. Riles, F.-W. Yang, Y. Zhao, Searching for Dark Photon Dark Matter in LIGO O1 Data, Commun. Phys. 2 (2019) 155.  arXiv:1905.04316, doi:10.1038/s42005-019-0255-0.
  * [454] A. L. Miller, et al., Probing new light gauge bosons with gravitational-wave interferometers using an adapted semicoherent method, Phys. Rev. D 103 (10) (2021) 103002.  arXiv:2010.01925, doi:10.1103/PhysRevD.103.103002.
  * [455] Y. Michimura, T. Fujita, J. Kume, S. Morisaki, K. Nagano, H. Nakatsuka, A. Nishizawa, I. Obata, Ultralight dark matter searches with KAGRA gravitational wave telescope, in: 17th International Conference on Topics in Astroparticle and Underground Physics, 2021.  arXiv:2111.00420.
  * [456] C.-R. Chen, C. S. Nugroho, Detection Prospects of Dark Matter in Einstein Telescope (11 2021).  arXiv:2111.11014.
  * [457] K. Nomura, A. Ito, J. Soda, Pulsar timing residual induced by ultralight vector dark matter, Eur. Phys. J. C 80 (5) (2020) 419.  arXiv:1912.10210, doi:10.1140/epjc/s10052-020-7990-y.
  * [458] X. Xue, et al., High-precision search for dark photon dark matter with the Parkes Pulsar Timing Array (12 2021).  arXiv:2112.07687.
  * [459] G. Bertone, et al., Gravitational wave probes of dark matter: challenges and opportunities, SciPost Phys. Core 3 (2020) 007.  arXiv:1907.10610, doi:10.21468/SciPostPhysCore.3.2.007.
  * [460] S. Ascenzi, G. Oganesyan, M. Branchesi, R. Ciolfi, Electromagnetic counterparts of compact binary mergers, J. Plasma Phys. 87 (1) (2021) 845870102\.  arXiv:2011.04001, doi:10.1017/S0022377820001646.
  * [461] M. Diamond, G. Marques-Tavares, Gamma-rays flashes from dark photons in neutron star mergers (6 2021).  arXiv:2106.03879.
  * [462] B. D. Metzger, Kilonovae, Living Rev. Rel. 23 (1) (2020) 1.  arXiv:1910.01617, doi:10.1007/s41114-019-0024-0.
  * [463] I. Goldman, S. Nussinov, Weakly Interacting Massive Particles and Neutron Stars, Phys. Rev. D 40 (1989) 3221–3230.  doi:10.1103/PhysRevD.40.3221.
  * [464] C. Kouvaris, WIMP Annihilation and Cooling of Neutron Stars, Phys. Rev. D 77 (2008) 023006.  arXiv:0708.2362, doi:10.1103/PhysRevD.77.023006.
  * [465] C. Kouvaris, P. Tinyakov, Excluding Light Asymmetric Bosonic Dark Matter, Phys. Rev. Lett. 107 (2011) 091301.  arXiv:1104.0382, doi:10.1103/PhysRevLett.107.091301.
  * [466] N. F. Bell, A. Melatos, K. Petraki, Realistic neutron star constraints on bosonic asymmetric dark matter, Phys. Rev. D 87 (12) (2013) 123507.  arXiv:1301.6811, doi:10.1103/PhysRevD.87.123507.
  * [467] C. Kouvaris, N. G. Nielsen, Asymmetric Dark Matter Stars, Phys. Rev. D 92 (6) (2015) 063526.  arXiv:1507.00959, doi:10.1103/PhysRevD.92.063526.
  * [468] J. Ellis, A. Hektor, G. Hütsi, K. Kannike, L. Marzola, M. Raidal, V. Vaskonen, Search for Dark Matter Effects on Gravitational Signals from Neutron Star Mergers, Phys. Lett. B 781 (2018) 607–610.  arXiv:1710.05540, doi:10.1016/j.physletb.2018.04.048.
  * [469] A. Nelson, S. Reddy, D. Zhou, Dark halos around neutron stars and gravitational waves, JCAP 07 (2019) 012.  arXiv:1803.03266, doi:10.1088/1475-7516/2019/07/012.
  * [470] D. Radice, A. Perego, F. Zappa, S. Bernuzzi, GW170817: Joint Constraint on the Neutron Star Equation of State from Multimessenger Observations, Astrophys. J. Lett. 852 (2) (2018) L29.  arXiv:1711.03647, doi:10.3847/2041-8213/aaa402.
  * [471] S. Reddy, D. Zhou, Dark Lepton Superfluid in Proto-Neutron Stars (7 2021).  arXiv:2107.06279.
  * [472] J. F. Beacom, M. R. Vagins, GADZOOKS! Anti-neutrino spectroscopy with large water Cherenkov detectors, Phys. Rev. Lett. 93 (2004) 171101.  arXiv:hep-ph/0309300, doi:10.1103/PhysRevLett.93.171101.
  * [473] S. Horiuchi, T. Kinugawa, T. Takiwaki, K. Takahashi, K. Kotake, Impact of binary interactions on the diffuse supernova neutrino background, Phys. Rev. D 103 (4) (2021) 043003.  arXiv:2012.08524, doi:10.1103/PhysRevD.103.043003.
  * [474] I. Andreoni, et al., Target of Opportunity Observations of Gravitational Wave Events with Vera C. Rubin Observatory (11 2021).  arXiv:2111.01945.


A comparative study of $2d$ Ising model at different boundary conditions using
non-deterministic Hexagonal Cellular Automata
†Jahangir Mohammed111 Present address: Department of Physics, Nabarangpur
College, Nabarangpur-764063, Odisha, India. and ∗,‡Swapna Mahapatra

Department of Physics, Utkal University,
Bhubaneswar, Odisha 751004, India
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>

Abstract

> The spin system of the $2d$ Ising model having a hexagonal-lattice is
> simulated using non-deterministic Cellular Automata. The method to implement
> this program is outlined and our results show a good approximation to the
> exact analytic solution. The phase transition in $2d$ Ising model is studied
> with a $40\times 40$ hexagonal-lattice with five different boundary
> conditions (bcs) i.e., adiabatic, periodic, reflexive, fixed $+1$ and fixed
> $-1$ with random orientation of spins as initial conditions in the absence
> of an external applied magnetic field. The critical temperature below which
> the spontaneous magnetization appears as well as other physical quantities
> such as the magnetisation, energy, specific heat, susceptibility and entropy
> with each of the bcs are calculated. The phase transition occurs around
> $T^{H}_{c}$ = 1.5 which approximates well with the result obtained from
> exact analytic solution by Wannier and Houtappel. We compare the behaviour
> of magnetisation per cell for five different types of bcs by calculating the
> number of points close to the line of zero magnetisation for $T>T^{H}_{c}$.
> We find that the periodic, adiabatic and reflexive bcs give closer
> approximation to the value of $T^{H}_{c}$ than fixed $+1$ and fixed $-1$ bcs
> with all three initial conditions for lattice size less than $50\times 50$.
> However, for lattice size between $50\times 50$ and $200\times 200$, fixed
> $+1$ bc and fixed $-1$ bc give closer approximation to the $T^{H}_{c}$ with
> initial conditions in which all spins are in down configuration and all
> spins are in up configuration respectively.

##  1 Introduction

The simple Hamiltonian for the Ising model and its variants are studied
extensively for their ability to describe critical phenomena. This can be
applied to understand magnetism, models for high-temperature superconductivity
and phase diagrams, disordered and non-equilibrium systems. The model also
provides means for testing algorithms. All the interesting phenomena of
statistical mechanics and phase transitions are found to have counterparts in
this model. It is an example of a lattice model for modeling simple
interacting many-particle systems in statistical mechanics where a set of
spins $s_{i}=\pm 1$ is assigned to each lattice site.

The spins interact with both their nearest neighbors and an external magnetic
field via a Hamiltonian of the form

$H(s)=-\sum\limits_{<ij>}J_{ij}s_{i}s_{j}-\mu\sum\limits_{k=1}^{N}h_{k}s_{k}$
(1.1)

where $s$ is an arbitrary spin configuration, the notation $<ij>$ indicates a
sum over nearest neighbor lattice points, $J_{ij}$($>0$) is a parameter
describing inter-particle interactions, $\mu$ is the magnetic moment and
$h_{k}$ is an external magnetic field at $k^{th}$ spin.

The model is best known for describing the emergence of ferromagnetism in
crystals of atoms that interact via spin-spin coupling. The Ising model is
useful even in cases for which it cannot be solved exactly. On both square-
lattice and hexagonal-lattice geometries where the Ising model can be exactly
solved and also for those where no exact solution has been found, numerical
simulations can provide insight into critical behavior of these models. The
most common method for simulating Ising systems is the Metropolis-Hastings
algorithm, originally developed for use in molecular dynamics simulations.

The $1d$ model introduced by Lenz was exactly solved by Ising, which showed
that the one-dimensional case does not exhibit a phase transition [1]. This
led to the general belief that the Ising model in higher dimensions would also
be ineffective at describing systems with critical points. However, Onsager in
1944 exactly solved the Ising model for a two-dimensional square lattice
configuration with periodic boundary condition by analytical method [2]. He
demonstrated that, for an infinite square lattice, the ferromagnetic phase
transition occurs at a critical temperature, $T^{S}_{c}$, of

${T^{S}_{c}}=\frac{2}{\ln{(1+\sqrt{2})}}\approx 2.269$ (1.2)

where $J=1$, $k_{B}=1$ and $S$ stands for square-lattice. The square and
hexagonal lattices, both, have exactly computable critical temperatures. The
hexagonal (honeycomb) lattice is the next simplest two dimensional lattice.
The pure Ising model on the honeycomb lattice has been studied by Wannier [3]
and Houtappel [4] and have been exactly solved. The most straightforward way
to find these critical temperatures is via duality between high-temperature
and low-temperature behaviour. Duality is the hidden symmetry found by Kramers
and Wannier that relates the partition function and the free energy of the
Ising model at low and high temperatures for the two dimensional square-
lattice [5, 6]. The system is mathematically modeled in two ways each of which
is independent of the other and both of which are valid to describe the
physical properties of the system. If one assumes that the free energy is
singular at the critical temperature, and that this singularity is unique,
then this leads to the determination of the critical temperature. For the two
dimensional hexagon system the critical temperature ($T^{H}_{c}$) was found to
be

${T^{H}_{c}}=\frac{2}{\ln{(2+\sqrt{3})}}\approx 1.519$ (1.3)

where $J=1$, $k_{B}=1$ and $H$ stands for hexagonal-lattice. The above results
are for lattice of infinite size with periodic boundary condition. Above this
temperature the average magnetization is zero in the absence of external
magnetic field. We performed the simulation taking lattices of finite size and
in the absence of external magnetic field so as to compare the result obtained
by the above analytic method as well as to those for our five different bcs.

A Cellular Automaton (CA) is a mathematical model, modeling a set of cells
which interact with their neighbors. In this model, each cell have values
known as states, all the cells update their states simultaneously at discrete
time steps, and the new state of a cell is determined by current state of its
neighbors according to a local function called rule of the CA. Hexagonal CA is
a tessellation of the plane by regular hexagons which provide for higher
packing density of cells and the unit cells of hexagonal grids are uniformly
connected in the sense that the distance from a given cell to any adjacent
cell is the same. To simulate the Ising model with Hexagonal CA, we can design
a two states CA, for spin up state ($+1$) and spin down state ($-1$) and
number of neighbors is six i.e, north, south, east, west, top and bottom. A
definite rule is designed for which states of the cells are either all in up
states or all in down states below $T^{H}_{c}$ and above $T^{H}_{c}$ on the
average half in $+1$ spin states and half in $-1$ spin states. The $2d$ Ising
model using non-deterministic hexagonal CA (HCA) has not been studied for five
different boundary conditions. The CA used in this work differs from that of
Metropolis-Hastings algorithm where the algorithm specifies that the
transitions must be made for one site at a time whereas we consider
transitions of many sites simultaneously. The work of Eltinge [7] uses
Metropolis-Hastings algorithm to study numerically the case of periodic
boundary condition.

Our computations will be on finite-size lattices with four different boundary
conditions, as well as, with the periodic boundary condition that effectively
simulates accurately infinite-size lattice results in many respects. Apart
from the Ising model, quantum spin models such as the Kitaev model [8] or
Kitaev-Heisenberg model [9, 10] on the honeycomb lattice have recently
received a lot of theoretical and experimental attention [11, 12, 13, 14, 15].
Hexagonal lattice has also been discussed in the context of high temperature
and low temperature susceptibility series and extracting the scaling function
[16]. Such a growing interest in the hexagonal-lattice systems motivates us to
revisit the Ising model on the hexagonal lattice. Attempts have been made for
mapping Ising models in different lattice geometry using CA. A deterministic
CA (DCA) proposed by Domany and Kinzel [17], the Q2R CA [18, 19, 20, 21] and
the Creutz CA [22, 23, 24, 25] are mostly used in analysing square-lattice
Ising model. All these CA models are deterministic and the computation can be
performed fast. It has been demonstrated that the probabilistic model of the
CA like Metropolis algorithm [26] is more realistic for description of the
Ising model even though the random number generation makes it slower.
Probabilistic CA model under five different boundary conditions has been
studied in the context of a square-lattice Ising model [27]. However, Ising
model using two dimensional hexagonal CA (HCA) under different boundary
conditions other than periodic boundary condition has not yet been studied.

The paper is organised as follows: in _section 2_ , we discuss the basic
theory to treat a $2d$ HCA and how to implement it in the Ising model with
five different bcs. The results obtained out of the simulations are given and
are analysed in _section 3_ and we have also compared Ising model with
different boundary conditions by considering their converging points. Our
conclusion and future perspective are discussed in _section 4_.

##  2 Implementation of Isotropic $2d$ Ising Model by Hexagonal CA

Two dimensional CA is described by finite states of cells ($s$), neighborhood
cells ($n$) and its distance among neighbourhood ($r$), boundary conditions
and transition functions or rules ($f$). In our $2d$ HCA model,
$s=\\{s_{i,j},\,s_{i,j}\in-1/+1\\}$, number of nearest neighbour cells $n=6$,
$r=1$ and we consider five different bcs.

Figure 1: Hexagonal-lattice with row and column vectors.

Figure 2: $4\times 4$ hexagonal-lattice.

Neighbourhoods of extreme cells are taken care of by boundary condition (bc).
If the extreme cells are adjacent to each other then it is called periodic bc
(pbc). In adiabatic bc (abc), the extreme cells replicate their state and in
reflexive bc (rbc), mirror position states replace the extreme cells. In fixed
bc, the extreme cells are connected to $+1$ or $-1$ state. If it is connected
to $+1$ state, it is called fixed $+1$ bc ($f+1bc$) and if it is connected to
$-1$ state, then it is called fixed -1 bc ($f-1bc$).

If the same rule is applied to all the elements of the matrix ($s$), then it
is called _uniform CA_ and if different rules are applied to individual
elements of the matrix or block of elements then it is called _nonuniform CA_.
At different time intervals, if different rules are applied to the matrix then
it is called varying CA e.g., _probabilistic CA_. With the application of
these rules, elements (states) of the matrix change at successive intervals as
shown in the following equation.

$\displaystyle s_{L\times L}^{t+1}=f_{L\times L}^{t}\times s_{L\times L}^{t}$
(2.1)

where $f$ is a time varying rule or transition matrix.

Consider an isotropic $2d$ Ising model in the form of hexagonal lattice ($s$)
with $L$ rows and $L$ columns. Here we consider each hexagonal cell as a
lattice point as shown in figure 2. The lattice has then $L^{2}=N$ sites. The
odd rows have different column positions i.e., $1,3,5,\dots$ etc. and even
rows have different column positions i.e., $2,4,6,\dots$ etc. At each of the
sites $s_{i,j},\,i=1,\,\dots,\,L;j=1,\,\dots,\,2L$, the spins are arranged in
such a way that $j$ increases as $1,3,\,\dots,\,2L-1$ from left to right in
case of odd rows, $j$ increases as $2,4,\,\dots,\,2L$ from left to right in
case of even rows and $i$ increases from bottom to top, and has one of the
$\pm 1$ spin, which are the two states in CA. So, there are $2^{L^{2}}$ spin
configurations. We consider only the nearest neighbor interactions, so the
number of neighbor cells are $6$. We include the five different bcs as
follows.

_For odd rows :_

  1. 1.

pbc : $s_{i,2L+1}=s_{i,1}$, $s_{L+1,j+1}=s_{1,j}$,
$s_{i,-1}=s_{i,2L-1}$ and $s_{-1,j}=s_{L,j+1}$.

  2. 2.

abc : $s_{i,2L+1}=s_{i,2L-1}$, $s_{L+1,j+1}=s_{L,j}$,
$s_{i,-1}=s_{i,1}$ and $s_{-1,j}=s_{1,j-1}$.

  3. 3.

rbc : $s_{i,2L+1}=s_{i,2L-3}$, $s_{L+1,j+1}=s_{L-1,j+1}$,
$s_{i,-1}=s_{i,3}$ and $s_{-1,j}=s_{2,j}$.

  4. 4.

f+1bc : $s_{i,2L+1}=+1$, $s_{L+1,j+1}=+1$,
$s_{i,-1}=+1$ and $s_{-1,j}=+1$.

  5. 5.

f-1bc : $s_{i,2L+1}=-1$, $s_{L+1,j+1}=-1$,
$s_{i,-1}=-1$ and $s_{-1,j}=-1$.

_For even rows :_

  1. 1.

pbc : $s_{i,2L+2}=s_{i,2}$, $s_{L+1,j}=s_{1,j-1}$,
$s_{i,0}=s_{i,2L}$ and $s_{0,j}=s_{L,j+1}$.

  2. 2.

abc : $s_{i,2L+2}=s_{i,2L}$, $s_{L+1,j}=s_{L,j+1}$,
$s_{i,0}=s_{i,2}$ and $s_{0,j}=s_{1,j-1}$.

  3. 3.

rbc : $s_{i,2L+2}=s_{i,2L-2}$, $s_{L+1,j}=s_{L-1,j}$,
$s_{i,0}=s_{i,4}$ and $s_{0,j}=s_{2,j}$.

  4. 4.

f+1bc : $s_{i,2L+2}=+1$, $s_{L+1,j}=+1$,
$s_{i,0}=+1$ and $s_{0,j}=+1$.

  5. 5.

f-1bc : $s_{i,2L+2}=-1$, $s_{L+1,j}=-1$,
$s_{i,0}=-1$ and $s_{0,j}=-1$.

The average magnetization for the configuration is defined as,

$\left\langle
M\right\rangle=\sum\limits_{i=1}^{L}{\sum\limits_{\genfrac{}{}{0.0pt}{}{j=1}{i\%2=j\%2}}^{2L}}s_{i,j}$
(2.2)

and the average magnetization per spin is given by

$\left\langle m\right\rangle=\frac{\left\langle M\right\rangle}{N}$ (2.3)

Energy for the configuration $s$ is defined as

$E(s)=-\frac{J}{6}\sum\limits_{i=1}^{L}{\sum\limits_{\genfrac{}{}{0.0pt}{}{j=1}{i\%2=j\%2}}^{2L}}s_{i,j}\times(s_{i,j-2}+s_{i,j+2}+s_{i-1,j-1}+s_{i-1,j+1}+s_{i+1,j-1}+s_{i+1,j+1})$
(2.4)

Here, the factor of $1/6$ has been put to remove the sextuple counting of
energy, otherwise the interacting energy will be computed six times.
$J_{ij}=J$ (isotropic) for $6$ neighbour cells, or else, $J_{ij}=0$.

The configuration energy per spin is

$\left\langle e\right\rangle=\frac{E(s)}{N}$ (2.5)

For updating the lattice in next iteration, we use the probabilistic approach
by constructing a probabilistic CA. We use the following procedure.

We calculate the change in energy i.e., the energy difference at successive
time intervals is $\Delta E(s^{t})=E(s^{t})-E(s^{t-1})$. Here, $\Delta E\geq
0$ or $\Delta E<0$ but we have considered the case $\Delta E\geq 0$ i.e.,
$E(s^{t})\geq E(s^{t-1})$. Next we calculate the probability of each site in
the spin configuration $s$ at time $t$ (number of iterations) by using the
Boltzmann factor

$p_{t}=\frac{p(E(s^{t}))}{p(E(s^{t-1}))}=e^{-\frac{\Delta E(s^{t})}{k_{B}T}}$
(2.6)

With the above probability for each site, we construct a probability weighted
matrix (or transition matrix). This matrix leads to our probabilistic CA
matrix ($PCA^{t}$) by comparing it with a random matrix and multiplying by a
factor of $0.1$ to normalise the $PCA^{t}$.

Successive spin configurations are obtained from

$[s_{i,j}^{t+1}]_{L\times L}=[PCA_{i,j}^{t}]_{L\times L}[s_{i,j}^{t}]_{L\times
L}.$ (2.7)

Here, we consider $s_{i,j}^{t+1}$ to be a new configuration. If in this new
configuration the calculated energy is less than the initial configuration,
then it is allowed to proceed forward; otherwise the new configuration is
flipped backward. After a finite number of iterations we calculate the average
energy of the system per cell ($e$), magnetisation per cell ($m$), and obtain
the susceptibility per cell ($\chi$), specific heat per cell ($C_{v}$) and the
entropy per cell ($S$).

$\chi=\frac{N}{k_{B}T}(\left\langle m^{2}\right\rangle-\left\langle
m\right\rangle^{2})$ (2.8)

$C_{v}=\frac{N}{k_{B}T^{2}}(\left\langle E^{2}\right\rangle-\left\langle
E\right\rangle^{2})$ (2.9)

$S=-k_{B}(r_{1}P_{1}\log_{2}P_{1}-r_{2}P_{2}\log_{2}P_{2})$ (2.10)

where $r_{1}$ is the total number of spin up states, $r_{2}$ is the total
number of spin down states, $P_{1}$ is the probability of spin up states and
$P_{2}$ is the probability of spin down states in the lattice $s$. The
Metropolis algorithm specifies that transitions must be made for one site at a
time. But here we have considered transitions of many sites simultaneously for
which we go to a regime of oscillating behavior. For instance, nearly every
site that is flipped in the direction of higher energy becomes unflipped on
the next iteration. CA is desirable to perform transitions on every lattice
site simultaneously. For instance such simultaneous transitions take advantage
of matrix operations for increased computational efficiency. This algorithm
checks the time complexity better than the Metropolis algorithm [26].

##  3 Simulation, results and discussions

In this work, we have considered hexagonal-lattice of size $40\times 40$ for
the determination of critical temperature at which a phase transition occurs.
Here we have considered lattice size of $40\times 40$ for our convenience in
computation. In this simulation we take $J=1$, $k_{B}=1$ and the external
magnetic field to be zero i.e., $h=0$. We have considered the temperature
range from $0.1$ to $4.0$ with increment of $0.05$ unit in each step. At each
temperature, we have considered $10^{5}$ runs. Each of the simulation starts
with random spin configurations. We have considered five different bcs i.e.,
periodic, adiabatic, reflexive, fixed $+1$ and fixed $-1$.

The magnetisation per cell ($m$) versus temperature results are shown for all
five bcs in figure 3 . One can observe that in the case of fixed $\pm 1$ bcs,
the curves obtained provide better result than either of the other boundary
conditions, i.e., pbc, abc or rbc at low temperature regions ($T<T^{H}_{c}$).

Figure 3: $m$ versus $T$ for all five bcs with random initial spin
configuration in $25\times 25$ hexagonal-lattice.

###  3.1 Phase transition with pbc, abc, rbc, f1bc and f-1bc

To study the phase transition we have calculated energy per cell ($e$),
magnetisation per cell ($m$), susceptibility per cell ($\chi$), specific heat
per cell ($C_{v}$) and entropy of the configuration ($S$) by considering
average of ten simulation. Here each simulation starts with random spin
configuration in hexagonal-lattice of size $40\times 40$ and we take
temperature range from $1.0$ to $4.0$ with increment of $0.05$ unit since we
have to find the phase transition corresponding to all five different bcs and
hence the $T^{H}_{c}$.

In figure 4, we have plotted $e$ vs $T$; $m$ vs $T$; $m$ vs $e$; $\chi$ vs
$T$; $C_{v}$ vs $T$ and $S$ vs $T$ with initial random spin configuration with
pbc, abc, rbc, f+1bc and f-1bc respectively. One finds that around temperature
$T^{H}_{c}=1.5$, the magnetisation approaches the value zero, the energy
gradually increases, as well as the entropy gradually increases, as shown in
figures 4(a), 4(b) and 4(f). The susceptibility and specific heat also change,
initially they increase with a sudden rise at $T=1.5$ and then start
decreasing as shown in figures 4(d) and 4(e) respectively. So, a phase
transition is clearly visible around $T^{H}_{c}=1.5$ as the second order
derivative of energy and magnetisation is extremely large with all five bcs.

In the $m$ vs $e$ graphs shown in figure 4(c), the regions with higher density
of points indicate three states. We find three states, two of them are low
temperature ground states around ($m=\pm 1$, $e=-6$) with random initial
conditions in the cases of abc, pbc and rbc, and the third one is the high
temperature phase which is centered at ($m=0$, $e=0$). But in the case of
f+1bc there is one low temperature ground state around ($m=1$, $e=-6$) and one
high temperature phase is centered at ($m=0$, $e=0$). Also in the case of
f-1bc there is one low temperature ground state around ($m=-1$, $e=-6$) and
one high temperature phase is centered at ($m=0$, $e=0$).

|
---|---
(a) $m$ versus $T$ | (b) $e$ versus $T$
|
(c) $m$ versus $e$ | (d) $\chi$ versus $T$
|
(e) $C_{v}$ versus $T$ | (f) $S$ versus $T$
Figure 4: Random initial spin configuration in $40\times 40$ hexagonal-lattice
with five different bcs.

##  4 Comparison of $2d$ Ising model simulation among boundary conditions

In figure 3, one can observe that the magnetization approach the zero line
after $T>T^{H}_{c}$ in different manner depending on the bcs. The analytic
solution gives a zero magnetization above $T^{H}_{c}$ and our simulation for
different bcs also show that it converge toward the zero value at increased
temperatures. With one simulation for all bcs, it is not possible to predict
which bc is closer to $T^{H}_{c}$. So, we analyze the points for magnetization
in the range $-0.01\leq m\leq+0.01$, $-0.1\leq m\leq+0.1$ and $-0.2\leq
m\leq+0.2$ which are close to the zero line of magnetization after
$T>T^{H}_{c}$. We call such points as converging points.

In this simulation, we have taken various lattice sizes ranging from $5\times
5$ to $200\times 200$ in this manner. We consider an increment of lattice size
$5$ in each step while going from lattice size from $5\times 5$ to $60\times
60$ and subsequently for $60\times 60$ to $100\times 100$ an increment of $10$
for each step and from $100\times 100$ to $200\times 200$ with increment of
$50$ for each step. The temperature ranges from $1.0$ to $4.0$ with small
increments of $0.05$ units. We have iterated $1000$ times for each simulation
for lattice sizes $\leq 10\times 10$, $10000$ times for lattice sizes in
between $10\times 10$ and $40\times 40$ and $100000$ times for lattice sizes
above $40\times 40$. Figure 5 shows the converging points in the above
mentioned range of magnetization. We have counted the number of converging
points as defined above by taking average result of ten simulations for each
bc with three different initial conditions. We carry out all the simulation
with three initial conditions and with all five bcs. Here, three initial
conditions are (i) all up (or most of spins up), (ii) all down (or most of the
spins down) and (iii) random (or randomly oriented spins up/down).

Figure 5: Converging points of bcs for $50\times 50$ lattice size after
$T>T^{H}_{c}$. Red lines are for magnetization $|m|\leq 0.2$, blue lines are
for magnetization $|m|\leq 0.1$ and yellow lines are for magnetization
$|m|\leq 0.01$ with the initial condition of random spin configuration.
Magenta line represents a parallel line to magnetization per cell at
$T=T^{H}_{c}$.

For lattice $sizes\leq 50\times 50$ with pbc, abc and rbc, one finds more
convergent points in both cases $-0.1\leq m\leq+0.1$ and $-0.2\leq m\leq+0.2$
than that for f+1bc and f-1bc for all up and all down spin initial conditions.

With lattice sizes in between $50\times 50$ and $200\times 200$ with f-1bc,
one finds more convergent points in both the cases corresponding to $-0.1\leq
m\leq+0.1$ and $-0.2\leq m\leq+0.2$ with initial condition of all up spins.
For f+1bc, one finds more convergent points in both the cases with initial
condition of all spins down.

In the case of lattice $sizes\leq 30\times 30$ with pbc and lattice
$size>30\times 30$ with abc, one finds more convergent points both cases
$-0.1\leq m\leq+0.1$ and $-0.2\leq m\leq 0.2$ among all bcs with random spin
initial conditions. For all initial cases with $-0.01\leq m\leq+0.01$, the
result is better for pbc, abc and rbc and converging points increase as
lattice sizes are increased.

It is observed that lattice $size=200\times 200$ with all five bcs, one can
observes same convergent points in both cases $-0.1\leq m\leq+0.1$ and
$-0.2\leq m\leq+0.2$ but pbc, abc and rbc show more convergent points in
$-0.01\leq m\leq+0.01$ than f+1bc and f-1bc.

##  5 Conclusion

In the present work, we have studied $2d$ Ising model in hexagonal-lattice
with five different boundary conditions using non-deterministic Cellular
Automata algorithm in which transitions of many sites takes place
simultaneously in contrast to transitions being made for one site at a time.
We have observed a phase transition that occurs at the critical temperature
$T^{H}_{c}\approx 1.5$ for each of the five bcs. The result of
$T^{H}_{c}\approx 1.5$ agrees with the exact solution presented in equation
$1.3$. The findings for fixed $+1/-1$ bcs also show a smoother curve at low
temperatures than pbc, abc and rbc with random initial condition. This implies
that with the different initial conditions on different lattice sizes $\leq
50\times 50$, one can take care of boundary spins by not only pbc but also by
abc and rbc. Further in our analysis, abc shows more converging points than
pbc and rbc for lattice size $>30\times 30$ with random initial condition. For
lattice size greater than $50\times 50$, f+1bc and f-1bc are better suited to
use than other bcs in case of initial spin configuration with all down or all
up spins respectively. It is also observed that for lattice $size=200\times
200$, in all cases of initial spin configuration, in all bcs give almost same
converging points. From the simulation point of view, our algorithm flipped
most of the spins at a single iteration. It is computationally more efficient
than Metropolis algorithm [26] which flips single spins at a single iteration.
The analysis done in this paper will help us in finding the values of critical
exponents more accurately for hexagonal-lattice, which we plan to study next.
The present work is expected to facilitate the use of non-deterministic CA in
the study of phase transition in case of anti-ferromagnetic materials, binary
alloys and spin glasses etc with different boundary conditions.

### Acknowledgement

We thank K. Maharana and S. Pattanayak for useful discussions. JM would like
to thank Institute of Physics, Bhubaneswar for hospitality where part of this
work has been done.

## References

  * [1] E. Ising, Beitrag zur Theorie des Ferromagnetismus,Z. Phys. 31:253 (1925).
  * [2] L. Onsager, Crystal statistics. I. A two dimensional model with an order-disorder transition, Phys. Rev., 65, 117–149 (1944).
  * [3] G. H. Wannier, Antiferromagnetism. The triangular Ising net., Phys. Rev., 79, 357–364 (1950).
  * [4] R. Houtappel, Order-disorder in hexagonal lattices, Physica, 16, 425–455 (1950).
  * [5] H. A. Kramers and G. H. Wannier, Statisics of the Two-Dimensional Ferromagnet. Part I, Phys. Rev., 60, 252-262 (1941).
  * [6] H. A. Kramers and G. H. Wannier, Statisics of the Two-Dimensional Ferromagnet. Part II, Phys. Rev., 60, 263-276 (1941).
  * [7] http://web.mit.edu/8.334/www/grades/projects/projects15/EltingeStephen.pdf
  * [8] A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys., 321, 2–111 (2006).
  * [9] J. Chaloupka, G. Jackeli and G. Khaliullin, Kitaev-Heisenberg model on a honeycomb lattice: Possible exotic phases in iridium oxides A2IrO3, Phys. Rev. Lett., 105, 027204 (2010).
  * [10] G. Jackeli and G. Khaliullin, Mott insulators in the strong spin-orbit coupling limit: From Heisenberg to a quantum compass and Kitaev models, Phys. Rev. Lett., 102, 017205 (2009).
  * [11] J. Chaloupka, G. Jackeli and G. Khaliullin, Zigzag magnetic order in the iridium oxide Na2IrO3, Phys. Rev. Lett., 110, 097204 (2013).
  * [12] S.K. Choi, R. Coldea, A.N. Kolmogorov, T. Lancaster, I.I. Mazin, S.J. Blundell, P.G. Radaelli, Y. Singh, P. Gegenwart, K.R. Choi, S.W. Cheong, P.J. Baker, C. Stock and J. Taylor, Spin waves and revised crystal structure of honeycomb iridate Na2IrO3, Phys. Rev. Lett., 108, 127204 (2012).
  * [13] J. Reuther, R. Thomale and S. Trebst, Finite-temperature phase diagram of the Heisenberg-Kitaev model, Phys. Rev. B, 84, 100406 (2011).
  * [14] Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst and P. Gegenwart, Relevance of the Heisenberg-Kitaev model for the honeycomb lattice iridates A2IrO3, Phys. Rev. Lett., 108, 127203 (2012).
  * [15] F. Ye, S. Chi, H. Cao, B.C. Chakoumakos, J.A. Fernandez-Baca, R. Custelcean, T.F. Qi, O.B. Korneta and G. Cao, Direct evidence of a zigzag spin-chain structure in the honeycomb lattice: A neutron and X-ray diffraction investigation of single-crystal Na2IrO3, Phys. Rev. B, 85, 180403 (2012).
  * [16] Y. Chan, A. J. Guttmann, B. G. Nickel and J. H. H. Perk, The Ising susceptibility scaling function, arXiv:1012.5272v3, [cond-mat.stat-mech] (2011).
  * [17] E. Domany and W. Kinzel, Equivalence of Cellular Automata to Ising Models and Directed Percolation, Phys. Rev. Lett., 53, 4 (1984).
  * [18] G. Y. Vichniac, Simulating physics with cellular automata, Physica D., 10, 96-116 (1984).
  * [19] J. G. Zabolitzky and H. J. Herrmann, Multitasking case study on the Cray-2: The Q2R cellular automaton, J. Comp. Phys., 76, 426-447 (1988).
  * [20] C. Moukarzel and N. Parga, On the evaluation of magnetisation fluctuations with Q2R cellular automata, J. Phys. A, 22, 943 (1989).
  * [21] D. Stauffer, Critical 2D and 3D Dynamics of Q2R Cellular Automata, Comp. Phys. Comm., 127, 113-119 (2000).
  * [22] M. Creutz, Deterministic Ising dynamics, Ann. Phys., 167(1), 62–72 (1986).
  * [23] N. Aktekin, Simulation of the three-dimensional Ising model on the Creutz cellular automaton, Physica A, 219, 3-4, 436–446 (1995).
  * [24] N. Aktekin, Simulation of the four-dimensional Ising model on the Creutz cellular automaton, Physica A, 232, 1-2, 397–407 (1996).
  * [25] N. Aktekin, Simulation of the Eight-Dimensional Ising Model on the Creutz Cellular Automaton, Int. J. Mod. Phys. C, 8(2), 287-292 (1997).
  * [26] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, Equation of state calculations by fast computing machines, The Journal of Chemical Physics, 21(6), 1087–1092 (1953).
  * [27] J. Mohammed and S. Mahapatra, A comparative study of 2d Ising model at different boundary conditions using Cellular Automata, Int. J. Mod. Phys. C, 29(8), 1850066 (2018).