# Collinear Laser Spectroscopy of $2\,{}^{3}\\!S_{1}\rightarrow
2\,{}^{3}\\!P_{\\!J}$ transitions in helium-like ${}^{12}\mathrm{C}^{4+}$

P. Imgram<EMAIL_ADDRESS>Institut für Kernphysik, Departement of
Physics, Technische Universität Darmstadt, Schlossgartenstraße 9, 64289
Darmstadt, Germany  K. König Institut für Kernphysik, Departement of Physics,
Technische Universität Darmstadt, Schlossgartenstraße 9, 64289 Darmstadt,
Germany  Helmholtz Research Academy Hesse for FAIR, Campus Darmstadt,
Schlossgartenstr. 9, 64289 Darmstadt  B. Maaß current address: Physics
Division, Argonne National Laboratory, 9700 S Cass Ave, IL 60439 Lemont, USA
Institut für Kernphysik, Departement of Physics, Technische Universität
Darmstadt, Schlossgartenstraße 9, 64289 Darmstadt, Germany  P. Müller Institut
für Kernphysik, Departement of Physics, Technische Universität Darmstadt,
Schlossgartenstraße 9, 64289 Darmstadt, Germany  W. Nörtershäuser Institut für
Kernphysik, Departement of Physics, Technische Universität Darmstadt,
Schlossgartenstraße 9, 64289 Darmstadt, Germany  Helmholtz Research Academy
Hesse for FAIR, Campus Darmstadt, Schlossgartenstr. 9, 64289 Darmstadt

###### Abstract

Transition frequencies and fine-structure splittings of the
$2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{\\!J}$ transitions in helium-like
${}^{12}\mathrm{C}^{4+}$ were measured by collinear laser spectroscopy on a
1-ppb level. Accuracy is increased by more than three orders of magnitude with
respect to previous measurements, enabling tests of recent non-relativistic
QED calculations including terms up to $m\alpha^{7}$. Deviations between the
theoretical and experimental values are within theoretical uncertainties and
are ascribed to $m\alpha^{8}$ and higher-order contributions in the series
expansion of the NR-QED calculations. Finally, prospects for an all-optical
charge radius determination of light isotopes are evaluated.

## Introduction. –

In the last decade, advances in experimental precision and sensitivity in
laser spectroscopy paired with ab-initio atomic structure calculations have
enabled the determination of nuclear charge radii in an all-optical approach
purely from optical transition frequencies and related QED calculations. Laser
spectroscopy of muonic systems in particular achieved remarkable results in
$\mu$H Pohl  _et al._ (2010), $\mu$D Pohl  _et al._ (2016) and $\mu$He Krauth
_et al._ (2021); The CREMA Collaboration  _et al._ (2023). The former
suggested a proton radius about 4% smaller than previously determined from
elastic electron scattering and atomic spectroscopy and led to the famous
proton-radius puzzle and questioned the lepton universality proposed by the
Standard Model of fundamental interactions (for a recent review see, e.g., Gao
and Vanderhaeghen (2022)). Although there is convincing evidence for the
smaller proton radius as observed in muonic hydrogen from improved
measurements of atomic hydrogen Beyer  _et al._ (2017), a new lamb-shift
measurement Bezginov  _et al._ (2019) as well as elastic electron scattering
at very small forward-scattering angles using a windowless target Xiong  _et
al._ (2019), some results are still under discussion and atomic hydrogen
experiments that support the larger proton radius are also reported Fleurbaey
_et al._ (2018); Brandt  _et al._ (2022). Therefore, an extended comparison
between muonic and electronic systems is demanded. However, only the nuclear
charge radius of H was studied so far in both systems Pohl  _et al._ (2010);
Udem  _et al._ (1997).
In order to expand the all-optical approach towards heavier atomic systems
than hydrogen, precise atomic structure calculations for two-electron systems
are needed since no laser-addressable narrow electric-dipole transition exists
in H-like systems beyond hydrogen. Towards this goal, non-relativistic QED
(NR-QED) calculations in He-like systems recently made significant progress
Yerokhin and Pachucki (2010); Patkóš  _et al._ (2021); Yerokhin  _et al._
(2022a, b, 2023). This approach is based on a perturbative expansion series of
the level energy in orders $n$ of $m\alpha^{n}$ with the electron mass $m$ and
the fine-structure constant $\alpha$. While the calculated
$1s2s\,^{3}S_{1}\rightarrow 1s2p\,^{3}P_{J}$ transition energies (in the
following text abbreviated as $2\,^{3}S_{1}$ and $2\,^{3}P_{J}$) are very
consistent with experimental data in He, the ionization energies of the
individual states as well as the transition energies to the $3\,^{3}\\!D$
states differ by up to $10\sigma$ Patkóš  _et al._ (2021); Clausen  _et al._
(2021); Yerokhin  _et al._ (2023). This could be caused by unknown theoretical
contribution shifting the $2\,^{3}\\!S$ and $2\,^{3}\\!P$ states by about the
same amount Patkóš  _et al._ (2021). Since QED contributions scale as $\sim
Z^{4}$, measurements in He-like systems of higher $Z$ have an increased
sensitivity and the determination of $2\,{}^{3}\\!S_{1}\rightarrow
2\,{}^{3}\\!P_{\\!J}$ transition frequencies and $2\,^{3}\\!P_{J}$ fine-
structure splittings in these ions might provide a hint for the origin of
these inconsistencies. So far, only transition frequency measurements in Li+
Riis  _et al._ (1994) and Be2+ Scholl  _et al._ (1993) reached the necessary
accuracy to test NR-QED calculations Yerokhin  _et al._ (2022b). While theory
and experiment agree for Li+, there is a significant discrepancy in Be2+
Scholl  _et al._ (1993); Pachucki and Yerokhin (2010); Yerokhin  _et al._
(2023). Helium-like carbon 12C4+ is an excellent candidate for such a test. It
is the first stable isotope beyond He that has no nuclear spin and is
therefore not plagued by hyperfine-induced fine-structure mixing. Moreover,
its nuclear charge radius is accurately known from elastic electron scattering
Cardman  _et al._ (1980); Sick (1982); Reuter  _et al._ (1982); Offermann  _et
al._ (1991) as well as muonic atom spectroscopy Schaller  _et al._ (1982);
Ruckstuhl  _et al._ (1984), yielding consistent results. Thus, spectroscopy of
the $2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{\\!J}$ transitions in He-like
carbon provides the ultimate testing ground for higher-order terms in NR-QED
calculations as well as for the evaluation of the prospects to extract all-
optical nuclear charge radii from He-like ions beyond He.

Figure 1: (a) Level scheme of 12C4+. (b) Experimental setup of the COALA
beamline. Highly charged 12C4+ ions are produced in an electron beam ion
source (EBIS) on a 10.5 kV starting potential. After acceleration of the ions
into the beamline, they are superposed with a collinear and an anticollinear
laser beam to perform quasi-simultaneous collinear and anticollinear laser
spectroscopy in the optical detection region where the fluorescence light is
collected and detected by two elliptical mirrors and photomultiplier tubes.

In this work, we present $2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{\\!J}$
transition frequency measurements in 12C4+ with 1.3 part-per-billion (ppb)
accuracy which represents an improvement by three orders of magnitude compared
to previous experiments. Our results surpass the theoretical accuracy by
roughly two orders of magnitude and therefore become sensitive to higher-order
QED terms in NR-QED beyond He.

## Method. –

The level structure of He-like ions can be subdivided in singlet ($S=0$) and
triplet ($S=1$) states where the electron spins are aligned anti-parallel or
parallel, respectively. As illustrated in Fig. 1(a), in 12C4+ no laser-
accessible transition exists from the $1s^{2}\,{}^{1}\\!S_{0}$ ground state.
However, the metastable $2\,^{3}\\!S_{1}$ triplet state can serve as the lower
state for laser excitation into the $2\,^{3}\\!P_{J}$ states at a laser
wavelength of $\lambda\approx 227$ nm if it is sufficiently populated in the
ion source. Here, we use an electron beam ion source (EBIS) which is well
known to produce highly charged ions and to populate metastable states. But
since spectroscopy inside an EBIS is limited by the very high temperature and
the corresponding Doppler broadening even if forced evaporative cooling is
applied Mäckel  _et al._ (2011), it is necessary to transport the ions into an
environment that is more appropriate for high-resolution spectroscopy. This
can boost accuracy by several orders of magnitude as it has been demonstrated
for Ar13+ ions in a Penning Egl  _et al._ (2019) or a Paul trap Micke  _et
al._ (2020); King  _et al._ (2022). However, the short lifetime of the
$2\,^{3}\\!S_{1}$ state in He-like C4+ ions ($\tau\approx 21$ ms) requires a
fast technique like collinear laser spectroscopy, which had originally been
developed to perform spectroscopy of short-lived isotopes Schinzler  _et al._
(1978). This was realized by coupling an EBIS (DREEBIT, EBIS-A) to the
Collinear Apparatus for Laser Spectroscopy and Applied Physics (COALA) König
_et al._ (2020) situated at the Institute for Nuclear Physics of TU Darmstadt
as depicted in Fig. 1(b). At COALA, we achieved accuracies of the order of 100
kHz in previous studies on singly charged ions Imgram  _et al._ (2019); Müller
_et al._ (2020); König  _et al._ (2020), comparable to ion-trap measurements
of allowed dipole transitions Herrmann  _et al._ (2009); Batteiger  _et al._
(2009); Gebert  _et al._ (2015); Shi  _et al._ (2016). The He-like ions are
produced through electron impact ionization and a significant fraction
($\approx 40\%$) of the ions is transferred into the metastable state almost
exclusively through charge exchange from C5+ to C4+. The initial energy spread
strongly depends on the production parameters as detailed in Imgram  _et al._
(2023) and was $1.1\,\mathrm{eV}$ in the applied continous-beam mode. Then,
the ions are accelerated from their starting potential
$U_{\mathrm{start}}\approx 10.5$ kV towards ground potential into the
beamline. The choice of $U_{\mathrm{start}}$ was a compromise of having a high
starting potential for maximum velocity compression while limiting recurrences
of discharges inside the source. Through subsequent electrostatic ion optics
and deflectors, the ion beam is superposed with two laser beams, of which one
is copropagating with the ion beam and the other one counterpropagating. By
alternately blocking one of the two laser beams, we perform frequency‑comb
referenced quasi-simultaneous collinear and anticollinear laser spectroscopy
as it was previously applied for short-lived Be isotopes Nörtershäuser  _et
al._ (2009); Krieger  _et al._ (2016). The resonance fluorescence signal is
recorded with photomultiplier tubes (PMT) attached to the optical detection
region that can be floated to a variable potential to realize Doppler tuning.
Hereby, the Doppler-shifted laboratory-frame transition frequencies
$\nu_{\mathrm{c/a}}=\nu_{0}\gamma(1\pm\beta)$ in collinear
($\nu_{\mathrm{c}}$) and anticollinear ($\nu_{\mathrm{a}}$) geometry are
measured in fast iterations. The rest-frame transition frequency $\nu_{0}$ is
extracted as the geometrical average

$\nu_{\mathrm{c}}\nu_{\mathrm{a}}=\nu_{0}^{2}\gamma^{2}(1+\beta)(1-\beta)=\nu_{0}^{2}$
(1)

without the need of the knowledge of the precise ion velocity
$\beta=\upsilon/c$ and the corresponding Lorentz factor
$\gamma=1/\sqrt{1-\beta^{2}}$. This reduces the systematic uncertainties,
which are then dominated by the uncertainty stemming from the alignment of the
laser beams as discussed in detail in Imgram  _et al._ (2023).
The laser system König  _et al._ (2020) consists of two continuous-wave
titanium:sapphire (Ti:Sa) lasers (Matisse 2), each pumped by a frequency-
doubled Nd:YAG laser (Millennia eV20). The emitted 908 nm light is frequency
quadrupled and transported in free space to the beamline. Lens telescopes are
used to achieve a collimated laser beam in both directions with a beam
diameter of about 1 mm inside the optical detection region. The laser power of
both lasers was held constant at 0.5 mW during the measurements. A GPS-
disciplined quartz oscillator provided the 10 MHz reference for the Menlo-
Systems FC1500-250-WG frequency comb used to measure and stabilize the
fundamental laser frequency. The resulting laser linewidth of the fundamental
laser light was approximately 200 kHz. The laser light had linear polarization
to suppress potential Zeeman shifts induced by external magnetic fields Imgram
_et al._ (2023).

## Transition frequencies. –

A typical resonance spectrum depicted in Fig. 2(a) is obtained by plotting the
PMT counts as a function of the laser frequency. To visualize the signal-to-
background ratio, the photon counts are normalized to the background rate. In
order to extract the center frequency $\nu_{\mathrm{c/a}}$ from the spectrum,
a Voigt profile, which is the convolution of a Gaussian and Lorentzian
profile, was fitted to the data. The resulting full width at half maximum
(FWHM) of roughly 170 MHz is dominated by the Gaussian contribution which
originates from the energy spread of the ions Imgram  _et al._ (2023). The
choice of the fitting lineshape does not shift the result significantly as
long as the same profile is used for both directions. Each pair consisting of
a collinear and an anticollinear measurement yields a rest-frame frequency
through Eq. (1). The final frequencies $\nu_{0}(\,^{3}\\!P_{J})$ of the
$2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{\\!J}$ transitions were
determined by taking the average over all pairs weighted by their statistical
fitting uncertainties. Results are listed in Tab. 1 and compared with previous
experimental values and theoretical predictions.

Table 1: Comparison between literature values and this work of the rest-frame frequencies $\nu_{0}(\,^{3}\\!P_{J})$ in the $2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{\\!J}$ transitions of 12C4+. The values in brackets denote the full $1\sigma$ uncertainty. All values are in GHz. Label Fig. 2 | $\nu_{0}(\,^{3}\\!P_{2})$ | $\nu_{0}(\,^{3}\\!P_{1})$ | $\nu_{0}(\,^{3}\\!P_{0})$ | Ref.
---|---|---|---|---
(A) | $1\,319\,747.\,(29)$ | $1\,315\,669.\,(29)$ | $1\,316\,056.\,(29)$ | [Edlen and Lofstrand, 1970] (Exp.)
(B) | $1\,319\,744.6\,(4.8)$ | $1\,315\,076.2\,(4.8)$ | $1\,316\,057.7\,(4.8)$ | [Drake, 1988] (Theory)
(C) | $1\,319\,753.0\,(3.6)$ | $1\,315\,679.7\,(4.2)$ | $1\,316\,056.2\,(4.5)$ | [Ozawa  _et al._ , 2001] (Exp.)
(D) | $1\,319\,748.6\,(1.0)$ | $1\,315\,677.5\,(1.0)$ | $1\,316\,052.5\,(1.0)$ | [Yerokhin and Pachucki, 2010] (Theory)
(E) | $1\,319\,748.55\,(13)$ | $1\,315\,677.60\,(75)$ | $1\,316\,052.03\,(27)$ | [Yerokhin  _et al._ , 2022b] (Theory)
| $1\,319\,748.571\,4\,(17)$ | $1\,315\,677.192\,8\,(17)$ | $1\,316\,052.219\,3\,(19)$ | This work

Statistical uncertainties are estimated as the standard error of the mean for
108 (${}^{3}\\!P_{2}$), 68 (${}^{3}\\!P_{1}$), and 28 (${}^{3}\\!P_{0}$)
measurements that were taken over several weeks. Our total systematic
uncertainty of the rest-frame transition frequency is $\sim$1.7 MHz, dominated
by the remaining ion beam divergence in combination with the superposition of
the two laser beams that might lead to probing slightly different velocity
classes in both directions (for details see Imgram  _et al._ (2023)). The
systematic uncertainty is added in quadrature to the statistical uncertainty
and a combined 1$\sigma$ uncertainty of less than 2 MHz is obtained for all
transitions. This represents an improvement of more than three orders of
magnitude compared to previous experimental results in 12C4+ Ozawa  _et al._
(2001) as illustrated in Fig. 2(b). Here, the relative accuracy
$\Delta\nu_{0}/\nu_{0}$ in experiment (blue) and theory (red) is plotted on a
logarithmic axis as a function of time. A steady improvement in accuracy is
visible with a substantial gain through our work. Our experimental results
agree very well with recent NR-QED calculations Yerokhin  _et al._ (2022b)
whose uncertainties are, however, two orders of magnitude larger than our
experimental uncertainties.

## Fine-structure splitting. –

Another important test of theory is the fine-structure splitting of the
$2\,^{3}\\!P_{J}$ states. While the theoretical accuracy of the
$2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{\\!J}$ transition frequencies in
12C4+ is currently limited to $\sim$130 MHz Yerokhin  _et al._ (2022b), the
fine-structure splittings in 12C4+ were calculated on the $\sim$10 MHz level
Pachucki and Yerokhin (2010) including $m\alpha^{7}$ terms, but could never be
tested hitherto because of the lack of experimental data. The fine structure
splittings can be directly obtained from the differences of our measured
transition frequencies and the values are listed in Tab. 2 together with the
theoretical predictions by Pachucki and Yerokhin Pachucki and Yerokhin (2010).
All values agree within the stated theoretical uncertainties. The difference
between the experimental and the theoretical value represents the sum of all
higher-order contributions and are, thus, an approximation of the dominant
$m\alpha^{8}$ terms. Their contribution was estimated by scaling the
calculated $m\alpha^{6}$ contribution by the factor $(Z\alpha)^{2}$ and used
to represent the theoretical uncertainty Yerokhin  _et al._ (2022b). In atomic
helium the corresponding experimental values Smiciklas and Shiner (2010);
Zheng  _et al._ (2017) were in excellent agreement with the theoretical result
and, thus, the $m\alpha^{8}$ contribution to the splitting between the unmixed
$J=0$ and $J=2$ levels were expected to be small. Confirming the
overestimation of the $m\alpha^{8}$ contributions in the more sensitive case
of 12C4+ would have allowed to revise the uncertainty estimation and, thus, to
improve the uncertainty of the fine-structure constant $\alpha$ Pachucki and
Yerokhin (2010). However, our experimental result indicates a rather large
contribution of higher orders. Instead, our result can now be used as a
benchmark for calculations that approximate the $m\alpha^{8}$ contributions to
identify the dominant terms, which finally can also lead to an improved He
fine-structure constant.

Finally, we note the larger difference between experiment and theory for the
${}^{3}\\!P_{0}-{{}^{3}\\!P_{1}}$ interval compared with the
${}^{3}\\!P_{0}-{{}^{3}\\!P_{2}}$ interval. This is most probably due to the
singlet-triplet configuration mixing of the $2\,^{3}\\!P_{1}$ and the
$2\,^{1}\\!P_{1}$ states having the same total angular momentum and parity and
are close in energy Yerokhin  _et al._ (2022b). This mixing is enhanced with
increasing $Z$ as compared to other $m\alpha^{8}$ effects which is also
reflected in the larger uncertainty from theory.

Table 2: Fine-structure splitting of the $2\,^{3}\\!P_{J}$ levels in 12C4+ acquired from this work and theory Pachucki and Yerokhin (2010). All values are in MHz. $(J,J^{\prime})$ | $(0,1)$ | $(0,2)$
---|---|---
This work | $-375\,026.5(2.5)$ | $3\,696\,352.1(2.5)$
Theory Pachucki and Yerokhin (2010) | $-374\,996.3\,(48.0)$ | $3\,696\,343.5\,(10.0)$
Difference | $-30.2\,(48.0)$ | $8.6\,(10.0)$
Figure 2: (a) Background-normalized anticollinear resonance spectrum of the
$2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{2}$ transition in 12C4+ with a
Voigt profile fit (red line) and the respective residuals. (b) Relative
precision $\Delta\nu_{0}/\nu_{0}$ of literature values (A) - (E) from Tab. 1
for the $2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{2}$ transition as a
function of time including our result. (c) Literature values (F) - (K) of the
nuclear charge radius of 12C from Tab. 3 in comparison to the all-optical
nuclear charge radius extracted here.

## Charge Radii. –

By confirming the NR-QED calculations in two-electron systems, these can be
combined with the experimental data to determine absolute nuclear charge radii
in an all-optical approach. The measured transition frequency $\nu_{0}$ of the
$2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{2}$ transition can be written as

$\nu_{0}=\nu_{\mathrm{point}}+F\cdot{R_{\mathrm{c}}^{2}}$ (2)

where $\nu_{\mathrm{point}}$ is the calculated transition frequency assuming a
point-like but finite-mass nucleus and $F$ the calculated field-shift factor
of the transition. Both values are calculated by Yerokhin et al. but not
explicitly tabulated in Yerokhin  _et al._ (2022b). Instead, we obtained
$\nu_{\mathrm{point}}=1\,319\,749.83(13)$ GHz as the difference of the
ionization energies of the $2\,^{3}\\!S_{1}$ and the $2\,^{3}\\!P_{2}$ levels
after subtracting the nuclear size contribution given also in the tables.
Dividing the latter by $R_{\mathrm{c}}^{2}$ using
$R_{\mathrm{c}}(^{12}\mathrm{C})=2.4702(22)$ fm as applied in Yerokhin  _et
al._ (2022b) provides $F=0.2115$ GHz/fm2. Using these values gives
$R_{\mathrm{c}}(^{12}\mathrm{C})=2.45\,(12)\,\mathrm{fm}$ in excellent
agreement with previous results as listed in Tab. 3. In case that
$\Delta\nu_{\mathrm{point}}$ can be reduced to the size of
$\Delta\nu_{0,\mathrm{exp}}$, the charge radius could be directly obtained
with an unprecedented accuracy of $\pm 0.0016$ fm, which is even more accurate
than the determination from muonic atom spectroscopy Ruckstuhl  _et al._
(1984). While this is a daunting task for theory, we note that for the stable
isotopes of boron, already an improvement of the theoretical uncertainty by a
factor of 2–3 would allow for an improved charge radius determination since
they are only poorly known from elastic electron scattering and muonic atom
spectroscopy Maaß  _et al._ (2019). A more reliable charge radius of stable
boron isotopes is of particular importance for the ongoing campaign to
determine the charge radius of the proton-halo candidate 8B Maaß  _et al._
(2017).

Table 3: Nuclear charge radius of 12C. All values are in fm. | Rc(12C) | Ref.
---|---|---
(F) | 2.472(15) | [Cardman  _et al._ , 1980]
(G) | 2.471(6) | [Sick, 1982]
(H) | 2.464(12) | [Reuter  _et al._ , 1982]
(I) | 2.478(9) | [Offermann  _et al._ , 1991]
(J) | 2.472(16) | [Schaller  _et al._ , 1982]
(K) | 2.4829(19) | [Ruckstuhl  _et al._ , 1984]
| 2.45(12) | This work

## Summary. –

We measured the transition frequencies of all fine-structure components in the
$2\,{}^{3}\\!S_{1}\rightarrow 2\,{}^{3}\\!P_{\\!J}$ transitions in He-like
12C4+ ions with 1.3 ppb accuracy and confirmed NR-QED calculations of these
transitions within their uncertainties. We combined the production of He-like
ions in the metastable $2\,^{3}\\!S_{1}$ state in an EBIS with collinear laser
spectroscopy and demonstrated frequency determinations at the 1-MHz level. The
experimentally determined fine-structure splittings of the ${}^{3}\\!P_{J}$
states pave the way to further decrease the theoretical uncertainty of the
fine-structure constant, required for constraining the fine-structure constant
in He. Based on the NR-QED calculations, the charge radius of 12C was
extracted directly from the optical transition frequency with an uncertainty
of $\sim$5%, only limited by the theoretical accuracy. The technique will be
applied next to He-like 13C ions. Measuring the hyperfine splitting and
isotope shift with respect to 12C, a very accurate charge radius of 13C can be
obtained based on conventional mass-shift calculations as they have been used
for determinations of charge radii of He, Li, Be Lu  _et al._ (2013) and B
isotopes Maaß  _et al._ (2019).
Finally, measurements on He-like and Li-like boron ions and neutral boron Maaß
_et al._ (2019) will provide a very precise differential charge radius between
the B isotopes allowing to test the consistency of NR-QED calculations in
two-, three-, and five-electron systems.

###### Acknowledgements.

## Acknowledgment. –

We thank J. Krämer for his contributions in the early stage of the project and
K. Pachucki & V. Yerokhin for many fruitful discussions. We acknowledge
funding by the Deutsche Forschungsgemeinschaft (DFG) Project No. 279384907 SFB
1245, as well as under Grant INST No. 163/392-1 FUGG, and from the German
Federal Ministry for Education and Research (BMBF) under Contract No.
05P21RDFN1. P.I. and P.M. acknowledge support from HGS-HIRE.

## References

  * Pohl  _et al._ (2010) R. Pohl, A. Antognini, F. Nez, F. D. Amaro, F. Biraben, J. M. R. Cardoso, D. S. Covita, A. Dax, S. Dhawan, L. M. P. Fernandes, A. Giesen, T. Graf, T. W. Hänsch, P. Indelicato, L. Julien, C.-Y. Kao, P. Knowles, E.-O. Le Bigot, Y.-W. Liu, J. A. M. Lopes, L. Ludhova, C. M. B. Monteiro, F. Mulhauser, T. Nebel, P. Rabinowitz, J. M. F. dos Santos, L. A. Schaller, K. Schuhmann, C. Schwob, D. Taqqu, J. F. C. A. Veloso, and F. Kottmann, Nature 466, 213 (2010).
  * Pohl  _et al._ (2016) R. Pohl, F. Nez, L. M. P. Fernandes, F. D. Amaro, F. Biraben, J. M. R. Cardoso, D. S. Covita, A. Dax, S. Dhawan, M. Diepold, A. Giesen, A. L. Gouvea, T. Graf, T. W. Hänsch, P. Indelicato, L. Julien, P. Knowles, F. Kottmann, E.-O. L. Bigot, Y.-W. Liu, J. A. M. Lopes, L. Ludhova, C. M. B. Monteiro, F. Mulhauser, T. Nebel, P. Rabinowitz, J. M. F. dos Santos, L. A. Schaller, K. Schuhmann, C. Schwob, D. Taqqu, J. F. C. A. Veloso, and A. Antognini, Science 353, 669 (2016).
  * Krauth  _et al._ (2021) J. J. Krauth, K. Schuhmann, M. A. Ahmed, F. D. Amaro, P. Amaro, F. Biraben, T.-L. Chen, D. S. Covita, A. J. Dax, M. Diepold, L. M. P. Fernandes, B. Franke, S. Galtier, A. L. Gouvea, J. Götzfried, T. Graf, T. W. Hänsch, J. Hartmann, M. Hildebrandt, P. Indelicato, L. Julien, K. Kirch, A. Knecht, Y.-W. Liu, J. Machado, C. M. B. Monteiro, F. Mulhauser, B. Naar, T. Nebel, F. Nez, J. M. F. dos Santos, J. P. Santos, C. I. Szabo, D. Taqqu, J. F. C. A. Veloso, J. Vogelsang, A. Voss, B. Weichelt, R. Pohl, A. Antognini, and F. Kottmann, Nature 589, 527 (2021).
  * The CREMA Collaboration  _et al._ (2023) The CREMA Collaboration, K. Schuhmann, L. M. P. Fernandes, F. Nez, M. A. Ahmed, F. D. Amaro, P. Amaro, F. Biraben, T.-L. Chen, D. S. Covita, A. J. Dax, M. Diepold, B. Franke, S. Galtier, A. L. Gouvea, J. Götzfried, T. Graf, T. W. Hänsch, M. Hildebrandt, P. Indelicato, L. Julien, K. Kirch, A. Knecht, F. Kottmann, J. J. Krauth, Y.-W. Liu, J. Machado, C. M. B. Monteiro, F. Mulhauser, B. Naar, T. Nebel, J. M. F. dos Santos, J. P. Santos, C. I. Szabo, D. Taqqu, J. F. C. A. Veloso, A. Voss, B. Weichelt, A. Antognini, and R. Pohl, “The helion charge radius from laser spectroscopy of muonic helium-3 ions,” (2023), arXiv:2305.11679 [physics.atom-ph] .
  * Gao and Vanderhaeghen (2022) H. Gao and M. Vanderhaeghen, Rev. Mod. Phys. 94, 015002 (2022).
  * Beyer  _et al._ (2017) A. Beyer, L. Maisenbacher, A. Matveev, R. Pohl, K. Khabarova, A. Grinin, T. Lamour, D. C. Yost, T. W. Hänsch, N. Kolachevsky, and T. Udem, Science 358, 79 (2017).
  * Bezginov  _et al._ (2019) N. Bezginov, T. Valdez, M. Horbatsch, A. Marsman, A. C. Vutha, and E. A. Hessels, Science 365, 1007 (2019).
  * Xiong  _et al._ (2019) W. Xiong, A. Gasparian, H. Gao, D. Dutta, M. Khandaker, N. Liyanage, E. Pasyuk, C. Peng, X. Bai, L. Ye, K. Gnanvo, C. Gu, M. Levillain, X. Yan, D. W. Higinbotham, M. Meziane, Z. Ye, K. Adhikari, B. Aljawrneh, H. Bhatt, D. Bhetuwal, J. Brock, V. Burkert, C. Carlin, A. Deur, D. Di, J. Dunne, P. Ekanayaka, L. El-Fassi, B. Emmich, L. Gan, O. Glamazdin, M. L. Kabir, A. Karki, C. Keith, S. Kowalski, V. Lagerquist, I. Larin, T. Liu, A. Liyanage, J. Maxwell, D. Meekins, S. J. Nazeer, V. Nelyubin, H. Nguyen, R. Pedroni, C. Perdrisat, J. Pierce, V. Punjabi, M. Shabestari, A. Shahinyan, R. Silwal, S. Stepanyan, A. Subedi, V. V. Tarasov, N. Ton, Y. Zhang, and Z. W. Zhao, Nature 575, 147 (2019).
  * Fleurbaey  _et al._ (2018) H. Fleurbaey, S. Galtier, S. Thomas, M. Bonnaud, L. Julien, F. Biraben, F. Nez, M. Abgrall, and J. Guéna, Phys. Rev. Lett. 120, 183001 (2018).
  * Brandt  _et al._ (2022) A. D. Brandt, S. F. Cooper, C. Rasor, Z. Burkley, A. Matveev, and D. C. Yost, Phys. Rev. Lett. 128, 023001 (2022).
  * Udem  _et al._ (1997) T. Udem, A. Huber, B. Gross, J. Reichert, M. Prevedelli, M. Weitz, and T. W. Hänsch, Phys. Rev. Lett. 79, 2646 (1997).
  * Yerokhin and Pachucki (2010) V. A. Yerokhin and K. Pachucki, Phys. Rev. A 81, 022507 (2010).
  * Patkóš  _et al._ (2021) V. Patkóš, V. A. Yerokhin, and K. Pachucki, Phys. Rev. A 103, 042809 (2021).
  * Yerokhin  _et al._ (2022a) V. A. Yerokhin, V. Patkóš, and K. Pachucki, The European Physical Journal D 76, 142 (2022a).
  * Yerokhin  _et al._ (2022b) V. A. Yerokhin, V. Patkóš, and K. Pachucki, Phys. Rev. A 106, 022815 (2022b).
  * Yerokhin  _et al._ (2023) V. A. Yerokhin, V. Patkóš, and K. Pachucki, Phys. Rev. A 107, 012810 (2023).
  * Clausen  _et al._ (2021) G. Clausen, P. Jansen, S. Scheidegger, J. A. Agner, H. Schmutz, and F. Merkt, Phys. Rev. Lett. 127, 093001 (2021).
  * Riis  _et al._ (1994) E. Riis, A. G. Sinclair, O. Poulsen, G. W. F. Drake, W. R. C. Rowley, and A. P. Levick, Phys. Rev. A 49, 207 (1994).
  * Scholl  _et al._ (1993) T. J. Scholl, R. Cameron, S. D. Rosner, L. Zhang, R. A. Holt, C. J. Sansonetti, and J. D. Gillaspy, Phys. Rev. Lett. 71, 2188 (1993).
  * Pachucki and Yerokhin (2010) K. Pachucki and V. A. Yerokhin, Phys. Rev. Lett. 104, 070403 (2010).
  * Cardman  _et al._ (1980) L. Cardman, J. Lightbody, S. Penner, S. Fivozinsky, X. Maruyama, W. Trower, and S. Williamson, Physics Letters B 91, 203 (1980).
  * Sick (1982) I. Sick, Physics Letters B 116, 212 (1982).
  * Reuter  _et al._ (1982) W. Reuter, G. Fricke, K. Merle, and H. Miska, Phys. Rev. C 26, 806 (1982).
  * Offermann  _et al._ (1991) E. A. J. M. Offermann, L. S. Cardman, C. W. de Jager, H. Miska, C. de Vries, and H. de Vries, Phys. Rev. C 44, 1096 (1991).
  * Schaller  _et al._ (1982) L. Schaller, L. Schellenberg, T. Phan, G. Piller, A. Ruetschi, and H. Schneuwly, Nuclear Physics A 379, 523 (1982).
  * Ruckstuhl  _et al._ (1984) W. Ruckstuhl, B. Aas, W. Beer, I. Beltrami, K. Bos, P. Goudsmit, H. Leisi, G. Strassner, A. Vacchi, F. De Boer, U. Kiebele, and R. Weber, Nuclear Physics A 430, 685 (1984).
  * Mäckel  _et al._ (2011) V. Mäckel, R. Klawitter, G. Brenner, J. R. Crespo López-Urrutia, and J. Ullrich, Phys. Rev. Lett. 107, 143002 (2011).
  * Egl  _et al._ (2019) A. Egl, I. Arapoglou, M. Höcker, K. König, T. Ratajczyk, T. Sailer, B. Tu, A. Weigel, K. Blaum, W. Nörtershäuser, and S. Sturm, Phys. Rev. Lett. 123, 123001 (2019).
  * Micke  _et al._ (2020) P. Micke, T. Leopold, S. A. King, E. Benkler, L. J. Spieß, L. Schmöger, M. Schwarz, J. R. Crespo López-Urrutia, and P. O. Schmidt, Nature 578, 60 (2020).
  * King  _et al._ (2022) S. A. King, L. J. Spieß, P. Micke, A. Wilzewski, T. Leopold, E. Benkler, R. Lange, N. Huntemann, A. Surzhykov, V. A. Yerokhin, J. R. Crespo López-Urrutia, and P. O. Schmidt, Nature 611, 43 (2022).
  * Schinzler  _et al._ (1978) B. Schinzler, W. Klempt, S. Kaufman, H. Lochmann, G. Moruzzi, R. Neugart, E.-W. Otten, J. Bonn, L. Von Reisky, K. Spath, J. Steinacher, and D. Weskott, Physics Letters B 79, 209 (1978).
  * König  _et al._ (2020) K. König, J. Krämer, C. Geppert, P. Imgram, B. Maaß, T. Ratajczyk, and W. Nörtershäuser, Review of Scientific Instruments 91, 081301 (2020).
  * Imgram  _et al._ (2019) P. Imgram, K. König, J. Krämer, T. Ratajczyk, R. A. Müller, A. Surzhykov, and W. Nörtershäuser, Phys. Rev. A 99, 012511 (2019).
  * Müller  _et al._ (2020) P. Müller, K. König, P. Imgram, J. Krämer, and W. Nörtershäuser, Phys. Rev. Research 2, 043351 (2020).
  * König  _et al._ (2020) K. König, J. Krämer, P. Imgram, B. Maaß, W. Nörtershäuser, and T. Ratajczyk, Phys. Rev. A 102, 042802 (2020).
  * Herrmann  _et al._ (2009) M. Herrmann, V. Batteiger, S. Knünz, G. Saathoff, T. Udem, and T. W. Hänsch, Phys. Rev. Lett. 102, 013006 (2009).
  * Batteiger  _et al._ (2009) V. Batteiger, S. Knünz, M. Herrmann, G. Saathoff, H. A. Schüssler, B. Bernhardt, T. Wilken, R. Holzwarth, T. W. Hänsch, and T. Udem, Phys. Rev. A 80, 022503 (2009).
  * Gebert  _et al._ (2015) F. Gebert, Y. Wan, F. Wolf, C. N. Angstmann, J. C. Berengut, and P. O. Schmidt, Phys. Rev. Lett. 115, 053003 (2015).
  * Shi  _et al._ (2016) C. Shi, F. Gebert, C. Gorges, S. Kaufmann, W. Nörtershäuser, B. K. Sahoo, A. Surzhykov, V. A. Yerokhin, J. C. Berengut, F. Wolf, J. C. Heip, and P. O. Schmidt, Applied Physics B 123, 2 (2016).
  * Imgram  _et al._ (2023) P. Imgram, K. König, B. Maaß, P. Müller, and W. Nörtershäuser, Phys. Rev. A , in print (2023).
  * Nörtershäuser  _et al._ (2009) W. Nörtershäuser, D. Tiedemann, M. Žáková, Z. Andjelkovic, K. Blaum, M. L. Bissell, R. Cazan, G. W. F. Drake, C. Geppert, M. Kowalska, J. Krämer, A. Krieger, R. Neugart, R. Sánchez, F. Schmidt-Kaler, Z.-C. Yan, D. T. Yordanov, and C. Zimmermann, Phys. Rev. Lett. 102, 062503 (2009).
  * Krieger  _et al._ (2016) A. Krieger, W. Nörtershäuser, C. Geppert, K. Blaum, M. L. Bissell, N. Frömmgen, M. Hammen, K. Kreim, M. Kowalska, J. Krämer, R. Neugart, G. Neyens, R. Sánchez, D. Tiedemann, D. T. Yordanov, and M. Zakova, Applied Physics B 123, 15 (2016).
  * Edlen and Lofstrand (1970) B. Edlen and B. Lofstrand, Journal of Physics B: Atomic and Molecular Physics 3, 1380 (1970).
  * Drake (1988) G. W. Drake, Canadian Journal of Physics 66, 586 (1988).
  * Ozawa  _et al._ (2001) S. Ozawa, T. Ariga, N. Inabe, M. Kase, I. Tanihata, M. Wakasugi, and Y. Yano, Physica Scripta T92, 195 (2001).
  * Smiciklas and Shiner (2010) M. Smiciklas and D. Shiner, Phys. Rev. Lett. 105, 123001 (2010).
  * Zheng  _et al._ (2017) X. Zheng, Y. R. Sun, J.-J. Chen, W. Jiang, K. Pachucki, and S.-M. Hu, Phys. Rev. Lett. 119, 263002 (2017).
  * Maaß  _et al._ (2019) B. Maaß, T. Hüther, K. König, J. Krämer, J. Krause, A. Lovato, P. Müller, K. Pachucki, M. Puchalski, R. Roth, R. Sánchez, F. Sommer, R. B. Wiringa, and W. Nörtershäuser, Phys. Rev. Lett. 122, 182501 (2019).
  * Maaß  _et al._ (2017) B. Maaß, P. Müller, W. Nörtershäuser, J. Clark, C. Gorges, S. Kaufmann, K. König, J. Krämer, A. Levand, R. Orford, R. Sánchez, G. Savard, and F. Sommer, Hyperfine Interactions 238, 25 (2017).
  * Lu  _et al._ (2013) Z.-T. Lu, P. Mueller, G. W. F. Drake, W. Nörtershäuser, S. C. Pieper, and Z.-C. Yan, Rev. Mod. Phys. 85, 1383 (2013).


# Stability and Bifurcation Analysis of a Fractional Order Delay Differential
Equation Involving Cubic Nonlinearity

Sachin Bhalekar111Corresponding Author Email<EMAIL_ADDRESS>Deepa Gupta  School of Mathematics and Statistics, University of Hyderabad,
Hyderabad, 500046 India

###### Abstract

Fractional derivative and delay are important tools in modeling memory
properties in the natural system. This work deals with the stability analysis
of a fractional order delay differential equation

$D^{\alpha}x(t)=\delta x(t-\tau)-\epsilon x(t-\tau)^{3}-px(t)^{2}+qx(t).$

We provide linearization of this system in a neighbourhood of equilibrium
points and propose linearized stability conditions. To discuss the stability
of equilibrium points, we propose various conditions on the parameters
$\delta$, $\epsilon$, $p$, $q$ and $\tau$. Even though there are five
parameters involved in the system, we are able to provide the stable region
sketch in the $q\delta-$plane for any positive $\epsilon$ and $p$. This
provides the complete analysis of stability of the system. Further, we
investigate chaos in the proposed model. This system exhibits chaos for a wide
range of delay parameter.

##  1 Introduction

The non-local operators and the delay are the crucial tools in modeling memory
properties in the natural system [1, 2, 3]. The non-local operator viz.
fractional order derivative(FD) is widely analyzed and applied by Scientists
and Engineers [4, 5]. The flexible order (integer, real, complex numbers as
well as functions) is yet another reason to employ the FD in the systems which
show an intermediate behaviour e.g. viscoelasticity [6].
Fractional calculus(FC) is used to model diffusion by Mainardi [7, 8], Wyss
[9], Luchko [10], Daftardar- Gejji and coworkers [11, 12].
Magin [13, 14] presented ample number of applications of FC in bioengineering
and related areas.
The FD is proved useful in designing robust controllers and other engineering
applications [5, 15, 16, 17, 18].
Existence and uniqueness of solutions of fractional differential equations is
discussed in [19, 20, 21].
Matignon proposed stability results of FDEs in his seminal work [22].
Since the past values of state are also included in the model, the delay
differential equations become an infinite dimensional dynamical systems [23,
24].
The delay models are observed in various phenomena [23, 25, 26].
The fractional order delay differential equations(FDDE) contain FD as well as
the delay. The stability analysis of FDDE is presented by Bhalekar in [27, 28,
29].
Exact and discritized stability of linear FDDEs is discussed in [30] by Kaslik
and Sivasundaram.
Stabilization problem of neutral FDDEs is given in [31]. Various numerical
schemes [32, 33, 34, 35] are designed by the researchers.
Some issues related with the initialization of FDDEs are examined in [36].
In this work, we propose the stability results of a FDDE involving a cubic
nonlinearity. The Section 2 deals with the preliminaries. Stability results
are proposed in Section 3. Section 4 provide stable region for an equilibrium
point. We analyze the chaos in proposed system in the Section 5. Finally,
conclusions are given in Section 6.

##  2 Preliminaries

In this section, we present definitions available in the literature[23, 24,
37, 38, 39].

######  Definition 1 (Fractional Integral).

For any $f\in\mathcal{L}^{1}(0,b)$ the Riemann-Liouville fractional integral
of order $\upmu>0$, is given by

$\textit{I}^{\upmu}f(t)=\dfrac{1}{\Gamma(\upmu)}\int_{0}^{t}(t-\tau)^{\upmu-1}f(\tau)d\tau,\quad
0<t<b.$

######  Definition 2 (Caputo Fractional Derivative).

For $f\in\mathcal{L}^{1}(0,b)$, $0<t<b$ and $m-1<\upmu\leq m$,
$m\in\mathbb{N}$, the Caputo fractional derivative of function $f$ of order
$\upmu$ is defined by,

$\textit{D}^{\upmu}f(t)=\begin{cases}\frac{d^{m}}{dt^{m}}f(t),\textit{ if
}\quad\upmu=m\\\ \textit{I}^{m-\upmu}\dfrac{d^{m}f(t)}{dt^{m}},\textit{ if
}\quad m-1<\upmu<m.\end{cases}$

Note that for $m-1<\upmu\leq m$, $m\in\mathbb{N},$

$\textit{I}^{\upmu}\textit{D}^{\upmu}f(t)=f(t)-\sum_{k=0}^{m-1}\dfrac{d^{k}f(0)}{dt^{k}}\dfrac{t^{k}}{k!}.$

######  Definition 3 (Equilibrium Point).

Consider the generalized delay differential equation

$D^{\alpha}{x}(t)=f(x(t),x(t-\tau)),\quad 0<\alpha\leq 1,$ (1)

where $\tau\textgreater 0$, $f:E\rightarrow\mathbb{R}$,
$E\subseteq\mathbb{R}^{2}$ is open and $f\in C^{1}(E)$.

A steady state solution of equation (1) is called an equilibrium point.

Note that $x^{*}$ is an equilibrium point if and only if

$f(x^{*},x^{*})=0.$ (2)

Consider the initial-value problem for the nonautonoumous delay differential
equation (1) with the initial data

$x(t)=\phi(t),-\tau\leq t\leq 0,\textit{ where
}\phi:[-\tau,0]\rightarrow\mathbb{R}.$ (3)

######  Notation 1.

The solution of delay differential equation (1) with initial data (3) is
denoted by $x(t,\phi).$

The norm of $\phi$ is given by

$||\phi||=\sup_{-\tau\leq t\leq 0}|\phi(t)|.$

######  Definition 4.

An equilibrium point $x^{*}$ of equation (1) is stable if for any given
$\epsilon\textgreater 0$, there exist $\delta\textgreater 0$ such that
$||\phi-x^{*}||\textless\delta\Rightarrow|x(t,\phi)-x_{*}|\textless\epsilon,\quad
t\geq 0.$

######  Definition 5.

An equilibrium point $x^{*}$ is asymptotically stable if it is stable and
there exists $b_{0}\textgreater 0$ such that $||\phi-x^{*}||\textless
b_{0}\Rightarrow\lim_{t\longrightarrow\infty}x(t,\phi)=x^{*}$.

######  Definition 6.

Equilibrium point which is not stable is called unstable.

######  Theorem 1.

[27] Suppose $x^{*}$ is an equilibrium solution of the fractional order delay
differential equation

$D^{\alpha}x(t)=ax(t)+bx(t-\tau).$

Case 1 If b $\in(-\infty,-|a|)$ then the stability region of $x^{*}$ in
$(\tau,a,b)$ parameter space is located between the plane $\tau=0$ and

$\tau_{*}=\dfrac{\arccos\Bigg{(}\dfrac{\Bigg{(}acos\Big{(}\dfrac{\alpha\pi}{2}\Big{)}+\sqrt{b^{2}-a^{2}\sin^{2}\Big{(}\dfrac{\alpha\pi}{2}}\Big{)}\Bigg{)}\cos\dfrac{\alpha\pi}{2}-a}{b}\Bigg{)}}{\Bigg{(}a\cos\Big{(}\dfrac{\alpha\pi}{2}\Big{)}\pm\sqrt{b^{2}-a^{2}\sin^{2}\Big{(}\dfrac{\alpha\pi}{2}\Big{)}}\Bigg{)}^{1/\alpha}}.$
(4)

The equation undergoes Hopf bifurcation at this value.
Case 2 If $b\in(-a,\infty)$ then $x^{*}$ is unstable for any $\tau\geq 0.$
Case 3 If $b\in(a,-a)$ and $a<0$ then $x^{*}$ is stable for any $\tau\geq 0.$

Note:  In Case 1, we say that $x^{*}$ is delay dependent stable.

###  2.1 Linearization near equilibrium [27]

Let $x(t)$ be a solution of the generalized fractional delay differential
equation (1) perturbed infinitesimally from the equilibrium solution. Let
$\xi(t)=x(t)-x^{*}$. Then by using first order Taylor’s approximation, we get
a linearized equation of (1) as

$\displaystyle D^{\alpha}\xi=D^{\alpha}x$ $\displaystyle\hskip
42.67912pt=f(x(t),x(t-\tau))$ $\displaystyle\hskip
76.82234pt=f(x^{*}+\xi(t),x^{*}+\xi(t-\tau))$ $\displaystyle\hskip
130.88284pt=f(x^{*},x^{*})+\partial_{1}f(x^{*},x^{*})\xi(t)+\partial_{2}f(x^{*},x^{*})\xi(t-\tau).$
$\therefore D^{\alpha}\xi=a\xi+b\xi(t-\tau),$ (5)

where $a=\partial_{1}f$, b=$\partial_{2}f$ are partial derivative of $f$ with
respect to the first and second variables evaluated at $(x^{*},x^{*})$,
respectively. Equation (5) is local linearization of equation (1) near
$x^{*}$. The trajectories of the generalized fractional order delay
differential equation (1) in the neighbourhood of an equilibrium point have
the same form as the trajectories of equation (5) [24, 38].

##  3 Main Results

We propose the following model

$D^{\alpha}x(t)=\delta x(t-\tau)-\epsilon(x(t-\tau))^{3}-p(x(t))^{2}+qx(t),$
(6)

where $\delta$, $\epsilon$, $p$ and $q$ are all real numbers.
The non-linear terms in this equation are $(x(t-\tau))^{3}$ and $(x(t))^{2}.$
If $p=0$ and $q=0$ then this is an Ucar system [40]. Such models occur in many
physical models [41, 42].
In this case $f(x(t),x(t-\tau))=\delta
x(t-\tau)-\epsilon(x(t-\tau))^{3}-p(x(t))^{2}+qx(t).$ Therefore $a=-2px+q$ and
$b=\delta-3\epsilon(x-\tau)^{2}$. The corresponding equilibrium points are
$x_{1}^{*}=0$ and
$x_{2,3}^{*}=\dfrac{-p\pm\sqrt{p^{2}+4\epsilon(\delta+q)}}{2\epsilon}.$ Note
that for the existence of equilibrium points $x_{2,3}^{*}$ we need
$p^{2}+4\epsilon(\delta+q)\geq 0.$

###  3.1 Stability and bifurcation analysis of equilibrium point $x_{1}^{*}$.

For $x_{1}^{*}=0$, we have $a=q$ and $b=\delta$.

######  Theorem 2.

If $\delta+q>0$, then the equilibrium point $x_{1}^{*}$ is unstable for all
$\tau\geq 0.$

###### Proof.

For $\delta+q>0$, we have $\delta>-q$.
Therefore, $\delta\in(-q,\infty)$ which implies that $b\in(-a,\infty)$.
So, by Theorem (1) Case (2), $x_{1}^{*}$ is unstable for all $\tau\geq 0$
which completes the proof. ∎

The illustration of this Theorem is given in Figure (1) by setting $\delta=2$,
$\epsilon=1$, $p=1$, $\tau=0.5$ and $q=1$.

######  Theorem 3.

If $\delta+q<0$ and $\delta\geqslant q$, then $x_{1}^{*}$ is asymptotically
stable $\forall\tau\geqslant 0$.

###### Proof.

The conditions $\delta+q<0$ and $\delta\geq q$ hold only when $q$ is negative
and $\delta\in(q,-q)$. Therefore by Theorem (1) Case (3), the equilibrium
point $x_{1}^{*}$ is stable for all $\tau\geq 0.$ This completes the proof. ∎

This result is verified in Figure (2) by putting $\delta=2$, $\epsilon=1$,
$p=1$, $q=-3$ and $\tau=0.8.$

######  Theorem 4.

If $\delta+q<0$, $\delta<q$, then there exists

$\displaystyle\tau_{*}=\dfrac{\arccos\Bigg{(}\dfrac{\Bigg{(}qcos\Big{(}\dfrac{\alpha\pi}{2}\Big{)}+\sqrt{\delta^{2}-q^{2}\sin^{2}\Big{(}\dfrac{\alpha\pi}{2}}\Big{)}\Bigg{)}\cos\dfrac{\alpha\pi}{2}-q}{\delta}\Bigg{)}}{\Bigg{(}q\cos\Big{(}\dfrac{\alpha\pi}{2}\Big{)}\pm\sqrt{\delta^{2}-q^{2}\sin^{2}\Big{(}\dfrac{\alpha\pi}{2}\Big{)}}\Bigg{)}^{1/\alpha}}$
(7)

such that the equilibrium point $x_{1}^{*}$ is asymptotically stable for
$0<\tau<\tau_{*}$ and unstable for $\tau>\tau_{*}$.

###### Proof.

If $q+\delta<0$ and $\delta<q$ then $\delta$ is always negative. Therefore,
$\delta\in(-\infty,-|q|)$ and hence the condition in the Case (1) of Theorem
(1) is satisfied. Therefore there exists the critical value of delay
$\tau_{*}$ given by the equation (7) and $x_{1}^{*}$ is delay dependent
stable. ∎

To verify this result, we take $\delta=-3,\quad\epsilon=1,\quad p=1\quad
and\quad q=-2$ which satisfy the conditions given in the Theorem 4. In this
case, the critical value of delay is $\tau_{*}\approx 1.0690$. So, for
$\tau=0.8$, we get stable solution (cf. Figure (3)) whereas for $\tau=1.1$ we
get unstable solution (cf. Figure (4)).

Note: The stability of equilibrium point $x_{1}^{*}$ is independent of values
of the parameters $p$ and $q.$
We also summarise the Theorems 2, 3 and 4 in Figure (5).

Figure 1:  $x_{1}^{*}$ is unstable for $\tau=0.5$, $\delta=2$, $\epsilon=1$,
$p=1$ and $q=1$ Figure 2: For $\tau=0.8$, $\delta=2$, $\epsilon=1$, $p=1$ and
$q=-3$, $x_{1}^{*}$ is stable Figure 3:  $x_{1}^{*}$ is stable for $\tau=0.8$,
$\delta=-3$, $\epsilon=1$, $p=1$ and $q=-2$ Figure 4: Unstable equilibrium
point $x_{1}^{*}$ for $\tau=1.1$, $\delta=-3$, $\epsilon=1$, $p=1$ and $q=-2$
Figure 5: Stability regions for $x_{1}^{*}$ in $\delta$ and $q$ plane

###  3.2 Stability and bifurcation analysis of $x_{2}^{*}$

We have

$a=\dfrac{p^{2}}{\epsilon}-\dfrac{p\sqrt{p^{2}+4\epsilon(\delta+q)}}{\epsilon}+q$
(8)

and

$b=\delta-\dfrac{3[2p^{2}+4\epsilon(\delta+q)-2p\sqrt{p^{2}+4\epsilon(\delta+q)}]}{4\epsilon}$
(9)

at $x_{2}^{*}$. Hence, we get

$a+b=\dfrac{\sqrt{p^{2}+4\epsilon(\delta+q)}[p-\sqrt{p^{2}+4\epsilon(\delta+q)}]}{2\epsilon}.$
(10)

######  Theorem 5.

If $\epsilon>0$, $p>0$ and $0<-q<\delta<-2q$, then the equilibrium point
$x_{2}^{*}$ is asymptotically stable for all $\tau\geq 0$.

###### Proof.

Since, $\epsilon>0$ and $-q<\delta$ we have $\delta+q>0.$
Also, $4\epsilon(\delta+q)>0.$
$\Rightarrow p^{2}+4\epsilon(\delta+q)>p^{2}$
$\Rightarrow\sqrt{p^{2}+4\epsilon(\delta+q)}-p>0$
$\Rightarrow-\sqrt{p^{2}+4\epsilon(\delta+q)}+p<0$
$\Rightarrow\sqrt{p^{2}+4\epsilon(\delta+q)}[-\sqrt{p^{2}+4\epsilon(\delta+q)}+p]<0.$
So, from equation (10) and $\epsilon>0$ we have $a+b<0$.
Further, we get
$a-b=4q+2\delta+\dfrac{5p^{2}}{2\epsilon}-\dfrac{5p\sqrt{p^{2}+4\epsilon(\delta+q)}}{2\epsilon}.$
Using, $\delta<-2q$ and $\sqrt{p^{2}+4\epsilon(\delta+q)}>p$, we have $a-b<0.$
$\Rightarrow a<b.$
Therefore, by using Case (3) of Theorem (1) we conclude that the equilibrium
point $x_{2}^{*}$ is asymptotically stable for all $\tau\geq 0$. ∎

By choosing $\delta=3$, $\epsilon=1$, $p=1$ and $q=-2$ in the equation (6) we
get stable solution for all $\tau\geq 0$ by the Theorem (5). We verified this
result by taking $\tau=\left\\{0.5k|k=1,2,\ldots,20\right\\}$. Figure (6)
shows stable orbit at $\tau=0.5$.

Figure 6:  $x_{2}^{*}$ is stable for $\tau=0.5$, $\delta=3$, $\epsilon=1$,
$p=1$ and $q=-2$

######  Theorem 6.

If $\epsilon>0$, $p>0$, $\delta<\dfrac{-p^{2}}{32\epsilon}<0$ and
$(q+\delta)>0$ then there exists $\tau_{*}$ as given in equation (4) (with $a$
and $b$ are defined in equations (8) and (9) respectively) such that
$x_{2}^{*}$ is asymptotically stable for $0<\tau<\tau_{*}$ and unstable for
$\tau>\tau_{*}.$

###### Proof.

Step 1: Since $q+\delta>0$, $\epsilon>0$ and $p>0$ we have,
$p^{2}+4\epsilon(\delta+q)>p^{2}.$
$\Rightarrow p-\sqrt{p^{2}+4\epsilon(\delta+q)}<0.$
Therefore by equation (10), $(a+b)<0.$
Step 2: Since, $\delta<\dfrac{-p^{2}}{32\epsilon},\textit{ we have
}q+\delta<q-\dfrac{p^{2}}{32\epsilon}$
$\Rightarrow 4\epsilon(q+\delta)<4\epsilon(q-\dfrac{p^{2}}{32\epsilon})$
$\Rightarrow 0<p^{2}+4\epsilon(\delta+q)<p^{2}+4\epsilon q-\dfrac{p^{2}}{8}\\\
\Rightarrow\sqrt{p^{2}+4\epsilon(\delta+q)}<\sqrt{\dfrac{7p^{2}}{8}+4\epsilon
q}\\\
\Rightarrow\dfrac{-p}{\epsilon}\sqrt{p^{2}+4\epsilon(\delta+q)}>-\dfrac{p}{\epsilon}\sqrt{\dfrac{7p^{2}}{8}+4\epsilon
q}\\\ \Rightarrow
a>q+\dfrac{p^{2}}{\epsilon}-\dfrac{p}{\epsilon}\sqrt{\dfrac{7p^{2}}{8}+4\epsilon
q}\\\ \Rightarrow
a>-\delta+\dfrac{p^{2}}{\epsilon}-\dfrac{p}{\epsilon}\sqrt{\dfrac{7p^{2}}{8}+4\epsilon
q}$.
Since, $32\epsilon q>-32\epsilon\delta>p^{2}$
we have $a>-\delta+\dfrac{p^{2}}{\epsilon}-\dfrac{p^{2}}{\epsilon}$
$\Rightarrow a>-\delta$
$\Rightarrow a>0$.
Step 3: Further, $b$ can also be written as
$b=\delta-\dfrac{3\epsilon(a-q)^{2}}{4p^{2}}.$
$\Rightarrow b<\dfrac{-p^{2}}{32\epsilon}-\dfrac{3\epsilon(a-q)^{2}}{4p^{2}}$
which shows that $b<0.$
So, $a+b<0$ and $b<a.$
$\Rightarrow b\in(-\infty,-|a|).$
Therefore, by using Case (1) of Theorem (1), we get the required critical
value $\tau_{*}$ of delay.

∎

By setting $\delta=-1/2$, $p=4$, $\epsilon=2$ and $q=1$ we get
$\tau_{*}\approx 2.6521$. The convergent solution for $\tau=1.8$ is given in
Figure (8), whereas divergent solution for $\tau=2.9$ is given in Figure (9).

######  Theorem 7.

If $\epsilon>0$, $p>0$ and $(q+\delta)<0$ then $x_{2}^{*}$ is unstable for all
$\tau\geq 0.$

###### Proof.

Since $\epsilon>0$ and $(q+\delta)<0,$
we have $0<p^{2}+4\epsilon(\delta+q)<p^{2}.$
Therefore, $\sqrt{p^{2}+4\epsilon(\delta+q)}<p.$
So, from (10) we get $a+b>0.$
Hence, by Case (2) of Theorem (1), we get the required result. ∎

We verify the Theorem (7) by putting $\delta=-2$, $\epsilon=1$, $p=3$, $q=1$
and $\tau=0.5$ in equation (6) for which we get the unstable curve which is
shown in Figure (7).

Figure 7:  $x_{2}^{*}$ is unstable for $\tau=0.5$, $\delta=-2$, $\epsilon=1$,
$p=3$ and $q=1$ Figure 8: $x_{2}^{*}$ is asymptotically stable for
$\delta=-1/2$, $\epsilon=2$, $p=4$, $q=1$ and $\tau=1.8$ Figure 9: $x_{2}^{*}$
is unstable for $\delta=-1/2$, $\epsilon=2$, $p=4$, $q=1$ and $\tau=2.9$

We propose few more stability results in the Theorems 8, 9 and 10.

######  Theorem 8.

If $\epsilon<0$, $p>0$ and $q+\delta<0$ then $x_{2}^{*}$ is unstable for all
$\tau\geq 0.$

###### Proof.

For, $\epsilon<0$ and $q+\delta<0$ we have $4\epsilon(q+\delta)>0.$
So, $p^{2}+4\epsilon(\delta+q)>p^{2}>0.$
$\Rightarrow p-\sqrt{p^{2}+4\epsilon(\delta+q)}<0$
$\Rightarrow\sqrt{p^{2}+4\epsilon(\delta+q)}[p-\sqrt{p^{2}+4\epsilon(\delta+q)}]<0$
Therefore, by equation (10), we conclude that $a+b>0$ because $\epsilon<0$.
Hence $b\in(-a,\infty).$
Therefore, we get the required result by employing Case (2) of Theorem (1). ∎

We verified Theorem (8) by setting $\delta=1$, $\epsilon=-1$, $p=1$, $q=-2$
and various values of $\tau\in(0,100)$. Figure (10) shows unbounded solution
in this case, with $\tau=0.3.$

Figure 10: For $\delta=1$, $\epsilon=-1$, $p=1$, $q=-2$ and $\tau=0.3$ we get
unstable solution for $x_{2}^{*}$

######  Theorem 9.

If $\epsilon>0$, $p<0$, $\delta\leq\dfrac{3p^{2}}{4\epsilon}$ and
$q>\dfrac{-p^{2}-4\delta\epsilon}{4\epsilon}$ then $\exists$ $\tau_{*}$ as
given in (4), where $a$ and $b$ are defined by (8) and (9) respectively, such
that $x_{2}^{*}$ is asymptotically stable for $0<\tau<\tau_{*}$ and unstable
for $\tau>\tau_{*}.$

###### Proof.

We can conclude that $a+b<0$ by using the equation (10) and the conditions
$\epsilon>0$ and $p<0$. Therefore, to prove this Theorem we need either $b<a$
or $a>0.$
However,
$a>\dfrac{-p^{2}-4\delta\epsilon}{4\epsilon}+\dfrac{p^{2}}{\epsilon}-\dfrac{p\sqrt{p^{2}+4\epsilon(\delta+q)}}{\epsilon}$
because $q>\dfrac{-p^{2}-4\delta\epsilon}{4\epsilon}$.
Also, we have $\delta\leq\dfrac{3p^{2}}{4\epsilon}$.
$\Rightarrow-\delta\geq\dfrac{-3p^{2}}{4\epsilon}.$
Using these inequalities we have,
$a>\dfrac{-p\sqrt{p^{2}+4\epsilon(\delta+q)}}{\epsilon}.$
Further, $p<0$ and $\epsilon>0$ so, we have $a>0.$
Hence, $b\in(-\infty,-|a|)$. Therefore, the proof follows by using the Case
(1) of Theorem (1). ∎

Illustration of this theorem is given in Figure (11) by setting the parameters
as $\delta=\dfrac{3}{8}$, $\epsilon=2$, $p=-1$ and $q=1$ in equation (6). For
this set of parameters, we get $x_{2}^{*}=1.11603$ and $\tau_{*}=0.157185.$
Hence $\tau=0.08$ gives stable solution (cf. Figure (11)) and $\tau=0.36$
gives unstable solution(cf. Figure (12)).

Figure 11:  $x_{2}^{*}$ is stable for $\tau=0.08,\delta=3/8,\epsilon=2$,$p=-1$
and $q=1$ Figure 12:  $x_{2}^{*}$ is unstable for
$\tau=0.36,\delta=3/8,\epsilon=2$, $p=-1$ and $q=1$

######  Theorem 10.

If $\epsilon<0$, $p<0$ and $\delta+q<0$ then $x_{2}^{*}$ is unstable for all
$\tau\geq 0$.

###### Proof.

We can write $a+b$ as
$\dfrac{-[4\epsilon(\delta+q)+p^{2}]+p\sqrt{p^{2}+4\epsilon(\delta+q)}}{2\epsilon}.$
Since $p<0$ and $\epsilon$ is also negative quantity so, $a+b>0$. Therefore by
Case (2) Theorem (1) we have the equilibrium point $x_{2}^{*}$ is unstable for
$\tau\geq 0$. ∎

Verification of this Theorem (10) is given in Figure (13) by setting
$\delta=-2$, $p=-1$, $q=1$, $\epsilon=-1$ and $\tau=0.2$.

Figure 13:  Theorem (10) is verified for the parameter values $\delta=-2$,
$\epsilon=-1$, p=-1 and q=1 and $\tau=0.2$

##  4 Stable region for $x_{2}^{*}$

In this section, we sketch the stable region for $x_{2}^{*}$ with $\epsilon>0$
and $p>0.$
From Theorem (1), it is clear that the curves $a-b=0$ (with $a\leq 0$) and
$a+b=0$ are bifurcation curves.
In this case, $a+b=0$ gives $\delta=-q.$
Further, $a-b=0$ with $a\leq 0$ gives, $\delta=g_{1}(p,q,\epsilon)$ and
$\delta=g_{2}(p,q,\epsilon),$ where
$g_{1}(p,q,\epsilon)=\dfrac{15p^{2}-16q\epsilon+5\sqrt{9p^{4}-16p^{2}q\epsilon}}{8\epsilon},$
$-\infty<q\leq q_{2}$,
$g_{2}(p,q,\epsilon)=\dfrac{15p^{2}-16q\epsilon-5\sqrt{9p^{4}-16p^{2}q\epsilon}}{8\epsilon},$
$q_{3}\leq q\leq q_{2}$,
$q_{2}=\dfrac{9p^{2}}{16\epsilon}$ and $q_{3}=\dfrac{-p^{2}}{\epsilon}.$
The branch $\delta=g_{2}(p,q,\epsilon)$ will be valid for
$q\in\Big{[}q_{3},q_{2}\Big{]}.$ If $q<q_{3}$ then either $a$ and $b$ become
complex numbers or $a\neq b$, along the curve $\delta=g_{2}(p,q,\epsilon)$. If
$q>q_{2}$ then $g_{2}(p,q,\epsilon)\notin\mathbb{R}.$
We have following observations:

  * (1)

$g_{1}(p,q,\epsilon)$ is decreasing in the interval $-\infty<q\leq q_{2}$
because

$\dfrac{\partial g_{1}}{\partial
q}=-2-5\dfrac{p}{\sqrt{9p^{2}-16q\epsilon}}<0$

,
for $p>0$ and $\epsilon>0$.

  * (2)

$\lim_{q\to q_{2}}g_{1}(p,q,\epsilon)=\lim_{q\to
q_{2}}g_{2}(p,q,\epsilon)=\dfrac{3p^{2}}{4\epsilon}.$ Furthermore,
$\lim_{q\to-\infty}g_{1}(p,q,\epsilon)=\infty.$

  * (3)

Nature of $g_{2}(p,q,\epsilon):$ •$g_{2}(p,q,\epsilon)$ is monotonically
decreasing in $q_{3}\leq q\leq q_{0}$ where
$q_{0}=\dfrac{11p^{2}}{64\epsilon}:$
We have $\dfrac{\partial g_{1}}{\partial
q}=-2-5\dfrac{p}{\sqrt{9p^{2}-16q\epsilon}}$ and $q_{3}<q<q_{0}.$
$\Rightarrow\dfrac{25p^{2}}{4}<9p^{4}-16q\epsilon<25p^{2}$
$\Rightarrow\dfrac{1}{5p}<\dfrac{1}{\sqrt{9p^{2}-16q\epsilon}}<\dfrac{2}{5p}$
$\Rightarrow 1<\dfrac{5p}{\sqrt{9p^{2}-16q\epsilon}}<2$
$\Rightarrow\dfrac{\partial g_{1}}{\partial q}<0$ $\forall
q\in\Big{[}q_{3},q_{0}\Big{]}.$
•$g_{2}(p,q,\epsilon)$ is monotonic increasing in $q_{0}\leq q\leq q_{2}:$
If $q\in\Big{[}q_{0},q_{2}\Big{]}$, then
$\dfrac{11p^{2}}{4}<16q\epsilon<9p^{2}.$
$\Rightarrow 0<9p^{2}-16q\epsilon<\dfrac{25p^{2}}{\epsilon}$
$\Rightarrow 2<\dfrac{5p}{\sqrt{9p^{2}-16q\epsilon}}.$
So, $\dfrac{\partial g_{2}}{\partial q}>0$ $\forall
q\in\Big{[}q_{0},q_{2}\Big{]}$.
•Local minima of $g_{2}(p,q,\epsilon)$ is at $q=q_{0}$ with minimum value
$\delta_{1}=\dfrac{-p^{2}}{32\epsilon}.$
Since $|q_{0}|>|\delta_{1}|$, $g_{2}(p,q,\epsilon)$ lies above the curve
$\delta=-q$ for $0\leq q\leq q_{2}.$

  * (4)

Intersection of $\delta=g_{2}(p,q,\epsilon)$ and $\delta=-q$ is $(0,0)$.
Further, $\delta=g_{2}(p,q,\epsilon)$ intersects $q$ axis at $(0,0)$ and at
$(q_{1},0)$ where $q_{1}=\dfrac{5p^{2}}{16\epsilon}$.
$\lim_{q\to
q_{3}}g_{2}(p,q,\epsilon)=\dfrac{3p^{2}}{4\epsilon}<\dfrac{p^{2}}{\epsilon}.$
So, in the interval $q_{3}\leq q\leq 0,$ the curve
$\delta=g_{2}(p,q,\epsilon)$ lies below the curve $\delta=-q.$

  * (5)

The curve $\delta=g_{1}(p,q,\epsilon)$ is always above the curve $\delta=-q:$
Suppose $q<0$.
Since,
$g_{1}(p,q,\epsilon)=-2q+\dfrac{15p^{2}+5p\sqrt{9p^{2}-16q\epsilon}}{8\epsilon}$
and $p$ and $\epsilon$ both are positive, $g_{1}(p,q,\epsilon)>-2q>-q.$
Further for any $q$,
$g_{1}(p,q,\epsilon)$ =
$\dfrac{6p^{2}+\sqrt{9p^{2}-16q\epsilon}(\sqrt{9p^{2}-16q\epsilon}+5p)}{8\epsilon}>0$
because $p$ and $\epsilon$ both are positive.
Therefore, $g_{1}(p,q,\epsilon)>-q$, $q\in\Big{(}-\infty,q_{2}\Big{)}.$

  * (6)

Intersection of $\delta=g_{1}(p,q,\epsilon)$ with $\delta$ axis is
$\delta_{0}=g_{1}(p,0,\epsilon)=\dfrac{30p^{2}}{8\epsilon}>0.$
Using these observations, we sketch the stability regions for $x_{2}^{*}$ in
Figure (14).
We have-

Figure 14:  Stability region for $x_{2}^{*}$ in $q\delta-$plane

  * (A)

If $q>q_{1}$ and $\delta>-q$, then $x_{2}^{*}$ is delay dependent stable.

  * (B)

If $-q<\delta<g_{1}(p,q,\epsilon)$ and $q<0$ then $x_{2}^{*}$ is
asymptotically stable $\forall\tau\geqslant 0.$

  * (C)

If $0\leq q\leq q_{1}$ and-

    * (i)

$\delta\in\Big{(}-q,g_{2}(p,q,\epsilon)\Big{)}\cup\Big{(}g_{1}(p,q,\epsilon),\infty\Big{)}$
then $\exists$ $\tau_{*}$ as given in (4) and $x_{2}^{*}$ is delay dependent
stable.

    * (ii)

$g_{2}(p,q,\epsilon)<\delta<g_{1}(p,q,\epsilon)$ then $x_{2}^{*}$ is
asymptotically stable, $\forall\tau\geqslant 0$.

##  5 Chaos

We observed chaotic oscillations in system (6) for some parameter values. The
Figure (15) shows the bifurcation diagram for the parameter set $\delta=5$,
$\epsilon=2$, $p=0.01$, $q=-2$ and $\alpha=0.95$. The horizontal axis is the
delay $\tau$. Equilibrium points in this case are $x_{1}^{*}=0$,
$x_{2}^{*}=1.22725$ and $x_{3}^{*}=-1.22725.$
The equilibrium point $x_{1}^{*}$ is unstable for all $\tau\geq 0$. For
$x_{2}^{*}$, $b<a<0$. Therefore, there exists $\tau_{*}=0.6427665$ such that
$x_{2}^{*}$ is asymptotically stable for $0<\tau<\tau_{*}$.
Similarly, for $x_{3}^{*}$, $b<a<0$ and $\tau_{*}=0.620447$. Thus, the system
is unstable for $\tau>0.642766.$
We observed periodic limit cycles for $0.65\leq\tau\leq 2.2$.
Figures (16) and (17) show periodic limit cycles for $\tau=1.4$ and
$\tau=1.95$ respectively.
Chaos is observed for $\tau>2.2$. Figures (18) and (19) show chaotic
attractors for $\tau=2.3$ and $\tau=2.5$ respectively. The chaos is confirmed
with bifurcation diagram (cf. Figure (15)) and the positive values of maximum
Lyapunov exponents (Table 1).

$\tau$ |  Maximum Lyapunov Exponents |  Behaviour of System
---|---|---
0.6 |  -0.912265 |  Stable
1.6 |  -0.002104 |  Limit cycle
1.8 |  -0.000428 |  Limit cycle
2.3 |  0.546279 |  Chaotic oscillations
2.5 |  1.083852 |  Chaotic oscillations
Table 1:  Maximum Lyapunov Exponents

We used the algorithm described by Kodba et al [43] which is based on the time
series analysis techniques and the work by Wolf et al [44].

Figure 15:  Bifurcation diagram Figure 16:  Graph of $x(t)$ versus $x(t-\tau)$
for $\tau=1.4$ Figure 17:  Periodic limit cycle for $\tau=1.95$ Figure 18:
Chaotic attractor for $\tau=2.3$ Figure 19:  Chaotic attractor for $\tau=2.5$

##  6 Conclusion

In this work, we considered a fractional order delay differential equation

$D^{\alpha}x(t)=\delta x(t-\tau)-\epsilon(x(t-\tau))^{3}-p(x(t))^{2}+qx(t).$

For some values of parameters, there are three equilibrium points viz.
$x_{1}^{*}$, $x_{2}^{*}$ and $x_{3}^{*}$. We provided explicit stability
conditions for equilibrium points $x_{1}^{*}$ and $x_{2}^{*}$. We proposed
delay-dependent as well as delay-independent stability conditions. The results
are verified by setting particular values to parameters. The key finding is to
sketch the stable regions in the $q\delta-$plane which are valid for any
$\epsilon>0$ and $p>0.$
It is observed that the system shows chaotic oscillations for some range of
parameters. We provided the bifurcation diagram and the values of maximum
Lyapunov exponents to confirm the chaos in this system.
The stability of $x_{3}^{*}$ can be done as a future work.

##  7 Acknowledgments

S. Bhalekar acknowledges the University of Hyderabad for Institute of
Eminence-Professional Development Fund (IoE-PDF) by MHRD (F11/9/2019-U3(A)).
D. Gupta thanks University Grants Commission for financial support
(Ref.No.:201610026200).

## References

  * [1] D. Baleanu, R. L. Magin, S. Bhalekar, V. Daftardar-Gejji, Chaos in the fractional order nonlinear bloch equation with delay, Communications in Nonlinear Science and Numerical Simulation 25 (1-3) (2015) 41–49.
  * [2] S. Bhalekar, V. Daftardar-Gejji, Fractional ordered liu system with time-delay, Communications in Nonlinear Science and Numerical Simulation 15 (8) (2010) 2178–2191.
  * [3] S. Bhalekar, V. Daftardar-Gejji, D. Baleanu, R. Magin, Fractional bloch equation with delay, Computers and Mathematics with Applications 61 (5) (2011) 1355–1365.
  * [4] L. Debnath, Recent applications of fractional calculus to science and engineering, International Journal of Mathematics and Mathematical Sciences 2003 (54) (2003) 3413–3442.
  * [5] J. Lai, S. Mao, J. Qiu, H. Fan, Q. Zhang, Z. Hu, J. Chen, Investigation progresses and applications of fractional derivative model in geotechnical engineering, Mathematical Problems in Engineering 2016 (2016).
  * [6] P. J. Torvik, R. L. Bagley, On the appearance of the fractional derivative in the behavior of real materials, Journal of Applied Mechanics 51 (2) (1984) 298–298.
  * [7] F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena, Chaos, Solitons & Fractals 7 (9) (1996) 1461–1477.
  * [8] F. Mainardi, The fundamental solutions for the fractional diffusion-wave equation, Applied Mathematics Letters 9 (6) (1996) 23–28.
  * [9] W. Wyss, The fractional diffusion equation, Journal of Mathematical Physics 27 (11) (1986) 2782–2785.
  * [10] Y. Luchko, Maximum principle for the generalized time-fractional diffusion equation, Journal of Mathematical Analysis and Applications 351 (1) (2009) 218–223.
  * [11] H. Jafari, V. Daftardar-Gejji, Solving linear and nonlinear fractional diffusion and wave equations by adomian decomposition, Applied Mathematics and Computation 180 (2) (2006) 488–497.
  * [12] V. Daftardar-Gejji, S. Bhalekar, Solving fractional diffusion-wave equations using a new iterative method, Fractional Calculus and Applied Analysis 11 (2) (2008) 193–202.
  * [13] R. Magin, Fractional calculus in bioengineering, part 1, Critical Reviews in Biomedical Engineering 32 (1) (2004).
  * [14] R. L. Magin, Fractional calculus models of complex dynamics in biological tissues, Computers & Mathematics with Applications 59 (5) (2010) 1586–1593.
  * [15] C. A. Monje, Y. Chen, B. M. Vinagre, D. Xue, V. Feliu-Batlle, Fractional-order systems and controls: fundamentals and applications, Springer Science & Business Media, 2010.
  * [16] Y. Chen, I. Petras, D. Xue, Fractional order control-a tutorial, in: 2009 American control conference, IEEE, 2009, pp. 1397–1411.
  * [17] Q. Yang, D. Chen, T. Zhao, Y. Chen, Fractional calculus in image processing: a review, Fractional Calculus and Applied Analysis 19 (5) (2016) 1222–1249.
  * [18] W. Chen, H. Sun, X. Li, et al., Fractional derivative modeling in mechanics and engineering, Springer, 2022.
  * [19] D. Delbosco, L. Rodino, Existence and uniqueness for a nonlinear fractional differential equation, Journal of Mathematical Analysis and Applications 204 (2) (1996) 609–625.
  * [20] Y. Zhou, Existence and uniqueness of solutions for a system of fractional differential equations, Fractional Calculus and Applied Analysis 12 (2) (2009) 195–204.
  * [21] V. Daftardar-Gejji, A. Babakhani, Analysis of a system of fractional differential equations, Journal of Mathematical Analysis and Applications 293 (2) (2004) 511–522.
  * [22] D. Matignon, Stability results for fractional differential equations with applications to control processing, in: Computational engineering in systems applications, Vol. 2, Citeseer, 1996, pp. 963–968.
  * [23] H. L. Smith, An introduction to delay differential equations with applications to the life sciences, Vol. 57, Springer New York, 2011.
  * [24] M. Lakshmanan, D. V. Senthilkumar, Dynamics of nonlinear time-delay systems, Springer Science & Business Media, 2011.
  * [25] A. Namajūnas, K. Pyragas, A. Tamaševičius, Stabilization of an unstable steady state in a mackey-glass system, Physics Letters A 204 (3-4) (1995) 255–262.
  * [26] G. A. Bocharov, F. A. Rihan, Numerical modelling in biosciences using delay differential equations, Journal of Computational and Applied Mathematics 125 (1-2) (2000) 183–199.
  * [27] S. Bhalekar, Stability and bifurcation analysis of a generalized scalar delay differential equation, Chaos: An Interdisciplinary Journal of Nonlinear Science 26 (8) (2016) 084306.
  * [28] V. Daftardar-Gejji, Y. Sukale, S. Bhalekar, A new predictor–corrector method for fractional differential equations, Applied Mathematics and Computation 244 (2014) 158–182.
  * [29] S. B. Bhalekar, Stability analysis of a class of fractional delay differential equations, Pramana 81 (2) (2013) 215–224.
  * [30] E. Kaslik, S. Sivasundaram, Analytical and numerical methods for the stability analysis of linear fractional delay differential equations, Journal of Computational and Applied Mathematics 236 (16) (2012) 4027–4041.
  * [31] C. Bonnet, J. R. Partington, Stabilization of some fractional delay systems of neutral type, Automatica 43 (12) (2007) 2047–2053.
  * [32] V. Daftardar-Gejji, Y. Sukale, S. Bhalekar, Solving fractional delay differential equations: a new approach, Fractional Calculus and Applied Analysis 18 (2) (2015) 400–418.
  * [33] S. Bhalekar, V. Daftardar-Gejji, A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order, Journal of Fractional Calculus and Applications 1 (5) (2011) 1–9.
  * [34] L. Shi, Z. Chen, X. Ding, Q. Ma, A new stable collocation method for solving a class of nonlinear fractional delay differential equations, Numerical Algorithms 85 (4) (2020) 1123–1153.
  * [35] B. Yuttanan, M. Razzaghi, T. N. Vo, Legendre wavelet method for fractional delay differential equations, Applied Numerical Mathematics 168 (2021) 127–142.
  * [36] R. Garrappa, E. Kaslik, On initial conditions for fractional delay differential equations, Communications in Nonlinear Science and Numerical Simulation 90 (2020) 105359.
  * [37] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Elsevier, 1998.
  * [38] Z. Vukic, Nonlinear control systems, CRC Press, 2003.
  * [39] K. Diethelm, N. J. Ford, Analysis of fractional differential equations, Journal of Mathematical Analysis and Applications 265 (2) (2002) 229–248.
  * [40] A. Uçar, A prototype model for chaos studies, International Journal of Engineering Science 40 (3) (2002) 251–258.
  * [41] S. Bhalekar, Dynamical analysis of fractional order uçar prototype delayed system, Signal, Image and Video Processing 6 (3) (2012) 513–519.
  * [42] S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering, CRC Press, 2018.
  * [43] S. Kodba, M. Perc, M. Marhl, Detecting chaos from a time series, European Journal of Physics 26 (1) (2004) 205.
  * [44] A. Wolf, J. B. Swift, H. L. Swinney, J. A. Vastano, Determining lyapunov exponents from a time series, Physica D: Nonlinear Phenomena 16 (3) (1985) 285–317.
  * [45] Y. Luchko, R. Gorenflo, An operational method for solving fractional differential equations with the caputo derivatives, Acta Math. Vietnam 24 (2) (1999) 207–233.
  * [46] R. Magin, M. Ovadia, Modeling the cardiac tissue electrode interface using fractional calculus, Journal of Vibration and Control 14 (9-10) (2008) 1431–1442.
  * [47] J. Čermák, L. Nechvátal, On exact and discretized stability of a linear fractional delay differential equation, Applied Mathematics Letters 105 (2020) 106296.


# Reflexive polytopes and discrete polymatroids

Jürgen Herzog and Takayuki Hibi  Jürgen Herzog, Fachbereich Mathematik,
Universität Duisburg-Essen, Campus Essen, 45117 Essen, Germany
<EMAIL_ADDRESS>Takayuki Hibi, Department of Pure and Applied
Mathematics, Graduate School of Information Science and Technology, Osaka
University, Suita, Osaka 565-0871, Japan<EMAIL_ADDRESS>
###### Abstract.

A classification of discrete polymatroids whose independence polytopes are
reflexive will be presented.

###### 2010 Mathematics Subject Classification:

Primary 52B20; Secondary 05E40

The second author was supported by JSPS KAKENHI 19H00637.

## Introduction

The discrete polymatroid is introduced in [1]. In the present paper, as a
supplement to [1], a classification of discrete polymatroids whose
independence polytopes are reflexive will be presented. We refer the reader to
[1] and [2] for fundamental materials on discrete polymatroids.

##  1\. Reflexive polytopes

A convex polytope ${\mathcal{P}}\subset{\mathbb{R}}^{d}$ of dimension $d$ is
called a lattice polytope if each of its vertices belongs to
${\mathbb{Z}}^{d}$. A reflexive polytope is a lattice polytope
${\mathcal{P}}\subset{\mathbb{R}}^{d}$ of dimension $d$ for which the origin
of ${\mathbb{R}}^{d}$ belongs to the interior of ${\mathcal{P}}$ and the dual
polytope
${\mathcal{P}}^{\vee}=\\{{\mathbf{x}}\in{\mathbb{R}}^{d}:\langle{\mathbf{x}},{\mathbf{y}}\rangle\leq
1,\forall{\mathbf{y}}\in{\mathcal{P}}\\}$ of ${\mathcal{P}}$ is a lattice
polytope, where $\langle{\mathbf{x}},{\mathbf{y}}\rangle$ stands for the
canonical inner product of ${\mathbb{R}}^{d}$. A lattice polytope which can be
a reflexive polytope by parallel shift is also called reflexive.

Let ${\mathbf{e}}_{1},\ldots,{\mathbf{e}}_{d}$ denote the canonical basis
vectors of ${\mathbb{R}}^{d}$. Let $P\subset{\mathbb{Z}}^{d}_{+}$ be a
discrete polymatroid [1, Definition 2.1] on the ground set $[d]$. In what
follows one assumes that each ${\mathbf{e}}_{i}$ belongs to $P$. Let
${\mathcal{P}}={\mathcal{P}}_{P}\subset{\mathbb{R}}^{d}$ denote the lattice
polytope which is the convex hull of $P$ in ${\mathbb{R}}^{d}$. We call
${\mathcal{P}}$ the independence polytope of $P$. One has
$\dim{\mathcal{P}}=d$. Let $\rho=\rho_{P}$ denote the ground set rank function
[1, pp. 243] of ${\mathcal{P}}$. It follows from [1, Theorem 7.3] that

######  Lemma 1.1.

The independence polytope ${\mathcal{P}}$ is reflexive if and only if, for
each subset $X\subset[d]$ which is $\rho$-closed and $\rho$-inseparable [1,
pp. 257–258], one has $\rho(X)=|X|+1$.

A sublattice of $2^{[d]}$ is a collection ${\mathcal{L}}$ of subsets of $[d]$
with $\emptyset\in{\mathcal{L}}$ and $[d]\in{\mathcal{L}}$ such that, for all
$A$ and $B$ belonging to ${\mathcal{L}}$, one has $A\cap B\in{\mathcal{L}}$
and $A\cup B\in{\mathcal{L}}$.

######  Theorem 1.2.

(a) Let $P$ be a discrete polymatroid on the ground set $[d]$ and
$\rho=\rho_{P}$ the ground set rank function of ${\mathcal{P}}$. Let
${\mathcal{A}}$ be the set of $\rho$-closed and $\rho$-inseparable subsets of
${\mathcal{P}}$. If ${\mathcal{P}}$ is reflexive, then
${\mathcal{A}}\cup\\{\emptyset\\}$ is a sublattice of $2^{[d]}$.

(b) Conversely, given a sublattice ${\mathcal{L}}$ of $2^{[d]}$, there exists
a unique discrete polymatroid $P$ on the ground set $[d]$ for which
${\mathcal{L}}$ is the set of $\rho$-closed and $\rho$-inseparable subsets of
${\mathcal{P}}$ and ${\mathcal{P}}$ is reflexive.

###### Proof.

(a) If the independence polytope ${\mathcal{P}}$ of $P$ is reflexive, then
Lemma 1.1 says that $\rho(A)=|A|+1$ for each $A\in{\mathcal{A}}$. It follows
from [1, Proposition 7.2] that ${\mathcal{P}}$ consists of those
$(x_{1},\ldots,x_{d})\in{\mathbb{R}}^{d}$ for which

$x_{i}\geq 0,\,\,\,\,\,i=1,2,\ldots,d,$

and

(1) $\displaystyle\sum_{i\in A}x_{i}\leq|A|+1,\,\,\,\,\,A\in{\mathcal{A}}.$

Since each ${\mathbf{e}}_{i}\in{\mathcal{P}}$ and ${\mathcal{P}}$ is compact,
it follows that

$\bigcup_{A\in{\mathcal{A}}}A=[d].$

Furthermore, if $X\not\in{\mathcal{A}}$, then $\rho(X)>|X|+1$.

In fact, if $|X|$ = 1 and $X=\\{i\\}$, then
$3{\mathbf{e}}_{i}\in{\mathcal{P}}$ and $\rho(X)>2$. In general, if $|X|=q\geq
2$ and $X=\\{i_{1},\ldots,i_{q}\\}$ with $i_{1}<\cdots<i_{q}$, then, one has

(2)
$\displaystyle{\mathbf{v}}=\frac{q}{\,q-1\,}\sum_{j=1}^{q}{\mathbf{e}}_{i_{j}}\in{\mathcal{P}}$

and

$\rho(X)\geq|{\mathbf{v}}|=\frac{q^{2}}{\,q-1\,}>q+1.$

To see why (2) holds, one shows that ${\mathbf{v}}$ satisfies each of the
inequalities (1). Let $A\in{\mathcal{A}}$ with $X\subsetneq A$, then

$\frac{q^{2}}{\,q-1\,}\leq q+2=|X|+2\leq|A|+1.$

Let $A\in{\mathcal{A}}$ with $|X\cap A|=k<q$, then

$k\frac{q}{\,q-1\,}\leq k+1\leq|A|+1.$

One claims that ${\mathcal{A}}\cup\\{\emptyset\\}$ is a sublattice of
$2^{[d]}$. Let $A,B\in{\mathcal{A}}$ and suppose that either $A\cup
B\not\in{\mathcal{A}}\cup\\{\emptyset\\}$ or $A\cap
B\not\in{\mathcal{A}}\cup\\{\emptyset\\}$. Then

$\rho(A)+\rho(B)=|A|+|B|+2=|A\cup B|+|A\cap B|+2<\rho(A\cup B)+\rho(A\cup B),$

which contradict the fact that $\rho$ is submodular. Furthermore, since
$\bigcup_{A\in{\mathcal{A}}}A=[d]$, one has $[d]\in{\mathcal{A}}$, as desired.

(b) By virtue of [1, Theorem 9.1] one introduces the nondecreasing submodular
function $\rho:2^{[n]}\rightarrow{\mathbb{Z}}_{+}$ by setting

$\rho(X)=\min\\{|A|+1:X\subseteq
A,A\in{\mathcal{L}}\\},\,\,\,\,\,\emptyset\neq X\subset[d]$

together with $\rho(\emptyset)=0$. Let $P$ be the discrete polymatroid on the
ground set $[n]$ and $\rho$ the ground set rank function of ${\mathcal{P}}$.
Then ${\mathcal{L}}$ is the set of $\rho$-closed and $\rho$-inseparable
subsets of $[d]$. Furthermore, Lemma 1.1 guarantees that the independence
polytope ${\mathcal{P}}$ of $P$ is reflexive. On the other hand, suppose that
$P^{\prime}$ is a discrete polymatroid on the ground set $[d]$ and
$\rho^{\prime}$ its ground set rank function of the independence polytope
${\mathcal{P}}^{\prime}$ of $P^{\prime}$ for which ${\mathcal{L}}$ is the set
of $\rho^{\prime}$-closed and $\rho^{\prime}$-inseparable subsets of
$P^{\prime}$ and for which ${\mathcal{P}}^{\prime}$ is reflexive. Then by
using Lemma 1.1 again one has $\rho^{\prime}(A)=|A|+1$ for each
$A\in{\mathcal{L}}$. Hence ${\mathcal{P}}={\mathcal{P}}^{\prime}$ ([1,
Proposition 7.2]). Thus $P=P^{\prime}$ ([1, Theorem 3.4]), as desired.

##  2\. Examples

######  Example 2.1.

Let ${\mathcal{P}}\subset{\mathbb{R}}^{3}$ be the convex polytope whose facets
are each $x_{i}=0$ together with

$x_{1}+x_{2}=3,\,\,\,x_{2}+x_{3}=3,\,\,\,x_{1}+x_{2}+x_{3}=4.$

It can be checked that ${\mathcal{P}}$ is reflexive. However, ${\mathcal{P}}$
cannot be the independence polytope of a discrete polymatroid on the ground
set $[3]$. In fact, if ${\mathcal{P}}$ is the independence polytope of a
discrete polymatroid $P$ on the ground set $[3]$, then both $u=(0,3,0)$ and
$v=(1,2,1)$ belong to the set of bases [1, p. 245] of $P$. One has $|u|<|v|$,
which contradicts [1, Theorem 2.3].

######  Example 2.2.

Let ${\mathcal{L}}$ be a chain of length $d$ of $2^{[d]}$, say,

${\mathcal{L}}=\\{\emptyset,\\{d\\},\\{d-1,d\\},\ldots,\\{1,\ldots,d\\}\\}.$

Let $P$ denote the discrete polymatroid constructed in Theorem 1.2 (b). Let

${\mathcal{B}}=\\{[d],[d],[d-1],\ldots,[2],[1]\\}.$

Let $P^{\prime}$ denote the transversal polymatroid [1, p. 267] presented by
${\mathcal{B}}$. If $X\subset[d]$ and $i=\min(X)$, then it follows from the
proof of Theorem 1.2 (b) that

$\rho_{P}(X)=|\\{i,i+1,\ldots,d\\}|+1=(d-(i-1))+1=d+2-i$

On the other hand, by the definition of the ground set rank function of a
transversal polymatroid, one has

$\rho_{P^{\prime}}(X)=(d+1)-(i-1)=d+2-i.$

Hence $\rho_{P}=\rho_{P^{\prime}}$. Thus $P=P^{\prime}$.

It would be of interest for which sublattice ${\mathcal{L}}$ of $2^{[d]}$ the
discrete polymatroid constructed in Theorem 1.2 (b) can be a transversal
polymatroid.

## References

  * [1] J. Herzog and T. Hibi, Discrete polymatroids, J. Algebraic Combin. 16 (2020), 239–268.
  * [2] J. Herzog and T. Hibi, “Monomial Ideals,” GTM 260, Springer, 2011.