11institutetext:  1Unité de Recherche de Matière Condensée, d’Électronique et
de traitement de signal (URMACETS),
Department of physics, University of Dschang, PO Box: 67 Dschang, Cameroon.
22institutetext:
2Centre d’Études et de Recherches en Agronomie et en Biodiversité, Faculté
d’Agronomie et des Sciences Agricoles, University of Dschang, P.O. Box: 222
Dschang, Cameroon.
∗22email<EMAIL_ADDRESS>
# Magneto-electric cooling rate of a multiferroic antiferromagnetic quantum
spin system: the cumulative influence of the site-dependent magnetic and
electric fields

William Degaulle Waladi Gueagni1 Lionel Tenemeza Kenfack1,∗ Martin Tchoffo1,2
Georges Collince Fouokeng1 Lukong Cornelius Fai1

(Received: / Accepted: xxx)

###### Abstract

The magneto-electrocaloric effect which can be defined as the coupling between
magnetocaloric and electrocaloric effects attracts currently considerable
attention due to the advantages provided by the caloric effect in designing
solid-state refrigeration technologies. The magneto-electrocaloric effect of a
multiferroic antiferromagnetic spin system with the Dzyaloshinskii Moriya (DM)
interaction is investigated in this paper. The DM interaction is assimilated
to a coupling between an external site-dependent electric field and a local
electric polarization. The external magnetic is also considered as a site-
dependent magnetic field. The spin-wave theory is used as a diagonalization
method and through the canonical partition function, some thermodynamic
properties such as the Boltzmann entropy and the specific heat capacity are
obtained. Then, the adiabatic magnetic, electric, and magnetoelectric cooling
rates are also derived. The graphs obtained for the adiabatic magnetic,
electric, and magnetoelectric cooling rate show a characteristic behavior of
the caloric or multi-caloric effect which is in good agreement with some
experimental and theoretical works. Besides, the entropy response due to the
variation of the external site-dependent magnetic field exhibits two anomalous
entropy peaks indicating the existence of an intermediate phase tuneable by
the magnetic and electric site-dependent parameters. Overall, it is
demonstrated that the cumulative influence of the site-dependent magnetic and
electric fields allows us not only to reveal quantum critical points hidden in
a multiferroic quantum spin system but also to control the caloric or multi-
caloric effect essential in the construction of solid-state refrigeration
devices.

###### Keywords:

magneto-electrocaloric effect adiabatic magnetic/electric/magnetoelectric
cooling ratemultiferroic antiferromagnetic entropy specific heat capacity
site-dependent magnetic/electric field

##  1 Introduction

The exposition of a given material to a suitable varying external field
instigates its heating or its cooling. This is the caloric effect. Thus, the
magnetocaloric effect, the electrocaloric, and elastocaloric effect occurs
when a material is exposed to a suitable magnetic field, electric field, and
mechanical elastic field respectively. The multiferroic materials which
exhibit two or more ferroics orders offer the possibility to obtain multi-
caloric effects 1 ; 2 ; 3 . Note that the ferroic material is a material which
chooses spontaneously a preferable switchable alignment (anti/ferromagnet is
witched by an external magnetic field, ferroelectric material is switched by
an external electric field, etc.). The multiferroic materials are then found
to be suitable for the technologies made by multi-caloric effects. The
magneto-electrocaloric effect is an example of multi-caloric effects resulting
from the coupling between magnetocaloric and electrocaloric effects. The
caloric phenomena have been already observed by Joule in the mid of 19th
century as the thermo-elastic properties result of solids 4 but the
magnetocaloric effect terminology appeared the first times in the literature
in the mid of 20th century from the works of Weiss and Picard on nickel (Ni) 5
. In recent decades and to date, the caloric effect received much effort from
both experimental and theoretician physicists due to their promising
application in cooling device technologies 6 ; 7 ; 8 ; 9 ; 10 . The
magnetocaloric effect which is the most caloric property studied today, allows
the creation of magnetic refrigeration recognized to be a viable alternative
to conventional vapor-compression refrigeration technology (see the review in
ref 11 and references therein for more detail). There are several
investigations done concerning the electrocaloric effect in the literature 12
; 13 ; 14 ; 15 ; 16 ; 17 ; 18 but so far, there are still few works done for
the magneto-electrocaloric effect 19 ; 20 ; 21 .

Moreover, the magneto-electrocaloric effect implies the suggestion of the
system under the cumulative influence of the magnetic and electric fields. The
coexistence of ferroelectricity and magnetism is more explained by the
Dzyaloshinskii Moriya (DM) interaction (the combination of the exchange
interaction and spin-orbit coupling) which plays a crucial role in the field
of study of multiferroic materials. Two theories describe the mechanism of
magnetism-driven ferroelectricity in the sense of the DM interaction 22 , the
theory of Sergienko and Dagotto ( called Sergienko-Dagotto model) based on the
ionic displacement 23 and the theory proposed by Katsura, Nagaosa, and
Balatsky (the so-called KNB theory) based on the electronic bias 24 . The main
prediction of the KNB theory exploited in this paper is the expression of
electric polarization. Based on the definition of the DM interaction and the
electric polarization, we have already demonstrated that the DM interaction
intrinsically induces ferroelectricity and then ensures the multiferroicity of
the system 25 . That idea is better understood with the works of Pradeep
Thakur and P. Durganandini 26 and Jozef Sznajd 27 , in which based on the KNB
theory they interpreted the Dzyaloshinskii Moriya interaction as the coupling
of an external electric field and a local electric polarization. Indeed, the
antisymmetric DM interaction, which is the unique spin-spin interaction
expressed by using a cross product is defined as
$H_{DM}=D\sum_{i}\hat{S}_{i}\times\hat{S}_{i+1}$ where $D$ is the strength of
DM interaction and $\hat{S}_{i}$ the spin operator of the atom located on-site
$i$. When the strength $D$ is taken in the z-direction, the DM interaction can
be assimilated to the coupling of an electric field taken in the y-direction
and a local polarization
$P\propto\sum_{i}e\times\hat{S}_{i}\times\hat{S}_{i+1}$ by taking the unit
vector$e$ along the x-axis and pointing from the site $i$ to the site $i+1$.
That is the best way to take into account the interplay between magnetism and
electricity in our system with DM interaction.

Furthermore, it is important to note that there are two kinds of DM
interactions such as the uniform and the site-dependent DM interaction also
called staggered or modulated DM interaction. A lot of systems have been
studied with uniform DM interaction (see 28 and references therein). For
instance, Shuling Wang et al. investigated the one-dimensional charge transfer
magnets and found that the staggered DM interaction allow to control the
transition temperature and the physical properties (electric polarization and
magnetic properties) of the such system 29 , Avalishvili et al. studied the
ground state phase diagram of a spin XXZ Heisenberg chain with modulated DM
interaction and theirs results reveal the formation of four states with two
types of phase transitions such as Benezinskii-Kosterlitz-Thoules type and
Ising type 30 , the study of the effect of the staggered DM interaction on a
quasi-two-dimensional Shastry-Sutherland lattice by Tianqi Chen et al. show
that the such interaction significantly modify the magnetization plateaux 31 .
The effect of staggered DM interaction allows us to explore the phase diagram
and to tune the phase transitions phenomena in multiferroics.

Indeed, the influence of both site-dependent magnetic and electric fields on
the quantum phase transition occurring in multiferroic antiferromagnets is one
of the main interests of this work. The caloric effect which is intimately
related to the quantum critical points appears to be the best issue to address
that investigation. One of the quantifiers of the caloric effect is the
adiabatic cooling rate, a physical parameter very helpful in the study of
quantum phase transitions 32 . The adiabatic cooling rate, defined as an
adiabatic temperature change of a system in response to a suitable applied
external field is closely related to the generalized Grüneisen ratio 33 .
Thus, the adiabatic magnetic and electric cooling rate interpreted as a
particular case of the generalized Grüneisen ratio change theirs sign when the
magnetic and electric fields (which are the parameters governing the zero-
temperature quantum phase transition) cross their critical values.

In this work, both the site-dependent DM interaction and the site-dependent
magnetic field is taking into account. Since the DM interaction is assimilated
to the coupling of electric field and electric polarization, the electric
field will be staggered. Throughout this work, the terminology “site-dependent
field” will be used to denote the “staggered or random field” because the
fields considered are closely related to the site parity of the atoms. Notice
that there exist in the literature systems dealing with the site-dependent
magnetic field 34 ; 35 . So far, to the best of our knowledge, the magneto-
electrocaloric effect in a multiferroic antiferromagnet has not yet been
investigated with both the site-dependent magnetic field and site-dependent
electric field.

Today, the multiferroic antiferromagnetic system attracts considerable
attention because it is susceptible to offer the new opportunity in multi-
caloric technologies 36 ; 37 ; 38 ; 39 , as also demonstrated in this work.
Indeed, because of their multifunctional physical properties, robustness
against external perturbation and very fast operation, multiferroic
antiferromagnetic system tends to replace ferromagnetic materials 40 .

The main aim of this work is to investigate in detail the cumulative influence
of the site-dependent magnetic field and site-dependent electric field on the
magneto-electrocaloric effect in a multiferroic antiferromagnetic quantum spin
system. The quantum phase transition phenomena are also addressed.

The rest of the paper is structured as follows, the Hamiltonian and spin-wave
theory are presented in Sec. II, the entropy and specific heat capacity
properties are derived in Sec. III, the adiabatic cooling rates are studied in
Sec. IV and the paper is closed in Sec. V with a conclusion.

##  2 The Hamiltonian and spin-wave theory

The model considered for a multiferroic antiferromagnetic is an anisotropic
Heisenberg model with Dzyaloshinskii Moriya interaction and on-site easy-axis
anisotropic interaction 41 ; 42 as shown in the previous work 25 . More
precisely, it is an antiferromagnetic lattice divided into two
interpenetrating sublattices a and b, and the full Hamiltonian is given as:

$\begin{split}&H=-g\mu_{B}\sum_{i}\left[B_{0}\left(1+(-1)^{i}\lambda\right)+B_{A}\right]\hat{S}_{a,i}^{z}-g\mu_{B}\sum_{j}\left[B_{0}\left(1+(-1)^{j}\lambda\right)-B_{A}\right]\hat{S}_{b,j}^{z}\\\
&+\sum_{i,\delta}J\left(\hat{S}_{a,i}^{x}\hat{S}_{b,i+\delta}^{x}+\hat{S}_{a,i}^{y}\hat{S}_{b,i+\delta}^{y}+\Delta\hat{S}_{a,i}^{z}\hat{S}_{b,i+\delta}^{z}\right)+\sum_{j,\delta}J\left(\hat{S}_{b,j}^{x}\hat{S}_{a,j+\delta}^{x}+\hat{S}_{b,j}^{y}\hat{S}_{a,j+\delta}^{y}+\Delta\hat{S}_{b,j}^{z}\hat{S}_{a,j+\delta}^{z}\right)\\\
&+\left[\sum_{i,\delta}E_{0}\left(1+(-1)^{i}\tau\right)\left(\hat{S}_{a,i}^{x}\hat{S}_{b,i+\delta}^{y}-\hat{S}_{a,i}^{y}\hat{S}_{b,i+\delta}^{x}\right)-\sum_{j,\delta}E_{0}\left(1+(-1)^{j}\tau\right)\left(\hat{S}_{b,j}^{x}\hat{S}_{a,j+\delta}^{y}-\hat{S}_{b,j}^{y}\hat{S}_{a,j+\delta}^{x}\right)\right]\\\
&-D^{\prime}\left[\sum_{i}\left(\hat{S}_{a,i}^{z}\right)^{2}+\sum_{j}\left(\hat{S}_{b,j}^{z}\right)^{2}\right]\end{split}$
(1)

where $g$ is the gyromagnetic factor, $\mu_{B}$ is the Bohr magneton, $J$ is
the exchange interaction, $\Delta$ the Heisenberg anisotropy parameter,
$\hat{S}_{a,i}$ and $\hat{S}_{b,j}$ represent the spin operators of the
$ith\;(jth)$ atom on sublattice $a$ and $b$, respectively, with each
sublattice containing $N$ atoms. The indices $i$ and $j$ label the $N$ atoms
in sublattice $a$ and $b$, respectively, whereas the vector $\delta$ connects
atom $i(j)$ with its nearest neighbors. $B_{0}\left(1+(-1)^{i}\lambda\right)$
is the applied external site-dependent magnetic field in the z-direction.
$B_{A}$ is the anisotropy field, assumed to be positive which approximates the
effect of the crystal anisotropic energy with the property of turning for a
positive magnetic moment $\mu_{B}$ to align the spin on sublattice $a$ in the
positive z-direction and spins on sublattice $b$ in the negative z-direction.
$E_{0}\left(1+(-1)^{i}\tau\right)$ is the strength of the z-component of the
staggered DM interaction assimilated in this work to the external site-
dependent electric field in the y-direction coupling to the local electric
polarization by taking the spin chain along the x-direction. The case with
$B_{0}=0\left(E_{0}=0\right)$ corresponds to the absence of the external
magnetic field (electric field) while the case with
$\lambda=0\left(\tau=0\right)$ corresponds to the uniform magnetic field
(electric field). So $\lambda$ and $\tau$ are the magnetic site-dependent
parameter and the electric site-dependent parameter respectively. $D^{\prime}$
is the on-site easy-axis anisotropy parameter for anisotropy interaction also
looking as potential energy term, assumed to be positive. Also, for sake of
simplicity, we have considered only the nearest neighbor interaction.

The Holstein-Primakoff transformations allow us to express the spin operators
of the system to bosonic operators,

$\hat{S}_{a,i}^{+}=\sqrt{2S-a_{i}^{+}a_{i}}a_{i},\,\,\,\hat{S}_{a,i}^{-}=a_{i}^{+}\sqrt{2S-a_{i}^{+}a_{i}},\,\,\,\hat{S}_{a,i}^{z}=S-a_{i}^{+}a_{i},$
(2)
$\hat{S}_{b,j}^{+}=b_{j}^{+}\sqrt{2S-b_{i}^{+}b_{i}},\,\,\,\hat{S}_{b,j}^{-}=\sqrt{2S-b_{i}^{+}b_{i}}b_{j},\,\,\,\hat{S}_{b,j}^{z}=b_{j}^{+}b_{j}-S,$
(3)

then consider the situation that our system is in the low-temperature and low-
excitation limit R11 , such that the thermal averages
$\mathrm{<}a_{i}^{+}a_{i}\mathrm{>}$ and $\mathrm{<}b_{j}^{+}b_{j}\mathrm{>}$
are very small compared to $2S$. By considering only the first term of the
expansion in power series of the square root part of Eqs. ((2) and (3)) in
terms of $a_{i}^{+}a_{i}/2S$ and $b_{j}^{+}b_{j}/2S$,

$\begin{split}&\left(1-\frac{1}{2S}a_{i}^{+}a_{i}\right)^{\frac{1}{2}}=1-\frac{1}{4S}a_{i}^{+}a_{i}+\ldots,\\\
&\left(1-\frac{1}{2S}b_{j}^{+}b_{j}\right)^{\frac{1}{2}}=1-\frac{1}{4S}b_{j}^{+}b_{j}+\ldots,\end{split}$

the spin ladder operators can be approximated as
$S_{a,i}^{+}\approx\sqrt{2S}a_{i}$ and $S_{b,j}^{+}\approx\sqrt{2S}b_{j}^{+}$.
Note that in this approximation the products containing more than three
bosonic operators are neglected. Thus, we neglect the product of four bosonic
operators and denote by $M$ the numbers of nearest neighbor, the Hamiltonian
in Eq. (1) is rewritten in spin-wave approximation R11 as,

$\begin{split}&H=H_{0}+\left[g\mu_{B}(B_{+}+B_{A})+2SD^{\prime}+J\Delta
MS\right]\sum_{i}a_{2i}^{+}a_{2i}+\left[-g\mu_{B}(B_{+}-B_{A})+2SD^{\prime}+J\Delta
MS\right]\sum_{j}b_{2j}^{+}b_{2j}\\\
&+\left[g\mu_{B}\left(B_{-}+B_{A}\right)+2SD^{\prime}+J\Delta
MS\right]\sum_{i}a_{2i+1}^{+}a_{2i+1}+\left[-g\mu_{B}\left(B_{-}-B_{A}\right)+2SD^{\prime}+J\Delta
MS\right]\sum_{j}b_{2j+1}^{+}b_{2j+1}\\\ &+J\Delta
S\left(\sum_{i,\delta}a_{2i+\delta}^{+}a_{2i+\delta}+\sum_{j,\delta}b_{2j+\delta}^{+}b_{2j+\delta}\right)+J\Delta
S\left(\sum_{i,\delta}a_{2i+1+\delta}^{+}a_{2i+1+\delta}+\sum_{j,\delta}b_{2j+1+\delta}^{+}b_{2j+1+\delta}\right)\\\
&+SJ\left[\sum_{i,\delta}\left(a_{2i}b_{2i+\delta}+a_{2i}^{+}b_{2i+\delta}^{+}\right)+\sum_{j,\delta}\left(b_{2j}^{+}a_{2j+\delta}^{+}+b_{2j}a_{2j+\delta}\right)\right]\\\
&+SJ\left[\sum_{i,\delta}\left(a_{2i+1}b_{2i+1+\delta}+a_{2i+1}^{+}b_{2i+1+\delta}^{+}\right)+\sum_{j,\delta}\left(b_{2j+1}^{+}a_{2j+1+\delta}^{+}+b_{2j+1}a_{2j+1+\delta}\right)\right]\\\
&-iSE_{+}\left[\sum_{i,\delta}\left(a_{2i}^{+}b_{2i+\delta}^{+}-a_{2i}b_{2i+\delta}\right)-\sum_{j,\delta}\left(b_{2j}a_{2j+\delta}-b_{2j}^{+}a_{2j+\delta}^{+}\right)\right]\\\
&-iSE_{-}\left[\sum_{i,\delta}\left(a_{2i+1}^{+}b_{2i+1+\delta}^{+}-a_{2i+1}b_{2i+1+\delta}\right)-\sum_{j,\delta}\left(b_{2j+1}a_{2j+1+\delta}-b_{2j+1}^{+}a_{2j+1+\delta}^{+}\right)\right],\end{split}$
(4)

where $H_{0}=-2g\mu_{B}SNB_{A}-2N(MJ\Delta+D^{\prime})S^{2}$,
$B_{\pm}=B_{0}(1\pm\lambda)$ and $E_{\pm}=E_{0}(1\pm\tau)$ are the site-
dependent magnetic field and electric field respectively, in which
$B_{+}(E_{+})$ and $B_{-}(E_{-})$ denote the site-dependent magnetic field
(electric field) for the even sites $2i$ or $2j$ and for the odd sites $2i+1$
or $2j+1$ respectively. It is worth mentioning that the site-dependent field
can be modeled by using the geometric configuration. The staggered DM
interaction have been introduced in the XY model through the zigzag geometry
see Ref.R12 .

Using the inverse of the Fourier transformation 43 of the bosonic operators
for even and odd atom sites of each sublattice given by

$\begin{split}&a_{l:even}=\sqrt{\frac{2}{N}}\sum_{k}e^{-ik.l}a_{1k},\,\,\,a_{l:odd}=\sqrt{\frac{2}{N}}\sum_{k}e^{-ik.l}a_{2k}\\\
&b_{l:even}=\sqrt{\frac{2}{N}}\sum_{k}e^{ik.l}b_{1k},\,\,\,b_{l:odd}=\sqrt{\frac{2}{N}}\sum_{k}e^{ik.l}b_{2k}\end{split}$
(5)

where the wave vector $k$ in each expression takes $N/2$ values and belong to
the reduced Brillouin zone
$BZ^{\prime}=\left[-\frac{\pi}{2};\frac{\pi}{2}\right]$, the Hamiltonian given
in Eq. (4) can be rewritten in momentum space and under the complex quadratic
form as

$\begin{split}&H=H_{0}+\sum_{k}\varepsilon\left(a_{1k}^{+}a_{1k}+b_{1k}^{+}b_{1k}\right)+\sum_{k}\left(\lambda_{1k}^{*}a_{1k}b_{1k}+\lambda_{1k}a_{1k}^{+}b_{1k}^{+}\right)+g\mu_{B}B_{+}\sum_{k}\left(a_{1k}^{+}a_{1k}-b_{1k}^{+}b_{1k}\right)\\\
&+\sum_{k}\varepsilon\left(a_{2k}^{+}a_{2k}+b_{2k}^{+}b_{2k}\right)+\sum_{k}\left(\lambda_{2k}^{*}a_{2k}b_{2k}+\lambda_{2k}a_{2k}^{+}b_{2k}^{+}\right)+g\mu_{B}B_{-}\sum_{k}\left(a_{2k}^{+}a_{2k}-b_{2k}^{+}b_{2k}\right)\end{split}$
(6)

where $\gamma_{k}=\frac{1}{M}\sum_{k}e^{ik\delta}$ is the Fourier transform
coupling constant, $\varepsilon=g\mu_{B}B_{A}+2SD^{\prime}+2\Delta JMS$,
$\lambda_{1k}=\gamma_{k}\left(2JMS-2E_{+}MSi\right)$,
$\lambda_{2k}=\gamma_{k}\left(2JMS-2E_{-}MSi\right)$ ,$\lambda_{1k}^{*}$ and
$\lambda_{2k}^{*}$ are the complex conjugates of $\lambda_{1k}$ and
$\lambda_{2k}$ respectively.

From the complex quadratic form, the matrix form of the Hamiltonian is
obtained as,

$H=H_{0}+\sum_{k}\left(\frac{1}{2}X^{+}TX-2\varepsilon\right)-g\mu_{B}\sum_{k}\frac{1}{2}X^{+}T^{\prime}X$
(7)

where

$X=\left(X^{+}\right)^{+}=\left(\begin{array}[]{c}{\begin{array}[]{c}{a_{1k}}\\\
{b_{1k}^{+}}\end{array}}\\\ {\begin{array}[]{c}{b_{1k}}\\\
{a_{1k}^{+}}\end{array}}\\\ {\begin{array}[]{c}{a_{2k}}\\\
{b_{2k}^{+}}\end{array}}\\\ {\begin{array}[]{c}{b_{2k}}\\\
{a_{2k}^{+}}\end{array}}\end{array}\right),\,\,\,T=\left(\begin{array}[]{cccc}{\begin{array}[]{cc}{\varepsilon}&{\lambda_{1k}}\\\
{\lambda_{1k}^{*}}&{\varepsilon}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{\varepsilon}&{\lambda_{1k}}\\\
{\lambda_{1k}^{*}}&{\varepsilon}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{\varepsilon}&{\lambda_{2k}}\\\
{\lambda_{2k}^{*}}&{\varepsilon}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}\\\ {\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{\varepsilon}&{\lambda_{2k}}\\\
{\lambda_{2k}^{*}}&{\varepsilon}\end{array}}\end{array}\right),$

and

$T^{\prime}=\left(\begin{array}[]{cccc}{\begin{array}[]{cc}{-B_{+}}&{0}\\\
{0}&{B_{+}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{B_{+}}&{0}\\\
{0}&{-B_{+}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{-B_{-}}&{0}\\\
{0}&{B_{-}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{B_{-}}&{0}\\\
{0}&{-B_{-}}\end{array}}\end{array}\right).$

Finally, the Hamiltonian in Eq. (7) can be diagonalized by following the

Bogoliubov transformations 44 given by,

$\left\\{\begin{split}&a_{mk}=u_{mk}\alpha_{mk}+v_{mk}\beta_{mk}^{+}\\\
&b_{mk}=v_{mk}\beta_{mk}+u_{mk}\alpha_{mk}^{+}\end{split}\right.$ (8)

where the coefficients $u_{mk}$ and $v_{mk}$ are complex numbers and
throughout this work $m=\left\\{1;2\right\\}$ are used to differ the
parameters denoting the even atom sites from those denoting the odd atoms
sites. It is important to note that the new operators $\alpha_{mk}$ and
$\beta_{mk}$ obey the boson commutation rules, thus, leading to the following
constraint $\left|u_{mk}\right|^{2}-\left|v_{mk}\right|^{2}=1$. Therefore,
$u_{mk}$ and$v_{mk}$ can be parameterized as $u_{mk}=\cosh\theta_{mk}$ and
$v_{mk}=e^{i\varphi_{m}}\sinh\theta_{mk}$ with
$\varphi_{m}=\varphi^{v_{mk}}-\varphi^{u_{mk}}$. From the above
transformations, the vectors $X$ and $X^{+}$ are given by
$X=\vartheta\cdot\phi$ and $X^{+}=\phi^{+}\cdot\vartheta^{+}$ where
$\phi=\left(\phi^{+}\right)^{+}=\left(\begin{array}[]{c}{\begin{array}[]{c}{\alpha_{1k}}\\\
{\beta_{1k}^{+}}\end{array}}\\\ {\begin{array}[]{c}{\beta_{1k}}\\\
{\alpha_{1k}^{+}}\end{array}}\\\ {\begin{array}[]{c}{\alpha_{2k}}\\\
{\beta_{2k}^{+}}\end{array}}\\\ {\begin{array}[]{c}{\beta_{2k}}\\\
{\alpha_{2k}^{+}}\end{array}}\end{array}\right)$ and

$\vartheta=\left(\vartheta^{+}\right)^{+}=\left(\begin{array}[]{cccc}{\begin{array}[]{cc}{u_{1k}}&{v_{1k}}\\\
{v_{1k}^{*}}&{u_{1k}^{*}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{u_{1k}}&{v_{1k}}\\\
{v_{1k}^{*}}&{u_{1k}^{*}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{u_{2k}}&{v_{2k}}\\\
{v_{2k}^{*}}&{u_{2k}^{*}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}\\\ {\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{u_{2k}}&{v_{2k}}\\\
{v_{2k}^{*}}&{u_{2k}^{*}}\end{array}}\end{array}\right).$

Substituting these vectors in Eq. (7) and once the calculations are performed,
the resulting Hamiltonian is written as

$H=H_{0}+\sum_{k}\left(\frac{1}{2}\phi^{+}D\phi-2\varepsilon\right),$ (9)

Where D is a diagonal matrix

$D=\left(\begin{array}[]{cccc}{\begin{array}[]{cc}{\omega_{1k}^{\left(+\right)}}&{0}\\\
{0}&{\omega_{1k}^{\left(-\right)}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{\omega_{1k}^{\left(-\right)}}&{0}\\\
{0}&{\omega_{1k}^{\left(+\right)}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\ {0}&{0}\end{array}}\\\
{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{\omega_{2k}^{\left(+\right)}}&{0}\\\
{0}&{\omega_{2k}^{\left(-\right)}}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}\\\ {\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{0}&{0}\\\
{0}&{0}\end{array}}&{\begin{array}[]{cc}{\omega_{2k}^{\left(-\right)}}&{0}\\\
{0}&{\omega_{2k}^{\left(+\right)}}\end{array}}\end{array}\right).$

Thus, we obtain the so-called diagonalized form for the Hamiltonian as

$\begin{split}&H=H^{\prime}_{0}+\sum_{k}\left[\omega_{1k}^{\left(+\right)}\left(\alpha_{1k}^{+}\alpha_{1k}+\frac{1}{2}\right)+\omega_{1k}^{\left(-\right)}\left(\beta_{1k}^{+}\beta_{1k}+\frac{1}{2}\right)\right]\\\
&+\sum_{k}\left[\omega_{2k}^{\left(+\right)}\left(\alpha_{2k}^{+}\alpha_{2k}+\frac{1}{2}\right)+\omega_{2k}^{\left(-\right)}\left(\beta_{2k}^{+}\beta_{2k}+\frac{1}{2}\right)\right]\end{split}$
(10)

where $H^{\prime}_{0}=H_{0}-2\varepsilon N$. From the above Hamiltonian, each
frequency of the magnon at the symmetric positions obtained in the case of
uniform fields 25 is divided into two pairs of frequencies in this case with
site-dependent magnetic and electric field, which can be easily derived as

$\omega_{1k}^{\left(\pm\right)}=\tilde{\varepsilon}_{1}\pm
g\mu_{B}B_{+},\,\,\,\omega_{2k}^{\left(\pm\right)}=\tilde{\varepsilon}_{2}\pm
g\mu_{B}B_{-}$ (11)

with
$\tilde{\varepsilon}_{m}=\sqrt{\varepsilon^{2}-\left|\lambda_{m}\right|^{2}}$.
It is important to remark that the equilibrium state of the system is ensured
by the following constraint $\varepsilon>\left|\lambda_{m}\right|$. In
contrast to the case of uniform magnetic field where only one equilibrium
critical point has been observed 25 , here there are two equilibrium critical
points corresponding to the two pairs of frequencies and which can be obtained
by taking $\omega_{mk}^{\left(-\right)}=0$ and $k=0$,

$B_{c1}=\frac{2JMS}{g\mu_{B}}\sqrt{\left(\frac{D^{\prime}}{JM}+\Delta+\frac{g\mu_{B}B_{A}}{2JMS}\right)^{2}-\left(\frac{E_{+}^{2}}{J^{2}}+1\right)},$
(12)
$B_{c2}=\frac{2JMS}{g\mu_{B}}\sqrt{\left(\frac{D^{\prime}}{JM}+\Delta+\frac{g\mu_{B}B_{A}}{2JMS}\right)^{2}-\left(\frac{E_{-}^{2}}{J^{2}}+1\right)}.$
(13)

Once the Hamiltonian of the system is diagonalized, the statistical sum is
found as

$Z=\operatorname{Tr}\left(e^{-\beta H}\right)$ (14)

Where $\beta=1/k_{B}T$ with $k_{B}$ as the Boltzmann constant, $T$ as the
absolute temperature and $H$ the diagonalized Hamiltonian in Eq. (10). The
Helmholtz free energy of the system is related to the statistical sum by the
following relation:

$F=-T\ln Z.$ (15)

Thus, after calculation, the free energy is given by,

$\begin{split}&F=H^{\prime}_{0}+\frac{1}{2}\sum_{k}\left(\omega_{1k}^{\left(+\right)}+\omega_{2k}^{\left(+\right)}\right)+T\sum_{k}\ln\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)+T\sum_{k}\ln\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)\\\
&+\frac{1}{2}\sum_{k}\left(\omega_{1k}^{\left(-\right)}+\omega_{2k}^{\left(-\right)}\right)+T\sum_{k}\ln\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)+T\sum_{k}\ln\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)\end{split}$
(16)

where the frequencies obtained in Eq. (11) are rewritten in terms of the
system parameters as follow,

$\omega_{1k}^{\left(\pm\right)}=2JMS\sqrt{\left(\frac{D^{\prime}}{JM}+\Delta+\frac{g\mu_{B}B_{A}}{2JMS}\right)^{2}+\left(\frac{E_{+}^{2}}{J^{2}}+1\right)\left(\frac{2k^{2}l^{2}}{M}-1\right)}\pm
g\mu_{B}B_{+},$
$\omega_{2k}^{\left(\pm\right)}=2JMS\sqrt{\left(\frac{D^{\prime}}{JM}+\Delta+\frac{g\mu_{B}B_{A}}{2JMS}\right)^{2}+\left(\frac{E_{-}^{2}}{J^{2}}+1\right)\left(\frac{2k^{2}l^{2}}{M}-1\right)}\pm
g\mu_{B}B_{-}$

with $l$ as the length sides of the cubic primitive cell of each sublattice.

From the free energy, some magnetoelectric properties 45 such as the
magnetization, the electric polarization, and the magnetoelectric
polarizability are derived R13 ,

$\begin{split}M&=-\mu_{0}\frac{\partial F}{\partial B}\\\
&=-g\mu_{0}\mu_{B}\left[\frac{N}{2}+\sum_{k}\left(\frac{e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{1-e^{-\frac{1}{T}\omega_{1\,k}^{\left(+\right)}}}+\frac{e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}\right)\right]\\\
&+g\mu_{0}\mu_{B}\left[\frac{N}{2}+\sum_{k}\left(\frac{e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{1-e^{-\frac{1}{T}\omega_{1\,k}^{\left(-\right)}}}+\frac{e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}\right)\right]\end{split}$
(17) $\begin{split}P&=-\frac{\partial F}{\partial E}\\\
&=\frac{2MS}{J}\sum_{k}\gamma_{k}^{2}\left(\frac{E_{+}}{\xi_{+}^{1/2}}+\frac{E_{-}}{\xi_{-}^{1/2}}\right)+\frac{2MS}{J}\sum_{k}\gamma_{k}^{2}\left[\frac{E_{+}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)}+\frac{E_{-}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)}\right]\\\
&+\frac{2MS}{J}\sum_{k}\gamma_{k}^{2}\left[\frac{E_{+}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)}+\frac{E_{-}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)}\right]\end{split}$
(18) $\begin{split}\alpha&=-\mu_{0}\frac{\partial F}{\partial B\partial E}\\\
&=-\frac{2g\mu_{0}\mu_{B}MS}{JT}\left\\{\begin{split}&\sum_{k}\left[\frac{E_{+}\gamma_{k}^{2}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)^{2}}+\frac{E_{-}\gamma_{k}^{2}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)^{2}}\right]\\\
&-\sum_{k}\left[\frac{E_{+}\gamma_{k}^{2}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)^{2}}+\frac{E_{-}\gamma_{k}^{2}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)^{2}}\right]\end{split}\right\\}\end{split}$
(19)

respectively, where

$\xi_{+}=\left(\frac{g\mu_{B}B_{A}}{2JMS}+\Delta+\frac{D^{\prime}}{JM}\right)^{2}+\left(\frac{E_{+}^{2}}{J^{2}}+1\right)\left(\frac{2k^{2}l^{2}}{M}-1\right),$
$\xi_{-}=\left(\frac{g\mu_{B}B_{A}}{2JMS}+\Delta+\frac{D^{\prime}}{JM}\right)^{2}+\left(\frac{E_{-}^{2}}{J^{2}}+1\right)\left(\frac{2k^{2}l^{2}}{M}-1\right).$

##  3 Entropy and specific heat capacity

###  3.1 Entropy

The Boltzmann entropy which quantified the rate of the disorder in the system
is derived from the free energy as,

$S=-\frac{\partial F}{\partial T}$ (20)

Thus, after performing calculation the entropy is given by,

$\begin{split}&S=-\left[\sum_{k}\ln\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)+\sum_{k}\ln\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)\right]\\\
&+\frac{1}{T}\left[\sum_{k}\frac{\omega_{1k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)}+\sum_{k}\frac{\omega_{2k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)}\right]+\frac{1}{T}\left[\sum_{k}\frac{\omega_{1k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)}+\sum_{k}\frac{\omega_{2k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)}\right]\end{split}$
(21)

First of all, we note that the numerical results presented in the current
section and next section are obtained using the following parameters:
$B_{A}=4.14T,D^{\prime}=5.6,\Delta=0.5,J=1$ and $M=6$.

The temperature dependence of the entropy for different values of the electric
(upper panels) and magnetic (lower panels) site-dependent parameters is
depicted in Fig. (1). By varying the electric site-dependent parameter, we
observe that when this parameter is greater than $0.5$ the curves exhibit the
anomalous oscillating-like behavior with a negative part at low temperature
($T\mathrm{<}2.5K$) but the amplitude of the negative part decreases with the
increasing of the magnetic site-dependent parameter and completely disappear
when $\lambda=0.9$. That negative part of the entropy can be interpreting as
evidence of the cooling of the system (inverse magneto/electrocaloric effect).
In addition, by varying the magnetic site-dependent parameter, it is observed
that the site-dependent magnetic field affects the entropy at the oscillating-
like part. Such an entropy behavior clearly demonstrates the violation of the
second law of thermodynamic at low temperatures.

Moreover, the plots in Fig. (2) highlight the magnetic response of the entropy
for different values of temperature. It is showed that one peak-like point
occurs in the case of the uniform magnetic field ($\lambda=0$) as observed in
the previous work 25 . However, two peaks-like points occur in the case of a
site-dependent magnetic field ($\lambda\mathrm{\neq}0$) and this confirms the
two values of the critical magnetic fields obtained in the Eqs. (12) and (13).
When the magnetic site-dependent parameter is equal to $0.5$, the values of
the first and second critical magnetic fields are approximately $B_{c1}=7.5T$
and $B_{c2}=22.5T$ respectively, both for a uniform and site-dependent
electric field. It is observed that the parameter $\tau$ considerable affects
the amplitudes of the peaks whereas the parameter $\lambda$ affects the
position of theses peaks (critical magnetic fields). Indeed, when $\lambda$
increase the first critical point is shifted into the left meanwhile and the
second critical point is shifted into the right. This implies that by
increasing the parameter $\lambda$ one can effectively decrease and increase
the values of the critical magnetic fields $B_{c1}$ and $B_{c2}$ respectively.
Notice that the two peaks observed when the system is under the influence of
the site-dependent magnetic field indicate the existence of the intermediate
phase between the order and the disordered phase as observed in reference 46 .

On the other hand, in Fig. (3), it is observed that the entropy increases with
the electric field and reaches a relative maximum point from which it shows an
oscillating-like behavior both in the case of uniform and site-dependent
fields. Also, increasing the temperature as well as the magnetic site-
dependent parameter enhances the value of the electric field from which
appears the first maximum point (critical electric field) whereas the opposite
situation occurs when the electric site-dependent parameter increases. Note
that the critical electric field is always approximately localized between
$1V/m$ and $2.5V/m$, no matter the values of the parameters of the system
considered. It is worth noticing that beyond the critical electric field, the
system violates the second law of thermodynamics. Indeed, according to the
second law of thermodynamics, the entropy of an isolated system always
increase but hereby the entropy displays an oscillatory behavior from the
critical electric field. In addition, that oscillatory behaviour of the
entropy with negative part is due to the heat dissipation for the system under
external electric field and thermal condition. Note that the site-dependent
electric field enhance such a behaviour of the entropy also observed in Fig.
(1). Furthermore, the phase diagram in term of the entropy obtained within the
$E_{0}-T$ plane, in Fig. (4), shows three regions as function of approximated
values of critical electric field and temperature: a region where the entropy
varies slightly with the electric field ($T\mathrm{<}1K$), a region where the
entropy increases with both the electric field and the temperature
($E_{0}\mathrm{<}1V/m$) and a region where the entropy oscillates both with
the electric field and the temperature ($E_{0}\mathrm{>}1V/m$ and
$T\mathrm{>}1K$ ). Note that the critical electric and the critical
temperature are tunable by the magnetic and electric site-dependent
parameters.

Figure 1: the upper panels show the evolution of the entropy of the system for
different values of the electric site-dependent parameters and for three
values of the magnetic site-dependent parameter namely 0 (a), 0.5 (b) and 0.9
(c). In the lower panels we plotted the evolution of the entropy for different
values of the magnetic site-dependent parameters and for three values of the
electric site-dependent parameter namely 0 (d), 0.5 (e) and 0.9 (f).

Figure 2: magnetic field dependence of entropy by varying the temperature with
the following magnetic and electric site-dependent parameters: $\lambda=0$,
$\tau=0.9$ (a); $\lambda=0.5$, $\tau=0$ (b); $\lambda=0.5$, $\tau=0.5$ (c);
$\lambda=0.5$, $\tau=0.9$ (d) ; $\lambda=0.75$, $\tau=0.5$ (e).

Figure 3: electric field dependence of entropy by varying the temperature with
the following magnetic and electric site-dependent parameters: $\lambda=0.9$,
$\tau=0$ (a); $\lambda=0$, $\tau=0.5$ (b); $\lambda=0.9$, $\tau=0.9$ (c);
$\lambda=0.9$, $\tau=0.1$ (d) ; $\lambda=0.1$, $\tau=0.9$ (e).

Figure 4: the surface plot of entropy against electric field and temperature
with the following magnetic field, magnetic and electric site-dependent
parameters :(a)$B_{0}=1.85T,\tau=\lambda=0$, (b)
$B_{0}=1.85T,\lambda=0.9,\tau=0.1$, (c) $B_{0}=1.85T,\lambda=0.9,\tau=0.9$.

###  3.2 specific heat capacity

From the entropy obtained in Eq. (21), the specific heat capacity at constant
magnetic and electric fields is obtained as

$C=T\left(\frac{\partial S}{\partial T}\right).$ (22)

Thus

$C=\frac{1}{T^{2}}\left[\begin{split}&\sum_{k}\frac{\omega_{1k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)^{2}}+\sum_{k}\frac{\omega_{1k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)^{2}}+\sum_{k}\frac{\omega_{2k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)^{2}}\\\
&+\sum_{k}\frac{\omega_{2k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)^{2}}\end{split}\right].$
(23)

The temperature dependence of the specific heat capacity is shown in Fig. (5).
It is observed that the specific heat capacity increases with the temperature
and reaches an asymptotic value. The value of the temperature beyond which the
heat capacity becomes constant (critical temperature) together with the
corresponding asymptotic value depend upon the values of the magnetic and
electric site-dependent parameters considered. Also, in the absence of site-
dependent-electric field ($\tau=0$), the specific heat capacity is only
slightly affected by the site-dependent magnetic field. However, when the
site-dependent-electric field switch on ($\tau\mathrm{\neq}0$), the impacts of
the site-dependent magnetic field becomes more prominent. Moreover, for
certain values of magnetic and electric site-dependent parameters, the
specific heat capacity exhibits an abrupt increase or decrease near-zero
temperature, indicating the occurrence of metamagnetic and metaelectric
transition in the system 47 .

In Fig. (6), the magnetic field dependence of the specific heat capacity with
different values of the temperature displays a quasi-oscillating-like behavior
accompanied by peak-like points which appear at one point in the case of zero
magnetic site-dependent parameter ($\lambda=0$) and two different points
corresponding to the critical magnetic fields in the case of nonzero magnetic
site-dependent parameter ($\lambda\mathrm{\neq}0$).

Furthermore, the electric dependence of the specific heat capacity highlighted
in Fig. (7) shows that the specific heat capacity firstly increases, then
exhibits a quasi-oscillating-like behavior associated with peak points and
finally freezes or becomes constant. This implies that the specific heat
capacity of the system becomes constant when the electric field is strong
enough, demonstrating that the system is in a thermal equilibrium state. It is
worth noting that similar behaviors of specific heat capacity have been
observed in Refs 48 ; 49 .

Figure 5: the upper panels show the evolution of the specific heat capacity of
the system for different values of the electric site-dependent parameters and
for three values of the magnetic site-dependent parameter namely 0 (a), 0.1
(b) and 0.75 (c). In the lower panels we plotted the evolution of the specific
heat capacity for different values of the magnetic site-dependent parameters
and for three values of the electric site-dependent parameter namely 0 (d),
0.5 (e) and 0.9 (f).

Figure 6: magnetic field dependence of the specific heat capacity by varying
the temperature with the following magnetic and electric site-dependent
parameters: $\lambda=0$, $\tau=0.1$ (a); $\lambda=0.5$, $\tau=0$ (b);
$\lambda=0.5$, $\tau=0.5$ (c); $\lambda=0.5$, $\tau=0.1$ (d) ; $\lambda=0.1$,
$\tau=0.5$ (e).

Figure 7: electric field dependence of the specific heat capacity by varying
the temperature with the following magnetic and electric site-dependent
parameters: $\lambda=0.1$, $\tau=0$ (a); $\lambda=0$, $\tau=0.9$ (b);
$\lambda=0.5$, $\tau=0.5$ (c); $\lambda=0.5$, $\tau=0.9$ (d) ; $\lambda=0.1$,
$\tau=0.5$ (e).

##  4 Adiabatic cooling rate

In this section, the adiabatic cooling rates such as the adiabatic magnetic
cooling rate, the adiabatic electric cooling rate, and adiabatic
magnetoelectric cooling rate are evaluated. These parameters help to quantify
the magnetocaloric effect, the electrocaloric effect, and the magneto-
electrocaloric effect respectively. In fact, the adiabatic cooling rate can be
defined clearly as an adiabatic temperature change of a system as a response
of a suitable applied external field. From the Maxwell equations and the basic
thermodynamic equations, the adiabatic cooling rate in a nonlinear media is
generalized as 50 ,

$\Gamma=-\frac{T}{C}\frac{\partial X_{i}}{\partial T}$ (24)

where $X_{i}$ is generalized displacement (magnetization, electric
polarization, strain, etc.). Thus, the adiabatic magnetic, electric, and
magnetoelectric cooling rates are derived as

$\begin{split}\Gamma^{m}&=\frac{\partial T}{\partial
H}=-\frac{T}{C}\frac{\partial M}{\partial T}\\\
&=\frac{g\mu_{0}\mu_{B}}{TC}\left[\begin{split}&\sum_{k}\frac{\omega_{1k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)^{2}}+\sum_{k}\frac{\omega_{2k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)^{2}}-\sum_{k}\frac{\omega_{1k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)^{2}}\\\
&-\sum_{k}\frac{\omega_{2k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)^{2}}\end{split}\right]\end{split}$
(25) $\begin{split}&\Gamma^{e}=\frac{\partial T}{\partial
E}=-\frac{T}{C}\frac{\partial P}{\partial T}\\\
&=-\frac{2MS}{JTC}\left[\begin{split}&\sum_{k}\frac{\gamma_{k}^{2}E_{+}\omega_{1k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)^{2}}+\sum_{k}\frac{\gamma_{k}^{2}E_{-}\omega_{2k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)^{2}}+\sum_{k}\frac{\gamma_{k}^{2}E_{+}\omega_{1k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)^{2}}\\\
&+\sum_{k}\frac{\gamma_{k}^{2}E_{-}\omega_{2k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)^{2}}\end{split}\right]\end{split}$
(26)

and

$\begin{split}&\Gamma^{me}=\frac{\partial T}{\partial H\partial
E}=-\frac{T}{C}\frac{\partial\alpha}{\partial T}\\\
&=-\frac{2g\mu_{0}\mu_{B}MS}{JTC}\left[\begin{split}&\sum_{k}\frac{\gamma_{k}^{2}E_{+}\omega_{1k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)^{2}}+\sum_{k}\frac{\gamma_{k}^{2}E_{-}\omega_{2k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)^{2}}\\\
&-\sum_{k}\frac{\gamma_{k}^{2}E_{+}\omega_{1k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)^{2}}-\sum_{k}\frac{\gamma_{k}^{2}E_{-}\omega_{2k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)^{2}}\end{split}\right]\\\
&+\frac{2g\mu_{0}\mu_{B}MS}{JT^{2}C}\left[\sum_{k}\frac{\gamma_{k}^{2}E_{+}\omega_{1k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\left(1+e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(+\right)}}\right)^{3}}+\sum_{k}\frac{\gamma_{k}^{2}E_{-}\omega_{2k}^{\left(+\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\left(1+e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(+\right)}}\right)^{3}}\right]\\\
&-\frac{2g\mu_{0}\mu_{B}MS}{JT^{2}C}\left[\sum_{k}\frac{\gamma_{k}^{2}E_{+}\omega_{1k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\left(1+e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)}{\xi_{+}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{1k}^{\left(-\right)}}\right)^{3}}+\sum_{k}\frac{\gamma_{k}^{2}E_{-}\omega_{2k}^{\left(-\right)}e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\left(1+e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)}{\xi_{-}^{1/2}\left(1-e^{-\frac{1}{T}\omega_{2k}^{\left(-\right)}}\right)^{3}}\right]\end{split}$
(27)

respectively, where $C$ is the specific heat capacity at constant fields (Eq.
(23)), $M$ the magnetization (Eq. (17)), $P$ the electric polarization (Eq.
(18)), $\alpha$ is the magnetoelectric polarizability (Eq. (19)).

###  4.1 Adiabatic magnetic and electric cooling rates

####  4.1.1 Adiabatic magnetic cooling rate

In Fig. (8), the temperature response of the adiabatic magnetic cooling rate
highlights that $\Gamma^{m}$ exhibits alternating negative and positive peak
points (in the region of low temperature) and then a linear increase or
decrease with temperature (in the region of high temperature). It is observed
that by varying the magnetic site-dependent parameter the peak points almost
appear at the same values of temperature while they almost appear at different
values of temperature when the electric site-dependent parameter varies. In
the particular case in which the electric site-dependent parameter is equal to
$0.9$ with different values of the magnetic site-dependent parameter, all the
peaks with the same sign appear exactly at the same points, and the change
points were $\Gamma^{m}=0$ also occur at the same temperature see( Fig. (8)
panel (c)). Besides, for a certain combination of the magnetic and electric
site-dependent parameters, there is a formation of two successive negative
peaks (Fig. (8) panel (b)). These behaviors are in good accordance with the
characteristic behaviors of magnetocaloric effects as obtained in the
following experimental works 51 . So it is important to notice that the
positive peaks of the magnetic cooling rate demonstrates the heating of the
system (normal magnetocaloric effect) whereas the negative ones signals the
cooling of the system (inverse magneto caloric effect).

Moreover, the magnetic field dependence of the adiabatic magnetic cooling rate
is plotted in Fig. (9) for different values of the temperature. It is observed
that $\Gamma^{m}$ exhibits a peak point near to zero magnetic field followed
in some cases by a damped oscillating peaks behavior. Indeed, in the
particular case in which the electric site-dependent parameter is equal to
$0.9$ the adiabatic magnetic cooling rate does not exhibit any oscillatory
part. Furthermore, the electric field response of the adiabatic magnetic
cooling rate as plotted in Fig. (10) for different values of temperature shows
similar behavior as its temperature dependence (Fig. (8)). However, in
contrary to Fig. (8) where the magnetic cooling rate saturates at values
different from zero, here it is observed that the magnetic cooling rate always
cancels (saturate at zero) after the alternating peaks part.

Furthermore, the phase diagram in term of the magnetization obtained within
the $B_{0}-T$ plane, in Fig. (11), shows three regions (see panel (c)) as
function of approximated values of critical magnetic field and temperature: a
region where the magnetization varies slightly with the temperature
($T\mathrm{<}1K$), an intermediate region where the magnetization oscillates
with the temperature ($1K\mathrm{<}T\mathrm{<}6K$) and a region where the
magnetization increases both with the magnetic field and the temperature
($T\mathrm{>}6K$ ). Note that the critical magnetic field and the critical
temperature are tunable by the magnetic and electric site-dependent
parameters. Thus, the phase diagram in Fig. (11) summarizes the phase
transition observed by interpreting the magnetic cooling rate graphs.

Figure 8: the upper panels show the evolution of adiabatic magnetic cooling
rate of the system for different values of the magnetic site-dependent
parameters and for three values of the electric site-dependent parameter
namely 0 (a), 0.1 (b) and 0.9 (c). In the lower panels we plotted the
evolution of adiabatic magnetic cooling rate for different values of the
electric site-dependent parameters and for three values of the magnetic site-
dependent parameter namely 0 (d), 0.5 (e) and 0.9 (f).

Figure 9: magnetic field dependence of the adiabatic magnetic cooling rate by
varying the temperature with the following magnetic and electric site-
dependent parameters: $\lambda=0$, $\tau=0.9$ (a); $\lambda=0.9$, $\tau=0$
(b); $\lambda=0.1$, $\tau=0.1$ (c); $\lambda=0.1$, $\tau=0.9$ (d) ;
$\lambda=0.9$, $\tau=0.1$ (e).

Figure 10: electric field dependence of adiabatic magnetic cooling rate by
varying the temperature with the following magnetic and electric site-
dependent parameters: $\lambda=0.9$, $\tau=0$ (a); $\lambda=0$, $\tau=0.9$
(b); $\lambda=0.5$, $\tau=0.5$ (c); $\lambda=0.1$, $\tau=0.9$ (d) ;
$\lambda=0.9$, $\tau=0.1$ (e).

Figure 11: the surface plot of magnetization against magnetic field and
temperature with the following electric field, magnetic and electric site-
dependent parameters :(a)$E_{0}=0V/m,\lambda=0.1$, (b)
$E_{0}=3.16V/m,\lambda=\tau=0$, (c) $E_{0}=3.16V/m,\lambda=0.1,\tau=0.9$.

####  4.1.2 Adiabatic electric cooling rate

In Fig. (12), the adiabatic electric cooling rate is plotted as a function of
temperature and for different values of $\lambda$ and $\tau$. We observed that
the site-dependent magnetic field affects the temperature response of the
adiabatic electric cooling rate only for temperatures approximately less than
$3.7K$, both under uniform and the site-dependent electric field. The evidence
of that influence is demonstrated by the presence of the peak points at the
different or same temperature for different values of the magnetic site-
dependent parameter. When $T\mathrm{>}3.7K$, both the site dependent magnetic
and electric parameters have no more effect on the adiabatic electric cooling
rate. However, the site-dependent electric field affects the temperature
dependence of the electric cooling rate from $T\mathrm{>}3.5K$ which is the
temperature from which the adiabatic electric cooling rate exhibits
alternating negative and positive sign maxima (down panels).

The plots in Fig. (13) show the magnetic dependence of the adiabatic electric
cooling rate. It observed that when the electric field is uniform ($\tau=0$)
and the magnetic field site-dependent ($\lambda\mathrm{\neq}0$), the adiabatic
electric cooling rate can be increased by decreasing the temperature of the
system. The situation is completely reversed when the electric field becomes
site-dependent ($\tau\mathrm{\neq}0$) and the magnetic field uniform
($\lambda=0$). Overall, it is observed that the adiabatic electric cooling
rate can increase or decrease linearly or nonlinearly with the magnetic field,
depending on the temperature of the system. Indeed, when the magnetic and
electric site-dependent parameters are combined as follow $\lambda=0.1$,
$\tau=0.9$ and vice-versa, the nonlinear behavior of the electric cooling rate
is strongly attenuated (panels (d) and (e)).

Moreover, the curves in Fig. (14) highlighting the electric dependence of the
adiabatic electric cooling rate show that it displays alternating maxima at
weak electric field ($E_{0}\mathrm{<}6V/m$), and finally freezes at zero with
the presence of peak points for certain values of temperature. Overall, the
results obtained here show that the adiabatic electric cooling rate displays
the characteristic behavior of the electrocaloric effect as observed in
experimental results 12 ; 14 . Notice that the coexistence of positive and
negative values of the electric cooling indicates the heating (positive
electrocaloric effect) and cooling (negative electrocaloric effect) of the
system, respectively.

Furthermore, the phase diagram in term of the electric polarization obtained
within the $E_{0}-T$ plane, in Fig. (15), shows three regions as function of
approximated values of critical electric field and temperature: a region where
the electric polarization varies slightly with the electric field
($E_{0}\mathrm{<}1V/m$), a region where the electric polarization increases
with the electric field but oscillates with the temperature
($1V/m\mathrm{<}E_{0}\mathrm{<}2.1V/m$) and a region where the electric
polarization oscillates both with the electric field and the temperature
($E_{0}\mathrm{>}2.1V/m$ ). Note that the critical electric field and the
critical temperature are tunable by the magnetic and electric site-dependent
parameters.

Figure 12: the upper panels show the evolution of adiabatic electric cooling
rate of the system for different values of the magnetic site-dependent
parameters and for three values of the electric site-dependent parameter
namely 0 (a), 0.01 (b) and 0.025 (c). In the lower panels we plotted the
evolution of adiabatic electric cooling rate for different values of the
electric site-dependent parameters and for three values of the magnetic site-
dependent parameter namely 0 (d) 0.05 (e) and 0.075 (f).

Figure 13: magnetic field dependence of the adiabatic electric cooling rate by
varying the temperature with the following magnetic and electric site-
dependent parameters: $\lambda=0.5$, $\tau=0$ (a); $\lambda=0$, $\tau=0.25$
(b); $\lambda=0.1$, $\tau=0.1$ (c); $\lambda=0.1$, $\tau=0.9$ (d) ;
$\lambda=0.9$, $\tau=0.1$ (e).

Figure 14: electric field dependence of adiabatic electric cooling rate by
varying the temperature with the following magnetic and electric site-
dependent parameters: $\lambda=0$, $\tau=0.25$ (a); $\lambda=0.75$, $\tau=0$
(b); $\lambda=0.1$, $\tau=0.1$ (c); $\lambda=0.1$, $\tau=0.9$ (d) ;
$\lambda=0.9$, $\tau=0.1$ (e).

Figure 15: the surface plot of electric polarization against electric field
and temperature with the following magnetic field, magnetic and electric site-
dependent parameters :(a)$B_{0}=0T,\tau=0.1$, (b)
$B_{0}=1.5mT,\lambda=\tau=0$, (c) $B_{0}=1.5mT,\lambda=0.05,\tau=0.1$.

###  4.2 Adiabatic magnetoelectric cooling rate

Fig. (16) clearly shows in panel (a) that the site-dependent magnetic field
doesn’t affect the temperature response of the magnetoelectric cooling rate
under the uniform electric field ($\tau=0$). However, by varying the magnetic
site-dependent parameter under the site-dependent electric field, alternating
negative and positive maxima appear from $T\mathrm{>}1.9K$. This thus shows
the influence of both site-dependent fields (panels b and c). Note that such a
behavior of the magnetoelectric cooling indicates the magneto-electrocaloric
effect. In addition, the maximum value of the magneto-electric cooling
increase with the magnetic site-dependent parameter. Thus, the site-dependent
magnetic field enhance the magneto-electrocaloric effect. Moreover, by varying
the electric site-dependent parameter it is observed that the temperature
dependence of the magnetoelectric cooling rate exhibits an oscillating
behavior both when the system is suggested to a uniform magnetic field
($\lambda=0$) and site-dependent magnetic field ($\lambda\mathrm{\neq}0$)
(lower panels). Indeed, contrary to the case of uniform electric field where
the site-dependent magnetic field has no effect on the evolution of the
magnetoelectric cooling rate, it is observed that the magnetoelectric cooling
rate is affected by the site-dependent electric field when the magnetic field
is uniform.

Furthermore, the magnetic field dependence of the magnetoelectric cooling rate
is depicted in Fig. (17). When $\lambda=0$ and $\tau\mathrm{\neq}0$ and vice-
versa, it is observed an appearance of the negative peak points from certain
values of the magnetic field (up panels). However, under the influence of both
the site-dependent magnetic and electric fields (that is $\tau\mathrm{\neq}0$
and $\lambda\mathrm{\neq}0$), the curves display almost the same qualitative
behavior as the magnetic response of the electric cooling rate.

The electric field response of the magnetoelectric cooling rate as depicted in
Fig. (18) shows that it displays three successive maxima with a positive
middle one as observed in ref.52 . Besides, the curve is smooth in the case of
the site-dependent electric field than in the case of site-dependent magnetic
field and these peaks appear between $1V/m$ and $6V/m$. Note that by analogy
to the magnetic and electric cooling rates, the negative magnetoelectric
cooling rate demonstrates the cooling of the system ( inverse magneto-electro
caloric effect) whereas the positive magnetoelectric cooling rate signals the
heating of the system (magneto-electro caloric effect).

Furthermore, the phase diagram in term of the magnetoelectric polarizability
obtained within the $E_{0}-T$ plane, in Fig. (19), shows three regions as
function of approximated values of critical electric field and temperature: a
region where the magnetoelectric polarizability varies slightly with the
electric field ($E_{0}\mathrm{<}1.15V/m$), a region where the magnetoelectric
polarizability exhibits a peak-like points
($E_{0}\mathrm{>}1.15V/m,T\mathrm{<}1K$) and a region where the
magnetoelectric polarizability oscillates both with the electric field and the
temperature ($E_{0}\mathrm{>}1.15V/m,T\mathrm{>}1K$ ). Note that the critical
electric field and the critical temperature are tunable by the magnetic and
electric site-dependent parameters.

Figure 16: the upper panels show the evolution of adiabatic magnetoelectric
cooling rate of the system for different values of the magnetic site-dependent
parameters and for three values of the electric site-dependent parameter
namely 0 (a), 0.5 (b) and 0.9 (c). In the lower panels we plotted the
evolution of adiabatic magnetoelectric cooling rate for different values of
the electric site-dependent parameters and for three values of the magnetic
site-dependent parameter namely 0 (d) 0.05 (e) and 0.09 (f).

Figure 17: magnetic field dependence of the adiabatic magnetoelectric cooling
rate by varying the temperature with the following magnetic and electric site-
dependent parameters: $\lambda=0.5$, $\tau=0$ (a); $\lambda=0$, $\tau=0.25$
(b); $\lambda=0.1$, $\tau=0.1$ (c); $\lambda=0.1$, $\tau=0.9$ (d) ;
$\lambda=0.9$, $\tau=0.1$ (e).

Figure 18: electric field dependence of adiabatic magnetoelectric cooling rate
by varying the temperature with the following magnetic and electric site-
dependent parameters: $\lambda=0$, $\tau=0.75$ (a); $\lambda=0.25$, $\tau=0$
(b); $\lambda=0.1$, $\tau=0.1$ (c); $\lambda=0.1$, $\tau=0.9$ (d) ;
$\lambda=0.9$, $\tau=0.1$ (e).

Figure 19: the surface plot of magnetoelectric polarizability against electric
field and temperature with the following magnetic field, magnetic and electric
site-dependent parameters :(a)$B_{0}=0.15mT,\lambda=\tau=0$, (b)
$B_{0}=0.15mT,\lambda=0.1,\tau=0.5$, (c) $B_{0}=0.15mT,\lambda=\tau=0.5$.

##  5 Conclusion

In this work, the influence of both the site-dependent magnetic and electric
fields on the magneto-electrocaloric effect in a multiferroic
antiferromagnetic quantum spin system has been investigated. The Hamiltonian
is diagonalized by the help of the spin-wave theory and the free energy is
derived via the statistical physic. As analytical results, the physical
parameters such as the entropy, the specific heat capacity, the adiabatic
magnetic, electric, and magnetoelectric cooling rate are calculated. The
results displayed by the graphs obtained demonstrate that the site-dependent
magnetic and electric fields strongly affect the caloric effect and quantum
phase transitions operate in multiferroic antiferromagnets.

Indeed, by interpreting the magnetic field dependence of entropy, It is showed
that one peak-like point occurs in the case of the uniform magnetic field
($\lambda=0$) as observed in the previous work 25 . However, two peaks-like
points occur in the case of a site-dependent magnetic field
($\lambda\mathrm{\neq}0$) indicating the formation of the intermediate phase
between the order and the disordered phase. It is observed that the parameter
$\tau$ considerable affects the amplitudes of the peaks whereas the parameter
$\lambda$ affects the position of theses peaks (critical magnetic fields). In
fact, when $\lambda$ increase the first critical point is shifted into the
left meanwhile and the second critical point is shifted into the right. This
implies that by increasing the parameter $\lambda$ one can effectively
decrease and increase the values of the critical magnetic fields. That result
is confirmed by the magnetic field dependence of the specific heat capacity.
In addition, the electric field dependence of the entropy displays an
oscillatory behaviour from certain value of the electric field (critical
electric field).

On the other hand, the numerical results obtained for the adiabatic magnetic,
electric, and magnetoelectric cooling rate show a characteristic behavior of
the caloric or multi-caloric effect as observed by the experimentalist and in
other theoretical works. Indeed, the result display by the graphs exhibits the
alternating negative and positive peaks occur at corresponding change points
(critical temperature, critical magnetic or electric fields) tuneable by the
electric or/and magnetic fields. Note that the site-dependent magnetic field
enhance the magneto-electrocaloric effect.

Overall, it is demonstrated that the cumulative influence of the site-
dependent magnetic and electric fields allows us not only to reveal quantum
critical points hide in a multiferroic quantum spin system but also to control
the caloric or multi-caloric effect essential in the construction of the
solid-state refrigeration devices.

## Acknowledgement

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

## References

  * (1) Stern-Taulats E., Castán T., Mañosa L., Planes A., Mathur N., and Moya X., MRS Bulletin 43, 295-299 (2018).
  * (2) Shashwat An and and Umesh V Waghmare,Mater. Res. Express 1, 045503 (2014).
  * (3) Ivan A. Starkov, Abdulkarim A. Amirov, and Alexander S. Starkov, Materials Research Proceedings 9, 131-135 (2018).
  * (4) Joule J P,Phil. Trans. R. Soc. 91, 149 91 (1859).
  * (5) Pierre Weiss, Auguste Piccard, J. Phys. Theor. Appl. 7 103-109 (1917).
  * (6) Angelo Maiorino, et al., Energies 12, 1871 (2019)
  * (7) M. E. Zhitomirsky, Phys. Rev. B 67, 104421 (2003).
  * (8) Kurt Engelbrecht,J. Phys.: Energy 1, 021001 (2019).
  * (9) A. S. Starkov and I. A. Starkovi, Journal of Experimental and Theoretical Physics 119, 258–263 (2014).
  * (10) B. P. Alho et al., Journal of Applied Physics 109, 083942 (2011).
  * (11) Lluís Mañosa and Antoni Planes, J. Phys. D: Appl. Phys. 51, 070201 (2018).
  * (12) B. Li et al., EPL 102, 47004 (2013).
  * (13) Xavier Moya et al., Adv. Mat. 25, 1360-1365 (2013).
  * (14) I. Ponomareva and S. Lisenkov, Phys. Rev. Lett. 108, 167604 (2012).
  * (15) Z. Kutnjak, B. Rozic, R Pirc,materials science, (2015),10.1002/047134608X.W8244.
  * (16) Anna Grgnebohm et al. ,Energy Technol. 6, 1491–1511 (2018).
  * (17) Jie Wang et al., Journal of Applied Physics 115,164102 (2014).
  * (18) Yang-Bin Ma et al., Phys. Rev. B 94, 094113 (2016).
  * (19) Gang Bai, et al. , Ferroelectrics 531, 186-195 (2018).
  * (20) Ursic, H., Bobnar, V., Malic, B. et al., Sci Rep 6, 26629 (2016).
  * (21) A. S. Starkov and I. A. Starkov , Journal of Experimental and Theoretical Physics 119, 258–263 (2014).
  * (22) T. Moriya, Phys. Rev. 120, 91 (1960).
  * (23) I. A. Sergienko and E. Dagotto, Phys. Rev. B 73, 094434 (2006).
  * (24) Hosho Katsura, Naoto Nagaosa, and Alexander V. Balatsky, Phys. Rev. Lett. 95, 057205( 2005).
  * (25) M. Tchoffo, International Journal of Modern Physics B 33, 1950051 (2019).
  * (26) Jozef Sznajd, phys. Rev. B 97, 214410 (2018).
  * (27) Pradeep Thakur and P. Durganandini, phys. Rev. B, Vol. 97, 064413 (2018).
  * (28) Zhong Ming, Xu Hui, Liu Xiao-Xian, and Tong Pei-Qing, Chin. Phys. B, Vol. 9, 090313 (2013).
  * (29) Shuling Wang, Ruixue Li, Hua-Hua Fu, Linjie Ding, and Kailun Yao, Journal of Applied Physics, Vol. 116, 014101 (2014).
  * (30) N. Avalishvili, G. I. Japaridze, and G. L. Rossini., phys. Rev. B 99, 205159 (2019).
  * (31) Tianqi Chen, Chu Guo, Pinaki Sengupta, and Dario Poletti.,Phys. Rev. B Vol. 101,064417 (2020).
  * (32) Christian Trippe et al.,Phys. Rev. B 81, 054402 (2010).
  * (33) Topilko M., Krokhmalskii T., Derzhko O. et al., Eur. Phys. J. B 85, 278 (2012).
  * (34) Moradmard, H., Shahri Naseri, M. and Mahdavifar, S., J Supercond Nov Magn 27, 1265–1271(2014).
  * (35) Tatsuhiko N. Ikeda, Condens. Matter 4, 92 (2019).
  * (36) G. Kozlowski and K. Rcgaciu, phys. stat. sol. (b) 114, K5 (1982).
  * (37) V. S. Borovi et al., J. Exp. Theor. Phys. 53, 816-821 (1981).
  * (38) E. Z. Valiev, The Physics of Metals and Metallography104, 8–12 (2007).
  * (39) P J von Ranke et al, J. Phys.: Condens. Matter 21, 056004 (2009).
  * (40) Liu et al., Reviews of Modern Physics, APS 90, 015005 (2018).
  * (41) C. J. Hamer,O. Rojas,and J. Oitma, Phys. Rev. B 81, 214424 (2010).
  * (42) Satoshi Ejima et al., SciPost Phys. 5, 059 (2018).
  * (43) X. Z. Yuan, H. S. Goan and K. D. Zhu, New Journal of Physics 9, 219 (2007).
  * (44) Baran et al. Phys. Rev. B 98, 064415 (2018).
  * (45) Phillip Weinberg and Marin Bukov,SciPost Phys 7, 020 (2019).
  * (46) Ming-wen Xiao, arXiv[math-ph], (2009), arXiv:0908.0787[math-ph],
  * (47) Raffaele Resta, J. Phys.: Condens. Matter 22, 123201 (2010).
  * (48) W. D. Waladi Gueagni et al., Journal of Magnetism and Magnetic Materials 534, 167987 (2021).
  * (49) F. Amiri et al, J. Low Temp. Phys. 177, 203-216 (2014).
  * (50) Sachio Horiuchi et al., Chem. Sci 11, 6183–6192 (2020).
  * (51) F. K. Fodouop et al., Journal of Magnetism and Magnetic Material 474, 456-461 (2019).
  * (52) J. Hasanzadeh, Z. Feiznejad,S. Mahdavifar, J Supercond Nov Magn 27, 595–600 (2014).
  * (53) Ivan A. Starkov, Alexander S. Starkov, International Journal of Solids and Structures 100–101, 187–194 (2016).
  * (54) Xavier Moya, Lluís Mañosa, and Antoni Planes, Phys. Rev. B 75, 184412 (2007).
  * (55) Anis Biswas et al., Phys. Rev. B 87, 134420 (2013).