# Probing the superconducting ground state of ZrIrSi: A $\mu$SR study

K. Panda  Department of Physics, Ramakrishna Mission Vivekananda Educational
and Research Institute, Howrah 711202, India  A. Bhattacharyya
<EMAIL_ADDRESS>Department of Physics, Ramakrishna Mission
Vivekananda Educational and Research Institute, Howrah 711202, India  D. T.
Adroja  ISIS Facility, Rutherford Appleton Laboratory, Chilton, Didcot, Oxon,
OX11 0QX, United Kingdom  Highly Correlated Matter Research Group, Physics
Department, University of Johannesburg, Auckland Park 2006, South Africa  N.
Kase  Department of Applied Physics, Tokyo University of Science, Tokyo,
125-8585, Japan  P. K. Biswas  ISIS Facility, Rutherford Appleton Laboratory,
Chilton, Didcot, Oxon, OX11 0QX, United Kingdom  Surabhi Saha  Department of
Physics, Indian Institute of Science, Bangalore, 560012, India  Tanmoy Das
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Science,
Bangalore, 560012, India  M. Lees  Department of Physics, University of
Warwick, Coventry CV4 7AL, United Kingdom  A. D. Hillier  ISIS Facility,
Rutherford Appleton Laboratory, Chilton, Didcot, Oxon, OX11 0QX, United
Kingdom

###### Abstract

The superconducting ground state of newly reported ZrIrSi is probed by means
of $\mu$SR technique along with resistivity measurement. The occurrence of
superconductivity at $T_{\mathrm{C}}$ = 1.7 K is confirmed by resistivity
measurement. ZF-$\mu$SR study revealed that below $T_{\mathrm{C}}$, there is
no spontaneous magnetic field in the superconducting state, indicates TRS is
preserved in case of ZrIrSi. From TF-$\mu$SR measurement, we have estimated
the superfluid density as a function of temperature, which is described by an
isotropic $s-$wave model with a superconducting gap
$2\Delta(0)/k_{\mathrm{B}}T_{\mathrm{C}}$ = 5.1, indicates the presence of
strong spin-orbit coupling. Ab-initio electronic structure calculation
indicates that there are four bands passing through the Fermi level, forming
four Fermi surface pockets. We find that the low-energy bands are dominated by
the $4d$-orbitals of transition metal Zr, with substantially lesser weight
from the $5d$-orbitals of the Ir-atoms.

###### pacs:

71.20.Be, 74.70.Dd, 76.75.+i

##  I Introduction

Recently, a number of studies have been carried out in ternary metal
phosphide, silicide and arsenide with general formula $TrT^{\prime}X$ ($Tr$
and $T^{\prime}$ are either 4$d$ or 3$d$ transition elements and whereas $X$
is either a group IV or V member) Barz1980 ; Muller1983 ; Seo1997 ; Ching2003
; Keiber1984 ; Shirotani2001 ; Shirotani1998 ; Shirotani2000 ; Shirotani1999 ;
Shirotani1995 . These systems have attracted considerable attention due to
their relatively high superconducting transition temperature
($T_{\mathrm{C}}$), for example 15.5 K for hexagonal $h-$MoNiP Shirotani2000 ,
13 K for $h-$ZrRuP Shirotani1993  and 12 K for $h-$ZrRuAs Meisner1983 . These
ternary equiatomic systems have provided a playground to investigate the role
of spin orbit (SO) coupling in superconductivity, which were not so well
studied in these systems. Compounds with Ir are often characterized with a
strong SO coupling effect, due to presence of the Ir $5d$-orbitals.
Superconductivity is observed in a number of Ir-based compounds such as
Li2IrSi3 ($T_{\mathrm{C}}$ = 4.2 K) Pyon2014 ; Lu2015 , IrGe ($T_{\mathrm{C}}$
= 4.7 K) Matthias1963 , RIr3 [$T_{\mathrm{C}}$ = 3.1 K (La), $T_{\mathrm{C}}$
= 3.3 K (Ce)] Hal ; Sato2018 ; BhattacharyyaLaIr3 , CaIr2($T_{\mathrm{C}}$ =
5.8 K) Haldolaarachchige2015 , HfIrSi ($T_{\mathrm{C}}$ = 3.1 K) Kase2016  and
ScIrP ($T_{\mathrm{C}}$ = 3.4 K) Pfannenschmidt . Cuamba et.al. Pyon2014 ;
Lu2015 ; Matthias1963 ; Cuamba  suggests that the presence strong spin orbit
(SO) coupling and a significant contribution to the total density of states
(DOS) comes from the Ir-atom in most of the Ir-based compounds. Recently, we
have reported time reversal symmetry (TRS) breaking superconductivity on the
transition metal based caged type R5Rh6Sn18 (R = Lu, Sc, and Y) Bhattacharyya1
; Bhattacharyya2 ; Bhattacharyya3  compounds due to strong spin orbit
coupling.

In ternary equiatomic compounds, superconductivity has only been found in two
types of crystal structures: the first one is the hexagonal Fe2P-type (space
group $P\bar{\mathrm{6}}m$2), Shirotani2001 ; Shirotani2000 ; Ching2003  and
the second one is the orthorhombic Co2Si-type (space group $Pnma$)
Shirotani1998 ; Shirotani2000 ; Ching2003 ; Muller1983 . It is interesting to
note that in these systems $T_{\mathrm{C}}$ is strongly associated with
crystal structure. Furthermore, $h-$Fe2P-type structure exhibit higher
$T_{\mathrm{C}}$ than $o-$Co2Si-type structure, for example: $h-$ZrRuP shows
$T_{\mathrm{C}}$ at 13.0 K whereas $o-$ZrRuP exhibit $T_{\mathrm{C}}$ at 3.5
K. In case of $h-$Fe2P-type structure each layer is fiiled up with either $Tr$
and $T^{\prime}$ or $X$ and $X$ elements. In case of $o-$ZrRuP, Shirotani et
al. Shirotani1999  reported the formation of two dimensional triangular Ru3
clusters and in the basal plane they are connected through Ru-P ionic bonds.
It also connected through Zr-Ru bonds where Zr atoms occupy the $z=1/2$ plane.
If phosphorus is replaced by the more electronegative silicon then nearest
neighbour Ru-Ru bond length enlarged to 2.87 Å. Surprisingly $o-$MoRuP shows
superconductivity at 15.5 K which is as high as isoelectronic $h-$ZrRuP
($T_{\mathrm{C}}$ = 13 K) and $h-$MoNiP ($T_{\mathrm{C}}$ = 13 K). Ching et.
al. Ching2003  shown that in $o-$MoRuP and $o-$ZrRuP higher value of DOS at
the FL is directly related to higher $T_{\mathrm{C}}$, as suggested for BCS
superconductor. In these systems, the density of states are governed by
Mo-$4d$ orbitals. Ching et. al. Ching2003  have calculated the values of
density of states which are 0.46 states per eV atom and 0.33 states per eV
atom for $o-$MoRuP and $o-$ZrRuP, respectively.

Figure 1: (a) The crystal structure of ZrIrSi which crystallize in the
orthorhombic structure with space group $Pnma$, where green, red and blue dark
yellow symbol represents Zr (big in size), Ir (medium), and Si atoms (small)
respectively. (b) Temperature variation of resistivity in zero field. Time
dependence TF-$\mu$SR asymmetry spectra for ZrIrSi recorded at (c) $T$ = 0.1 K
(d) $T$ = 1.9 K in presence of an applied field $H$ = 10 mT and at (e) $T$ =
0.1 K at (f) $T$ = 2.0 K in an applied field $H$ = 30 mT. The red solid line
shows the fittings to the data using Eqn. 1 described in the text.

Therefore, to investigate the superconducting pairing mechanism in ZrIrSi, we
have reported a systematic $\mu$SR studies. Zero field (ZF)-$\mu$SR is a
powerful technique to know the time reversal symmetry (TRS) breaking in the
superconducting state Sonier . ZF-$\mu$SR data indicates the absence of any
spontaneous magnetic fields below $T_{\mathrm{C}}$, thus implying that TRS is
not broken in the superconducting state. The superfluid density as a function
of temperature is determined from the depolarization rate of the transverse
field (TF)-$\mu$SR which is well described by an isotropic $s-$wave model.
These results are further supported by ab-initio electronic structure
calculation.

##  II Experimental Details

For this study, a polycrystalline sample of ZrIrSi was synthesized using
typical arc melting process on a water-cooled copper hearth using Zr (99.99%),
Ir (99.99%), and Si (99.999%) in a stoichiometric ratio. The arc melted ingot
was remelted several times to confirm the phase homogeneity. After that the
sample was annealed at 10000C for a week in a sealed vacuum quartz tube. X-ray
diffraction were carried out using Cu-$K_{\mathrm{\alpha}}$ radiation.
Electrical resistivity measurement was done in a standard dc-four probe
technique down to 0.5 K.

$\mu$SR experiments were performed at the ISIS pulsed neutron and muon source
of Rutherford Appleton Laboratory, UK using MuSR spectrometer with 64
detectors at transverse and longitudinal direction Lee1999 . 100% spin-
polarized muons were probed into the sample as a result positive muons
decomposed into a positron, preferably in the direction of muon spin vector
and two neutrinos with an average lifetime of 2.2 $\mu$s. These positrons are
detected by the detectors, placed either forward ($F$), or backward ($B$)
direction. The time dependence of $\mu$SR asymmetry spectra $A$ is calculated
as $A(t)=\frac{N_{\mathrm{F}}(t)-\alpha
N_{\mathrm{B}}(t)}{N_{\mathrm{F}}(t)+\alpha N_{\mathrm{B}}(t)}$, where
$N_{\mathrm{F}}(t)$ and $N_{\mathrm{B}}(t)$ are the number of positrons
counted in the forward and backward detectors respectively, and $\alpha$ is an
instrumental calibration factor. ZF-$\mu$SR is carried out with detectors in
the longitudinal configuration, where a correction coil is applied to
neutralize any stray magnetic fields up to 10-4 mT. ZF-$\mu$SR measurement is
crucial to understand the type of pairing symmetry in superconductors Sonier .
TF-$\mu$SR measurements were carried out in the vortex state in presence of
10, 20, 30 and 40 mT applied field, which is above the lower critical field
$H_{\mathrm{c1}}$ (= 0.7 mT), and below the upper critical field
$H_{\mathrm{c2}}$(= 0.6 T). The sample was mounted onto a high purity
(99.995%) silver sample holder using diluted GE-varnish and then wrapped with
thin silver foil. This was inserted in the sample chamber using a dilution
refrigerator that can go down to 50 mK. We have analyzed the $\mu$SR data
using WiMDA Pratt2000  software.

Figure 2: (a) The superconducting depolarization rate
$\sigma_{\mathrm{sc}}(T)$ as a function of temperature in presence of applied
field of 10 mT$\leq$ H$\leq$ 40 mT. (b) Magnetic field dependence of the muon
spin-depolarization rate is shown for a range of different temperatures. The
solid lines are the results of fitting the data using Brandt equation as
discussed in the text Eq. 2. (c) The inverse magnetic penetration depth
squared as a function of temperature is shown here. The lines show the fits
using $s-$wave (solid), and $d-$wave (dashed) gap functions.

##  III Results and discussion

###  III.1 Crystal structure and resistivity

X-ray powder diffraction data revealed that ZrIrSi crystallize in the
orthorhombic structure (space group $Pnma$) as displayed in Fig. 1(a). The
calculated lattice parameters are $a$ = 6.557 Å, $b$ = 3.942 Å and $c$ = 7.413
Å, which are in agreement with previous report Kase2016 . The temperature
($T$) variation of the electrical resistivity $\rho(T)$ in zero applied
magnetic field is presented in Fig. 1(b). The electrical resistivity data
reveals superconductivity at $T_{\mathrm{C}}$ = 1.7 K. Kase et. al. Kase2016
have estimated the Ginzburg Landau coherence length ($\zeta$) = 23.1 nm. It is
interesting to note that $T-$ dependence of upper critical field shows a
convex curvature Kase2016 , which might suggest the presence of SO coupling.
Similar curvature is also found in R3T4Sn13 (R = La, Sr; T = Rh, Ir) which is
a SO coupled superconductor Kase2011 .

###  III.2 TF-$\mu$SR analysis

To explore the pairing mechanism and gap structure of superconducting state of
ZrIrSi, TF-$\mu$SR measurements were performed down to 0.05 K. Fig. 1(c)-(f)
represents the TF-$\mu$SR asymmetry time spectra in presence of 10 mT and 30
mT applied magnetic field at above and below $T_{\mathrm{C}}$. Below
$T_{\mathrm{C}}$ the spectra depolarizes strongly because of the inhomogeneous
field distribution in the vortex state. TF-$\mu$SR data fitted using two
Gaussian oscillatory, functions Bhattacharyyarev ; AdrojaThFeAsN ;
BhattacharyyaThCoC2 :

$G_{\mathrm{TF}}(t)=\sum_{i=1}^{2}A_{\mathrm{i}}\cos(\omega_{\mathrm{i}}t+\phi)\exp(-\frac{\sigma_{\mathrm{i}}^{2}t^{2}}{2})$
(1)

where $A_{\mathrm{i}}$, $\sigma_{\mathrm{i}}$, $\omega_{\mathrm{i}}$, $\phi$
is the initial asymmetry, Gaussian relaxation rate, muon spin precession
frequency and the initial phase of the offset, respectively. In this fit
$\sigma_{\mathrm{i}}$ for the $2^{nd}$ part is equal to zero, which
corresponds to the background term. This term comes from those muons which
missed the sample and directly hit the silver sample holder and therefore the
depolarization of this oscillating term is zero, i.e. $\sigma_{\mathrm{2}}$ =
0 as silver has a minimal nuclear moment. $\sigma_{\mathrm{1}}$ can be
expressed as
$\sigma_{\mathrm{1}}=\sqrt{\sigma_{\mathrm{sc}}^{2}+\sigma_{\mathrm{n}}^{2}}$,
where $\sigma_{\mathrm{sc}}$ comes from superconducting part and
$\sigma_{\mathrm{n}}$ comes from nuclear magnetic dipolar moment which is
fixed in the entire temperature range, supported by the ZF-$\mu$SR later.

The temperature variation of $\sigma_{\mathrm{sc}}$ is depicted in Fig. 2(a).
As $H_{\mathrm{c2}}$ value is low in this sample, $\sigma_{\mathrm{sc}}$
depends on the applied field as displayed in Fig. 2(b). Brandt Brandt2003 ;
Brandt1988  has reported that for a superconductor with
$H_{\mathrm{ext}}/H_{\mathrm{c2}}$ $\leq$0.25, $\sigma_{\mathrm{sc}}$ is
associated with the London penetration depth [$\lambda(T)$] by the following
equation:

$\begin{split}\sigma_{\mathrm{sc}}[\mu s^{-1}]=4.83\times
10^{4}(1-H_{\mathrm{ext}}/H_{\mathrm{c2}})\\\
\times[1+1.21(1-\sqrt{{(H_{\mathrm{ext}}/H_{\mathrm{c2}}})})^{3}]\lambda^{-2}[nm]\end{split}$
(2)

This equation is a good approximation for $\kappa\geq 5$, which is valid for
our case as $\kappa=40.5$ for ZrIrSi Kase2016 . From this relation we have
determined the temperature dependence of $\lambda(T)$ and
$\mu_{0}H_{\mathrm{c2}}(T)$. Isothermal cuts perpendicular to the temperature
axis of $\sigma_{\mathrm{sc}}$ data sets were used to determine the
$H$-dependence of the depolarization rate $\sigma_{\mathrm{sc}}(H)$ as
displayed in Fig.2(b). We have estimated the London penetration depths
$\lambda$ = 254.4(3) nm, using $s$-wave model.

We have plotted the temperature variation of normalized
$\lambda^{-2}(T)/\lambda^{-2}(0)$, which is directly proportional to the
superfluid density. $\lambda^{-2}(T)/\lambda^{-2}(0)$ data were fitted using
the following equation Prozorov ; BhattacharyyaBiS2 ; AdrojaK2Cr3As3 ;
AdrojaCs2Cr3As3 ; DasLaPt2Si2 :

$\displaystyle\frac{\sigma_{sc}(T)}{\sigma_{sc}(0)}$ $\displaystyle=$
$\displaystyle\frac{\lambda^{-2}(T)}{\lambda^{-2}(0)}$ $\displaystyle=$
$\displaystyle
1+\frac{1}{\pi}\int_{0}^{2\pi}\int_{\Delta(T)}^{\infty}(\frac{\delta f}{\delta
E})\times\frac{EdEd\phi}{\sqrt{E^{2}-\Delta(T})^{2}}$

here $f$ is the Fermi function which can be expressed as
$f=[1+\exp(-E/k_{\mathrm{B}}T)]^{-1}$.
$\Delta(T,0)=\Delta_{0}\delta(T/T_{\mathrm{C}})g(\phi)$ whereas $g(\phi)$ is
the angular dependence of the gap function, $\phi$ is the azimuthal angel in
the direction of FS. The temperature variation of the superconducting gap is
approximated by the relation
$\delta(T/T_{\mathrm{C}})=\tanh\\{{1.82[1.018(T_{\mathrm{C}}/T-1)]^{0.51}}\\}$.
The spatial dependence $g(\phi$) is substituted by (a) 1 for $s-$wave gap, (b)
$|\cos(2\phi)|$ for $d-$wave gap with line nodes.

Fig. 2(c) represents the fits to the $\lambda^{-2}(T)/\lambda^{-2}(0)$ data of
ZrIrSi using a single gap $s-$wave and nodal $d-$wave models. It is clear that
the data can be well described by the isotropic $s-$wave model with a gap
value 0.374 meV. This model gives a gap to $T_{\mathrm{C}}$ ratio,
$2\Delta(0)/k_{\mathrm{B}}T_{\mathrm{C}}$ = 5.10. The higher value of gap
compare to BCS gap (3.53) suggest the presence of strong SO coupling. Similar
high gap value was obtained for Ir-based superconductors, for example: IrGe
[$2\Delta(0)/k_{\mathrm{B}}T_{\mathrm{C}}$ = 5.14] Cuamba ; Matthias1963 ,
CaIrSi3 [$2\Delta(0)/k_{\mathrm{B}}T_{\mathrm{C}}$ = 5.4] SinghCaTSi3 . On the
other hand $d-$wave model is clearly not suitable for this system as the
$\chi^{2}$ value increased significantly for this fit ($\chi^{2}$ = 6.82). As
ZrIrSi is a type II superconductor, supposing that approximately all the
normal states carriers ($n_{\mathrm{e}}$) contribute to the superconductivity
($n_{\mathrm{s}}\approx$ ne), superconducting carrier density
$n_{\mathrm{s}}$, the effective-mass enhancement $m^{*}$ have been estimated
to be $n_{\mathrm{s}}$ = 6.9(1) $\times$ 1026 carriers $m^{-3}$, and $m^{*}$ =
1.474(3) $m_{\mathrm{e}}$ respectively for ZrIrSi. Detail calculations can be
found in Ref.[42-44] Hillier1997 ; Adroja2005 ; AnandLaIrSi3 .

Figure 3: ZF-$\mu$SR asymmetry time spectra for ZrIrSi at 0.05 K (black
circle) and 2 K (dark yellow square) are shown together. The lines are the
least squares fit to the data using Eq. 4.

###  III.3 ZF-$\mu$SR analysis

In order to investigate the pairing mechanism in the superconducting ground
state, we used ZF-$\mu$SR study. The time evolution of asymmetry spectra is
shown in Fig. 3 for $T$ = 0.05 K $<T_{\mathrm{C}}$ and $T$ = 2 K
$>T_{\mathrm{C}}$. The spectra below and above $T_{\mathrm{C}}$ are found
identical, ruling out the presence of any magnetic field. This reveals that
the TRS is preserved in the superconducting state of ZrIrSi. This ZF data were
fitted by a Lorentzian function:

$G_{\mathrm{ZF}}(t)=A_{\mathrm{0}}(t)\exp{(-\lambda t)}+A_{\mathrm{bg}}$ (4)

$A_{\mathrm{0}}$ and $A_{\mathrm{bg}}$ are the sample and background asymmetry
respectively, which are nearly temperature independent. $\lambda$ is the
relaxation rate which comes from nuclear moments. The red and blue lines in
Fig. 3 indicate the fits to the ZF-$\mu$SR data. The fitting parameters of the
ZF-$\mu$SR asymmetry data are as follows: $\lambda$ = 0.030(9)
$\mu\mathrm{s}^{-1}$ at 0.05 K and $\lambda$ = 0.026(3) $\mu\mathrm{s}^{-1}$
at 2 K. The change is the relaxation rate is within the error bar, indicates
no clear evidence of TRS breaking in ZrIrSi.

###  III.4 Theoretical Calculations

ZrIrSi unit-cell has a $\it{mmm}$ point group symmetry and it belongs to the
$\it{Pnma(62)}$ space group (orthorhombic crystal structure). We have used the
Vienna Ab-initio Simulation Package(VASP) for ab-initio electronic structure
calculation. The projector augmented wave (PAW) pseudo-potentials are used to
describe the core electrons and for the exchange-correlation functional
Perdew-Burke-Ernzerhof (PBE) form is used. We have used a local density
application (LDA) functional with a cut-off energy for the plane wave basis
set of 500 eV. The Monkhorst-Pack $k$-mesh is set to $14\times 14\times 14$ in
the Brillouin zone for the self-consistent calculation. The relaxed lattice
parameters were found to be as follows $a$ = 3.9643 Å, $b$ = 6.5893 Å, and $c$
= 7.4070 Å and $\alpha=\beta=\gamma=90\degree$. To deal with the strong
correlation effect of the $d$-electrons of the Ir atoms, we employed the LDA+U
method with $U$ = 2.8 eV. For the Fermi surface calculation, we have used a
larger $k$-mesh of $31\times 31\times 31$.

In Fig. 4, we show the band structure and Fermi surface plots. We find that
there are four bands passing through the Fermi level, forming four Fermi
surface pockets. Two Fermi pockets are centered around the $\Gamma$-point,
while the other two pockets are centered around the X-point. Unlike the multi-
gap superconductivity in MgB2 MgB2  and Mo8Ga41Mo8Ga41  which are driven by
the presence of multiple Fermi surfaces, we do not find any evidence of
multigap superconductivity in ZrIrSi, which is in agreement with the
TF-$\mu$SR results. This is presumably because of the absence of
$E_{2g}$-phonon mode which could enable inter-band scattering, while in the
present case phonon modes cause intra-band electron-phonon coupling. Moreover,
we observe substantial three-dimensionality in all four Fermi surfaces. This
substantially weakens Fermi surface nesting strength. Thus the possibility of
an inter-band nesting driven unconventional $s^{\pm}$-pairing symmetry is
suppressed as compared to a two-dimensional iron-pnictide family with similar
Fermi surface topology (see, e.g., Ref. FeAs ).

Figure 4: (a) Computed DFT band structure is plotted along the standard high-
symmetric directions for the orthorhombic crystal structure of ZrIrSi. The
bands are colored with the difference in the $d$-orbital weight of the Zr and
Ir-atoms, where red color gives stronger Zr-$d$ orbital weight while blue
color dictates Ir-$d$ orbital weight. Four bands which pass through the Fermi
level are indicated by 1,2,3,4 numbers. (b) Corresponding Fermi surfaces are
plotted in the full three-dimensional Brillouin zone. We find that Fermi
surface 1 and 2 form pockets around the X-point, while Fermi surface 3 and 4
constitute pockets centering the $\Gamma$-point. The colors on the Fermi
surface only aid visualization.

Finally, we study the orbitals’ contributions to the low-energy electronic
states. We find that the low-energy bands are dominated by the $4d$-orbitals
of transition metal Zr, with substantially lesser weight from the
$5d$-orbitals of the Ir-atoms. Due to the presence of $d$-orbitals in the low-
energy states, it is natural to anticipate the involvement of strong
correlation in these systems and hence strong coupling superconductivity which
changes from typical electron-phonon to quasiparticle-phonon mechanism within
the Eliashberg theoryQPCoupling . However, to our surprise, we find a
substantially lower effective mass $1.5m_{e}$ (where $m_{e}$ is the bare
electron’s mass) which is in remarkable agreement with the results from the
TF-$\mu$SR measurements. This is well captured within our LDA+U calculation
without essentially including dynamical correlations. We repeated the
calculations for isostructural compounds TiIrSi and HfIrSi and find that the
essential Fermi surface topology and three-dimensionality remain the same in
all three materials (not shown). Therefore, we conclude that the
superconductivity in ZrIrSi and its isostructural materials (such as TiIrSi
and HfIrSi) can be well understood within the conventional BCS theory.
Although our estimates of the BCS ratio of 5.1 is slightly higher than the BCS
estimate of 3.5, however, we believe this slight increment is caused by the
spin-orbit coupling of the Ir-atoms, and the Fermi surface anisotropy.

##

##  IV Conclusions

In conclusion, we have performed ZF and TF-$\mu$SR measuremnts in the mixed
state of ZrIrSi. Using Brandt’s equation we have determine the temperature
dependence of the magnetic penetration depth. The superfluid density
$n_{\mathrm{s}}\propto 1/\lambda^{2}$ well described by an isotropic $s-$wave
model. The obtained gap value is $2\Delta(0)/k_{\mathrm{B}}T_{\mathrm{C}}$ =
5.1, which suggest ZrIrSi to be a strongly coupled BCS superconductor. Ab-
initio electronic structure calculation indicates BCS superconductivity, which
supports our experimental results. The low-energy bands are dominated by the
$4d-$orbitals of the transition metal Zr, with a substantially lesser weight
from the $5d-$orbitals of the Ir-atoms. ZF-$\mu$SR reveals, there is no
spontaneous magnetic field below $T_{\mathrm{C}}$, which suggest the absence
of TRS breaking. The present results pave the way to develop a realistic
theoretical model to interpret the origin of superconductivity in ternary
systems.

### Acknowledgements

KP acknowledge the financial support from DST India, for Inspire Fellowship
(IF170620). AB would like to acknowledge DST India, for Inspire Faculty
Research Grant (DST/INSPIRE/04/2015/000169), and UK-India Newton funding for
funding support. DTA and ADH would like to thank the Royal Society of London
for the UK-China Newton funding and CMPC-STFC, grant number CMPC-09108, for
financial support. TD acknowledges the financial support from the Science and
Engineering Research Board (SERB), Department of Science & Technology (DST),
Govt. of India for the Start-Up Research Grant (Young Scientist).

## References

  * (1) H. Barz, H. C. Ku, G. P. Meisner, Z. Fisk and B. T. Matthias Proc. Natl. Acad. Sci. 77, 3132 (1980).
  * (2) R. Müller, R. N. Shelton, J. W. Richardson, and R. A. Jacobson J. Less-Common Met. 42, 177 (1983).
  * (3) D. K. Seo, J. Ren, M. H. Whangbo and E. Canadell Inorg. Chem. 36, 6058-6063 (1997).
  * (4) W. Y. Ching, Y. N. Xu, L. Ouyang, and W. Wong-Ng J. Appl. Phys. 93, 8209 (2003).
  * (5) H. Keiber, H. W$\ddot{u}$hl, G. P. Meisner and G. R. Stewart J. Low Temp. Phys. 55, 11 (1984).
  * (6) I. Shirotani, M. Takaya, I. Kaneko, C. Sekine and T. Yagi Physica C 357-360, 329 (2001).
  * (7) I. Shirotani, Y. Konno, Y. Okada, C. Sekine, S. Tado and T. Yagi Solid State Commun. 108, 967 (1998).
  * (8) I. Shirotani, M. Takeya, I. Kaneko, C. Sekine, and T. Yagi Solid State Comm. 116, 683 (2000).
  * (9) I. Shirotani, K. Tachi, Y. Konnon, S. Todo, and T. Yagi Solid State Commun. 79, 767 (1999).
  * (10) I. Shirotani, K. Tachi, K. Takeda, S. Todo, T. Yagi and K. Kanoda Phys. Rev. B 52, 6197 (1995).
  * (11) I. Shirotani, N. Ichihashi, K. Nozawa, M. Kinoshita, T. Yagi, K. Suzuki, and T. Enoki Jpn. J. Appl. Phys. 32-3, 695 (1993).
  * (12) G. P. Meisner, H. C. Ku, and H. Barz Matter. Res. Bull. 18, 983 (1983).
  * (13) S. Pyon, K. Kudo, J. Matsumura, H. Ishii, G. Matsuo, M. Nohara, H. Hojo, K. Oka, M. Azuma, O. Garlea, K. Kodama, and S. Shamoto, J. Phys. Soc. Jpn. 83, 093706 (2014).
  * (14) H. Y. Lu, N. N. Wang, L. Geng, S. Chen, Y. Yang, W. J. Lu, W. S. Wang, and J. Sun, Europhys. Lett. 110, 17003 (2015).
  * (15) B. T. Matthias, T. H. Gkballk and V. B. Compton Rev. Mod. Phys. 35, 1 (1963).
  * (16) N. Haldolaarachchige1, L. Schoop, M. A. Khan, W. Huang, H. Ji, K. Hettiarachchilage and D. P. Young J. Phys.: Condens. Matter 29, 475602 (2017).
  * (17) Y. J. Sato, A. Nakamura, Y. Shimizu, A. Maurya, Y. Homma, D. Li, F. Honda, and D. Aoki J. Phys. Soc. Jpn. 87, 053704 (2018).
  * (18) A. Bhattacharyya, D. T. Adroja, P. K. Biswas, Y. J. Sato, M. R. Lees, D. Aoki, and A. D. Hillier arXiv: 1901.04796(2019).
  * (19) N. Haldolaarachchige, Q. Gibson, L. M. Schoop, H. Luo, and R. J. Cava, J. Phys. Condens. Matter 27, 185701 (2015).
  * (20) N. Kase, H. Suzuki, T. Nakano, and N. Takeda Supercond. Sci. Technol 29, 035011 (2016).
  * (21) U. Pfannenschmidt, U. C. Rodewald, and R. P. Ottgen, Z. Naturforsch. B 66, 205 (2011).
  * (22) A. S. Cuamba, H.-Y. Lu, and C. S. Ting Phys. Rev. B 94, 094513 (2016).
  * (23) A. Bhattacharyya, D. T. Adroja, J. Quintanilla, A. D. Hillier, N. Kase, A. M. Strydom and J. Akimitsu Phys. Rev. B 91, 060503 (2015).
  * (24) A. Bhattacharyya, D. T. Adroja, N. Kase, A. Hillier, J. Akimitsu and A. M. Strydom Sci. Rep. 5, 12926 (2015).
  * (25) A. Bhattacharyya, D. T. Adroja, N. Kase, A. D. Hillier, A. M. Strydom and J. Akimitsu Phys. Rev. B 98, 024511 (2018).
  * (26) J. E. Sonier, J. H. Brewer, and R. F. Kiefl Rev. Mod. Phys. 72, 769 (2000).
  * (27) S. L. Lee, S. H. Kilcoyne, and R. Cywinski, Muon Science: Muons in Physics, Chemistry, and Materials (SUSSP and IOP, Bristol, UK, 1999).
  * (28) F. L. Pratt, Physica B 289-290, 710 (2000).
  * (29) N. Kase, H. Hayamizu, and J. Akimitsu Phys. Rev. B 83, 184509 (2011).
  * (30) A. Bhattacharyya, D. T. Adroja, M. Smidman, and V. K. Anand, Sci. China-Phys. Mech. Astron. 61, 127402 (2018).
  * (31) D. T. Adroja, A. Bhattacharyya, P. K. Biswas, M. Smidmann, A. Hillier, H. Mao, H. Luo, G. H. Cao, Z. Wang and C. Wang Phys. Rev. B 96, 144502 (2017).
  * (32) A. Bhattacharyya, D. T. Adroja, K. Panda, S. Saha, T. Das, A. J. S. Machado, T. W. Grant, Z. Fisk, A. D. Hillier and P. Manfrinetti arXiv:1810.06072 (2018).
  * (33) E. H. Brandt, J. Low Temp. Phys. 73, 355 (1988).
  * (34) E. H. Brandt, Phys. Rev. B 68, 054506 (2003).
  * (35) R. Prozorov, and R. W. Giannetta, Supercond. Sci. Technol. 19, R41 (2006).
  * (36) D. Adroja, A. Bhattacharyya, M. Telling, Y. Feng, M. Smidman, B. Pan, J. Zhao, A. Hillier, F. Pratt, and A Strydom Phys. Rev. B 92, 134505 (2015).
  * (37) D. T. Adroja, A. Bhattacharyya, M. Smidman, A. D. Hillier, Y. Feng, B. Pan, J. Zhao, M. R. Lees, A. M. Strydom, and P. K. Biswas J. Phys. Soc. Jpn. 86, 03470 (2017).
  * (38) D. Das, R. Gupta, A. Bhattacharyya, P. K. Biswas, D. T. Adroja, and Z. Hossain Phys. Rev. B 97, 184509 (2018).
  * (39) A. Bhattacharyya, D. T. Adroja, A. D. Hillier, R. Jha, V. P. S. Awana, and A. M. Strydom J. Phys: Cond. Matt. 29, 265602 (2017).
  * (40) G. Pang, M. Smidman, W. Jiang, J. Bao, Z. Weng, Y. Wang, L. Jiao, J. Zhang, G. Cao, and H. Yuan Phys. Rev. B 91, 220502 (2015).
  * (41) R. P. Singh, A. D. Hillier, D. Chowdhury, J. A. T. Barker, D. McK. Paul, M. R. Lees, and G. Balakrishnan Phys. Rev. B 90, 104504 (2014).
  * (42) A. Hillier and R. Cywinski Appl. Magn. Resonan. 13 95–109 (1997).
  * (43) D. Adroja, A. Hillier, J. G. Park, E. Goremychkin, K. McEwen, N. Takeda, R. Osborn, B. Rainford, and R. Ibberson Phys. Rev. B 72, 184503 (2005).
  * (44) V. K. Anand, D. Britz, A. Bhattacharyya, D. T. Adroja, A. D. Hillier, A. M. Strydom, W. Kockelmann, B. D. Rainford and K. A. McEwen Phys. Rev. B 90 (1), 014513 (2014).
  * (45) Amy Y. Liu, I. I. Mazin, and Jens Kortus, Phys. Rev. Lett. 87, 087005 (2001).
  * (46) Anshu Sirohi, Surabhi Saha, Prakriti Neha, Shekhar Das, Satyabrata Patnaik, Tanmoy Das and Goutam Sheet Phys. Rev. B 99, 054503 (2019).
  * (47) P. J. Hirschfeld, M. M. Korshunov, and I. I. Mazin, Rep. Prog. Phys. 74, 124508 (2011); T. Das, EPJ Web of Conferences 23, 00014 (2012).
  * (48) T. Das and K. Dolui, Phys. Rev. B 91, 094510 (2015).


# Complex dynamics near extinction in a predator-prey model with ratio
dependence and Holling type III functional response

K.B. Blyuss, Y.N. Kyrychko

Department of Mathematics, University of Sussex, Falmer,
Brighton, BN1 9QH, United Kingdom

Corresponding author. Email<EMAIL_ADDRESS>O.B. Blyuss

Wolfson Institute of Population Health, Queen Mary University of London,
London, EC1M 6BQ, United Kingdom

###### Abstract

In this paper, we analyse a recently proposed predator–prey model with ratio
dependence and Holling type III functional response, with particular emphasis
on the dynamics close to extinction. By using Briot-Bouquet transformation we
transform the model into a system, where the extinction steady state is
represented by up to three distinct steady states, whose existence is
determined by the values of appropriate Lambert W functions. We investigate
how stability of extinction and coexistence steady states is affected by the
rate of predation, predator fecundity, and the parameter characterising the
strength of functional response. The results suggest that the extinction
steady state can be stable for sufficiently high predation rate and for
sufficiently small predator fecundity. Moreover, in certain parameter regimes,
a stable extinction steady state can coexist with a stable prey-only
equilibrium or with a stable coexistence equilibrium, and it is rather the
initial conditions that determine whether prey and predator populations will
be maintained at some steady level, or both of them will become extinct.
Another possibility is for coexistence steady state to be unstable, in which
case sustained periodic oscillations around it are observed. Numerical
simulations are performed to illustrate the behaviour for all dynamical
regimes, and in each case a corresponding phase plane of the transformed
system is presented to show a correspondence with stable and unstable
extinction steady state.

##  1 Introduction

A starting point for many ecological models analysing predator-prey
interactions is provided by models of Gause-Kolmogorov type [1, 2, 3]

$\begin{array}[]{l}\dot{u}=uf(u)-vg\left(u,v\right),\\\ \\\
\dot{v}=bvg\left(u,v\right)-dv,\end{array}$

where $u(t)$ and $v(t)$ are population densities or abundances of prey and
predator, $d$ is the natural death rate of predators, and functions $f(u)$ and
$g(u,v)$ describe, respectively, natural per capita growth rate of prey, and
the trophic function or functional response of predators [4], which quantifies
how efficiently predators are consuming prey, and how prey consumption
increases predator reproduction. Function $f(\cdot)$ is often chosen in the
form of either a constant (so-called Malthusian growth), or as a monotonically
decreasing linear function $f(u)=r(1-u/K)$, which describes intra-specific
competition of prey for available resources, and is known as the logistic
growth of prey with a linear growth rate $r$ and the carrying capacity $K$. In
the simplest case, the function $g(u,v)$ would only depend on the abundance of
prey $g(u,v)=g(u)$, in which case it is known as prey-dependent functional
response. Perhaps, the simplest example of prey-dependent functional response
is when it is proportional to the number of prey $g(u)=au$, which yields the
classical Lotka-Volterra model [5, 6]. A more realistic representation of
interactions between prey and predators is given by a functional response
suggested by Holling [7, 8], which accounts for the time predators spend
searching for prey, as well as handling (i.e. chasing, killing and digesting)
the prey. Holling proposed three types of functional response $g(u)$, which
all satisfy $g(0)=0$ and approach some constant value for large $u$. The
difference between different types of Holling functional response is in the
behaviour for small prey numbers/densities. Type I response is linearly
increasing for small prey densities, whereas for large prey numbers it
saturates at some constant value; type II and type III are functions that are
also saturating at high prey numbers, and are, respectively, concave and
sigmoidal. Besides direct impact of predation, prey behaviour can also be
affected by the fear of predation, which has been explored in a series of
recent papers [9, 10, 11, 12, 13].

By the early 1990s, a number of ecological studies showed that when studying
predator-prey interactions across a variety of scales and biological species,
from insects to mammals [14, 15, 16, 17], it may be more realistic to consider
the functional response that depends on the ratio of prey to predators, which
became known as ratio-dependent functional response [18, 19]. Such functional
response can effectively model an observation that at higher predator
densities, predators would have to compete for and share some of the prey.
Formally, this would correspond to writing the functional response $g(u,v)$ in
the form $g(u/v)$, which would then yield the following model

$\begin{array}[]{l}\displaystyle{\dot{u}=uf(u)-vg\left(\frac{u}{v}\right),}\\\
\\\ \displaystyle{\dot{v}=bvg\left(\frac{u}{v}\right)-dv.}\end{array}$

A number of theoretical studies have studied the dynamics of models with ratio
dependence [20, 21, 22, 23, 24, 25, 26] and various forms of functional
response that have usually taken the form of $g(z)=a$,
$\displaystyle{g(z)=a_{1}z/(1+a_{2}z)}$, and
$\displaystyle{g(z)=a_{1}z^{n}/(1+a_{2}z^{n})}$ with $n>1$, where $z=u/v$, for
Holling type I, II and III responses, respectively [27, 28]. Some additional
ecological effects that have been studied in models with ratio dependence
include time delays [29, 30, 31], as well as spatial dependence [32, 33, 34].
One important observation to make here is that due to the very nature of
ratio-dependent functional response, the origin, i.e. the point $(u,v)=(0,0)$
that is characterised by the extinction of both species, may create certain
mathematical challenges for analysis from the perspective that either the
vector field itself is not defined at that point, or it is defined, but the
linearisation is not, which hinders standard stability analysis. To overcome
these challenges, one approach that is often used consists in rescaling time
with the denominator of functional response, often in combination with Briot-
Bouquet transformation to remove the singularity [22, 24, 26].

In this paper we are interested in Holling type III functional response with
ratio dependence that has been observed in a number of experimental settings
[35, 36, 37, 38]. Rather than using a functional form
$\displaystyle{g(z)=a_{1}z^{n}/(1+a_{2}z^{n})}$ given by a ratio of two
polynomials, where $z=u/v$ and $n>1$, we consider an alternative form that can
be written as follows,

$g(z)=ae^{-\alpha/z},$ (1)

which satisfies the conditions of $g(0)=0$, is monotonically increasing, and
settles at a constant value as $z\to\infty$. Biological motivation for this
form of functional response comes from the studies of plant-parasite
interactions, in which plants serve as hosts for insect predators, and hence,
current plant population represents a carrying capacity for insect population,
in which case the interaction term has the form $g(u,v)=a(1-v/u)$ [39, 40].
Formally, trophic function (1) is reminiscent of the Ivlev trophic function
$g(z)=a\left(1-e^{-\alpha z}\right)$ [41] that represents Holling type II
response, as well as of the Ricker model $N_{t+1}=N_{t}\exp[r(1-N_{t})]$ for
single-species populations [42]. Making ratio dependence in (1) explicit, we
have the functional response $g(u,v)$ being given by

$g(u,v)=ae^{-\alpha v/u}.$ (2)

A recent work has used a delayed version of this trophic function to analyse
the dynamics of vector-plant interactions in the context of modelling plant
mosaic disease [43], while the role of stochastic effects in the time-delayed
model has been studied in [44]. Of particular interest to us is the analysis
of model behaviour in the neighbourhood of extinction $(u,v)=(0,0)$, and the
effects this has on overall dynamics. While earlier numerical simulations have
suggested that extinction is indeed possible in some parameter regimes [44],
the question of when exactly this happens, and whether extinction can coexist
with other states has not been explored.

The remainder of this paper is organized as follows. In the next section we
identify different steady states of the model and show that the origin that
corresponds to extinction of both species is a well-defined steady state.
Applying Briot-Bouquet transformation, we will then obtain another version of
our model, where the original extinction steady state is unfolded as an entire
axis in the phase space. Whereas the origin of the modified model is always
unstable, there are two more steady states that also correspond to the
extinction steady state of the original model, and whose locations are given
by a Lambert W function. We will derive analytical conditions for stability of
these steady states and will illustrate regions of their feasibility and
stability depending on model parameters. Model dynamics is further explored by
numerically computing regions of feasibility and stability for different
steady states and periodic orbits of the model depending on predation rate,
predator growth rate, and the parameter characterising functional response. We
also demonstrate phase plane of the modified model in each scenario, clearly
indicating different steady states that correspond to extinction and
coexistence, and complement these by numerical solutions of the original
model, which can exhibit such distinct types of behaviour, as extinction of
both species; survival of prey only, when predator growth rate is not
sufficiently high; a regime of bi-stability between these two scenarios, where
for the same values of parameters it is the initial values that determine
whether both species or only predators will go extinct. Other possibilities
include a stable coexistence steady state, where both prey and predators are
maintained at some constant level, suggesting that predation is compensated or
balanced by the prey growth, a periodic solution around the coexistence steady
state, which is reminiscent of oscillations in the standard Lotka-Volterra
model, and a regime of bi-stability between a stable coexistence and
extinction. The paper concludes with a discussion of results.

##  2 Methods

As discussed in the Introduction, we consider a predator-prey model with
logistic growth of prey and a Holling type III functional response with ratio
dependence, which has the form

$\begin{array}[]{l}\displaystyle{\dot{u}=ru\left(1-\frac{u}{K}\right)-ave^{-\alpha
v/u},}\\\ \\\ \displaystyle{\dot{v}=bve^{-\alpha v/u}-dv.}\end{array}$ (3)

In order to simplify the model and to reduce the number of free parameters, we
rescale the variables and parameters as follows,

$\begin{array}[]{l}u=K\widehat{u},\quad v=K\widehat{v},\quad
rt=\widehat{t},\\\ \\\
\displaystyle{\frac{d}{r}=\widehat{d},\quad\frac{a}{r}=\widehat{a},\quad\frac{b}{r}=\widehat{b}.}\end{array}$

The model (3) can then be rewritten in the form

$\begin{array}[]{l}\displaystyle{\dot{u}=u(1-u)-ave^{-\alpha v/u},}\\\ \\\
\displaystyle{\dot{v}=bve^{-\alpha v/u}-dv,}\end{array}$ (4)

where we have dropped hats for notational convenience. It is straightforward
to show that this system is well-posed in that for arbitrary non-negative
initial conditions, its solutions will remain non-negative and bounded for all
$t\geq 0$. For any values of parameters, the system (4) has a prey-only steady
state $E=(1,0)$, which is stable for $b<d$ and unstable for $b>d$ [44]. If
$b>d$ and $\alpha b-ad\ln(b/d)>0$, this system also has a coexistence steady
state $E^{*}=(u^{*},v^{*})$, with

$\displaystyle{u^{*}=\frac{1}{\alpha b}\left[\alpha
b-ad\ln\left(\frac{b}{d}\right)\right],\quad v^{*}=\frac{1}{\alpha
b}\ln\left(\frac{b}{d}\right)\left[\alpha
b-ad\ln\left(\frac{b}{d}\right)\right].}$

This steady state is stable, provided the following condition is satisfied
[44]

$-1+d\ln\left(\frac{b}{d}\right)\left[\frac{a}{\alpha
b}\left(2-\ln\left(\frac{b}{d}\right)\right)-1\right]<0.$

In an earlier paper [44], we focused on the analysis of coexistence steady
state $E^{*}=(u^{*},v^{*})$ under effects of maturation time delay in
predators, and also investigated the role of stochasticity.

Figure 1: Lambert functions $W_{0}(x)$ and $W_{-1}(x)$.

In the neighbourhood of the point $(u,v)=(0,0)$, in the first quadrant the
term $ve^{-\alpha v/u}$ is well-defined and positive, and in the limit
$(u,v)\to(0,0)$, we have

$\displaystyle{0\leq\left|ve^{-\alpha
v/u}\right|\leq|v|\xrightarrow[(u,v)\to(0,0)]{}0,}$

which shows that $(0,0)$ is indeed a steady state of the system (4), which
biologically describes extinction of both species. Since the system (4) is not
differentiable at the point $(0,0)$, to explore its dynamics close to this
steady state, we employ the Briot-Bouquet transformation [45] that in our case
consists in introducing a new variable $z$, such that $v=zu$, and the system
(4) then transforms into

$\begin{array}[]{l}\displaystyle{\dot{u}=u(1-u)-azue^{-\alpha z},}\\\ \\\
\displaystyle{\dot{z}=z(b+az)e^{-\alpha z}-z(d+1-u).}\end{array}$ (5)

The Briot-Bouquet transformation blows up the origin into an entire $z$-axis
[46]. Hence, we should investigate equilibria of the system (5) that lie on
the $z$-axis.

For any parameter values, the transformed system (5) has a steady state
$E_{0}=(0,0)$. Besides this steady state, as long as the following condition
is satisfied

$\displaystyle{\frac{\alpha(d+1)e^{-b\alpha/a}}{a}\leq\frac{1}{e},}$ (6)

the system (5) can have up to two more steady states $E_{1}=(0,z_{1})$ and
$E_{2}=(0,z_{2})$, where

$\displaystyle{z_{1}=-\frac{b}{a}-\frac{1}{\alpha}W_{0}\left[-\frac{\alpha(d+1)e^{-b\alpha/a}}{a}\right],\quad
z_{2}=-\frac{b}{a}-\frac{1}{\alpha}W_{-1}\left[-\frac{\alpha(d+1)e^{-b\alpha/a}}{a}\right],}$

and $W_{0}(\cdot)$ and $W_{-1}(\cdot)$ are two branches of Lambert W function.
Figure 1 illustrates the dependence of these two functions on their argument,
and we note that the smallest possible value of $W_{0}(x)$ is attained at
$x=-1/e$ and is equal to $W_{0}(-1/e)=-1$, which suggests that the steady
state $E_{1}$ is only biologically feasible when $0<\alpha b/a<1$.
Furthermore, since $W_{-1}(x)\leq W_{0}(x)<0$ for $x<0$, this implies that
whenever both $E_{1}$ and $E_{2}$ exist, we have $z_{1}\leq z_{2}$.

Linearisation of the system (5) near any steady state
$\widehat{E}=(\widehat{u},\widehat{z})$ gives the Jacobian

$J=\left|\begin{array}[]{cc}1-2\widehat{u}-a\widehat{z}e^{-\alpha
z}&-a\widehat{u}e^{-\alpha\widehat{z}}(1-\alpha\widehat{z})\\\
\widehat{z}&(b+a\widehat{z})(1-\alpha\widehat{z})e^{-\alpha\widehat{z}}+aze^{-\alpha\widehat{z}}-(d+1-\widehat{u})\end{array}\right|.$
(7)

At $E_{0}$, the two eigenvalues are $1$ and $b-d+1$, suggesting that this
steady state is always unstable and is either a saddle, or an unstable node.
At the steady states $E_{1}$ and $E_{2}$, we have the Routh-Hurwitz conditions
for stability in the form

$1<a\widehat{z}e^{-\alpha\widehat{z}}<\alpha(d+1)\widehat{z}.$

Figure 2: Regions of feasibility and stability of steady states $E_{1}$ (a,c)
and $E_{2}$ (b,d). Steady states are feasible in a coloured region, and
biologically infeasible in the white region. Colour code denotes the largest
real part of characteristic eigenvalues. Parameter values are $a=d=1$ (a,b),
$\alpha=d=1$ (c,d).

##  3 Results

In order to illustrate the effects of different parameters on feasibility and
stability of the two steady states $E_{1}$ and $E_{2}$ that correspond to
extinction in the original model, we plot in Fig. 2 maximum real part of the
characteristic eigenvalues as determined by the Jacobian (7). We observe that
the steady state $E_{1}$ only exists for sufficiently small levels of predator
growth rate $b$, and either low enough predation rate $a$, or sufficiently
small rate $\alpha$. For any parameter combination, where this steady state is
biologically feasible, it is unstable. In contrast, starting with rather small
values of $b$, the steady state $E_{2}$ is biologically feasible for an entire
range of values of $a$ or $\alpha$ and is stable for smaller values of $a$, or
for larger values of $\alpha$, and the region of stability shrinks with
increasing predator growth rate $b$.

Figure 3: Regions of feasibility and stability of different steady states of
the model (4). The prey-only steady state $E=(1,0)$ is stable to the left of
the dashed black line at $b=1$, and unstable to the right of that line. The
extinction steady state $E_{0}$ is stable above red line in plot (a) and below
red line in plot (b). In region 1, there is bi-stability between a stable
steady state $E$ and the extinction steady state $E_{0}$; in region 2 only the
steady state $E$ is stable; in region 3, only the coexistence steady state
$E^{*}$ is stable; in region 4, the coexistence steady state $E^{*}$ is
unstable, and there is periodic orbit around this state; in region 5, there is
bi-stability between the extinction steady state $E_{0}$ and the coexistence
steady state $E^{*}$; in region 6, only extinction steady state $E_{0}$ is
stable. Parameter values are $\alpha=d=1$ (a), $a=d=1$ (b).

Under the inverse Brio-Bouquet transformation, we can interpret regions of
stability of the steady state $E_{2}$ of the transformed system (5) as
parameter regions, where the extinction steady state $E_{0}=(0,0)$ of the
original system (4) is stable. Combining these results with those on stability
of the prey-only equilibrium $E$ and the coexistence equilibrium $E^{*}$, we
now illustrate in Fig. 3 an entire range of dynamical scenarios that can be
exhibited in the model.

For values of $b$ smaller than $1$, which biologically describes a scenario,
where the rate of growth in the population of predators is too small to be
able to maintain their population as based on the available prey, the
coexistence steady state $E^{*}$ is not biologically feasible, while the prey-
only steady state $E$ is stable. In this case, the system either goes to this
prey-only steady state as the only stable equilibrium for parameter values in
region 1, or it exhibits a bi-stability between this steady state, and a
stable extinction steady state $E_{0}$ in region 2. This is shown in detail in
Fig. 4, where, for convenience, for each parameter combination, in the right
column we plot numerical solution of the model (4), and in the left column we
plot the corresponding phase plane of the transformed system (5) to
demonstrate how changes in stability of the extinction steady state $E_{0}$
affect global dynamics. Since in region 1 there is bi-stability between two
steady states, both of which are characterised by the absence of predators,
this suggests that it is rather the initial conditions of the system that
determine whether prey will survive (though, initially, prey population will
decline due to predation), or will also be driven to extinction.
Interestingly, even though in this parameter regime, predator population goes
to zero, this does not automatically guarantee the survival of the prey, which
is affected not only by the values of parameters, but also by the initial
presence of sufficiently many prey to maintain its stable population.

Figure 4: Phase planes (left column) of the transformed system (5) and the
corresponding numerical solution (right column) of the rescaled original
system (4), illustrating different dynamical regimes listed in Fig. 3(a).
Blue/red dots in phase plane plots indicate stable/unstable steady states.
Parameter values are $\alpha=d=1$. (a)-(b) Case 1, $a=6$, $b=0.5$, bi-
stability between prey-only steady state $E$ and the extinction steady state
$E_{0}$. (c)-(d) Case 2, $a=3$, $b=0.5$, prey-only steady state $E$ is the
only stable state. (e)-(f) Case 3, $a=1.5$, $b=1.4$, only the coexistence
steady state $E^{*}$ is stable.

As the value of $b$ exceeds 1, the prey-only equilibrium $E$ becomes unstable,
while the coexistence steady state $E^{*}$ becomes feasible and stable, as
shown in Fig. 4(e)-(f) (region 3). As the value of predation rate $a$
increases, as long it stays below the boundary of stability of the extinction
steady state $E_{0}$, the coexistence steady state $E^{*}$ loses its stability
through a supercritical Hopf bifurcation, giving rise to stable periodic
solutions around this steady state, as illustrated in Fig. 5(a)-(b) (region
4). Once the stability boundary of the extinction steady state is crossed, for
smaller values of $b$ and sufficiently high values of predation rate $a$,
stable coexistence steady state $E^{*}$ co-exists with a stable extinction
steady state $E_{0}$ shown in Fig. 5(c)-(d) (region 5). Biologically, this
means that for the same parameter values, whether or not the system evolves
towards coexistence of prey and predators, or they drive each other to
extinction, is again determined by the choice of initial conditions. In region
6, which corresponds to high values of both predation rate and the growth rate
of predators, the only stable state of the system is that of extinction, as
illustrated in Fig. 5(e)-(f).

Figure 5: Phase planes (left column) of the transformed system (5) and the
corresponding numerical solution (right column) of the rescaled original
system (4), illustrating different dynamical regimes listed in Fig. 3(a).
Blue/red dots in phase plane plots indicate stable/unstable steady states.
Parameter values are $\alpha=d=1$. (a)-(b) Case 4, $a=3.3$, $b=1.5$, periodic
solution around the unstable coexistence steady state $E^{*}$. (c)-(d) Case 5,
$a=6$, $b=1.05$, bi-stability between the extinction steady state $E_{0}$ and
the coexistence steady state $E^{*}$. (e)-(f) Case 6, $a=6$, $b=1.4$, only the
extinction steady state $E_{0}$ is stable.

##  4 Discussion

In this article, we have studied the dynamics of a predator-prey model with
Holling type III functional response and ratio dependence. Earlier numerical
results for the same model with maturation delay suggested the possibility of
extinction of both prey and predator populations that was, however, not
investigated at that time [44]. Whereas extinction steady state is a feasible
steady state of the model for any values of parameters, due to ratio
dependence in the functional response, the model could not be linearised in
the neighbourhood of extinction steady state in the standard manner. To
overcome this difficulty, we have transformed the system using a Briot-Bouquet
transformation, which resulted in the unfolding of the extinction steady state
into three distinct steady states, whose stability could now be studied. Two
of those steady states are unstable for any values of parameters, while the
third is stable for sufficiently large predation rate to have an effect on
prey who are reproducing logistically, and for sufficiently small values of
the parameter $\alpha$ characterising functional response. In both cases, we
note that the region of stability of extinction steady state is shrinking for
higher rates of predator fecundity, which can be explained by the fact that if
for the same predation rate predators are better able to reproduce, this
increases the possibility of a stable coexistence, where both populations are
maintained at some steady level. Regions of feasibility of two steady states
of the transformed model that correspond to extinction in the original model
are determined explicitly in terms of system parameters, since they are given
in the form of principal and $n=-1$ branches of the, in general, complex-
valued Lambert W function.

By combining the results on stability of extinction steady state with
conditions for stability of prey-only and coexistence equilibria, we are able
to obtain a full picture of system dynamics in different parameter regions,
identifying situations where only one of those steady states is stable, as
well as cases of bi-stability, where extinction coexists with either a prey-
only, or coexistence equilibria. In the parameter region, where the system
exhibits sustained periodic oscillations around an unstable coexistence
equilibrium, resulting from a supercritical Hopf bifurcation of this steady
state, the extinction steady state is unstable, hence, no bi-stability is
observed. While coexistence equilibrium is only biologically feasible, when
the prey-only steady state is unstable, in contrast, the extinction steady
state exists for all values of parameters and can coexist with either of those
steady states. In the case of coexistence of steady states, for the same
parameter values, depending on the initial conditions, the system will
approach one of the two coexisting steady states. This observation is
important from the perspective of analysing real ecological data, as it
provides clues about the possibility of different types of dynamics for the
same or very similar values of parameters, some of which may be difficult to
accurately measure. Earlier papers that focused largely on ratio dependence in
models with Holling type II functional response [22, 24, 26] primarily showed
either extinction of both species, or their stable coexistence. In contrast,
our model with Holling type III response exhibits significant parameter
regions, where prey-only equilibrium is stable either by itself, or in
combination with extinction steady state. On the other hand, in the case when
coexistence steady state is unstable, and there is a periodic solution around
it, it is the only attractor in the model.

There are several directions, in which the work presented in this paper could
be extended. In an earlier paper [44], we studied the effects of stochasticity
and maturation delays for the same model, and some of the simulations there
also indicated the possibility of observing extinction in the time-delayed
model. Unfortunately, the standard Briot-Bouquet transformation, as used in
this paper, would not work for the analysis of stability of extinction steady
state in the time-delayed model, because even after the transformation, the
singularity at the point $(0,0)$ would remain, and some alternative approach
for studying stability of extinction equilibrium in a model with time delay
needs to be developed. Another feature of many ecological models is the so-
called stochastic extinction, where for sufficiently small population
densities, species can go extinct, even though deterministically they could
recover to some sustained levels at a later stage. It would be interesting and
important to investigate how (in)stability of extinction steady state is
affected by stochasticity, and how it impacts the regimes of bi-stability
between extinction and other equilibria. This could be achieved through
numerical exploration of basins of attraction of different steady states in
the stochastic model in a manner similar to how it was done by Fatehi et al.
[47, 48] in the context of modelling autoimmunity arising from immune response
to a viral infection.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

All authors listed have made a substantial, direct, and intellectual
contribution to the work and approved it for publication.

## Funding

OB has been supported by Cancer Research UK (EDDCPJT\100022).

## References

  * [1] H. I. Freedman, Deterministic Mathematical Models in Population Ecology.  New York: Marcel Dekker, Inc., 1980.
  * [2] G. F. Gause, The Struggle for Existence.  Baltimore: Williams and Wilkins, 1934.
  * [3] A. N. Kolmogorov, “A qualitative study of mathematical models of population dynamics,” (In Russian) Problemy Kibernetiki, vol. 25, p. 100, 1972.
  * [4] M. E. Solomon, “The natural control of animal populations,” Journal of Animal Ecology, vol. 18, pp. 1–35, 1949.
  * [5] A. J. Lotka, Elements of physical biology.  Baltimore: Williams and Wilkins, 1925.
  * [6] V. Volterra, “Fluctuations in the abundance of a species considered mathematically,” Nature, vol. 188, pp. 558–560, 1926.
  * [7] C. S. Holling, “The components of predation as revealed by a study of small-mammal predation of the european pine sawfly,” The Canadian Entomologist, vol. 91, pp. 293–320, 1959.
  * [8] C. S. Holling, “Some characteristics of simple types of predation and parasitism,” The Canadian Entomologist, vol. 91, pp. 385–398, 1959.
  * [9] B. Mondal, S. Roy, U. Ghosh, and P. K. Tiwari, “A systematic study of autonomous and nonautonomous predator–prey models for the combined effects of fear, refuge, cooperation and harvesting,” The European Physical Journal Plus, vol. 137, p. 724, 2022.
  * [10] S. Roy, P. K. Tiwari, H. Nayak, and M. Martcheva, “Effects of fear, refuge and hunting cooperation in a seasonally forced eco-epidemic model with selective predation,” The European Physical Journal Plus, vol. 137, p. 528, 2022\.
  * [11] S. Roy, P. K. Tiwari, H. Nayak, and M. Martcheva, “An ecoepidemic seasonally forced model for the combined effects of fear, additional foods and selective predation,” Journal of Biological Systems, vol. 30, pp. 285–321, 2022\.
  * [12] M. Hossain, S. Pal, P. K. Tiwari, and N. Pal, “Bifurcations, chaos, and multistability in a nonautonomous predator–prey model with fear,” Chaos, vol. 31, p. 123134, 2021.
  * [13] P. K. Tiwari, M. Verma, S. Pal, Y. Kang, and A. K. Misra, “A delay nonautonomous predator–prey model for the effects of fear, refuge and hunting cooperation,” Journal of Biological Systems, vol. 29, pp. 927–969, 2021.
  * [14] H. R. Akcakaya, “Population cycles of mammals: evidence for a ratio-dependent predation hypothesis,” Ecology, vol. 62, pp. 119–142, 1992.
  * [15] R. Arditi, L. R. Ginzburg, and H. R. Akcakaya, “Variation in plankton densities among lakes: a case for ratio-dependent predation models,” American Naturalist, vol. 138, pp. 1287–1296, 1991.
  * [16] A. P. Gutierrez, “Physiological basis of ratio-dependent predator-prey theory: The metabolic pool model as a paradigm,” Ecology, vol. 73, pp. 1552–1563, 1992.
  * [17] R. Arditi, N. Perrin, and H. Saïah, “Functional responses and heterogeneities: an experimental test with cladocerans,” Oikos, vol. 60, pp. 69–75, 1991.
  * [18] R. Arditi and L. R. Ginzburg, “Coupling in predator-prey dynamics: Ratio-dependence,” Journal of Theoretical Biology, vol. 139, pp. 311–326, 1989.
  * [19] R. Arditi and L. R. Ginzburg, How species interact: altering the standard view on trophic ecology.  Oxford: Oxford University Press, 2012.
  * [20] E. Beretta and Y. Kuang, “Global analyses in some delayed ratio-dependent predator-prey systems,” Nonlinear Analysis: Theory, Methods & Applications, vol. 32, pp. 381–408, 1998.
  * [21] C. Jost, O. Arino, and R. Arditi, “About deterministic extinction in ratio-dependent predator–prey models,” Bulletin of Mathematical Biology, vol. 61, pp. 19–32, 1999.
  * [22] D. Xiao and S. Ruan, “Global dynamics of a ratio-dependent predator-prey system,” Journal of Mathematical Biology, vol. 43, pp. 268–290, 2001.
  * [23] M. Haque, “Ratio-dependent predator-prey models of interacting populations,” Bulletin of Mathematical Biology, vol. 71, pp. 430–452, 2009.
  * [24] C. Arancibia-Ibarra, P. Aguirre, J. Flores, and P. Heijster, “Bifurcation analysis of a predator-prey model with predator intraspecific interactions and ratio-dependent functional response,” Applied Mathematics and Computation, vol. 402, p. 126152, 2021.
  * [25] P. Aguirre, J. D. Flores, and E. González-Olivares, “Bifurcations and global dynamics in a predator–prey model with a strong allee effect on the prey, and a ratio-dependent functional response,” Nonlinear Analysis: Real World Applications, vol. 16, pp. 235–249, 2014.
  * [26] P. J. Pal and T. Saha, “Qualitative analysis of a predator–prey system with double allee effect in prey,” Chaos, Solitons & Fractals, vol. 73, pp. 36–63, 2015.
  * [27] Y. V. Tyutyunov and L. I. Titova, “From Lotka–Volterra to Arditi–Ginzburg: 90 years of evolving trophic functions,” Biology Bulletin Reviews, vol. 10, pp. 167–185, 2020.
  * [28] J. H. P. Dawes and M. O. Souza, “A derivation of Holling’s type I, II and III functional responses in predator–prey systems,” Journal of Theoretical Biology, vol. 327, pp. 11–22, 2013.
  * [29] Y.-H. Fan and W.-T. Li, “Permanence for a delayed discrete ratio-dependent predator–prey system with holling type functional response,” Journal of Mathematical Analysis and Applications, vol. 299, pp. 357–374, 2004.
  * [30] P. J. Pal, P. K. Mandal, and K. K. Lahiri, “A delayed ratio-dependent predator–prey model of interacting populations with Holling type III functional response,” Nonlinear Dynamics, vol. 76, pp. 201–220, 2014.
  * [31] X. Wang, M. Peng, and X. Lu, “Stability and Hopf bifurcation analysis of a ratio-dependent predator–prey model with two time delays and Holling type III functional response,” Applied Mathematics and Computation, vol. 268, pp. 496–508, 2015.
  * [32] W.-T. Li and S.-L. Wu, “Traveling waves in a diffusive predator–prey model with Holling type-III functional response,” Chaos, Solitons & Fractals, vol. 37, pp. 476–486, 2008.
  * [33] N. Apreutesei and G. Dimitriu, “On a prey–predator reaction–diffusion system with Holling type III functional response,” Journal of Computational and Applied Mathematics, vol. 235, pp. 366–379, 2010.
  * [34] L. N. Guin and P. K. Mandal, “Spatial pattern in a diffusive predator–prey model with sigmoid ratio-dependent functional response,” International Journal of Biomathematics, vol. 7, p. 1450047, 2014.
  * [35] A. Morozov, “Emergence of Holling type III zooplankton functional response: Bringing together field evidence and mathematical modelling,” Journal of Theoretical Biology, vol. 265, pp. 45–54, 2010.
  * [36] D. Schenk, L.-F. Bersies, and S. Bacher, “An experimental test of the nature of predation: neither prey- nor ratio-dependent,” Journal of Animal Ecology, vol. 74, pp. 86–91, 2005.
  * [37] O. Sarnelle and A. E. Wilson, “Type III functional response in daphnia,” Ecology, vol. 89, pp. 1723–1732, 2008.
  * [38] P. Kratina, M. Vos, A. Bateman, and B. R. Anholt, “Functional responses modified by predator density,” Oecologia, vol. 159, pp. 425–433, 2009\.
  * [39] J. Holt, M. J. Jeger, J. M. Thresh, and G. W. Otim-Nape, “An epidemiological model incorporating vector population dynamics applied to african cassava mosaic virus disease,” Journal of Applied Ecology, vol. 34, pp. 793–806, 1997.
  * [40] K. B. Blyuss, F. A. Basir, V. A. Tsygankova, L. O. Biliavska, G. O. Iutynska, and et al., “Control of mosaic disease using microbial biostimulants: insights from mathematical modelling,” Ricerche di Matematica, vol. 69, pp. 437–455, 2020.
  * [41] V. S. Ivlev, Experimental ecology of the feeding of Fishes.  New Haven, CT: Yale University Press, 1947.
  * [42] W. E. Ricker, “Stock and recruitment,” Journal of the Fisheries Research Board of Canada, vol. 11, pp. 559–623, 1954.
  * [43] F. A. Basir, Y. N. Kyrychko, K. B. Blyuss, and S. Ray, “Effects of vector maturation time on the dynamics of Cassava Mosaic Disease,” Bulletin of Mathematical Biology, vol. 83, p. 87, 2021.
  * [44] K. B. Blyuss, S. N. Kyrychko, and Y. N. Kyrychko, “Time-delayed and stochastic effects in a predator-prey model with ratio dependence and Holling type III functional response,” Chaos, vol. 31, p. 073141, 2021.
  * [45] T. Zhang, T. Ding, W. Huang, and Z. Dong, Qualitative Theory of Differential Equations.  Providence: American Mathematical Society, 1992.
  * [46] F. Dumortier, J. Llibre, and J. Artés, Qualitative Theory of Planar Differential Systems.  Berlin: Springer-Verlag, 2006.
  * [47] F. Fatehi, S. Kyrychko, A. Ross, Y. N. Kyrychko, and K. B. Blyuss, “Stochastic effects in autoimmune dynamics,” Frontiers in Physiology, vol. 9, p. 45, 2018.
  * [48] F. Fatehi, Y. N. Kyrychko, and K. B. Blyuss, “Stochastic dynamics in a time-delayed model for autoimmunity,” Mathematical Biosciences, vol. 322, p. 108323, 2020.