# Rhodium based half-Heusler alloys as possible optoelectronic and
thermoelectric materials

Dhurba R. Jaishi  Sujit Bati  Nileema Sharma  Bishnu Karki  Bishnu P. Belbase
Madhav Prasad Ghimire<EMAIL_ADDRESS>Central Department of
Physics, Tribhuvan University, Kirtipur 44613, Kathmandu, Nepal & Condensed
Matter Physics Research Center (CMPRC), Butwal 32907, Rupandehi, Nepal

###### Abstract

On the basis of density functional theory and semi-classical Boltzmann theory,
we have investigated the structural, elastic, electronic, optical and
thermoelectric properties of 18–valence electron count rhodium based half-
Heusler alloys focusing on RhTiP, RhTiAs, RhTiSb, and RhTiBi. The absence of
imaginary frequencies in the phonon dispersion curve for these system verifies
that they are structurally stable. RhTiP is ductile in nature, while others
are brittle. The alloys are found to be semiconducting with indirect band gaps
ranging from 0.94 to 1.01 eV. Our calculations suggest these materials to have
high absorption coefficient and optical conductivity in the ultraviolet as
well as visible region. While considering thermoelectricity, we found that
$p$–type doping is more favorable in improving the thermoelectric properties.
The calculated values of power factor with $p$-type doping are comparable to
some of the reported half-Heusler materials. The optimum figure of merit $ZT$
is $\sim 1$ for RhTiBi suggesting it as a promising candidate for
thermoelectric applications while RhTiP, RhTiAs, and RhTiSb with optimum $ZT$
values between 0.38 to 0.67 are possible candidates for use in thermoelectric
devices.

##  I Introduction

In recent years, research on materials science are found to focus in search
for novel materials suitable for applications in optoelectronics, spintronics,
thermoelectrics, etc. For these, various compounds including oxides,
perovskites, Heusler alloys, chalcogenides, and many other groups were
explored to tune their properties.Hu (2012); Leskelä  _et al._ (2019); Ghimire
_et al._ (2010, 2015); Kumar  _et al._ (2020); Ghimire and Richter (2017);
Roknuzzaman  _et al._ (2017) The growing demand for energy and over reliance
on fossil fuels have alerted the scientific community to seek alternative
sources of energy. Another crucial aspect to look into for energy conservation
is to reduce power loss. A significant portion of total energy is wasted in
form of unusable heat every year.Fitriani  _et al._ (2016) Thus, technology to
recuperate some of the wasted energy can be very helpful in solving the
problem of energy shortage and escalating energy costs. This can also help in
reducing the emission of greenhouse gases caused by the excessive use of
fossil fuels. One of the techniques used to generate electricity from waste-
heat generated from vehicles, factories, etc. is the thermoelectrics (TE).
Rowe (2018); Snyder and Toberer (2008); Sootsman  _et al._ (2009); Tritt
(2000)

The suitability of materials for TE application is measured in terms of a
dimensionless figure of merit, $ZT=\sigma\alpha^{2}T/\kappa$, where $\sigma$
is the electrical conductivity, $\alpha$ is the Seebeck coefficient, $T$ is
the temperature, and $\kappa$ is the thermal conductivity. The thermal
conductivity consists of the contribution from the lattice thermal
conductivity ($\kappa_{l}$) and the electronic thermal conductivity
($\kappa_{e}$). Enhancement of TE efficiency of the system can be done by
tailoring these interrelated parameters. The value of $ZT$ can be improved
through reduction of the lattice thermal conductivity and increasing the power
factor, and this can be achieved through various means like doping, band
engineering, and reducing dimension of the materials.Poudel  _et al._ (2008);
Mao  _et al._ (2017); Zhao  _et al._ (2016); Pei  _et al._ (2012); Neophytou
and Kosina (2011)

Among others, half-Heusler (hH) compounds are found to be mechanically stable
and possess magnificent thermal stability and large power factor.He  _et al._
(2016); Rogl  _et al._ (2016); Rausch  _et al._ (2015) These features make hH
alloys suitable for high temperature TE applications. Remarkably, many doped
hH alloys are found to have a good TE properties with $ZT$ values that are
comparable to the bismuth, lead, silicon-germanium containing conventional TE
materials such as Bi2Te3, BiSbTe, PbTe, BiSb, and SiGe.Goldsmid and Douglas
(1954); Poudel  _et al._ (2008); Gelbstein  _et al._ (2005); Poudeu  _et al._
(2006); Lee  _et al._ (2012) However, high values of lattice thermal
conductivity limits the $ZT$ values for undoped hH alloys. Recently, $p$–type
FeNbSb,Fu  _et al._ (2015), ZrCoBi,Zhu  _et al._ (2018) and TaFeSbZhu  _et
al._ (2019)–based hH compounds have been explored as TE materials with high
$ZT$ $\sim$ 1.5. This high value of $ZT$ is due to the reduction of thermal
conductivity with the help of dopants which cause more phonon scattering,
while maintaining high carrier mobility. Many other experimentally and
theoretically identified hH alloys show interesting TE properties.Kim  _et
al._ (2007); Yu  _et al._ (2009); Chen  _et al._ (2013); Fu  _et al._ (2014);
Rai  _et al._ (2015); Zeeshan  _et al._ (2017a, b, 2018); Singh  _et al._
(2019); Ning  _et al._ (2020); Jaishi  _et al._ (2021); Hori  _et al._ (2020)
Although doping has been the most widely used method to improve $ZT$ and power
factor of hH alloys, application of pressure has also been shown to be
fruitful in improving the TE properties.Kumarasinghe and Neophytou (2019);
Ning  _et al._ (2020) At present, it is well known that hH alloys are
excellent candidates for use in TE devices, but some of them require expensive
dopants to simultaneously tune the carrier concentration and suppress the
lattice thermal conductivity. Therefore, further research is required to
identify efficient, non-toxic, low-cost, and naturally abundant
multifunctional materials.

Apart from TE, the hH alloys have drawn significant attention due to the
possibility of tunable band gap,Kumarasinghe and Neophytou (2019) high optical
absorption, optical conductivity in the broad range of incident spectrum
etc.Wei and Wang (2018); Kieven  _et al._ (2010) The noticeable band gap
becomes useful in optoelectronic applications, and they further act as a
potential source in replacement of Pb containing alloys.Eperon  _et al._
(2015); Ramasamy  _et al._ (2016); Akkerman  _et al._ (2016)

From our literature survey, we noticed very few theoretical works on rhodium
based hH alloys.Anissa  _et al._ (2019); Bamgbose  _et al._ (2021); Azin-
Sanjabod  _et al._ (2021) The reported work is preliminary and needs extensive
investigation to understand the electronic, optical, mechanical, and
thermoelectric properties of these materials. With this motivation, we explore
RhTiZ (where Z are P, As, Sb, and Bi) which consists of two transition
elements - Rh and Ti. The localized $d$ orbitals of these elements are
expected to contribute heavy bands near the Fermi level ($E_{\textrm{F}}$)Yang
_et al._ (2008) which can be beneficial for TE properties as their Seebeck
coefficient usually increases with increase in the effective mass of the
charge carriers.

##  II COMPUTATIONAL DETAILS

The hH compounds RhTiZ (Z = P, As, Sb, Bi) belong to the cubic MgAgAs$-$type
structure with space group $F\bar{4}3m$. The crystal structure of RhTiZ is
shown in Figure 1. In this structure Ti and $Z$ form rock salt structure, with
$Z$ and either of the Rh or Ti forming zinc-blende structure. The Rh, Ti, and
$Z$ atoms occupy the Wyckoff position 4c (1/4, 1/4, 1/4), 4b (1/2,1/2,1/2),
and 4a (0,0,0) respectively.

Figure 1: The crystal structure of cubic hH RhTiZ (Z = P, As, Sb, Bi). The
balls in red, blue, and black color represent Rh, Ti, and $Z$ atoms
respectively. The structure was drawn using VESTA.Momma and Izumi (2011)

The structure optimization, electronic structure, mechanical stability, and
optical property calculations are performed using the linearized augmented
plane-wave (LAPW) method based on density functional (DF) theory as
implemented in the WIEN2k code.Blaha  _et al._ (2001) We used the basis set
cut-off parameter $RK_{max}$ = 7, and the muffin-tin radii (RMTs) were taken
in the ranges 2.0–2.5 Bohr for all the atoms. The standard generalized
gradient approximation (GGA)Perdew  _et al._ (1996) along with the modified
Becke$-$Johnson (mBJ) potentialTran and Blaha (2009) were used. A well
converged $k-$mesh of $12\times 12\times 12$ was used for the self-consistent
calculations and structure optimization. A denser $k-$mesh of $50\times
50\times 50$ was used for the optical and transport properties calculations.
To check the phonon stability, we performed the phonon band structure
calculation using Quantum ESPRESSO (QE) packageGiannozzi  _et al._ (2009)
based on DF perturbation theory on a 4$\times$4$\times$4 q–mesh with the norm-
conserving Perdew-Zunger exchange-correlation functional. The electronic
transport properties are calculated by solving the semi-classic Boltzmann
transport equation (BTE). The calculation is performed by using the BoltzTraP2
code.Madsen  _et al._ (2018) Though preliminary work was reported
earlier,Azin-Sanjabod  _et al._ (2021) $\kappa_{l}$ has not been considered
properly involving lattice vibrations and phonon dispersion. In order to
determine the accuracy of $ZT$, values of $\kappa$ is essential. Thus, the
lattice thermal conductivity, $\kappa_{l}$, is evaluated here by solving the
linearized BTE within the single-mode relaxation time approximation (SMA)
using thermal2 code.Fugallo  _et al._ (2013, 2014); Cepellotti  _et al._

##  III RESULTS AND DISCUSSION

###  I Structure optimization, mechanical and phonon stability

The lattice constants obtained by optimizing the volume and corresponding band
gaps using GGA and mBJ functionals are listed in Table 1. The energy vs.
volume (E–V) plot for the RhTiP, RhTiAs, RhTiSb, and RhTiBi are shown in
Figure 2. The optimized volumes of the unit cells are obtained by using
Murnaghan equation of state.Murnaghan (1944) The optimized lattice constants
are found consistent with the previous theoretical results.Ma  _et al._ (2017)

Table 1: The optimized lattice constant $a$ and the band gap $E_{g}$ within GGA and mBJ for the cubic hH alloys RhTiZ (Z = P, As, Sb, Bi). |  |  |  |  |  |  | GGA |  | mBJ |
---|---|---|---|---|---|---|---|---|---|---
System |  |  |  | a (Å) |  |  | Eg(eV) |  | Eg(eV) |
RhTiP |  |  |  | 5.76 |  |  | 0.84 |  | 1.01 |
RhTiAs |  |  |  | 5.90 |  |  | 0.85 |  | 1.06 |
RhTiSb |  |  |  | 6.14 |  |  | 0.77 |  | 0.94 |
RhTiBi |  |  |  | 6.27 |  |  | 0.73 |  | 0.94 |
Figure 2: E–V diagram for RhTiP, RhTiAs, RhTiSb, and RhTiBi. Here E0 is the
ground state energy corresponding to optimized volume.

Mechanical stability of the crystal is checked by calculating the elastic
constants on the basis of the well-known Born stability criteriaBorn (1940):
C11-C12 $>$ 0, C11+2C12 $>$ 0, C44 $>$0\. The independent elastic constants
C11, C12, and C44 for cubic hH compounds calculated by using the finite strain
theoryMurnaghan (1937) and are presented in the Table 2. The calculated values
are positive and thus shows that RhTiZ is mechanically stable. Using the
second order elastic constants, bulk modulus ($B$), shear modulus ($G$),
Young’s modulus ($Y$), Pugh’s ratio (i.e., the ratio of the bulk and shear
modulus) ($B/G$), and Poisson ratio ($\nu$) of the hH alloys are calculated as
tabulated in Table 2.

Table 2: Calculated elastic constants Cij, bulk modulus (B), shear modulus (G), Pugh’s ratio (B/G), Young’s modulus (Y), Poisson ratio ($\nu$), hardness (H), volume (V), density ($\rho$), transverse sound velocity ($\nu_{t}$), longitudinal sound velocity ($\nu_{l}$), average sound velocity ($\nu_{m}$), and Debye temperature ($\Theta_{D}$) for the cubic hH alloys RhTiZ (Z = P, As, Sb, Bi). The unit of $B$, $G$, $Y$, and H is GPa. |  | RhTiP | RhTiAs | RhTiSb | RhTiBi |
---|---|---|---|---|---|---
$C_{11}$ |  | 201.41 | 305.85 | 243.93 | 250.72 |
$C_{12}$ |  | 160.87 | 86.71 | 96.67 | 66.08 |
$C_{44}$ |  | 107.62 | 176.03 | 137.78 | 123.68 |
$B$ |  | 174.38 | 159.76 | 145.75 | 127.63 |
$G$ |  | 56.09 | 145.55 | 107.15 | 110.01 |
$B/G$ |  | 3.12 | 1.09 | 1.36 | 1.16 |
$Y$ |  | 151.99 | 334.94 | 258.18 | 156.37 |
$\nu$ |  | 0.35 | 0.15 | 0.20 | 0.17 |
$H$ |  | 5.629 | 33.979 | 21.515 | 14.701 |
$V$ (Å3) |  | 191.55 | 205.51 | 232.15 | 247.15 |
$\rho$(g/cm3) |  | 6.30 | 7.29 | 7.79 | 9.67 |
$\nu_{t}$ (m/s) |  | 2983.51 | 4467.21 | 3707.23 | 3373.41 |
$\nu_{l}$ (m/s) |  | 6288.25 | 6965.06 | 6084.34 | 5326.89 |
$\nu_{m}$ (m/s) |  | 3356.55 | 5906.76 | 4094.82 | 3710.74 |
$\Theta_{D}$ (K) |  | 397.07 | 567.00 | 454.34 | 403.22 |
Figure 3: Variation of the Pugh’s ratio of the cubic RhTiZ (Z = P, As, Sb, Bi)
hH materials. The horizontal dashed line separates the ductile materials from
brittle.

The bulk modulus is a measure of degree of stiffness of the material. The
calculated values of bulk modulus of hH alloys lies between 127.63 GPa (for
RhTiBi) to 174.38 GPa (for RhTiP). Upon substitution of different atoms (say
lighter P to heavier Bi) at the $Z-$sites, the value of bulk modulus is found
to decrease. This indicates an increase in flexibility of RhTiZ as $Z$ is
replaced by P, As, Sb, and Bi respectively. The shear and Young’s modulii are
maximum for RhTiAs. The Pugh’s ratio and Poisson ratio determines the
ductility and brittleness of the material. The critical value of Pugh’s ratio
and Poisson ratio are 1.75 and 0.26, Frantsevich  _et al._ respectively, that
separates the ductile material from the brittle one. If the value of Pugh’s
ratio is greater than 1.75 or when the Poisson ratio is greater than 0.26, the
material is ductile or else the material is brittle in nature. Under this
criteria, RhTiP is a ductile material while the other three hH alloys are
brittle in nature. This property not only measures the strength in response to
mechanical shocks but also ease in manufacturing. The enhancement of
mechanical strength makes it more durable. The plot of Pugh’s ratio for these
four compounds are shown in Figure 3. We further consider the hardness ($H$)
which can be derived by using Poisson’s ratio and Young’s modulus. As seen
from Table 2, RhTiAs has the highest hardness while RhTiP has the least. Here,
we also calculate the longitudinal and transverse acoustic wave velocities,
average sound velocity and Debye temperatureSchreiber  _et al._ (1973);
Anderson (1963) shown in Table 2. This information is important for the
synthesis and design of these materials. The $\Theta_{D}$ obtained at 0 K are
397.07, 567.00, 454.34, and 403.22 for RhTiP, RhTiAs, RhTiSb, and RhTiBi,
respectively. The calculated melting temperature Fine  _et al._ (1984) of the
system ranges from 1740 – 2360 ($\pm$ 300) K. The obtained values of
$\nu_{m}$, $\Theta_{D}$, and high value of melting temperature shows that
these hH alloys will have lower lattice thermal conductivity and high
temperature stability. These features indicates that RhTiZ can be a suitable
for TE applications.

Figure 4: Phonon dispersion curves for (a) RhTiP, (b) RhTiAs, (c) RhTiSb, and
(d) RhTiBi.

The calculated phonon dispersion of hH materials are shown in Figure 4. This
measures the dynamical stability of the materials. The absence of phonon modes
with imaginary frequencies suggests that these alloys are dynamically stable.
As hH alloys are composed of three types of atoms, it comprises a total of
nine phonon modes: three acoustic modes at lower-frequency region and six
optical modes at the high-frequency region. The maximum acoustic phonon mode
of RhTiP lies around 200 cm-1 which decreases as the element P is replaced by
As, Sb, or Bi. In case of RhTiBi the acoustic modes lie within frequency range
$\sim$ 100 cm-1. This indicates that phonon group velocity are in the order:
RhTiP $>$ RhTiAs $>$ RhTiSb $>$ RhTiBi.

###  II Electronic structure

The calculated band structure of RhTiZ (Z = P, As, Sb, Bi) along the high
symmetry directions in the first Brillouin zone are shown in Figure 5. The
materials are found to be semiconducting with indirect band gaps ranging from
$0.94-1.01$ eV within mBJ calculations. The values are comparable with the
earlier works. Anissa  _et al._ (2019); Bamgbose  _et al._ (2021); Azin-
Sanjabod  _et al._ (2021) The band gap variation is small between the valence
band maximum (VBM) and conduction band minimum (CBM). To note is, the band
next to CBM moves close to it when P is replaced by As, Sb, or Bi (see Figure
5). The reason for this is due to the extended nature of the band from the
heavier elements.

Our calculation reveals that these compounds are indirect band gap
semiconductors having CBM at $X$\- high symmetry point and valence band
maximum (VBM) at $\Gamma$. The VBM at $\Gamma$ comprises of light and heavy
band with 3-fold degeneracy shown in Figure 5. For the TE properties, this
feature plays significant role. Here, we also calculated the effective mass of
the charge carriers and the calculated values of effective mass are shown in
the Table 3. To know more about the nature of Fermi surface in RhTiZ, we
computed the iso-energy surfaces around VBM and CBM as shown in Figure
5(e)-(g). Figure 5(e) represents the iso-energy surface for RhTiBi at 4 meV
below the VBM. Similar features are observed also for RhTiP, RhTiAs, and
RhTiSb. Figure 5(f) represents the iso-energy surface plot for RhTiAs at 4 meV
above the CBM, whose features are found similar to RhTiP and RhTiSb. Unlike
others, the iso-energy surface of RhTiBi near CBM is noticeably different.
This results the band shifts from $E_{\textrm{F}}$ significantly towards the
CBM in RhTiBi. The additional energy surface observed in RhTiBi increasess the
electronic density of states (DOS) (see Figure 6(a)). This is expected to
increase the thermopower i.e. the magnitude of Seebeck coefficient. Thus, we
expect to get comparatively better values of Seebeck coefficient for RhTiBi.

Figure 5: Calculated electronic band structure for (a) RhTiP, (b) RhTiAs, (c)
RhTiSb, and (d) RhTiBi within mBJ and iso-energy surfaces (e) at 4 meV below
VBM for RhTiBi, (f) at 4 meV above CBM for RhTiAs, and (g) at 4 meV above CBM
for RhTiBi. In (f) and (g), blue and violet colors are the inner surface,
whereas yellow and green colors are the outer surface. Figure 6: (a) Total DOS
for RhTiZ (Z = P, As, Sb, Bi) and (b)-(d) Partial DOS of Rh-4$d$, Ti-3$d$, and
P-3$p$ orbitals in RhTiP within mBJ. Vertical dotted line represent $E_{\rm
F}$. Table 3: The effective mass of charge carriers in valence band and
conduction band of RhTiZ (Z = P, As, Sb, Bi) in unit of electronic rest mass
(me). The notation ($\Gamma$-K$\Gamma$) represents the effective mass of
charge carriers at $\Gamma$ point along $\Gamma$K direction.

Effective mass (m∗/me) High-symmetry point RhTiP RhTiAs RhTiSb RhTiBi Valence
Band ($\Gamma$-K$\Gamma$) HH1 5.29 5.04 6.62 4.87 HH2 0.66 0.72 0.86 0.85 LH
0.31 0.30 0.35 0.30 Conduction band (X-$\Gamma$X) Electron 1.03 1.21 1.49 1.82
Conduction band (X-XL) Electron 0.62 0.66 0.68 0.61

Figure 6 represents the total and partial DOS. The region below the
$E_{\textrm{F}}$ is found to be contributed mainly by the Rh$-4d$ ( $d-t_{2g}$
and $d-e_{g}$ orbitals), hybridizing with the Ti$-3d$ and Z$-p$ states, while
above $E_{\textrm{F}}$ the major contributions are from the Ti$-3d$ and Z$-p$
states. This suggests that the carrier concentration can be tuned near
$E_{\textrm{F}}$ by doping on to the Z–site.

###  III Optical Properties

To explore the optical properties, we calculated the real and imaginary parts
of dielectric function, refractive index, reflectivity, absorption spectra and
optical conductivity of the hH systems up to the photon energy of 30 eV to
point out the response of the materials in solar and high energy radiations.
The real and imaginary parts of the dielectric function, represented by
$\epsilon_{1}$ and $\epsilon_{2}$ are shown in Figure 7. For all the four
compounds, there are two major peaks within the photon energy 0 – 5 eV. The
magnitude of the peaks are comparable in all four cases. The obtained static
dielectric functions $\epsilon_{1}(0)$ for RhTiP, RhTiAs, RhTiSb, and RhTiBi
are 17.3, 18.0, 18.8, and 21.2, respectively. The imaginary part of the
dielectric function shows similar trend for all the alloys. For photon
energies higher than 25 eV, the imaginary part of dielectric function tends to
zero. This indicates that the radiations with higher energies transmit through
the materials without significant loss of energy.

Figure 7: Calculated real and imaginary parts of dielectric function
($\epsilon_{1}$ and $\epsilon_{2}$) as function of incident photon energies
for (a) RhTiP, (b) RhTiAs, (c) RhTiSb, and (d) RhTiBi within mBJ. Figure 8:
Calculated (a) refractive index, (b) optical absorption, (c) reflectivity, and
(d) optical conductivity as function of incident photon energies for RhTiZ (Z
= P, As, Sb, Bi ) within mBJ.

Figure 8 shows the optical response of RhTiZ. The real part of the refractive
index $n(\omega)$ is found to show similar trend as that of dielectric
function. The static refractive index $n(0)$ for the studied hHs are found to
lie within $4.2-4.6$, while the maximum values of the energy dependent
refractive index obtained are 6.35 at photon energy 2.4 eV for RhTiP, 6.36 at
2.1 eV for RhTiAs, 6.40 at 1.8 eV for RhTiSb, and 6.83 at 1.7 eV for RhTiBi,
respectively.

It is well-known that the optical absorption coefficient signifies the
thickness of the material up to which the incident light can penetrate before
being completely absorbed. It gives information on the solar energy conversion
efficiency which is important for application as optoelectronic devices. The
calculated absorption spectra is shown in Figure 8 (b). The absorption peaks
are found around 2 eV. Despite the replacement of P by other elements, the
effect on the magnitude and position of the absorption spectra are negligible.

In the visible region, RhTiP has the smallest reflectivity out of the four
compounds under study, and is found to increase as P is replaced by As, Sb,
and Bi (see Figure 8 (c)). These alloys have the maximum reflectivity for
photon energies within 3–5 eV. In addtion, the real part of the optical
conductivity have been calculated (see Figure 8 (d)). At low energy region
(i.e., $<$ 8 eV), high optical conductivity has been noticed. The values are
comparable with other hH alloys.Anissa  _et al._ (2019); Dey  _et al._ (2020)

Figure 9: Variations of Seebeck coefficient with chemical potential for (a)
RhTiP, (b) RhTiAs, (c) RhTiSb, and (d) RhTiBi at different temperatures. The
vertical dotted lines are the respective middle of band gap, and the solid
vertical lines represent the chemical potential at 300 K.

###  IV Transport Properties

The suitability of a material for TE purpose depends on properties like
Seebeck coefficient ($\alpha$), electrical conductivity ($\sigma$) and thermal
conductivity ($\kappa$). A factor that improves one of these properties can
adversely affect another. For instance, a heavy conduction or valence band is
good for Seebeck coefficient but usually results in smaller electrical
conductivity. Therefore, it is necessary to optimize all these properties to
maximize the figure of merit $ZT$. In this work, we have studied $\alpha$,
$\sigma$, and $\kappa$ of considered hH compounds using constant relaxation
time approximation (CRTA) and rigid band approximation (RBA).

The magnitude of $\alpha$ is dependent upon the density of state effective
mass, which is directly proportional to the energy valley degeneracy in the
band structure and the effective mass of the charge carriers. To estimate the
values of $\alpha$ for the RhTiZ compounds, we calculated the effective mass
of charge carriers using the electronic band structure. The effective mass is
decided by the curvature of the energy band. The flat bands in the band
structure corresponds to larger effective mass. In case of RhTiZ, VBM is
threefold degenerate and lies at the $\Gamma$. The CBM lies at the high
symmetric point $X$. An interesting observation is that, as we move down the
group from P to Bi, the band above the CBM shifts down, and in case of RhTiBi,
the conduction band is also twofold degenerate at $X$ point. We have estimated
the effective mass of holes at the $\Gamma$ along $K\Gamma$, and the effective
mass of electrons at $X$ along $XL$ and $\Gamma X$. The calculated effective
masses are listed in Table 3. Here, the notation $\Gamma-K\Gamma$ implies that
the effective mass was calculated at $\Gamma$ along the $K\Gamma$. The valence
bands are labeled as HH1, HH2 and LH as shown in Figure 5(a). The effective
mass of holes in the heaviest valence band HH1 is much larger than the
effective mass of electrons in the conduction band. The effective mass for LH
bands are comparable among four hH alloys. The maximum value of effective mass
for holes in HH1 at $\Gamma$ along K$\Gamma$ is 6.62 me for RhTiSb followed by
RhTiP (5.29 me ), RhTiAs (5.04 me), and RhTiBi (4.87 me), respectively.
Furthermore, the maximum value of effective mass for electron at $X$ along
$\Gamma X$ is 1.82 me for RhTiBi, whereas the lowest value achieved at $X$
along $\Gamma X$ is 1.03 me for RhTiP.

From the above calculation on the effective mass, the maximum values noted are
for holes. This suggests that holes will contribute largely to thermopower.
Thus, these materials are expected to have positive value of $\alpha$. It
should however be noted that the chemical potential $\mu$ for an intrinsic
semiconductor depends on the temperature and the effective masses of holes and
electrons following the relation:Kittel (2004)

$\mu=\epsilon_{\nu}+\frac{1}{2}E_{g}+\frac{3}{4}\kappa_{B}Tln(m_{v}/m_{c})$
(1)

where $\epsilon_{\nu}$ is the energy of the electron at the VBM, $m_{v}$ and
$m_{c}$ are the effective mass of the charge carriers in the valence band and
conduction band respectively, and $E_{g}$ is the energy band gap.

Figure 10: Variations of Seebeck coefficient with temperature for RhTiP,
RhTiAs, RhTiSb, and RhTiBi. Figure 11: Variations of power factor with
chemical potential for (a) RhTiP, (b) RhTiAs, (c) RhTiSb, and (d) RhTiBi at
different temperatures. The vertical dotted lines at positive (negative) $\mu$
corresponds to the optimal value for $ZT$ of $n$–type ($p$–type) doping cases,
and the vertical line at $\mu$=0 represent middle of the band gap.

Using the effective mass for HH1 band and the effective mass of electron along
$\Gamma X$ direction from Table 3 for $m_{v}$ and $m_{c}$ respectively, we
calculated the value of chemical potential $\mu$ at 300 K. In all the cases,
the chemical potential is shifted towards the conduction band. The chemical
potential at 300 K was 31 meV, 28 meV, 29 meV, and 19 meV above the middle of
the band gap for RhTiP, RhTiAs, RhTiSb, and RhTiBi, respectively. The
variation of $\alpha$ with $\mu$ at different temperatures are shown in Figure
9. The dashed (solid) vertical lines represents the middle of band gap
(chemical potential at 300 K). The calculated values of $\alpha$ in the
temperature range 400 – 1200 K for a fixed $\mu$ for 300 K are shown in Figure
10. All the hH alloys have negative Seebeck coefficients that decrease in
magnitude as temperature is increased. The largest values of $\alpha$ are
observed for RhTiBi, which can be attributed to the two-fold degeneracy of the
conduction band at the CBM in the bandstructure of RhTiBi. At 400 K, the
values of $\alpha$ are -727 $\mu$VK-1, -896 $\mu$VK-1, -1030 $\mu$VK-1, and
-1000 $\mu$VK-1for RhTiP. RhTiAs, RhTiSb, and RhTiBi, respectively. The
decrease in the magnitude of $\alpha$ with increase in temperature could be
due to the rise in the carrier concentration.

To determine the optimum doping level for enhancing the TE performance, the
power factor (PF) has been calculated for all four systems. Due to the
unavailability of experimental value of relaxation time of charge carriers,
$\tau$, we have used the CRTA and fixed the value of $\tau$ at 2$\times$10-15
s. This value was chosen based on our previous work,Jaishi  _et al._ (2021)
where we estimated $\tau$ for NiTiSn by using the experimental values of
electrical conductivity from Kim $et$ $al$.Kim  _et al._ (2007) The variations
in the PF with respect to $\mu$ for different temperatures are shown in Figure
11. Maxima are noted in the values of PF for both positive and negative values
of $\mu$. Positive region of $\mu$ represents the $n$–type doping, while the
negative represents $p$–type doping. Up to the temperature of 600 K, the
maximum values of PF for $n$–type and $p$–type are comparable. The maximum
values of PF obtained at 600 K for $p$ ($n$)–type are 15.16(13.20), 15.08
(16.87), 14.48 (15.03), and 14.27 (16.65) $\mu$Wcm-1K-2 for RhTiP, RhTiAs,
RhTiSb, and RhTiBi, respectively.

These values are comparable to the reported PF for some cobalt–based hH
materials. Zeeshan  _et al._ (2017b) For temperatures larger than 600 K, the
maximum PF values achievable by $p$–type doping are significantly larger than
those obtainable by $n$–type doping. Hence, for high temperature region (T >
600 K), $p$–type doping is more preferable for enhancing the value of PF.

Figure 12: (a) Calculated electronic thermal conductivity ($\kappa_{e}$), (b) Calculated lattice thermal conductivity ($\kappa_{l}$), and the total thermal conductivity ($\kappa_{total}$) as a function of temperature for RhTiP, RhTiAs, RhTiSb, and RhTiBi. Figure 13: Variations of $ZT$ with chemical potential for (a) RhTiP, (b) RhTiAs, (c) RhTiSb, and (d) RhTiBi at 1200 K. The upper horizontal axis denotes the number of electrons per unit cell (uc), and the negative value corresponds to hole doping. Table 4: Calculated $n$–type and $p$–type doping levels and the corresponding Seebeck coefficient, electrical conductivity, PF, thermal conductivity, and $ZT$ of RhTiZ ( Z= P, As, Sb, Bi) in cubic symmetry $F\bar{4}3m$ at 1200 K. Here, negative values of ne correspond to $p$–type doping, and the positive ne correspond to $n$–type doping. |  | RhTiP | RhTiAs | RhTiSb | RhTiBi
---|---|---|---|---|---
$n$–type | n (e/uc) | 0.15 | 0.20 | 0.21 | 0.24
$\alpha$ ($\mu$VK-1) | -183.25 | -200.40 | -202.20 | -228.61
$\sigma$$\times$104 ($\Omega^{-1}$m-1) | 8.22 | 7.53 | 6.00 | 4.32
$\alpha^{2}\sigma$ ($\mu$Wcm-1K-2) | 27.60 | 30.24 | 24.51 | 22.60
$\kappa$ (Wm-1K-1) | 8.81 | 7.85 | 5.19 | 2.96
$ZT$ | 0.38 | 0.46 | 0.58 | 0.91
$p$–type | n (e/uc) | -0.10 | -0.10 | -0.10 | -0.074
$\alpha$ ($\mu$VK-1) | 188.08 | 191.54 | 210.58 | 232.71
$\sigma$ ($\Omega^{-1}$m-1) | 9.57 | 9.31 | 7.27 | 5.24
$\alpha^{2}\sigma$ ($\mu$Wcm-1K-2) | 33.85 | 34.15 | 32.22 | 28.39
$\kappa$ (Wm-1K-1) | 9.32 | 8.69 | 5.76 | 3.48
$ZT$ | 0.44 | 0.47 | 0.67 | 0.98
Figure 14: The temperature-dependent variation of $ZT$ values for $p$–type and
$n$–type RhTiP, RhTiAs, RhTiSb, and RhTiBi.

The optimum $p$–type ($n$–type) doping levels for RhTiP, RhTiAs, RhTiSb, and
RhTiBi are 0.30 (0.28), 0.28 (0.31), 0.29 (0.34), and 0.27 (0.40) holes
(electrons) per unit cell, respectively.

In all cases, the optimum PF values for $p$–type doping are significantly
larger than that for $n$–type doping. The additional advantage of $p$–type
doping is that, the better value of PF is achievable by smaller doping
concentrations. However, the doping levels with optimum PF do not necessarily
correspond to optimum $ZT$ values.

The thermal conductivity of the materials at different temperatures are shown
in Figure 12.

The calculated values of the electronic thermal conductivity $\kappa_{e}$
(lattice thermal conductivity $\kappa_{l}$) are 2.30 (7.02), 2.17 (6.51), 1.60
(4.16), and 1.15 (2.31) Wm-1K-1 for RhTiP, RhTiAs, RhTiSb, and RhTiBi
respectively at 1200 K (Figure 12(a and b)). The values of $\kappa_{e}$
reported here corresponds to the optimum $ZT$ with $p$–type doping. As can be
seen in Figure 12(b), $\kappa_{l}$ is found to decrease with increase in
temperature for all cases. Note that the values reported here for $\kappa_{l}$
are calculated for undoped RhTiZ. These values were used for the calculation
of $ZT$ for different doping levels. Among the studied alloys, the
$\kappa_{l}$ value for RhTiP is highest, and is found to decrease as P is
replaced by As, Sb, and Bi. This might be due to the scattering of phonons
from the heavier atoms.Fu  _et al._ (2015, 2016) Due to very small value of
$\kappa_{l}$, RhTiBi has the largest $ZT$ value among the studied alloys.

Figure 13 shows the variation of $ZT$ with the chemical potential and doping
levels at 1200 K. The doping level corresponding to optimum PF (see Figure11)
are different when compared to those corresponding to optimum $ZT$ values (see
Figure 13). The different values of TE parameters corresponding to $n$ and
$p$–type doping at 1200 K are summarized in the Table 4. The variation of $ZT$
with temperature for $n$ and $p$–type doping are shown in Figure 14. The
optimum $ZT$ values of $p$–type ($n$–type) RhTiP, RhTiAs, RhTiSb, and RhTiBi
are 0.44 (0.38), 0.47 (0.46), 0.67 (0.58), and 0.98 (0.91), respectively at
1200 K. These values of $ZT$ are lower than those for the state-of-the-art TE
materials,Snyder and Toberer (2008); Fu  _et al._ (2015); Zhu  _et al._ (2019)
but are comparable to several other works done on Heusler alloys.Kim  _et al._
(2007); Jaishi  _et al._ (2021); Hu  _et al._ (2020) Thus, these Rh-based hH
alloys are suitable candidates for use as both $p$–type and $n$–type elements
in TE devices.

##  IV CONCLUSIONS

We performed the density functional theory calculations to explore the
structural, electronic, mechanical, optical, and thermoelectric properties of
rhodium based hH compounds RhTiP, RhTiAs, RhTiSb and RhTiBi. All these alloys
are found to be semiconducting in nature. In case of RhTiBi, an extra electron
pocket is observed at the CBM. The Pugh’s ratio calculation indicates that
RhTiP is ductile while the remaining alloys are brittle in nature. Optical
study reveals that RhTiZ (where Z are P, As, Sb, and Bi) have high absorption
coefficient and optical conductivity, due to which these alloys are predicted
to be suitable for optoelectronic applications. We calculated the power factor
which suggests that $p$–type doping is more favorable to enhance the TE
properties. The obtained optimum $ZT$ values are 0.38, 0.46, 0.58 and 0.91 for
$n$–type and 0.44, 0.47, 0.67, and 0.98 for $p$–type RhTiP, RhTiAs, RhTiSb and
RhTiBi, respectively at 1200 K. This indicates that the alloys are good
candidates for use in TE devices. Brittleness can be an issue during synthesis
and use in TE devices in few systems. To address this, doping and annealing
may be possible ways to improve ductility. Furthermore, the value of $\tau$
used in this study might be larger than the actual value of the relaxation
time. In that case, the $ZT$ values might be smaller than the values
calculated here. Hence, an experimental study of the transport properties of
these materials can be a very interesting prospect for experimentalists.

##  V acknowledgments

D. R. Jaishi thanks Rajendra Adhikari for fruitful discussion and suggestions.
Most of the computational works was performed at the Advanced Materials
Research Laboratory (AMRL), Central Department of Physics, Kirtipur, Tribhuvan
University supported by Alexander von Humboldt Foundation, Germany as
equipment grants.

## References

  * Hu (2012) X. Hu, Adv. Mater. 24, 294 (2012).
  * Leskelä  _et al._ (2019) M. Leskelä, M. Mattinen, and M. Ritala, J. Vac. Sci. Technol. B 37, 030801 (2019).
  * Ghimire  _et al._ (2010) M. P. Ghimire, Sandeep, and R. K. Thapa, Mod. Phys. Lett. B 24, 2187 (2010).
  * Ghimire  _et al._ (2015) M. P. Ghimire, R. K. Thapa, D. P. Rai, Sandeep, T. P. Sinha, and X. Hu, J. Appl. Phys. 117, 063903 (2015).
  * Kumar  _et al._ (2020) N. Kumar, S. N. Guin, K. Manna, C. Shekhar, and C. Felser, Chem. Rev. 121, 2780 (2020).
  * Ghimire and Richter (2017) M. P. Ghimire and M. Richter, Nano Lett. 17, 6303 (2017).
  * Roknuzzaman  _et al._ (2017) M. Roknuzzaman, K. Ostrikov, H. Wang, A. Du, and T. Tesfamichael, Sci. Rep. 7, 1 (2017).
  * Fitriani  _et al._ (2016) Fitriani, R. Ovik, B. D. Long, M. C. Barma, M. Riaz, M. F. M. Sabri, S. M. Said, and R. Saidur, Renew. Sustain. Energy Rev. 64, 635 (2016).
  * Rowe (2018) D. M. Rowe,  _Thermoelectrics Handbook_ (CRC Press, 2018).
  * Snyder and Toberer (2008) G. J. Snyder and E. S. Toberer, Nat. Mater. 7, 105 (2008).
  * Sootsman  _et al._ (2009) J. R. Sootsman, D. Y. Chung, and M. G. Kanatzidis, Angew. Chem. Int. Ed. 48, 8616 (2009).
  * Tritt (2000) T. Tritt,  _Recent trends in thermoelectric materials research, part two_ (Academic Press, 2000).
  * Poudel  _et al._ (2008) B. Poudel, Q. Hao, Y. Ma, Y. Lan, A. Minnich, B. Yu, X. Yan, D. Wang, A. Muto, D. Vashaee, X. Chen, J. Liu, M. S. Dresselhaus, G. Chen, and Z. Ren, Science 320, 634 (2008).
  * Mao  _et al._ (2017) J. Mao, Y. Wang, Z. Liu, B. Ge, and Z. Ren, Nano Energy 32, 174 (2017).
  * Zhao  _et al._ (2016) L.-D. Zhao, X. Zhang, H. Wu, G. Tan, Y. Pei, Y. Xiao, C. Chang, D. Wu, H. Chi, L. Zheng, S. Gong, C. Uher, J. He, and M. G. Kanatzidis, J. Am. Chem. Soc. 138, 2366 (2016).
  * Pei  _et al._ (2012) Y. Pei, H. Wang, and G. J. Snyder, Adv. Mater. 24, 6125 (2012).
  * Neophytou and Kosina (2011) N. Neophytou and H. Kosina, Phys. Rev. B 83, 245305 (2011).
  * He  _et al._ (2016) R. He, D. Kraemer, J. Mao, L. Zeng, Q. Jie, Y. Lan, C. Li, J. Shuai, H. S. Kim, Y. Liu, D. Broido, C.-W. Chu, G. Chen, and Z. Ren, Proc. Natl. Acad. Sci. U.S.A 113, 13576 (2016).
  * Rogl  _et al._ (2016) G. Rogl, A. Grytsiv, M. Gürth, A. Tavassoli, C. Ebner, A. Wünschek, S. Puchegger, V. Soprunyuk, W. Schranz, E. Bauer, H. Müller, M. Zehetbauer, and P. Rogl, Acta Mater. 107, 178 (2016).
  * Rausch  _et al._ (2015) E. Rausch, B. Balke, S. Ouardi, and C. Felser, Energy Technol. 3, 1217 (2015).
  * Goldsmid and Douglas (1954) H. J. Goldsmid and R. W. Douglas, Br. J. Appl. Phys. 5, 386 (1954).
  * Gelbstein  _et al._ (2005) Y. Gelbstein, Z. Dashevsky, and M. P. Dariel, Physica B 363, 196 (2005).
  * Poudeu  _et al._ (2006) P. F. P. Poudeu, J. D’Angelo, H. Kong, A. Downey, J. L. Short, R. Pcionek, T. P. Hogan, C. Uher, and M. G. Kanatzidis, J. Am. Chem. Soc. 128, 14347 (2006).
  * Lee  _et al._ (2012) E. K. Lee, L. Yin, Y. Lee, J. W. Lee, S. J. Lee, J. Lee, S. N. Cha, D. Whang, G. S. Hwang, K. Hippalgaonkar, A. Majumdar, C. Yu, B. L. Choi, J. M. Kim, and K. Kim, Nano Lett. 12, 2918 (2012).
  * Fu  _et al._ (2015) C. Fu, S. Bai, Y. Liu, Y. Tang, L. Chen, X. Zhao, and T. Zhu, Nat. Commun. 6, 1 (2015).
  * Zhu  _et al._ (2018) H. Zhu, R. He, J. Mao, Q. Zhu, C. Li, J. Sun, W. Ren, Y. Wang, Z. Liu, Z. Tang, A. Sotnikov, Z. Wang, D. Broido, D. J. Singh, G. Chen, K. Nielsch, and Z. Ren, Nat. Commun. 9, 1 (2018).
  * Zhu  _et al._ (2019) H. Zhu, J. Mao, Y. Li, J. Sun, Y. Wang, Q. Zhu, G. Li, Q. Song, J. Zhou, Y. Fu, R. He, T. Tong, Z. Liu, W. Ren, L. You, Z. Wang, J. Luo, A. Sotnikov, J. Bao, K. Nielsch, G. Chen, D. J. Singh, and Z. Ren, Nat. Commun. 10, 1 (2019).
  * Kim  _et al._ (2007) S.-W. Kim, Y. Kimura, and Y. Mishima, Intermetallics 15, 349 (2007).
  * Yu  _et al._ (2009) C. Yu, T.-J. Zhu, R.-Z. Shi, Y. Zhang, X.-B. Zhao, and J. He, Acta Mater. 57, 2757 (2009).
  * Chen  _et al._ (2013) S. Chen, K. C. Lukas, W. Liu, C. P. Opeil, G. Chen, and Z. Ren, Adv. Energy Mater. 3, 1210 (2013).
  * Fu  _et al._ (2014) C. Fu, T. Zhu, Y. Pei, H. Xie, H. Wang, G. J. Snyder, Y. Liu, Y. Liu, and X. Zhao, Adv. Energy Mater. 4, 1400600 (2014).
  * Rai  _et al._ (2015) D. P. Rai, A. Shankar, Sandeep, M. P. Ghimire, R. Khenata, and R. K. Thapa, RSC Adv. 5, 95353 (2015).
  * Zeeshan  _et al._ (2017a) M. Zeeshan, J. van den Brink, and H. C. Kandpal, Phys. Rev. Mater. 1, 074401 (2017a).
  * Zeeshan  _et al._ (2017b) M. Zeeshan, H. K. Singh, J. van den Brink, and H. C. Kandpal, Phys. Rev. Mater. 1, 075407 (2017b).
  * Zeeshan  _et al._ (2018) M. Zeeshan, T. Nautiyal, J. van den Brink, and H. C. Kandpal, Phys. Rev. Mater. 2, 065407 (2018).
  * Singh  _et al._ (2019) S. Singh, M. Zeeshan, U. Singh, J. van den Brink, and H. C. Kandpal, J. Phys. Condens. Matter 32, 055505 (2019).
  * Ning  _et al._ (2020) S. Ning, S. Huang, Z. Zhang, R. Zhang, N. Qi, and Z. Chen, Phys. Chem. Chem. Phys. 22, 14621 (2020).
  * Jaishi  _et al._ (2021) D. R. Jaishi, N. Sharma, B. Karki, B. P. Belbase, R. P. Adhikari, and M. P. Ghimire, AIP Adv. 11, 025304 (2021).
  * Hori  _et al._ (2020) A. Hori, S. Minami, M. Saito, and F. Ishii, Appl. Phys. Lett. 116, 242408 (2020).
  * Kumarasinghe and Neophytou (2019) C. Kumarasinghe and N. Neophytou, Phys. Rev. B 99, 195202 (2019).
  * Wei and Wang (2018) J. Wei and G. Wang, J. Alloy. Compd. 757, 118 (2018).
  * Kieven  _et al._ (2010) D. Kieven, R. Klenk, S. Naghavi, C. Felser, and T. Gruhn, Phys. Rev. B 81, 075208 (2010).
  * Eperon  _et al._ (2015) G. E. Eperon, G. M. Paternò, R. J. Sutton, A. Zampetti, A. A. Haghighirad, F. Cacialli, and H. J. Snaith, J. Mater. Chem. A 3, 19688 (2015).
  * Ramasamy  _et al._ (2016) P. Ramasamy, D.-H. Lim, B. Kim, S.-H. Lee, M.-S. Lee, and J.-S. Lee, Chem. Commun. 52, 2067 (2016).
  * Akkerman  _et al._ (2016) Q. A. Akkerman, M. Gandini, F. Di Stasio, P. Rastogi, F. Palazon, G. Bertoni, J. M. Ball, M. Prato, A. Petrozza, and L. Manna, Nat. Energy 2, 1 (2016).
  * Anissa  _et al._ (2019) B. Anissa, D. Radouan, and B. Benaouda, Int. J. Mod. Phys. B 33, 1950247 (2019).
  * Bamgbose  _et al._ (2021) M. Bamgbose, P. Adebambo, G. Solola, and G. Adebayo, Mater. Sci. Eng., B 264, 114987 (2021).
  * Azin-Sanjabod  _et al._ (2021) M. R. Azin-Sanjabod, B. Nedaee-Shakarab, Y. Azizian-Kalandaragh, and H. Akbari, Int. Nano Lett. , 1 (2021).
  * Yang  _et al._ (2008) J. Yang, H. Li, T. Wu, W. Zhang, L. Chen, and J. Yang, Adv. Funct. Mater. 18, 2880 (2008).
  * Momma and Izumi (2011) K. Momma and F. Izumi, J. Appl. Cryst. 44, 1272 (2011).
  * Blaha  _et al._ (2001) P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, J. Luitz, R. Laskowsk, T. Fabien, L. Marks, and L. Marks, An augmented plane wave+ local orbitals program for calculating crystal properties (2001).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Tran and Blaha (2009) F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Condens. Matter 21, 395502 (2009).
  * Madsen  _et al._ (2018) G. K. H. Madsen, J. Carrete, and M. J. Verstraete, Comput. Phys. Commun. 231, 140 (2018).
  * Fugallo  _et al._ (2013) G. Fugallo, M. Lazzeri, L. Paulatto, and F. Mauri, Phys. Rev. B 88, 045430 (2013).
  * Fugallo  _et al._ (2014) G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, and F. Mauri, Nano Lett. 14, 6109 (2014).
  * (58) A. Cepellotti, G. Fugallo, L. Paulatto, M. Lazzeri, F. Mauri, and N. Marzari, Nat. Commun. 6.
  * Murnaghan (1944) F. Murnaghan, Proc. Natl. Acad. Sci. U.S.A 30, 244 (1944).
  * Ma  _et al._ (2017) J. Ma, V. I. Hegde, K. Munira, Y. Xie, S. Keshavarz, D. T. Mildebrath, C. Wolverton, A. W. Ghosh, and W. H. Butler, Phys. Rev. B 95, 024411 (2017).
  * Born (1940) M. Born, in  _Mathematical Proceedings of the Cambridge Philosophical Society_ , Vol. 36 (Cambridge University Press, 1940) pp. 160–172.
  * Murnaghan (1937) F. D. Murnaghan, Am. J. Math. 59, 235 (1937).
  * (63) I. N. Frantsevich, F. F. Voronov, and S. A. Bakuta, Kiev, Izdatel’stvo Naukova Dumka , 288.
  * Schreiber  _et al._ (1973) E. Schreiber, O. L. Anderson, N. Soga, _et al._ ,  _Elastic constants and their measurement_ , Vol. 6 (McGraw-Hill New York, 1973).
  * Anderson (1963) O. L. Anderson, J. Phys. Chem. Solids 24, 909 (1963).
  * Fine  _et al._ (1984) M. E. Fine, L. D. Brown, and H. L. Marcus, Scr. Metall. Mater. 18, 951 (1984).
  * Dey  _et al._ (2020) A. Dey, R. Sharma, and S. A. Dar, Mater. Today Commun. 25, 101647 (2020).
  * Kittel (2004) C. Kittel,  _Introduction to Solid State Physics_ , 8th ed. (Wiley, 2004).
  * Fu  _et al._ (2016) C. Fu, H. Wu, Y. Liu, J. He, X. Zhao, and T. Zhu, Adv. Sci. 3, 1600035 (2016).
  * Hu  _et al._ (2020) Y. Hu, Y. Jin, G. Zhang, and Y. Yan, RSC Adv. 10, 28501 (2020).


# Steady-state moments under resetting to a distribution

Kristian Stølevik Olsen  Nordita, Royal Institute of Technology and Stockholm
University, Hannes Alfvéns väg 12, 23, SE-106 91 Stockholm, Sweden
Institut für Theoretische Physik II - Weiche Materie, Heinrich-Heine-
Universität Düsseldorf, D-40225 Düsseldorf, Germany

###### Abstract

The non-equilibrium steady states emerging from stochastic resetting to a
distribution is studied. We show that for a range of processes, the steady-
state moments can be expressed as a linear combination of the moments of the
distribution of resetting positions. The coefficients of this series are
universal in the sense that they do not depend on the resetting distribution,
only underlying dynamics. We consider the case of a Brownian particle and a
run-and-tumble particle confined in a harmonic potential, where we derive
explicit closed-form expressions for all moments for any resetting
distribution. Numerical simulations are used to verify the results, showing
excellent agreement.

###### pacs:

Valid pACS appear here

††preprint: ApS/123-QED

##  I Introduction

Large fluctuations are present in almost all systems in Nature, and are
typically unavoidable in small systems. Instantaneous and large jumps in a
systems state can have drastic implications for the systems dynamics, and
potentially drive it out of thermal equilibrium. Stochastic resetting is one
example of large and sudden fluctuations, whereby a systems state is typically
brought back to its initial state at a constant rate [1, 2]. The study of
(non-equilibrium) steady states emerging from stochastic resetting has gained
significant attention in the past decade due to both its broad relevance in
diverse scientific disciplines, and the analytically tractable steady states
that emerge. Most famous is perhaps the application of resetting in search
processes, where it has been shown that resets, when optimized, can expedite a
search process [1, 3, 4, 5, 6, 7]. In biology resetting can be found across
multiple scales, ranging from protein-bound search processes [8] to the
migratory and foraging patterns of animals [9, 10]. In physics, fundamental
aspects have been in regards to the non-equilibrium nature of resetting, such
as relaxation dynamics to non-equilibrium steady states [11, 12], and
stochastic thermodynamics [13, 14, 15, 16, 17, 18, 19, 20].

The paradigmatic example of stochastic resetting in physics is that of a reset
Brownian motion [1, 21, 2, 22, 23]. Since the first studies of resetting in
its modern form more than a decade ago, a myriad of generalizations have
surfaced, including Brownian motion in potentials [24, 25], processes where
the diffusivity switches between multiple possible values [26], fractional
Brownian processes [27, 28], and in active matter models like active Brownian
and run-and-tumble particles [29, 30, 31, 32]. Steady states have also been
studied under a variety of resetting schemes, for example time-dependent
resetting rates [33], non-instantaneous resets [34], non-Poissonian waiting
times [35, 36, 34], or resetting mediated by an external trap [37] to mention
a few. For a recent review of stochastic resetting and its applications, see
[2].

Another resetting scheme of high degree of relevance to both theory and
experiments is one where the resetting position is random and drawn from a
distribution $p_{R}(x)$. Recently, it has been argued that only for resetting
distributions that are not delta-peaked does the standard path-based framework
of stochastic thermodynamics make sense for resetting systems [18]. It is also
highly natural to assume that resets cannot be performed with perfect
precision, neither in experiments or in natural systems. In addition, from a
theoretical perspective, the resetting distribution $p_{R}(x)$ introduces new
lengthscales into the problem, which can give rise to interesting phenomena.

Figure 1: Sketch of the system under consideration. A particle moving in a
potential $V(x)$ undergoes resetting with constant rate $r$ to a resetting
distribution $p_{R}(x)$. Multiple lengthscales are present in the problem,
competing to produce interesting steady-state properties.

In this paper we characterize the steady-states of resetting processes where
the reset position is drawn from a distribution $p_{R}(x_{0})$. See Fig. (1)
for a sketch of the system considered. We consider Poissonian resets, where
the waiting time between resets is distributed exponentially with constant
rate $r$. The exact steady-state distribution can in general be rather arduous
to compute for any interesting choice of $p_{R}(x_{0})$, and using it to
calculate moments and cumulants can be tricky from a practical standpoint.
Here we provide an alternate approach, where we find exact closed-form
formulas for all the moments in the steady state, assuming that the moments of
the resetting distribution $\langle x_{0}^{n}\rangle_{R}$ are known _a
priori_. The main results take the form

$\langle x^{n}\rangle_{*}=\sum_{j=0}^{n}C_{nj}\langle x_{0}^{j}\rangle_{R}$
(1)

where the asterisk denotes steady states and where the coefficients
$\\{C_{nj}\\}_{j=1}^{n}$ depend only on resetting rate and underlying system
parameters. The fact that the coefficients will not depend on the shape of the
resetting distribution makes the results rather universal, and can be used for
gain insight into the steady states under any choice of resetting distribution
once the coefficients $C_{nj}$ are known. The presented results could be
verified experimentally using optical tweezers, and be of relevance to further
theoretical studies of distributed resetting for example in thermodynamical
settings [18].

This paper is organized as follows. Section II introduces the (last) renewal
equation, and uses it to prove the main result of the paper, Eq. (1), for the
steady state moments for a class of processes. Section III derives the series
coefficients $C_{nj}$ for Brownian motion in a harmonic potential, i.e. for an
Ornstein-Uhlenbeck process. In section IV we consider a run-and-tumble
particles in a harmonic potential, which may be though of as a non-Markovian
version of the Ornstein-Uhlenbeck process, driven by a telegraphic noise.
Section V offers concluding remarks.

##  II The renewal equation approach

The standard route to find the propagator $p(x,t|x_{0})$ in the presence of
resetting is through the (last) renewal equation [2]

$\displaystyle p(x,t|x_{0})=$ $\displaystyle e^{-rt}p_{0}(x,t|x_{0})$
$\displaystyle+r\int_{0}^{t}d\tau e^{-r\tau}\int
dx_{0}p_{R}(x_{0})p_{0}(x,\tau|x_{0}).$ (2)

Here $p_{0}(x,t)$ is the propagator of the underlying system in absence of
resetting ($r=0$). The first term corresponds to trajectories where no
resetting took place (which happens with probability $e^{-rt}$). The second
term takes into account trajectories where resetting does happen, and makes
use of the fact that at every reset all dynamical memory is erased. Therefore
one only needs information regarding the last reset, which took place at time
$t-\tau$, before evolving freely for an exponentially distributed waiting time
$\tau$ from a reset position $x_{0}$ drawn from $p_{R}(x_{0})$.

In the late-time regime, the renewal equation predicts the following steady
state

$\displaystyle p_{*}(x)$ $\displaystyle=r\int_{0}^{\infty}d\tau e^{-r\tau}\int
dx_{0}p_{R}(x_{0})p_{0}(x,\tau|x_{0})$ (3) $\displaystyle=r\int
dx_{0}p_{R}(x_{0})\tilde{p}_{0}(x,r|x_{0})$ (4)

where $\tilde{p}_{0}(x,s|x_{0})$ is the Laplace transform of $p_{0}$. For many
interesting cases beyond free Brownian motion, performing the above integrals
analytically at best arduous and in the worst case impossible in closed form.
Even for simple case of freely diffusive Brownian particles with Gaussian
resetting distribution, the steady state density becomes rather involved [18].
Here we propose a method based on moments that is useful when the underlying
process is more involved, such as in the case of distributed resets for
diffusion in a potential.

The expression in Eq. (1) can be shown to follow from the renewal equation for
a certain class of processes. From the last renewal equation we know that the
steady-state moments can in general be written as

$\displaystyle\langle x^{n}\rangle_{*}$ $\displaystyle=r\int_{0}^{\infty}d\tau
e^{-r\tau}\int dx_{0}\int dxx^{n}p_{0}(x,\tau|x_{0})p_{R}(x_{0})$ (5)

To proceed, we consider a class of processes where the propagator satisfies
$p_{0}(x,\tau|x_{0})=p_{0}(x-x_{0}g(\tau),\tau)$, for some time-dependent
function $g(\tau)$ that may be specific to each process. For spatially
homogeneous processes, such as free Brownian motion, Levy flights etc, we have
$p_{0}(x,\tau|x_{0})=p_{0}(x-x_{0},\tau)$ and hence $g(\tau)=1$. For the
Ornstein-Uhlenbeck process with relaxation-time $\tau_{\text{rel}}$,
$g(\tau)=\exp(-\tau/\tau_{\text{rel}})$. Making a change of coordinates to
$y=x-x_{0}g(\tau)$ we find

$\langle x^{n}\rangle_{*}=r\int d\tau e^{-r\tau}\int dx_{0}\int
dy[y+x_{0}g(\tau)]^{n}p_{0}(y,\tau)p_{R}(x_{0})$ (6)

Using the binomial formula we can write

$[y+x_{0}g(\tau)]^{n}=\sum_{j=0}^{n}{n\choose j}y^{n-j}g^{j}(\tau)x_{0}^{j}$
(7)

Combining Eq. (6) and Eq. (7) one immediately finds Eq. (1), where the
expansion coefficients take the form

$C_{n,j}=r{n\choose j}\int dy\int_{0}^{\infty}d\tau
e^{-r\tau}y^{n-j}g^{j}(\tau)p_{0}(y,\tau)$ (8)

We see that the coefficients depends on the resetting rate, and in addition on
the propagator of the underlying process through $p_{0}(y,\tau)$ and
$g(\tau)$. In particular, they are related to the Laplace transform of the
moments of the underlying process without resetting. Importantly, the
coefficients do not depend on the resetting distribution, and therefore holds
universally for any choice of $p_{R}(x)$.

Figure 2: Steady state variance for a Brownian particle in a harmonic
potential, for Gaussian resetting distribution with mean $z$ and variance
$\sigma^{2}$. Panel a) shows resetting centered at the origin $z=0$, while
panel b) shows the case $z=3$. Dots show numerical simulations, while solid
curves correspond to Eq. (24). Parameters are set to $k=\mu=D=1$.

###  II.1 Spatially homogeneous processes

Before we move on to study cases where the above approach is useful, we
briefly mention some properties of spatially homogeneous systems. Homogeneity
sets $g(t)=1$, which will in some cases de-couple the space- and time-
integration in Eq. (8). For example, we have the following identities

$\displaystyle C_{nn}$ $\displaystyle=1$ (9) $\displaystyle C_{n-1,n-2}$
$\displaystyle=\frac{n-1}{n}C_{n,n-1}$ (10)

which follow immediately from Eq. (8). The latter identity is a consequence of
the fact that for homogeneous systems the integrand in Eq. (8) depends on $n$
and $j$ only through the combination $n-j$. This gives rise to the more
general identity

${n\choose n-k}C_{n-\ell,n-\ell-k}={n-\ell\choose n-\ell-k}C_{n,n-k}$ (11)

These properties of the series coefficients imply that the cumulants in the
steady state are linearly related to the corresponding cumulants of the
resetting distribution. For example, for the variance $\kappa_{2}$ and the
third cumulant $\kappa_{3}$, we have

$\displaystyle\kappa_{2}$
$\displaystyle=\kappa_{2}^{R}+C_{2,0}-\frac{C_{2,1}^{2}}{4}$ (12)
$\displaystyle\kappa_{3}$
$\displaystyle=\kappa_{3}^{R}+C_{3,0}-3C_{2,0}C_{1,0}+2C_{1,0}^{2}$ (13)

Hence the variance and the skew, as measured by $\kappa_{3}$, are directly
”inherited” from $p_{R}(x)$, plus terms that comes from the dynamics in
between resetting events. Another way to look at this is that the cumulants
(except the mean) are simply additive combinations of the corresponding
cumulant of the resetting distriution and the value one would obtain if the
resetting was to a fixed position $p_{R}(x)=\delta(x-x_{r})$. This is because
in homogeneous systems exactly where the particle is reset does not matter for
cumulants of higher order than $1$ due to spatial translational invariance. In
this sense, the case of spatially homogeneous processes are not all that
interesting to study under the effect of a distribution of resetting
positions, as most properties are transferred directly from $p_{R}(x)$ or the
corresponding resetting process when $p_{R}(x)=\delta(x-x_{r})$.

However, if the system is _not_ homogeneous there are bound to be some
lengthscales already present in the system. These can interfere with the
scales introduced by the resetting distribution, giving rise to more complex
steady-state properties. For more general processes, such as Brownian motion
in a harmonic potential, $g(t)\neq 1$ and the coefficients do not satisfy the
above relations Eqs. (9, 10, 11).

While Eq. (1) and (8) shows the general form of the steady-state moments,
these equations are not always convenient as a practical way to derive
$C_{n,j}$ for specific cases, as the integrals involved in Eq. (8) are not
always trivial. The approach we will take, relying on the knowledge that Eq.
(1) holds, will be to identify the coefficients $C_{nj}$ rather easily
directly from the master equation for a range of processes.

##  III Distributed resets in a harmonic potential

In this section we consider Brownian motion in a harmonic potential with
resetting to a distribution $p_{R}(x)$. During a timestep $dt$, the particle
has a probability $rdt$ to follow the underlying dynamics

$x(t+dt)=x(t)+\mu F(x)dt+\sqrt{2D}\Delta W(t),$ (14)

and a complimentary probability $1-rdt$ to be reset

$x(t+dt)=x_{r}.$ (15)

Here $x_{r}$ is the resetting position which we draw from a resetting
distribution $p_{R}(x)$ at every resetting event. In Eq. (14), $\Delta W(t)$
is the increment of a Wiener process, which has mean zero and variance
$\langle\Delta W(t)^{2}\rangle=\sqrt{dt}$. $F(x)=-V^{\prime}(x)$ is a
conservative force field originating from a potential $V(x)$, $\mu=1/\gamma$
is the mobility (inverse friction coefficient), and $D$ is the diffusivity.

The associated master equation for the density $p(x,t|x_{0},t)$ can be
obtained from noting that in a small step $dt$ there are two ways the density
can be updated [2]

$p(x,t+dt|x_{0},0)=(1-rdt)\left\langle
p(x-dx,t|x_{0})\right\rangle_{dx}+rdtp_{R}(x)$ (16)

Here the first term takes into account trajectories where the particle at time
$t$ was located at position $x^{\prime}=x-dx$ before taking a (random) step
$dx$ without resetting. This happens with probability $1-rdt$. With
probability $rdt$, the particle resets to a position drawn from $p_{R}(x)$.
Expanding to first order in $dt$, meaning second order in $dx$ due to the
properties of the Wiener increments, one finds in the $dt\to 0$ limit the
master equation

$\displaystyle\partial_{t}p(x,t)=$ $\displaystyle\partial_{x}\left[-\mu
F(x)p(x,t)+D\partial_{x}p(x,t)\right]$ $\displaystyle-rp(x,t)+rp_{R}(x)$ (17)

where we have written $p(x,t|x_{0},t_{0})=p(x,t)$ for simplicity.

The time-evolution of any generic observable $\langle f(x)\rangle$ can be
derived using this master equation. If the time evolution operator in the
absence of resetting is denoted $\hat{\mathcal{L}}$, one generally has

$\partial_{t}\langle f(x)\rangle=\int dxf(x)\hat{\mathcal{L}}p(x,t)-r\langle
f(x)\rangle+r\langle f(x)\rangle_{R}$ (18)

where the subscript $R$ denoted averages calculated using the resetting
distribution $p_{R}(x)$. If a steady state can be reached at late times
$\partial_{t}\langle f(x)\rangle_{*}=0$, one finds

$\langle f(x)\rangle_{*}=\langle
f(x)\rangle_{R}+r^{-1}\langle\hat{\mathcal{L}}^{\dagger}f(x)\rangle_{*}$ (19)

where the asterisk denotes steady states. For monomial functions $f(x)=x^{n}$
Eq. (19) together with Eq. (III) results in the following hierarchy for the
steady-state moments

$\langle x^{n}\rangle_{*}=\langle x^{n}\rangle_{R}-\frac{n\mu}{r}\langle
x^{n-1}\partial_{x}V(x)\rangle_{*}+\frac{n(n-1)D}{r}\langle
x^{n-2}\rangle_{*}$ (20)

For a harmonic potential $V(x)=\frac{1}{2}kx^{2}$ the hierarchy simplifies,
becoming

$\langle x^{n}\rangle_{*}=\frac{\langle x^{n}\rangle_{R}}{1+\frac{n\mu
k}{r}}+\frac{n(n-1)D}{r(1+\frac{n\mu k}{r})}\langle x^{n-2}\rangle_{*}$ (21)

Rather than calculating the integrals in Eq. (8), we can solve this second-
order recursive equation directly using standard methods, resulting in

$\langle x^{n}\rangle_{*}=\sum_{j\in 2\mathbf{N}}^{n}\left(\prod_{i\in
2\mathbf{N},i<j}\frac{(n-i)(n-i-1)}{1+\frac{(n-i)\mu
k}{r}}\right)\left(\frac{D}{r}\right)^{j/2}\frac{1}{1+\frac{({n-j})\mu
k}{r}}\>\langle x^{n-j}\rangle_{R}$ (22)

where the sum is over positive even numbers, $2\mathbf{N}$,(with zero
included) up to and including $n$, and the product is over positive even
numbers up to the largest even number smaller than $j$. Note that the order of
summation is opposite to that of Eq. (1), which can easily be reversed by
summing over $i=n-j$ in stead. This gives the expansion coefficients for
Brownian motion in a harmonic potential, valid for any resetting distribution
$p_{R}(x)$.

From Eq. (22) any moment or cumulant can be obtained directly. For example,
the mean takes the simple form

$\langle x\rangle_{*}=\frac{\langle x\rangle_{R}}{1+\mu k/r}$ (23)

where we see that due to the potential the steady-state mean is always less
than the mean of the resetting distribution. The variance takes the form

$\kappa_{2}=\frac{2D/r+\langle x^{2}\rangle_{R}}{1+2\mu k/r}-\frac{\langle
x\rangle_{R}^{2}}{(1+\mu k/r)^{2}}$ (24)

This variance can display non-monotonic behavior as a function of resetting
rate $r$ resetting with mean $\langle x\rangle_{R}$ non-zero. As a concrete
example, consider resetting to a Gaussian distribution
$p_{R}(x)=\mathcal{N}(x,\sigma)$ with mean $z$ and variance $\sigma^{2}$. Fig.
(2 a) shows the variance as a function of resetting rate for centered
resetting $z=0$, while (2 b) shows the same quantity for $z\neq 0$. While
resetting to a position with mean centered at the potential minimum gives rise
to monotonic increase of decrease in the variance, a non-centered resetting
$z\neq 0$ gives rise to non-monotonic variance with a global maximum at a
critical resetting rate. High values of the resetting width $\sigma$, this
non-monotonic behavior can be erased.

The skew of the steady-state also shows interesting behaviors. Using again Eq.
(22) for the third moments and combining them into the third cumulant, we find

$\displaystyle\kappa_{3}=$ $\displaystyle\frac{2r^{3}\langle
x\rangle_{R}^{3}}{(k\mu+r)^{3}}+\frac{6D\langle
x\rangle_{R}r}{(k\mu+r)(3k\mu+r)}$ $\displaystyle-\frac{3\langle
x\rangle_{R}r(2D+\langle
x^{2}\rangle_{R}r)}{(k\mu+r)(2k\mu+r)}+\frac{r}{3k\mu+r}\langle
x^{3}\rangle_{R}$ (25)

Just like the variance, the skew may also display non-monotonic behavior. More
interestingly, in some parameter regions the skew even changes sign as a
function of resetting rate; for $r<r_{c}$ the skew is positive, while for
$r>r_{c}$ the skew is negative. Returning again to the case of Gaussian
resetting $p_{R}\sim\mathcal{N}(z,\sigma)$ one can find the critical resetting
rate $r_{c}$ where the skew is zero to be

$\displaystyle r_{c}=$
$\displaystyle\frac{\sqrt{k^{3}\mu^{3}z^{2}\left(24D+k\mu\left(z^{2}-24\sigma^{2}\right)\right)}}{6\left(D-k\mu\sigma^{2}\right)}$
$\displaystyle-\frac{k\mu\left(6D+k\mu\left(z^{2}-6\sigma^{2}\right)\right)}{6\left(D-k\mu\sigma^{2}\right)}$
(26)

This zero only exists when the resetting width is sufficiently small, namely

$\sigma^{2}<\sigma_{c}^{2}\equiv\frac{D}{k\mu}+\frac{z^{2}}{24}$ (27)

where it is assumed that $z\neq 0$. For larger resetting widths
$\sigma>\sigma_{c}$ the third cumulant may still be non-monotonic, as seen in
Fig. (3 c), but will not change sign as resetting rate is varied.

Figure 3: Steady state third cumulant from Eq. (III) for an Ornstein-Uhlenbeck
process with Gaussian resetting $p_{R}\sim\mathcal{N}(z,\sigma)$, where
parameters used are $k=\mu=D=1$. For sharp resetting to a point high in the
potential landscape, $z=3$ in this case, the skew of the distribution may
change sign at a critical resetting rate given by Eq. (III). According to Eq.
(27), this will only happen for $\sigma<\sigma_{c}\approx 1.173$. Dots show
numerical simulations, while solid curves correspond to Eq. (III).

##  IV Resetting of run-and-tumble particles

This section derives the series coefficients entering into Eq. (1) for an
active run-and-tumble particle in a harmonic potential. In active particle
systems non-trivial steady states emerge when the particles are put into
confinement, either due to interactions with solid boundaries or smooth
potentials [38, 39]. Intriguing phenomena has been studied in this context in
the past, ranging from accumulation and trapping of individual particles at
confining boundaries [40, 41, 42, 43, 44, 45, 46, 47], to collective phenomena
only possible due to the presence of confinement [48, 49, 50]. Understanding
the interplay between the activity of the particles and the confinement is
crucial for elucidating the complex behavior of active matter in realistic
environments. On top of confinement, we here include the effects of
distributed stochastic resetting. In addition to a model of active matter, the
model considered here may also be interpreted as distributed resetting in a
non-Markovian Ornstein-Uhlenbeck process driven by a telegraphic rather than
Gaussian white noise.

###  IV.1 Stochastic equation and coupled master equations

We consider a RTP following the stochastic dynamics

$x_{t+dt}=\left\\{\begin{array}[]{ll}x_{t}+f(x_{t})dt+v_{0}\kappa(t)dt&\textnormal{with
probability }1-rdt\\\ x_{r}&\textnormal{with probability
}rdt\end{array}\right.$ (28)

where $\kappa(t)$ is a telegraph process that switches between values $+1$ and
$-1$ with a constant rate $\alpha$. Hence, $\kappa$ obeys

$\kappa(t+dt)=\left\\{\begin{array}[]{ll}\kappa(t)&\textnormal{with
probability }1-\alpha dt\\\ -\kappa(t)&\textnormal{with probability }\alpha
dt\end{array}\right.$ (29)

The value of $\kappa(t)$ indicates the particles direction of motion in which
it swims with speed $v_{0}$. In the above, $f(x_{t})$ a time-independent
external force which we will assume to be conservative
$f(x)=-\partial_{x}V(x)$ for some potential $V(x)$. The resetting position
$x_{r}$ is as before drawn at each reset from the resetting distribution
$p_{R}(x)$.

Upon resetting also the internal velocity state $\kappa(t)$ of the RTP is also
subject to resetting, and we will assume that motion in the positive
(negative) direction is chosen with probability $\rho_{+}$ (resp. $\rho_{-}$)
at each reset. Hence, in addition to the underlying telegraphic process where
$\kappa(t)$ changes sign with probability $\alpha dt$ in a timestep of size
$dt$, the process is at each resetting updated according to

$\kappa(t+dt)=\left\\{\begin{array}[]{ll}+1&\textnormal{with probability
}\rho_{+}\\\ -1&\textnormal{with probability }\rho_{-}\end{array}\right.$ (30)

We assume that the process is also initialized according to the probabilities
$\rho_{\pm}$. The state of the system is specified by the probability density
$p_{z}(x)$, where $z=\pm 1$ corresponds to the two swimming directions. The
master equation for this process can be conveniently expressed as

$\displaystyle\partial_{t}p_{+}=-v\partial_{x}p_{+}-\mu\partial_{x}[fp_{+}]-(\alpha+r)p_{+}+\alpha
p_{-}+rp_{R}\rho_{+}$ (31)
$\displaystyle\partial_{t}p_{-}=v\partial_{x}p_{-}-\mu\partial_{x}[fp_{-}]-(\alpha+r)p_{-}+\alpha
p_{+}+rp_{R}\rho_{-}$ (32)

In the case of a confining potential, it is known that the densities $p_{\pm}$
separately become stationary [51]. Since resetting effectively confines the
process by curbing large excursions, we will assume that
$\partial_{t}p_{\pm}=0$ in the steady-state, even if $V(x)=0$. In this model,
a non-equilibrium steady state is expected to be reached for two reasons; one
is the resetting dynamics which produces steady states even in passive
systems. In addition there is the effect of activity, which further pushes the
system away from thermodynamic equilibrium.

Figure 4: Comparisons between numerical simulations and theory for
harmonically confined run-and-tumble with resetting to a Gaussian distribution
$p_{R}(x)=\mathcal{N}(z,\sigma)$. a) Probability of finding the particle
moving to the right $(+)$ or to the left $(-)$ in the steady state, given by
Eq. (38), for $\rho_{+}=0.8$. b) Steady-state mean $\langle x\rangle_{*}$ as a
function of resetting rate for $z=0$. We see that the mean is generally not
monotonic for asymmetric velocity resetting $\rho_{\pm}\neq 0$. c) Steady-
state mean $\langle x\rangle_{*}$ as a function of resetting rate for $z=3$.
Points correspond to numerical simulations, while solid lines theory.
Parameters used are $\alpha=\mu=k=v_{0}=1$.

###  IV.2 Steady-state moments

In order to obtain moments for the RTP case, we proceed similarly to the
passive Brownian case. We define

$\displaystyle m_{n}$ $\displaystyle=\sum_{z=\pm}\int dxx^{n}p_{z}=\int
dxx^{n}(p_{+}+p_{-})$ (33) $\displaystyle d_{n}$
$\displaystyle=\sum_{z=\pm}\int dxx^{n}zp_{z}=\int dxx^{n}(p_{+}-p_{-})$ (34)

where $\\{m_{n}\\}$ are the full bare moments of the spatial density
$p=p_{+}+p_{-}$, while $\\{d_{n}\\}$ measures the difference in the $n$’th
moment between the right- and left-moving density. To emphasize the notation
of this section, Eq. (1) will now take the form

$m_{n}=\sum_{j=0}^{n}C_{nj}m^{R}_{j}$ (35)

where $m_{1}=\langle x\rangle_{*}$, $m_{2}=\langle x^{2}\rangle_{*}$ _et
cetera_. While $d_{n}$ does not have as clear of an interpretation as the bare
moments $m_{n}$, it does contain information about the steady state. For
example, $d_{0}$ contains information regarding the mean current in the
system, i.e. gives the probability of finding the particle in the right- or
left-moving states in the steady state through $\int dxp_{\pm}(x)=\frac{1\pm
d_{0}}{2}$, where we also used that by normalization $m_{0}=1$.

For a harmonic potential $V(x)=\frac{1}{2}kx^{2}$ the master equation gives
rise to the following steady-state hierarchy for the moments

$\displaystyle m_{n}$
$\displaystyle=\frac{r}{r+nk\mu}m_{n}^{R}+\frac{nv_{0}}{r+nk\mu}d_{n-1}$ (36)
$\displaystyle d_{n}$
$\displaystyle=\frac{r(\rho_{+}-\rho_{-})}{2\alpha+r+nk\mu}m_{n}^{R}+\frac{nv_{0}}{2\alpha+r+nk\mu}m_{n-1}$
(37)

Here we immediately see that $d_{0}=(r(\rho_{+}-\rho_{-}))/(2\alpha+r)$,
giving the probability of finding the particle swimming to the right (+) or
left(-) as

$\int dxp_{\pm}(x)=\frac{2\alpha+r\pm r(\rho_{+}-\rho_{-})}{4\alpha+2r}$ (38)

This is shown in Fig. (4 a). As $d_{n}$ in Eq. (34) is given entirely in terms
of the bare moments $m_{n-1}$, we can easily find the closed hierarchy

$\displaystyle m_{n}$
$\displaystyle=\frac{rm_{n}^{R}}{r+nk\mu}+\frac{nv_{0}}{r+nk\mu}\frac{r(\rho_{+}-\rho_{-})}{2\alpha+r+(n-1)k\mu}m_{n-1}^{R}$
(39)
$\displaystyle+\frac{n(n-1)v_{0}^{2}}{(r+nk\mu)(2\alpha+r+(n-1)k\mu)}m_{n-2}$

Although this recursive relation has somewhat complex coefficients, it is of
the general form $m_{n}=b_{n}+g_{n}m_{n-2}$ just as in the passive case, where
we defined coefficients

$\displaystyle b_{n}$
$\displaystyle=\frac{rm_{n}^{R}}{r+nk\mu}+\frac{nv_{0}}{r+nk\mu}\frac{r(\rho_{+}-\rho_{-})}{2\alpha+r+(n-1)k\mu}m_{n-1}^{R}$
(40) $\displaystyle g_{n}$
$\displaystyle=\frac{n(n-1)v_{0}^{2}}{(r+nk\mu)(2\alpha+r+(n-1)k\mu)}$ (41)

Here $m_{n}^{R}=\langle x^{n}\rangle_{R}$ is again the $n$’th moment of the
resetting distribution. In this notation, the solution reads

$m_{n}=\sum_{j\in 2\mathbf{N}}^{n}\left(\prod_{i\in
2\mathbf{N},i<j}g_{n-i}\right)b_{n-j}.$ (42)

Figure 5: Steady state variance of a run-and-tumble particle under resetting
to a Gaussian with mean $z=2$ and variance $\sigma^{2}$, in the presence of a
harmonic potential. Different colored curves correspond to different
asymmetries in the internal velocity states, controlled by $\rho_{+}$. For
high resetting width $\sigma$, the variance becomes monotonically growing,
with little distinction between different values of $\rho_{+}$. Parameters are
set to $k=\mu=\alpha=v_{0}=1$.

The mean of the run-and-tumble particle takes a form similar to the passive
Brownian case, with an additional contribution coming from the potential
asymmetry in how the internal velocity states $\sigma(t)$ are reset:

$\langle
x\rangle_{*}=\frac{rv_{0}(\rho_{+}-\rho_{-})}{(r+2\alpha)(r+k\mu)}+\frac{r}{r+k\mu}\langle
x\rangle_{R}.$ (43)

When the internal state is reset asymmetrically $\rho_{+}\neq\rho_{-}$, the
mean obtains a more complex behavior. Indeed, even for $\langle
x\rangle_{R}=0$, the first term alone gives rise to non-monotonic behavior as
a function of reset rate. This is due to the competition between the spatial
resetting and the asymmetric velocity resets. In this case, the mean is
maximized at $r=\sqrt{2k\mu\alpha}$, in agreement with the numerics in Fig (4
B), where again a Gaussian resetting distribution was considered.
Interestingly, this maximum of the mean is independent of the degree of
asymmetry $|\rho_{+}-\rho_{-}|$ between the velocity states. Also, note that
for symmetric velocity resets $\rho_{\pm}=1/2$ it is not possible to
distinguish an active RTP from a passive Brownian by looking solely at the
mean.

The fluctuations in the steady state are characterized by the variance, which
reads

$\displaystyle\kappa_{2}$
$\displaystyle=\frac{2v_{0}^{2}}{(r+2k\mu)(r+2\alpha+k\mu)}$
$\displaystyle+\frac{2rv_{0}(\rho_{+}-\rho_{-})\langle
x\rangle_{R}}{(r+2k\mu)(r+2\alpha+k\mu)}$
$\displaystyle-\frac{[rv_{0}(\rho_{+}-\rho_{-})+r(r+2\alpha)\langle
x\rangle_{R}]^{2}}{(r+2\alpha)^{2}(r+k\mu)^{2}}$
$\displaystyle+\frac{r}{r+2k\mu}\langle x^{2}\rangle_{R}$ (44)

Fig. (5) shows the variance as a function of the resetting rate for Gaussian
resetting distribution. As in the passive case, we see a non-monotonic
behavior at sharp resetting distributions for $z\neq 0$, while large resetting
widths $\sigma$ give rise to monotonic variance.

In all calculations for the run-and-tumble particle, the passive limit can be
taken to ensure consistency with the Brownian results. This limit can be taken
by letting $v_{0}\to\infty$, $\alpha\to\infty$, while keeping the ratio equal
to a constant proportional to the diffusion coefficient. This reproduced the
results of the previous sections.

##  V Conclusion

The steady-states under stochastic resetting to a distribution was studied
analytically and the results supported by numerical simulations. We have
proposed a method for calculating moments of any order in the steady state,
with results valid for arbitrary resetting distributions $p_{R}(x)$. We show
that for homogeneous systems with translational invariance, distributed
resetting will not introduce new phenomena, while for non-homogeneous ones the
presence of multiple scales may give rise to novel behavior. This includes a
strongly non-monotonic dependence of the cumulants on resetting rate in some
cases, an effect which can be washed away by introducing too broad resetting
distributions.

The presented results are useful when exact calculation of the steady state
probability density is not possible due to potentially complex nature of the
resetting distribution. Distributed resetting positions also closely mimics
resetting schemes used in experimental setups with optical tweezers. We
exemplified the framework by considering both passive Brownian particles and
active run-and-tumble particles confined by harmonic potentials in one
dimension. In both cases we calculate exact expressions for the moments, which
we further study for a Gaussian resetting distribution. Numerical simulations
show excellent agreement with the theoretical predictions.

In future works, it would be interesting to extend Eq. (1) and Eq. (8) to
other resetting schemes, including non-Poissonian waiting time distributions.
It would also be interesting to compare the presented predictions to
experiments on resetting colloids using optical tweezers, where resetting
distributions naturally appear.

###### Acknowledgements.

The author would like to thank Supriya Krishnamurthy, Francesco Mori,
Francesco Coghi and Prabal S. Negi for insightful discussions and
interactions. KSO acknowledges the Nordita fellowship program. Nordita is
partially supported by Nordforsk. KSO acknowledges support by the Deutsche
Forschungsgemeinschaft (DFG) within the project LO 418/29-1.

## References

  * Evans and Majumdar [2011a] M. R. Evans and S. N. Majumdar, Physical review letters 106, 160601 (2011a).
  * Evans  _et al._ [2020] M. R. Evans, S. N. Majumdar, and G. Schehr, Journal of Physics A: Mathematical and Theoretical 53, 193001 (2020).
  * Reuveni [2016] S. Reuveni, Physical review letters 116, 170601 (2016).
  * Pal and Reuveni [2017] A. Pal and S. Reuveni, Physical review letters 118, 030603 (2017).
  * Pal and Prasad [2019] A. Pal and V. Prasad, Physical Review E 99, 032123 (2019).
  * Ahmad  _et al._ [2019] S. Ahmad, I. Nayak, A. Bansal, A. Nandi, and D. Das, Physical Review E 99, 022130 (2019).
  * Tal-Friedman  _et al._ [2020] O. Tal-Friedman, A. Pal, A. Sekhon, S. Reuveni, and Y. Roichman, The journal of physical chemistry letters 11, 7350 (2020).
  * Cui and Joo [2019] T. J. Cui and C. Joo, RNA biology 16, 1093 (2019).
  * Grecian  _et al._ [2018] W. J. Grecian, J. V. Lane, T. Michelot, H. M. Wade, and K. C. Hamer, Journal of the Royal Society Interface 15, 20180084 (2018).
  * Pal  _et al._ [2020] A. Pal, Ł. Kuśmierz, and S. Reuveni, Physical Review Research 2, 043174 (2020).
  * Majumdar  _et al._ [2015] S. N. Majumdar, S. Sabhapandit, and G. Schehr, Physical Review E 91, 052131 (2015).
  * Singh  _et al._ [2020] R. Singh, R. Metzler, and T. Sandev, Journal of Physics A: Mathematical and Theoretical 53, 505003 (2020).
  * Fuchs  _et al._ [2016] J. Fuchs, S. Goldt, and U. Seifert, EPL (Europhysics Letters) 113, 60009 (2016).
  * Gupta  _et al._ [2020a] D. Gupta, C. A. Plata, and A. Pal, Physical review letters 124, 110608 (2020a).
  * Gupta and Plata [2022] D. Gupta and C. A. Plata, New Journal of Physics 24, 113034 (2022).
  * Pal and Rahav [2017] A. Pal and S. Rahav, Physical Review E 96, 062135 (2017).
  * Pal  _et al._ [2021] A. Pal, S. Reuveni, and S. Rahav, Physical Review Research 3, 013273 (2021).
  * Mori  _et al._ [2023] F. Mori, K. S. Olsen, and S. Krishnamurthy, Phys. Rev. Res. 5, 023103 (2023).
  * Sunil  _et al._ [2023] J. C. Sunil, R. A. Blythe, M. R. Evans, and S. N. Majumdar, arXiv preprint arXiv:2304.09348 (2023).
  * Goerlich  _et al._ [2023] R. Goerlich, M. Li, L. B. Pires, P.-A. Hervieux, G. Manfredi, and C. Genet, arXiv preprint arXiv:2306.09503 (2023).
  * Evans and Majumdar [2011b] M. R. Evans and S. N. Majumdar, J. Phys. A Math. Theor. 44, 435001 (2011b).
  * Sherman [1958] B. Sherman, The Annals of Mathematical Statistics 29, 267 (1958).
  * Feller [1954] W. Feller, Transactions of the American Mathematical Society 77, 1 (1954).
  * Pal [2015] A. Pal, Physical Review E 91, 012113 (2015).
  * Ray and Reuveni [2020] S. Ray and S. Reuveni, The Journal of chemical physics 152 (2020).
  * Bressloff [2020] P. C. Bressloff, Journal of Physics A: Mathematical and Theoretical 53, 275003 (2020).
  * Majumdar and Oshanin [2018] S. N. Majumdar and G. Oshanin, Journal of Physics A: Mathematical and Theoretical 51, 435001 (2018).
  * Wang  _et al._ [2021] W. Wang, A. G. Cherstvy, H. Kantz, R. Metzler, and I. M. Sokolov, Physical Review E 104, 024105 (2021).
  * Evans and Majumdar [2018] M. R. Evans and S. N. Majumdar, Journal of Physics A: Mathematical and Theoretical 51, 475003 (2018).
  * Kumar  _et al._ [2020] V. Kumar, O. Sadekar, and U. Basu, Physical Review E 102, 052129 (2020).
  * Tucci  _et al._ [2022] G. Tucci, A. Gambassi, S. N. Majumdar, and G. Schehr, Physical Review E 106, 044127 (2022).
  * Santra  _et al._ [2020] I. Santra, U. Basu, and S. Sabhapandit, Journal of Statistical Mechanics: Theory and Experiment 2020, 113206 (2020).
  * Pal  _et al._ [2016] A. Pal, A. Kundu, and M. R. Evans, Journal of Physics A: Mathematical and Theoretical 49, 225001 (2016).
  * Radice [2022] M. Radice, Journal of Physics A: Mathematical and Theoretical 55, 224002 (2022).
  * Eule and Metzger [2016] S. Eule and J. J. Metzger, New Journal of Physics 18, 033006 (2016).
  * Nagar and Gupta [2016] A. Nagar and S. Gupta, Physical Review E 93, 060102 (2016).
  * Gupta  _et al._ [2020b] D. Gupta, C. A. Plata, A. Kundu, and A. Pal, Journal of Physics A: Mathematical and Theoretical 54, 025003 (2020b).
  * Elgeti and Gompper [2015] J. Elgeti and G. Gompper, Europhysics Letters 109, 58003 (2015).
  * Gompper  _et al._ [2020] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, _et al._ , Journal of Physics: Condensed Matter 32, 193001 (2020).
  * Moen  _et al._ [2022] E. Q. Z. Moen, K. S. Olsen, J. Rønning, and L. Angheluta, Physical Review Research 4, 043012 (2022).
  * Bressloff [2023] P. C. Bressloff, arXiv preprint arXiv:2302.09105 (2023).
  * Debnath  _et al._ [2021] T. Debnath, P. Chaudhury, T. Mukherjee, D. Mondal, and P. K. Ghosh, The Journal of Chemical Physics 155, 194102 (2021).
  * Souzy  _et al._ [2022] M. Souzy, A. Allard, J.-F. Louf, M. Contino, I. Tuval, and M. Polin, Physical Review Research 4, L022029 (2022).
  * Volpe  _et al._ [2014] G. Volpe, S. Gigan, and G. Volpe, American Journal of Physics 82, 659 (2014).
  * Yang  _et al._ [2014] X. Yang, M. L. Manning, and M. C. Marchetti, Soft matter 10, 6477 (2014).
  * Caprini and Marconi [2018] L. Caprini and U. M. B. Marconi, Soft matter 14, 9044 (2018).
  * Caprini and Marconi [2019] L. Caprini and U. M. B. Marconi, Soft matter 15, 2627 (2019).
  * Olsen  _et al._ [2022] K. S. Olsen, L. Angheluta, and E. G. Flekkøy, Physical Review Research 4, 043017 (2022).
  * Kumar  _et al._ [2019] N. Kumar, R. K. Gupta, H. Soni, S. Ramaswamy, and A. Sood, Physical Review E 99, 032605 (2019).
  * Olsen  _et al._ [2020] K. S. Olsen, L. Angheluta, and E. G. Flekkøy, Physical Review Research 2, 043314 (2020).
  * Dhar  _et al._ [2019] A. Dhar, A. Kundu, S. N. Majumdar, S. Sabhapandit, and G. Schehr, Physical Review E 99, 032132 (2019).