# Study of the Ashkin Teller model with spins $S$ = $1$ and $\sigma$ = $3/2$
subjected to different crystal fields using the Monte-Carlo method

Z. Amimer${}^{\textsf{\footnotesize{\mbox{?}}}}$, S.
Bekhechi${}^{\textsf{\footnotesize{\mbox{?}}}}$, B.N.
Brahmi${}^{\textsf{\footnotesize{\mbox{?}}}}$, R.
Boudefla${}^{\textsf{\footnotesize{\mbox{?}}}}$, H. Ez-
Zahraouy${}^{\textsf{\footnotesize{\mbox{?}}}}$, A.
Rachadi${}^{\textsf{\footnotesize{\mbox{?}}}}$

(Received March 3, 2020, in final form June 18, 2020)

###### Abstract

Using the Monte-Carlo method, we study the magnetic properties of the Ashkin-
Teller model (ATM) under the effect of the crystal field with spins $S=1$ and
$\sigma=3/2$. First, we determine the most stable phases in the phase diagrams
at temperature $T=0$ using exact calculations. For higher temperatures, we use
the Monte-Carlo simulation. We have found rich phase diagrams with the ordered
phases: a Baxter $3/2$ and a Baxter $1/2$ phases in addition to a
$\left\langle\sigma S\right\rangle$ phase that does not show up either in ATM
spin 1 or in ATM spin $3/2$ and, lastly, a
$\left\langle\sigma\right\rangle=1/2$ phase with first and second order
transitions.

Key words: Ising model, Ashkin-Teller, spin-1, spin-3/2, Monte Carlo

###### Abstract

Âèêîðèñòîâóþè ìåòîä Ìîíòå Êàðëî, äîñëäæåíî ìàãíòí âëàñòèâîñò ìîäåë Àøêíà-
Òåëëåðà (ÌÀÒ) ç ñïíàìè $S$ = $1$ $\sigma$ = $3/2$ ïä äþ êðèñòàëíîãî ïîëÿ.
Ñïîàòêó âèçíàåíî íàéñòéêø ôàçè íà ôàçîâèõ äàãðàìàõ ïðè òåìïåðàòóð $T=0$ ç
äîïîìîãîþ òîíèõ îáèñëåíü. Ïðè âèùèõ òåìïåðàòóðàõ ìè âèêîðèñòîâóìî ìîäåëþâàííÿ
ìåòîäîì Ìîíòå Êàðëî. Çíàéäåíî áàãàò ôàçîâ äàãðàìè ç âïîðÿäêîâàíèìè ôàçàìè:
ôàçó Áàõòåðà $3/2$ ôàçó Áàõòåðà $1/2$ äîäàòêîâî äî ôàçè $\left\langle\sigma
S\right\rangle$, ÿêà íå ç’ÿâëÿòüñÿ í â ÌÀÒ ç ñïíîì $1$, í â ÌÀÒ ç ñïíîì $3/2$
, íàðåøò, ôàçó $\left\langle\sigma\right\rangle=1/2$ ç ïåðåõîäàìè ïåðøîãî
äðóãîãî ðîäó.

Ключов слова: ìîäåëü çíãà, Àøêí-Òåëëåð, ñïí-$1$, ñïí-$3/2$, Ìîíòå Êàðëî

##  1 Introduction

In recent years, the functioning of spins in different network structures has
been a magnetic manifestation. It also allowed one to verify the nature of the
phase transition as well as the critical behavior in the field of statistical
mechanics [1]. In addition, the properties of magnetic materials and their
technological applications such as thermomagnetic recording media and micro-
electromechanical systems [2]are characterized by the phenomenon of mixed
spins, which are well defined in the Ising model approach [3]. Studies of
magnetic materials of mixed spins have been extended to the Ising model in the
presence of a crystalline field and, specifically, are applied to the mixed
spin ($1,3/2$). The latter studies have shown some interesting behaviors using
an effective field theory. The results of the field theory study have shown
that the mixed spin system has first order transition lines as well as offers
tricritical and triple points. They also found out that the system is of the
types [4]. However, in the context of mixed spin ($1,3/2$), the Blume-Capel
Ising model was realized when a first order transition line was found
separating two ferromagnetic regions on a square cubic lattice [5]. Using
Monte Carlo simulations, it was shown that the interactions between the
nearest neighbors of the Ising model $J_{1}$, $J_{2}$ and $J_{3}$ with
frustrations are the main barriers to the transmission change in the increase
in temperature and also indicate an Ashkin-Teller behavior. This study
estimates the transition points at the critical point of Potts and confirms
the first order transition behavior in the stabilization state of $J_{3}$
antiferromagnetic [6]. In addition, they conducted studies on the nature of
the four-stage thermal phase transition degraded in a Monte Carlo simulation
and the finite size scaling. On the other hand, first-order behaviors are
noted under Potts’ critical points with four states, and thus, his work
indicates that the four-state transition in the Ising antimagnetic model
represents a similar transition [7]. In this context and to properly describe
the notion of phase transitions, Ashkin and Teller [8] developed a very
interesting model in these Ising systems and, thus, simplified the study of
statistical mechanics. In this model, one could introduce the cooperative
phenomena of quaternary alloys into a network, which is described by a
Hamiltonian in a form suitable for magnetic systems [9]. Kramers and Wannier
observed critical points of a particular case of the Ashkin-Teller model in
which three of the four components are degraded [10]. Their hypotheses
extended to the Ashkin-Teller model shown by Fan [11], and tended to be that
of Wegner who generally proved that the argument did not exist at a critical
point. Hence, it is interesting to study closely the problem of transition in
this model [12]. In addition to this, Wagner proved that the Ashkin-Teller
model was the equivalent of the alternating eight vertex model, which has not
been solved exactly; only one critical line in the phase diagram of the
isotropic Ashkin-Teller model is as precise as possible thanks to the duality
relation found by Fan [11]. One of the most interesting critical properties of
this model is the non-universality of critical behaviour on the self-doubling
lines in which the critical exponents evolve continuously [13]. The model is a
two-dimensional system in which two layers of Ising spins interact with each
other by a four-spin system. Within each of the models, there is an
interaction at two turns between the nearest neighbors, i.e., a coupling of
two Ising models that are characterized by the spins located in each cubic
network site with a four-spin interaction parameter [14]. In particular, the
numerical study of the Ashkin-Teller spin-1 model under the crystal field
effect was carried out by Badehdah et al. [15]. In addition, Wu and Lin found
out a diversity of Ising type phase transitions of the anisotropic Ashkin-
Teller model [16], and in this system, Bekhechi et al. analyzed the critical
behavior of the Ashkin-Teller model using the averaged field and the Monte
Carlo methods, with which they established that the
$\left\langle\sigma\right\rangle$ phase appears in the isotropic case when the
interactions are antiferromagnetic [14]. Recently, the development of the
Ashkin-Teller model has been achieved by the absorption of the selenium
compound on the Ni surface [17], and the phase diagram obtained from the
elastic DNA response [18]. This model also needs the study of thermodynamic
properties of superconducting capsules (CuO2 chips) [19]. Furthermore, the
oxygen estimate can also be adapted in YBa2Cu3Oz to the two-dimensional
Ashkin-Teller model [13]. This model has also been applied over the years in
other fields such as chemical reactions in metal alloys [20]. In addition and
because of the similarities to this model, which presents a complex and
important phase diagram, different theoretical and numerical methods have been
applied including Monte Carlo simulation [21], mean field approximation [22],
effective field theory [23], matrix transfer method [24] and renormalization
group theory [25]. This model can also describe phenomena [15]. In this paper,
we essentially study the isotropic spin ($S=1$, $\sigma=3/2$) of the Ashkin-
Teller model in the fundamental state at zero temperature ($T=0$) using the
Monte Carlo method. We have also examined this model at high temperatures
using Monte Carlo simulations. We also check the stable phases of this model
for different parameters $K_{4}$, $D_{1}$ and $D_{2}$. The paper is structured
as follows; after an introduction, we have determined the fundamental state of
the model and their basic diagrams. The subsequent section describes the Monte
Carlo simulation and its formalism with the presentation of the phase diagram
of the model. Then, we discuss the results of the simulations. Finally, we
present a concluding section.

##  2 Model and phase diagram of the fundamental state

In this work, we consider the Ashkin-Teller model in the case of mixed spins
$\sigma=3/2$ and $S=1$. We analyze this case under the effect of different
crystal fields. Thus, this model is described by the following Hamiltonian:

$H=-K_{2}\sum_{(\left\langle
i,j\right\rangle)}(\sigma_{i}\sigma_{j}+S_{i}S_{j})-K_{4}\sum_{(\left\langle
i,j\right\rangle)}\sigma_{i}\sigma_{j}S_{i}S_{j}-D_{1}\sum_{i}S_{i}^{2}-D_{2}\sum_{i}\sigma_{i}^{2}\,,$
(2.1)

where the variables $\sigma_{i}$ and $S_{i}$ take the values ($\pm 3/2$, $\pm
1/2$) and ($\pm 1,0$), respectively, and are located on the sites of a cubic
lattice, $\left\langle i,j\right\rangle$ refers to a pair of nearest neighbor
spins. The first term of equation (2.1) refers to the bilinear interactions
between the spins located at the sites $i$ and $j$ using the coupling
parameter $K_{2}$. Moreover, the second term refers to the interaction of the
four spins with the coupling constant $K_{4}$. The last term refers to the
existence of two ionic crystal fields $D_{1}$ and $D_{2}$. From the
contribution of a pair of $S_{1}$, $S_{2}$, $\sigma_{1}$ and $\sigma_{2}$, the
Hamiltonian is expressed as a sum of contributions of the nearest neighbors,
we obtain the pair energy as follows:

$E_{P}=-K_{2}\left[(\sigma_{1}\sigma_{2})+S_{1}S_{2}\right]-K_{4}(\sigma_{1}\sigma_{2}S_{1}S_{2})-\frac{D_{1}}{2}(S_{1}^{2}+S_{2}^{2})-\frac{D_{2}}{2}(\sigma_{1}^{2}+\sigma_{2}^{2}).$
(2.2)

According to the values containing the variables $S_{i}$ and $\sigma_{i}$, we
extract 144 $(3^{2}\times 4^{2})$ possible configurations for the ground state
at $T=0$. Using symmetry configurations, this number reduces to 24
configurations. For each set of parameters: $K_{2}$, $K_{4}$, $D_{1}$ and
$D_{2}$, we select the configuration with minimal energy $E_{p}$. This leads
to the phase diagram in the fundamental state ($T=0$). Different phases will
be given in the form ($S_{1},\sigma_{1},S_{2},\sigma_{2}$). In what follows,
we consider different situations by fixing one parameter and varying the
others (the latters will be normalized by $K_{2}$). In figure 1, we plot the
phase diagram by varying the parameter $K_{4}/K_{2}$ as a function of
$D_{2}/K_{2}$ (letting $D_{1}=0$):

Figure 1: Phase diagram of the fundamental state in the case of $D_{1}=0$.

Figure 2: Phase diagram of the fundamental state in the case $D_{1}=D_{2}=D$.

  * •

For $D_{2}/K_{2}<-1$, if $K_{4}/K_{2}<-1$ and $K_{4}/K_{2}>D_{2}/K_{2}+1/4$:
we observe that the $S_{i}$ spins are parallel such that $\left\langle
S_{i}\right\rangle=1$ and $\sigma_{i}$ spins are antiparallel. Consequently,
the stable phase obtained is the antiferromagnetic phase. For $K_{4}/K_{2}>-1$
and $K_{4}/K_{2}>-D_{2}/K_{2}-1$. We can distinguish that the spins $S_{i}$
and $\sigma_{i}$ are both aligned in the same direction, so
$\left\langle\sigma S\right\rangle=1/2$; this corresponds to the ferromagnetic
phase.

  * •

In the case: $D_{2}/K_{2}>-1$, if $K_{4}/K_{2}<-4/9$, we can observe that the
$\sigma_{i}$ spins are parallel such that
$\left\langle\sigma\right\rangle=3/2$ and $S_{i}$ spins are antiparalleled.
Consequently, the stable phase obtained is the antiferromagnetic phase.
Otherwise if $K_{4}/K_{2}>-1$ and $K_{4}/K_{2}>-D_{2}/K_{2}-1$. We can
distinguish that the spins $S_{i}$ and $\sigma_{i}$ are both aligned in the
same direction, including $\left\langle S\right\rangle=1$ and
$\left\langle\sigma\right\rangle=3/2$ so $\left\langle\sigma
S\right\rangle=3/2$ while we have the ferromagnetic phase $3/2$.

In the second situation, we obtain the figure which represents the variation
of $K_{4}/K_{2}$ as a function of $D_{1}/K_{2}$. Lastly, we put
$D_{1}=D_{2}=D$ and draw the diagram $K_{4}/K_{2}$ as a function of $D/K_{2}$
(figure 2). Figure 3 shows a phase diagram in the fundamental state and
parameter variation $K4$ as a function of crystalline field D1 ($D2=0$).

Figure 3: Phase diagram of the fundamental state in the case of $D_{2}=0$.

This figure shows a phase diagram in the fundamental state and parameter
variation $K_{4}$ as a function of crystalline field $D_{1}$ ($D_{2}=0$):

  * •

For $D_{1}/K_{2}<-0.1$ or $K_{4}/K_{2}>4/9D_{1}/K_{2}-4/9$ and
$K_{4}/K_{2}>-4/9D_{1}/K_{2}-4/9$: we have a stable phase called the phase
$\left\langle\sigma\right\rangle$ because we have the spins $S_{i}$ being
equal to zero such that $\left\langle S\right\rangle=\left\langle\sigma
S\right\rangle=0$ and $\left\langle\sigma\right\rangle=-3/2$.

  * •

For $D_{2}/K_{2}>-0.1$, if $K_{4}/K_{2}>-4/9$ and
$K_{4}/K_{2}>-D_{2}/K_{2}-1$: we observe that the $S_{i}$ spins are parallel
such that $\left\langle S\right\rangle=1$ so that the $\sigma_{i}$ spins are
antiparallel or $\left\langle\sigma\right\rangle=3/2$, consequently the stable
phase obtained is the antiferromagnetic phase. Otherwise if $K_{4}/K_{2}>-1$
and $K_{4}/K_{2}>-D_{2}/K_{2}-1$. We can observe that the spins $S_{i}$ and
$\sigma_{i}$ are both aligned in the same direction, including $\left\langle
S\right\rangle=1$ and $\left\langle\sigma\right\rangle=3/2$ so
$\left\langle\sigma S\right\rangle=3/2$ while we have the ferromagnetic phase
3/2 as the stable phase.

The third case was made using $D_{1}=D_{2}=D$ (figure 2). We draw the diagram
$K_{4}$ according to crystalline field $D$ such that:

  * •

For $K_{4}/K_{2}>-12/9D/K_{2}-12/9$ and $K_{4}/K_{2}>-0.4$; we have
$\left\langle\sigma\right\rangle=3/2$, $\left\langle S\right\rangle=1$ and
$\left\langle\sigma S\right\rangle=3/2$ so
$\left\langle\sigma\right\rangle=\left\langle\sigma S\right\rangle$ such that
the spins $\sigma_{i}$ and $S_{i}$ are both parallel, then we have the Baxter
3/2 phase called the ferromagnetic Baxter phase (the stable phase).

  * •

For $K_{4}/K_{2}<4/9D/K_{2}-4/9$ and $K_{4}/K_{2}<-0.4$; in this part of the
diagram we see that the spins $\sigma_{i}$ are parallel so that the spins
$S_{i}$ are antiparallel. This means that we have an antiferromagnetic Baxter
phase, which is always the Baxter phase (3/2).

  * •

For $K_{4}/K_{2}<-12/9D/K_{2}-12/9$ and $K_{4}/K_{2}>12/9D/K_{2}+4/9$, in this
region we have $\left\langle\sigma\right\rangle=1/2$ and $\left\langle
S\right\rangle=0=\left\langle\sigma S\right\rangle=0$, so the phase here is
the phase symbolized by $\left\langle\sigma\right\rangle$, because the spins
$S_{i}$ are equal to zero and the spins $\sigma_{i}$ while the parallel ones
designate phase 1/2.

  * •

For the zone that is noticed in the phase diagram and which is specificed by
the equations: $K_{4}/K_{2}>4/9D/K_{2}-4/9$ and $K_{4}/K_{2}<-4/9D/K_{2}-4/9$
and if $D=-1$, we have the spins $\left\langle S\right\rangle=0$ and
$\left\langle\sigma\right\rangle=3/2$, such that $S_{i}$ are equal to zero and
$\sigma_{i}$ are parallel. Finally, we obtain the phase
$\left\langle\sigma\right\rangle$ is the phase
$\left\langle\sigma\right\rangle=3/2$ which does not exist either in the case
of mixed spin $-1/2$ [14] or in the Ashkin-Teller model for spin-3/2 [26].

##  3 The Monte Carlo simulation

In our work, to determine the magnetic properties of the Ashkin-Teller model
for non-zero temperatures, we use Monte Carlo simulations implemented with the
Metropolis algorithm with periodic boundary conditions to update the lattice
configurations. We consider a 2d square lattice of $L\times L$ size which
contains $N=L^{2}$ sites. We performed the simulations for system size $L=30$.
We performed simulations for certain values of the parameters $K_{4}$, $D_{1}$
and $D_{2}$ using $P=100000$ Monte Carlo steps after discarding the first
$20000$ MCS for thermal equilibrium. The magnetization of the system is given
by the formula:

$|M_{\alpha}|=\left\langle|M_{\alpha}|\right\rangle=\frac{1}{N_{P}}\sum_{c}\sum_{i}\alpha_{i}(c)$
(3.1)

with $\alpha=\sigma,S,\sigma S$, where $i$ runs over the lattice sites and $c$
runs over the obtained system configurations obtained. The lattice is updated
by a sweep of the $N$ spins (the Monte Carlo step) after the system reaches
thermal equilibrium. The magnetic susceptibility relationship is given by:

$\chi_{\alpha}=N/(K_{\text{B}}T)\left(\left\langle
M_{\alpha}^{2}\right\rangle-\left\langle|M_{\alpha}\right\rangle^{2}\right)$
(3.2)

with $\alpha=\sigma,S,\sigma S$ and the Binder cumulant is given by:

$U_{\alpha}=1-\frac{\left\langle M_{\alpha}^{4}\right\rangle}{3\left\langle
M_{\alpha}^{2}\right\rangle^{2}}\,.$ (3.3)

Errors are deducted from the blocking method.

##  4 Results and discussions

We obtain the magnetization behavior as a function of temperature as well as
the susceptibilities of the studied system for different values of the
coupling parameters. As shown in figure 4 our MC results at low temperature
show a ferromagnetic Baxter phase $(S_{1}\sigma_{1}S_{2}\sigma_{2})$=(1 3/2 1
3/2) with $\left\langle\sigma S\right\rangle=3/2$ (figure 4a) and a
ferromagnetic Baxter phase (1 1/2 1 1/2) with $\left\langle\sigma
S\right\rangle=1/2$ (figure 4b) as expected from the $T=0$ phase diagram
(figure 1), where we find a new partially ordered phase $\left\langle\sigma
S\right\rangle$ identified by $\left\langle\sigma\right\rangle=\left\langle
S\right\rangle=0$, and $\left\langle\sigma S\right\rangle\neq 0$. For high
temperatures in both cases, the system becomes disordered. The critical
transition temperature is estimated from the maximum of the susceptibility
associated with the different magnetization. We found for the case (a)
$T_{c}=7.39$ and for case (b) $T_{c}=1.89$. In addition, the transition
between the phases mentioned is always of second order due to the continuity
of the order parameters across the transition line.

Figure 4: (Colour online) Phase diagram shows the magnetization (the
parameters $\sigma,\,\,S,\,\,\sigma S$) as a function of the temperature, with
the size $L=30$, with crystal field $D_{1}=0$; we take $K_{4}=1$ and for a)
$D_{2}=6$ and b) $D_{2}=-6$. Continuous behavior is observed signaling second
order transition.

In figure 5, the first case (a) at low temperature, we have
$\left\langle\sigma\right\rangle$ = 1/2, $\left\langle
S\right\rangle=\left\langle\sigma S\right\rangle=0$ corresponds to the phase
(0 1/2 0 1/2). However, at high temperature the system undergoes a transition
to the paramagnetic phase. For the second case (b) $D=-2$ and $K_{4}=3$, the
ground state is the ferromagnetic Baxter phase (1 3/2 1 3/2). The
susceptibility plot shows a peak corresponding to
$\left\langle\sigma\right\rangle$ and $\left\langle S\right\rangle$ at the
critical temperature $T_{c1}=11.09$; by contrast, the susceptibility
corresponding to $\left\langle\sigma S\right\rangle$ shows a distinct peak at
the transition temperature $T_{c2}=14.19$, clearly defining a partially
ordered phase $\left\langle\sigma S\right\rangle$ at high temperature
separating the disordered phase from the Baxter phase.

Figure 5: (Colour online) The magnetization (for the parameters
$\sigma,\,\,S,\,\,\sigma S$) as a function of the temperature, with system
size $L=30$, with crystal field $D_{1}=D_{2}=D$, in the two cases a)
$K_{4}=1$, $D=-6$ and b) $K_{4}=3$, $D=-2$ a partially ordered is observed at
high temperature.

The phase diagram in figure 6 shows the stable phases at different temperature
in the plane $(D_{2}/K_{2},\\\ T/K_{2})$ in the case $D_{1}=0$ for the
coupling parameter $K_{4}=1$, we found that for low values of $D_{2}/K_{2}$,
two phases are separated by a first order transition line. This is shown in
figure 8. The verification of the phase transition nature is determined from
the discontinuity or continuity of the order parameters [25].The two phases
are: ferromagnetic Baxter $1/2$ and ferromagnetic Baxter $3/2$. The former
phase $(\sigma=1/2)$ was neither found for this ATM model with spin-1/2 [9]
nor in the Ashkin-Teller model with spin-3/2 [26]. At high temperature, a
second order transition to the paramagnetic disorder phase takes place. The
Baxter ferromagnetic-3/2 is stable for large values of $D_{2}/K_{2}$.

Figure 6: (Colour online) Phase diagram in the plane $(D_{2}/K_{2},T/K_{2})$
for $K_{4}=1$ by MC simulation with $L=30$.

Figure 7: (Colour online) Phase diagram in the plane $(D/K_{2},T/K_{2})$ for
$K_{4}=1$ by MC simulation with $L=30$.

Figure 8: (Colour online) The magnetization (the parameters $\sigma$, S,
$\sigma S$) as a function of the temperature with system size $L=30$ with
$T=1$.

In figure 7 we plot the phase diagram in the $(T/K_{2},D/K_{2})$ plane. We
found out a similar form of the phase diagram as in figure 6 except that the
low temperature low $D$ phase is now $\left\langle\sigma\right\rangle=1/2$
phase with a first transition line to the ferromagnetic Baxter $3/2$ for large
values of $D$. We also note that the $\sigma-1/2$ phase was not found in the
Ashkin-Teller for spin-3/2 [26]. Moreover, when the coupling parameter values
are increased $K_{4}/K_{2}=3$ for the case of $D_{1}=D_{2}=D$ with the growth
of the values of $D/K_{2}$. In figure 9 we found out a partially ordered
$\left\langle\sigma S\right\rangle$ phase, figure 5b, between the ordered
phase of Baxter ferromagnetic 3/2 and the disordered paramagnetic phase at
high temperature. These phases were illustrated in previous cases(figure 6,
7). This new phase $\left\langle\sigma S\right\rangle$ was also found in the
model Ashkin-Teller mixed of spin-1/2 [14] but not in the Ashkin-Teller for
spin-3/2 [26]. By decreasing the crystal field $D/K_{2}$, we found a
transition from Baxter ferromagnetic $3/2$ to the paramagnetic phase which is
of first order type. We also observe the same phase
$\left\langle\sigma\right\rangle=1/2$ as in figure 7 which separates high
temperature paramagnetic phase in low values of $D/K_{2}$ by a second order
transition, and is separated with a ferromagnetic Baxter $3/2$ phase by a
first order transition.

Figure 9: In-plane phase diagram $(D/K_{2},T/K_{2})$ for $K_{4}=3$ by MC
simulation with $L=30$.

Figure 10: (Colour online) Dependance of the accumulate (a), the
susceptibility (b) and the magnetization (c) for different choices of the size
$L$ as a function of temperature in the case $K_{4}/K_{2}=1$ and
$D_{2}/K_{2}=3$ $(D_{1}=0)$ for the order parameter $\sigma$.

The phase transition points to a function of temperature at the fixed values
of the coupling $K_{4}/K_{2}$ as well as the crystalline field $D_{2}/K_{2}$,
which were pre-located from the points of intersection of the cumulative
curves of Binder specified by the equation (3.3). We show these cases for
$D_{2}/K_{2}=3$ $(D_{1}/K_{2}=0)$ for different size of
$L=(10,\,20,\,30,\,40,\,60)$ and the susceptibilities as a function of
temperature as well as the Binder cumulant as a function of temperature figure
10 plotted. It is noted in figure 10 that the peaks when $L$ increases, show
the transition point which means the critical temperature or the change of
phase transition. Nevertheless, the Binder accumulate curves whose figure 10
shows that there is an intersection point that defines the critical
temperature knowing that $T_{c}=6.55$.

##  5 Conclusion

In order to well describe the magnetic properties of Ising typical systems in
statistical mechanics, and within this framework, we analyzed the Ashkin-
Teller model with spins $(1,3/2)$ on a cubic lattice under the effect of
different crystalline fields and the coupling parameters which are defined in
equation (2.1). The first step of our study was the most important of the
stable phase in the fundamental state (zero temperature) in three cases of
crystalline field; the system undergoes a first-order phase transition between
these stable phases because we noticed that we have phases that did not exist
in ATM spin-1/2. On the other hand, when the temperature is non-zero, we have
processed the AT model by the Monte Carlo simulation, specifically, using the
Metropolis method. As a result, we found out that the coupling parameter
values were fixed and the crystal field varied with the temperature variation.
We also found out that we had the second-order phase diagram, which contained
stable phases such as the Baxter phase 3/2 as well as the paramagnetic phase
in different cases of crystalline field in the parameter space ($K_{4}/K_{2}$,
$D_{1}/K_{2}$, $D_{2}/K_{2}$, $D/K_{2}$, $T/K_{2}$) delimited by lines with
multicritical points. Crucially, we found a new phase in the phase diagram in
space ($K_{4}/K_{2}$, $D/K_{2}$, $T/K_{2}$). Finally, we verified the phase
transition nature of this model, which is of second order phase transition of
Ising type.

## References

  * [1] Ndizeye P., Hontinfinde F., Kounouhewa B., Bekhechi S., Cent. Eur. J. Phys., 2014, 12, No. 5, 375–382, doi:10.2478/s11534-014-0446-y.
  * [2] Zahir H., Bahlagui T., El Kenz A., El Bouziani M., Benyoussef A., Hasnaoui A., Sbiaai K., J. Supercond. Novel Magn., 2019, 32, 963–970, doi:10.1007/s10948-018-4768-7.
  * [3] Sidi Ahmed S., Bahmad L., El Yousfi A., Benyoussef A., El Kenz A., El Hachimi A.G. Superlattices Microstruct., 2017, 109, 841–851, doi:10.1016/j.spmi.2017.11.034.
  * [4] Sabri S., El Falaki M., El Yadaria M., Benyoussef A., El Kenz A., Physica A, 2016, 460, 210–221, doi:10.1016/j.physa.2016.04.012.
  * [5] Tucker J.W., J. Magn. Magn. Mater., 2001, 237, No. 2, 215–224, doi:10.1016/S0304-8853(01)00691-6.
  * [6] Liu R.M., Zhuo W.Z., Dong S., Lu X.B., Gao X.S., Qin M.H., Liu J.-M., Phys. Rev. E, 2016, 93, No. 3, 032114, doi:10.1103/PhysRevE.93.032114
  * [7] Liu R.M., Zhuo W.Z., Chen J., Qin M.H., Zeng M., Lu X.B., Gao X.S., Liu J.-M., Phys. Rev. E, 2017, 96, No. 1, 012103, doi:10.1103/PhysRevE.96.012103
  * [8] Ashkin J., Teller E., Phys. Rev., 1943, 64, 178, doi:10.1103/PhysRev.64.178
  * [9] Bezerra C.G., Mariz A.M., de Araujo J.M., da Costa F.A., Physica A, 2001, 292, No. 1–4, 429–436,
doi:10.1016/S0378-4371(00)00568-9.

  * [10] Kramers H.A., Wannier G.H., Phys. Rev., 1941, 60, No. 3, 252, doi:org/10.1103/PhysRev.60.252.
  * [11] Fan C., Phys. Rev. B, 1972, 6, No. 3, 902, doi:org/10.1103/PhysRevB.6.902.
  * [12] Kamieniarz G., Kozłowski P., Dekeyser R., Phys. Rev. E, 1997, 55, No. 3, 3724, doi:10.1103/PhysRevE.55.3724.
  * [13] Akinci Ü., Physica A, 2017, 469, 740–749, doi:10.1016/j.physa.2016.11.064.
  * [14] Bekhechi S., Benyoussef A., Ettaki B., Loulidi M., Eur. Phys. J. B, 2000, 18, 275–282,
doi:10.1007/s100510070058.

  * [15] Badehdah M., Bekhechi S., Benyoussef A., Ettaki B., Phys. Rev. B, 1999,  59, 6250,
doi:10.1103/PhysRevB.59.6250.

  * [16] Wu F.Y., Lin K.Y., J. Phys. C: Solid State Phys., 1974, 7, L181, doi:10.1088/0022-3719/7/9/002.
  * [17] Bak P., Kleban P., Unertl W.N., Ochab J., Akinci G., Bartelt N.C., Einstein T.L., Phys. Rev. Lett., 1985,  54, No. 14, 1539, doi:10.1103/PhysRevLett.54.1539.
  * [18] Zhe C., Ping W., Yinǵ-Honǵ Z., Commun. Theor. Phys., 2008, 49, 525, doi:10.1088/0253-6102/49/2/57.
  * [19] Grønsleth M.S., Nilssen T.B., Dahl E.K., Stiansen E.B., Varma C.M., Sudbo A., Phys. Rev. B, 2009, 79, 094506, doi:10.1103/PhysRevB.79.094506.
  * [20] Fransisco R.M., Carvalho J.J., Oyadomari G.Y., Santos J.P., Physica A, 2019, 531, 121819,
doi:10.1016/j.physa.2019.121819.

  * [21] Musiał G., Dębski L., Kamieniarz G., Phys. Rev. B, 2002, 66, 012407, doi:10.1103/PhysRevB.66.012407.
  * [22] Santos J.P., Barreto F.C.S., Physica A, 2015,  421, 316–329, doi:10.1016/j.physa.2014.11.043.
  * [23] Santos J.P., Rosa D.S., Barreto F.C.S., Phys. Lett. A, 2018, 382, 272–275, doi:10.1016/j.physleta.2017.11.029.
  * [24] Bahmed L., Benyoussef A., Ez-Zahraouy H., Phys. Status Solidi B, 2001, 226, 403–411, doi:10.1002/1521-3951(200108)226:2%3C403::AID-PSSB403%3E3.0.CO;2-3.
  * [25] De Oliveira P.M.C., Barreto F.C.S., J. Stat. Phys., 1989, 57, 53–63, doi:10.1007/BF01023634.
  * [26] Boudefla R., Bekhechi S., Hontinfinde F., Condens. Matter Phys., 2015, 18, No. 3, 33703,
doi:10.5488/CMP.18.33703.

Ukrainian 0 Âèâåííÿ ìîäåë Àøêíà-Òåëëåðà ç ñïíàìè $S$ = $1$ $\sigma$ = $3/2$ ïä
äþ ðçíèõ êðèòñàëíèõ ïîëâ ç çàñòîñóâàííÿì ìåòîäó Ìîíòå Êàðëî Ç.
Àììåð${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ñ.
Áåõåø${}^{\textsf{\footnotesize{\mbox{?}}}}$, Á.Í.
Áðàõì${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ð.
Áóäåôëà${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ã. Åç-Çàõðî¿
${}^{\textsf{\footnotesize{\mbox{?}}}}$,
A.Ðàøàä${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ Ëàáîðàòîðÿ àâòîìàòèêè, ôàêóëüòåò ïðèðîäíèèõ
íàóê, Óíâåðñèòåò ì. Òëåìñåí, BP119, Àëæèð ${}^{\textsf{\footnotesize 2}}$
Ëàáîðàòîðÿ òåîðåòèíî¿ ôçèêè, ôàêóëüòåò ïðèðîäíèèõ íàóê, Óíâåðñèòåò ì. Òëåìñåí,
BP119, Àëæèð ${}^{\textsf{\footnotesize 3}}$ Ëàáîðàòîðÿ êîíäåíñîâàíî¿ ìàòåð¿
ìæäèñöèïëíàðíèõ íàóê (LaMCScl), ôàêóëüòåò ïðèðîäíèèõ íàóê, Óíâåðñèòåò ì.
Ìîõàììåäà V, ì. Ðàáàò, Ìàðîêêî


# Mean-field Study of Charge, Spin, and Orbital Orderings
in Triangular-lattice Compounds $A$NiO2 ($A$=Na, Li, Ag)

Hiroshi Uchigaito E-mail address<EMAIL_ADDRESS>Masafumi
Udagawa  and Yukitoshi Motome Department of Applied PhysicsDepartment of
Applied Physics University of Tokyo University of Tokyo Hongo Hongo Bunkyo-ku
Bunkyo-ku Tokyo 113-8656 Tokyo 113-8656 Japan Japan

###### Abstract

We present our theoretical results on the ground states in layered triangular-
lattice compounds $A$NiO2 ($A$=Na, Li, Ag). To describe the interplay between
charge, spin, orbital, and lattice degrees of freedom in these materials, we
study a doubly-degenerate Hubbard model with electron-phonon couplings by the
Hartree-Fock approximation combined with the adiabatic approximation. In a
weakly-correlated region, we find a metallic state accompanied by
$\sqrt{3}\times\sqrt{3}$ charge ordering. On the other hand, we obtain an
insulating phase with spin-ferro and orbital-ferro ordering in a wide range
from intermediate to strong correlation. These phases share many
characteristics with the low-temperature states of AgNiO2 and NaNiO2,
respectively. The charge-ordered metallic phase is stabilized by a compromise
between Coulomb repulsions and effective attractive interactions originating
from the breathing-type electron-phonon coupling as well as the Hund’s-rule
coupling. The spin-orbital-ordered insulating phase is stabilized by the
cooperative effect of electron correlations and the Jahn-Teller coupling,
while the Hund’-rule coupling also plays a role in the competition with other
orbital-ordered phases. The results suggest a unified way of understanding a
variety of low-temperature phases in $A$NiO2. We also discuss a keen
competition among different spin-orbital-ordered phases in relation to a
puzzling behavior observed in LiNiO2.

multi-orbital Hubbard model, electron-phonon coupling, charge order, orbital
order, metal-insulator transition, triangular lattice, NaNiO2, LiNiO2, AgNiO2

##  1 Introduction

One of the most distinctive features of strongly-correlated electron systems
is diverse cooperative phenomena. [1, 2] Correlated electron systems generally
show intricate phase diagrams full of competing or coexisting states, and the
phase competition often leads to exotic many-body phenomena. Among many
aspects, there are two factors which promote such complexity, i.e., multiple
degrees of freedom [2, 3, 4] and geometrical frustration [5]. The multiple
degrees of freedom are composed of charge, spin, and orbital of electrons.
Strong electron correlations induce interplay among them, resulting in a
variety of phases. In addition, the interplay often causes exotic response to
external perturbations, such as the colossal magneto-resistance in perovskite
manganites [6]. On the other hand, the geometrical frustration promotes the
phase competition. In general, the geometrical frustration results in a huge
number of low-energy degenerate states, by suppressing conventional long-range
orders. The degeneracy yields nontrivial phenomena such as complicated
ordering, glassy behavior, and spin-liquid states. It is a long-standing
problem in condensed-matter physics to understand a variety of phenomena which
emerge from synergetic effects between the multiple degrees of freedom and the
geometrical frustration.

The family of compounds $A$NiO2 ($A$=Li, Na, Ag) is a typical example of such
geometrically-frustrated systems with multiple degrees of freedom. $A$NiO2
takes two different lattice structures depending on the cation $A$, that is,
the delafossite structure [7] (AgNiO2) and the ordered rock-salt structure [8]
(NaNiO2 and LiNiO2). Both structures are quasi-two-dimensional, composed of
stacking of Ni, O, and $A$ layers. The magnetic and transport properties are
dominated by Ni cations, which are surrounded by the octahedron of oxygens.
The NiO6 octahedra share their edges so that the Ni sites constitute the
frustrated triangular layers as shown in Fig. 1. Ni3+ cation has seven $3d$
electrons in the low-spin configuration: Six out of seven fully occupy the
lower $t_{2g}$ levels and the remaining one electron enters in the higher
$e_{g}$ levels. Hence the doubly-degenerate $e_{g}$ orbital degree of freedom
is active in these systems.

Figure 1: (Color online). Schematic picture of NiO6 layer in $A$NiO2. NiO6
octahedra share their edges with neighbors to form the triangular lattice of
Ni cations. $\gamma$ denotes the Ni-Ni bond directions. Two $e_{g}$ orbitals,
$3z^{2}-r^{2}$ and $x^{2}-y^{2}$, are shown.

$A$NiO2 shows a variety of behaviors depending on the cation $A$ in spite of
the similarity in lattice and electronic structure. NaNiO2 is a Mott insulator
with a gap of 0.24eV [9]. This compound shows a first-order structural
transition accompanied by cooperative Jahn-Teller distortions at 480K [8, 10]
and a second-order magnetic transition at 20K [10, 12, 11]. At the lowest
temperature($T$), the system exhibits orbital-ferro and A-type antiferro-spin
order (antiferromagnetic stacking of spin-ferro ordered layers)[13, 11].
LiNiO2 is also a Mott insulator with a gap of 0.2eV [14], however, it shows no
clear phase transition down to 1.4K in contrast to NaNiO2 [15]. There is
strong sample dependence in $T$ dependence of the magnetic susceptibility; the
field-cool and zero-field-cool bifurcation appears in different ways depending
on samples. This sample dependence strongly suggests the relevance of
extrinsic disorder [16]. The ground state as well as the finite-$T$ properties
remain controversial [16, 17, 18, 19, 20, 21, 22, 23, 24]. In addition to
these insulating materials, recently, a new metallic compound AgNiO2 was
synthesized [7, 25]. It shows a structural transition associated with a
$\sqrt{3}\times\sqrt{3}$ charge ordering at 365K and antiferromagnetic
transition at 20K [26, 27]. The system remains metallic down to the lowest
$T$. It was claimed that in the low-$T$ phase the system separates into rather
localized spins at 1/3 Ni sites and itinerant electrons at the remaining sites
[26, 27]. Magnetic properties at low $T$ were analyzed by considering the
competing nearest- and second-neighbor exchange couplings between localized
spins [28, 29].

So far, $A$NiO2 has been theoretically studied by the first-principle
calculations and the strong-coupling analyses. The orbital and spin ordering
observed in NaNiO2 was reproduced by the LSDA+$U$ first-principle calculations
[30, 31]. The metallicity as well as the charge ordering in AgNiO2 was also
reproduced by the first-principle calculations [25, 26]. On the other hand,
effective models in the limit of strong electron correlation, the so-called
Kugel-Khomskii models, were studied to understand the orbital and spin
ordering in NaNiO2 and the peculiar disordered state in LiNiO2 [20, 32, 33].
For AgNiO2, recently, the magnetic phase diagram was investigated by the
classical spin model with ignoring the itinerant electrons [34].

Despite the extensive studies so far, comprehensive understanding of the
ground states of $A$NiO2 has not been reached yet. Although the low-$T$ states
of NaNiO2 and AgNiO2 are reproduced by the first-principle calculations, the
mechanism of stabilizing these states is not fully clarified. In addition, the
effective-model approach does not fully succeed in reproducing the ground
states of NaNiO2 and LiNiO2 in an unified way. A possible way to explore the
comprehensive understanding is to investigate a model in a wide region of
interaction parameters systematically beyond the strong-coupling approach. In
fact, in both NaNiO2 and LiNiO2, the Mott gap is not large; the gap is
comparable to the transfer integrals. Furthermore, the newly-synthesized
AgNiO2 shows metallic behavior. These facts suggest an importance of charge
fluctuations in weakly- or intermediately-correlated regions. Electron-phonon
couplings may be another essential factors, which have not been considered
seriously in spite of the experimental facts that the structural transitions
are observed in NaNiO2 and AgNiO2.

In this study, aiming at a unified picture of these compounds $A$NiO2, we
investigate a multi-orbital Hubbard model with electron-phonon couplings on a
two-dimensional triangular lattice. Our purpose is to elucidate the
microscopic mechanism for the variety of phases in these compounds $A$NiO2. In
particular, we focus on electron-phonon couplings and charge degrees of
freedom, both of which have not been carefully examined in the previous
studies. We clarify that these two elements play an important role in the
phase competition in these complicated systems. We obtain a
$\sqrt{3}\times\sqrt{3}$ charge-ordered metallic (COM) phase in the weakly-
correlated region and an insulating phase with ferro-type spin and orbital
ordering in the intermediately- to strongly-correlated region. These two
phases reproduce many aspects of the low-$T$ states in AgNiO2 and NaNiO2,
respectively. We discuss the peculiar disordered state in LiNiO2 in relation
with a keen phase competition in the obtained phase diagram.

The organization of this paper is as follows. In Sec. 2, we describe our model
and method. We introduce the Hamiltonian term by term, and present the
approximations adopted in the calculations. We present our results in Sec. 3.
We discuss the parameter region and the mechanism to stabilize the
$\sqrt{3}\times\sqrt{3}$ COM phase and the spin-orbital-ordered insulating
phase. Finally, Sec. 4 is devoted to summary.

##  2 Model and Method

###  2.1 Model

In the present study, to elucidate the phase competition in $A$NiO2, we
investigate the ground state of the multi-orbital Hubbard model with electron-
phonon couplings. Among the five $3d$ orbitals, we consider only the twofold
degenerate $e_{g}$ orbitals, by taking account of the low-spin state of Ni3+
cations. Our Hamiltonian is written as

$\mathcal{H}=\mathcal{H}_{\rm{kin}}+\mathcal{H}_{\rm{int}}+\mathcal{H}_{\text{el-
ph}}+\mathcal{H}_{\rm{ph}},$ (1)

where $\mathcal{H}_{\rm{kin}}$, $\mathcal{H}_{\rm{int}}$,
$\mathcal{H}_{\text{el-ph}}$, and $\mathcal{H}_{\rm{ph}}$ represent the
kinetic term of electrons, the electron-electron interactions, the electron-
phonon couplings, and the elastic term of phonons, respectively. We describe
the detailed forms of each term in the following.

####  2.1.1 Kinetic term

Due to the spatial anisotropy of the $e_{g}$-orbital wave functions, transfer
integrals between Ni sites depend on the bond direction as well as orbital
types (see Fig. 1). The kinetic term in eq. (1) is written as

$\mathcal{H}_{\rm{kin}}=-\sum_{\langle
ij\rangle}\sum_{\alpha,\beta}\sum_{\sigma}t^{\gamma_{ij}}_{\alpha\beta}\bigl{(}c_{i\alpha\sigma}^{\dagger}c_{j\beta\sigma}+\rm{H.c}.\bigr{)}.$
(2)

Here, $i$ and $j$ denote the site indices, $\alpha$ and $\beta$ represent the
orbital indices with $\alpha=a(b)$ corresponding to the $3z^{2}-r^{2}$
($x^{2}-y^{2}$) orbital, $\sigma$ is the spin, and $\gamma_{ij}$ denotes the
direction of bond between the site $i$ and $j$, as shown in Fig. 1. The sum
over $\langle ij\rangle$ is taken for the nearest-neighbor sites on the
triangular lattice. The transfer integrals are given by the following
matrices;

$t^{\gamma=1}=\begin{pmatrix}t&0\\\ 0&t^{\prime}\end{pmatrix},\
t^{\gamma=2}=\begin{pmatrix}\displaystyle t_{2}&t_{3}\\\
t_{3}&t_{4}\end{pmatrix},\ t^{\gamma=3}=\begin{pmatrix}t_{2}&-t_{3}\\\
-t_{3}&t_{4}\end{pmatrix},$ (3)

for the two bases of $3z^{2}-r^{2}$ and $x^{2}-y^{2}$ orbitals. From the
symmetry of orbitals, we obtain the following relations:
$t_{2}=t/4+3t^{\prime}/4$, $t_{3}=\sqrt{3}(t-t^{\prime})/4$, and
$t_{4}=3t/4+t^{\prime}/4$, with two independent parameters, $t$ and
$t^{\prime}$. We set $t=1$ as an energy scale. The value of $t^{\prime}$
depends on both $d$-$d$ direct transfer integrals and $d$-$p$-$d$ indirect
ones in a complicated manner [32, 33]. In the following, we show the results
for $t^{\prime}=-1$ by noting that the orbital overlaps between atomic
orbitals at the neighboring sites lead to $t^{\prime}/t\sim-1$ when one
consider both contributions. An extended study in wider range of
$t^{\prime}/t$ for an effective model without phonon is found in Ref. Vernay.
The choice of $t$ and $t^{\prime}$ gives the non-interacting bandwidth $8t$.

####  2.1.2 Electron-electron interactions

Next we introduce the electron-electron interaction term
$\mathcal{H}_{\rm{int}}$ in eq. (1). We consider only the on-site Coulomb
interactions. For the doubly-degenerate $e_{g}$ orbital system,
$\mathcal{H}_{\rm{int}}$ is written as

$\displaystyle\mathcal{H}_{\rm{int}}=\mathcal{H}_{U}+\mathcal{H}_{U^{\prime}}+\mathcal{H}_{J_{\rm{H}}}+\mathcal{H}_{J_{\rm{H}}^{\prime}},$
(4)

where

$\displaystyle\mathcal{H}_{U}=$ $\displaystyle\
U\sum_{i}\sum_{\alpha}n_{i\alpha\uparrow}n_{i\alpha\downarrow},$ (5)
$\displaystyle\mathcal{H}_{U^{\prime}}=$ $\displaystyle\
U^{\prime}\sum_{i}\sum_{\sigma\sigma^{\prime}}n_{ia\sigma}n_{ib\sigma^{{}^{\prime}}},$
(6) $\displaystyle\mathcal{H}_{J_{\rm{H}}}=$ $\displaystyle\
J_{\rm{H}}\sum_{i}\sum_{\sigma\sigma^{\prime}}c^{\dagger}_{ia\sigma}c^{\dagger}_{ib\sigma^{\prime}}c_{ia\sigma^{\prime}}c_{ib\sigma},$
(7) $\displaystyle\mathcal{H}_{J_{\rm{H}}^{\prime}}=$ $\displaystyle\
J_{\rm{H}}^{\prime}\sum_{i}\sum_{\alpha\neq\alpha^{\prime}}c^{\dagger}_{i\alpha\uparrow}c^{\dagger}_{i\alpha\downarrow}c_{i\alpha^{\prime}\downarrow}c_{i\alpha^{\prime}\uparrow}.$
(8)

Here $n_{i\alpha\sigma}=c^{\dagger}_{i\alpha\sigma}c_{i\alpha\sigma}$, $U$ and
$U^{\prime}$ denote the intra- and inter-orbital Coulomb repulsions, and
$J_{\rm{H}}$ and $J_{\rm{H}}^{\prime}$ denote the exchange interaction and the
pair hopping, respectively. $\mathcal{H}_{J_{\rm{H}}}$ and
$\mathcal{H}_{J_{\rm{H}}^{\prime}}$ are called the Hund’s-rule couplings.
Hereafter we assume the relations $\ U^{\prime}=U-2J_{\rm{H}}$ and
$J_{\rm{H}}=J_{\rm{H}}^{\prime}$ to retain the rotational symmetry of the
Coulomb interaction.

####  2.1.3 Electron-phonon couplings

As to the electron-phonon couplings, we consider two relevant distortions of
NiO6 octahedra in $A$NiO2, namely, the $A_{1g}$ breathing mode and the $E_{g}$
Jahn-Teller modes. $\mathcal{H}_{\text{el-ph}}$ in eq. (1) is given by the sum
of these two contributions as

$\mathcal{H}_{\text{el-ph}}=\mathcal{H}_{\text{el-
ph}}^{\rm{br}}+\mathcal{H}_{\text{el-ph}}^{\rm{JT}}.$ (9)

The $A_{1g}$ mode corresponds to the isotropic expansion (contraction) of NiO6
octahedron [Fig. 2(a)], which lowers (raises) two $e_{g}$ energy levels
without lifting their degeneracy. Namely, the $A_{1g}$ mode couples to the
local charge on each Ni site, written in the form

$\mathcal{H}_{\text{el-ph}}^{\rm{br}}=-\gamma_{\rm{br}}\sum_{i}x_{{\rm
br},i}\left(n_{ia}+n_{ib}-1\right),$ (10)

where $n_{i\alpha}=\sum_{\sigma}n_{i\alpha\sigma}$, $x_{{\rm br},i}$ is the
amplitude of the $A_{1g}$ lattice distortion, and $\gamma_{\rm{br}}>0$ is the
corresponding coupling constant. A positive (negative) $x_{{\rm br},i}$
corresponds to an expansion (contraction).

The $E_{g}$ mode has two components, $E_{g,u}$ and $E_{g,v}$, as shown in
Figs. 2(b) and 2(c), respectively. The $E_{g,u}$ mode corresponds to the
$z$-axis elongation of NiO6 octahedron, which splits the energy levels of
$x^{2}-y^{2}$ and $3z^{2}-r^{2}$ orbitals, while the $E_{g,v}$ mode causes a
mixing of the two orbitals: The coupling to the $E_{g}$ modes is written as

$\displaystyle\mathcal{H}_{\text{el-ph}}^{\rm{JT}}=-\gamma_{\rm{JT}}$
$\displaystyle\sum_{i,\sigma}\Bigl{\\{}x_{{\rm
JT},i}\left(n_{ia\sigma}-n_{ib\sigma}\right)$
$\displaystyle\quad+\bar{x}_{{\rm
JT},i}(c_{ia\sigma}^{\dagger}c_{ib\sigma}+c_{ib\sigma}^{\dagger}c_{ia\sigma})\Bigr{\\}},$
(11)

where $x_{{\rm JT},i}$ and $\bar{x}_{{\rm JT},i}$ are the amplitudes of
$E_{g,u}$ and $E_{g,v}$ modes, respectively, and $\gamma_{\rm{JT}}$ represents
the common coupling constant.

Figure 2: Schematic pictures of the displacement of oxygens (open circles)
around a Ni cation at the origin in (a) $A_{1g}$ mode, (b) $E_{g,u}$ mode, and
(c) $E_{g,v}$ mode.

####  2.1.4 Phonon term

The phonon term $\mathcal{H}_{\rm{ph}}$ in eq. (1) consists of the on-site
term ($\mathcal{H}_{\rm{ph}}^{\rm{elastic}}$) and the inter-site term
($\mathcal{H}_{\rm{ph}}^{\rm{coop}}$) as,

$\displaystyle\mathcal{H}_{\rm{ph}}=$
$\displaystyle\mathcal{H}_{\rm{ph}}^{\rm{elastic}}+\mathcal{H}_{\rm{ph}}^{\rm{coop}}.$
(12)

Each term is given by the sum of contributions from the $A_{1g}$ and $E_{g}$
mode phonons. The first term $\mathcal{H}_{\rm{ph}}^{\rm{elastic}}$ is the
elastic energy of lattice distortions, which is given by the sum of the
following two terms;

$\displaystyle\mathcal{H}_{\rm{ph}}^{\rm{elastic,br}}=$
$\displaystyle\frac{1}{2}\sum_{i}x_{{\rm br},i}^{2},$ (13)
$\displaystyle\mathcal{H}_{\rm{ph}}^{\rm{elastic,JT}}=$
$\displaystyle\frac{1}{2}\sum_{i}\left(x_{{\rm JT},i}^{2}+\bar{x}_{{\rm
JT},i}^{2}\right).$ (14)

The elastic constants of $A_{1g}$ mode and $E_{g}$ modes are taken as unity
without losing generality, by normalizing the amplitudes of lattice
distortions ($x_{{\rm br},i}$, $x_{{\rm JT},i}$, and $\bar{x}_{{\rm JT},i}$)
and the coupling constants ($\gamma_{\rm{br}}$ and $\gamma_{\rm{JT}}$).

The second term $\mathcal{H}_{\rm{ph}}^{\rm{coop}}$ describes the cooperative
couplings of lattice distortions. The precise form of the couplings is, in
principle, determined by the phonon dispersion in the materials, but here, for
simplicity, we consider the nearest-neighbor couplings only. Then
$\mathcal{H}_{\rm{ph}}^{\rm{coop}}$ is defined by the sum of two terms;

$\displaystyle\mathcal{H}_{\rm{ph}}^{\rm{coop,br}}=$
$\displaystyle\lambda_{\rm{br}}\sum_{\langle ij\rangle}x_{{\rm br},i}x_{{\rm
br},j},$ (15) $\displaystyle\mathcal{H}_{\rm{ph}}^{\rm{coop,JT}}=$
$\displaystyle-\lambda_{\rm{JT}}\sum_{\langle ij\rangle}x_{{\rm JT},i}x_{{\rm
JT},j}-\bar{\lambda}_{\rm{JT}}\sum_{\langle ij\rangle}\bar{x}_{{\rm
JT},i}\bar{x}_{{\rm JT},j},$ (16)

with the coupling constants $\lambda_{\rm{br}}$, $\lambda_{\rm{JT}}$, and
$\bar{\lambda}_{\rm{JT}}$. Although the values of $\lambda_{\rm{JT}}$ and
$\bar{\lambda}_{\rm{JT}}$ are generally different, we take
$\lambda_{\rm{JT}}=\bar{\lambda}_{\rm{JT}}$ for simplicity. It is reasonable
to assume $\lambda_{\rm{br}}$ to be positive, i.e., the “antiferro”-type
coupling, since the $A_{1g}$-type expansion (contraction) of NiO6 octahedron
tends to shrink (expand) the neighboring octahedra due to the edge-sharing
network of octahedra. On the other hand, the tendency is opposite for the
$E_{g}$ modes; the Jahn-Teller distortion of an octahedra induces the same
distortion in the neighboring octahedra. Hence we consider the “ferro”-type
coupling $\lambda_{\rm{JT}}>0$ in the following study. In the following
calculations, the “antiferro”-type $A_{1g}$ coupling tends to stabilize a
charge ordering by differentiating charge density between the neighboring
sites, while the “ferro”-type Jahn-Teller coupling favors “ferro”-type orbital
ordering.

###  2.2 Method

In order to study the ground state of the model (1) in a wide range of
parameters, we adopt the Hartree-Fock approximation to decouple the electron-
electron interactions, and the adiabatic approximation to treat the electron-
phonon couplings. Within the Hartree-Fock approximation, the two-body
interaction terms in $\mathcal{H}_{\mathrm{int}}$ are decoupled by introducing
mean fields, $\langle
c_{i\alpha\sigma}^{\dagger}c_{i\alpha^{\prime}\sigma^{\prime}}\rangle$. The
amplitudes of lattice distortions are determined by the adiabatic
approximation. Within this approximation, the equilibrium values of $x_{{\rm
br},i}$, $x_{{\rm JT},i}$, and $\bar{x}_{{\rm JT},i}$ are determined by using
the Hellmann-Feynman theorem as

$\Big{\langle}\frac{\partial\mathcal{H}}{\partial x_{{\rm
br},i}}\Big{\rangle}=\Big{\langle}\frac{\partial\mathcal{H}}{\partial x_{{\rm
JT},i}}\Big{\rangle}=\Big{\langle}\frac{\partial\mathcal{H}}{\partial\bar{x}_{{\rm
JT},i}}\Big{\rangle}=0.$ (17)

These relations lead to the set of equations in the form

$\displaystyle x_{{\rm br},i}+\lambda_{\rm{br}}{\sum_{j}}^{\prime}x_{{\rm
br},j}=$ $\displaystyle\gamma_{\rm{br}}\left(\langle n_{ia}\rangle+\langle
n_{ib}\rangle-1\right),$ (18) $\displaystyle x_{{\rm
JT},i}-\lambda_{\rm{JT}}{\sum_{j}}^{\prime}x_{{\rm JT},j}=$
$\displaystyle\gamma_{\rm{JT}}\sum_{\sigma}\left(\langle
n_{ia\sigma}\rangle-\langle n_{ib\sigma}\rangle\right),$ (19)
$\displaystyle\bar{x}_{{\rm
JT},i}-\lambda_{\rm{JT}}{\sum_{j}}^{\prime}\bar{x}_{{\rm JT},j}=$
$\displaystyle\gamma_{\rm{JT}}\sum_{\sigma}(\langle
c_{ia\sigma}^{\dagger}c_{ib\sigma}\rangle+\langle
c_{ib\sigma}^{\dagger}c_{ia\sigma}\rangle),$ (20)

where the sum $\sum^{\prime}_{j}$ is taken over the nearest neighbors of the
site $i$.

We determine the mean fields and the lattice distortions in a self-consistent
way. For a given set of $\\{\langle
c_{i\alpha\sigma}^{\dagger}c_{i\alpha^{\prime}\sigma^{{}^{\prime}}}\rangle,x_{{\rm
br},i},x_{{\rm JT},i},\bar{x}_{{\rm JT},i}\\}$, we diagonalize the Hamiltonian
under the Hartree-Fock approximation and obtain one-particle eigenenergies and
eigenstates, which are used to calculate a new set of $\\{\langle
c_{i\alpha\sigma}^{\dagger}c_{i\alpha^{\prime}\sigma^{{}^{\prime}}}\rangle\\}$.
These new mean fields are substituted in Eqs. (18)-(20) to determine the new
set of $\\{x_{{\rm br},i}$, $x_{{\rm JT},i}$, $\bar{x}_{{\rm JT},i}\\}$. These
procedures are repeated until the convergence is reached within the precision
less than $10^{-4}$ for all the variables.

In the following calculations, we take the unit cell which includes six Ni3+
sites in the triangular lattice as shown in Fig. 3. To incorporate different
orbital orderings depending on the way of taking the $\gamma$ direction in
Fig. 1, we consider two different ways of embedding the unit cell as shown in
Fig. 3. These unit cells accommodate a $\sqrt{3}\times\sqrt{3}$ charge
ordering and a two-sublattice ordering such as a stripe-type antiferromagnetic
state. Note that the ordering patterns observed in NaNiO2 and AgNiO2 are both
included by taking the unit cells. For the initial state in the iteration, we
consider more than 30 different states with different symmetry in spin,
orbital, and charge sectors, which are relevant in the parameter space we
study. The different initial configurations are adopted for each type of the
unit cells in Fig. 3.

The integration over the wave number in the calculations of the mean fields is
replaced by the sum over $24\times 24$ grids in the Brillouin zone for the
supercell. Hereafter, we focus on the quarter-filling case, i.e., one electron
per site on average, corresponding to one $e_{g}$ electron in the low-spin
state of Ni3+.

Figure 3:  Two different ways of taking the unit cell with six Ni sites on the
triangular lattice used in the calculations. The axis $\gamma$ for the
orbitals is shown (see also Fig. 1).

##  3 Results and Discussion

In this section, we show the results obtained for the ground state of the
model (1). In particular, we focus on the roles of $U$, $J_{\rm H}/U$, and
$\gamma_{\rm{JT}}$, and discuss how these parameters affect the ground state.
We take $\gamma_{\rm{br}}=1.6$ and $\lambda_{\rm{br}}=\lambda_{\rm{JT}}=0.05$,
which result in reasonable energy gain in forming CO, as we will see below.

As a result, we find that the tendency to charge ordering becomes pronounced
by a compromise among $U$, $U^{\prime}$, $J_{\rm H}$, and $\gamma_{\rm br}$.
On-site repulsions $U$ and $U^{\prime}$ suppress charge disproportionation,
while the Hund’s-rule coupling $J_{\rm H}$ as well as the breathing-type
coupling $\gamma_{\rm br}$ works as an inter-orbital effective attractive
interaction and promotes charge disproportionation. On the other hand, for
larger $U$ and $U^{\prime}$, the system becomes insulating, and the Jahn-
Teller coupling $\gamma_{\rm{JT}}$ becomes important and enhances the tendency
to orbital ordering, concomitant with magnetic ordering. In order to
characterize the orbital-ordered phases, we introduce the pseudospin operators
in the orbital sector, defined as

$\boldsymbol{\tau}_{i}\equiv\sum_{\sigma}c^{{\dagger}}_{i\alpha\sigma}\boldsymbol{\sigma}_{\alpha\beta}c_{i\beta\sigma},$
(21)

where $\boldsymbol{\sigma}_{\alpha\beta}$ denotes the pauli matrix. For
example, $\tau_{z}$-OF means a ferro-type order of $\tau_{iz}$, which is the
$z$ component of $\boldsymbol{\tau}_{i}$.

The main result is summarized as the phase diagrams shown in Fig. 4. In the
following, we will focus on the $\sqrt{3}\times\sqrt{3}$-type charger-ordered
metallic (COM) phase found in the weak-coupling region [Fig. 4(a)], and the
orbital-ferro spin-ferro insulating ($\tau_{z}$-OF SF I) phase obtained in the
intermediately- to strongly-correlated region [Figs. 4(b) and 4(c)]. These two
phases are candidates for the low-$T$ states of AgNiO2 and NaNiO2,
respectively. We will identify the parameter range for these phases, and
discuss the origin and the nature of them in the following.

Figure 4: (Color online). Ground-state phase diagrams for the model (1) in the
plane of $J_{\rm H}/U$ and $\gamma_{\rm JT}$ at (a) $U=4$, (b) $U=8$, and (c)
$U=20$. We take $t^{\prime}=-1$, $\gamma_{\rm br}=1.6$, and $\lambda_{\rm
br}=\lambda_{\rm JT}=0.05$. The ordering patterns of each phase are shown in
Fig. 5.  Figure 5:  Schematic pictures of charge, spin, and orbital ordering
patterns in the six-site unit cell for the phases in Fig. 4. The size of the
circles schematically indicates the charge density at each Ni site. The length
of the arrows denotes the magnitude of the spin moment at each site. The open
circles denote the orbital para state, while the filled or shaded circles show
orbitally polarized states. The different patterns denote different orbitally
polarized states. See the text for details of (a)-(d). Both $\tau_{z}$-OF SF I
and $\tau_{y}$-OF SF are represented by Fig. 5(b) because these phases have
the same symmetry (OF SF). Both OF SAF I and $\tau_{x}$-OF SAF M (I) are
represented by Fig. 5(e). Typical values of of the charge disproportionation
and orbital polarization in (e)-(i) are as follows. (e) OF SAF I:
$\langle\tau_{x}\rangle\simeq\langle\tau_{z}\rangle\simeq 0.6$ at all sites.
$\tau_{x}$-OF SAF M: $\langle\tau_{x}\rangle\simeq 0.2$ and
$\langle\tau_{z}\rangle\simeq-0.06$ at all sites. (f): $\langle n\rangle\simeq
3.0$ and $\langle n\rangle\simeq 1.7$ at the two charge-rich sites and
$\langle n\rangle\simeq 0.3$ at the other charge-poor sites. (g):
$\langle\tau_{x}\rangle\simeq 0.2$, $\langle\tau_{z}\rangle\simeq 0.7$, and
$\langle n\rangle\simeq 1.2$ at the charge rich sites;
$\langle\tau_{x}\rangle\simeq 0.1$, $\langle\tau_{z}\rangle\simeq-0.2$, and
$\langle n\rangle\simeq 0.8$ at the charge poor sites. (h):
$\langle\tau_{x}\rangle\simeq 1.4$, $\langle\tau_{z}\rangle\simeq 0.6$, and
$\langle n\rangle\simeq 1.7$ at the charge rich sites and
$\langle\tau_{x}\rangle\simeq 0.2$, $\langle\tau_{z}\rangle\simeq 0.0$, and
$\langle n\rangle\simeq 0.3$ at the charge poor sites. (i): $\langle
n\rangle\simeq 1.02$ at the charge rich sites (1 and 4 in Fig. 3), $\langle
n\rangle\simeq 0.98$ at the charge poor sites (3 and 6), and $\langle
n\rangle\simeq 1.00$ at the other charge-neutral sites (2 and 5).

###  3.1  Weakly-correlated region

####  3.1.1 Phase diagram

A tendency toward charge ordering is widely observed in a weakly-correlated
region. Figure 4(a) shows the phase diagram at $U=4$. Among the competing
phases, a COM phase with $\sqrt{3}\times\sqrt{3}$-type charge ordering is
stabilized in a wide range of $\gamma_{\rm JT}$ at $J_{\rm H}/U\sim 0.25$. The
$\sqrt{3}\times\sqrt{3}$ charge ordering is a three-sublattice order, in which
one is charge rich (the density is almost two) and the other two charge poor
(the density is almost 0.5 per site). This $\sqrt{3}\times\sqrt{3}$ COM state,
shown in Fig. 5(a), is remarkable, since it has the same charge ordering
pattern as the low-$T$ state of AgNiO2 [26, 27]. The spin state is also
interesting; large moments appear at charge-rich sites ($S\simeq 0.6$), while
the moments are suppressed at charge-poor sites ($S\simeq 0.05$). A similar
differentiation was proposed for the low-$T$ state of AgNiO2 [26, 27],
although the calculated spin pattern does not fully agree with the
experimental result. We note that importance of interlayer coupling is
experimentally suggested for the magnetic ordering [28, 26, 27], which is not
taken into account in our model. We will discuss the nature of this COM phase
in Sec. 3.1.3 in detail.

Around $J_{\rm{H}}/U=0.24$ and $\gamma_{\rm{JT}}=1.1$, we find another COM
phase, i.e., the sixfold COM phase. Although this phase is stabilized in a
narrow region in the phase diagram, it remains to be metastable in a wide
parameter range, with a slightly higher energy than the ground state, as we
will discuss in the next section 3.1.2. This phase has the consistent ordering
structure with the low-$T$ state of AgNiO2 in terms of both charge and spin,
as shown in Fig. 5(c). Strictly speaking, the charge pattern of this phase has
lower symmetry compared with AgNiO2, due to the superposition of stripe-type
charge modulation onto the $\sqrt{3}\times\sqrt{3}$ charge ordering. However,
the magnitude of this modulation is very small: For charge rich sites (sites 1
and 4 in Fig. 3), we have $\langle n_{1}\rangle=1.5040$ and $\langle
n_{4}\rangle=1.5027$, while for charge poor sites (sites 2, 3, 5, and 6 in
Fig. 3), we have $\langle n_{3}\rangle=\langle n_{5}\rangle=0.7479$ and
$\langle n_{2}\rangle=\langle n_{6}\rangle=0.7487$, at $J_{\rm H}/U=0.24$ and
$\gamma_{\rm{JT}}=1.1$. The modulation gives very small charge
disproportionation of the order of $\sim 0.001$ within charge rich and poor
sites. [The small disproportionations are exaggerated in the schematic picture
in Fig. 5(c).]

In addition to the two COM states, we obtain a variety of ordered phases in
the weakly-correlated region. Among them, we focus on two phases which compete
with the COM states; the sixfold charge-ordered insulating (sixfold COI) phase
in the large $J_{\rm H}/U$ region, and the spin-ferro metallic phase with a
weak charge ordering (SF-COM) stabilized for smaller $J_{\rm H}/U$ [Fig.
4(a)]. We argue the stability of the $\sqrt{3}\times\sqrt{3}$ COM as well as
the sixfold COM in comparison with the two competing phases in the next
section.

####  3.1.2 Stability of the $\sqrt{3}\times\sqrt{3}$ COM phase

In order to clarify the competition among the $\sqrt{3}\times\sqrt{3}$ COM,
sixfold COM, sixfold COI, and SF-COM phases, we investigate the internal
energy in detail by comparing the contributions from different terms in the
Hamiltonian; $E_{\rm{kin}}\equiv\langle\mathcal{H}_{\rm{kin}}\rangle$,
$E_{U}\equiv\langle\mathcal{H}_{U}\rangle$,
$E_{U^{\prime}}\equiv\langle\mathcal{H}_{U^{\prime}}\rangle$,
$E_{J_{\rm{H}}}\equiv\langle\mathcal{H}_{J_{\rm{H}}}+\mathcal{H}_{J_{\rm{H}}^{\prime}}\rangle$,
$E_{\rm{br}}\equiv\langle\mathcal{H}_{\text{el-
ph}}^{\rm{br}}+\mathcal{H}_{\rm{ph}}^{\rm{elestic,br}}+\mathcal{H}_{\rm{ph}}^{\rm{coop,br}}\rangle$,
and the total energy $E_{\rm{tot}}\equiv\langle\mathcal{H}\rangle$. We show
the comparison as a function of $U$ at $J_{\rm H}/U=0.27$ and $\gamma_{\rm
JT}=0.5$ in Fig. 6.

For $U\lesssim 3.9$, the sixfold COI state has the lowest energy. As shown in
Fig. 5(f), this phase has a polaronic nature, namely, one site is almost fully
occupied (the local density is almost 4), and another one site accommodates
almost two electrons. At the latter site, spins of two electrons are aligned
parallel by the Hund’s-rule coupling. This phase is stabilized in a region
where the repulsive Coulomb interactions are compensated by effective
attractive interactions originating in the breathing-type electron-phonon
coupling as well as the inter-orbital Hund’s-rule coupling. In fact, it is
clearly observed in Figs. 6(d) and 6(e) that the energy gain in
$E_{J_{\rm{H}}}$ and $E_{\rm{br}}$ contributes to the stabilization of the
sixfold COI phase.

On the other hand, for $U\gtrsim 4.9$, the SF-COM state is most stabilized. As
schematically shown in Fig. 5(d), the charge ordering in this phase is a
stripe type, but the charge disproportionation is very small ($\langle
n\rangle\sim 1.03-1.09$ at charge rich sites, while $\langle n\rangle\sim
0.91-0.97$ at charge poor sites): the main feature is the spin ferromagnetic
ordering. The origin of this phase can be attributed to the Stoner mechanism
[35]. As shown in the inset of Fig. 7, the non-interacting Fermi level is
located in the vicinity of the steep peak of the density of states (DOS).
Consequently, a ferromagnetic instability is caused at a relatively small
$U\simeq 3.4$, according to the Stoner’s criterion. The characteristics of
Stoner ferromagnet are observed in the energy comparison in Fig. 6(c); $E_{U}$
becomes smallest among the competing phases.

The $\sqrt{3}\times\sqrt{3}$ COM state intervenes these two, and has the
lowest energy for $4.0\lesssim U\lesssim 4.8$. In the same parameter range,
the sixfold COM state appears as a metastable state and stays very close to
the ground state, as shown in the inset of Fig. 6(a) (the energy difference is
less than 0.02). These two COM states are stabilized by a compromise between
the different stabilization mechanisms for the sixfold COI and the SF-COM
phases. According to Fig. 6, the COM phases have higher (lower)
$E_{J_{\rm{H}}}$ and $E_{\rm{br}}$, while they have lower (higher) $E_{U}$ and
$E_{U^{\prime}}$, compared with the sixfold COI phase (the SF-COM phase).
Namely, the two COM phases are stabilized in a subtle balance between
repulsive Coulomb interactions and effective attractive interactions due to
the Hund’s-rule coupling and the breathing-type electron-phonon coupling.

Since the COM phases are stabilized in a delicate compromise, it is important
to consider their stability against the elements which are ignored in our
current analysis, such as fluctuations beyond the mean-field level and the
long-range part of electron interactions. First, we note that the Stoner
ferromagnetism is fragile when considering the electron correlation effect
beyond the mean-field approximation[36]. Therefore, the COM phases are
expected to extend to larger $U$ or smaller $J_{\rm H}/U$. Second, the
amplitudes of breathing-type distortions are fairly large in the sixfold COI
phase compared to those in the other phases. Hence this phase will be
suppressed by considering more realistic contributions from phonons, e.g.,
anharmonic terms of phonons. This may give a chance for the COM phases to
become wider also in smaller $U$ or larger $J_{\rm H}/U$ regions. Finally, the
long-range part of electron interactions generally works in favor of the
charge ordering, in particular, the $\sqrt{3}\times\sqrt{3}$ type and the
sixfold COM, as is discussed in several transition metal compounds and organic
materials [37]. Therefore, we expect that the COM phases relevant to AgNiO2
become more stable in a wider parameter range when extending the analyses
beyond the present model and method. Although the COM phases are robust in
this parameter region, the energy difference between the sixfold COM phase and
the $\sqrt{3}\times\sqrt{3}$ COM phase is very small, implying that the
magnetic ordering pattern might be affected by small perturbations, such as
inter-layer coupling. More accurate studies beyond the mean-field
approximation are necessary for fully determining the spin pattern of the
ground state.

Figure 6:  (Color online). Energy comparisons among the
$\sqrt{3}\times\sqrt{3}$ COM (circle), sixfold COM (triangle), sixfold COI
(square), and SF-COM states (diamond): $U$ dependences of (a) the total
energy, and the contribution from (b) kinetic term, (c) Coulomb repulsions
[closed (open) symbols denote $E_{U}$ ($E_{U^{\prime}}$)], (d) Hund’-rule
coupling, and (e) breathing-type electron-phonon coupling. The inset of (a)
shows the energy differences between the $\sqrt{3}\times\sqrt{3}$ COM phase
and the other competing phases. The parameters are $t^{\prime}=-1$,
$J_{\rm{H}}/U=0.27$, $\gamma_{\rm{br}}=1.6$, $\gamma_{\rm{JT}}=0.5$, and
$\lambda_{\rm{br}}=\lambda_{\rm{JT}}=0.05$.

####  3.1.3 Nature of the $\sqrt{3}\times\sqrt{3}$ COM phase

Reflecting the subtle balance between the attractive and repulsive
interactions, the $\sqrt{3}\times\sqrt{3}$ COM phase shows peculiar electronic
properties. The density of states (DOS) in the $\sqrt{3}\times\sqrt{3}$ COM
phase is shown in Fig. 7. The site- and spin-resolved DOS indicates that the
system exhibits a half-metallic nature: Up-spin electrons are localized at
charge-rich sites, showing a gap at the Fermi level, on the other hand, down-
spin electrons remain conductive, with a finite DOS at the Fermi level.
Electron correlations dominantly affect up-spin electrons; down-spin
conductive electrons preserve the non-interacting band structure. For
comparison, we show DOS for the non-interacting case in the inset of Fig. 7,
which is quite similar to DOS of down-spin conductive electrons. DOS in the
$\sqrt{3}\times\sqrt{3}$ COM phase reproduces several aspects of the result
obtained by the first-principle band calculation [26]: The electrons at charge
rich site tend to localize and the electronic structure at the other two
charge poor sites resembles each other.

This peculiar electronic state can be attributed to the delicate balance
between the attractive and repulsive interactions. We plot the effective one-
body potential in Fig. 8, which is defined as the sum of the terms in
$\mathcal{H}_{\rm{int}}$ and $\mathcal{H}_{\text{el-ph}}$, which couple to the
density operator $n_{i}$ at each site under the Hartree-Fock approximation.
Figure 8 shows that the charge-rich (-poor) sites bear attractive (repulsive)
potentials for up-spin electrons. In contrast, the cancellation between the
breathing-type electron-phonon coupling and the repulsive Coulomb interactions
leads to an almost flat potential for down-spin electrons. Consequently, these
interactions only work as a shift of chemical potential, and the down-spin
electrons retain the non-interacting band structure.

To conclude the discussions for the weakly-correlated region, the
$\sqrt{3}\times\sqrt{3}$ COM phase is stabilized by a compromise between
repulsive Coulomb interactions and attractive interactions originating from
the breathing-type electron-phonon coupling as well as the Hund’s-rule
coupling. This phase shows a half-metallic behavior with large magnetic
moments almost localized at charge-rich sites and conduction electrons moving
almost freely in the entire lattice. We note that this phase is distinguished
from the so-called pinball liquid state, in which the electrons at the charge-
rich sites exclude the conduction electrons as hard core potentials and
confine them to the honeycomb network of charge-poor sites, as discussed in a
spinless tight-binding model with intersite Coulomb repulsion on the
triangular lattice [38].

Figure 7:  (Color online). DOS per site for the $\sqrt{3}\times\sqrt{3}$ COM
phase. The Fermi level is set to be the origin of energy. As to the
$\sqrt{3}\times\sqrt{3}$ COM phase, the up-spin (down-spin) DOS is drawn on
the upper (lower) side. The parameters are chosen as $U=4$,
$J_{\rm{H}}/U=0.27$, $\gamma_{\rm{br}}=1.6$, $\gamma_{\rm{JT}}=0.5$, and
$\lambda_{\rm{br}}=\lambda_{\rm{JT}}=0.05$. The inset shows DOS for the non-
interacting case.  Figure 8:  (Color online). (a) Effective one-body
potentials at each site for up-spin and down-spin electrons (see the text for
details). The horizontal axis denotes the site indices in the unit cell, shown
in Fig. 3. The charge-rich (-poor) sites are the site 1 and 4 (2, 3, 5, and
6). The parameters are chosen as $U=4$, $J_{\rm{H}}/U=0.27$,
$\gamma_{\rm{br}}=1.6$, $\gamma_{\rm{JT}}=0.5$, and
$\lambda_{\rm{br}}=\lambda_{\rm{br}}=0.05$, consistent with Fig. 7.

###  3.2 Intermediately- to strongly-correlated region

Next, let us consider the intermediately- to strongly-correlated region.
Representative phase diagrams are shown in Figs. 4(b) and 4(c).

####  3.2.1 Phase diagram

We show the ground-state phase diagram at $U=8$ in Fig. 4(b). Note that the
value of $U$ is comparable with the non-interacting bandwidth., i.e., the
system is in the intermediately-correlated region. In this region, we find
three dominant orbital-ordered phases; $\tau_{y}$-ordered spin-ferro insulator
($\tau_{y}$-OF SF I), $\tau_{x}$-ordered spin-antiferro insulator
($\tau_{x}$-OF SAF I), and $\tau_{z}$-ordered spin-ferro insulator
($\tau_{z}$-OF SF I). The ordering pattern of $\tau_{y}$\- or $\tau_{z}$-OF SF
I ($\tau_{x}$-OF SAF I) is schematically shown in Fig. 5(b) [Fig. 5(e)]. For
the small $\gamma_{\rm JT}$ and $J_{\rm H}/U$ region, $\tau_{y}$-OF SF I is
stabilized, while it is replaced by $\tau_{x}$-OF SAF I for larger
$\gamma_{\rm JT}$ or by a charge-ordered state for larger $J_{\rm H}/U$.
Meanwhile, when both $\gamma_{\rm JT}$ and $J_{\rm H}/U$ become large,
$\tau_{z}$-OF SF I is stabilized. Among these phases, the $\tau_{z}$-OF SF I
phase deserves attention, since the spin and orbital pattern of this phase is
consistent with the low-$T$ phase of NaNiO2.

Remarkably, the three OF phases remain stable in a wide range of $U$ toward
the strongly-correlated regime. Figure 4(c) shows an example of the phase
diagram at large $U$. The result indicates that the three phases remain in
similar parameter regions of $\gamma_{\rm JT}$ and $J_{\rm H}/U$ compared to
the intermediate-$U$ case. This implies a possibility to understand the origin
of these phases from the strong-coupling analysis, i.e., by starting from the
Mott insulating state at $U=\infty$.

In fact, the $\tau_{y}$-OF SF I state was obtained for an effective spin-
orbital model in the strong-coupling limit (F3 phase of Fig. 8 in Ref.
Vernay). Our result is consistent with the previous study. Meanwhile, the
competition between the $\tau_{x}$-OF SAF I phase and the $\tau_{z}$-OF SF I
phase is obtained for the first time by explicitly taking account of electron-
phonon couplings. In the following, we will consider the mechanism of
stabilization of these phases through the detailed study of energetics.

####  3.2.2 Stability of the $\tau_{z}$-OF SF insulating phase

In order to understand the stability condition, it is instructive to rewrite
$\mathcal{H}^{\rm{JT}}_{\text{el-ph}}$, $\mathcal{H}_{J_{\rm{H}}}$, and
$\mathcal{H}_{J^{\prime}_{\rm{H}}}$ by using the pseudospin operators in eq.
(21) as

$\mathcal{H}_{\text{el-ph}}^{\rm{JT}}=-\gamma_{\rm{JT}}\sum_{i}(x_{{\rm
JT},i}\tau_{iz}+\bar{x}_{{\rm JT},i}\tau_{ix}),$ (22)
$\mathcal{H}_{J_{\rm{H}}}+\mathcal{H}_{J_{\rm{H}}^{\prime}}=\frac{J_{\rm
H}}{2}\sum\limits_{i}\\{\tau_{ix}^{2}-(n_{ia}+n_{ib})\\}.$ (23)

These equations clearly show that the Jahn-Teller coupling stabilizes the
$\tau_{x}$ and $\tau_{z}$ orbital orderings, while the Hund’s-rule coupling
destabilizes the $\tau_{x}$ orbital ordering. In Fig. 9, we compare the energy
contributions including these terms,
$E_{\rm{JT}}\equiv\langle\mathcal{H}_{\text{el-
ph}}^{\rm{JT}}+\mathcal{H}_{\rm{ph}}^{\rm{elastic,JT}}+\mathcal{H}_{\rm{ph}}^{\rm{coop,JT}}\rangle$
and $E_{J_{\rm H}}$, together with other relevant energy contributions, for
the three orbital-ordered phases.

Figure 9(b) shows that the stability of the $\tau_{y}$-OF SF I phase in the
small $\gamma_{\rm{JT}}$ region is attributed to the energy gain in the
kinetic energy $E_{\rm{kin}}$. This is consistent with the result of strong-
coupling analysis, where the kinetic energy gain through the spin-orbital
superexchange interactions is claimed to be the origin of this phase. [33] It
is also observed that all energy contributions in this phase are fairly
insensitive to $\gamma_{\rm JT}$, as shown in Fig. 9. This is expected from
the absence of coupling between $\tau_{y}$ and Jahn-Teller distortions, as is
clear from eq. (22).

On the other hand, the $\tau_{x}$-OF SAF I and $\tau_{z}$-OF SF I states lower
their energy through the coupling to the Jahn-Teller distortions [Fig. 9(e)],
as expected from eq. (22); thus they replace the $\tau_{y}$-OF SF I phase and
become the ground state for larger $\gamma_{\rm JT}$, as shown in Fig. 9(a).
The Hund’s-rule coupling plays an important role in the relative stability
between the $\tau_{x}$-OF SAF I and the $\tau_{z}$-OF SF I phases. As is
evident from eq. (23), the Hund’s-rule coupling affects only the $\tau_{x}$
ordering, and destabilizes it [Fig. 9(d)]. Furthermore, the $\tau_{z}$-OF SF I
phase is stabilized by the kinetic energy gain from the interorbital hoppings,
compared with the $\tau_{x}$-OF SAF I phase as shown in Fig. 9(b). In fact,
according to the second-order perturbation from the strong coupling limit
$U\to\infty$, the $\tau_{z}$-OF SF I phase has lower energy than $\tau_{x}$-OF
SAF I phase for $J_{\rm H}/U>\frac{8-\sqrt{10}}{18}\simeq 0.27$. The phase
boundary between the two phases in Fig. 4(c) is roughly located around this
critical value, which indicates that the phase competition is essentially
understood from the strong coupling picture.

Figure 9:  (Color online). $\gamma_{\rm JT}$ dependences of (a) the total
energy, and the contribution from (b) kinetic term, (c) Coulomb repulsions
[closed (open) symbols denote $E_{U}$ ($E_{U^{\prime}}$)], (d) Hund’s-rule
coupling, and (e) Jahn-Teller coupling. Comparison is made for the
$\tau_{z}$-OF SF I (circle), $\tau_{x}$-OF SAF I (square), and $\tau_{y}$-OF
SF I states (triangle). The parameters are $t^{\prime}=-1$, $U=20$,
$J_{\rm{H}}/U=0.25$, $\gamma_{\rm{br}}=1.6$, and
$\lambda_{\rm{br}}=\lambda_{\rm{JT}}=0.05$.

####  3.2.3 Nature of the $\tau_{z}$-OF SF insulating phase

Due to the strong electron repulsion, an excitation gap opens at the Fermi
level in the $\tau_{z}$-OF SF I phase. Hence, this state is insulating,
consistent with the low-$T$ insulating phase of NaNiO2. Figure 10 shows DOS in
the $\tau_{z}$-OF SF I phase. DOS is composed of four sectors (two in each
spin component), and the total weight of each sector is equal to one electron
per site, as indicated in the integrated DOS in the figure. This structure can
be qualitatively understood by considering the atomic limit with ignoring the
hoppings, $t$ and $t^{\prime}$. Let us assume the perfect $\tau_{z}$-OF SF
order, i.e., $\langle n_{i,3z^{2}-r^{2},\uparrow}\rangle=1$ and let other mean
fields to be zero. Then the excitation energies in the atomic limit are
estimated as

$\displaystyle E_{x^{2}-y^{2},\uparrow}$
$\displaystyle=U^{\prime}-J_{\rm{H}}-2E^{*}_{\rm{JT}},$ (24) $\displaystyle
E_{x^{2}-y^{2},\downarrow}$ $\displaystyle=U^{\prime}-2E^{*}_{\rm{JT}},$ (25)
$\displaystyle E_{3z^{2}-r^{2},\downarrow}$
$\displaystyle=U+2E^{*}_{\rm{JT}},$ (26)

where $E_{\alpha\sigma}$ signifies the energy necessary to add one electron
with orbital $\alpha$ and spin $\sigma$ to the $(3z^{2}-r^{2},\uparrow)$
ground state. $E^{*}_{\rm{JT}}$ is the energy gain from the Jahn-Teller
distortion per site, estimated as

$E^{*}_{\rm{JT}}=-\frac{\gamma_{\rm{JT}}^{2}}{2(1-6\lambda_{\rm{JT}})}.$ (27)

Substituting the parameters used in Fig. 10 ($U=20$, $U^{\prime}=10$, $J_{{\rm
H}}=5$, $\gamma_{\rm JT}=1.0$, and $\lambda_{\rm JT}=0.05$) into these
equations, we obtain $E_{x^{2}-y^{2},\uparrow}\simeq 6.4$,
$E_{x^{2}-y^{2},\downarrow}\simeq 11$, and $E_{3z^{2}-r^{2},\downarrow}\simeq
19$. These values well correspond to the mean energy of each sector of DOS
measured from that for the lowest one in Fig. 10.

Eq. (24) gives an estimate of the energy gap $\Delta$ in the atomic limit.
This atomic value is reduced for finite $t$ and $t^{\prime}$, since each
atomic level is broadened by a renormalized bandwidth $\tilde{W}$: The
estimate of energy gap is modified as $\Delta\sim
E_{x^{2}-y^{2},\uparrow}-\tilde{W}$. We use this simple estimate with
replacing $\tilde{W}$ by the bare bandwidth $W=8t$ for a comparison with the
Hartree-Fock solutions. As shown in Fig. 11, our simple estimate from the
atomic limit is qualitatively consistent with Hartree-Fock results. This fact
also supports that the electronic spectrum of the $\tau_{z}$-OF SF I phase is
adiabatically continued from the strong-coupling limit.

Figure 10:  (Color online). DOS for the $\tau_{z}$-OF SF I state. The
parameters are $t^{\prime}=-1.0$, $U=20$, $J_{\rm{H}}/U=0.25$,
$\gamma_{\rm{JT}}=1.0$, $\gamma_{\rm{br}}=1.6$, and
$\lambda_{\rm{br}}=\lambda_{\rm{JT}}=0.05$. The vertical line denotes the
Fermi level. The up-spin (down-spin) DOS is represented on the upper (lower)
side. The integrated DOS is also shown.  Figure 11:  (Color online). Energy
gap in the $\tau_{z}$-OF SF I phase as a function of $U$ at $t^{\prime}=-1.0$,
$J_{\rm{H}}/U=0.25$, $\gamma_{\rm{JT}}=1.0$, $\gamma_{\rm{br}}=1.6$, and
$\lambda_{\rm{br}}=\lambda_{\rm{br}}=0.05$. The line is the simple estimate
from the strong-coupling analysis. See the text for details.

To conclude this part, three orbital-ordered insulating phases appear
dominantly in the intermediately- to strongly-correlated region. Among them,
the $\tau_{z}$-OF SF I phase, which is relevant to NaNiO2, becomes stable in
the region where both the Jahn-Teller coupling and the Hund’-rule coupling are
substantial. This phase is understood by the strong-coupling picture under the
Jahn-Teller type electron-phonon couplings.

###  3.3 Comparison to experiments

In our calculations, we successfully reproduce the $\sqrt{3}\times\sqrt{3}$
COM phase and the $\tau_{z}$-OF SF I phase, whose ordering patterns are
consistent with the low-$T$ phases in AgNiO2 and NaNiO2, respectively. These
two phases appear in close parameter regions of $\gamma_{\rm JT}$ and $J_{\rm
H}/U$, but for different range of the on-site repulsion $U$. The
$\sqrt{3}\times\sqrt{3}$ COM phase is stabilized in the weakly-correlated
region, where $U$ is smaller than the bare bandwidth, while the $\tau_{z}$-OF
SF I phase is stabilized in the intermediately- to strongly-correlated region.

The difference in $U$ may be attributed to the structure of cation bands in
these compounds. The magnitude of the effective on-site repulsion for Ni $3d$
electrons is not solely determined by its atomic value, but it is considerably
affected by the screening effect brought about by electrons in the $A$ cation
and oxygen $p$ bands. According to the first-principle calculations, Ag bands
in AgNiO2 reside in the vicinity of the Fermi level [25, 26], while the Na
bands in NaNiO2 are located about 4eV above the Fermi level [30, 31].
Consequently, a larger screening effect is expected for AgNiO2, which reduces
the magnitude of $U$ considerably, compared to that for NaNiO2. Our results
are consistent with this trend. In the first-principle calculations [30, 31,
25, 26], the bandwidth of Ni 3$d$ bands was roughly estimated to be $2\sim
3$eV, leading to a rough estimate of $t=0.25\sim 0.4$eV. In the LSDA+$U$
calculation for NaNiO2 [30, 31], the value of $U$ was taken to be 5eV to
reproduce the correct size of excitation gap. On the other hand, a cluster-
model analysis of the photoemission spectra gave an estimate of $U$=7.0eV
[39]. These studies suggest that $U\simeq 10-30t$ is reasonable, consistent
with our results. It is noteworthy that the gap in our calculation at $U=20t$
corresponds to $0.5\sim 0.8$eV, which is in the same order of magnitude as the
experimental value $\sim 0.24$eV in NaNiO2 [9]. We also note that the CO
stabilization energy, which is estimated from the energy difference between
the $\sqrt{3}\times\sqrt{3}$ COM phase and a para phase, is $\sim$ $0.05t$ for
the present parameters $\gamma_{\rm br}$ and $\lambda_{\rm br}$: This result
leads a rough estimate that the CO stabilization energy is $\sim$ 0.01 - 0.02
eV, which is in the same order of magnitude as the CO transition temperature
observed in AgNiO2 (365K [26, 27]).

Finally, we make a brief comment on the absence of any explicit ordering in
LiNiO2. Since LiNiO2 is also a Mott insulator with a gap of 0.2eV [14], we
suppose that the compound is in the strongly-correlated region similar to
NaNiO2. In our results, there exists phase competition among three insulating
phases with different spin and orbital patterns, the $\tau_{z}$-OF SF,
$\tau_{x}$-OF SAF, and $\tau_{y}$-OF SF orderings. The competition brings
about a frustration in the spin and orbital sectors in the vicinity of the
phase boundaries. To argue the consequence of such frustration, we need to go
beyond the present mean-field analysis; however, we can expect severe
suppression of the orderings and some liquid-like or glassy behavior in the
spin-orbital coupled system. Hence, one possibility is that LiNiO2 is located
in such competing regime. It is noteworthy that the competition is brought
about by explicitly taking account of the electron-phonon couplings, which
have not been considered in the previous effective model approaches [20, 32,
33]. In addition to this intrinsic phase competition, extrinsic defects on Li
sites may play an important role in the glassy behavior. Furthermore, a long-
range strain effect might also play a role through the frustrating orbital and
lattice sectors [24]. It is interesting to take account of these factors
explicitly, by extending our model and analysis. We leave this problem for a
future study.

##  4 Summary

We have investigated the ground state of the multi-orbital Hubbard model with
electron-phonon couplings by the Hartree-Fock approximation and the adiabatic
approximation, in order to elucidate the origin of various phases observed in
$A$NiO2 in a unified way. We found the $\sqrt{3}\times\sqrt{3}$ charge-ordered
metallic phase in the weakly-correlated region and the orbital-ferro spin-
ferro ordered insulating phase in the strongly-correlated region. The
$\sqrt{3}\times\sqrt{3}$ charge-ordered metallic phase is stabilized by a
compromise between Coulomb repulsions and effective attractive interactions
from the breathing-type electron-phonon coupling as well as the Hund’s-rule
coupling. The electronic state is half metallic; up-spin electrons are
localized at the charge-rich sites, but down-spin electrons are extended and
almost free. On the other hand, the orbital-ferro spin-ferro ordered
insulating phase is stabilized by the Jahn-Teller coupling under strong
electron correlation, with a help by the Hund’s-rule coupling in the
competition with other orbital-ordered phases. These two phases are promising
candidates for the low-$T$ phases in AgNiO2 and NaNiO2, respectively. A
possible origin of the quite different electron repulsion between AgNiO2 and
NaNiO2 might be a screening effect from the cation and oxygen $p$ bands. The
puzzling glassy behavior in LiNiO2 might be ascribed to the competition among
different spin and orbital ordered states in the strongly-correlated region,
which occurs under a substantial Jahn-Teller type electron-phonon coupling.

## Acknowledgements

The authors thank M. Imada, S. Watanabe, and Y. Yamaji for fruitful
discussions. This work was supported by Grants-in-Aid for Scientific Research
(No. 17071003, 17740244, 19014020, and 19052008), Global COE Program “the
Physical Sciences Frontier”, the Next Generation Super Computing Project, and
Nanoscience Program, from MEXT, Japan.

## References

  * [1] M. Imada, A. Fujimori, and Y. Tokura: Rev. Mod. Phys. 70 (1998) 1039.
  * [2] Y. Tokura and N. Nagaosa: Science 288 (2000) 462.
  * [3] K. I. Kugel and D. I. Khomskii: Sov. Phys. JETP 37 (1973) 725.
  * [4] K. I. Kugel and D. I. Khomskii: Sov. Phys. Solid State 17 (1975) 285.
  * [5] “Frustrated Spin Systems”, edited by H. T. Diep (World Scientific, Singapore, 2005).
  * [6] “Colossal Magnetoresitive Oxides”, edited. by Y. Tokura (G & B Science Pub, 2000).
  * [7] T. Sorgel and M. Jansen: Z. Anorg. Allg. Chem. 631 (2005) 2970.
  * [8] E. Chappel, M. D. Nunez-Regueiro, G. Chouteau, O. Isnard, and C. Darie: Eur. Phys. J. B 17 (2000) 615.
  * [9] J. Molenda and A. Stoklosa: Solid State Ionics 38 (1990) 1.
  * [10] P. F. Bongers and U. Enz: Solid State Comm. 4 (1966) 153.
  * [11] C. Darie, P. Bordet, S. de Brion, M. Holzapfel, O. Isnard, A. Lechi, J. E. Lorenzo, and E. Suard: Eur. Phys. J. B 43 (2005) 159.
  * [12] P. J. Baker, T. Lancaster, S.J. Blundell, M. L. Brooks, W. Hayes, D. Prabhakaran, and F. L. Pratt: Physica B 374-375 (2006) 47-50.
  * [13] M. J. Lewis, B. D. Gaulin, L. Filion, C. Kallin, A. J. Berlinsky, H. A. Dabkowska, Y. Qiu, and J. R. D. Copley: Phys. Rev. B 72 (2005) 014408.
  * [14] R. Kanno, H. Kubo, Y. Kawamoto, T. Kamiyama, F. Izumi, Y. Takeda, M. Takano: J. Solid State Chem. 110 (1994) 216.
  * [15] K. Hirakawa, H. Kadowaki, and K. Ubukoshi: J. Phys. Soc. Jpn. 54 (1985) 3526.
  * [16] J. N. Reimers, J. R. Dahn, I. Davidson, and U. Von Sacken: J. Solid State Chem. 102 (1993) 542.
  * [17] K. Yamaura, M. Takano, A. Hirano, and R. Kanno: J. Solid State Chem. 127 (1996) 109.
  * [18] Y. Q. Li, M. Ma, D. N. Shi, and F. C. Zhang: Phys. Rev. Lett. 81 (1998) 3527.
  * [19] Y. Kitaoka, T. Kobayashi, A. Koda, H. Wakabayashi, Y. Niino, H. Yamakage, S. Taguchi, K. Amaya, K. Yamaura, M. Takano, A. Hirano, and R. Kanno: J. Phys. Soc. Jpn. 67 (1998) 3703.
  * [20] M. V. Mostovoy and D. I. Khomskii: Phys. Rev. Lett. 89 (2002) 227203.
  * [21] H. Yoshizawa, H. Mori, K. Hirota, and M. Ishikawa: J. Phys. Soc. Jpn. 59 (1990) 2631.
  * [22] J. Sugiyama, K. Mukai, Y. Ikedo, P. L. Russo, H. Nozaki, D. Andreica, A. Amato, K. Ariyoshi, and T. Ohzuku: Phys. Rev. B 78 (2008) 144412.
  * [23] A. Rougier, C. Delmas, A. V. Chadwick: Solid. State. Comm. 94 (1995) 123.
  * [24] J. -H. Chung, Th. Proffen, S. Shamoto, A. M. Ghorayeb, L. Croguennec, W. Tian, B. C. Sales, R. Jin, D. Mandrus, and T. Egami: Phys. Rev. B 71 (2005) 064410.
  * [25] T. Sorgel and M. Jansen: J. Solid State Chem. 180 (2007) 8.
  * [26] E. Wawrzynska, R. Coldea, E. M. Wheeler, I. I. Mazin, M. D. Johannes, T. Sorgel, M. Jansen, R. M. Ibberson, and P. G. Radaelli: Phys. Rev. Lett. 99 (2007) 157204.
  * [27] E. Wawrzynska , R. Coldea, E. M. Wheeler, T. Sorgel, M. Jansen, R. M. Ibberson, P. G. Radaelli, and M. M. Koza: Phys. Rev. B 77 (2008) 094439.
  * [28] E. M. Wheeler, R. Coldea, E. Wawrzynska, T. Sorgel, M. Jansen, M. M. Koza, J. Taylor, P. Adroguer, and N. Shannon: Phys. Rev. B 79 (2009) 104421.
  * [29] A. I. Coldea, A. Carrington, R. Coldea, L. Malone, A.F. Bangura, M. D. Johannes, I. I. Mazin, E.A. Yelland, J. G. Analytis, J.A.A.J. Perenboom, C. Jaudet, D. Vignolles, T. Sorgel, M. Jansen: preprint (arXiv:0908.4169).
  * [30] H. Meskine and S. Satpathy: Phys. Rev. B 72 (2005) 224423.
  * [31] H. Meskine and S. Satpathy: J. Appl. Phys. 97 (2005) 10A314.
  * [32] A.-M. Dare, R. Hayn, and J.-L. Richard: Europhys. Lett. 61 (2003) 803.
  * [33] F. Vernay, K. Penc, P. Fazekas, and F. Mila: Phys. Rev. B 70 (2004) 014428.
  * [34] L. Seabra and N. Shannon: Phys. Rev. Lett. 104 (2010) 237205.
  * [35] E. C. Stoner: Rep. Prog. Phys. 11 (1946) 43.
  * [36] J. Kanamori: Prog. Theor. Phys. 30 (1963) 275.
  * [37] H. Seo, J. Phys. Soc. Jpn. 75 (2006) 051009.
  * [38] C. Hotta and N. Furukawa: Phys. Rev. B 74 (2006) 193107.
  * [39] T. Mizokawa and A. Fujimori: Phys. Rev. B 54 (1996) 5368.