Also at ]Dokuz Eylül University, Graduate School of Natural and Applied
Sciences, Turkey

# Critical behavior and phase diagrams of a spin-1 Blume-Capel model with
random crystal field interactions: An effective field theory analysis

Yusuf Yüksel  [  Ümit Akıncı  Department of Physics, Dokuz Eylül University,
TR-35160 Izmir, Turkey  Hamza Polat<EMAIL_ADDRESS>Department of
Physics, Dokuz Eylül University, TR-35160 Izmir, Turkey

###### Abstract

A spin-1 Blume-Capel model with dilute and random crystal fields is examined
for honeycomb and square lattices by introducing an effective-field
approximation that takes into account the correlations between different spins
that emerge when expanding the identities. For dilute crystal fields, we have
given a detailed exploration of the global phase diagrams of the system in
$k_{B}T_{c}/J-D/J$ plane with the second and first order transitions, as well
as tricritical points. We have also investigated the effect of the random
crystal field distribution characterized by two crystal field parameters $D/J$
and $\triangle/J$ on the phase diagrams of the system. The system exhibits
clear distinctions in qualitative manner with coordination number $q$ for
random crystal fields with $\triangle/J,D/J\neq 0$. We have also found that,
under certain conditions, the system may exhibit a number of interesting and
unusual phenomena, such as reentrant behavior of first and second order, as
well as a double reentrance with three successive phase transitions.

PACS numbers

75.10.Dg, 75.10.Hk, 75.30.Kz

††preprint: APS/123-QED

###### Contents

  1. I Introduction
  2. II Formulation
  3. III Results and Discussion
    1. III.1 Phase diagrams of the system with dilute crystal field
    2. III.2 Phase diagrams of the system with random crystal field
  4. IV Conclusions
  5. A Fundamental correlation functions of the system for a square lattice
  6. B The complete set of twenty one linear equations of a honeycomb lattice

##  I Introduction

Spin-1 Blume-Capel (BC) model blume ; capel  is one of the most extensively
studied models in statistical mechanics and condensed matter physics. The
model exhibits a variety of multicritical phenomena such as a phase diagram
with ordered ferromagnetic and disordered paramagnetic phases separated by a
transition line that changes from a continuous phase transition to a first-
order transition at a tricritical point. On the other hand, as an extension of
the model, BC model with a random crystal field represents the critical
behavior of $\mathrm{{}^{3}He}-\mathrm{{}^{4}He}$ mixtures in a random media,
i.e., aerogel where $S=0$ and $S=\pm 1$ states represent $\mathrm{{}^{3}He}$
and $\mathrm{{}^{4}He}$ atoms, respectively aerogel1 ; aerogel2 . From the
theoretical point of view, BC model with a random crystal field (RCF) has been
studied by a variety of techniques such as cluster variational method (CVM)
aerogel2 , Bethe lattice approximation (BLA) albayrak , effective field theory
(EFT) kaneyoshi1 ; kaneyoshi2 ; kaneyoshi3 ; kaneyoshi4 ; yan , finite cluster
approximation (FCA) benyoussef1 ; ilkovic , mean field theory (MFT)
benyoussef2 ; borelli ; carneiro2 ; boccara ; carneiro1 ; bahmad , Monte Carlo
(MC) simulations puha , pair approximation (PA) lara , and renormalization
group (RG) method branco . Among these studies, EFT and MFT have been widely
used to investigate the thermal and magnetic properties of BC model with a RCF
distribution. For example, Kaneyoshi and Mielnicki kaneyoshi3  investigated
the phase diagram of the system for a honeycomb lattice by using EFT with
correlations and they found some important differences from the results
obtained by the standard MFT. Similarly, in a recent paper, Yan and Deng yan
considered the same model within the framework of EFT, and they derived the
expressions of magnetizations for honeycomb and square lattices. On the other
hand, in several studies based on MFT, benyoussef2 ; carneiro2  the authors
paid attention to the effects of crystal field dilution on the phase diagrams
of the system and they observed that the system may exhibit a reentrant
behavior, as well as first order phase transitions.

However, EFT and MFT studies mentioned above have some unsatisfactory results.
Namely, the results obtained by EFT are limited to second-order phase
transitions and tricritical points, and a detailed description of first-order
transitions has not been reported. Other than this, it is well known that
magnetic systems with dilute crystal fields exhibit qualitatively similar
characteristics when compared to site dilution problem of magnetic atoms. From
this point of view, for a BC model with diluted crystal fields, MFT predicts
that the phase transition temperature of the system will remain at a finite
value until zero concentration is reached. In this context, we believe that BC
model with RCF still deserves particular attention for investigating the
proper phase diagrams, especially the first-order transition lines that
include reentrant phase transition regions. Conventional EFT approximations
include spin-spin correlations resulting from the usage of the Van der Waerden
identities, and provide results that are superior to those obtained within the
traditional MFT. However, these conventional EFT approximations are not
sufficient enough to improve the results, due to the usage of a decoupling
approximation (DA) that neglects the correlations between different spins that
emerge when expanding the identities. Therefore, taking these correlations
into consideration will improve the results of conventional EFT
approximations. In order to overcome this point, recently we proposed an
approximation that takes into account the correlations between different spins
in the cluster of a considered lattice ak nc . Namely, an advantage of the
approximation method proposed by this study is that no decoupling procedure is
used for the higher-order correlation functions.

In this paper, we intent to investigate the effects of RCF distributions on
the phase diagrams of spin-1 BC model on 2D lattices, namely honeycomb $(q=3)$
and square $(q=4)$ lattices. For this purpose, the paper is organized as
follows: In Sec. II we briefly present the formulations. The results and
discussions are presented in Sec. III, and finally Sec. IV contains our
conclusions.

##  II Formulation

In this section, we give the formulation of the present study for a 2D lattice
which has $N$ identical spins arranged. We define a cluster on the lattice
which consists of a central spin labeled $S_{0}$, and $q$ perimeter spins
being the nearest neighbors of the central spin. The cluster consists of
$(q+1)$ spins being independent from the spin operator $\hat{S}$. The nearest-
neighbor spins are in an effective field produced by the outer spins, which
can be determined by the condition that the thermal average of the central
spin is equal to that of its nearest-neighbor spins. The Hamiltonian
describing our model is

$H=-J\sum_{\langle
i,j\rangle}S_{i}^{z}S_{j}^{z}-\sum_{i}D_{i}(S_{i}^{z})^{2},$ (1)

where the first term is a summation over the nearest-neighbor spins with
$S_{i}^{z}=\pm 1,0$ and the term $D_{i}$ on the second summation represents a
random crystal field, distributed according to a given probability
distribution. In this paper, we primarily deal with two kinds of probability
distributions, namely, a quenched diluted crystal field distribution and a
double peaked delta distribution which are given by Eqs. (2) and (3),
respectively as follows

$\displaystyle P(D_{i})$ $\displaystyle=$ $\displaystyle
p\delta(D_{i}-D)+(1-p)\delta(D_{i}),$ (2) $\displaystyle P(D_{i})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left\\{\delta[D_{i}-(D-\triangle)]+\delta[D_{i}-(D+\triangle)]\right\\}.$
(3)

where $p$ denotes the concentration of the spins on the lattice which are
influenced by a crystal field $D$.

We can construct the mathematical background of our model by using the
approximated spin correlation identities sa_barreto_et_al  by taking into
account random configurational averages

$\langle\langle\\{f_{i}\\}S_{i}^{z}\rangle\rangle_{r}=\left\langle\left\langle\\{f_{i}\\}\frac{\mathrm{Tr}_{i}\left(S_{i}^{z}\right)\exp{(-\beta
H_{i})}}{\mathrm{Tr}_{i}\exp{(-\beta H_{i})}}\right\rangle\right\rangle_{r},$
(4)
$\langle\langle\\{f_{i}\\}(S_{i}^{z})^{2}\rangle\rangle_{r}=\left\langle\left\langle\\{f_{i}\\}\frac{\mathrm{Tr}_{i}(S_{i}^{z})^{2}\exp{(-\beta
H_{i})}}{\mathrm{Tr}_{i}\exp{(-\beta H_{i})}}\right\rangle\right\rangle_{r},$
(5)

where $\beta=1/k_{B}T$, $\\{f_{i}\\}$ is an arbitrary function which is
independent of the spin variable $S_{i}$ and the inner $\langle...\rangle$ and
the outer $\langle...\rangle_{r}$ brackets represents the thermal and random
configurational averages, respectively.

In order to apply the differential operator technique kaneyoshi_honmura ;
kaneyoshi5 , we should separate the Hamiltonian (1) into two parts as
$H=H_{i}+H^{{}^{\prime}}$. Here, the effective Hamiltonian $H_{i}$ includes
all the contributions associated with the site $i$, and the other part
$H^{{}^{\prime}}$ does not depend on the site $i$.

$-H_{i}=ES_{i}^{z}+D_{i}\left(S_{i}^{z}\right)^{2},$ (6)

where $E=J\sum_{j}S_{j}^{z}$ is the local field on the site $i$. If we use the
matrix representations of the operators $S_{i}^{z}$ and $(S_{i}^{z})^{2}$ for
the spin-1 system then we can obtain the matrix form of Eq. (6)

$-H_{i}=\left(\begin{array}[]{ccc}E+D&0&0\\\ 0&0&0\\\ 0&0&-E+D\\\
\end{array}\right).$ (7)

Hereafter, we apply the differential operator technique in Eqs. (4) and (5)
with $\\{f_{i}\\}=1$. From Eq. (4) we obtain the following spin identity with
thermal and configurational averages of a central spin for a lattice with a
coordination number $q$ as

$\langle\langle
S_{0}^{z}\rangle\rangle_{r}=\left\langle\left\langle\prod_{j=1}^{q}\left[1+S_{j}^{z}\mathrm{sinh}(J\nabla)+(S_{j}^{z})^{2}\\{\mathrm{cosh}(J\nabla)-1\\}\right]\right\rangle\right\rangle_{r}F(x)|_{x=0}.$
(8)

The function $F(x)$ in Eq. (8) is defined by

$F(x)=\int dD_{i}P(D_{i})f(x,D_{i}),$ (9)

where

$\displaystyle f(x,D_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{\sum_{n=1}^{3}\exp(\beta\lambda_{n})}\sum_{n=1}^{3}\langle\varphi_{n}|S_{i}^{z}|\varphi_{n}\rangle\exp(\beta\lambda_{n}),$
$\displaystyle=$ $\displaystyle\frac{2\sinh(\beta x)}{2\cosh(\beta
x)+e^{-\beta D_{i}}}.$

In Eq. (II), $\lambda_{n}$ denotes the eigenvalues of $-H_{i}$ matrix in Eq.
(7), and $\varphi_{n}$ represents the eigenvectors corresponding to the
eigenvalues $\lambda_{n}$ of $-H_{i}$ matrix. With the help of Eq. (II), and
by using the distribution functions defined in Eqs. (2) and (3), the function
$F(x)$ in Eq. (9) can be easily calculated by numerical integration.
Hereafter, we will focus our attention on the construction of the correlation
functions, as well as magnetization and quadrupole moment identities of a
honeycomb lattice with $q=3$. A brief formulation of the fundamental spin
identities for a square lattice with $q=4$ can be found in Appendix A.

By expanding the right-hand side of Eq. (8) for a honeycomb lattice with
$q=3$, we get the longitudinal magnetization as

$\displaystyle m_{z}=\langle\langle S_{0}^{z}\rangle\rangle_{r}$
$\displaystyle=$ $\displaystyle l_{0}+3k_{1}\langle\langle
S_{1}\rangle\rangle_{r}+3(l_{1}-l_{0})\langle\langle
S_{1}^{2}\rangle\rangle_{r}+3l_{2}\langle\langle S_{1}S_{2}\rangle\rangle_{r}$
(11) $\displaystyle+6(k_{2}-k_{1})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}+3(l_{0}-2l_{1}+l_{3})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}$ $\displaystyle+k_{3}\langle\langle
S_{1}S_{2}S_{3}\rangle\rangle_{r}+3(l_{4}-l_{2})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+3(k_{1}-2k_{2}+k_{4})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-l_{0}+3l_{1}-3l_{3}+l_{5})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r},$

We note that, for the sake of simplicity, the superscript $z$ is omitted from
the correlation functions on the right-hand side of Eq. (11). The coefficients
in Eq. (11) are defined as follows:

$\displaystyle l_{0}=F(0),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ $\displaystyle
l_{1}=\mathrm{cosh}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ $ $\displaystyle k_{1}=\mathrm{sinh}(J\nabla)F(x)|_{x=0},$
$\displaystyle l_{2}=\mathrm{sinh}^{2}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ $ $\displaystyle
k_{2}=\mathrm{cosh}(J\nabla)\mathrm{sinh}(J\nabla)F(x)|_{x=0},$ $\displaystyle
l_{3}=\mathrm{cosh}^{2}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ $ $\displaystyle k_{3}=\mathrm{sinh}^{3}(J\nabla)F(x)|_{x=0},$
$\displaystyle
l_{4}=\mathrm{cosh}(J\nabla)\mathrm{sinh}^{2}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \
\ \ \ \ $ $\displaystyle
k_{4}=\mathrm{cosh}^{2}(J\nabla)\mathrm{sinh}(J\nabla)F(x)|_{x=0},$
$\displaystyle l_{5}=\mathrm{cosh}^{3}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ $ (12)

Next, the average value of the perimeter spin in the system can be written as
follows, and it is found as

$m_{1}=\langle\langle S_{\delta}^{z}\rangle\rangle_{r}=\langle\langle
1+S_{0}^{z}\mathrm{sinh}(J\nabla)+(S_{0}^{z})^{2}\\{\mathrm{cosh}(J\nabla)-1\\}\rangle\rangle_{r}F(x+\gamma)|_{x=0},$
(13) $\langle\langle
S_{1}\rangle\rangle_{r}=a_{1}\left(1-\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}\right)+a_{2}\langle\langle
S_{0}^{z}\rangle\rangle_{r}+a_{3}\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r},$
(14)

with the coefficients

$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle F(\gamma),$
$\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle\mathrm{sinh}(J\nabla)F(x+\gamma)|_{x=0},$ $\displaystyle a_{3}$
$\displaystyle=$ $\displaystyle\mathrm{cosh}(J\nabla)F(x+\gamma)|_{x=0},$ (15)

where $\gamma=(q-1)A$ is the effective field produced by the $(q-1)$ spins
outside the system and $A$ is an unknown parameter to be determined self-
consistently. In the effective-field approximation, the number of independent
spin variables describes the considered system. This number is given by the
relation $\nu=\langle\langle(S_{i}^{z})^{2S}\rangle\rangle_{r}$. As an example
for the spin-1 system, $2S=2$ which means that we have to introduce the
additional parameters $\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}$ and
$\langle\langle(S_{\delta}^{z})^{2}\rangle\rangle_{r}$ resulting from the
usage of the Van der Waerden identity for the spin-1 Ising system. With the
help of Eq. (5), quadrupolar moment of the central spin can be obtained as
follows

$\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}=\left\langle\left\langle\prod_{j=1}^{q}\left[1+S_{j}^{z}\mathrm{sinh}(J\nabla)+(S_{j}^{z})^{2}\\{\mathrm{cosh}(J\nabla)-1\\}\right]\right\rangle\right\rangle_{r}G(x)|_{x=0},$
(16)

where the function $G(x)$ is defined as

$G(x)=\int dD_{i}P(D_{i})g(x,D_{i}).$ (17)

Definition of the function $g(x,D_{i})$ in Eq. (17) is given as follows and
the expression in Eq. (II) can be evaluated by using the eigenvalues and
corresponding eigenvectors of the effective Hamiltonian matrix in Eq. (7).

$\displaystyle g(x,D_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{\sum_{n=1}^{3}\exp(\beta\lambda_{n})}\sum_{n=1}^{3}\langle\varphi_{n}|\left(S_{i}^{z}\right)^{2}|\varphi_{n}\rangle\exp(\beta\lambda_{n}),$
$\displaystyle=$ $\displaystyle\frac{2\cosh(\beta x)}{2\cosh(\beta
x)+e^{-\beta D_{i}}}.$

Hence, we get the quadrupolar moment by expanding the right-hand side of Eq.
(16)

$\displaystyle\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}$
$\displaystyle=$ $\displaystyle r_{0}+3n_{1}\langle\langle
S_{1}\rangle\rangle_{r}+3(r_{1}-r_{0})\langle\langle
S_{1}^{2}\rangle\rangle_{r}+3r_{2}\langle\langle S_{1}S_{2}\rangle\rangle_{r}$
(19) $\displaystyle+6(n_{2}-n_{1})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}+3(r_{0}-2r_{1}+r_{3})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+n_{3}\langle\langle
S_{1}S_{2}S_{3}\rangle\rangle_{r}$ $\displaystyle+3(r_{4}-r_{2})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}+3(n_{1}-2n_{2}+n_{4})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-r_{0}+3r_{1}-3r_{3}+r_{5})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r},$

with

$\displaystyle r_{0}=G(0),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ $\displaystyle
r_{1}=\mathrm{cosh}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ $ $\displaystyle n_{1}=\mathrm{sinh}(J\nabla)G(x)|_{x=0},$
$\displaystyle r_{2}=\mathrm{sinh}^{2}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ $ $\displaystyle
n_{2}=\mathrm{cosh}(J\nabla)\mathrm{sinh}(J\nabla)G(x)|_{x=0},$ $\displaystyle
r_{3}=\mathrm{cosh}^{2}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ $ $\displaystyle n_{3}=\mathrm{sinh}^{3}(J\nabla)G(x)|_{x=0},$
$\displaystyle
r_{4}=\mathrm{cosh}(J\nabla)\mathrm{sinh}^{2}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \
\ \ \ \ $ $\displaystyle
n_{4}=\mathrm{cosh}^{2}(J\nabla)\mathrm{sinh}(J\nabla)G(x)|_{x=0},$
$\displaystyle r_{5}=\mathrm{cosh}^{3}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ $ (20)

Corresponding to Eq. (13),

$\langle\langle(S_{\delta}^{z})^{2}\rangle\rangle_{r}=\langle\langle
1+S_{0}^{z}\mathrm{sinh}(J\nabla)+(S_{0}^{z})^{2}\\{\mathrm{cosh}(J\nabla)-1\\}\rangle\rangle_{r}G(x+\gamma),$
(21) $\langle\langle
S_{1}^{2}\rangle\rangle_{r}=b_{1}\left(1-\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}\right)+b_{2}\langle\langle
S_{0}^{z}\rangle\rangle_{r}+b_{3}\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}.$
(22)

where

$\displaystyle b_{1}$ $\displaystyle=$ $\displaystyle G(\gamma),$
$\displaystyle b_{2}$ $\displaystyle=$
$\displaystyle\mathrm{sinh}(J\nabla)G(x+\gamma)|_{x=0},$ $\displaystyle b_{3}$
$\displaystyle=$ $\displaystyle\mathrm{cosh}(J\nabla)\rangle
G(x+\gamma)|_{x=0}.$ (23)

Eqs. (11), (14), (19) and (22) are the fundamental spin identities of the
system. When the right-hand sides of Eqs. (8) and (16) are expanded, the
multispin correlation functions appear. The simplest approximation, and one of
the most frequently adopted is to decouple these correlations according to

$\left\langle\left\langle
S_{i}^{z}(S_{j}^{z})^{2}...S_{l}^{z}\right\rangle\right\rangle_{r}\cong\left\langle\left\langle
S_{i}^{z}\right\rangle\right\rangle_{r}\left\langle\left\langle(S_{j}^{z})^{2}\right\rangle\right\rangle_{r}...\left\langle\left\langle
S_{l}^{z}\right\rangle\right\rangle_{r},$ (24)

for $i\neq j\neq...\neq l$ tamura_kaneyoshi . The main difference of the
method used in this study from the other approximations in the literature
emerges in comparison with any decoupling approximation (DA) when expanding
the right-hand sides of Eqs. (8) and (16). In other words, one advantage of
the approximation method used in this study is that no uncontrolled decoupling
procedure is used for the higher-order correlation functions.

For spin-1 Ising system with $q=3$, taking Eqs. (11), (14), (19) and (22) as a
basis, we derive a set of linear equations of the spin identities. At this
point, we assume that (i) the correlations depend only on the distance between
the spins, (ii) the average values of a central spin and its nearest-neighbor
spin (it is labeled as the perimeter spin) are equal to each other with the
fact that, in the matrix representations of spin operator $\hat{S}$, the
spin-1 system has the properties $(S_{\delta}^{z})^{3}=S_{\delta}^{z}$ and
$(S_{\delta}^{z})^{4}=(S_{\delta}^{z})^{2}$. Thus, the number of the set of
linear equations obtained for the spin-1 Ising system with $q=3$ reduces to
twenty one, and the complete set is given in Appendix B.

If Eq. (B) is written in the form of a $21\times 21$ matrix and solved in
terms of the variables $x_{i}[(i=1,2,...,21)(e.g.,x_{1}=\langle\langle
S_{0}^{z}\rangle\rangle_{r},x_{2}=\langle\langle
S_{1}S_{0}\rangle\rangle_{r},...)]$ of the linear equations, all of the spin
correlation functions, as well as magnetizations and quadrupolar moments can
be easily determined as functions of the temperature and Hamiltonian
parameters. Since the thermal and configurational average of the central spin
is equal to that of its nearest-neighbor spins within the present method then
the unknown parameter $A$ can be numerically determined by the relation

$\langle\langle S_{0}^{z}\rangle\rangle_{r}=\langle\langle
S_{1}\rangle\rangle_{r}\qquad{\rm{or}}\qquad x_{1}=x_{4}.$ (25)

By solving Eq. (25) numerically at a given fixed set of Hamiltonian parameters
we obtain the parameter $A$. Then we use the numerical values of $A$ to obtain
the spin correlation functions which can be found from Eq. (B). Note that
$A=0$ is always the root of Eq. (25) corresponding to the disordered state of
the system. The nonzero root of $A$ in Eq. (25) corresponds to the long-range
ordered state of the system. Once the spin identities have been evaluated then
we can give the numerical results for the thermal and magnetic properties of
the system. Since the effective field $\gamma$ is very small in the vicinity
of $k_{B}T_{c}/J$, we can obtain the critical temperature for the fixed set of
Hamiltonian parameters by solving Eq. (25) in the limit of $\gamma\rightarrow
0$ then we can construct the whole phase diagrams of the system. Depending on
the values of Hamiltonian and crystal field distribution parameters, there may
be two solutions [i.e., two critical temperature values which satisfy Eq.
(25)] corresponding to the first (or second) and second-order phase-transition
points, respectively. We determine the type of the transition by looking at
the temperature dependence of magnetization for selected values of system
parameters.

##  III Results and Discussion

In this section, we will discuss the effect of the crystal field distributions
defined in Eqs. (2) and (3) on the global phase diagrams of the system where
the second and first order transitions are shown by solid and dashed curves,
respectively with tricritical points (shown by hollow circles) for honeycomb
$(q=3)$ and square $(q=4)$ lattices. Also, in order to clarify the type of the
transitions in the system, we will give the temperature dependence of the
order parameter.

###  III.1 Phase diagrams of the system with dilute crystal field

In this section, we illustrate the phase diagrams and magnetization curves of
the system with a dilute crystal field distribution defined in Eq.(2) where
crystal field $D$ is turned on, or turned off with probabilities $p$ and
$(1-p)$ on the lattice sites, respectively. In Figs. (1a) and (1c), we plot
the phase diagrams of the system in $(k_{B}T_{c}/J-D/J)$ plane for honeycomb
and square lattices with coordination numbers $q=3$ and $q=4$, respectively.
As seen in Figs. (1a) and (1c), phase diagrams of the system can be divided
into three parts with different concentration values $p$. For the curves in
the first group with $p<p^{*}$, the system always exhibits a second order
phase transition with a finite critical temperature $k_{B}T_{c}/J$ which
extent to $D/J\rightarrow-\infty$. If the concentration $p$ reaches its
critical value $p^{*}$ then the critical temperature depresses to zero.
Physical reason underlying this behavior can be explained as follows; when we
select sufficiently large negative crystal field values
$(\mathrm{i.e.}D/J\rightarrow-\infty)$, all of the spins in the system tend to
align in $S=0$ state. As $p$ increases starting from zero, the ratio of spins
which aligned in $S=0$ state increases, and therefore, magnetization weakens,
and accordingly, critical temperature of the system decreases. According to
our numerical results, the critical concentration value is obtained as
$p^{*}=0.3795$ for $q=3$ and $p^{*}=0.5875$ for $q=4$. In the second group of
the phase diagrams in Figs. (1a) and (1c), the system exhibits a reentrant
behavior of second order, whereas the curves in the third group, exhibit a
reentrant behavior of first order with a tricritical point at which a first
order transition line turns into a second order transition line. Besides, the
curves which exhibit a reentrant behavior of first (or second) order, depress
to zero at three successive values of crystal field $D/J=-3.0,-2.0,-1.0$.
Moreover, in $D/J\rightarrow\infty$ limit, the system behaves like spin-$1/2$
for $p=1.0$. In the case of $p\neq 0$, the ratio of spins that behave like
$S=\pm 1$ increases as $p$ increases. Therefore, for $0\leq p\leq 1.0$, all
transition lines have finite critical temperatures which increase with
increasing $p$ values for $D/J\rightarrow\infty$. At this point, we also note
that if we select $D/J=0$ in Eq. (2), all lattice sites expose to a crystal
field $D_{i}/J=0$ independent from $p$. Hence, all transition lines intersect
each other on the point $(D/J,k_{B}T_{c}/J)=(0,1.3022)$ for $q=3$, and
$(D/J,k_{B}T_{c}/J)=(0,1.9643)$ for $q=4$. Meanwhile, previous studies based
on EFT are not capable of obtaining first order transition lines of the
system. From this point of view, we see that our method improves the results
of the other EFT works and we take the conventional EFT method one step
forward by investigating the global phase diagrams, especially the first-order
transition lines that include reentrant phase transition regions.

Table 1: Critical concentration $p^{*}$ obtained by several methods and the
present work for honeycomb $(q=3)$ and square $(q=4)$ lattices.
---
Lattice | EFT-I kaneyoshi3 ; kaneyoshi4  |  EFT-II yan  |  MFTbenyoussef2 ; carneiro2  |  PA lara  |  Present Work
$q=3$ | 0.484 |  0.492 |  1.0 |  0.5 |  0.3795
$q=4$ | 0.604 |  0.610 |  1.0 |  0.667 |  0.5875

On the other hand, Figs. (1b) and (1d) shows the phase boundary in
$(k_{B}T_{c}/J-p)$ plane which separates the ferromagnetic and paramagnetic
phases with $D/J\rightarrow-\infty$. According to this figure, critical
temperature $k_{B}T_{c}/J$ of system decreases gradually, and ferromagnetic
region gets narrower as $p$ increases, and $k_{B}T_{c}/J$ value depresses to
zero at $p=p^{*}$. Such a behavior is an expected fact in dilution problems.
Numerical value of critical concentration $p^{*}$ for honeycomb $(q=3)$ and
square $(q=4)$ lattices is given in Table 1, and compared with the other works
in the literature. As seen in Table 1, numerical values of $p_{c}$ for $q=3$
and $q=4$ are new results in literature. Furthermore, MFT benyoussef2 ;
carneiro2  predicts that the system always has a finite critical temperature
and exists in a ferromagnetic state at lower temperatures in
$D/J\rightarrow-\infty$ limit, except that $p=1.0$. This artificial result can
be regarded as a failure of the MFT.

In Fig.(2), we plot the temperature dependencies of magnetization curves
corresponding to the phase diagrams depicted in Fig. (1) for $q=3$. As seen in
Fig. (2), as $p$ increases then critical temperature $k_{B}T_{c}/J$ values
decrease for $D/J<0$, except the reentrant phase transition temperatures which
occur at low temperatures. On the other hand, effect of increasing $p$ values
on the shape of magnetization curves depends on value of $D/J$. Namely, in
Figs. (2a) and (2b) we see that ground state saturation values of
magnetization curves decreases as $p$ increases for $D/J=-10.0$ and $-3.1$.
Moreover, for $D/J=-3.1$, magnetization curves of the system exhibit a broad
maximum at low temperatures for $p=0.37$, and a reentrant behavior of second
order for $p=0.4$. If we select $D/J=-2.5$ as in Fig. (2c), saturation values
of magnetization curves remain unchanged for $p=0,0.2,0.3$ and tend to
decrease for $p>0.3$. If $p$ increases further, a reentrant behavior of second
order appears for $p=0.53$, and we see a broad maximum at low temperatures for
$p=0.517$ and $0.5$ which tends to depress as $p$ decreases. This broad
maximum behavior of magnetization curves originates from the increase in the
number of spins directed parallel to the z-direction with increasing
temperature, due to the thermal agitation. For $D/J=-2.0$ in Fig. (2d),
magnetization curves saturates at $m=1$ at the ground state and reentrant
behavior disappears. If we increase $D/J$ further, for example for $D/J=-1.5$
(Fig. (2e)), another type of reentrant behavior occurs in the system in which
a first order transition is followed by a second order transition. Finally,
for sufficiently large positive values of crystal field, magnetization curves
always saturate at $m=1$ and the system always undergoes a second order phase
transition from a ferromagnetic phase to a paramagnetic phase with increasing
temperature, which can be seen in Fig. (2f) with $D/J=10.0$. As a common
property of the curves in Fig. (2), we see that effect of $p$ on the
saturation values, as well as temperature dependence of magnetization curves
strictly depend on the strength of $D/J$. Hence, according to us, the presence
of dilute crystal fields on the system should produce a competition effect on
the phase diagrams of the system. We also note that, although it has not been
shown in the present work, magnetization curves for $q=4$ corresponding to the
phase diagrams depicted in Fig. (1c) exhibit qualitatively similar behavior
with those of Fig. (2) with $q=3$.

As seen in Fig. (1), for a dilute crystal field distribution defined in Eq.
(2), the global phase diagrams which are plotted in $(k_{B}T_{c}/J-D/J)$
plane, as well as the phase boundaries in $(k_{B}T_{c}/J-p)$ plane for $q=3$
exhibit qualitatively similar characteristics when compared with those for
$q=4$. Hence, in order to examine the phase diagrams which are plotted in
$(k_{B}T_{c}/J-D/J)$ plane in Figs. (1a) and (1c) in detail, we plot the
evolution of the global phase diagrams in Fig. (3) only for $q=4$. From this
point of view, Fig. (3a) shows how the phase diagrams in Fig. (1c) evolve when
the concentration $p$ changes from 0.5 to 0.6. As seen in Fig. (3a), we
observe a second order phase transition line with a finite critical
temperature $k_{B}T_{c}/J$ which extent to $D/J\rightarrow-\infty$ for
$p=0.575$. If $p$ increases, namely for $p=0.580$ and $0.583$, we see that a
low temperature transition line arises between $-4.0<D/J<-3.0$, as well as a
high temperature phase boundary which extents to $D/J\rightarrow-\infty$. If
$p$ increases further, such as for $p=0.584$, $0.585$ and $0.587$, high
temperature phase boundary is gradually connected to the transition line which
arises between $-4.0<D/J<-3.0$, and the phase diagrams exhibit a bulge on the
right hand side of $(k_{B}T_{c}/J-D/J)$ plane, whereas another transition line
emerges within the range of $-\infty<D/J<-4.0$, which disappears for
$p>0.587$. Similarly, evolution of the phase diagrams in Fig. (1c) when the
concentration $p$ changes from 0.6 to 0.7 can be seen in Fig. (3b). As seen in
this figure, the curves for $p=0.60$, $0.62$, $0.64$, $0.66$ exhibit a
reentrant behavior of second order, while for $p=0.68$ reentrance disappears
and for $p=0.70$ and $0.71$ double reentrance with three successive second
order phase transitions occurs in a very narrow region of $D/J$. On the other
hand, increasing values of $p$ generates first order phase transitions with
tricritical points, as well as reentrant behavior of first order. This
phenomena is illustrated in Fig. (3c). From Fig. (3c), we see that, the second
order transition temperatures decrease as absolute value of $D/J$ increases,
and turn into first order transition lines at tricritical points. Evidently,
the phase diagrams change abruptly for $p\geq 0.7164$. Hence, the behavior of
the $p=0.7163$ curve is completely different from that of $p=0.7164$. Namely,
first order transition temperatures of the system for $p\leq 0.7163$ and
$p\geq 0.7164$ depress to zero at $D/J=-3.0$ and $D/J=-2.0$, respectively. In
order to investigate the phase transition features of the system further, we
should continue increasing the value of $p$. In Fig. (3d), we see that the
curves for $p=0.72$, $0.74$, $0.76$ and $0.78$ exhibit a reentrant behavior of
first order, whereas the curves with $p=0.80$, $0.82$, and $0.84$ exhibit
double reentrance with two first order and a second order transition
temperature.

###  III.2 Phase diagrams of the system with random crystal field

Next, in order to investigate the effect of the random crystal fields defined
in Eq.(3) on the thermal and magnetic properties of the system, we represent
the phase diagrams and corresponding magnetization curves for honeycomb
$(q=3)$ and square lattices $(q=4)$ throughout Figs. (4) and (8).

We note that random crystal field distribution given in Eq. (3) with $D/J=0$
corresponds to a bimodal distribution function, while for $\triangle/J=0$, we
obtain pure BC model with homogenous crystal field $D/J$. In Fig. (4), phase
diagrams of the system corresponding to the bimodal distribution function are
shown in $(k_{B}T_{c}/J-\triangle/J)$ plane. For a bimodal distribution, the
phase diagrams have symmetric shape with respect to $\triangle/J$ which comes
from the fact that $p=1/2$, and as seen in Fig. (4), transition temperatures
are second order, and it is clear that the system exhibit different
characteristic features depending on the coordination number $q$. Namely, for
$q=3$, transition temperature decreases with increasing $\triangle/J$ and
exhibits double reentrance with three second order phase transition
temperatures, then falls to zero at $\triangle/J=3.0$ (left panel in Fig.
(4)). On the other hand, as seen on the right panel in Fig. (4), as
$\triangle/J$ increases then the transition temperature of the system for
$q=4$ decreases and remains at a finite value for
$\triangle/J\rightarrow\infty$ which means that ferromagnetic exchange
interactions for $q=3$ are insufficient for the system to keep its
ferromagnetic order for $\triangle/J>3.0$, while for $q=4$ these interactions
are dominant in the system, and the presence of a disorder in the crystal
fields cannot destruct the ferromagnetic order.

At the same time, in order to see the effect of the random crystal fields with
$\triangle/J,D/J\neq 0$ on the phase diagrams and magnetization curves of the
system for $q=3$ and $4$, we plot the phase diagrams in $(k_{B}T_{c}/J-D/J)$
plane in Fig. (5) and variation of the corresponding magnetization curves with
temperature in Figs. (6) and (7), respectively. At first sight, it is obvious
that the phase diagrams in Fig.(5) represent evident differences in
qualitative manner with coordination number $q$. That is, as seen in Fig.
(5a), the curve corresponding to $\triangle/J=0$ represents the phase diagram
of pure BC model for a honeycomb lattice which exhibits a reentrant behavior
of first order with first and second order transition lines, as well as a
tricritical point. From Fig. (5a), we see that as $\triangle/J$ increases then
the tricritical point and first order transitions disappear, and the first
order reentrance turns into double reentrance with three transition
temperatures of second order, and phase transition lines shift to positive
crystal field direction without changing their shapes. On the other hand, the
situation is very different for a square lattice. Namely, as seen in Figs.
(5a) and (5b), $\triangle/J=0$ curves for $q=3$ and $q=4$ are qualitatively
identical to each other. However, as seen in Fig. (5b), for $\triangle/J\neq
0$, first order transition lines and tricritical points do not disappear from
the system for $q=4$, but shift to negative crystal field values. Besides, the
system does not exhibit double reentrance for $q=4$. Furthermore, for
$\triangle/J\geq 1.5$ in Fig. (5a), and $\triangle/J\geq 0$ in Fig. (5b), all
phase diagrams exhibit similar behavior as $D/J$ varies. Namely, critical
temperature $k_{B}T_{c}/J$ in Fig. (5a) reduces to zero at
$D/J=\triangle/J-3.0$. On the other hand, first order transition temperatures
in Fig. (5b), reduce to zero at $D/J=-\triangle/J$.

It is important to note that these observations are consistent with the
results shown in Figs. (1a) and (1c). In other words, the distribution
function given in Eq. (2) with $p=0.5$ and $D/J=2D_{0}/J$ is identical to Eq.
(3) for $\triangle/J=\pm D_{0}/J$ and $D/J=D_{0}/J$. For example, according to
Eq. (2), if we select $D_{0}/J=4.0$ with $p=0.5$, it means that half of the
spins on the lattice sites expose to a crystal field $D/J=0$, while a crystal
field given by $D/J=8.0$ acts on the other half of the spins. On the other
hand, if we select $\triangle/J=\pm D_{0}/J$ and $D/J=D_{0}/J$ by using Eq.
(3), we generate the same distribution again. Hence, we expect to get the same
results in Figs. (1) and (5) for these system parameters. For instance, for
$D_{0}/J=4.0$ in Fig. (1a), we get $D/J=8.0$, and the system exhibits a
ferromagnetic order in the ground state, which can also be seen in Fig. (5a)
with a critical temperature $k_{B}T_{c}/J=1.4395$. These conditions are also
valid for $q=4$, and for the whole temperature region on the phase diagrams.
Therefore, the state (para-or ferro), as well as thermal and magnetic
properties of a selected $(k_{B}T/J,D/J)$ point with respect to $p=0.5$ curves
in Figs. (1a) and (1c) is identical to the state of a point $(k_{B}T/J,D/2J)$
in Figs. (5a) and (5b) with respect to the curve $\triangle/J=\pm D/J$,
respectively. Moreover, the qualitative differences between Figs. (5a) and
(5b) mentioned above are strongly related to the percolation threshold value
of the lattice. Namely, distribution function Eq.(3) is valid only for
$p=0.5$. However, as seen in Table 1, we obtain $p_{c}<0.5$ for $q=3$, and
$p_{c}>0.5$ for $q=4$.

In Fig. (6), we examine the temperature dependence of magnetization curves for
$q=3$, corresponding to the phase diagrams shown in Fig. (5a) with $D/J=-1.0$.
Fig. (6a), shows how the temperature dependence of magnetization curves evolve
when $\triangle/J$ changes. According to Fig. (6a), magnetization curves
saturate at a partially ordered state at low temperatures. Besides, for
$\triangle/J=1.68,1.74$ and $1.80$, the system undergoes three successive
phase transitions of second order, which confirms the existence of double
reentrance. Similarly, Fig. (6b) shows how the shape of the magnetization
curves change as $D/J$ changes for constant $\triangle/J=6.0$. As seen in Fig.
(6b), magnetization curves exhibit a second order phase transition from a
ferromagnetic (fully ordered) to a paramagnetic phase at certain values of
crystal field, namely at $D/J=4.0,4.1$ and $4.4$, whereas for
$D/J=3.3,3.5,3.7$ and $3.9$ the system can only achieve a partially ordered
phase. In addition, the curves corresponding to $D/J=3.5,3.7$ and $3.9$
exhibit a broad maximum at low temperatures, and then decrease as the
temperature increases, whereas for $D/J=3.3$, we observe double reentrance.
Fig. (7) shows the magnetization curves for $q=4$, corresponding to the phase
diagrams shown in Fig. (5b). In Fig. (7a), we see that magnetization curves
exhibit a second order phase transition from a paramagnetic phase to a fully
ordered ferromagnetic phase for $\triangle/J=0$ and $5.0$. On the other hand,
the curves corresponding to $D/J=5.5$, $5.8$ and $6.0$, saturate at a
partially ordered state at low temperatures, and exhibit a broad maximum with
increasing temperature which depresses gradually as $\triangle/J$ increases,
then fall rapidly at a second order phase transition temperature. The broad
maximum behavior observed in these curves disappears for $\triangle/J=8.0$.
Additionally, Fig. (7b) represents the magnetization versus temperature curves
for $q=4$ with with $D/J=-4.0$. In Fig. (7b), it is clearly evident that, at
low temperatures, the system saturates at a partially ordered phase for
$\triangle/J=4.0$, $5.0$, $6.0$ and $7.0$, while for $\triangle/J=3.7$ and
$3.8$, a reentrant behavior of first order occurs. Again we see that, there is
a competition between ferromagnetic exchange interactions and disorder effects
in crystal fields which determines the saturation values and temperature
dependence of magnetization curves of the system.

Finally, dependence of magnetization of the system on the crystal field
$\triangle/J$ for fixed temperature values $k_{B}T/J=0.01,0.05,0.1$ and $0.2$
with $D/J=2.0$ is shown in Fig. (8) for $q=3$ and $4$, respectively. We see
that at sufficiently low temperatures such as $k_{B}T/J=0.01$, magnetization
curves exhibit three phases for $q=3$. A first order transition is
characterized by a gap in this figure. On the left panel in Fig. (8), which is
plotted for $q=3$, we observe two successive first order phase transitions for
$k_{B}T/J=0.01$. The first one is from the fully ordered ferromagnetic phase
($m=1.0$) to the partly ordered phase ($m=0.47$), and the other is from partly
ordered phase to disordered phase ($m=0.0$). On the other hand, according to
the right panel in Fig. (8), the system can not reach a paramagnetic phase at
the ground state for $q=4$.Hence, we observe two phases. Namely, for
$k_{B}T/J=0.01$, a first order phase transition from a fully ordered phase
($m=1.0$) to a partly ordered phase $(m=0.59)$. Then, saturation magnetization
of the partly ordered phase reduces continuously to ($m=0.527$). Moreover, the
first order transitions disappear with increasing temperatures, both for $q=3$
and $4$. Since in Figs. (5a) and (5b), all phase diagrams exhibit similar
behavior as $D/J$ varies, behavior of magnetization curves in Fig. (8) should
be the same for different $D/J$ values. Namely, the left panel of Fig. (8) can
be regarded as the variation of magnetization with $\triangle/J$ within the
range $D/J+1.0<\triangle/J<D/J+4.0$ for $q=3$. for a selected value of $D/J$.
It is possible to observe the similar behavior on the right panel in Fig. (8).
This general behavior is an expected result, since the behavior of the system
depends on the relation between $\triangle/J$ and $D/J$ parameters for a given
temperature, not their values independently.

##  IV Conclusions

In this work, we have studied the phase diagrams of a spin-1 Blume-Capel model
with diluted and random crystal field interactions on two dimensional
lattices. We have introduced an effective-field approximation that takes into
account the correlations between different spins in the cluster of a
considered lattice and examined the phase diagrams as well as magnetization
curves of the system for different types of crystal field distributions,
namely, dilute crystal fields and a double peaked delta distribution, given by
Eqs. (2) and (3), respectively.

For dilute crystal fields, we have given a detailed exploration of the global
phase diagrams of the system in $k_{B}T_{c}/J-D/J$ plane with the second and
first order transitions, as well as tricritical points. We have also shown
that the system with dilute crystal fields exhibits a percolation threshold
value $p_{c}$ which can not be predicted by standard MFA. In addition, we have
observed multi-reentrant phase transitions for specific set of system
parameters.

On the other hand, we have investigated the effect of the random crystal field
distribution characterized by two crystal field parameters $D/J$ and
$\triangle/J$ on the phase diagrams of the system. As a limited case, we have
also focused on a bimodal distribution with $D/J=0$. Particulary, we have
reported the following observations for a bimodal distribution: It has been
found that the phase diagrams have symmetric shape with respect to
$\triangle/J$ which comes from the fact that $p=1/2$. The transition
temperatures are of second order, and the system exhibit different
characteristic features depending on the coordination number $q$. Besides, we
have realized that the system may exhibit clear distinctions in qualitative
manner with coordination number $q$ for random crystal fields with
$\triangle/J,D/J\neq 0$. Moreover, we have discussed a competition effect
which arises from the presence of dilution, as well as random crystal fields,
and we have observed that saturation values of the magnetization curves are
strongly related to these effects.

As a result, we can conclude that all of the points mentioned above show that
our method improves the conventional EFT methods based on decoupling
approximation. Therefore, we hope that the results obtained in this work may
be beneficial from both theoretical and experimental points of view.

## Acknowledgements

One of the authors (Y.Y.) would like to thank the Scientific and Technological
Research Council of Turkey (TÜBİTAK) for partial financial support. This work
has been completed at Dokuz Eylül University, Graduate School of Natural and
Applied Sciences. Partial financial support from SRF (Scientific Research
Fund) of Dokuz Eylül University (2009.KB.FEN.077) (H.P.) is also acknowledged.

##  Appendix A Fundamental correlation functions of the system for a square
lattice

Magnetization of the central spin for a square lattice is given as follows

$\displaystyle\left\langle\langle S_{0}^{z}\right\rangle\rangle$
$\displaystyle=$ $\displaystyle\mu_{0}+4c_{1}\langle\langle
S_{1}\rangle\rangle_{r}+4(\mu_{2}-\mu_{0})\langle\langle
S_{1}^{2}\rangle\rangle_{r}$ (26) $\displaystyle+6\mu_{1}\langle\langle
S_{1}S_{2}\rangle\rangle_{r}+12(c_{2}-c_{1})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}$
$\displaystyle+6(\mu_{0}-2\mu_{2}+\mu_{3})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+4c_{3}\langle\langle
S_{1}S_{2}S_{3}\rangle\rangle_{r}$
$\displaystyle+12(\mu_{4}-\mu_{1})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}+12(c_{4}-2c_{2}+c_{1})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+4(\mu_{5}-3\mu_{3}+3\mu_{2}-\mu_{0})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}+\mu_{8}\langle\langle
S_{1}S_{2}S_{3}S_{4}\rangle\rangle_{r}$
$\displaystyle+4(c_{5}-c_{3})\langle\langle
S_{1}S_{2}S_{3}S_{4}^{2}\rangle\rangle_{r}+6(\mu_{1}-2\mu_{4}+\mu_{6})\langle\langle
S_{1}S_{2}S_{3}^{2}S_{4}^{2}\rangle\rangle_{r}$
$\displaystyle+4(c_{6}-3c_{4}+3c_{2}-c_{1})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}S_{4}^{2}\rangle\rangle_{r}$
$\displaystyle+(\mu_{0}-4\mu_{2}+6\mu_{3}-4\mu_{5}+\mu_{7})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}S_{4}^{2}\rangle\rangle_{r},$

where the coefficients are given by

$\displaystyle\mu_{0}=F(0),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ $ $\displaystyle\mu_{1}=\sinh^{2}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \
\ \ \ \ \ \ \ $ $\displaystyle c_{1}=\sinh(J\nabla)F(x)_{x=0},$
$\displaystyle\mu_{2}=\cosh(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \ \ $
$\displaystyle c_{2}=\sinh(J\nabla)\cosh(J\nabla)F(x)_{x=0},$
$\displaystyle\mu_{3}=\cosh^{2}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \
$ $\displaystyle c_{3}=\sinh^{3}(J\nabla)F(x)_{x=0},$
$\displaystyle\mu_{4}=\sinh^{2}(J\nabla)\cosh(J\nabla)F(x)|_{x=0},$
$\displaystyle c_{4}=\cosh^{2}(J\nabla)\sinh(J\nabla)F(x)_{x=0},$
$\displaystyle\mu_{5}=\cosh^{3}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \
$ $\displaystyle c_{5}=\sinh^{3}(J\nabla)\cosh(J\nabla)F(x)_{x=0},$
$\displaystyle\mu_{6}=\sinh^{2}(J\nabla)\cosh^{2}(J\nabla)F(x)_{x=0},$
$\displaystyle c_{6}=\cosh^{3}(J\nabla)\sinh(J\nabla)F(x)_{x=0},$
$\displaystyle\mu_{7}=\cosh^{4}(J\nabla)F(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \
$ $\displaystyle\mu_{8}=\sinh^{4}(J\nabla)F(x)_{x=0}.\ \ \ \ \ \ \ \ \ \ \ \ \
\ $

Quadrupolar moment corresponding to equation (B) defined as

$\displaystyle\left\langle\langle(S_{0}^{z})^{2}\right\rangle\rangle$
$\displaystyle=$ $\displaystyle\rho_{0}+4\eta_{1}\langle\langle
S_{1}\rangle\rangle_{r}+4(\rho_{2}-\rho_{0})\langle\langle
S_{1}^{2}\rangle\rangle_{r}$ (27) $\displaystyle+6\rho_{1}\langle\langle
S_{1}S_{2}\rangle\rangle_{r}+12(\eta_{2}-\eta_{1})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}$
$\displaystyle+6(\rho_{0}-2\rho_{2}+\rho_{3})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+4\eta_{3}\langle\langle
S_{1}S_{2}S_{3}\rangle\rangle_{r}$
$\displaystyle+12(\rho_{4}-\rho_{1})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}+12(\eta_{4}-2\eta_{2}+\eta_{1})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+4(\rho_{5}-3\rho_{3}+3\rho_{2}-\rho_{0})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}+\rho_{8}\langle\langle
S_{1}S_{2}S_{3}S_{4}\rangle\rangle_{r}$
$\displaystyle+4(\eta_{5}-\eta_{3})\langle\langle
S_{1}S_{2}S_{3}S_{4}^{2}\rangle\rangle_{r}+6(\rho_{1}-2\rho_{4}+\rho_{6})\langle\langle
S_{1}S_{2}S_{3}^{2}S_{4}^{2}\rangle\rangle_{r}$
$\displaystyle+4(\eta_{6}-3\eta_{4}+3\eta_{2}-\eta_{1})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}S_{4}^{2}\rangle\rangle_{r}$
$\displaystyle+(\rho_{0}-4\rho_{2}+6\rho_{3}-4\rho_{5}+\rho_{7})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}S_{4}^{2}\rangle\rangle_{r},$

where

$\displaystyle\rho_{0}=G(0),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ $ $\displaystyle\rho_{1}=\sinh^{2}(J\nabla)G(x)|_{x=0},\ \ \ \ \
\ \ \ \ \ \ \ \ $ $\displaystyle\eta_{1}=\sinh(J\nabla)G(x)_{x=0},$
$\displaystyle\rho_{2}=\cosh(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \ \ $
$\displaystyle\eta_{2}=\sinh(J\nabla)\cosh(J\nabla)G(x)_{x=0},$
$\displaystyle\rho_{3}=\cosh^{2}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \
$ $\displaystyle\eta_{3}=\sinh^{3}(J\nabla)G(x)_{x=0},$
$\displaystyle\rho_{4}=\sinh^{2}(J\nabla)\cosh(J\nabla)G(x)|_{x=0},$
$\displaystyle\eta_{4}=\cosh^{2}(J\nabla)\sinh(J\nabla)G(x)_{x=0},$
$\displaystyle\rho_{5}=\cosh^{3}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \
$ $\displaystyle\eta_{5}=\sinh^{3}(J\nabla)\cosh(J\nabla)G(x)_{x=0},$
$\displaystyle\rho_{6}=\sinh^{2}(J\nabla)\cosh^{2}(J\nabla)G(x)_{x=0},$
$\displaystyle\eta_{6}=\cosh^{3}(J\nabla)\sinh(J\nabla)G(x)_{x=0},$
$\displaystyle\rho_{7}=\cosh^{4}(J\nabla)G(x)|_{x=0},\ \ \ \ \ \ \ \ \ \ \ \ \
$ $\displaystyle\rho_{8}=\sinh^{4}(J\nabla)G(x)_{x=0}.\ \ \ \ \ \ \ \ \ \ \ \
\ \ $

Finally, perimeter spin identities are as follows

$\displaystyle\langle\langle S_{1}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle\alpha_{1}(1-\langle\langle(S_{0})^{2}\rangle\rangle_{r})+\alpha_{2}\langle\langle
S_{0}\rangle\rangle_{r}+\alpha_{3}\langle\langle(S_{0})^{2}\rangle\rangle_{r},$
(28) $\displaystyle\langle\langle S_{1}^{2}\rangle\rangle_{r}$
$\displaystyle=$ $\displaystyle\omega_{1}+\omega_{2}\langle\langle
S_{0}\rangle\rangle_{r}+(\omega_{3}-\omega_{1})\langle\langle
S_{0}^{2}\rangle\rangle_{r}$ (29)

with the coefficients

$\displaystyle\alpha_{1}=F(\gamma)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $
$\displaystyle\omega_{1}=G(\gamma)|_{x=0}$
$\displaystyle\alpha_{2}=\sinh(J\nabla)F(x+\gamma)$
$\displaystyle\omega_{2}=\sinh(J\nabla)G(x+\gamma)|_{x=0}$
$\displaystyle\alpha_{3}=\cosh(J\nabla)F(x+\gamma)$
$\displaystyle\omega_{3}=\cosh(J\nabla)G(x+\gamma)|_{x=0}$

where $\gamma=(q-1)A$ with $q=4$, and the functions $F(x)$ and $G(x)$ are
defined in equations (9) and (17).

##  Appendix B The complete set of twenty one linear equations of a honeycomb
lattice

$\displaystyle\langle\langle S_{0}^{z}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle l_{0}+3k_{1}\langle\langle
S_{1}\rangle\rangle_{r}+3(l_{1}-l_{0})\langle\langle
S_{1}^{2}\rangle\rangle_{r}+3l_{2}\langle\langle S_{1}S_{2}\rangle\rangle_{r}$
$\displaystyle+6(k_{2}-k_{1})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}+3(l_{0}-2l_{1}+l_{3})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}$ $\displaystyle+k_{3}\langle\langle
S_{1}S_{2}S_{3}\rangle\rangle_{r}+3(l_{4}-l_{2})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+3(k_{1}-2k_{2}+k_{4})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-l_{0}+3l_{1}-3l_{3}+l_{5})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{0}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle(3l_{1}-2l_{0})\langle\langle
S_{1}\rangle\rangle_{r}+3k_{1}\langle\langle
S_{1}^{2}\rangle\rangle_{r}+3(l_{0}-2l_{1}+l_{2}+l_{3})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}$ $\displaystyle+6(k_{2}-k_{1})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+k_{3}\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-l_{0}+3l_{1}-3l_{2}-3l_{3}+3l_{4}+l_{5})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+3(k_{1}-2k_{2}+k_{4})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{2}S_{0}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle({l_{0}}-3{l_{1}}+3{l_{2}}+3{l_{3}})\langle\langle{S_{1}}{S_{2}}\rangle\rangle_{r}+(6{k_{2}}-3{k_{1}})\langle\langle{S_{1}}S_{2}^{2}\rangle\rangle_{r}$
$\displaystyle+(-{l_{0}}+3{l_{1}}-3{l_{2}}-3{l_{3}}+3{l_{4}}+{l_{5}})\langle\langle{S_{1}}{S_{2}}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(3{k_{1}}-6{k_{2}}+{k_{3}}+3{k_{4}})\langle\langle{S_{1}}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle\langle\langle S_{1}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle
a_{1}(1-\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r})+a_{2}\langle\langle
S_{0}^{z}\rangle\rangle_{r}+a_{3}\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}$
$\displaystyle\langle\langle S_{1}S_{2}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle a_{1}\langle\langle S_{1}\rangle\rangle_{r}+a_{2}\langle\langle
S_{0}S_{1}\rangle\rangle_{r}+(a_{3}-a_{1})\langle\langle
S_{1}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{2}S_{3}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
a_{1}\langle\langle S_{1}S_{2}\rangle\rangle_{r}+a_{2}\langle\langle
S_{0}S_{1}S_{2}\rangle\rangle_{r}+(a_{3}-a_{1})\langle\langle
S_{1}S_{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{1}(1-\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r})+b_{2}\langle\langle
S_{0}^{z}\rangle\rangle_{r}+b_{3}\langle\langle(S_{0}^{z})^{2}\rangle\rangle_{r}$
$\displaystyle\langle\langle S_{1}S_{2}^{2}\rangle\rangle_{r}$
$\displaystyle=$ $\displaystyle b_{1}\langle\langle
S_{1}\rangle\rangle_{r}+b_{2}\langle\langle
S_{0}S_{1}\rangle\rangle_{r}+(b_{3}-b_{1})\langle\langle
S_{1}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{1}\langle\langle S_{1}^{2}\rangle\rangle_{r}+b_{2}\langle\langle
S_{0}S_{1}^{2}\rangle\rangle_{r}+(b_{3}-b_{1})\langle\langle
S_{1}^{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{0}S_{1}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{3}\langle\langle S_{0}\rangle\rangle_{r}+b_{2}\langle\langle
S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{0}S_{1}S_{2}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{3}\langle\langle S_{0}S_{1}\rangle\rangle_{r}+b_{2}\langle\langle
S_{1}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{0}S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{3}\langle\langle S_{0}S_{1}^{2}\rangle\rangle_{r}+b_{2}\langle\langle
S_{1}^{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{1}\langle\langle S_{1}S_{2}\rangle\rangle_{r}+b_{2}\langle\langle
S_{0}S_{1}S_{2}\rangle\rangle_{r}+(b_{3}-b_{1})\langle\langle
S_{1}S_{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{1}\langle\langle S_{1}S_{2}^{2}\rangle\rangle_{r}+b_{2}\langle\langle
S_{0}S_{1}S_{2}^{2}\rangle\rangle_{r}+(b_{3}-b_{1})\langle\langle
S_{1}S_{2}^{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
b_{1}\langle\langle S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+b_{2}\langle\langle
S_{0}S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+(b_{3}-b_{1})\langle\langle
S_{1}^{2}S_{2}^{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle=$ $\displaystyle
r_{0}+3n_{1}\langle\langle
S_{1}\rangle\rangle_{r}+3(r_{1}-r_{0})\langle\langle
S_{1}^{2}\rangle\rangle_{r}+3r_{2}\langle\langle S_{1}S_{2}\rangle\rangle_{r}$
$\displaystyle+6(n_{2}-n_{1})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}+3(r_{0}-2r_{1}+r_{3})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+n_{3}\langle\langle
S_{1}S_{2}S_{3}\rangle\rangle_{r}$ $\displaystyle+3(r_{4}-r_{2})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}+3(n_{1}-2n_{2}+n_{4})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-r_{0}+3r_{1}-3r_{3}+r_{5})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle(3r_{1}-2r_{0})\langle\langle
S_{1}\rangle\rangle_{r}+3n_{1}\langle\langle
S_{1}^{2}\rangle\rangle_{r}+(3r_{2}+3r_{0}-6r_{1}+3r_{3})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}$ $\displaystyle+6(n_{2}-n_{1})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+n_{3}\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-r_{0}+3r_{1}-3r_{2}-3r_{3}+3r_{4}+r_{5})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+3(n_{1}-2n_{2}+n_{4})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}^{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle(3r_{1}-2r_{0})\langle\langle
S_{1}^{2}\rangle\rangle_{r}+3n_{1}\langle\langle
S_{1}\rangle\rangle_{r}+(3r_{2}+3r_{0}-6r_{1}+3r_{3})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}$
$\displaystyle+6(n_{2}-n_{1})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}+(3n_{1}-6n_{2}+n_{3}+3n_{4})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-r_{0}+3r_{1}-3r_{2}-3r_{3}+3r_{4}+r_{5})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle(r_{0}-3r_{1}+3r_{2}+3r_{3})\langle\langle
S_{1}S_{2}\rangle\rangle_{r}+(-3n_{1}+6n_{2})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}$
$\displaystyle(-r_{0}+3r_{1}-3r_{2}-3r_{3}+3r_{4}+r_{5})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(3n_{1}-6n_{2}+n_{3}+3n_{4})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}S_{2}^{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle(r_{0}-3r_{1}+3r_{2}+3r_{3})\langle\langle
S_{1}S_{2}^{2}\rangle\rangle_{r}+(-3n_{1}+6n_{2})\langle\langle
S_{1}S_{2}\rangle\rangle_{r}$
$\displaystyle(-r_{0}+3r_{1}-3r_{2}-3r_{3}+3r_{4}+r_{5})\langle\langle
S_{1}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(3n_{1}-6n_{2}+n_{3}+3n_{4})\langle\langle
S_{1}S_{2}S_{3}^{2}\rangle\rangle_{r}$ $\displaystyle\langle\langle
S_{1}^{2}S_{2}^{2}S_{0}^{2}\rangle\rangle_{r}$ $\displaystyle=$
$\displaystyle({r_{0}}-3{r_{1}}+3{r_{2}}+3{r_{3}})\langle\langle
S_{1}^{2}S_{2}^{2}\rangle\rangle_{r}+(-3{n_{1}}+6{n_{2}})\langle\langle{S_{1}}S_{2}^{2}\rangle\rangle_{r}$
$\displaystyle+(3{n_{1}}-6{n_{2}}+{n_{3}}+3{n_{4}})\langle\langle{S_{1}}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}$
$\displaystyle+(-{r_{0}}+3{r_{1}}-3{r_{2}}-3{r_{3}}+3{r_{4}}+{r_{5}})\langle\langle
S_{1}^{2}S_{2}^{2}S_{3}^{2}\rangle\rangle_{r}.$

## References

  * (1)
  * (2) M. Blume, Phys. Rev. 141, 517 (1966).
  * (3) H. W. Capel, Physica 32, 966 (1966).
  * (4) A. Maritan, M. Cieplak, M. R. Swift, F. Toigo, J.R. Banavar, Phys. Rev. Lett. 69, 221 (1992).
  * (5) C. Buzano, A. Maritan, A. Pelizzola, J. Phys. Condens. Matter 6, 327 (1994).
  * (6) E. Albayrak, Physica A 390, 1529 (2011).
  * (7) T. Kaneyoshi, J. Phys. C 19, L557 (1986).
  * (8) T. Kaneyoshi, J. Phys. C 21, L679 (1988).
  * (9) T. Kaneyoshi and J. Mielnicki, J. Phys. Condens. Matter 2, 8773 (1990).
  * (10) T. Kaneyoshi, Phys. Status Solidi B 170, 313 (1992).
  * (11) S. L. Yan and L. L. Deng, Physica A 308, 301 (2002).
  * (12) A. Benyoussef and H. Ez-Zahraouy, J. Phys. Condens. Matter 6, 3411 (1994).
  * (13) V. Ilkovic, Phys. Status Solidi B 192, K7 (1995).
  * (14) A. Benyoussef, T. Biaz, M. Saber, and M. Touzani, J. Phys. C 20, 5349 (1987).
  * (15) M. E. S. Borelli and C. E. I Carneiro, Physica A 230, 249 (1996).
  * (16) C. E. I Carneiro, V. B. Henriques, and S. R. Salinas, J. Phys. Condens. Matter 1, 3687 (1989).
  * (17) N. Boccara, A. El Kenz, and M. Saber, J. Phys. Condens. Matter 1, 5721 (1989).
  * (18) C. E. I Carneiro, V. B. Henriques, and S. R. Salinas, J. Phys. A Math. Gen. 23, 3383 (1990).
  * (19) L. Bahmad, A. Benyoussef, and A. El Kenz, J. Magn. Magn. Mater 320, 397 (2008).
  * (20) I. Puha and H. T. Diep, J. Magn. Magn. Mater. 224, 85 (2001).
  * (21) D. P Lara and J. A. Plascak, Physica A 260, 443 (1998).
  * (22) N. S. Branco and B. M. Boechat, Phys. Rev. B 56, 11673 (1997).
  * (23) Ü. Akıncı, Y. Yüksel, and H. Polat, Phys. Rev. E 83, 061103 (2011).
  * (24) F. C. SáBarreto, I. P. Fittipaldi, B. Zeks, Ferroelectrics 39, 1103 (1981).
  * (25) R. Honmura and T. Kaneyoshi, J. Phys. C 12, 3979 (1979).
  * (26) T. Kaneyoshi, Acta Phys. Pol. A 83, 703 (1993).
  * (27) I. Tamura, T. Kaneyoshi, Prog. Theor. Phys. 66, 1892 (1981).

Figure 1: (a) Phase diagrams of the system for $q=3$ in a $(k_{B}T_{c}/J-D/J)$
plane corresponding to dilute crystal field distribution defined in Eq. (2).
The solid and dashed lines correspond to second- and first-order phase
transitions, respectively. The open circles denote the tricritical points, and
the numbers on each curve represent the value of concentration $p$. (b) Phase
diagrams of the system for $q=3$ in a $(k_{B}T_{c}/J-p)$ plane with a selected
value of the crystal field $D/J=-10.0$. (c) Phase diagrams of the system for
$q=4$ in a $(k_{B}T_{c}/J-D/J)$ plane corresponding to dilute crystal field
distribution defined in Eq. (2). The solid and dashed lines correspond to
second- and first-order phase transition, respectively. The open circles refer
to the tricritical points, and the numbers on each curve represent the value
of concentration $p$. (d) Phase diagrams of the system for $q=4$ in a
$(k_{B}T_{c}/J-p)$ plane with a selected value of the crystal field
$D/J=-10.0$.
Figure 2: (Color online) Temperature dependence of magnetization corresponding
to Fig. (1a) with some selected values of crystal field. (a) $D/J=-10.0$, (b)
$D/J=-3.1$, (c) $D/J=-2.5$, (d) $D/J=-2.0$, (e) $D/J=-1.5$, and (f)
$D/J=10.0$. The numbers on each curve denote the value of concentration $p$.
The solid and dashed lines correspond to second- and first-order phase
transitions, respectively.

Figure 3: (Color online) Evolution of the phase diagrams corresponding to Fig.
(1c). The numbers on each curve denote the value of concentration $p$. The
solid and dashed lines correspond to second- and first-order phase
transitions, respectively. The open circles indicate the tricritical points.
Figure 4: Phase diagrams of the system in a $(k_{B}T_{c}/J-\triangle/J)$ plane
for a bimodal crystal field distribution corresponding to Eq. (3) with
$D/J=0.0$. Left and right-hand side panels are plotted for $q=3$ and $q=4$,
respectively.

Figure 5: (Color online) Phase diagrams of the system in a
$(k_{B}T_{c}/J-D/J)$ plane corresponding to random crystal field distribution
defined in Eq. (3) for (a) $q=3$, (b) $q=4$. The numbers on each curve denote
the value of $\triangle/J$. The open circles represent the tricritical points,
and the solid and dashed lines correspond to second- and first-order phase
transitions, respectively.

Figure 6: (Color online) (a) Temperature dependence of magnetization curves
corresponding to Fig. (5a) for $q=3$ (a) with $D/J=-1.0$ and for some selected
values of $\triangle/J$, and (b) with $\triangle/J=6.0$ and for some selected
values of $D/J$.

Figure 7: (Color online) (a) Temperature dependence of magnetization curves
for $q=4$ corresponding to Fig. (5b) for (a) $D/J=2.0$ and (b) $D/J=-4.0$ with
some selected values of $\triangle/J$.
Figure 8: (Color online) Variation of the magnetization curves as a function
of $\triangle/J$ with $D/J=2.0$. Left and right panels are plotted for $q=3$
and $4$, respectively. Each curve is plotted for different temperature values,
namely $k_{B}T/J=0.01,0.05,0.1$ and $0.2$.