# Stretched-exponential relaxation in weakly-confined Brownian systems through
large deviation theory

Lucianno Defaveri  Department of Physics, Bar-Ilan University, Ramat-Gan
52900, Israel  Eli Barkai  Department of Physics, Institute of Nanotechnology
and Advanced Materials, Bar-Ilan University, Ramat-Gan 52900, Israel  David A.
Kessler  Department of Physics, Bar-Ilan University, Ramat-Gan 52900, Israel

###### Abstract

Stretched-exponential relaxation is a widely observed phenomenon found in
glassy systems. It was previously modeled with non-Markovian dynamics
reflecting a memory effect. Here, we study a Brownian particle under the
influence of a confining, albeit weak, potential field that grows with
distance as a sub-linear power law. We find that for this memoryless model,
observables display stretched-exponential relaxation. The probability density
function of the system is studied using a rate function ansatz. We obtain
analytically the stretched-exponential exponent along with an anomalous power-
law scaling of length with time. The rate function exhibits a point of
nonanalyticity, indicating a dynamical phase transition. In particular, the
rate function is double-valued both to the left and right of this point,
leading to four different rate functions, depending on the choice of initial
conditions and symmetry.

Introduction.– Anomalous relaxation, characterized by non-exponential decay,
is observed in a wide range of physical systems [1]. One class of such
behavior is stretched-exponential relaxation [2, 3, 4, 5, 6]. This has been
seen, for example, in disordered or heterogeneous systems, where the complex
interplay of interactions leads to a broad distribution of relaxation times
[7]. In the context of diffusion in disordered media, the heterogeneous nature
of the medium gives rise to a stretched-exponential pattern in particle
movement, as observed in disordered dielectrics [8, 9, 10]. Moreover, the
relaxation of systems with glassy dynamics often exhibits a stretched-
exponential behavior [11, 12].

However, the possibility of simple Brownian models exhibiting such relaxation
properties remains largely unexplored. In this Letter, we present a solution
to the Fokker-Planck equation with an external potential that grows with
distance as a sub-linear power law. Remarkably, we observe that the relaxation
of the various moments of the position to their equilibrium values follows an
anomalous stretched exponential. Formally, the Fokker-Planck equation (see Eq.
(2)) for $P(x,t)$, the probability density function (PDF), can be solved via
an eigenfunction expansion, and yields stretched-exponential relaxation.
However, it turns out that there is a simpler, more direct, approach, which
also yields additional physical insight. We find that in the long-time limit,
$P(x,t)$ takes the form of the exponential of a rate function which possesses
a scaling form, as usually results from a large-deviation formalism [13],
which looks at the far tails of the distribution of an observable. Large
deviation theory is a subject of much active interest in statistical physics,
[14, 15, 16, 17, 18, 19, 20, 21, 22, 23], in particular, due to the recent
discovery of dynamical phase transitions in the large-deviation behavior of
some model systems [24, 25, 26, 27, 28, 29, 30, 31]. These are nonanalytic
points of the rate function, and are so-called due to the analogy of the rate
function to an equilibrium free energy [17, 26]. We show the relationship
between the stretched-exponential relaxation and the presence of a dynamical
phase transition.

The rate function is defined as the logarithm of the PDF $P(x,t)$, divided by
a power of the time [17, 18, 26]

$\displaystyle\mathcal{I}(z)\equiv\lim_{\stackrel{{\scriptstyle
t,x\to\infty}}{{z=x/t^{\gamma}\textrm{fixed}}}}\frac{-\ln P(x,t)}{t^{\nu}}\,,$
(1)

with the anomalous time-exponent $\nu\neq 1$ and where $\mathcal{I}$ is a
function of the scaling variable $z\equiv x/t^{\gamma}$. The stretched-
exponential relaxation of observables will be governed by the same anomalous
exponent $\nu$. It should be noted that in the previously identified cases
with anomalous temporal scaling, the observable in question was nonlocal in
time, whereas here it is the PDF of $x$ itself that exhibits the anomalous
scaling. The appearance of a rate function in our problem implies the
surprising result that the anomalous stretched-exponential relaxation we
observe is a result of large deviations, i.e. the dynamics at large $x$.

We find that the rate function is multivalued, and possesses a critical-point
$z_{c}$. The rate function has two possible branches for $z<z_{c}$ and two for
$z>z_{c}$, all meeting at $z_{c}$. All four possible combinations of branches
below and above $z_{c}$ have different interpretations, corresponding to
different classes of initial conditions and parity. Two of the combinations
have a jump discontinuity in $\mathcal{I}^{\prime\prime}(z)$, marking a
dynamical phase transition. Such a multi-valued rate function appears not to
have previously been encountered. Using the appropriate branches of the rate
function, we obtain the time scales of the stretched-exponential relaxation of
the even and odd moments of $x$.

Model.– We study non-interacting Brownian particles, in contact with a heat
bath at temperature $T$, that are also subject to a binding potential field
$V(x)$. The spatial spreading of the cloud of particles can be described via
the PDF $P(x,t)$, obtained from the Fokker-Planck equation

$\displaystyle\frac{\partial}{\partial
t}P(x,t)=D\left[\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial}{\partial
x}\left(\frac{V^{\prime}(x)}{k_{B}T}\right)\right]P(x,t)\,,$ (2)

where $D$ is the diffusion coefficient. We consider herein even potentials,
satisfying $V(x)=V(-x)$, which for large $x$ grow as a sublinear power-law,
$V(x)\propto x^{\alpha}$, where $0<\alpha<1$. This means that the force,
$F(x)=-V^{\prime}(x)$, will be negligible for large $x$, as $F(x)\propto
x^{\alpha-1}$. At long times, the particles will reach the stationary
equilibrium Boltzmann-Gibbs state $P_{\mathrm{BG}}(x)=\exp[-V(x)/k_{B}T]/Z$,
with $Z$ being the normalizing partition function. For the more commonly
considered superlinear growth of the potential, $\alpha>1$, $F(x)$ grows with
distance and everything is standard. The system exponentially relaxes to the
Boltzmann-Gibbs state, at a rate given by the first nonzero eigenvalue of the
Fokker-Planck operator, which has a discrete spectrum starting at 0. This
discreteness follows from the fact that under a similarity transformation, the
Fokker-Planck equation becomes a Schrödinger equation [32] with an effective
potential $V_{S}(x)=F(x)^{2}/4k_{B}T+F^{\prime}(x)/2k_{B}T$, which grows
without bound as $x\to\infty$. However, for $\alpha<1$, $V_{S}(x)$ decays at
large $x$ as $x^{-2(1-\alpha)}$ and so the spectrum goes continuously down to
0. Potentials of this form have already been studied in the context of
resetting processes [33] and active processes [34], in particular, the
$\alpha=1$ limiting case [35, 36, 37]. For $\alpha\to 0$ (and assuming
$V(x)\sim x^{\alpha}/\alpha$), we have that for large $x$, $V(x)\sim\ln x$. In
that limit, the relaxation has been shown to be governed by a power law [38,
39, 40, 41]. The question is then what happens for $0<\alpha<1$?

For our numerical examples, we use the family of potentials

$\displaystyle V(x)=V_{0}\left(1+\frac{x^{2}}{\ell^{2}}\right)^{\alpha/2}\,.$
(3)

For $x\gg\ell$, with $\ell$ being the lengthscale of the center region, the
potential exhibits the desired power law growth. Throughout this Letter, we
shall scale the position variable by $\ell$ and correspondingly the time by
$\ell^{2}/D$.

As mentioned above, the $P(x,t)$ can be decomposed in a sum of eigenfunctions,
each decaying as $e^{-\lambda t}$ with a continuous eigenvalue spectrum
$\lambda$ starting at $0^{+}$, together with the Boltmann-Gibbs bound state at
$\lambda=0$. It turns out that the dominant continuum contribution to $P(x,t)$
for large times comes from the vicinity of a particular finite eigenvalue
$\lambda_{*}$. This dominant eigenvalue scales as a negative power of the
time, $\lambda_{*}(t)\sim t^{\nu-1}$, $0\leq\nu<1$, so that
$e^{-\lambda_{*}(t)t}~{}\sim e^{-t^{\nu}}$, a stretched-exponential relaxation
to equilibrium. This stretched-exponential relaxation is easily demonstrated
numerically (see below). The eigenvalue calculation will be sketched in the
Supplemental Material (SM). We turn now instead to the rate-function
calculation. Most strikingly, we shall see that the power of $t$ in the
stretched exponential can be obtained immediately from our scaling ansatz for
the rate function.

Figure 1: All possible rate functions, defined by Eq. (1), which are solutions
of the differential equation in Eq. (6), versus the scaled position $z\equiv
x/t^{\gamma}$ (the scaled exponents defined in Eq. (5)) for $\alpha=1/2$. The
rate function has two low-$z$ and two high-$z$ branches, leading to four
different possible forms. The region where the discriminant of Eq. (8) is
negative is highlighted in gray.

Large deviation formalism.– From Eq. (1) we can write the PDF, up to pre-
exponential factors, as $P(x,t)\sim e^{-t^{\nu}\mathcal{I}(z)}$. Inserting
this ansatz into the Fokker-Planck equation (2), we find that in the long-time
limit

$\displaystyle\frac{z\gamma\mathcal{I}^{\prime}(z)-\nu\mathcal{I}(z)}{t^{1-\nu}}=\frac{\mathcal{I}^{\prime}(z)^{2}}{t^{2\gamma-2\nu}}-\frac{V_{0}}{k_{B}T}\frac{\alpha}{z^{1-\alpha}}\frac{\mathcal{I}^{\prime}(z)}{t^{(2-\alpha)\gamma-\nu}}.$
(4)

We impose that the time dependence of all terms is the same, namely,
$1-\nu=2\gamma-2\nu=(2-\alpha)\gamma-\nu$, and find that the scaling exponents
are

$\displaystyle\nu=\frac{\alpha}{2-\alpha}~{}\mathrm{and}~{}\gamma=\frac{1}{2-\alpha}\,.$
(5)

With these exponents, both the exponent of the Boltzmann-Gibbs solution,
$V(x)/k_{B}T=V_{0}/k_{B}Tx^{\alpha}=t^{\nu}V_{0}/k_{B}Tz^{\alpha}$, and of
free-diffusion, $x^{2}/4Dt=t^{\nu}z^{2}/4D$, are compatible with our scaling.

The rate function can be found by solving the non-linear differential equation
(4), which now reads

$\displaystyle\mathcal{I}^{\prime}(z)^{2}-\left(\frac{V_{0}}{k_{B}T}\frac{\alpha}{z^{1-\alpha}}+\frac{z}{2-\alpha}\right)\mathcal{I}^{\prime}(z)+\frac{\alpha\mathcal{I}(z)}{2-\alpha}=0.$
(6)

When we consider the large-$z$ limit, the rate function solution will assume
the form $\mathcal{I}(z)\approx\xi_{0}z^{\mu_{0}}$, where $\mu_{0}$ and
$\xi_{0}$ must be determined. Plugging $\mathcal{I}(z)$ into Eq. (6) in the
limit of large $z$, we have

$\displaystyle\xi_{0}\mu_{0}^{2}z^{\mu_{0}-2}-\frac{V_{0}}{k_{B}T}\alpha\mu_{0}z^{\alpha-2}\sim\left(\frac{\mu_{0}-\alpha}{2-\alpha}\right)\,.$
(7)

There are two possible solutions for this equation. First, we have
$\mu_{0}=\alpha$, leading to $\mathcal{I}(z)=V_{0}z^{\alpha}/k_{B}T$, i.e. the
mentioned Boltzmann-Gibbs state, which is a solution of Eq. (6) for all $z$.
Second, we have $\mu_{0}=2$ and $\xi_{0}=1/4$, which is equivalent to large
$z$ diffusive behavior. We emphasize that these are the only possible
asymptotic solutions of Eq. (6). Both these behaviors are necessary. The
Boltzmann-Gibbs state is a possible time-independent solution of the problem,
when the initial state is the equilibrium state. However, under any initial
conditions that decay with $x$ faster than the Boltzmann-Gibbs state and for
any finite time $t$, we cannot expect the Boltzmann-Gibbs state to describe
the whole PDF. The particles cannot spread faster than what is permitted by
diffusion. Since $F(x)\to 0$ as $x\to\infty$, we expect diffusive behavior at
large $|x|$, consistent with the second asymptotic behavior.

The next step is to solve Eq. (6) globally. Since Eq. (6) is a second-order
polynomial in ${\cal I}^{\prime}(z)$, we have in fact two different options
for the ODE:

$\displaystyle\mathcal{I}^{\prime}(z)$
$\displaystyle=\frac{1}{2}\left(\frac{V_{0}}{k_{B}T}\frac{\alpha}{z^{1-\alpha}}+\frac{z}{2-\alpha}\right)$
$\displaystyle{}\pm\frac{1}{2}\sqrt{\left(\frac{V_{0}}{k_{B}T}\frac{\alpha}{z^{1-\alpha}}+\frac{z}{2-\alpha}\right)^{2}-\frac{4\alpha\mathcal{I}(z)}{2-\alpha}}\,.$
(8)

These two ODEs give rise to two smooth solutions that cross at a critical
point $z_{c}$, at which the square-root vanishes. One solution has the
positive sign of the square root for $z<z_{c}$ and the negative sign for
$z>z_{c}$ and the other with the opposite choices. Interestingly, two other
solutions are also possible, where the sign does not switch, and which have a
jump in ${\cal I}^{\prime\prime}$, leaving us with four possible descriptions,
two for $z<z_{c}$ and two for $z>z_{c}$, as shown in Fig. 1. The next task is
to uncover the physical content of these rate-function branches.

Boltzmann-Gibbs thermal initial condition.– As noted above, if we were to
start at time $t=0$ with the thermal state in all of space, the state would
remain unchanged for all time. In Fig. 1, this is shown as the curve $ACE$.
The square-root in Eq. (8) vanishes at a critical point $z_{c}$,

$\displaystyle
z_{c}=\left(\frac{\alpha(2-\alpha)V_{0}}{k_{B}T}\right)^{1/(2-\alpha)}\,.$ (9)

The pure thermal state is obtained when we switch from positive ($z<z_{c}$) to
negative ($z>z_{c}$) in Eq. (8).

Localized initial condition.– Our goal is to associate the solutions shown in
Fig. 1 with classes of initial conditions and with parity. For an initially
localized packet of particles, we expect the large $z$ behavior to be
diffusive rather than Boltzmann-Gibbs. However, for small $z$, we expect the
behavior to match Boltzmann-Gibbs, so that the Boltzmann-Gibbs regime in $x$
expands as $t^{\gamma}$. Keeping the positive sign of the square-root for
$z>z_{c}$ leads to the desired diffusive-like large $z$ behavior. The
resulting singularity in ${\cal I}^{\prime\prime}$ is precisely the dynamical
phase transition. The resulting curve, $ACD$ in Fig. 1, describes the
localized initial condition.

Figure 2: Comparison between theoretical and numerical (at finite times) rate
functions, Eq. (1), as a function of the scaled position $z\equiv
x/t^{\gamma}$. The solid lines represent numerical results obtained
integrating Eq. (2) at different times, indicated in the legend (panel (a))
for all panels. In panel (a) we show the clear convergence to the rate
function $\mathcal{I}_{ACD}(z)$ and highlight (white-filled circle) the
critical transition point $z_{c}$. In panel (b) we show the derivative
$\mathcal{I}^{\prime}(z)$, and the non-analytical behavior becomes clear. The
solution is equivalent to Boltzmann-Gibbs up until the critical point, where
the system changes to the free-particle solution for larger values of $z$. In
panel (c) we have the odd rate functions, where there is no dynamical phase
transition. The numerical solutions demonstrate clear convergence towards the
expected rate functions. We have used $\alpha=3/4$, $\ell=1$, $D=1$ and
$V_{0}/k_{B}T=1$.

We can write odd/even PDFs using distinct rate functions as

$\displaystyle P_{\mathrm{even}}\sim
e^{-t^{\nu}\mathcal{I}_{\mathrm{even}}(z)}~{}\mathrm{and}~{}P_{\mathrm{odd}}\sim
e^{-t^{\nu}\mathcal{I}_{\mathrm{odd}}(z)}\,.$ (10)

If the particles start at the origin ($x_{0}=0$), by symmetry,
$P_{\mathrm{odd}}(x,t)=0$, since $V(x)$ is even and so all odd moments vanish,
as they do in equilibrium. On the other hand, for $x_{0}\neq 0$, the odd part
is present, $P_{\mathrm{odd}}=(P(x,t)-P(-x,t))/2$, and is described by the
curve $BCD$. We then write that
$\mathcal{I}_{\mathrm{odd}}(z)\equiv\mathcal{I}_{BCD}(z)$, which tends to a
constant value $\mathcal{I}_{\mathrm{odd}}(0)$ as $z$ approaches zero, see
Fig. 1. Therefore, we have that $P_{\mathrm{odd}}(x,t)$ has an upper bound
that decays as $e^{-\mathcal{I}_{\mathrm{odd}}(0)t^{\nu}}$, indicating that
the odd contributions will decay as a stretched exponential. Curve $BCE$
describes the contribution from the continuum modes. This contribution is
obtained by subtracting from $P(x,t)$ the Boltzmann-Gibbs solution,
$P^{*}\equiv P-P_{\mathrm{BG}}$. At large $x$, this results in the negative of
the Boltzmann-Gibbs solution, as the full solution is small. The stretched-
exponential relaxation of the even moments of $x$ to their equilibrium values
is controlled by $P^{*}$. The rate function, $BCE$ in Fig. 1, also displays a
dynamical phase transition at the critical point. Finally, the even part,
$P_{\mathrm{even}}=(P(x,t)+P(-x,t))/2$, is described by the $ACD$ rate
function, $\mathcal{I}_{\mathrm{even}}(z)\equiv\mathcal{I}_{ACD}(z)$. This
solution also captures the leading behavior of the density $P(x,t)$ for
localized initial conditions.

We presented the four different scenarios for the rate function represented by
branches in Fig. 1. Which scenario is the relevant one is based on the choice
of the initial condition and the parity. We now use finite time numerical
integration of the FPE to show that, using Eq. (1) (adapted for finite times),
the results converge to the expected rate functions. In Fig. 2(a), for
localized initial conditions, we show the numerical convergence to the rate
function $ACD$, while in Fig. 2(b) we show the derivative of the curve $ACD$.
The derivative clearly shows the non-analytical behavior at the point $z_{c}$.
In Fig. 2(c), we compare our long-time prediction for
$\mathcal{I}_{\mathrm{odd}}(z)$ ($I_{BCD}(z)$) with numerical results
($t^{-\nu}\log P_{\mathrm{odd}}$, with $P_{\mathrm{odd}}$ obtained
numerically), showing clear convergence.

In summary, we find that two of the rate functions are completely analytical.
Those are the Boltzmann-Gibbs (curve $ACE$) and the odd rate functions (curve
$BCD$). The other two possibilities, the localized initial condition rate
function (curve $ACD$) and the continuum modes rate function (curve $BCE$)
have a non-analytical behavior, and therefore a dynamical phase transition, at
the critical point $z_{c}$.

Remark on the notation and prefactors.– We highlight that the rate functions
defined in Eq. (10) satisfy
$\mathcal{I}_{\mathrm{even}/\mathrm{odd}}(z)=\mathcal{I}_{\mathrm{even}/\mathrm{odd}}(-z)$
as these functions are even. Further, they do not depend explicitly on the
initial conditions. The parity of $P_{\mathrm{even}/\mathrm{odd}}$ is
determined by the pre-exponential factors. The exponential prefactors of
$P_{\mathrm{odd}}$ will depend on the initial condition, since for symmetric
initial conditions, the whole odd part must vanish. The prefactor is obtained
using the WKB method in the SM.

Figure 3: The log of the numerical ensemble average of the position for a
system initialized with localized condition $P(x,t=0)=\delta(x-x_{0})$ for
different values of $\alpha$ (shown in legend). The stretched-exponential
behavior is clearly shown for long times, where we compare with
$\mathcal{I}_{\mathrm{odd}}(0)t^{\nu}$, where the anomalous time scale
$\mathcal{I}_{\mathrm{odd}}(0)$ is obtained in Eq. (12), and
$\nu=\alpha/(2-\alpha)$ (dashed black lines). On the inset, the complete
theoretical prediction (dotted black lines), described by Eq. (14), is
compared with the numerical simulations (colored lines) for the same exponents
$\alpha$ as in the main label. We have used $V_{0}/k_{B}T=1$, $\ell=1$ and
$x_{0}=0.04$.

Anomalous relaxation.– Here we study the relaxation properties of the system,
focusing on the first moment of the position $x$, starting with an asymmetric
initial condition ($x_{0}\neq 0$). Because this observable, $x$, is odd, the
only non-zero contribution comes from the asymmetric part of the PDF. The odd
part of the PDF can be written using $\mathcal{I}_{\mathrm{odd}}$, up to pre-
exponential factors, as shown in Eq. (10). We obtain, up to pre-exponential
factors, the stretched-exponential characteristic of the mean as

$\displaystyle\langle
x\rangle=\int_{-\infty}^{\infty}x\,P_{\mathrm{odd}}(x,t)dx\sim
e^{-\mathcal{I}_{\mathrm{odd}}(0)t^{\nu}}\,.$ (11)

The value of $\mathcal{I}_{\mathrm{odd}}(0)$ governs the anomalous time-scale
$\tau$ of the relaxation. It is possible to obtain this value numerically by
integrating Eq. (8) in the correct branch (positive sign). We have obtained an
analytical expression through our eigenfunction calculations (see SM)

$\displaystyle\frac{1}{\tau^{\nu}}=\mathcal{I}_{\mathrm{odd}}(0)=\frac{\left(\sqrt{\pi}\left(\frac{\alpha
V_{0}}{2k_{B}T}\right)^{\frac{1}{1-\alpha}}\frac{\Gamma\left(\frac{\alpha}{2-2\alpha}\right)}{\Gamma\left(\frac{1}{2-2\alpha}\right)}\right)^{1-\nu}}{\nu^{\nu}(1-\nu)^{1-\nu}}.$
(12)

It is remarkable that a rate function controls the relaxation, in the sense
that it cannot be considered describing a rare event, nor is it particularly
hard to measure it. In Fig. 3, the stretched-exponential behavior is shown
numerically in the long-time limit. The observed time scale matches our
prediction in Eq. (12) and is independent of the odd moment under
consideration. The same can be extended for even moments $x^{2n}$ (see SM). In
order to obtain a complete expression for the mean in Eq. (11), the rate
function is not enough. We must account for the pre-exponential factor,
$\mathcal{A}(t)$, that is,
$P_{\mathrm{odd}}\approx\mathcal{A}(t)\,x^{1-\alpha}e^{-\mathcal{I}_{\mathrm{odd}}(z)t^{\nu}}$.

In the long-time limit, the main contribution to the integral in Eq. (11)
arises from the region where $x$ is much smaller than the critical $x_{c}(t)$,
corresponding to the small $z$ region [42]. For small $z$, we have
$\mathcal{I}_{\mathrm{odd}}(z)\approx\mathcal{I}_{\mathrm{odd}}(0)(1+k_{B}Tz^{2-\alpha}/V_{0}(2-\alpha)^{2})$,
and the pre-exponential factor is (see SM)

$\displaystyle\mathcal{A}(t)=x_{0}\frac{k_{B}T}{(2-\alpha)V_{0}}\frac{\mathcal{I}_{\mathrm{odd}}(0)e^{\frac{V(0)}{k_{B}T}}}{t^{3/2-\nu}}\sqrt{\frac{\gamma-\nu}{4\pi}}\,.$
(13)

Note that for the transition value $\alpha=1$, $\gamma=\nu=1$, and
$\mathcal{A}(t)$ will be null. With the contribution of the prefactor, we
obtain, from Eq. (11), the long-time expression for the relaxation,

$\displaystyle\frac{\langle x\rangle}{x_{0}}$ $\displaystyle\sim$
$\displaystyle\frac{\mathcal{C}_{1}\left(\frac{V_{0}}{k_{B}T},\alpha\right)}{t^{\nu^{2}/2}}e^{-\mathcal{I}_{\mathrm{odd}}(0)t^{\nu}}\,,$
(14)

where the definition of $\mathcal{C}_{1}\left({V_{0}}/{k_{B}T},\alpha\right)$
is found in the SM. Thus, the time relaxation of the system to equilibrium is
through a stretched exponential (multiplied by a power law in time). We show
the excellent agreement of Eq. (14) with the numerical results in the inset of
Fig. 3.

Discussion.– The results we obtained in this Letter are quite general, but
nevertheless, many extensions to this work are possible. As a first step in
this direction, we studied a case in dimensions higher than one, showing that
the main results remain valid. The characteristics of time-averaged
observables and the potential link between different branches and
singularities in the cumulant generating function warrant attention [43, 44].
It is likely that the multivalued nature of the rate function under study,
which depends on symmetry and initial condition, is an important feature for
other systems. Thus, the appearance of the multivalued rate function in more
systems and its relationship with dynamical phase transitions requires further
study.

Acknoledgements: The support of Israel Science Foundation’s grant 1614/21 is
acknowledged. LD and EB thank Naftali Smith for the insightful conversations.

## References

  * Uchaikin and Sibatov [2013] V. Uchaikin and R. Sibatov,  _Fractional Kinetics in Solids_ (WORLD SCIENTIFIC, 2013).
  * Kohlrausch [1854] R. Kohlrausch, Theorie des elektrischen rückstandes in der leidener flasche, Annalen der Physik und Chemie 91, 179 (1854).
  * Klafter and Shlesinger [1986] J. Klafter and M. F. Shlesinger, On the relationship among three theories of relaxation in disordered systems, Proc. Nat’l. Acad. Sci. 83, 848 (1986).
  * Kisslinger [1993] C. Kisslinger, The stretched exponential function as an alternative model for aftershock decay rate, J. Geophys. Res. Solid Earth 98, 1913 (1993).
  * Flomenbom  _et al._ [2005] O. Flomenbom, K. Velonia, D. Loos, S. Masuo, M. Cotlet, Y. Engelborghs, J. Hofkens, A. E. Rowan, R. J. M. Nolte, M. V. der Auweraerand F. C. de Schryver, and J. Klafter, Stretched exponential decay and correlations in the catalytic activity of fluctuating single lipase molecules, Proc. Natl. Acad. Sci. 102, 2368 (2005).
  * Mukherjee  _et al._ [2023] S. Mukherjee, P. Pareek, M. Barma, and S. K. Nandi, From stretched exponential to power-law: crossover of relaxation in a simple model,  (2023), arXiv:2307.01801.
  * Malacarne  _et al._ [2001] L. C. Malacarne, R. S. Mendes, I. T. Pedron, and E. K. Lenzi, Nonlinear equation for anomalous diffusion: Unified power-law and stretched exponential exact solution, Phys. Rev. E 63, 030101(R) (2001).
  * Milovanov  _et al._ [2007] A. V. Milovanov, K. Rypdal, and J. J. Rasmussen, Stretched exponential relaxation and ac universality in disordered dielectrics, Phys. Rev. B 76, 104201 (2007).
  * Metzler and Klafter [2002] R. Metzler and J. Klafter, From stretched exponential to inverse power-law: Fractional dynamics, cole-cole relaxation processes, and beyond, Journal of Non-Crystalline Solids 305, 81 (2002).
  * Magdziarz  _et al._ [2005] M. Magdziarz, H. Steinhaus, and K. Weron, Anomalous diffusion schemes underlying the stretched exponential relaxation. the role of subordinators, Acta Physica Polonica, Series B 36, 1823 (2005).
  * Kakalios  _et al._ [1987] J. Kakalios, R. A. Street, and W. B. Jackson, Stretched-exponential relaxation arising from dispersive diffusion of hydrogen in amorphous silicon, Phys. Rev. Lett. 59, 1037 (1987).
  * Wu  _et al._ [2018] Z. W. Wu, W. Kob, W.-H. Wang, and L. Xu, Stretched and compressed exponentials in the relaxation dynamics of a metallic glass-forming melt, Nat Commun 9, 5334 (2018).
  * Freidlin and Wentzell [1998] M. I. Freidlin and A. D. Wentzell,  _Random perturbations of dynamical systems_ (Springer-Verlag, New York, 1998).
  * Touchette [2009] H. Touchette, The large deviation approach to statistical mechanics, Physics Reports 478, 1 (2009).
  * Baek and Kafri [2015] Y. Baek and Y. Kafri, Singularities in large deviation functions, Journal of Statistical Mechanics: Theory and Experiment 2015, P08026 (2015).
  * Baek  _et al._ [2018] Y. Baek, Y. Kafri, and V. Lecomte, Dynamical phase transitions in the current distribution of driven diffusive channels, Journal of Physics A: Mathematical and Theoretical 51, 105001 (2018).
  * Nickelsen and Touchette [2018] D. Nickelsen and H. Touchette, Anomalous scaling of dynamical large deviations, Phys. Rev. Lett. 121, 090602 (2018).
  * Meerson [2019] B. Meerson, Anomalous scaling of dynamical large deviations of stationary gaussian processes, Phys. Rev. E 100, 042135 (2019).
  * Nyawo and Touchette [2016] P. T. Nyawo and H. Touchette, A minimal model of dynamical phase transition, Europhysics Lett. 116, 50009 (2016).
  * Touchette [2018] H. Touchette, Introduction to dynamical large deviations of markov processes, Physica A 504, 5–19 (2018).
  * Jack and Harris [2020] R. L. Jack and R. J. Harris, Giant leaps and long excursions: Fluctuation mechanisms in systems with long-range memory, Phys. Rev. E 102, 012154 (2020).
  * Hurtado-Gutiérrez  _et al._ [2023] R. Hurtado-Gutiérrez, P. I. Hurtado, , and C. Pérez-Espigares, Spectral signatures of symmetry-breaking dynamical phase transitions, Phys. Rev. E 108, 014107 (2023).
  * Teza  _et al._ [2023] G. Teza, R. Yaacoby, and O. Raz, Eigenvalue crossing as a phase transition in relaxation dynamics, Phys. Rev. Lett. 130, 207103 (2023).
  * Janas  _et al._ [2016] M. Janas, A. Kamenev, and B. Meerson, Dynamical phase transition in large-deviation statistics of the Kardar-Parisi-Zhang equation, Phys. Rev. E 94, 032133 (2016).
  * Nemoto  _et al._ [2017] T. Nemoto, R. L. Jack, and V. Lecomte, Finite-size scaling of a first-order dynamical phase transition: Adaptive population dynamics and an effective model, Phys. Rev. Lett. 118, 115702 (2017).
  * Smith [2022a] N. R. Smith, Large deviations in chaotic systems: Exact results and dynamical phase transition, Phys. Rev. E 106, L042202 (2022a).
  * Smith [2022b] N. R. Smith, Exact short-time height distribution and dynamical phase transition in the relaxation of a Kardar-Parisi-Zhang interface with random initial condition, Phys. Rev. E 106, 044111 (2022b).
  * Smith [2022c] N. R. Smith, Anomalous scaling and first-order dynamical phase transition in large deviations of the ornstein-uhlenbeck process, Phys. Rev. E 105, 014120 (2022c).
  * Mukherjee and Smith [2023] S. Mukherjee and N. R. Smith, Dynamical phase transition in the occupation fraction statistics for non-crossing brownian particles (2023), arXiv:2303.17250.
  * Gutiérrez  _et al._ [2023] R. Gutiérrez, A. Canella-Ortiz, and C. Pérez-Espigares, Making rare events typical in chaotic maps (2023), arXiv:2304.13754.
  * Agranov  _et al._ [2023] T. Agranov, M. E. Cates, and R. L. Jack, Tricritical behavior in dynamical phase transitions, Phys. Rev. Lett. 131, 017102 (2023).
  * Risken [1996] H. Risken,  _The Fokker-Planck Equation_, Vol. 18 (Springer Berlin Heidelberg, 1996).
  * Capała and Dybiec [2023] K. Capała and B. Dybiec, Optimization of escape kinetics by reflecting and resetting,  (2023), arXiv:2305.14832.
  * Dhar  _et al._ [2019] A. Dhar, A. Kundu, S. N. Majumdar, S. Sabhapandit, and G. Schehr, Run-and-tumble particle in one-dimensional confining potentials: Steady-state, relaxation, and first-passage properties, Phys. Rev. E 99, 032132 (2019).
  * Chen  _et al._ [2013] Y. Chen, A. Baule, H. Touchette, , and W. Just, Weak-noise limit of a piecewise-smooth stochastic differential equation, Phys. Rev. E 88, 052103 (2013).
  * Sabhapandit and Majumdar [2020] S. Sabhapandit and S. N. Majumdar, Freezing transition in the barrier crossing rate of a diffusing particle, Phys. Rev. Lett. 125, 200601 (2020).
  * Zarfaty  _et al._ [2022] L. Zarfaty, E. Barkai, and D. A. Kessler, Discrete sampling of extreme events modifies their statistics, Phys. Rev. Lett. 129, 094101 (2022).
  * Kessler and Barkai [2010] D. A. Kessler and E. Barkai, Infinite covariant density for diffusion in logarithmic potentials and optical lattices, Phys. Rev. Lett. 105, 120602 (2010).
  * Dechant  _et al._ [2011] A. Dechant, E. Lutz, E. Barkai, and D. A. Kessler, Solution of the Fokker-Planck Equation with a logarithmic potential, J. Stat. Phys. 145, 1524–1545 (2011).
  * Hirschberg  _et al._ [2011] O. Hirschberg, D. Mukamel, and G. M. Schütz, Approach to equilibrium of diffusion in a logarithmic potential, Phys. Rev. E 84, 041111 (2011).
  * Hirschberg  _et al._ [2012] O. Hirschberg, D. Mukamel, and G. M. Schütz, Diffusion in a logarithmic potential: scaling and selection in the approach to equilibrium, Journal of Statistical Mechanics: Theory and Experiment 2012, P02001 (2012).
  * [42] The maximum contribution to the integral in Eq. (11) is centered at $x^{*}(t)$. In the long-time limit, $x^{*}(t)/x_{c}(t)\equiv z^{*}(t)/z_{c}(t)\sim t^{-\gamma\nu}\to 0$, and therefore, the important region is the $z\to 0$ region.
  * Stella  _et al._ [2023a] A. L. Stella, A. Chechkin, and G. Teza, Anomalous dynamical scaling determines universal critical singularities, Phys. Rev. Lett. 130, 207104 (2023a).
  * Stella  _et al._ [2023b] A. L. Stella, A. Chechkin, and G. Teza, Universal singularities of anomalous diffusion in the Richardson class, Phys. Rev. E 107, 054118 (2023b).


# Magnetocaloric effect for a $Q$-clock type system

Michel Aguilera<EMAIL_ADDRESS>Instituto de Física, Pontificia
Universidad Católica de Valparaíso, Casilla 4950, Valparaíso 2373223, Chile
Departamento de Física, Universidad Técnica Federico Santa María, Valparaíso
2390123, Chile  Sergio Pino-Alarcón  Departamento de Física, Universidad
Técnica Federico Santa María, Valparaíso 2390123, Chile  Francisco J. Peña
<EMAIL_ADDRESS>Departamento de Física, Universidad Técnica Federico
Santa María, Valparaíso 2390123, Chile  Millennium Nucleus in NanoBioPhysics
(NNBP), Av. España 1680, Valparaíso 2390123, Chile  Eugenio E. Vogel
Departamento de Ciencias Físicas, Universidad de La Frontera, Casilla 54-D,
Temuco, Chile  Center for the Development of Nanoscience and Nanotechnology,
Santiago 9170124, Chile  School of Engineering, Central University of Chile,
Santiago, Chile  Patricio Vargas  Departamento de Física, CEDENNA, Universidad
Técnica Federico Santa María, Av. España 1680, Valparaíso 2390123, Chile

###### Abstract

In this work, we study the magnetocaloric effect applied to a magnetic working
substance corresponding to a square lattice of spins with $Q$ possible
orientations known as the “$Q$-state clock model” where for $Q\geq 5$, the
systems present the famous Berezinskii–Kosterlitz–Thouless phase (BKT).
Thermodynamic quantities are obtained in exact form for a small lattice size
of $L\times L$ with $L=3$ and by the mean-field approximation and Monte Carlo
simulations for $Q$ pairs between 2 and 8 with $L=3,8,16,32$ with free
boundary conditions, and magnetic fields varying between $B=0$ and $1$ in
natural units of the system. By obtaining the entropy, it is possible to
quantify the caloric effect through an isothermal process in which the
external magnetic field on the spin system is varied. In particular, we find
the values of $Q$ that maximize the effect depending on the lattice size and
the magnetic phase transitions related to maximizing the caloric phenomena.
These indicate that in a small lattice (up to $\sim 7\times 7$), when $Q\geq
5$, the transition that maximizes the effect is related to ferromagnetic to
BKT type. In contrast, transitioning from BKT to paramagnetic type increases
the system’s caloric response when we work with a larger lattice size.

##  I Introduction

The caloric phenomena constitute a fundamental axis of study in the physics of
materials to find replacement candidates for the toxic gases used in
conventional refrigeration systems [1, 2]. To this end, considerable
temperature variations have been achieved for different types of materials
after an adiabatic variation of the value of the applied hydrostatic pressure
(barocaloric effect) [3, 4, 5, 6], mechanical stress (elastocaloric effect)
[7, 8, 9, 10, 11], electric field (electrocaloric effect) [12, 13, 14, 15] and
external magnetic field (magneto-caloric effect, MCE) [16, 17, 18, 19].
Although each phenomenon’s nature and physical interpretation are different,
they comply with a basic concept that characterizes them: the aim is to
maximize the variation of entropy in the face of changes in the control
parameter that governs the thermodynamic process that describes the phenomenon
[20, 21, 22, 23, 24, 25, 26]. In particular, the systems that present phase
transitions maximize entropy variation near the transition [27, 28, 29, 30].
One of many corresponds to magnetic materials, characterized by a determined
magnetic order or ground state. Remarkably, such an order depends on the
energy contributions in the material [31, 32, 33]. Within these transitions,
we will focus on those magnetic systems that present as the temperature
increases and thus change their magnetic phase. An interesting system to
explore corresponds to the $Q$-state clock model that constitutes the discrete
version of the famous 2D XY model [34, 35, 36, 37], which is probably the most
extensively studied example showing the Berezinskii–Kosterlitz–Thouless (BKT)
transition in the presence of a frustrated quenched disordered phase [38, 39,
40]. The $Q-$state clock model is one of many magnetic models to mimic the
thermodynamics of some materials, and it can be viewed as a classical
Heisenberg spins model with very strong planar anisotropy [41, 42, 43].

One way to characterize the phase transitions of the $Q$-clock state model is
through the maxima obtained in the specific heat as a function of temperature.
Each location of a maximum of the specific heat on the temperature axis will
represent a value for a so-called critical temperature. It has been shown [44,
45, 46, 47, 48, 49] (in the absence of an external magnetic field) that for
the $Q$-clock state model, values $Q\geq 5$ (where $Q$ represents the number
of possible orientations that the spins can take), the specific heat presents
two maxima. The first maximum corresponds to transitioning from a
ferromagnetic phase (FP) to a BKT phase. In contrast, the second maximum
transitions from BKT to a paramagnetic disordered phase (PP) [36].

In this sense, this model, which presents these two phase transitions
mentioned above, is very interesting to exploit in studying caloric phenomena.
In addition, applying this spin model to this effect raises several questions:
which value of $Q$ is the optimum, i.e., maximizes the thermodynamic response?
How does the estimation of the magnitude of the effect through mean-field
approximation compare to the exact calculations and Monte Carlo simulations?
Is it possible to know which magnetic phases maximize the effect? Is there
convergence in the results when we increase the size of the lattice?

One of the technical difficulties with this model is the scaling of the
accessible states, which scales as $\sim Q^{L^{2}}$ where $L\times L$
represents the size of the lattice you work. This means that when solving
thermodynamics exactly via a canonical partition function, the computational
costs result in long computation times. With this limitation, exact model
computations can only be applied to lattices of small size. That is why two
relevant mechanisms are employed to see what happens with larger lattices:
Mean field approximation and Monte Carlo simulations.

In this work, we study the MCE in the $Q$-clock model, in which the
thermodynamic quantities of the model are obtained exactly for small lattices
and using a mean-field approximation and Monte Carlo simulations for larger
lattice sizes. The $Q$ values studied correspond to the even values $Q=2,4,6$
and $8$. The exact calculation for a small lattice of $L\times L$ with $L=3$
is presented to review the MCE concerning the magnetic phase and to quantify
the caloric effect per spin of small-size lattices. In completeness, the
results are analyzed using the mean-field approximation of this lattice to
verify if the estimates in the orders of magnitude of the effect are
consistent with this approximation. Finally, our analysis is extended to
$16\times 16$ and $32\times 32$ lattices employing a mean-field and Monte
Carlo simulations.

Our work is structured as follows: Section II presents the $Q$-clock model,
Section III revisits the mean-field approximation, Section IV analyzes the
thermodynamic quantities involved in the calculations, Section V gives an
analysis of the caloric phenomena, Section VI shows the results, and finally,
Section VII presents our work’s conclusions.

##  II Spin Model

###  II.1 $Q$-state clock model

Figure 1: (a) $6\times 6$ spin lattice and (b) $Q$ clock model for $Q=2,4,6$
and $8$.

The working substance under study corresponds to the $Q$-states clock model on
a two-dimensional (2D) square lattice of dimensions $L\times L=N$, where the
local magnetic moment or “spin” $S_{i}$ at site $i$ can point in any of $Q$
directions in a given plane. $S_{i}$ is then a 2D vector, i.e.
$S_{i}=(\cos(\frac{2\pi}{q}k),\sin(\frac{2\pi}{q}k))$, where $k=0,1,...Q-1$,
with equal probability. The magnitude of $S_{i}$ is chosen as the unity.
The isotropic Hamiltonian for such a system can be written as:

$\mathcal{H}=-\sum_{\langle
i,j\rangle}\mathcal{J}(\vec{S}_{i}\cdot\vec{S}_{j})-\sum_{i}\vec{B}\cdot\vec{S}_{i}\,,$
(1)

where $\mathcal{J}>0$ is the ferromagnetic exchange interaction to nearest
neighbors; the sum runs over all pairs of nearest neighbors $(i,j)$, which is
indicated by the symbol $\langle i,j\rangle$ under the summation symbol.
$\vec{B}$ is an external field applied along one direction in the plane. In
this work, we used arbitrary units choosing $\mathcal{J}=1$, making all the
calculated quantities for the exchange energy constant. From this Hamiltonian
model, we must extract the accessible states of the system to obtain the
entropy of the system.

To quantify the caloric effect, it is necessary to obtain the entropy from the
thermodynamic study of the model. We will first analyze two cases where this
calculation comes from a canonical partition function. The first corresponds
to an exact calculation of all the accessible microstates of the system for a
$3\times 3$ lattice for $Q=2,4,6$ and $8$. For the same values of $Q$ studied,
the mean-field theory will be employed for lattices of size up to $32\times
32$. On the other hand, Monte Carlo simulations offer us the possibility to
obtain thermal averages with a sample of the configurations and to define an
entropy via integration of the specific heat.

In the next section, we will review the mean-field approximation to emphasize
the calculation of the approximate partition function and the definition taken
from there for the entropy coming from the Helmholtz free energy temperature
derivative, which also applies to the exact case in small lattices.

##  III Mean Field Approximation

The mean-field theory assumes that the fluctuations around the average value
of the order parameter, in this case, the magnetization per site
$\vec{m}=\vec{M}/N$, are so small that they can be neglected. The first term
of the Hamiltonian of Eq. (1) that corresponds to the interaction term between
the spin of the lattice in different sites is modified by performing the
following approximations.

We can write the spin term as follows:

$\vec{S}_{j}=\vec{m}+\delta\vec{S}_{j},$ (2)

where $\vec{m}$ is the average thermodynamic spin, the same for all sites in
the lattice. Therefore, we have

$\delta\vec{S}_{j}=\vec{S}_{j}-\vec{m}_{j}.$ (3)

Thus, the spin-spin interaction term can be written as

$\vec{S}_{i}\cdot\vec{S}_{j}=-m^{2}+\vec{m}\cdot\left(\vec{S}_{i}+\vec{S}_{j}\right),$
(4)

where we have neglected the square terms of the fluctuation
($O(\delta\vec{S})^{2}$). Therefore, the interaction term of the Hamiltonian
of Eq. (1) (that we call $\mathcal{H}_{J}$) can take the form

$\displaystyle\mathcal{H}_{\mathcal{J}}$ $\displaystyle=$
$\displaystyle-\mathcal{J}\sum_{\langle
i,j\rangle}\left(-m^{2}+\vec{m}\cdot\left(\vec{S}_{i}+\vec{S}_{j}\right)\right)$
(5) $\displaystyle=$
$\displaystyle\sum_{i}(2\mathcal{J}m^{2}-\mathcal{J}z\vec{m}\cdot\vec{S}_{i}),$

where now the sum runs for each site in the lattice, and $z$ are the effective
nearest neighbors of the model. Consequently, we can define a Hamiltonian per
site given by the structure

$h_{i}=\frac{z}{2}\mathcal{J}m^{2}-\mathcal{J}z\vec{m}_{i}\cdot\vec{S}_{i}-\vec{B}\cdot\vec{S}_{i}.$
(6)

Now, we can calculate the partition function per site, which will depend on
$Q$, $B$, $m$ and $T$ given by

$Z(q,B,m,T)=\sum_{q}e^{-\frac{\xi(q,B,m)}{T}},$ (7)

where $\xi(q,B,m)$ is the energy per site coming from the Hamiltonian of Eq.
(6).

To extrapolate the results of mean field approximation of the thermodynamic
quantities in smaller lattices, in the Ref. [43], we discuss the idea of
finding an effective nearest neighbor ($z_{eff}$) through a simple
optimization protocol. This protocol matches the internal energies obtained
for a lattice of $L\times L$ with $L=3$ employing exact and mean-field
calculations. In this way, by releasing the number of neighbors,
$z_{eff}\in\,\mathbb{R}^{+}$, and minimizing the internal energy difference
(between the exact and approximate case via mean-field), it was found that for
$L=3$, the optimal number of neighbors was $z_{eff}=2.67$ (for all values of
$Q$).

Amplifying the above qualitatively, we propose an expression for the number of
near neighbors effective that adjusts according to the weighting of the effect
of non-interacting edges in the system when the lattice has a generic
resolution $L\times L$. The effective neighbors for a central spin for the
mean-field are determined by the following expression (independent of $Q$) for
the square lattice with up, down, left, and right near neighbors [43]

$z_{eff}=\frac{4\times(L-2)^{2}+3\times 4\times(L-2)+2\times 4}{L^{2}}.$ (8)

##  IV Thermodynamic Quantities

Once the partition function of the system has been obtained by either of the
two approaches discussed above, it is possible to compute all the
thermodynamical observables in a general way through the expressions ( with
$k_{B}=1$ )

$F=-T\ln Z,$ (9)

$U=T^{2}\frac{\partial\ln Z}{\partial T},$ (10)

and

$C=\frac{\partial U}{\partial T},$ (11)

where $F$ is the Helmholtz free energy, $U$ is the internal energy and
$\mathcal{C}$ is the specific heat at constant magnetic field. In addition,
with the differential expression of Helmholtz free energy given by $dF=-SdT-
MdB$, we can obtain the entropy and the magnetization of the system given by

$S=-\frac{\partial F}{\partial T};\quad M=-\frac{\partial F}{\partial B}.$
(12)

Figure 2: Specific heat as a function of temperature for even $Q$ values
between $Q=2$ and $Q=8$ for (a) $B=0.05$ and (b) $B=0.1$ for a $3\times 3$
lattice. We observe the peak displacements to the right at a higher magnetic
field.

We will begin by reviewing the thermodynamic results of the model, focusing
first on the analysis of the specific heat for the exact results of a $3\times
3$ lattice. This thermodynamic quantity for a field $B=0.05$ and field $B=1$
is shown in Fig. 2(a) and Fig. 2(b), respectively for $Q=2,4,6,8$. Here, we
observe two peaks for $q\geq 6$, indicative of a double-phase transition. It
is, therefore, the BKT phase that appears in these cases. It is also
noticeable that the critical temperature of each transition increases with the
increasing external field, as we can observe from Fig. 2(b) compared with the
critical temperatures obtained for (a) panel of the same figure. This is
because the external field favors ordered phases (FM and BKT) over disordered
ones. A different behavior is obtained by fixing the field and increasing $Q$,
where we observe that the high-temperature peak remains at the same
temperature, whereas the low-temperature peak tends toward lower temperatures
(eventually to a zero temperature) as $Q$ increases, tending to infinity.
Therefore, for $Q\rightarrow\infty$, we expect the temperature of the first
transition (from FM to BKT phase) to go to zero. The case of internal energy
is shown in Fig. 3(a), where we observe that we obtain lower values for the
case with a magnetic field than for the case without an external magnetic
field. This is due to the Zeeman term present in the Hamiltonian given by Eq.
(1), which decreases the ground state energy compared to the case in the
absence of a magnetic field. Continuing with the thermodynamic description of
our system of study, Fig. 3(b) shows the magnetization of the system for
$Q=2,4,6,$ and $8$ for high fields $B=1$ (blue curves) and $B=0.05$ (purple
curves). First, we notice that all systems start saturated (in a ferromagnetic
state) and then decrease their magnetization as the temperature increases. The
case with a higher external magnetic field shows higher magnetization than the
case with a lower magnetic field as the temperature increases for all $Q$
displayed. It is observed that the change in magnetization as a function of
temperature happens much faster for larger and larger q values. This is a
direct consequence of the degrees of freedom of the system. For example, we
can think of $Q=2$, where two possible orientations are up and down. In terms
of energy, generating a flip of such a configuration is much more difficult.
Therefore, the magnetization at the same temperature will be higher than in
the other cases for higher values of $Q$.

Figure 3: Internal energy (a), magnetization (b), and entropy (c) as a
function of temperature for different values of $Q$ pairs between $Q=2$ and
$Q=8$ for fields $B=0.05$ (purple) and $B=0.1$ (blue) for a $3\times 3$
lattice.

Next, we analyze the entropy in Fig. 3(c). We observe that for a given $Q$,
the entropy is lower for the case with field $B=1$ than for the case of
$B=0.05$. This is consistent with the interpretation that the magnetic field
tends to order the spin system. Also, it is noticeable that for more spin
degrees of freedom (as $Q$ increases), the entropy is higher as a function of
temperature consistent with the greater number of accessible states the system
acquires. These graphs show that the entropy fulfills that $\lim_{T\rightarrow
0}S(T,B\neq 0)=0$ because, with even a very low field value, the system at low
temperature tends to have only one preferred state due to the presence of the
external field applied to the system. It is important to remember that in the
case of a null field, the entropy at low temperature has the value of
$\lim_{T\rightarrow 0}S(T,0)=\ln Q$.

Figure 4: In panel (a), we show the spin magnetization of a $3\times 3$
lattice with an external magnetic field equal to unity for different values of
$Q=2,4,6,8$ calculated exactly (orange color) and approximated via mean field
(blue color). In panel (b), we exhibit the spin magnetization of a $3\times 3$
lattice with an external magnetic field close to zero ($B=0.05$ in
$\mathcal{J}$ units) for different values of $Q=2,4,6,8$ calculated exactly
(purple color) and approximated via mean field (green color).

In Fig. 4(a), we note the saturation in the magnetization for low
temperatures, finding for all the cases studied the unity for the limit of
$T\rightarrow 0$. Additionally, we observe that the temperature increase
increases the difference between the mean field and the exact calculation for
the Ising case $(Q=2)$. In contrast, for $Q\geq 4$, the difference between
both methods decreases as the temperature increases. In Fig. 4(b), we find the
spin magnetization of a $3\times 3$ lattice with an external magnetic field
close to zero ($B=0.05$ in $\mathcal{J}$ units) for different values of
$Q=2,4,6,8$ calculated in an exact way (purple color) and approximated via
mean field (green color). We note a difference in the shape of the curves as
the temperature increases, although both methods start to become saturated at
low temperatures. The shape of the exact calculated curve shows that the
magnetization behavior for low non-zero external fields varies with increasing
field and not so for the medium field case.

##  V Caloric Phenomena

To analyze this effect, we start by treating the total entropy of a material
as the sum of three entropies: electronic $S_{e}$, magnetic $S_{m}$, and
lattice $S_{l}$.

$S(T,B)=S_{e}(T,B)+S_{m}(T,B)+S_{l}(T,B)$ (13)

Magnetic entropy depends strongly on the magnetic field, whereas it is very
common to find materials with electronic and lattice entropy with negligible
dependence on the magnetic field. However, this is not general, and we
highlight that in a low temperatures regime ($\approx<$ 10 K), some reports
indicate that the electronic entropy can present nonlinear dependence on the
magnetic field. In this work, let us assume that we are dealing with a more
usual system, where we can establish null dependence of the magnetic field in
the lattice and electronic entropies so that Eq. (13) is of the form:

$S(T,B)\approx S_{e}(T)+S_{m}(T,B)+S_{l}(T)$ (14)

The quantification of the caloric phenomenon is associated with a
thermodynamic refrigerator cycle, where we can take two paths to quantify the
effect: i) an adiabatic trajectory and ii) an isothermal trajectory. The
temperature variation that suffers the systems along the isentropic process is
called $\Delta T_{ad}$ and is given by the expression

$\Delta T_{ad}=-\int_{B_{i}}^{B_{f}}\frac{T}{C_{B}}\left(\frac{\partial
S}{\partial B}\right)_{T}dB,$ (15)

where we use that $C_{B}=\left(\frac{\partial U}{\partial
T}\right)_{B}=T\left(\frac{\partial S}{\partial T}\right)_{B}$ corresponds to
specific heat at constant $B$ parameter.

In the case of quantifying the effect employing an isothermal trajectory, we
obtain an entropy variation at constant temperature ($\Delta S$) given by

$\Delta S=\int_{B_{1}}^{B_{2}}\left(\frac{\partial S}{\partial
B}\right)_{T}dB.$ (16)

The quantity $\left(\frac{\partial S}{\partial B}\right)_{T}$ that appears in
the last two equations can be rewritten in terms of the magnetization of the
system via Maxwell’s relationship $\left(\frac{\partial S}{\partial
B}\right)_{T}=\left(\frac{\partial M}{\partial T}\right)_{B}$.

If we look at Eq. (15) and Eq. (16), we can find a relationship between these
quantities concerning each other. It is found that $-\Delta S\propto\Delta
T_{ad}$. Consequently, it is essential to note that when you have a case in
which $-\Delta S>0$, let’s call this kind of response of direct type. The
system will heat up, while when the response is of the $-\Delta S<0$ type, we
call this inverse type response, and consequently, the system will cool down.
Therefore, we expect for direct response a $\Delta T_{ad}>0$, and for the case
of an inverse response, we expect a $\Delta T_{ad}<0$ for the final result of
the caloric phenomena.

For the case of the entropy variation at constant temperature ($T$), the
expression can be given by the difference of the entropy at the initial and
final point of the process in the following form:

$-\Delta S=-\Delta S_{m}=S_{m}(T,B_{0})-S_{m}(T,B),$ (17)

where the contributions of $S_{l}$ and $S_{e}$ have been canceled due to its
dependence only on temperature, and consequently, in an isothermal process, it
will have no associated variations. The last equation generates a graph for a
given $B_{0}$ and final $B$ as a function of $T_{0}$ that quantifies the
caloric phenomena.

##  VI Results

Figure 5: (a) Entropy as a function of temperature for the case $Q=2$ for the
case of $B=0.05$ (purple curve), $B=0.5$ (green curve) and $B=1$ (blue curve).
(b) Entropy difference $-\Delta S=S(T,B_{0}=0.05)-S(T,B)$ as a function of
temperature for $B=0.5$ (green curve), $B=1$ (blue curve) and $B=0.05$ (purple
curve).

To start the discussion of the caloric phenomena, we must fix the zero value
for temperature and the external magnetic field for our calculations. The zero
temperature for our thermodynamic quantities will be $T=0.01$ while the
smallest field we will work with will be $B=0.05$. These two quantities
combined produce zero entropy, as theoretically expected. In addition, the
selection of $B=0.05$ as the lowest field is supported to avoid the problem
generated when we have $B\rightarrow 0$ in the mean-field self-consistency
equation. Regarding the maxima used in our calculations, the maximum
temperature corresponds to $T=3.01$ and $B=1$ in units of $\mathcal{J}$, both
quantities.

###  VI.1 Small Lattices

Let us start the analysis with a simple example of the exact case ($3\times 3$
lattice) for $Q=2$ to understand the caloric response of the system. For this
purpose, we plotted the entropy for three different field values of the system
in Fig. 5(a), where $B=0.05$ is the purple curve, $B=0.5$ is the green curve,
and B=1 is the blue curve. The entropy difference that quantifies the effect
is given by $-\Delta S=S(T,0.05)-S(T,B)$, where we have selected the initial
point with an external field strength $B=0.05$. These entropy differences can
be plotted for different field values on a curve as a function of the
temperature to which the process is taken, and they are presented in Fig.
5(b). From this last figure, we note that the response of the effect is
positive, which indicates a direct type of behavior for both cases presented.
It is also observed that below 0.9 temperature (approximately), the thermal
response of the effect for $B=0.5$ and $B=1$ are equal. After this temperature
value, it is observed that the case with the higher field is notably different
from the one with a lower field. The explanation of this phenomenon is
relatively simple. The system prefers to be in a ferromagnetic configuration
in an external magnetic field, i.e., in a single configuration state. To
remove it from this state requires a considerable increase in temperature in
terms of energy. This is reflected in the entropy for cases $B=0.5$ (green
curve) and $B=1$ (blue curve) in Fig. 5(a) where both entropies are equal up
to $T\sim 0.9$ and then differentiate with $B=0.5$ being higher than $B=1.0$.

Figure 6: Caloric response with an initial $B$ of $0.05$ for the case of a
final field of (a) $B=0.1$ and (b) $B=1$ for exact calculations for a $3\times
3$ lattice and using the mean-field approximation.

To compare the effects for the exact calculations and mean field, we have
analyzed $-\Delta S$ for the case of a final field of $B=0.1$ and $B=1$ and
the initial fixed value of field $B=0.05$. From Fig. 6(a), we note a large
discrepancy in the effect’s location as a temperature function for the case
with a final $B$ of $0.1$. As we observe in Fig. 6(b), this discrepancy
decreases as we increase the final field. It is remarkable from these two
comparative plots that if we only focus on estimating the maximum effect at
the value associated only with $-\Delta S$ and not the location of the
temperature, both calculations do not differ too much. We will return to this
point later in the discussion when we deal with other $Q$ cases.

Figure 7: Caloric response for a $3\times 3$ lattice with an initial field of
B=0.05 and a final field of (a) $B=0.1$ and (b) $B=1$ for the exact and mean
field cases for $Q=2,4,6$ and $8$.

We will now analyze the MCE for different values of $Q$ of the $3\times 3$
lattice, both in the exact case and using the mean field. This is presented in
Fig. 7, where the cases for a final field of $B=0.1$ (panel (a)) and the case
of $B=1$ (panel (b)) with an initial field of $B=0.05$ are shown. We first
notice that for the value of $Q=4$, what was observed for $Q=2$, associated
with the location of the maximum temperature effect for both methods, is also
true. As the magnetic field increases, the maximum temperature locations
become closer. However, a considerable error is obtained when extrapolating
mean-field solutions for small lattices, and it has to do with the location of
the maximum (or maxima) at $-\Delta S$ when we have the cases $Q=6$ and $Q=8$
as we can appreciate in Fig. 7(b). From this last figure, it can be seen that
both $Q=6$ and $Q=8$ present two maxima in $-\Delta S$, vestiges of the double
transition present for those values of $Q$ in the model, both for the mean
field and the exact calculations. We note that the exact calculations indicate
that the peak at the lowest temperature is where the greatest effect occurs,
while the mean field indicates the opposite; the second peak indicates the
maximum value of the effect. Our previous work [41, 42] shows this phenomenon
in the specific heat up to lattice of approximately $7\times 7$. After that
value, the second peak observed in the specific heat becomes larger, following
the trend shown and observed in the mean-field approximation. In addition,
from Fig. 7(a), we observe that in the exact case for the same parameter and
lattice configuration, the maximum of the effect is slightly higher for $Q=6$
and $8$ than for $Q=4$. Meanwhile, the mean-field indicates that the ordering
of the effect from highest to lowest is $Q=4$, $Q=6$, $Q=8$, and finally,
$Q=2$. In contrast, from Fig. 7(b), we can see that the largest effect is
obtained for $Q=4$ independent of the calculation method. This is remarkable
since previous studies found that when applying thermodynamic cycles using the
$Q$-clock model as the working substance, the maximum efficiency and useful
work are obtained for the parameter $Q=4$ [43]. These two cases, therefore,
indicate that there is a field value at which this inversion observed in the
effect maxima occurs. For this purpose, we plot in Fig. 8 the value of the
maximum $-\Delta S_{max}$ as a function of the field independent of the
location in temperature at which it occurs both for the mean field and for the
exact calculation. From this plot, we can notice that in the exact case, $Q=6$
and $Q=8$, responses higher than $Q=4$ and $Q=2$ are observed up to a field of
approximately $\sim 0.26$. After that value, the effect is ordered in the same
way as the mean-field results, i.e., in the form:

$\displaystyle-\Delta S_{max}^{Q=4}>-\Delta S_{max}^{Q=6}>-\Delta
S_{max}^{Q=8}>-\Delta S_{max}^{Q=2},$ (18)

equation valid only for a final $B>0.26$.

Figure 8: Maximum caloric response ($-\Delta S_{max}$) as a function of the
applied external field between $0.1$ and $1.0$ for a $3\times 3$ lattice with
an initial magnetic field $B=0.05$ for $Q=2,4,6$ and $8$ employing exact
calculations (solid lines) and mean field calculations (dashed lines).

From Fig. 8, we can also confirm what happens when estimating the maximum
effect between the exact and mean-field calculations. Although the mean-field
and exact differ quite a bit in the temperature location of the maximum point,
if we focus only on the numerical value of effect maximization at $-\Delta S$,
the mean-field does not show significant variations compared to the exact
calculation even for a small lattice, as shown in this study. The dotted
curves in Fig. 8 representing mean-field calculations differ very little from
the solid lines of the exact case, with the largest difference observed for
$Q=4$ and the smallest differences for $Q=6$ and $Q=8$. From this tendency, we
can conclude that for values larger than $Q=4$, the numerical estimate of
$-\Delta S$ will be closer between the two methods presented. We can say,
then, that the mean field is useful in this context even for small lattices in
the numerical estimation of the maximum effect but fails logically in the
location of the temperature where it occurs because the nature of the
approximation used is regularly employed for lattices of larger sizes.

Figure 9: $T$ vs. $B$ phase diagram for (a) $Q=2$ and $Q=4$ and for (b) $Q=6$
and $Q=8$. The solid lines for (a) represent the critical temperatures of the
FP-PP type transition of the final state of the spin system $S(T,B)$, while
for (b), the solid lines represent the critical temperatures of the FP-BKT
type transition of the final state of the spin system $S(T,B)$. For (a), the
light blue triangles indicate the maximum obtained from $-\Delta
S=S(T,B_{0}=0.05)-S(T,B)$ with $Q=2$ and the green squares indicate the same
but for $Q=4$. The same applies to $Q=6$ represented by magenta circles and
$Q=8$ with blue circles. The horizontal dotted lines indicate the critical
temperature of the $S(T,B_{0}=0.05)$ state for $Q=2$ (light blue), $Q=4$
(green), $Q=6$ (magenta) and $Q=8$ (blue). The black vertical dotted lines(for
panels (a) and (b)) represent the location where the systems maximize $-\Delta
S$ passing through an effective phase change.

An interesting point to discuss is the location in temperature where the
maximum value of the effect is given compared to the critical temperature
$(T_{c})$ of the different phases present in the model. This is shown in Fig.
9, which represents a plot of $T$ vs. $B$ for (a) $Q=2$ and $Q=4$ and for (b)
$Q=6$ and $Q=8$. To understand these plots, we will start with the cases $Q=2$
and $Q=4$ presented in Fig. 9(a). We must think that the variation of the
effect is given by $-\Delta S=S(T,B_{0}=0.05)-S(T,B)$ and hence is described
by a spin system starting with an entropy of the form $S(T,B_{0}=0.05)$. This
quantity has associated a specific heat $C(T,B_{0}={0.05})$ that defines a
critical temperature represented by horizontal lines in Fig 9(a). On the other
hand, the continuous lines represent the final entropy state $S(T,B)$, which
has an associated specific heat $C(T,B)$, which will deliver as the field
changes different critical temperature values, indicating the respective phase
transition for these cases. The graphs with geometrical figures (triangular
for $Q=2$ and square for $Q=4$) are the points where, for a given field,
$-\Delta S$ is the maximum. Precisely, the intersection between the curve with
geometric figures and the dotted horizontal one generates, for $Q=2$ and
$Q=4$, two possible regions where the system maximizes its effect. Under a
field value given by the vertical black line, the system maximizes its effect
in a region where the system starts in a ferromagnetic phase $(FP_{B_{0}})$
and ends in a ferromagnetic phase ($FP_{B}$) while past that critical field
value found the system maximizes the caloric response starting from a
paramagnetic phase ($PP_{B_{0}}$) and ending in a ferromagnetic phase
($FP_{B}$).

The same analysis can be made if we recall that for the exact case of $Q=6$
and $Q=8$, the peak of the effect was concentrated in the first peak of Fig. 7
and, therefore, more oriented to transitions relating to the BKT and FP
phases. This is shown in Fig. 9(b), where the maxima are concentrated when the
system remains in a ferromagnetic phase or involves transitions linking the
BKT and FP phases. We note that the maxima of the effect for $Q=8$ are much
closer to the $T_{c}$ of the first phase transition related to an FP-BKT type
transition of the final state of the spin system than what occurs for the case
of $Q=6$.

###  VI.2 Large Lattices

Including mean-field calculations allows for studying the caloric response for
larger lattice sizes and analyzing the convergence of results. Fig. 10(a)
presents the spin caloric response as a function of temperature for $Q=4$, and
in Fig.(b) for the case of $Q=8$ both for lattice sizes of $3\times 3$,
$8\times 8$, $16\times 16$ and $32\times 32$ with an initial field of $B=0.05$
and a final field of $B=1$. As discussed, the $Q=4$ case only presents one
transition of type FP to PP, and we therefore see only a peak in the entropy
difference, while $Q=8$ presents two phase transitions where we see that
$-\Delta S$ is larger when we go from BKT to PP over that from FP to BKT. From
both graphs, we notice that there is a convergence of the caloric response
already for a $32\times 32$ lattice because it presents the same results as a
$16\times 16$ lattice both for the case of $Q=4$ and $Q=8$. It can also be
observed by comparing the two graphs that it is still $Q=4$ that presents the
highest $-\Delta S/N$ evaluated on the system.

Figure 10: Mean field calculations for (a) $Q=4$ and (b) $Q=8$ for $3\times
3$, $8\times 8$, $16\times 16$ and $32\times 32$ lattice.

Figure 11: Magneto-caloric per spin effect for $16\times 16$ and $32\times 32$
lattices using mean-field calculations (orange curves) and using the Monte
Carlo method (black-grey curves) for the case (a) $Q=4$ and (b) $Q=8$.

For the completeness of our calculations, we performed a large lattice scan
where the thermodynamic properties of the magnetic substance are obtained
through the Monte Carlo method and state selection by the Metropolis
algorithm, thus obtaining the internal energy and magnetization for fixed
external fields. Fixing the external magnetic field of the simulation at
$0.05$ and $1.00$ (in units of J) and calculating the observables with a
thermalization of $50000$ flips and a set of $100000$ iterations for the
calculation of thermal averages, the curve of the magnetocaloric effect is
obtained for $L=16$ and $L=32$. In Fig. 11, we observe the magneto-caloric per
spin effect for $16\times 16$ and $32\times 32$ lattices using mean-field
calculations (orange curves) and using the Monte Carlo method (black-grey
curves) for the case (a) $Q=4$ and (b) $Q=8$. We can observe from there that
the Monte Carlo calculation concerning the mean-field calculation coincides
with each other in the maximum value of the effect but not in the temperature
at which the maximum occurs. It can also be seen that for the case of the
Monte Carlo analysis, the $16\times 16$ and $32\times 32$ effects still show
differences with each other (small), while, as we mentioned earlier in the
discussions, the mean-field case already converges for a $32\times 32$
lattice. However, it is observed that for Monte Carlo simulations up to the
location of the first maximum (or the only one for $Q<5$) corresponding to an
FP to PP transition, $-\Delta S$ presents almost the same values for $16\times
16$ and $32\times 32$. However, there is a small difference between the values
obtained for these networks near the single transition for Q=4 and from the
second transition for $Q=8$.

Finally, returning to the discussion of the phases and their role in
quantifying the effect, Fig. 11(a)-(b) shows that for large networks, the peak
of the effect is much closer to the BKT-PP type transition, and this is
corroborated by both mean-field and Monte Carlo simulations. However, we must
reiterate that the opposite occurs in small lattices, where exact calculations
have been performed since we observe that the maximum $-\Delta S$ is close to
the FP-BKT type transitions.

##  VII Conclusions

In this work, we have studied the magneto-caloric effect for a spin system
with $Q$ orientations known as the $Q$-clock model. We reviewed the
thermodynamic quantities describing the model and calculated for $Q$ pairs
between $2$ and $8$ using three methods: exact, mean-field, and Monte Carlo
simulations. For small lattices ($3\times 3$) by exact calculations, we have
found that there is a region of field where the effects of $Q=6$ and $Q=8$ are
higher than those presented for the cases of $Q=4$ and $Q=2$, and then in
higher regimes of the external field, they are ordered in descending form as
$-\Delta S_{max}^{Q=4}>-\Delta S_{max}^{Q=6}>-\Delta S_{max}^{Q=8}>-\Delta
S_{max}^{Q=2}$. Forcing the mean-field approximation to explore these lattice
sizes, we note that mean-field allows us to quantify the maximum (in
intensity) fairly close to the exact calculation; however, for reasons of
applying it to small lattices, it fails considerably in the temperature
location where it occurs.

Related to the magnetic phases involved in maximizing the caloric phenomenon,
we can report that for small lattices for $Q\geq 5$, the FP-BKT transition
plays a fundamental role. This changes abruptly when we explore large
lattices, where Monte Carlo simulations are employed, which indicate that the
transition where the effect is maximized is closer to a BKT-PP type. In
addition, there are two considerable effects for large-size lattices: First,
there is a convergence of results, evidenced by Monte Carlo and mean-field
simulations for lattices as large as $32\times 32$ in the results, and the
superior effect always occurs for $Q=4$.

##  VIII Acknowledgments

M. A., F. J. P., E. V., and PV acknowledge the financial support of CEDENNA
(Chile) under contract ANID AFB220001 and Fondecyt (Chile) under contract
1230055. F. J. P. and P.V. acknowledge support from ANID Fondecyt grant under
contract 1240582. F.J.P. acknowledges financial support from “Millennium
Nucleus in NanoBioPhysics” project NNBP NCN2021 _021 and UTFSM DGIIE.

## References

  * Reis and Ma [2020] M. S. Reis and N. Ma, Caloric effects of quantum materials: An outlook, Physics Open 4, 100028 (2020).
  * Ram  _et al._ [2018] N. R. Ram, M. Prakash, U. Naresh, N. S. Kumar, T. S. Sarmash, T. Subbarao, R. J. Kumar, G. R. Kumar, and K. C. B. Naidu, Review on magnetocaloric effect and materials, Journal of Superconductivity and Novel Magnetism 31, 1971 (2018).
  * Cirillo  _et al._ [2022] L. Cirillo, A. Greco, and C. Masselli, Cooling through barocaloric effect: A review of the state of the art up to 2022, Thermal Science and Engineering Progress 33, 101380 (2022).
  * Li  _et al._ [2019] B. Li, Y. Kawakita, S. Ohira-Kawamura, T. Sugahara, H. Wang, J. Wang, Y. Chen, S. I. Kawaguchi, S. Kawaguchi, K. Ohara, _et al._ , Colossal barocaloric effects in plastic crystals, Nature 567, 506 (2019).
  * Matsunami  _et al._ [2015] D. Matsunami, A. Fujita, K. Takenaka, and M. Kano, Giant barocaloric effect enhanced by the frustration of the antiferromagnetic phase in mn 3 gan, Nature materials 14, 73 (2015).
  * Ma and Reis [2017] N. Ma and M. S. Reis, Barocaloric effect on graphene, Scientific reports 7, 13257 (2017).
  * Xiao  _et al._ [2015] F. Xiao, M. Jin, J. Liu, and X. Jin, Elastocaloric effect in ni50fe19ga27co4 single crystals, Acta Materialia 96, 292 (2015).
  * Chen  _et al._ [2021] J. Chen, L. Lei, and G. Fang, Elastocaloric cooling of shape memory alloys: A review, Materials Today Communications 28, 102706 (2021).
  * Cong  _et al._ [2019] D. Cong, W. Xiong, A. Planes, Y. Ren, L. Mañosa, P. Cao, Z. Nie, X. Sun, Z. Yang, X. Hong, _et al._ , Colossal elastocaloric effect in ferroelastic ni-mn-ti alloys, Physical review letters 122, 255703 (2019).
  * Chauhan  _et al._ [2015a] A. Chauhan, S. Patel, and R. Vaish, Elastocaloric effect in ferroelectric ceramics, Applied Physics Letters 106 (2015a).
  * Chauhan  _et al._ [2015b] A. Chauhan, S. Patel, R. Vaish, and C. R. Bowen, A review and analysis of the elasto-caloric effect for solid-state refrigeration devices: Challenges and opportunities, MRS Energy & Sustainability 2, E16 (2015b).
  * Mischenko  _et al._ [2006] A. Mischenko, Q. Zhang, J. Scott, R. Whatmore, and N. Mathur, Giant electrocaloric effect in thin-film pbzr0. 95ti0. 05o3, Science 311, 1270 (2006).
  * Liu  _et al._ [2016] Y. Liu, J. F. Scott, and B. Dkhil, Direct and indirect measurements on electrocaloric effect: Recent developments and perspectives, Applied Physics Reviews 3 (2016).
  * Correia and Zhang [2013] T. Correia and Q. Zhang, Electrocaloric effect: an introduction, in  _Electrocaloric Materials: New Generation of Coolers_ (Springer, 2013) pp. 1–15.
  * Jia and Sungtaek Ju [2012] Y. Jia and Y. Sungtaek Ju, A solid-state refrigerator based on the electrocaloric effect, Applied Physics Letters 100 (2012).
  * Tishin and Spichkin [2016] A. M. Tishin and Y. I. Spichkin,  _The magnetocaloric effect and its applications_ (CRC Press, 2016).
  * De Oliveira and von RANKE [2010] N. De Oliveira and P. J. von RANKE, Theoretical aspects of the magnetocaloric effect, Physics Reports 489, 89 (2010).
  * Pecharsky and Gschneidner Jr [1999] V. K. Pecharsky and K. A. Gschneidner Jr, Magnetocaloric effect and magnetic refrigeration, Journal of magnetism and magnetic materials 200, 44 (1999).
  * Franco  _et al._ [2018] V. Franco, J. Blázquez, J. Ipus, J. Law, L. Moreno-Ramírez, and A. Conde, Magnetocaloric effect: From materials research to refrigeration devices, Progress in Materials Science 93, 112 (2018).
  * Reis [2020] M. S. Reis, Magnetocaloric and barocaloric effects of metal complexes for solid state cooling: Review, trends and perspectives, Coordination Chemistry Reviews 417, 213357 (2020).
  * Reis [2012] M. Reis, Oscillating magnetocaloric effect on graphenes, Applied Physics Letters 101 (2012).
  * Reis [2011] M. S. Reis, Oscillating magnetocaloric effect, Applied Physics Letters 99 (2011).
  * Reis [2014] M. Reis, Diamagnetic magnetocaloric effect due to a transversal oscillating magnetic field, Physics Letters A 378, 1903 (2014).
  * Martínez-Rojas  _et al._ [2023] P. Martínez-Rojas, M. E. Benavides-Vergara, F. J. Peña, and P. Vargas, Caloric effect due to the aharonov–bohm flux in an antidot, Nanomaterials 13, 2714 (2023).
  * Negrete  _et al._ [2018a] O. A. Negrete, F. J. Peña, and P. Vargas, Magnetocaloric effect in an antidot: the effect of the aharonov-bohm flux and antidot radius, Entropy 20, 888 (2018a).
  * Negrete  _et al._ [2018b] O. A. Negrete, F. J. Peña, J. M. Florez, and P. Vargas, Magnetocaloric effect in non-interactive electron systems:“the landau problem” and its extension to quantum dots, Entropy 20, 557 (2018b).
  * Tegus  _et al._ [2002] O. Tegus, E. Brück, L. Zhang, K. Buschow, F. De Boer, _et al._ , Magnetic-phase transitions and magnetocaloric effects, Physica B: Condensed Matter 319, 174 (2002).
  * Tishin  _et al._ [1999] A. Tishin, K. Gschneidner Jr, and V. Pecharsky, Magnetocaloric effect and heat capacity in the phase-transition region, Physical Review B 59, 503 (1999).
  * De Oliveira and Von Ranke [2008] N. De Oliveira and P. Von Ranke, Magnetocaloric effect around a magnetic phase transition, Physical Review B 77, 214439 (2008).
  * Franco  _et al._ [2009] V. Franco, A. Conde, M. Kuz’Min, and J. Romero-Enrique, The magnetocaloric effect in materials with a second order phase transition: Are tc and tpeak necessarily coincident?, Journal of Applied Physics 105 (2009).
  * Halder  _et al._ [2010] M. Halder, S. Yusuf, M. Mukadam, and K. Shashikala, Magnetocaloric effect and critical behavior near the paramagnetic to ferrimagnetic phase transition temperature in tbco 2- x fe x, Physical Review B 81, 174402 (2010).
  * Shamberger and Ohuchi [2009] P. J. Shamberger and F. Ohuchi, Hysteresis of the martensitic phase transition in magnetocaloric-effect ni-mn-sn alloys, Physical Review B 79, 144407 (2009).
  * Zhang  _et al._ [2008] C. Zhang, D. Wang, Q. Cao, Z. Han, H. Xuan, and Y. Du, Magnetostructural phase transition and magnetocaloric effect in off-stoichiometric mn1. 9- xnixge alloys, Applied Physics Letters 93 (2008).
  * Tobochnik [1982] J. Tobochnik, Properties of the q-state clock model for q= 4, 5, and 6, Physical Review B 26, 6201 (1982).
  * Li  _et al._ [2020] Z.-Q. Li, L.-P. Yang, Z. Xie, H.-H. Tu, H.-J. Liao, and T. Xiang, Critical properties of the two-dimensional q-state clock model, Physical Review E 101, 060105 (2020).
  * Miyajima  _et al._ [2021] Y. Miyajima, Y. Murata, Y. Tanaka, and M. Mochizuki, Machine learning detection of berezinskii-kosterlitz-thouless transitions in q-state clock models, Physical Review B 104, 075114 (2021).
  * Lupo and Ricci-Tersenghi [2017] C. Lupo and F. Ricci-Tersenghi, Approximating the xy model on a random graph with a q-state clock model, Physical Review B 95, 054433 (2017).
  * Ilker and Berker [2014] E. Ilker and A. N. Berker, Odd q-state clock spin-glass models in three dimensions, asymmetric phase diagrams, and multiple algebraically ordered phases, Physical Review E 90, 062112 (2014).
  * Chatterjee  _et al._ [2020] S. Chatterjee, S. Sutradhar, S. Puri, and R. Paul, Ordering kinetics in a q-state random-bond clock model: Role of vortices and interfaces, Physical Review E 101, 032128 (2020).
  * Ali  _et al._ [2024] A. Ali, H. Xu, W. Bernoudy, A. Nocera, A. D. King, and A. Banerjee, Quantum quench dynamics of geometrically frustrated ising models, arXiv preprint arXiv:2403.00091 (2024).
  * Negrete  _et al._ [2021] O. A. Negrete, P. Vargas, F. J. Peña, G. Saravia, and E. E. Vogel, Short-range berezinskii-kosterlitz-thouless phase characterization for the q-state clock model, Entropy 23, 1019 (2021).
  * Negrete  _et al._ [2018c] O. A. Negrete, P. Vargas, F. J. Peña, G. Saravia, and E. E. Vogel, Entropy and mutability for the q-state clock model in small systems, Entropy 20, 933 (2018c).
  * Aguilera  _et al._ [2022] M. A. Aguilera, F. J. Peña, O. A. Negrete, and P. Vargas, Otto engine for the q-state clock model, Entropy 24, 268 (2022).
  * Elitzur  _et al._ [1979] S. Elitzur, R. Pearson, and J. Shigemitsu, Phase structure of discrete abelian spin and gauge systems, Physical Review D 19, 3698 (1979).
  * Cardy [1980] J. Cardy, General discrete planar models in two dimensions: Duality properties and phase diagrams, Journal of Physics A: Mathematical and General 13, 1507 (1980).
  * Fröhlich and Spencer [1981] J. Fröhlich and T. Spencer, The kosterlitz-thouless transition in two-dimensional abelian spin systems and the coulomb gas, Communications in Mathematical Physics 81, 527 (1981).
  * Ortiz  _et al._ [2012] G. Ortiz, E. Cobanera, and Z. Nussinov, Dualities and the phase diagram of the p-clock model, Nuclear Physics B 854, 780 (2012).
  * Kumano  _et al._ [2013] Y. Kumano, K. Hukushima, Y. Tomita, and M. Oshikawa, Response to a twist in systems with z p symmetry: The two-dimensional p-state clock model, Physical Review B 88, 104427 (2013).
  * Kim [2017] D.-H. Kim, Partition function zeros of the p-state clock model in the complex temperature plane, Physical Review E 96, 052130 (2017).