# Negative/positive electrocaloric effect in antiferroelectric squaric acid

A.P. Moina  Institute for Condensed Matter Physics, 1 Svientsitskii St., 79011
Lviv, Ukraine
<EMAIL_ADDRESS>
###### Abstract

Using the previously developed model we study the electrocaloric effect (ECE)
in antiferroelectric crystals of squaric acid. At low temperatures, the
polarization reorientation by the electric field applied along the
crystallographic $a$ axis of these crystals is predicted to be a two-stage
process, with one of the sublattice polarizations switched twice, by 90∘ each
time. The $T$-$E$ landscapes of the model polarization and entropy are
explored. The ECE, characterized by the introduced electric Grüneisen
parameter, is found to be negative in the antiferroelectric phase. In the
intermediate ferrielectric phase with perpendicular sublattice polarizations
it is positive at low fields, but becomes slightly negative below the
transition to the ferroelectric phase. In the ferroelectric phase the ECE is
positive at all temperatures and fields. Maximal EC temperature shift
magnitude is predicted to be just less than $-3$ K at 200 kV/cm. The
supercritical behavior of the Grüneisen parameter in the crossover region
between two bicritical end points is explored.

##  I Introduction

The electrocaloric effect (ECE) is a change of temperature or entropy of a
sample by varying external electric field. The ECE observed in ferroics is
believed to be promising for potential applications in solid-state
refrigeration, as the experimentally obtained temperature shifts in these
materials can be quite significant (see the reviews [1, 2, 3]).

Normally, the external field applied to a dielectric orders the electric
dipoles in the system, thereby decreasing its entropy, or, if the process is
adiabatic, increasing its temperature. This is the conventional or positive
ECE ($dT/dE>0$). However, if the initial ordering of the dipoles is
antiferroelectric (AFE), then the electric field tends to destroy it and
increases the system entropy. In this case the ECE is negative or inverse, and
the sample temperature is decreased ($dT/dE<0$). This and other circumstances,
when the ECE in ferroics is negative, are substantially discussed in the
review [4]. Commonly, the negative ECE in antiferroelectrics coexists with the
positive one, observed in other field or temperature domains. Those are, for
instance, the temperatures above the Néel temperature and the fields above the
field-induced transition from the antiferroelectric to ferroelectric (FE)
phase, at which one of the sublattice polarizations is switched by
$180^{\circ}$, from the antiparallel position to the parallel one.

The magnitude of the EC temperature shift $\Delta T$ is expected to increase
in the vicinities of phase transitions, where the variation of entropy and
polarization should be most prominent. If a line of the phase transitions in
the $T$-$E$ plane exhibits a multicritical point, the largest $\Delta T$ often
occurs near that point [5, 6, 3]. In antiferroelectrics, the line of the AFE-
FE transition is of a special interest in this context. As follows from the
analysis of the simplest Kittel model, this line, terminating at the Néel
temperature at zero field, does have a tricritical point at the used in
calculations values of the Landau coefficients, and $\Delta T$ is indeed the
largest in its vicinity [5, 3].

Squaric acid, H2C4O4, is a classical two-dimensional antiferroelectric.
However, some unique peculiarities of the process of the AFE-FE switching in
it render the ECE in squaric acid particularly interesting to explore.

The high-temperature phase of squaric acid is tetragonal ($I4/m$), whereas the
low-temperature structure, formally monoclinic ($P2_{1}/m$), is
pseudotetragonal [7]. Parallel to the $ac$ plane sheets of C4O4 groups, linked
by hydrogen bonds along two perpendicular directions as shown in fig. 1, are
stacked along the unique $b$ axis. Below the transition temperature, 373 K, a
ferroelectric polarization arises in each sheet, and the neighboring ones
polarize alternatively in the opposite directions.

Figure 1: The pseudotetragonal fully ordered AFE structure of a sheet of
H-bonded C4O4 groups in the squaric acid, drawn using the neutron diffraction
data [8] (proton off-center displacements are exaggerated). Two adjacent
protons near each group are localized in the vicinal to this group sites on
the H-bonds, and the two others in the distal sites (a lateral configuration).
Note a correlation between the positions of the double bonds and of the
vicinal protons. One of four ground state configurations is shown; the three
others can be obtained by rotating the entire picture by a multiple of 90∘.
Direction of the sublattice polarization $\textbf{P}^{+}$ for the shown
configuration and the crystallographic axes $a$ and $c$ are indicated.

Protons in squaric acid move in double-well potentials between two off-center
cites on the H-bonds. The sublattice polarizations are formed [9] by the
protonic displacements and, most importantly, by switchable $\pi$-bond
dipoles, positions of which in a given C4O4 group are determined by the proton
arrangement around it (see fig. 1). Because of the pseudotetragonal symmetry
of the hydrogen bond network, the resulting ground state sublattice
polarizations can be oriented along one of two perpendicular axes, depending
on which two adjacent vicinal H-sites are occupied by the ordered protons
(fig. 2). These polarization axes, most likely, run close to the diagonals of
the $ac$ plane [9, 11].

Figure 2: The allowed directions of the sublattice polarization in fully
ordered squaric acid and the schemes of the associated lateral proton
configurations. Squares are the C4O4 groups, and black circles are the vicinal
H-sites, occupied by protons.

It has been suggested [9, 10] that, as two polarization axes are present, the
external electric field applied within the $ac$ plane can switch the
sublattice polarization not only by 180∘ at once, but also by two consecutive
(at different field magnitudes) steps of 90∘. In this case, between the
initial AFE and field-induced collinear FE phases, the intermediate
noncollinear phase with perpendicular sublattice polarizations (NC90) should
appear [12]. Switching of the fully ordered sublattice polarization by 90∘
requires relocation of two protons around each C4O4 group to the other
potential wells within the respective H-bonds, causing also switching of two
$\pi$-bonds to other C and O atoms of the group [9, 12]. Switching by 180∘
requires all four protons around each group to be relocated and three
$\pi$-bonds to be switched simultaneously, which is apparently more energy
consuming. Due to thermal fluctuations and spontaneous switching of dipole
moments of some H2C4O4 groups in presence of the bias field, the resulting
angle between the sublattice polarizations, as statistically averaged
quantities, may, in fact, be arbitrary at the elevated temperatures, and not
just 0/90/180∘.

The deformable two-sublattice proton ordering model [15, 12, 13, 14] has been
proposed in order to explore the electric field effects in squaric acid. The
model calculations predict [12, 13] that the polarization switching at low
temperatures is the two-step process for all orientations of the external
field within the $ac$ plane, except for some particular ones: i) when the
field is directed along the axis of spontaneous sublattice polarizations, the
intermediate noncollinear phase is absent, and the switching occurs as in
usual uniaxial antiferroelectrics by 180∘ at once, and ii) when the field is
at 45∘ to this axis, the field of the transition to the ferroelectric phase
tends to infinity.

A typical $T$-$E$ phase diagram of squaric acid, exhibiting the two-step
polarization reorientation, is presented in fig. 3 for the case of the field
$E_{1}$, directed along the crystallographic axis $a$. These diagrams have
been constructed for different field orientations and extensively discussed in
Refs. 12, 13, 14. There are three phases in the diagram: AFE*, non-collinear
antiferrielectric with nearly antiparallel and partially uncompensated
sublattice polarizations; FE, collinear field-induced ferroelectric 111The
AFE* phase is nothing but a tweaked by the external field initial fully
compensated collinear AFE phase. At low fields and high temperatures (above
$T_{{\rm N}0}$) the FE phase is the initial paraelectric phase, in which the
external field induces small polarization; the sublattice polarization vectors
are equal and parallel. With increasing field and lowering temperature,
circumventing lines IV and III, it gradually transforms into the fully ordered
ferroelectric phase. The term “polar” could be also appropriate for this
phase., and NC90, noncollinear ferrielectric, where the angle between the
sublattice polarizations mostly remains close to 90∘. The phases are separated
by the lines of the first order phase transitions I, II, and III, and of the
second order phase transitions IV, each terminating at some critical end point
(bicritical end points BCE, tricritical point TCP, or critical end point CEP).
At low temperatures, when the field is increased, one of the sublattice
polarizations is switched twice by approximately 90∘ each time, first at line
II, to become perpendicular to the other sublattice polarization, and then at
line III to become parallel to it. In the region between the critical points
BCE1 and BCE2, a crossover from the AFE* to NC90 phase occurs, at which the
sublattice polarization rotates continuously with increasing field.

The theoretical values of the AFE*-NC90 switching fields and the
$P_{1}(E_{1})$ polarization curves and the experimental data [11] are in a
reasonable agreement, which improves significantly when the experimental
samples are of a better quality and of the increased dielectric strength [14].
The NC90-FE transition has not been detected experimentally yet, as the
maximal magnitudes of the electric field that can be safely applied to the
best quality crystals of squaric acid [11] are still below 230 kV/cm. Also,
the calculated fields of the NC90-FE transition strongly depend on the choice
of the model parameter values [12, 13, 14]. Therefore, the position of line
III in the phase diagram (fig. 3) is nominal and subject to change, should new
experimental data become available, allowing us to ascertain the values of the
fitting parameters.

The motivation to explore the ECE in squaric acid is thus two-fold. First,
because of the existence of the unique intermediate noncollinear ferrielectric
phase NC90. It is interesting to determine the sign of the ECE in this phase
and to study its field and temperature variations, as they are likely to be
different from those in the AFE* and FE phases. Second, because of the
presence of two bicritical end points in the field range below the limit of
the experimental dielectric strength of the crystal. The ECE is expected to be
strongly enhanced in the vicinity of the critical end points.

In the next Section we present some basic formulae, used for description of
the ECE. Then we proceed to explore the $T$-$E$ dependences of the entropy and
polarization of squaric acid, and then to the ECE itself. The expressions for
the entropy and polarization, obtained using the proton ordering model [12,
15], are given in the Appendix.

##  II The calculations

By definition, the adiabatic EC temperature shift is the change of the system
temperature from $T^{init}$ to $T^{end}=T^{init}+\Delta T$ at changing the
external applied field from $\textbf{E}^{init}$ to $\textbf{E}^{end}$, while
the (molar) entropy of the system is kept constant

$S(T^{init}+\Delta T,\textbf{E}^{end})=S(T^{init},\textbf{E}^{init}).$ (1)

If $S={}$const, and external pressure is assumed to be zero, then from

$dS=\left(\frac{\partial
S}{\partial\textbf{E}}\right)_{T}d\textbf{E}+\left(\frac{\partial S}{\partial
T}\right)_{\textbf{E}}dT=0,$

it follows that

$dT=T\bm{\Gamma}_{E}d\textbf{E}.$ (2)

We introduce here the electric Grüneisen parameter

$\bm{\Gamma}_{E}=-\frac{1}{T}\left(\frac{\partial
S}{\partial\textbf{E}}\right)_{T}\left(\frac{\partial S}{\partial
T}\right)_{E}^{-1}=-\frac{V_{m}}{C_{E}}\left(\frac{\partial\textbf{P}}{\partial
T}\right)_{E}$ (3)

as a differential characteristic of the ECE, as opposed to $\Delta T$ being
the integral one. Here $V_{m}$ is the molar volume;
$C_{E}=T\left(\frac{\partial S}{\partial T}\right)_{E}$ is the molar heat
capacity at constant electric field, and the Maxwell relation

$\left(\frac{\partial
S}{\partial\textbf{E}}\right)_{T}=V_{m}\left(\frac{\partial\textbf{P}}{\partial
T}\right)_{E}$

was used. Efficiency of the EC heating/cooling is sometimes characterized by
the ECE strength or EC responsivity

$\frac{\Delta T}{\Delta
E}=\frac{T^{end}-T^{init}}{|\textbf{E}^{end}-\textbf{E}^{init}|}.$ (4)

We shall consider the case of the electric field $E_{1}$ applied along the
$a$-axis. We start from analyzing the $(T,E)$ landscapes of the molar entropy
$S$ and polarization $P_{1}$ of squaric acid, calculated using the previously
developed model [12]. A brief description of the model and the obtained
expressions for the entropy, Eq. (9), and polarization, Eq. (10), are given in
Appendix. All the relevant details, discussions, formulae, fitting procedure,
and values of the model parameters can be found in Refs. 15, 12, 13, 14. The
entropy is differentiated numerically with respect to temperature and field to
find the component of the Grüneisen parameter vector $\Gamma_{1E}$, hereafter
referred to as simply the Grüneisen parameter. To determine the EC temperature
shift we solve Eq. (1) with respect to $\Delta T$ at different $T_{init}$ and
$E_{1}^{end}$ and fixed $E_{1}^{init}=0$.

The contour plot of the molar entropy is presented in fig. 3. Bending of the
isentropic contour lines towards lower temperatures indicate the entropy
accumulation at the transition lines I and II and at the crossover between the
AFE* and NC90 phases. The entropy increases with the field in the AFE* phase
and decreases in FE and most of NC90 phase. According to Eq. (3) this means
negative and positive Grüneisen parameter $\Gamma_{1E}$, respectively, or, in
other words, lowering temperature down by the increasing field $E_{1}$ in the
AFE* phase and raising it up in the FE and NC90 phases. A better insight can
be obtained from the analysis of the field dependence of the order parameters,
illustrated in fig. 4.

Figure 3:  The $T$-$E_{1}$ phase diagram of squaric acid taken from Ref. 12,
overlapping the contour color plot of the molar entropy $S$. Solid and dashed
black lines indicate the first and second order phase transitions,
respectively. The open square $\square$, star ☆, and full circles $\bullet$
indicate the critical end point (CEP), tricritical point (TCP), and bicritical
end points (BCE), respectively. The white lines are the contour lines of
constant entropy.

As described in Appendix, the system behaviour is governed by four order
parameters $\eta_{1\pm}$ and $\eta_{2\pm}$, two for each sublattice, ‘+‘ and
‘$-$‘. The system is most ordered and has the lowest entropy, when the
absolute values of the order parameters are close to 1. This can be confirmed
by plotting the entropy, Eq. (9), as a function of, say, $\eta_{1-}$ and
$\eta_{2-}$ at constant temperature and electric field. As seen in the bottom
panel of fig. 4, at approaching the first order transition at lines I and II
from below, the absolute values of all order parameters decrease, which
results in a marked increase of the entropy (negative $\Gamma_{1E}$). At the
transition $\eta_{2-}$ jumps to a positive value, indicating the rotation of
the ‘$-$‘ sublattice polarization. Above the transition, the perfect
ferrielectric ordering is approached, as all four $|\eta_{f\pm}|$ increase.
The entropy decreases, and $\Gamma_{1E}$ becomes positive. The entropy
behavior in the vicinity of the NC90-FE phase transition at line III is
qualitatively similar, which means that $\Gamma_{1E}$ again becomes negative
at approaching this line from below and positive above it. However, at fields
this high, the system is almost completely ordered in both phases, especially
at low temperatures. The entropy, therefore, remains low and barely changes,
and the magnitudes of $\Gamma_{1E}$ in this region are expected to be low too.

Figure 4:  The field dependences of the molar entropy $S$ (top) and order
parameters $\eta_{f\pm}$ (bottom) at 300 and 371 K. The first order
transitions at lines I, II, and III from the $T$-$E_{1}$ phase diagram (fig.
3) are indicated by arrows. Vertical dotted lines indicate the jumps of the
order parameters $\eta_{1-}$ and $\eta_{2-}$ at the first order transitions.
The jump of entropy at line II is too small to be discernible.

The field variation of the order parameters and entropy at 371 K is stronger
than at 300 K. There is, however, a more important distinction. At 300 K this
variation is nearly symmetric/antisymmetric at approaching the transition
field at line II from below and from above. At 371 K, on the other hand, it is
highly asymmetric below and above line I. Note, for instance, that $|\partial
S/\partial E_{1}|$ and $|\partial\eta_{f-}/\partial E_{1}|$ nearly diverge at
approaching line II from both sides and line I from below (and do diverge at
the endpoints BCE2 and BCE1 of these lines), but remain moderate as line I is
approached from above. One should expect different behavior of the Grüneisen
parameter in the vicinities of lines I and II and a strong asymmetry of its
field dependence near line I.

The entropy and ECE behavior of squaric acid under external electric field
$E_{1}$ can then be summarized as follows. In the AFE* phase the field spoils
the perfect compensated AFE ordering and increases the entropy
($\Gamma_{1E}<0$). In the NC90 phase just above the switching, the field
increases the noncollinear ferrielectric ordering and decreases the entropy
($\Gamma_{1E}>0$). Much higher fields begin to destroy the perfect
noncollinear ordering in favor of the collinear FE one and increase the
entropy again ($\Gamma_{1E}<0$). In the FE phase the field only increases the
ordering and decreases the entropy ($\Gamma_{1E}>0$).

The temperature curves of the polarization $P_{1}$, presented in fig. 5, lead
to the same conclusions. The polarization increases with temperature in the
AFE* phase (negative $\Gamma_{1E}$ expected) and decreases in FE and most of
the NC90 phase (positive $\Gamma_{1E}$). At very high fields just below line
III, the polarization is almost independent of the field, but begins again to
slightly increase with temperature, implying $\Gamma_{1E}\lesssim 0$ here.

Figure 5:  The temperature dependences of polarization $P_{1}$ at different
values of the electric field $E_{1}$. Open circles show the phase transitions
at lines I-IV. The rounded peak at 150 kV/cm corresponds to the crossover
between the AFE* and NC90 phases.

The $T$-$E_{1}$ landscape of the Grüneisen parameter $\Gamma_{1E}$ in fig. 6
agrees with the qualitative predictions, following from the analysis of the
entropy and polarization. As one can see, the ECE is indeed negative
($\Gamma_{1E}<0$) within the AFE* phase. At lines I and II of the first order
phase transition to the FE and noncollinear NC90 phases, the Grüneisen
parameter becomes positive. In the NC90 phase $\Gamma_{1E}$ mostly stays
positive, but decreases with increasing field. Within the narrow stripe along
the line III it becomes negative again. In the FE phase $\Gamma_{1E}>0$ at any
temperature or field.

In the crossover region between the AFE* and NC90 phases the Grüneisen
parameter changes its sign at line V. As a continuation of the transition
lines I and II beyond the end points BCE1 and BCE2, this line represents yet
another way to mark the crossover between the two phases. Frequently, the
supercritical lines are drawn using the loci of the extrema of the response
functions – second derivatives of the thermodynamic potential (susceptibility,
compressibility, etc). In this case, the condition $\Gamma_{1E}=0$, according
to Eq. (3), means that $(\partial S/\partial E)_{T}=0$ and $(\partial
P/\partial T)_{E}=0$. In other words, the maxima of the observables, i.e. the
first derivatives of the thermodynamic potential: entropy isotherms or
polarization isofields, are used.

Figure 6:  The $T$-$E_{1}$ contour plot of the Grüneisen parameter. The solid
and dashed lines, as well as the symbols are the same as in fig. 3. Dash-
dotted lines indicate the change of sign of $\Gamma_{E}$.

The largest absolute values of $\Gamma_{1E}$ are predicted to occur in the
vicinities of the bicritical end points BCE1 and BCE2 in the AFE* phase, and
also next to BCE2 in the NC90 phase, where $\partial S/\partial
E_{1}\to\infty$, but not next to BCE1 in the NC90 phase, where $\partial
S/\partial E_{1}$ remains finite (see fig. 4). This confirms the predicted
earlier different behavior of the Grüneisen parameter near the two endpoints
BCE1 and BCE2 and its asymmetry below and above the switching field near BCE1.

Both findings are further illustrated by the shown in fig. 7 field dependences
of the Grüneisen parameter $\Gamma_{1E}$ in the crossover region. Near BCE2
(the left panel) all curves, corresponding to different temperatures, cross at
approximately the same “isosbestic” point, where $\Gamma_{1E}=0$. It reflects
the fact that line V near BCE2 is almost parallel to the temperature axis,
i.e. the crossover field is barely temperature dependent. Near BCE1 (the right
panel), the crossover field strongly increases with decreasing temperature.
The field dependences of the Grüneisen parameter near BCE2 are nearly
antisymmetric. The sharp negative and positive peaks below and above the
crossover point are smeared out with increasing temperature, as the system
moves away from BCE2. On the other hand, the field dependences of
$\Gamma_{1E}$ near BCE1 are highly asymmetric. Negative peaks are observed
below the crossover field, similarly smeared out, as the system moves away
from BCE1. On the other hand, above the crossover field $\Gamma_{1E}$
monotonically increases with $E_{1}$ into a rounded shoulder.

Figure 7: Dependences of the Grüneisen parameter on the field $E_{1}$ at
different temperatures in the vicinities of the BCE2 (left) and BCE1 (right)
points. Lines: the theory. Symbols: the power law fits according to Eq. (5)
with $G=3.6\cdot 10^{-5}$ cm/kV ($\square$), $2.6\cdot 10^{-5}$ cm/kV
($\circ$), and $3.1\cdot 10^{-5}$ cm/kV ($\bigtriangleup$).

The tails of all three peaks in fig. 7 at $T$ closest to $T_{\textrm{BCE}}$
are well approximated by the power law dependences

$\Gamma_{1E}=\text{sgn}(e-1)G|e-1|^{-\alpha},$ (5)

where $e$ is the applied field normalized per the coordinate of the respective
end point BCE1 or BCE2

$e=\frac{E_{1}}{E_{\textrm{BCE}}}.$

Those coordinates are quite sensitive to the choice of the model parameters
values. At the used in the present calculation ones [12, 13] they are:
$T_{\textrm{BCE2}}=307.4$ K, $E_{\textrm{BCE2}}=217.4$ kV/cm,
$T_{\textrm{BCE1}}=370.1$ K, $E_{\textrm{BCE1}}=111$ kV/cm. The exponent
$\alpha$ was found to be about 0.42 both below and above BCE2 and 0.48 below
BCE1. The coefficients $G$ are different for each branch and given in the
caption to fig. 5. The $G_{1E}(E_{1})$ curves at larger $|T-T_{\textrm{BCE}}|$
approach the power law dependences (5) asymptotically with increasing $|e-1|$.

It has been discussed earlier [12] that the noncollinear ferrielectric phase
within the narrow wedge between lines I and IV above CEP, may, in fact, be
different from the NC90 phase above line II. The differences between the field
behaviour of the order parameters, entropy, and the Grüneisen parameter above
line I and above line II seem to give further evidence for this assumption.
There is also a possibility that the system of the critical end points BCE1
and CEP is an artefact of the mean field approximation [17] used for the long-
range interactions, and in a more advanced approximation these end points
should merge into a single tricritical point. In that case, the wedge between
lines I and IV disappears, and the system behavior in this region will be
different. Additional experimental measurements and Monte-Carlo calculations
could shed some light on this problem.

Figure 8:  The $T^{\rm{init}}$-$E^{\rm{end}}$ contour plot of the adiabatic
ECE temperature shift by the field $E_{1}$ (right). $E_{1}^{\rm{init}}=0$. The
solid and dashed lines, as well as the symbols are the same as in fig. 3 and
correspond to $T^{init}$.

Finally, the $T^{init}$-$E_{1}^{end}$ contour plot of the EC temperature shift
$\Delta T$ is given in fig. 8. Since it is calculated for $E_{1}^{init}=0$,
the entire landscape bears the imprint of the system entropy behavior at zero
field, revealed as the jumps in $\Delta T$ at $T^{init}=T_{{\rm N}0}$ at all
fields. Note that other anomalies of $\Delta T$, visualized, for instance, as
bends of the contour lines, occur at $T^{end}$ and do not coincide with the
phase transition lines, drawn for $T^{init}$.

In the AFE* phase the EC temperature shift $\Delta T$ is negative, and its
magnitude increases with the end field $E_{1}^{end}$. At $T^{init}>T_{{\rm
N}0}$, it is positive at all values of $E_{1}^{end}$, both in the FE and in
NC90 phases. Below $T_{{\rm N}0}$ in the NC90 phase, the temperature shift
first remains negative despite the change of sign of $\Gamma_{E}$ at lines I,
II, and V. However, its amplitude decreases with increasing field
$E_{1}^{end}$, until the effect of positive $\Gamma_{E}$ compensates the
negative shift accumulated in the AFE* phase. At $E_{1}^{end}$ above 400 kV/cm
$\Delta T$ is positive at all temperatures.

The largest negative EC temperature shifts $\Delta T$, up to $-2.7$ K, are
predicted to be achieved in the vicinity of the crossover line V at around 200
kV/cm, i.e. with the experimentally attainable fields [11]. The respective
maximal ECE strength $\Delta T/\Delta E$ is about 0.013 K cm/kV. Positive EC
temperature shifts do not exceed 1 K at $E_{1}^{end}$ below 250 kV/cm. It is
predicted that $\Delta T$ can reach 8 K and more, however, the fields of the
order of 800 kV/cm are required for that, which is well beyond the obtained so
far largest dielectric strength of squaric acid crystals [11].

The field dependence of the adiabatic EC temperature shifts at different
$T^{init}$ is shown in fig. 9. At low fields $\Delta T$ can be well
approximated by quadratic functions $\Delta T=k(E_{1}^{end})^{2}$ for
$T^{init}$ both below and above $T_{{\rm N}0}$. The temperature dependent
coefficient $k$ is negative for $T^{init}<T_{{\rm N}0}$ and positive for
$T^{init}>T_{{\rm N}0}$. This is in agreement with the predicted scaling law
for the EC temperature shift in antiferroelectrics [16]. Just above the
transition (crossover) to the NC90 phase, $\Delta T$ is a concave down
function of $E_{1}^{end}$ at 300 K, concave up at 350 K, and a linear function
at 400 K. Small jumps of $\Delta T$ are observed at crossing the first order
transition lines II and III. The jump directions are opposite to the
directions of the respective polarization $P_{1}$ jumps. In the high-field
region of the FE phase (above line III of the phase diagram) $\Delta T$ is a
weakly nonlinear, concave down function of $E_{1}^{end}$.

Figure 9: The field dependence of the adiabatic EC temperature shifts $\Delta
T$ at $E_{1}^{init}=0$ and different $T^{init}$. Lines: the theory. Symbols:
quadratic fits $\Delta T=k(E_{1}^{end})^{2}$, where $k=-1.7\cdot 10^{-5}$ K
cm2/kV2 for $T^{init}=300$ K ($\bigcirc$), $-6.7\cdot 10^{-5}$ for 350 K
($\bigtriangleup$), and $1.7\cdot 10^{-5}$ for 400 K ($\square$). Anomalies of
$\Delta T$ in the vicinities of the phase transition lines II, III, IV and of
the crossover line V are indicated by the arrows.

##  III Concluding remarks

In the present paper we study the electrocaloric effect in squaric acid, two-
dimensional hydrogen-bonded pseudotetragonal antiferroelectric. The case of
the electric field $E_{1}$ applied along the crystallographic $a$ axis of the
crystal is considered. The $T$-$E_{1}$ phase diagram of this crystal is quite
rich [12], exhibiting three phases, including the noncollinear ferrielectric
one (NC90) with perpendicular sublattice polarizations, as well the first and
second order transitions, critical end points, etc. The field-induced
polarization switching is expected to be an unusual two-stage process, with
the sublattice polarization vector rotated by $90^{\circ}$ at each stage, and
the NC90 phase being the intermediate one between the AFE* and FE phases.

The calculations are performed using the previously elaborated [15, 12, 13,
14] deformable proton-ordering model, which has proved to be successful for
description of the squaric acid behavior at ambient conditions, under
hydrostatic pressure, and in external electric fields, applied within the
plane of the hydrogen bonds. The $T$-$E_{1}$ landscapes of the system entropy
and polarization are explored. The electric Grüneisen parameter $\Gamma_{1E}$
and ECE temperature shift $\Delta T$, being the local in temperature/field and
the integral characteristics of the ECE, respectively, are calculated.

It has been found that the ECE, as expected, is negative ($\Gamma_{1E}<0$) in
the AFE* phase, where the field destroys the perfect AFE order and increases
the entropy. The ECE is positive ($\Gamma_{1E}>0$) in the FE phase, where the
field enhances the order and decreases the entropy. In the intermediate
noncollinear phase NC90, the ECE is positive at fields just above the
switching from the AFE* phase, where the field increases the noncollinear
ferrielectric order. With increasing field $\Gamma_{1E}$ gradually decreases
and becomes negative again at fields below the switching to the FE phase,
where the field begins to destroy the perfect noncollinear ordering in favor
of the collinear FE one.

At low fields, the EC temperature shift is proportional to $(E_{1}^{end})^{2}$
at all temperatures, both in the AFE* and in the FE phases. The largest
calculated EC temperature shift below 230 kV/cm, which is the best obtained to
date dielectric strength of squaric acid samples [11], is expected to occur
along this line, reaching $-2.7$ K at 200 kV/cm with the EC responsivity of
$-0.013$ K cm/kV.

Most interesting are the peculiarities of $\Gamma_{1E}$ caused by the presence
of two bicritical end points BCE1 and BCE2 in the phase diagram. The line,
connecting BCE1 and BCE2, where $\Gamma_{1E}=0$ (line V in fig. 3), is
suggested to be used to mark the crossover between the AFE and NC90 phases.
Behavior of $\Gamma_{1E}$ is found to be qualitatively different in the
vicinities of BCE1 and BCE2: practically antisymmetric near BCE2 and
completely asymmetric near BCE1. It indicates that the noncollinear
ferrielectric phase near BCE1 might be different from the NC90 one, if only
the BCE1 endpoint is not the artefact of the used mean field approximation.
Further experimental studies and Monte-Carlo calculations are required to
clarify this situation.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Appendix A

The used deformable proton ordering model [15, 12] considers the motion of
protons in two alternating plane sublattices, to be denoted hereafter as ’+’
and ’$-$’, placed in the thermally expanding host lattice of heavy ions that
forms the double-well potentials for each proton. This motion is described by
Ising pseudospins, whose two eigenvalues $\sigma=\pm 1$ are assigned to two
equilibrium positions of each proton. The arrangement of pseudospins in each
sublattice corresponds to the structure presented in fig. 1.

It is assumed that electric fields $E_{1}$ and $E_{3}$ are applied to the
crystal along its $a$ and $c$ axes. These fields break the equivalence of the
hydrogen bonds linking the C4O4 groups along these axes, as well as the
equivalence of the two sublattices. That leaves us with four independent order
parameters [12] $\eta=\langle\sigma\rangle$, two for each sublattice:
$\eta_{1\pm}$ and $\eta_{2\pm}$. In absence of the electric fields we have
$\eta_{1+}=\eta_{2+}=-\eta_{1-}=-\eta_{2-}$. Polarization rotation by 90∘ in a
fully ordered system is signalled by the change of the sign of one of the
order parameters in the ‘$-$‘ sublattice.

The system Hamiltonian [12, 13] includes intra- and interplane long-range
interactions between the pseudospins, ensuring ferroelectric ordering within
each separate plane and alternation of polarizations in the stacked planes, as
well as the short-range interactions, ensuring localization of protons in the
lateral configurations (see fig. 1). The mean field and four-particle cluster
approximations are used for the long-range and short-range interactions,
respectively.

In the short-range Hamiltonian the usual Slater-Takagi type scheme [18, 19,
15, 12] of 16 degenerate levels of lateral/diagonal/single-ionized/double-
ionized proton configurations is exploited. The degeneracy of the energy
levels of the lateral and single-ionized configurations, possessing dipole
moments in the $ac$ plane, is removed by the electric fields $E_{1}$ and
$E_{3}$. Assignment of the dipole moments to these configurations is a crucial
point of the model, and it has been described earlier [12, 13] in great
detail. For instance, the vector of the dipole moment, assigned to the H2C4O4
group with protons in the lateral configuration shown in fig. 1 is [12]
$(2\mu^{H}+2\mu_{\parallel}^{\pi},0,-2\mu_{\perp}^{\pi})$. This is based on
the results of the Berry phase calculations [9], which have shown that the
ground-state sublattice polarization in this crystal is formed directly by
displacements of protons along the hydrogen bonds and, mostly, by the
electronic contributions of switchable $\pi$-bond dipoles, and have indicated
the directions for those contributions.

From the presented in Refs. 12, 13 thermodynamic potential of the model
system, its molar entropy is obtained straightforwardly

$\displaystyle\frac{\Delta
S}{R}=\frac{\mathrm{v}}{2R}\sum_{ij=1}^{3}c_{ij}^{(0)}\alpha_{i}^{(0)}u_{j}+\frac{1}{2}\left(\ln
D_{+}+\ln
D_{-}\right)+\frac{1}{4}\left[\ln(1-\eta^{2}_{1+})+\ln(1-\eta^{2}_{2+})+\ln(1-\eta^{2}_{1-})+\ln(1-\eta^{2}_{2-})\right]$
$\displaystyle{}-(\beta\nu-\frac{2\delta_{T}\lambda_{\delta}\nu_{0}}{k_{\rm
B}})\left[\frac{(\eta_{1+}-\eta_{1-})^{2}}{4}+\frac{(\eta_{2+}-\eta_{2-})^{2}}{4}\right]-(\beta\nu^{\prime}-\frac{2\delta_{T}\lambda_{\delta}\nu^{\prime}_{0}}{k_{\rm
B}})\left[\frac{(\eta_{1+}+\eta_{1-})^{2}}{4}+\frac{(\eta_{2+}+\eta_{2-})^{2}}{4}\right]$
$\displaystyle+\left(\beta-\frac{2\delta_{T}}{k_{\rm
B}\lambda_{\delta}}\right)\left[\frac{a\varepsilon+2bw(\cosh{z_{1+}}+\cosh{z_{2+}})}{2D_{+}}+\frac{a\varepsilon+2bw(\cosh{z_{1-}}+\cosh{z_{2-}})}{2D_{-}}\right]$
(6)
$\displaystyle-\beta\sum_{f=1}^{2}\frac{\eta_{f+}+\eta_{f-}}{4}(\mu_{f1}E_{1}+\mu_{f3}E_{3}),$

where $R$ is the gas constant, $\beta=(k_{\textrm{B}}T)^{-1}$, and

$\displaystyle D_{\pm}=a+\cosh(z_{1\pm}+z_{2\pm})+2b(\cosh z_{1\pm}+\cosh
z_{2\pm})+\cosh(z_{1\pm}-z_{2\pm}),$ $\displaystyle
z_{f\pm}=\frac{1}{2}\ln\frac{1+\eta_{f\pm}}{1-\eta_{f\pm}}\pm\beta\nu\frac{\eta_{f+}-\eta_{f-}}{2}+\beta\nu^{\prime}\frac{\eta_{f+}+\eta_{f-}}{2}+\beta\frac{\mu_{f1}E_{1}}{2}+\beta\frac{\mu_{f3}E_{3}}{2},\quad(f=1,2);$
(7) $\displaystyle a=\exp(-\beta\varepsilon),\quad b=\exp(-\beta w),$
$\displaystyle\mu_{11}=\mu_{23}=\mu^{H}+\mu^{\pi}_{\parallel}-\mu^{\pi}_{\perp},\quad-\mu_{13}=\mu_{21}=\mu^{H}+\mu^{\pi}_{\parallel}+\mu^{\pi}_{\perp}.$

The first term in (A) comes from the “seed” energy of thermal expansion [15]
of the host lattice; $u_{1}$, $u_{2}$, and $u_{3}$ are the uniform diagonal
lattice strains; $c_{ij}^{(0)}$ and $\alpha_{i}^{(0)}$ are the “seed” elastic
constants and thermal expansion coefficients; $\mathrm{v}$ is the unit cell
volume just above the AFE transition temperature $T_{{\rm N}0}$ in absence of
the electric fields.

The interaction parameters of the deformable model are temperature dependent.
Thus, the short-range Slater-Takagi type energy parameters $\varepsilon$ and
$w$ are considered [15] to be quadratic functions of the H-site distance
$\delta$. In its turn, the distance $\delta$ is taken to vary according to its
experimentally observed [8] above the transition linear temperature dependence

$\delta=\delta_{0}[1+\delta_{T}(T-T_{\textrm{N0}})]\equiv\delta_{0}\lambda_{\delta},$
(8)

where $T_{\textrm{N0}}$ is the transition temperature at zero electric field.
It yields

$\varepsilon=\varepsilon^{\prime}_{0}\lambda_{\delta}^{2},\quad
w=w_{0}\lambda_{\delta}^{2}.$

The long-range interaction parameters $\nu$ and $\nu^{\prime}$ from Eq. (A)
are taken to depend on $\delta$ and on the lattice strains as [15, 12]

$\nu=\nu_{0}\lambda_{\delta}^{2}+\sum_{i=1}^{3}\psi_{i}u_{i},\quad\nu^{\prime}=\nu^{\prime}_{0}\lambda_{\delta}^{2}+\sum_{i=1}^{3}\psi_{i}^{\prime}u_{i}.$

Finally, the dipole moments $\mu_{f1}$ and $\mu_{f3}$ are assumed to be
temperature or strain independent.

Differentiating Eq. (A) numerically, we obtain the contribution to the molar
specific heat of the crystal from its subsystem, described by the model
Hamiltonian

$\Delta C_{E}=T\left(\frac{\partial\Delta S}{\partial T}\right)_{E}.$

To get the total molar specific heat of the crystal that can be compared to
experimental data, the background regular contribution from the host lattice
has to be added to $\Delta C_{E}$

$C_{E}=\Delta C_{E}+C_{add}.$

A sufficiently good agreement with the experiment, as illustrated in fig. 10,
is obtained if the added regular part is a linear function of temperature

$C_{add}=k_{0}+k_{1}T,$

with $k_{0}=45.62$ J/(mol K), $k_{1}=0.217$ J/(mol K2). One can also see, that
$\Delta C_{E}$ does not represent a purely anomalous part of the total
specific heat, but has a significant regular part as well, which, actually, is
larger than $C_{add}$.

Figure 10:  The temperature dependence of the specific heat of squaric acid at
zero electric field. Symbols: experimental points, taken from Refs. 20
($\bigcirc$) and 21 ($\triangle$). Dashed line: the added regular part of the
specific heat, $C_{add}$.

Integrating $C_{add}/T$ over temperature we find the added contribution to the
entropy

$S_{add}=k_{1}(T-T_{{\rm N}0})+k_{0}\ln(T/T_{{\rm N}0})+\textrm{const}.$

The total entropy, to be used to explore the ECE in squaric acid, is

$S=\Delta S+S_{add}.$ (9)

Note that $S_{add}$ and $S$ are determined within a constant, which in the
calculations is chosen to yield $S_{add}(T_{{\rm N}0})=0$.

The net polarization $\mathbf{P}$ of squaric acid has been obtained earlier
[12, 14]. It is a sum of the sublattice polarization vectors
$\mathbf{P}_{\pm}$

$\mathbf{P}=\mathbf{P}_{+}+\mathbf{P}_{-},\quad\mathbf{P}_{\pm}=(P_{1\pm},0,{P_{3\pm}}),$
(10)

where [12]

$\displaystyle
P_{1\pm}=\frac{\mu^{H}+\mu^{\pi}_{\parallel}}{2{\textrm{v}}}(\eta_{1\pm}+\eta_{2\pm})+\frac{\mu_{\perp}^{\pi}}{2{\textrm{v}}}(\eta_{2\pm}-\eta_{1\pm}),$
$\displaystyle
P_{3\pm}=\frac{\mu^{H}+\mu^{\pi}_{\parallel}}{2{\textrm{v}}}(\eta_{2\pm}-\eta_{1\pm})-\frac{\mu_{\perp}^{\pi}}{2{\textrm{v}}}(\eta_{2\pm}+\eta_{1\pm}).$

Naturally, in absence of the electric fields $\mathbf{P}=0$.

The order parameters $\eta_{f\pm}$ and lattice strains $u_{i}$ are found by
numerical minimization of the presented in Refs. 12, 13 thermodynamic
potential of the system, using the given therein values of the model
parameters. The fitting procedure has been discussed earlier [15, 12, 14] very
extensively.

## References

  * [1] Y. Liu, J.F. Scott, and B. Dkhil, Appl. Phys. Rev. 3, 031102 (2016) 10.1063/1.4958327.
  * [2] X. Moya, S. Kar-Narayan, and N. D. Mathur, Nature Materials 13, 439 (2014) 10.1038/nmat3951.
  * [3] Z. Kutnjak, B. Rožič, and R. Pirc,“Electrocaloric Effect: Theory, Measurements, and Applications”, in Wiley Encyclopedia of Electrical and Electronics Engineering, edited by J. G. Webster (Wiley, Hoboken, 2015), p. 1. 10.1002/047134608X.W8244.
  * [4] A. Grünebohm, Y.-B. Ma, M. Marathe, B.-X. Xu et al, Energy Technol. 6, 1491 (2018) 10.1002/ente.v6.8.
  * [5] R. Pirc, B. Rožič, J. Koruza, B. Malič, and Z. Kutnjak, EPL, 107, 17002 (2014) 10.1209/0295-5075/107/17002
  * [6] L. Squillante, I.F. Mello, A.C. Seridonio, and M. de Souza, Materials Research Bulletin 142, 111413 (2021) 10.1016/j.materresbull.2021.111413
  * [7] D. Semmingsen, F. J. Hollander, and T. F. Koetzle, J. Chem. Phys. 66, 4405 (1977). 10.1063/1.433745
  * [8] D. Semmingsen, Z. Tun, R.J. Nelmes, R. K. McMullan, and T. F. Koetzle, Zeit. für Kristallographie 210, 934 (1995) 10.1524/zkri.1995.210.12.934.
  * [9] S. Horiuchi, R. Kumai, and S. Ishibashi, Chem. Sci. 9, 425 (2018) 10.1039/C7SC03859C
  * [10] S. Ishibashi, S. Horiuchi, and R. Kumai, Phys. Rev. B 97, 184102 (2018) 10.1103/PhysRevB.97.184102
  * [11] S. Horiuchi, S. Ishibashi, Chem. Phys. 12, 14198 (2021) 10.1039/d1sc02729h
  * [12] A. P. Moina, Phys. Rev. B 103, 214104 (2021) 10.1103/PhysRevB.103.214104
  * [13] A. P. Moina, Condens. Matter Phys. 24, 43703 (2021) 10.5488/CMP.24.43703
  * [14] A. P. Moina, Condens. Matter Phys. 25, to be published (2022).
  * [15] A. P. Moina, Condens. Matter Phys. 23, 33704 (2020) 10.5488/CMP.23.33704
  * [16] S. Lisenkov, B.K. Mani, E. Glazkova, C. W. Miller, and I. Ponomareva, Sci. Rep. 6, 19590 (2016) 10.1038/srep19590.
  * [17] J. A. Plascak, D. P. Landau, Phys. Rev. E 67, 015103(R) (2003) 10.1103/PhysRevE.67.015103.
  * [18] E. Matsushita and T. Matsubara, Progr. Theor. Phys. 64, 1176 (1980).
  * [19] E. Matsushita and T. Matsubara, Progr. Theor. Phys. 68, 1811 (1982).
  * [20] E. Barth, J. Helwig, H.-D. Maier, H.E. Müser, and J. Petersson, Z. Physik B 34, 393 (1979) 10.1007/BF01325204.
  * [21] F. Gronvold, Pure and Appl. Chem. 47, 251 (1976).


#  Random-search efficiency in a bounded interval with spatially heterogeneous
diffusion coefficient

L. Menon Jr.1, M. A. F. dos Santos1 and C. Anteneodo1,2 1 Department of
Physics, PUC-Rio, Rua Marquês de São Vicente 225, 22451-900, Rio de Janeiro,
RJ, Brazil  2 Institute of Science and Technology for Complex Systems, INCT-CS
Brazil

###### Abstract

We consider random walkers searching for a target in a bounded one-dimensional
heterogeneous environment, in the interval $[0,L]$, where diffusion is
described by a space-dependent diffusion coefficient $D(x)$. Boundary
conditions are absorbing at the position of the target (set at $x=0$) and
reflecting at the border $x=L$. We calculate and compare the estimates of
efficiency $\varepsilon_{1}=\langle 1/t\rangle$ and $\varepsilon_{2}=1/\langle
t\rangle$. For the Stratonovich framework of the multiplicative random
process, both measures are analytically calculated for arbitrary $D(x)$. For
other interpretations of the stochastic integrals (e.g., Itô and anti-Itô), we
get general results for $\varepsilon_{2}$, while $\varepsilon_{1}$ is obtained
for particular forms of $D(x)$. The impact of the diffusivity profile on these
measures of efficiency is discussed. Symmetries and peculiar properties arise
when the search starts at the border ($x_{0}=L$), in particular, heterogeneity
spoils the efficiency of the search within the Stratonovich framework, while
for other interpretations the searcher can perform better in certain
heterogeneous diffusivity profiles.

Keywords: random search, finite interval, heterogeneous diffusion, search
efficiency, first-arrival times

##  1 Introduction

The study of random searches is valuable in several areas of science and
technology [1], offering effective strategies for exploring the accessible
space and finding targets. Random searches can be found from the molecular
scale, as in problems of protein-DNA binding [2, 3, 4], to ecological scales,
where random searches are crucial in foraging for resources or viable
habitats, necessary for survival and reproduction [5, 6, 7, 8, 9, 10]. We can
also find technological applications, as in robotic design [11] or in the
optimization of hyperparameters in neural networks [12].

We address the situation in which a random search occurs in an heterogeneous
confined environment. We assume that the search dynamics follows a Brownian
(or diffusive) motion, such that a non-uniform, i.e., position-dependent,
diffusion coefficient reflects the heterogeneity of the environment. State-
dependent diffusion appears in the description of diverse physical systems
that exhibit spatial complexity, producing anomalous diffusion, e.g.,
subdiffusion [13] and turbulent diffusion [14], as well as ergodicity breaking
[15]. It also appears in the modeling of biological [16, 17], socioeconomic
[18], and ecological [19] systems, among other examples.

Random searches in heterogeneous media have been previously studied in non-
confined domains, exploring different diffusivity profiles within Stratonovich
scheme [20], for general schemes of the stochastic integrals [21], or with
stochastic resetting [22]. However, in real systems random searches occur in
bounded domains. As examples, reaction-diffusion processes take place in
limited plates [23], animal search for food occurs in geographically delimited
habitats [5], the visual search for an object in an image is bounded by image
size [24]. In such bounded domains, the homogeneous problem has been
extensively studied. A review about random searches in homogeneous bounded
media can be found in Ref. [25], after which, extensions such as searches
under stochastic resetting [26], and partially-reactive targets [27] have been
addressed. Studies about search processes in confined heterogeneous
environments, involving position-dependent diffusion coefficient, are
relatively recent, and address particular diffusivity profiles [28, 29] or
particular interpretations of the stochastic integrals, e.g., within Hänggi-
Klimontovich (or anti-Itô interpretation [30]). In this work, we aim to obtain
analytical results for the search problem in bounded one-dimensional
heterogeneous environments with general forms of the state-dependent
diffusivity and for general interpretation of the stochastic integrals.

For this purpose, we consider the stochastic equation

$\displaystyle\dot{x}=\sqrt{2D(x^{*})}\,\eta(t),$ (1)

where $\eta(x)$ represents a delta-correlated Gaussian noise and the
(definite-positive) function $D(x)$ denotes the state-dependent diffusion
coefficient. Importantly, since the white noise is multiplicative, it is
necessary to choose the specific instant at which $x^{*}$ is calculated. We
consider $x^{*}=[(2-A)\,x(t+dt)+A\,x(t)]/2$ [29], where $A\in[0,2]$, focusing
on three special interpretations: Itô ($A=2$) [31], Stratonovich ($A=1$) [32],
anti-Itô ($A=0$) [33, 34]. Alternatively, Eq. (1) can be cast in the Itô form
($x^{*}=x(t)$)

$\displaystyle\dot{x}=(1-A/2)D^{\prime}(x)+\sqrt{2D(x)}\,\eta(t),$ (2)

where $D^{\prime}=dD/dx$ and the first term is the noise-induced drift [35].
The Fokker-Planck equation associated to Eq. (2) is given by [36]

$\displaystyle\frac{\partial\ }{\partial t}p(x,t|x_{0})=\frac{\partial\
}{\partial x}\left\\{D(x)^{1-\frac{A}{2}}\frac{\partial}{\partial
x}[D(x)^{\frac{A}{2}}p(x,t|x_{0})]\right\\},$ (3)

where $x_{0}\in\Omega\equiv(0,L]$ is the initial position and the boundary at
$L$ is modelled by a reflecting boundary condition. Additionally, without loss
of generality we assume that the target is located at $x=0$, where we set an
absorbing boundary condition, meaning that the searcher is removed when the
search ends, upon reaching the target.

It is interesting to note that similar forms of spatial heterogeneity can be
introduced through the potential of a field acting on the Brownian particle
[37, 38]. Moreover, as can be seen in Eq. (2), non-uniform diffusivity
produces a spurious drift for interpretations other than Itô.

To analytically grasp random searches, it is central to determine the first-
passage time distribution (FPTD)

$\displaystyle\wp(t)=-\frac{d\ }{dt}Q(x_{0},t)\,,$ (4)

where $Q(x_{0},t)=\int_{\Omega}p(x,t|x_{0})dx$, is the survival probability at
time $t$. The FPTD represents the probability density function of first-
arrival times when the particles hit the target for the first time, and after
which are removed from the system [39, 40]. The FPTD can be used to reveal
efficient strategies that minimize the characteristic time to encounter a
target, or optimize other search criteria, which can be crucial in diverse
contexts and scales [39, 41, 42]. As a single-value measure, the mean first-
passage time (MFPT), $\left\langle
t\right\rangle=\int_{0}^{\infty}t\,\wp(t)dt$, is an important quantity to be
calculated, and in some cases may be enough to characterize a random search,
as well as in other first-passage problems, when it is finite and the
corresponding standard deviation relatively small. In the context of search
processes, associated to the so called cruise motion [43] in which the
searcher keeps exploring all the points along the trajectory, the search
efficiency is defined as number of visited targets over average number of
steps, which in the case of a single target can be interpreted as the inverse
of the MFPT [44],

$\displaystyle\varepsilon_{1}=\left\langle
t\right\rangle^{-1}=\left(\int_{0}^{\infty}t\,\wp(t)dt\right)^{-1}.$ (5)

However, this quantity vanishes in cases where the MFPT is divergent.
Alternatively, a different measure of efficiency has been proposed [45],
namely,

$\displaystyle\varepsilon_{2}=\left\langle
t^{-1}\right\rangle=\int_{0}^{\infty}t^{-1}\,\wp(t)dt.$ (6)

This measure is adequate for systems where the MFPT diverges [45], typically
when $\wp$ has long tails, while these tails do not contribute to the first-
order negative moment, which dismisses trajectories that take very long times
to reach the target, and preferentially weights the contribution of short
arrival times. Therefore, it is useful to consider additionally the so-called
reliability, which is the complement of the survival probability,
$1-Q(x_{0},t\to\infty)$, at asymptotically long times, providing the fraction
of particles that reached the target. Eq. (6) has been used to characterize
the performance of Lévy searches, facing multiple targets [46], under external
bias [44], comb structures [47] and asymmetric Lévy flights [48]. Let us note
that other negative moments might also be considered.

Our focus is to obtain analytically and compare the measures of search
performance defined in Eqs. (5) and (6), for general types of bounded
heterogeneous environments, characterized by different forms of $D(x)$ and
interpretations of the heterogeneous diffusion process (HDP). In particular,
within Stratonovich framework, we will be able to find general results for
arbitrary $D(x)$.

The remaining of the paper is organized as follows. In Sec. 2, under the
Stratonovich interpretation, we manage to calculate both $1/\langle t\rangle$
and $\langle 1/t\rangle$, for arbitrary profiles $D(x)$. In Sec. 3, we
consider arbitrary interpretations ($A\in[0,2]$). In this case, for general
$D(x)$, we obtain and discuss a closed expression for the MFPT, hence for
$\varepsilon_{1}$, while for the efficiency given by the first negative moment
$\varepsilon_{2}$, analytical expressions are obtained for particular shapes
of $D(x)$. From the analytical results, the impact of the diffusivity profile
and interpretation of the stochastic integrals on search performance is
discussed. Final remarks are presented in Sec. 4.

##  2 Search efficiency under Stratonovich interpretation

The backward FP equation associated to Eq. (3), under Stratonovich
interpretation ($A=1$), is given by [40]

$\displaystyle\frac{\partial\ }{\partial
t}Q(x_{0},t)=\sqrt{D(x_{0})}\frac{\partial\ }{\partial
x_{0}}\left\\{\sqrt{D(x_{0})}\frac{\partial\ }{\partial
x_{0}}Q(x_{0},t)\right\\},$ (7)

where $Q(x_{0},t)=\int_{\Omega}p(x,t|x_{0})dx$ is the survival probability.
The position of the target, located at $x=0$, can be represented by an
absorbing boundary condition, while the confinement of the accessible domain
is given by a reflecting boundary at $x=L$, namely,

$\displaystyle Q(x_{0}=0,t)$ $\displaystyle=$ $\displaystyle 0,$ (8)
$\displaystyle\frac{\partial\ }{\partial x_{0}}Q(x_{0}=L,t)$ $\displaystyle=$
$\displaystyle 0,$ (9)

for all $t$, where $0\leq x_{0}\leq L$.

To solve Eq. (7), we apply the Laplace transform in the temporal variable,
using
$\mathcal{L}\\{f(t)\\}\equiv\tilde{f}(s)=\int_{0}^{\infty}f(t)e^{-ts}ds$, and
additionally use the change of variables

$\displaystyle y(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{x}D(x^{\prime})^{-\frac{1}{2}}dx^{\prime},$ (10)

which, to be well-defined, requires that the integrand not grow faster than
$x^{2}$ at the origin. In such case, $y$ is a monotonically increasing
function of $x$. The integral $y(x)$ will play a central role, as it contains
the information about the profile $D(x)$, and, the fact that its integrand
depends only on $D(x)$ will print important properties to the heterogeneous
search process.

After this change of variables, the solution in Laplace space is (for details,
see A):

$\displaystyle\tilde{Q}(y_{0},s)$ $\displaystyle=$
$\displaystyle\frac{1}{s}\left(1-\frac{e^{2y_{L}\sqrt{s}-y_{0}\sqrt{s}}}{e^{2y_{L}\sqrt{s}}+1}-\frac{e^{y_{0}\sqrt{s}}}{e^{2y_{L\sqrt{s}}}+1}\right)$
(11) $\displaystyle=$
$\displaystyle\frac{1}{s}\left(1-\frac{\cosh([y_{L}-y_{0}]\sqrt{s})}{\cosh(y_{L}\sqrt{s})}\right),$

where we defined $y_{L}\equiv y(L)$ and $y_{0}\equiv y(x_{0})$. The Laplace
transformed FPTD associated with the survival probability in Eq. (11) is given
by

$\displaystyle\tilde{\wp}(y_{0},s)=1-s\tilde{Q}(y_{0},s)=\frac{\cosh([y_{L}-y_{0}]\sqrt{s})}{\cosh(y_{L}\sqrt{s})}.$
(12)

Using Eqs. (11) and (12), we will be able to determine the measures $\langle
1/t\rangle$ and $1/\langle t\rangle$ for an arbitrary $D(x)$, implicitly
embodied in $y$.

###  2.1  Obtaining $1/\langle t\rangle$

We first calculate the MFPT using $\wp(t)=-\frac{d\ }{dt}Q(x_{0},t)$, namely,

$\displaystyle\langle t\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}t\,\frac{d\
}{dt}Q(x_{0},t)\,dt=\tilde{Q}(y_{0},s\to 0).$ (13)

Thereby, using the form of $Q(x_{0},t)$ in Eq. (11), we have

$\displaystyle\langle t\rangle$ $\displaystyle=$ $\displaystyle
y_{L}y_{0}-\frac{y_{0}^{2}}{2},$ (14)

hence,

$\displaystyle\langle t\rangle^{-1}=\frac{2}{y_{0}(2y_{L}-y_{0})},$ (15)

recalling that $y_{L}=y(L)$, where $y(x)$ was defined in Eq. (10).

On the other hand, if $D(x)$ does not grow faster than $x^{2}$ for large $x$,
then, $L\to\infty$ implies $y_{L}\to\infty$. In such case we obtain that
$\lim_{L\to\infty}\langle t\rangle^{-1}=0$, because when exploring the semi-
infinite region, $\langle t\rangle$ increases indefinitely [39].

###  2.2  Obtaining $\langle 1/t\rangle$

The average $\langle 1/t\rangle$ defined in Eq. (6) can be rewritten in terms
of the FPTD (12) in Laplace space, as follows

$\displaystyle\langle t^{-1}\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}t^{-1}\,\wp(y_{0},t)\,dt=\int_{0}^{\infty}\tilde{\wp}(y_{0},s)\,ds$
(16) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{\cosh([y_{L}-y_{0}]\sqrt{s})}{\cosh(y_{L}\sqrt{s})}ds.$

Defining the new variable $z=e^{2y_{L}\sqrt{s}}$, Eq. (16) becomes

$\displaystyle\langle t^{-1}\rangle$ $\displaystyle=$
$\displaystyle\int_{1}^{\infty}\frac{z^{-\frac{1}{2}\frac{y_{0}}{y_{L}}}+z^{\frac{1}{2}\frac{y_{0}}{y_{L}}-1}}{2y_{L}^{2}}\frac{\log
z}{z+1}dz$ (17) $\displaystyle=$
$\displaystyle\frac{1}{8y_{L}^{2}}\left[\psi^{(1)}\left(\frac{y_{0}/y_{L}}{4}\right)-\psi^{(1)}\left(\frac{y_{0}/y_{L}+2}{4}\right)\right.$
$\displaystyle-\left.\psi^{(1)}\left(1-\frac{y_{0}/y_{L}}{4}\right)+\psi^{(1)}\left(\frac{1}{2}-\frac{y_{0}/y_{L}}{4}\right)\right],\;\;\;$

with $y_{0}/y_{L}\in[0,1]$, and where $\psi^{(1)}(z)$ is the polygamma
function [49].

###  2.3 Properties for $x_{0}=L$

Note that both $1/\langle t\rangle$ and $\langle 1/t\rangle$ depend on the
diffusivity profile through $y_{0}$ and $y_{L}$ only, which integrate a
function of the diffusivity from the target to $x_{0}$ and up to $L$
respectively. This implies that, curiously, shuffling the values of the
profiles within each one of the intervals $(0,x_{0})$ and $(x_{0},L)$ will not
alter the results.

The particular case $x_{0}=L$ means that the searcher is initially positioned
on the reflecting wall, which implies $y_{L}=y_{0}$. In such case, the results
are not altered by shuffling the values of the diffusivity in the whole
accessible region. Similar property was observed for the unbounded case
($L\to\infty$) [21]. In our case, the effect is illustrated in Fig. 1, where
we consider the family of profiles $D(x)=1+d\cos(n\pi x/L)$, with
$n=1,2,\ldots,10$, and $d=\pm 0.5$, which all contain the same values with
equal probability (as can be seen by stretching and unfolding). In fact, for
any integer $n$, the theoretical result, for each efficiency (large symbols),
remains invariant, in good agreement with the simulational data (small
symbols). Furthermore, we verified that for the shuffled profile shown in the
inset (after decomposition in 200 fragments), Langevin simulations also yield
the expected same level.

Let us remark that, changing the average level $D_{0}$ of the diffusivity
proportionally changes the efficiency measures, then, for fair comparisons in
all numerical examples, we set unitary average level.

Figure 1: Within Stratonovich framework, we plot $1/\langle t\rangle$ (large
triangles) and $\langle 1/t\rangle$ (large quadrangles), using Eqs. (15) and
(17), respectively, for the diffusivity profiles of the form
$D(x)=1+d\cos(n\pi x/L)$, with $d=\pm 0.5$ (illustrated in the inset for
$n=4$) as a function of integer values of $n$. The small symbols represent the
corresponding averages over $10^{6}$ realizations obtained from the
integration of the stochastic simulation of Eq. (2), using the Euler method
with a time step $10^{-5}\leq\Delta t\leq 10^{-3}$. In the inset, we included
a shuffled profile (after fragmentation in 200 segments), for which the same
invariant levels are obtained. The horizontal lines (solid for $1/\langle
t\rangle$ and dashed for $\langle 1/t\rangle$) represent the values for the
homogeneous case $d=0$, which performs better than the heterogeneous case,
according to both measures.

Moreover, for $x_{0}=L$, hence $y_{0}=y_{L}$, Eq. (15) becomes

$\displaystyle\langle t\rangle^{-1}\Big{|}_{y_{0}=y_{L}}=\frac{2}{y_{L}^{2}}.$
(18)

While, using Eq. (16) and the the change of variables $z=e^{y_{L}\sqrt{s}}$,
we obtain

$\displaystyle\langle t^{-1}\rangle\Big{|}_{y_{0}=y_{L}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{1}{\cosh(y_{L}\sqrt{s})}ds=\frac{4}{y_{L}^{2}}\int_{1}^{\infty}\frac{\log
z}{z^{2}+1}dz=\frac{4G}{y_{L}^{2}}$ (19) $\displaystyle\simeq$ $\displaystyle
1.83\times\frac{2}{y_{L}^{2}},$

where $G\simeq 0.915$ is the Catalan’s constant. Therefore, a general
proportionality exists between both measures independently of $D(x)$, when the
searcher is initially at the wall ($x_{0}=L$).

The limit $L\to\infty$, implying $y_{L}\to\infty$, when applied in Eq. (12),
leads to $\tilde{\wp}(y_{0},s)=e^{-y_{0}\sqrt{s}}$, which substituted in the
first line of Eq. (16) immediately yields

$\displaystyle\langle t^{-1}\rangle\Big{|}_{y_{L}\to\infty}$ $\displaystyle=$
$\displaystyle\frac{2}{y_{0}^{2}}.$ (20)

This result was already found in Ref. [21], where semi-infinite environment
was investigated. For such unbounded domain, particles far from the target
will likely reach it in so long times that will have a negligible contribution
to $\langle t^{-1}\rangle$. Here we remark the coincidence of Eq. (20) with
Eq. (18) when $x_{0}=L$. Results are summarized in Table 1.

| $x_{0}=L$ | $L\to\infty$
---|---|---
$1/\langle t\rangle$ | $2/y_{L}^{2}$ | 0
$\langle 1/t\rangle$ | $\sim 1.83\times 2/y_{L}^{2}$ | $2/y_{0}^{2}$
Table 1: Comparison of particular cases within Stratonovich framework.
Recalling that $(0,L]$ is the accessible region, $x_{0}$ the initial position
of the searcher for a target at $x=0$, and $y(x)$ is given by Eq. (10).

Given the inverse scaling of the efficiency with $y_{L}$, when $x_{0}=L$,
another important property emerges within the Stratonovich framework. On the
one hand, it can be shown [19] that

$y_{L}=\int_{0}^{L}D(x)^{-\frac{1}{2}}dx\geq LD_{0}^{-1/2},$ (21)

where $D_{0}$ is the average level of the diffusivity in $[0,L]$. Since the
efficiencies scale inversely with $y_{L}$, then the inequality (21) means that
in an heterogeneous (nonH) profile the searcher performs less efficiently than
in an homogeneous (H) environment with average level. That is, we find as a
general result, valid for the Stratonovich framework when $x_{0}=L$, that the
efficiency $\varepsilon$, measured either by $1/\langle t\rangle$ or
$1/\langle t\rangle$, verifies

$\varepsilon_{nonH}\leq\varepsilon_{H}.$ (22)

This effect is illustrated in Fig. 1, where the corresponding homogeneous
cases are represented by horizontal lines, above the heterogeneous values.

###  2.4 Properties for arbitray initial position $x_{0}\in(0,L]$

As paradigmatic example we first consider the linear profile
$D(x)=1+h(2x/L-1)$, which encompasses diffusivity profiles that increase
($0<h<1$) or decrease ($-1<h<0$) from the position of the target, as well as
the homogeneous case ($h=0$). Note that, if $h>0$ (resp. $<0$), the
diffusivity increases (resp. decreases) when moving away from the target).
Moreover this family has fixed (unitary) average level, which turns
comparisons fair. For these linear profiles, we plot, in Fig. 2, the
theoretical results for $\langle 1/t\rangle$ and $1/\langle t\rangle$ vs.
$x_{0}$ (a) and vs. $L$ (b), in good agreement with the simulations. In panel
(a), we can observe that, besides decaying with the distance $x_{0}$, as
expected, both measures approach each other when $x_{0}\to L$, reaching the
ratio $2G$ when $x_{0}=L$, as predicted by Eqs. (18) and (19).

These effects can be also observed (not shown) for the power-law family
$D(x)=x^{\alpha}$, where $\alpha>(<)0$ also represents profiles that increase
(decrease) from the target, although in this case the average level varies
with $\alpha$.

Figure 2:  $1/\langle t\rangle$ (solid lines) and $\langle 1/t\rangle$ (dashed
lines), using Eqs. (15) and (17), for the linear diffusivity profile
$D(x)=1+h(2x/L-1)$, for different values of $h$, as a function of the initial
position $x_{0}$, fixing $L=10$ (a) and as a function of the boundary position
$L$, fixing $x_{0}=0.1$ (b). The symbols represent the stochastic simulation
of Eq. (2), averaged over $10^{6}$ realizations.

With regard to the dependence on $L$, $1/\langle t\rangle$ decays with $L$,
since the MFPT grows with $L$, indicating that individual arrival times also
grow, however, these times do not have a significant contribution to $\langle
1/t\rangle$, which for large $L$ tends to a constant. This dichotomy reflects
the fact that for large $L$ both efficiencies measure different properties, to
the extent that long search times become more likely. For both measures, the
case $h<0$ (implying larger diffusivity near the target) allows a better
performance

##  3 Search efficiency under arbitrary interpretation of the HDP

The backward FP equation for all interpretations of the HDP is given by

$\displaystyle\frac{\partial\ }{\partial t}Q(x_{0},t)$ $\displaystyle=$
$\displaystyle D(x_{0})^{\frac{A}{2}}\frac{\partial\ }{\partial
x_{0}}\left\\{D(x_{0})^{1-\frac{A}{2}}\frac{\partial\ }{\partial
x_{0}}Q(x_{0},t)\right\\}\,,$ (23)

recalling that $Q(x_{0},t)$ is the survival probability and $x_{0}$ the
initial position.

###  3.1 MFPT for arbitrary $D(x)$

According to Eq. (23) the mean first-passage time, i.e., $\langle t\rangle$,
is given by

$\displaystyle\langle t\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{L}D(x^{\prime\prime})^{-\frac{A}{2}}dx^{\prime\prime}\int_{0}^{x_{0}}D(x^{\prime})^{-1+\frac{A}{2}}dx^{\prime}$
(24) $\displaystyle-$
$\displaystyle\int_{0}^{x_{0}}D(x^{\prime\prime})^{-1+\frac{A}{2}}\int_{0}^{x^{\prime\prime}}D(x^{\prime})^{-\frac{A}{2}}dx^{\prime}dx^{\prime\prime},$

whose derivation can be found in B. From Eq. (24), it is possible to calculate
$1/\langle t\rangle$ for all the HDP interpretations. In the case $A=1$
(Stratonovich), for $\langle t\rangle$, we recover Eq. (15).

With this general expression for the MFPT, we can study the implications of
Itô and anti-Itô interpretations.

Considering $x_{0}=L$ in Eq. (24), for $A=0$ and $A=2$, we have

$\displaystyle\langle t\rangle_{A=0}$ $\displaystyle=$
$\displaystyle\int_{0}^{L}\frac{L-x}{D(x)}dx=\int^{L}_{0}\frac{x^{\prime}}{D(L-x^{\prime})}dx^{\prime}$
(25)

and

$\displaystyle\langle t\rangle_{A=2}$ $\displaystyle=$
$\displaystyle\int_{0}^{L}\frac{x}{D(x)}dx.$ (26)

First we remark that, in contrast to Stratonovich framework ($A=1$), the
integrands in Eqs. (25) and (26) do not depend only on $D$. This means that
the insensitivity to the ordering of the heterogeneity, observed for the MFPT
when $A=1$, and $x_{0}=L$, is broken for other interpretations of the HDP,
implying that the shape of the profile is relevant and not only the
distribution of values within the accessible domain.

Second, note that Itô and anti-Itô interpretations produce the same outcomes
for profiles that are symmetric around $L/2$, when $x_{0}=L$. That is, the
reflection of the heterogeneity profile $D(x)\leftrightarrow D(L-x)$ emulates
the anticipating Itô ($A=0$) and non-anticipating anti-Itô ($A=2$) character
of the stochastic integration. This symmetry effect is illustrated in Fig.
3(a), using the linear profile $D(x)=1+h(2x/L-1)$, which is monotonic in the
accessible region $[0,L]$, and setting $x_{0}=L$. Notice, however that this
symmetry does not hold if $x_{0}<L$, as can be seen in Fig. 3(b).
Qualitatively similar results are obtained for nonlinear profiles, e.g., for
$D(x)=1+d\cos(\pi x/L)$ (not shown).

Figure 3:  Efficiency $\langle t\rangle^{-1}$, obtained from Eq. (24) (solid
lines) and from numerical simulations (symbols), for the diffusivity profile
$D(x)=1+h(2x/L-1)$ (depicted in the inset of panel (b)), for different values
of $A$, setting $L=1$, and (a) $x_{0}=1$ and (b) $x_{0}=0.6$. Notice in panel
(a) that since $x_{0}=L$, exchanging $x\leftrightarrow L-x$ is equivalent to
the change of interpretation parameter $A$: $0\leftrightarrow 2$, while this
symmetry is broken if $x_{0}<L$, as illustrated in panel (b).

Figure 4:  Efficiency $\langle t^{-1}\rangle$ obtained from Eq. (LABEL:eq:ps)
(solid lines) and from numerical simulations (symbols), with the same
diffusivity and parameters used in Fig. 3.

It is important to note, that this exchange symmetry is already present in
Langevin Eq. (2). Indeed, if we use Itô’s Lemma for the linear change of
variables $x=L-\bar{x}$, together with the exchange $A=2-\bar{A}$, the
Langevin equation remains the same. But the observed exchange symmetry emerges
only when the initial position is at the reflecting boundary $x_{0}=L$. Then,
similar behaviors are expected for $\langle 1/t\rangle$, as seen for the
linear profile in Fig. 4, which will be developed in the next section.

Another important implication of the results for $x_{0}=L$ is that while for
the Stratonovich interpretation the searcher performs better in homogeneous
media, for $A\neq 1$ there is a heterogeneous profile that optimizes the
search, although it is not the same for both efficiencies.

Finally, we observe that, when the profile increases with the distance from
the target $(h>0)$, the efficiency increases with larger $A$, while the
contrary occurs for a decreasing profile $(h<0)$. This effect has been
previously observed and explained for the unbounded domain in Ref. [21].
Differently, when $A=1$, the shape of the profile is not relevant, but only
its distribution of values, hence the results for given $h$ coincide with
those for $-h$.

###  3.2  Search efficiency $\langle 1/t\rangle$ for a special family of
profiles $D(x)$

Although the next calculations can be straightforwardly extended to the family
of profiles of the power-law form $D(x)=(ax+b)^{\alpha}$, such that $ax+b>0$
for $x\in[0,L]$ and $\alpha<2$, we will develop explicitly the linear case,
setting $\alpha=1$, which provide illustrative examples with less complicated
expressions.

Starting from Eq. (23), in C), we obtain the FPTD in Laplace space

$\displaystyle\tilde{\wp}(s)$ $\displaystyle=$
$\displaystyle\frac{(ax_{0}+b)^{\frac{\gamma}{2}}}{b^{\frac{\gamma}{2}}}\;\frac{I_{\gamma-1}\left[\xi(s,L)\right]\,K_{\gamma}\left[\xi(s,x_{0})\right]+K_{\gamma-1}\left[\xi(s,L)\right]\,I_{\gamma}\left[\xi(s,x_{0})\right]}{I_{\gamma-1}\left[\xi(s,L)\right]K_{\gamma}\left[\xi(s,0)\right]+K_{\gamma-1}\left[\xi(s,L)\right]I_{\gamma}\left[\xi(s,0)\right]},$

where $\gamma=A/2$ and $\xi(s,x_{0})=2\sqrt{s(ax_{0}+b)}/a$.

Then, we can calculate $\langle 1/t\rangle$ by numerical integration of
equation $\left\langle 1/t\right\rangle=\int_{0}^{\infty}\tilde{\wp}(s)ds$.
Results, compared to Langevin simulations, are displayed in Fig. 4 for the
linear profile with unitary average, $D(x)=1+h(2x/L-1)$. Note that the
symmetries and main features discussed for $1/\langle t\rangle$, when
$x_{0}=L$, also hold for $1/\langle t\rangle$.

##  4 Final remarks

We have calculated and compared two standard measures of efficiency for the
random search problem in a one-dimensional bounded interval where the
diffusion coefficient is space-dependent.

The first step was to consider a general profile $D(x)$ under the Stratonovich
interpretation ($A=1$), for which we discussed the similarities and
differences between both efficiencies. A relevant common feature is that they
do not depend on the sequence of values of $D(x)$, since the dependence on the
diffusivity profile occurs only through the integral $y$ of a function of
$D(x)$. Particular features emerge when $x_{0}=L$ (search initialized on the
reflecting wall), namely, we have shown that (i) heterogeneous profiles are
less efficient that the homogeneous one with same average diffusivity and (ii)
there is a proportionality between both measures, which is valid regardless of
the profile $D(x)$.

For general interpretation of stochastic integration, characterized by
parameter $A$, we obtained results for $1/\langle t\rangle$, valid for any
$D(x)$. From which, an important consequence is a symmetry property that
emerges when $x_{0}=L$, namely, the change $D(x)\leftrightarrow D(x-L)$ yields
the same results when changing the interpretation parameter $A$:
$0\leftrightarrow 2$. This symmetry also emerges for the measure of efficiency
$\langle 1/t\rangle$, for which we managed to obtain an analytical expression
for particular choices of $D(x)$, as the power-law diffusivity, among which we
developed the linear case, which embraces environments with increasing and
decreasing mobility versus the distance from the target, as well as the
homogeneous case. Furthermore, in contrast to the case $A=1$, for which the
homogeneous profile with equal average diffusivity allows a more efficient
search than heterogeneous ones, when $A\neq 1$, heterogeneous environments can
enhance the efficiency to reach the target. Another general feature is that,
increasing $A$ favors the search when the diffusivity increases with the
distance from the target and hinders the search otherwise. Moreover, this is
not unique to a particular shape of the diffusivity profile, or initial
position, but is determined by its monotonic character.

As perspectives of continuation, it would be interesting to extend the present
study to higher dimensions, consider colored instead of white noise, introduce
stochastic resetting, among other variants.

Acknowledgments: We all acknowledge partial financial support by the
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior \- Brazil (CAPES)
- Finance Code 001. C.A. also acknowledges partial support by Conselho
Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and Fundação de
Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ).

## References

## References

  * [1] Chupeau M, Bénichou O and Voituriez R 2015 Nature Phys. 11 844–847
  * [2] Mirny L, Slutsky M, Wunderlich Z, Tafvizi A, Leith J and Kosmrlj A 2009 J. Phys. A: Math. Theor. 42 434013
  * [3] Chen X, Cheng X, Kang Y and Duan J 2019 J. Stat. Mech. 2019 033501
  * [4] Bhattacherjee A and Levy Y 2014 Nucleic acids research 42 12404–12414
  * [5] Viswanathan G M, Buldyrev S V, Havlin S, Da Luz M, Raposo E and Stanley H E 1999 Nature 401 911–914
  * [6] O’brien W J, Browman H I and Evans B I 1990 American Scientist 78 152–160
  * [7] Bartumeus F, da Luz M G E, Viswanathan G M and Catalan J 2005 Ecology 86 3078–3087
  * [8] Viswanathan G M, Da Luz M G, Raposo E P and Stanley H E 2011 The physics of foraging: an introduction to random searches and biological encounters (Cambridge University Press)
  * [9] Martínez-García R, Calabrese J M, Mueller T, Olson K A and López C 2013 Phys. Rev. Lett. 110 248106
  * [10] Bhat U and Redner S 2022 J. Stat. Mech. 2022 033402
  * [11] Castello E, Yamamoto T, Libera F D, Liu W, Winfield A F, Nakamura Y and Ishiguro H 2016 Swarm Intelligence 10 1–31
  * [12] Bergstra J and Bengio Y 2012 J. Machine learning research 13
  * [13] Zodage A, Allen R J, Evans M R and Majumdar S N 2023 J. Stat. Mech. 2023 033211
  * [14] Stella A L, Chechkin A and Teza G 2023 Phys. Rev. E 107 054118
  * [15] Cherstvy A G, Chechkin A V and Metzler R 2013 New J. Phys. 15 083039
  * [16] English B P, Hauryliuk V, Sanamrad A, Tankov S, Dekker N H and Elf J 2011 Proc. Nat. Acad. Sci. 108 E365–E373
  * [17] Martinello M, Hidalgo J, Maritan A, Di Santo S, Plenz D and Muñoz M A 2017 Phys. Rev. X 7 041071
  * [18] Vieira A R and Anteneodo C 2018 Phys. Rev. E 97 052106
  * [19] Dos Santos M, Dornelas V, Colombo E and Anteneodo C 2020 Phys. Rev. E 102 042139
  * [20] Mutothya N M, Xu Y, Li Y, Metzler R and Mutua N M 2021 J. Phys.: Complexity 2 045012
  * [21] Dos Santos M, Menon Jr L and Anteneodo C 2022 Phys. Rev. E 106 044113
  * [22] Ray S 2020 J. Chem. Phys. 153
  * [23] Rice S A 1985 Diffusion-limited reactions (Elsevier)
  * [24] Wolfe J M and Horowitz T S 2004 Nature Rev. Neuro. 5 495–501
  * [25] Bénichou O and Voituriez R 2014 Phys. Reports 539 225–284 ISSN 0370-1573 from first-passage times of random walks in confinement to geometry-controlled kinetics
  * [26] Christou C and Schadschneider A 2015 J. Phys. A: Math. Theor. 48 285003
  * [27] Bressloff P C 2022 J. Phys. A: Math. Theor. 55 205001
  * [28] Li Y, Mei R, Xu Y, Kurths J, Duan J and Metzler R 2020 New J. Phys. 22 053016
  * [29] Vaccario G, Antoine C and Talbot J 2015 Phys. Rev. Lett. 115 240601
  * [30] Godec A and Metzler R 2015 Phys. Rev. E 91 052134
  * [31] Itô K 1944 Proc. Impe. Acad. 20 519–524
  * [32] Stratonovich R 1966 SIAM Journal on Control 4 362–371
  * [33] Hanggi P 1982 Phys. Rev. A 25 1130
  * [34] Klimontovich Y L 1990 Phys. A: Stat. Mech. App. 163 515–532
  * [35] Volpe G and Wehr J 2016 Rep. Prog. Phys. 79 053901
  * [36] Leibovich N and Barkai E 2019 Phys. Rev. E  99 042138
  * [37] Palyulin V V and Metzler R 2012 J. Stat. Mec.  2012 L03001
  * [38] Chupeau M, Gladrow J, Chepelianskii A, Keyser U F and Trizac E 2020 Proc. Nat. Acad. Sci. 117 1383-1388
  * [39] Redner S 2001 A guide to first-passage processes (Cambridge university press)
  * [40] Risken H and Risken H 1996 Fokker-planck equation (Springer)
  * [41] Bénichou O, Loverdo C, Moreau M and Voituriez R 2011 Rev. Mod. Phys. 83 81
  * [42] Zaburdaev V, Denisov S and Klafter J 2015 Rev. Mod. Phys. 87(2) 483–530
  * [43] John O’Brien W, Evans B I and Browman H I 1989 Oecologia 80 100–110
  * [44] Palyulin V V, Chechkin A V and Metzler R 2014 J. Stat. Mech. 2014 P11031
  * [45] Palyulin V V, Chechkin A V and Metzler R 2014 Proc. Nat. Acad. Sci. 111 2931–2936
  * [46] Palyulin V V, Mantsevich V N, Klages R, Metzler R and Chechkin A V 2017 Eur. Phys. J. B 90 1–16
  * [47] Sandev T, Iomin A and Kocarev L 2019 J. Phys. A: Math. Theor. 52 465001
  * [48] Padash A, Sandev T, Kantz H, Metzler R and Chechkin A V 2022 Fractal and Fractional 6 ISSN 2504-3110
  * [49] Abramowitz M and Stegun I A 1968 Handbook of mathematical functions with formulas, graphs, and mathematical tables vol 55 (US Government printing office)
  * [50] Bowman F 2012 Introduction to Bessel functions (Courier Corporation)

##  Appendix A FPTD in Laplace space for Stratonovich case

Using the change of variables defined in Eq. (10), for the new variable $y$,
Eq. (7) becomes

$\frac{\partial}{\partial t}Q(y_{0},t)=\frac{\partial^{2}}{\partial
y_{0}^{2}}Q(y_{0},t).$ (28)

Upon performing a Laplace transform with respect to time, the equation becomes

$s\tilde{Q}(y_{0},s)-1=\frac{\partial^{2}}{\partial
y_{0}^{2}}\tilde{Q}(y_{0},s).$ (29)

This leads to

$\tilde{Q}(y_{0},s)=\frac{1}{s}\left(1+B_{1}(s)e^{-y_{0}\sqrt{s}}+B_{2}(s)e^{y_{0}\sqrt{s}}\right).$
(30)

Applying the boundary conditions $Q(0,s)=0$ and
$\left.\partial_{y}Q\right|{y_{L}}=0$, we obtain

$\displaystyle B_{1}(s)$ $\displaystyle=$ $\displaystyle-
B_{2}(s)-1=e^{2y{L}\sqrt{s}}B_{2}(s)\,.$

Solving for $B_{1}(s)$ and $B_{2}(s)$, gives

$\displaystyle B_{1}(s)$ $\displaystyle=$
$\displaystyle\frac{-e^{2y_{L}\sqrt{s}}}{e^{2y_{L}\sqrt{s}}+1},$
$\displaystyle B_{2}(s)$ $\displaystyle=$
$\displaystyle\frac{-1}{e^{2y_{L}\sqrt{s}}+1}.$

Substituting these results back into Eq. (30), we get Eq. (11).

##  Appendix B MFPT for arbitrary $A$

The Laplace transform of Eq. (23) is

$\displaystyle s\tilde{Q}_{s}(x_{0})-1$ $\displaystyle=$ $\displaystyle
D(x_{0})^{\frac{A}{2}}\frac{\partial\ }{\partial
x_{0}}\left\\{D(x_{0})^{1-\frac{A}{2}}\frac{\partial\ }{\partial
x_{0}}\tilde{Q}_{s}(x_{0})\right\\}.$ (31)

By defining the MFPT $T(x_{0})=\langle t\rangle=\lim_{s\to
0}\tilde{Q}_{s}(x_{0})$, Eq. (31) implies

$\displaystyle-D(x_{0})^{-\frac{A}{2}}$ $\displaystyle=$
$\displaystyle\frac{\partial\ }{\partial
x_{0}}\left\\{D(x_{0})^{1-\frac{A}{2}}\frac{\partial\ }{\partial
x_{0}}T(x_{0})\right\\}.$

After the first integration in the variable $x_{0}$, we obtain

$\displaystyle\frac{\partial\ }{\partial x_{0}}T(x_{0})$ $\displaystyle=$
$\displaystyle
BD(x_{0})^{-1+\frac{A}{2}}-D(x_{0})^{-1+\frac{A}{2}}\int_{0}^{x_{0}}D(x^{\prime})^{-\frac{A}{2}}dx^{\prime}\,,$

where $B$ is a constant. After the second integration in $x_{0}$, we get

$\displaystyle T(x_{0})$ $\displaystyle=$ $\displaystyle
B\int_{0}^{x_{0}}D(x^{\prime})^{-1+\frac{A}{2}}dx^{\prime}$ (32)
$\displaystyle-$
$\displaystyle\int_{0}^{x_{0}}D(x^{\prime\prime})^{-1+\frac{A}{2}}\int_{0}^{x^{\prime\prime}}D(x^{\prime})^{-\frac{A}{2}}dx^{\prime}dx^{\prime\prime}.$

To find out the constant $B$, we consider the reflection condition (that is,
$\partial_{x_{0}}\tilde{Q}|_{x_{0}=L}=\partial_{x_{0}}T|_{x_{0}=L}=0$, for
$s\to 0$), then

$\displaystyle
BD(L)^{-1+\frac{A}{2}}-D(L)^{-1+\frac{A}{2}}\int_{0}^{L}D(x^{\prime})^{-\frac{A}{2}}dx^{\prime}=0,$

that leads to

$\displaystyle B$ $\displaystyle=$
$\displaystyle\int_{0}^{L}D(x^{\prime})^{-\frac{A}{2}}dx^{\prime}.$ (33)

Substituting $B$ into Eq. (32), we arrive to Eq. (24).

##  Appendix C Survival probability for $D(x)=ax+b$

Applying the transformation defined in Eq. (10) into Eq. (23), we obtain

$\displaystyle\frac{\partial\ }{\partial t}Q(y_{0},t)$ $\displaystyle=$
$\displaystyle\frac{1-A}{2\sqrt{D(x)}}\left.\frac{dD(x)}{dx}\right|_{y_{0}}\frac{\partial\
}{\partial y_{0}}Q(y_{0},t)+\frac{\partial^{2}\ }{\partial
y_{0}^{2}}Q(y_{0},t)\,,$ (34)

which for $A=1$ recovers Eq. (29).

For a linear profile $D(x)=ax+b$, Eq. (10) becomes
$y(x)=2(\sqrt{D(x)}-\sqrt{b})/a$. Then Eq. (34) can be rewritten as

$\displaystyle\frac{\partial\ }{\partial t}Q(y_{0},t)$ $\displaystyle=$
$\displaystyle\frac{1-A}{\mu+y_{0}}\frac{\partial\ }{\partial
y_{0}}Q(y_{0},t)+\frac{\partial^{2}\ }{\partial y_{0}^{2}}Q(y_{0},t)\,,$ (35)

where $\mu=2\sqrt{b}/a$.

The Laplace transform of Eq. (35) is given by

$\displaystyle\frac{\partial^{2}\ }{\partial
y_{0}^{2}}\tilde{Q}(y_{0},s)+\frac{1-A}{\mu+y_{0}}\frac{\partial\ }{\partial
y_{0}}\tilde{Q}(y_{0},s)-s\tilde{Q}(y_{0},s)+1=0\,,$ (36)

and using the ansatz

$\tilde{Q}(y_{0},s)=(y_{0}+\mu)^{\frac{A}{2}}R(y_{0},s)+1/s\,,$ (37)

we rewrite Eq. (36) as

$\frac{\partial^{2}}{\partial
y_{0}^{2}}R(y_{0},s)+\frac{1}{(y_{0}+\mu)}\frac{\partial}{\partial
y_{0}}R(y_{0},s)-\left[\left(\frac{A}{2(y_{0}+\mu)}\right)+s\right]R(y_{0},s)=0,$

which can be identified as a generalization of the Bessel equation of order
$A/2$ [50], whose solution is the linear combination

$R(y_{0},s)=c_{1}K_{\frac{A}{2}}\left[\sqrt{s}(y_{0}+\mu)\right]+c_{2}I_{\frac{A}{2}}\left[\sqrt{s}(y_{0}+\mu)\right].$
(38)

From Eq. (37), we have

$\tilde{Q}(y_{0},s)=(y_{0}+\mu)^{\frac{A}{2}}\left(c_{1}K_{\frac{A}{2}}\left[\sqrt{s}(y_{0}+\mu)\right]+c_{2}I_{\frac{A}{2}}\left[\sqrt{s}(y_{0}+\mu)\right]\right)+\frac{1}{s}\,.$
(39)

The boundary condition $\tilde{Q}(0,s)=0$ gives

$c_{1}=-\frac{\mu^{-\frac{A}{2}}}{sK_{\frac{A}{2}}(\mu\sqrt{s})}-c_{2}\,\frac{I_{\frac{A}{2}}(\mu\sqrt{s})}{K_{\frac{A}{2}}(\mu\sqrt{s})},$
(40)

and the reflection boundary condition, $\left.\frac{\partial
Q(y_{0},s)}{\partial y_{0}}\right|_{y_{0}=y_{L}}=0$, gives

$\displaystyle c_{2}$ $\displaystyle=$ $\displaystyle
c_{1}\frac{K_{\frac{A}{2}-1}[(y_{L}+\mu)\sqrt{s}]}{I_{\frac{A}{2}-1}[(y_{L}+\mu)\sqrt{s}]}.$

Combining Eqs. (40) and (C), we obtain

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{\mu^{A/2}s}\;\frac{I_{\frac{A}{2}-1}\left[(y_{L}+\mu)\sqrt{s}\right]}{I_{\frac{A}{2}-1}\left[(y_{L}+\mu)\sqrt{s}\right]K_{\frac{A}{2}}\left[\mu\sqrt{s}\right]+K_{\frac{A}{2}-1}\left[(y_{L}+\mu)\sqrt{s}\right]I_{\frac{A}{2}}\left[\mu\sqrt{s}\right]}.$

Finally, inserting Eq. (39) into the relation
$\tilde{\wp}(s)=1-s\tilde{Q}(s)$, we obtain the FPTD in Laplace space

$\displaystyle\wp(s)$ $\displaystyle=$
$\displaystyle\frac{(y_{0}+\mu)^{\frac{A}{2}}}{\mu^{A/2}}\frac{I_{\frac{A}{2}-1}\left[y_{L}^{\prime}\sqrt{s}\right]K_{\frac{A}{2}}\left[y_{0}^{\prime}\sqrt{s}\right]+K_{\frac{A}{2}-1}\left[y_{L}^{\prime}\sqrt{s}\right]I_{\frac{A}{2}}\left[y_{0}^{\prime}\sqrt{s}\right]}{I_{\frac{A}{2}-1}\left[y_{L}^{\prime}\sqrt{s}\right]K_{\frac{A}{2}}\left[\mu\sqrt{s}\right]+K_{\frac{A}{2}-1}\left[y_{L}^{\prime}\sqrt{s}\right]I_{\frac{A}{2}}\left[\mu\sqrt{s}\right]},$
(41)

where we defined $y_{0}^{\prime}\equiv y_{0}+\mu$ and $y_{L}^{\prime}\equiv
y_{L}+\mu$. Recalling the change of variables used in Eq. (35),
$y(x)=2\sqrt{ax+b}/a-\mu$, Eq. (41) immediately leads to Eq. (LABEL:eq:ps).