#  Accepted for publication in The Astrophysical Journal.
Infrared Emission of Specific Polycyclic Aromatic Hydrocarbon Molecules:
Cyanonaphthalenes
DRAFT:

Kaijun Li11affiliation: Department of Astronomy, Xiamen University, Xiamen,
Fujian 361005, China<EMAIL_ADDRESS>, Aigen Li22affiliation: Department of
Physics and Astronomy, University of Missouri, Columbia, MO 65211, USA;
<EMAIL_ADDRESS>, X.J. Yang33affiliation: Department of Physics, Xiangtan
University, 411105 Xiangtan, Hunan Province, China<EMAIL_ADDRESS>, and
Taotao Fang11affiliation: Department of Astronomy, Xiamen University, Xiamen,
Fujian 361005, China<EMAIL_ADDRESS>
###### Abstract

The unidentified infrared emission (UIE) features at 3.3, 6.2, 7.7, 8.6, 11.3
and 12.7$\,{\rm\mu m}$ are ubiquitously seen in a wide variety of
astrophysical regions and commonly attributed to polycyclic aromatic
hydrocarbon (PAH) molecules. However, the unambiguous identification of any
individual, specific PAH molecules has proven elusive until very recently two
isomers of cyanonapthalene, which consists of two fused benzene rings and
substitutes a nitrile (–-CN) group for a hydrogen atom, were discovered in the
Taurus Molecular Cloud based on their rotational transitions at radio
frequencies. To facilitate the James Webb Space Telescope (JWST) to search for
cyanonapthalenes in astrophysical regions, we model the vibrational excitation
of cyanonapthalenes and calculate their infrared emission spectra in a number
of representative astrophysical regions. The model emission spectra and
intensities will allow JWST to quantitatively determine or place an upper
limit on the abundances of cyanonapthalenes.

dust, extinction — ISM: lines and bands — ISM: molecules

##  1 Introduction

The “unidentified infrared emission” (UIE) bands, a distinct set of spectral
features at wavelengths of 3.3, 6.2, 7.7, 8.6, 11.3 and 12.7$\,{\rm\mu m}$,
dominate the infrared (IR) spectra of many bright astronomical objects. They
are ubiquitously seen in the interstellar medium (ISM) of our own galaxy and
star-forming galaxies, both near (e.g., the Local Group) and far (e.g.,
distant galaxies at redshifts $z>4$), and account for up to 20% of their total
infrared (IR) luminosity (see Li 2020). The identification of the UIE bands is
important as they are a useful probe of the cosmic star-formation history, and
their carriers are an essential player in galactic evolution. Although the
exact nature of the UIE bands remains unknown, a popular hypothesis is that
they arise from the C–H and C–C vibrational transitions of polycyclic aromatic
hydrocarbon (PAH) molecules (Léger & Puget 1984, Allamandola et al. 1985,
Tielens 2008).

Nevertheless, despite the popularity of the PAH hypothesis, unambiguous
astronomical identification of any individual, specific PAH molecules has long
proven elusive for nearly four decades ever since the PAH model was first
proposed in the 1980s. However, recent years have witnessed breakthroughs in
this aspect. Based on radio observations of the dark molecular cloud TMC-1
located within the Taurus Molecular Cloud made by the 100 m Green Bank
Telescope (GBT) in the frequency range of 8 to 34$\,{\rm GHz}$, McGuire et al.
(2021) reported the detection of the rotational transitions of
1-cyanonaphthalene (1-CNN) and 2-cyanonaphthalene (2-CNN), two isomers of
cyanonaphthalene (C10H7CN), in which a nitrile/cyano (–CN) group replaces one
of the hydrogen atoms of naphthalene (see Figure 1). This is the first
definitive identification of a specific PAH molecule in space. Similarly,
Burkhardt et al. (2021) conducted GBT radio observations of TMC-1 in the
frequency range of 2–12 and 18–34$\,{\rm GHz}$. They detected indene (C9H8), a
pure hydrocarbon PAH species composed of both a five- and six-membered ring
(see Figure 1).

Prior to these, benzene (C6H6) was detected in CRL 618, a protoplanetary
nebula, through a single weak absorption feature arising from its $\nu_{4}$
bending mode at $\sim\,$14.85$\,{\rm\mu m}$ (Cernicharo et al. 2001).
Recently, benzene has also been detected by the James Webb Space Telescope
(JWST) through the 14.85$\,{\rm\mu m}$ emission feature in a protoplanetary
disk around a low-mass star (Tabone et al. 2023). Benzonitrile (C6H5CN), a
single benzene ring with an attached CN group, was discovered also with GBT by
McGuire et al. (2018) in TMC-1. This discovery was also achieved by radio
observations of its rotational lines in the frequency range of 18 to 23$\,{\rm
GHz}$. By definition, benzene and benzonitrile are aromatic but not PAH
molecules. For illustration, Figure 1 lists all these specific aromatic
species spectroscopically identified in space.111Several small PAHs, including
naphthalene (C10H8), phenanthrene (C14H10), perylene (C16H10), pyrene
(C20H12), and their derivatives have been found in the Stardust samples
collected from comet Wild 2 (Sandford et al. 2006, Clemett et al. 2010), in
the Rosetta samples collected from comet 67P/Churyumov-Gerasimenko (Schuhmann
et al. 2019), in the Hayabusa 2 samples collected from the carbonaceous
asteroid Ryugu (Naraoka et al. 2023), and in interplanetary dust particles
possibly of cometary origin (Clemett et al. 1993). These PAH molecules were
identified through mass spectrometry, not through vibrational or rotational
spectroscopy.  These aromatic molecules may be precursors to more complex
PAHs. The identification of these specific aromatic molecules sheds light on
the composition of aromatic material within the ISM that will eventually be
incorporated into new stars and planets.

Nevertheless, identification of specific PAH molecules through IR spectroscopy
has so far been ineffective. The UIE bands are attributed to collective
emission from many different PAH species and thus do not allow one to
fingerprint individual PAH molecules. However, the IR spectra of different PAH
molecules are readily distinguishable in the laboratory from one species to
another. Therefore, in principle, IR spectroscopy could provide individual
identification if the molecule is sufficiently abundant and the telescope
instrument is sufficiently sensitive. With the advent of JWST, this may become
possible due to its unprecedented sensitivity. With an aim to offer
spectroscopic guidance for JWST to search for specific PAH molecules, we
perform a systematic exploration of the theoretical IR emission spectra of
various specific PAH species expected in different astrophysical environments.
In this work, we report the model IR emission spectra for cyanonaphthalenes
(i.e., 1- and 2-CNN), the first PAH molecules ever detected in the ISM. In §2
we synthesize the ultraviolet (UV) absorption cross sections which determine
how cyanonaphthalenes absorb starlight in space. The IR absorption cross
sections are also discussed in §2 which determine how cyanonaphthalenes emit
in the IR. The vibrational excitation and de-excitation of cyanonaphthalenes
are discussed in §3. We present in §4 the model IR emission spectra of
cyanonaphthalenes. The results are discussed in §5 and summarized in §6.

##  2 UV and IR Absorption Cross Sections

The UV absorption cross sections of cyanonaphthalenes determine how they
absorb starlight photons in space. Unfortunately, the UV absorption has only
been measured for 1-CNN over the wavelength range of $0.21\,{\rm\mu
m}<\lambda<0.33\,{\rm\mu m}$ (Perkampus 1992). Therefore, for $0.14\,{\rm\mu
m}<\lambda<0.21\,{\rm\mu m}$ we adopt the UV absorption of naphthalene (Robin
1975) and for $0.04\,{\rm\mu m}<\lambda<0.14\,{\rm\mu m}$ we take the vacuum
UV absorption measured by Koch et al. (1972) for naphthalene. Finally, as
shown in Figure 2, the UV absorption spectra of naphthalene are smoothly
scaled to join that of 1-CNN. As cyanonaphthalenes are too small to absorb
appreciably in the visible, for $\lambda>0.33\,{\rm\mu m}$ we therefore
extrapolate from that of Perkampus (1992) at $\lambda<0.33\,{\rm\mu m}$. Due
to lack of experimental data, we do not distinguish the UV absorption between
1-CNN and 2-CNN.

In the IR, Bauschlicher (1998) have computed the vibrational frequencies and
intensities of 1-CNN and 2-CNN, using the B3LYP density functional theory in
conjunction with the 4-31G basis set. We take the vibrational frequencies and
intensities of Bauschlicher (1998) which are available from the NASA Ames PAH
IR Spectroscopic Database (Boersma et al. 2014, Bauschlicher et al. 2018,
Mattioda et al. 2020). We represent each vibrational line by a Drude
function,222If we approximate the vibrational modes of PAHs as harmonic
oscillators, the absorption cross sections are then expected to be Drude
functions (see Li 2009).  characterized by the peak wavelength and intensity
of the vibrational transition. In addition, we assign a width of 30$\,{\rm
cm}^{-1}$ for each line, consistent with the natural line width expected from
a vibrationally excited PAH molecule (see Allamandola et al. 1999). The
resulting IR absorption sections of 1-CNN and 2-CNN are shown in Figure 3. It
is remarkable that the IR vibrational spectra of 1-CNN and 2-CNN are quite
similar at $\lambda<10\,{\rm\mu m}$ and even at $\lambda>10\,{\rm\mu m}$ they
do not differ substantially. This will be further discussed in §4.

##  3 Vibrational Excitation of Cyanonaphthalenes

Cyanonaphthalenes only have 19 atoms and 51 vibrational degrees of freedom.
Upon absorption of an UV stellar photon, cyanonaphthalenes will undergo
stochastic heating since their energy contents are often smaller than the
energy of a single stellar photon. We model the stochastic heating of
cyanonaphthalenes by employing the “exact-statistical” method of Draine & Li
(2001). We characterize the state of a cyanonaphthalene molecule (i.e., 1-CNN
or 2-CNN) by its vibrational energy $E$, and group its energy levels into
$(M+1)$ “bins”, where the $j$-th bin ($j$ = 0, …, $M$) is
$[E_{j,\min},E_{j,\max})$, with representative energy $E_{j}$ $\equiv$
$(E_{j,\min}$+$E_{j,\max})$/2, and width $\Delta E_{j}$ $\equiv$
$(E_{j,\max}$–$E_{j,\min})$. Let $P_{j}$ be the probability of finding 1-CNN
(or 2-CNN) in bin $j$ with energy $E_{j}$. The probability vector $P_{j}$
evolves according to

$dP_{i}/dt=\sum_{j\neq i}{\bf T}_{ij}P_{j}-\sum_{j\neq i}{\bf
T}_{ji}P_{i}~{}~{},~{}~{}i\,=\,0,...,M~{}~{},$ (1)

where the transition matrix element ${\bf T}_{ij}$ is the probability per unit
time for 1-CNN (or 2-CNN) in bin $j$ to make a transition to one of the levels
in bin $i$. We solve the steady state equations

$\sum_{j\neq i}{\bf T}_{ij}P_{j}=\sum_{j\neq i}{\bf
T}_{ji}P_{i}~{}~{},~{}~{}i\,=\,0,...,M~{}~{}$ (2)

to obtain the $M$+1 elements of $P_{j}$, and then calculate the resulting IR
emission spectrum (see eq. 55 of Draine & Li 2001).

In calculating the state-to-state transition rates ${\bf T}_{ji}$ for
transitions $i$$\rightarrow$$j$, we distinguish the excitation rates ${\bf
T}_{ul}$ (from $l$ to $u$, $l$ $<$ $u$) from the deexcitation rates ${\bf
T}_{lu}$ (from $u$ to $l$, $l$ $<$ $u$). For a given starlight energy density
$u_{E}$, the rates for upward transitions $l$$\rightarrow$$u$ (i.e., the
excitation rates) are just the photon absorption rates:

${\bf T}_{ul}\approx C_{{\rm abs}}(E)\,c\,u_{E}\Delta
E_{u}/(E_{u}-E_{l})~{}~{}.$ (3)

The rates for downward transitions $u$$\rightarrow$$l$ (i.e., the deexcitation
rates) can be determined from the detailed balance analysis of the Einstein
$A$ coefficient:

${\bf T}_{lu}\approx\frac{8\pi}{h^{3}c^{2}}\frac{g_{l}}{g_{u}}\frac{\Delta
E_{u}}{E_{u}-E_{l}}E^{3}\times C_{{\rm abs}}(E)\left[1+\frac{h^{3}c^{3}}{8\pi
E^{3}}u_{E}\right]~{}~{},$ (4)

where $h$ is the Planck constant, and the degeneracies $g_{u}$ and $g_{l}$ are
the numbers of energy states in bins $u$ and $l$, respectively:

$g_{j}\equiv
N(E_{j,\max})-N(E_{j,\min})\approx\left(dN/dE\right)_{E_{j}}\Delta
E_{j}~{}~{},$ (5)

where $\left(dN/dE\right)_{E_{j}}$ is the vibrational density of states of
1-CNN (or 2-CNN) at internal energy $E_{j}$. For a cyanonaphthalene molecule,
if we know the frequencies of all its 51 vibrational modes, we can employ the
Beyer-Swinehart numerical algorithm (Beyer & Swinehart 1973, Stein &
Rabinovitch 1973) to calculate the vibrational density of states
$\left(dN/dE\right)_{E_{j}}$ and therefore the degeneracies $g_{j}$ for each
vibrational energy bin. Meanwhile, if we know the oscillator strength of each
vibrational mode, we can obtain the IR absorption cross section $C_{\rm
abs}(E)$ by summing up all the vibrational transitions with each approximated
as a Drude profile. Without a prior knowledge of the width of each vibrational
transition, we assign a width of 30$\,{\rm cm}^{-1}$, consistent with the
natural line width expected from free-flying molecule (see Allamandola et al.
1999). This natural line width arises from intramolecular vibrational energy
redistribution.

As mentioned earlier (see §2), Bauschlicher (1998) have computed with
B3LYP/4-31G the frequency and intensity for each of the 51 vibrational
transitions of 1-CNN and 2-CNN. Therefore, we can calculate the vibrational
density of states and the degeneracy $g_{j}$ for each vibrational energy bin
as well as the IR absorption cross section $C_{\rm abs}(E)$ by adopting the
vibrational frequencies and intensities of 1-CNN and 2-CNN computed by
Bauschlicher (1998). For a given astrophysical environment characterized by
its starlight energy density $u_{E}$, we first calculate the excitation rates
${\bf T}_{ul}$ (from $l$ to $u$, $l$ $<$ $u$) according to eq. 3, using the UV
absorption cross setions described in §2. We then calculate the deexcitation
rates ${\bf T}_{lu}$ (from $u$ to $l$, $u$ $>$ $l$) according to eq. 4, using
the degeneracy $g_{j}$ and the IR absorption cross section $C_{\rm abs}(E)$
derived from the vibrational frequencies and intensities of Bauschlicher
(1998). With the state-to-state transition rates ${\bf T}_{ji}$ determined, we
solve the steady-state probability evolution equation (see eq. 2) to obtain
the steady-state energy probability distribution $P_{j}$ and finally calculate
the resulting IR emission spectrum, according to eq. 55 of Draine & Li (2001).
For computational convenience, we consider 500 energy bins (i.e., $M=500$).

##  4 Model Emission Spectra

We first consider cyanonaphthalenes in the diffuse ISM excited by the solar
neighborhood interstellar radiation field of Mathis et al. (1983; hereafter
MMP83). Figure 4a shows the IR emissivities (erg$\,{\rm s}^{-1}\,{\rm
sr}^{-1}\,{\rm cm}^{-1}$) per molecule for 1-CNN and 2-CNN. A first glance of
Figure 4a reveals that the IR emission spectra of 1-CNN and 2-CNN are rather
similar. Both exhibit a prominent C–H stretching band at 3.26$\,{\rm\mu m}$, a
prominent C–N stretching band at 4.69$\,{\rm\mu m}$, and a broad complex at
$\sim\,$6–9$\,{\rm\mu m}$ consisting of a number of sub-features arising from
C–C stretching and C–H in-plane bending vibrations. The C–H out-of-plane
bending bands at $\sim\,$10.5–13.5$\,{\rm\mu m}$ are also similar, except for
1-CNN they coalesce into a broad band at 12.7$\,{\rm\mu m}$, while for 2-CNN
several sub-features remain to be noticeable. The C-C-C skeletal bending band
occurs at a longer wavelength for 1-CNN ($\sim\,$22$\,{\rm\mu m}$) compared to
2-CNN ($\sim\,$21$\,{\rm\mu m}$), otherwise they are also similar. Because of
the close similarity between the emission spectra of 1-CNN and 2-CNN, in the
following we will only consider 1-CNN.

To verify the single-photon heating nature of cyanonaphthalenes, we calculate
the IR emission of 1-CNN excited by enhanced MMP83-type radiation fields.
Figure 4b shows the IR emission spectra of 1-CNN excited by the MMP83
radiation field, but enhanced by a factor of $U$ = 1000 and 105 ($U$ = 1
corresponds to the MMP83 radiation field). As illustrated in Figure 4b, when
scaled by $U$, the IR emission spectra are identical. This justifies the
stochastic heating treatment of the vibrational excitation of
cyanonaphthalenes. In the single-photon heating regime, the resulting IR
emission is only dependent on the energy of the illuminating starlight
photons, not on the starlight intensity.

We also examine how the IR emission of cyanonaphthalenes depends on the
starlight spectrum (i.e., the “hardness” of the exciting radiation field). To
this end, we consider 1-CNN exposed to stars with different effective
temperatures: $T_{\rm eff}$ = 40,000, 22,000, 8,000$\,{\rm K}$. The starlight
spectra are approximated by the Kurucz model atmospheric spectra. Figure 5a
shows the IR emission of 1-CNN illuminated by stars of $T_{\rm eff}$ =
40,000$\,{\rm K}$ (like O6V stars) with an intensity of $U$ = 104, where $U$
is defined as

$U=\frac{\int_{1\mu{\rm m}}^{912{\rm\,{\rm\AA}}}4\pi J_{\star}(\lambda,T_{\rm
eff})\,d\lambda}{\int_{1\mu{\rm m}}^{912{\rm\,{\rm\AA}}}4\pi J_{\rm
ISRF}(\lambda)\,d\lambda}~{}~{},$ (6)

where $J_{\rm ISRF}(\lambda)$ is the MMP83 radiation intensity, and
$J_{\star}(\lambda,T_{\rm eff})$ is the intensity of starlight approximated by
the Kurucz model atmospheric spectrum. Such a starlight spectrum and intensity
resemble that of the Orion Bar photodissociation region and the M17 star-
forming region. A comparison of Figure 5a with Figure 4b reveals that the IR
emission spectra are almost identical for 1-CNN excited by the MMP83 field and
by stars of $T_{\rm eff}$ = 40,000$\,{\rm K}$, except the emissivity level is
higher for the latter since the mean absorbed photon energy is higher for the
latter.

Figure 5b shows the IR emission of 1-CNN illuminated by stars of $T_{\rm eff}$
= 22,000$\,{\rm K}$ (like B1.5V stars) with an intensity of $U$ = 103. This
would be the case if 1-CNN is located in the reflection nebula NGC 2023. If
scaled by $U$ and the mean absorbed photon energy, the IR emission spectrum
for $T_{\rm eff}$ = 22,000$\,{\rm K}$ is essentially identical to that for
$T_{\rm eff}$ = 40,000$\,{\rm K}$.

Figure 5c shows the IR emission of 1-CNN illuminated by stars of $T_{\rm eff}$
= 8,000$\,{\rm K}$ (like A5V stars) with an intensity of $U$ = 105. The Red
Rectangle protoplanetary nebula is illuminated by such a starlight spectrum
and intensity.333It is worth noting that Witt et al. (2009) argued that the
Red Rectangle protoplanetary nebula is illuminated by a binary, and the UV
radiation in this nebula does not come from HD 44179 of which $T_{\rm
eff}\approx 8,000\,{\rm K}$, but rather from the accretion disk surrounding
the secondary star in this binary, with $T_{\rm eff}$ in the range of
$\sim\,$17,000–25,000$\,{\rm K}$.  Cyanonaphthalenes in the Red Rectangle
would be excited by stellar photons much softer than that in the Orion Bar. It
is therefore not unexpected that, as illustrated in Figure 5c, 1-CNN in the
Red Rectangle would emit somewhat more at longer wavelengths compared to that
in the Orion Bar. Nevertheless, the overall spectral shape does not differ
considerably between that of the Red Rectangle and that of the Orion Bar.

Finally, we also consider the TMC-1 molecular cloud where cyanonaphthalenes
were first detected through their rotational lines. The TMC-1 cloud is
externally illuminated by the general interstellar radiation field and has a
total visual extinction of $\sim\,$3.6$\,{\rm mag}$ (Whittet et al. 2004).
Therefore, to calculate the IR emission of cyanonaphthalenes in the TMC-1
cloud, we take the MMP83 interstellar radiation field, but attenuated by dust
extinction with $A_{V}$ = 1.8$\,{\rm mag}$ (from the cloud surface to the
cloud core) and a wavelength-dependence $A_{\lambda}/A_{V}$ like that of the
Galactic average extinction curve of $R_{V}=3.1$, where $R_{V}$ is the total-
to-selective extinction ratio (see Cardelli et al. 1989), i.e.,
cyanonaphthalenes in the TMC-1 cloud are assumed to be excited by starlight
with an intensity of $J_{\rm
ISRF}(\lambda)\times\exp\\{-\left(A_{V}/1.086\right)\times\left(A_{\lambda}/A_{V}\right)\\}$.
Figure 5d shows the IR emission of 1-CNN calculated for the TMC-1 cloud.
Compared to that excited by the MMP83 radiation field (see Figure 4b), the
overall IR emission spectrum of 1-CNN expected in the TMC-1 cloud is closely
similar, except the emissivity level is appreciably reduced and the
3.26$\,{\rm\mu m}$ C–H stretch and 4.69$\,{\rm\mu m}$ C–N stretch emit
slightly less.

##  5 Discussion

It is gratifying that, as illustrated in Figure 5a–d, the IR emission spectra
of cyanonaphthalenes are not sensitive to the illuminating starlight spectrum.
This is because, as shown in Figure 6, due to their small heat contents, the
energy probability distribution functions of cyanonaphthalenes closely
resemble each other, upon excited by starlight of different intensities and
different spectral shapes.

While observationally the identification of cyanonaphthalenes through the
$\sim\,$6–9$\,{\rm\mu m}$ complex and the 10.5–13.5$\,{\rm\mu m}$ complex
could be complicated by the 6.2, 7.7, 8.6, 11.3 and 12.7$\,{\rm\mu m}$ UIE
bands, the detection of the 3.26 and 4.69$\,{\rm\mu m}$ bands, and, to a less
degree, the 21 or 22$\,{\rm\mu m}$ band, could potentially allow one to
identify cyanonaphthalenes in space. Prior to JWST, the detection and
identification of the spectral features as possibly due to cyanonaphthalenes
would have been hampered by their small intensities which would put them at
the limit of modern observational techniques, including the Short Wavelength
Spectrometer (SWS) on board the Infrared Space Observatory (ISO) and the
Infrared Camera (IRC) on board AKARI. With the advent of JWST, this will
change. The Near InfraRed Spectrograph (NIRSpec) and Mid-Infrared Instrument
(MIRI) on JWST span the wavelength range of the characteristic vibrational
bands of cyanonaphthalenes. JWST’s unique high sensitivity and high resolution
capabilities will open up an IR window unexplored by Spitzer444The Infrared
Spectrograph (IRS) on Spitzer only operates longward of $\sim\,$5.2$\,{\rm\mu
m}$ and does not cover the characteristic 3.26$\,{\rm\mu m}$ C–H and
4.69$\,{\rm\mu m}$ C–N stretching bands of cyanonaphthalenes.  and unmatched
by ISO observations and thus could potentially place the detection of the IR
vibrational bands of individual PAH molecules on firm ground.

Admittedly, it is not clear if the detection of the 3.26$\,{\rm\mu m}$ C–H and
4.69$\,{\rm\mu m}$ C–N stretching bands could uniquely pinpoint the presence
of cyanonaphthalenes, as other cyano-substituted PAHs may also emit at similar
wavelengths. Note that cyano-benzene (benzonitrile) and cyano-indene have also
been detected in TMC-1 (McGuire et al. 2018, Sita et al. 2022). To this end, a
systematic calculation of the IR emission spectra of these molecules would be
crucial. One can imagine that different-sized cyano-containing aromatic
molecules would exhibit different C–H/C–N band ratios since, with different
energy contents, they are expected to be excited to different energy levels by
the same stellar photon.

Let $N_{\rm H}$ be the hydrogen column density along the line of sight to an
astrophysical region. Let $\left[{\rm C/H}\right]_{\rm CNN}$ be the number of
carbon atoms (per H nucleon) locked up in a cyanonaphthalene molecule. The
intensity of the IR emission (erg$\,{\rm s}^{-1}\,{\rm cm}^{-3}\,{\rm
sr}^{-1}$) expected from cyanonaphthalenes would be

$I_{\lambda}=j_{\lambda}\times\left(\frac{N_{\rm
H}}{11}\right)\times\left[\frac{\rm C}{\rm H}\right]_{\rm CNN}~{}~{},$ (7)

where the denominator “11” accounts for the fact that a cyanonaphthalene
molecule has 11 carbon atoms. The column density of cyanonaphthalenes is
simply

$N_{\rm CNN}=\left(\frac{N_{\rm H}}{11}\right)\times\left[\frac{\rm C}{\rm
H}\right]_{\rm CNN}~{}~{}.$ (8)

Therefore, by comparing the observed intensity $I_{\lambda}^{\rm obs}$ with
the model emissivity $j_{\lambda}$ (e.g., see Figure 5a–d), one can derive
$N_{\rm CNN}$. If $N_{\rm H}$ is known, then one can determine how much carbon
is locked up in cyanonaphthalenes. On the other hand, if $N_{\rm CNN}$ is
known (e.g., for 1-CNN or 2-CNN in TMC-1 from the measurements of their
rotational lines), one can predict the IR emission intensity $I_{\lambda}$ by
multiplying the model emissivity shown in Figure 5a–d with $N_{\rm CNN}$.

While cyanonaphthalenes are present in the UV-attenuated TMC-1 molecular
cloud, it is not clear if they can survive in hostile regions where the UV
radiation is intense (e.g., the diffuse ISM, the Orion Bar and the M17 star-
forming cloud). According to Jacovella et al. (2022), HCN removal from
protonated benzonitrile requires only $\sim\,$3.4$\,{\rm eV}$. As similar HCN
removal energies are expected for cyanonaphtalenes, the survivability of these
species in UV-irradiated environments such as the Orion Bar is questionable.
It has been argued that only PAHs with more than 20 carbon atoms may survive
in these regions (see Tielens 2008). However, Stockett et al. (2023) recently
found that cyanonaphthalene can be efficiently stabilized following
ionization, with the aid of the so-called “Recurrent Fluorescence” (also known
as Poincaré fluorescence, see Léger et al. 1988), a radiative relaxation
channel in which optical photons are emitted from thermally populated
electronically excited states. Iida et al. (2022) also found that the
Recurrent Fluorescence could efficiently stabilize small cationic carbon
clusters with as few as nine carbon atoms. The Recurrent Fluorescence of PAHs
excited by UV photons could reach high quantum efficiencies in photon
conversion and has been suggested as the source of the so-called “extended red
emission” (ERE; Witt & Lai 2020).

If the Recurrent Fluorescence does occur, the absorbed stellar photon energy
will not be completely released through vibrational transitions. Therefore,
the IR vibrational emissivity is expected to be somewhat lower than that
calculated in Figure 5a–d, roughly by a factor of $\langle h\nu\rangle_{\rm
RF}/\langle h\nu\rangle_{\rm abs}$, where $\langle h\nu\rangle_{\rm RF}$ is
the mean energy of the optical photons emitted due to Recurrent Fluorescence,
and $\langle h\nu\rangle_{\rm abs}$ is the mean enegy of the absorbed stellar
photons. For cyanonaphthalenes in the diffuse ISM, the reflection nebula NGC
7023, and the Orion Bar, $\langle h\nu\rangle_{\rm abs}$ is in the order of
8–9$\,{\rm eV}$. With $\langle h\nu\rangle_{\rm RF}\approx 2\,{\rm eV}$ (see
Stockett et al. 2023), the IR emissivity $j_{\lambda}$ will only be reduced by
$\sim\,$20% and the emission spectral profile will largely remain the same.

Finally, we note that the UV absorption cross sections of cyanonaphthalenes
(i.e., 1-CNN and 2-CNN) adopted in this work were “synthesized” from 1-CNN and
naphthalene (see §2). It is not clear how different the actual UV absorption
cross sections of 1-CNN would be compared with naphthalene and of 2-CNN
compared with 1-CNN. As a priori, one would imagine that the incorporation of
a CN group in naphthalene, the parental molecule of cyanonaphthalenes, would
cause significant changes to the UV absorption spectrum, as the $\pi$
electrons of the CN group are partially conjugated with the main aromatic
$\pi$ system. Nevertheless, as mentioned earlier (see Figure 5a–d), the IR
emission spectra of cyanonaphthalenes are not sensitive to the illuminating
starlight spectrum. This also implies that the IR emission spectral shape of
cyanonaphthalenes is not very sensitive to the exact UV absorption. In any
case, future experimental measurements of the absorption spectra of
cyanonaphthalenes from the optical to the far-UV will be helpful in examining
their vibrational excitation.

We also note that the IR vibrational frequencies and intensities of
cyanonaphthalenes adopted here were computed by Bauschlicher (1998), using a
DFT force field and intensitie in the harmonic limit. It is well recognized
that anharmonic interactions could appreciably affect the vibrational
frequencies and intensities. In addition, Fermi interactions can have a
significant impact on the C–H stretches. Moreover, quadratic dipole dependence
can also affect the vibrational intensities. While some of these effects could
be grossly accounted for by empirical scaling coefficients as applied by
Bauschlicher (1998), more accurate results may be obtained with the standard
second-order vibrational perturbation theory (VPT2). Note that the VPT2 theory
is feasible with DFT force fields for molecules as large as cyanonaphthalenes.
Future quantum chemical computations utilizing the VPT2 theory and even
better, experimental measurements of gas-phase cyanonaphthalenes will be
helpful in evaluating the accuracy of the IR emission spectra calculated here.

In calculating the IR emission sepctra of cyanonaphthalenes, we assume a
natural line width of $\gamma=30\,{\rm cm}^{-1}$ for each vibrational
transition (see §2) due to internal vibration redistribution (IVR). This may
be somewhat simplified. Molecules like cyanonaphthalenes are in a somewhat
intermediate size-range, where IVR is prevalent, but not necessarily ergodic
or completely statistical. For these molecules, spectral broadening may more
likely be dominated by anharmonic shifts of the thermally excited vibrational
levels. Indeed, Hewett et al. (1994) argued that the C–H fundamentals of
naphthalene are far from the IVR threshold. In view of these, we have also
considered $\gamma=10\,{\rm cm}^{-1}$ and $50\,{\rm cm}^{-1}$. As shown in
Figure 7, the model IR emission spectra are essentially the same, particularly
in the characteristic C–H and C–N stretching regions. This justifies the
choice of a natural line width of $\gamma=30\,{\rm cm}^{-1}$.

##  6 Summary

We have modeled the vibrational excitation of cyanonaphthalenes, the first
specific PAH species ever identified in space, and have calculated their IR
emission spectra for a number of representative astrophysical regions, from
the diffuse ISM to regions illuminated by stars of different effective
temperatures. Also calculated includes the TMC-1 dark molecular cloud where
cyanonaphthalenes were first detected through their rotational lines. It is
found that while the emissivity level varies from one region to another, the
overall IR emission spectra do not vary much with environments or between two
isomers of cyanonaphthalene (i.e., 1-CNN and 2-CNN). Both isomers exhibit
prominent bands at 3.26$\,{\rm\mu m}$ (C–H stretch), 4.69$\,{\rm\mu m}$ (C–N
stretch), 21 or 22$\,{\rm\mu m}$ (skeletal bend), and broad complexes at
$\sim\,$6–9$\,{\rm\mu m}$ and $\sim\,$10.5–13.5$\,{\rm\mu m}$. The NIRSpec and
MIRI instruments on board JWST are well suited for searching for
cyanonaphthalenes in space, particularly through the characteristic bands at
3.26 and 4.69$\,{\rm\mu m}$.

We thank the anonymous referees for helpful comments and suggestions. We thank
B.M. Broderick, B.T. Draine, B.A. McGuire, and E.F. van Dishoeck for
stimulating discussions. KJL and TTF are supported by the National Key R&D
Program of China under No. 2017YFA0402600, and the NSFC grants 11890692,
12133008, and 12221003, as well as CMS-CSST-2021-A04. XJY is supported in part
by NSFC 12333005 and 12122302 and CMS-CSST-2021-A09.

## References

  * (1) Allamandola, L.J., Tielens, A.G.G.M., & Barker, J.R. 1985, ApJ, 290, L25
  * (2) Allamandola, L.J., Hudgins, D.M., & Sandford, S.A. 1999, ApJ, 511, 115
  * (3) Bauschlicher, C.W. 1998, Chem. Phys., 234, 87
  * (4) Bauschlicher, C. W., Ricca, A., Boersma, C., & Allamandola, L.J. 2018, ApJS, 234, 32
  * (5) Beyer, T., & Swinehart, D.F. 1973, Commun. Assoc. Comput. Machinery, 16, 379
  * (6) Boersma, C., Bauschlicher, C.W., Ricca, A., et al. 2014, ApJS, 211, 8
  * (7) Burkhardt, A. M., Lee, K.L.K., Changala, P.B., et al. 2021, ApJL, 913, L18
  * (8) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 354, 245
  * (9) Cernicharo, J., Heras, A. M., Tielens, A. G. G. M., et al. 2001, ApJL, 546, L123
  * (10) Clemett, S.J., Maechling, C.R., Zare, R.N., Swan, P.D., & Walker, R.M. 1993, Science, 262, 721
  * (11) Clemett, S. J., Sandford, S. A., Nakamura-Messenger, K., Hörz, F., & McKay, D.S. 2010, Meteorit. Planet. Sci., 45, 701
  * (12) Draine, B.T., & Li, A. 2001, ApJ, 551, 807
  * (13) Hewett, K. B., Shen, M., Brummel, C. L., et al. 1994, J. Chem. Phys., 100, 4077.
  * (14) Iida, S., Hu, W., Zhang, R., et al. 2022, MNRAS, 514, 844
  * (15) Jacovella, U., Noble, J. A., Guliani, A., et al. 2022, A&A, 657, A85
  * (16) Koch, E.E., Otto, A., & Radler, K. 1972, Chem. Phys. Lett., 16, 131
  * (17) Léger, A., & Puget, J. 1984, A&A, 137, L5
  * (18) Léger, A., Boissel, & d’Hendecourt, L.B. 1988, Phys. Rev. Lett., 60, 921
  * (19) Li, A. 2009, Optical Properties of Dust, in Small Bodies in Planetary Sciences, ed. I. Mann, A. Nakamura, & T. Mukai, Springer, 167
  * (20) Li, A. 2020, Nature Astronomy, 4, 339
  * (21) Mathis, J.S., Mezger, P.G., & Panagia, N. 1983, A&A, 128, 212
  * (22) Mattioda, A.L., Hudgins, D.M., Boersma, C., et al. 2020, ApJS, 251, 22
  * (23) McGuire, B. A., Burkhardt, A. M., Kalenskii, S., et al. 2018, Science, 359, 202
  * (24) McGuire, B. A., Loomis, R. A., Burkhardt, A. M., et al. 2021, Science, 371, 1265
  * (25) Naraoka, H., Takano, Y., Dworkin, J.P., et al. 2023, Science, 379, 9033
  * (26) Perkampus, H.H. 1992, UV-VIS Atlas of Organic Compounds (2nd ed.), VCH, Weinheim
  * (27) Robin, M.B. 1975, Higher Excited States of Polyatomic Molecules, Vol. II (London: Academic Press)
  * (28) Sandford, S.A., Aléon, J., Alexander, C.M.O.D., et al. 2006, Science, 314, 1720
  * (29) Schuhmann, M., Altwegg, K., Balsiger, H., et al. 2019, A&A, 630, A31
  * (30) Sita, M. L., Changala, P. B., Xue, C., et al. 2022, ApJL, 938, L12
  * (31) Stein, S.E., & Rabinovitch, B.S. 1973, J. Chem. Phys., 58, 2438
  * (32) Stockett, M. H., Bull, J. N., Cederquist, H., et al. 2023, Nature Communications, 14, 395
  * (33) Tabone, B., Bettoni, G., van Dishoeck, E.F., et al. 2023, Nature Astronomy, 7, 805
  * (34) Tielens, A. G. G. M. 2008, ARA&A, 46, 289
  * (35) Whittet, D. C. B., Shenoy, S. S., Clayton, G. C., & Gordon, K.D. 2004, ApJ, 602, 291
  * (36) Witt, A. N., Vijh, U. P., Hobbs, L. M., Aufdenberg, J.P., Thorburn, J.A., & York, D.G. 2009, ApJ, 693, 1946
  * (37) Witt, A. N., & Lai, T. S.-Y. 2020, Ap&SS, 365, 58

Figure 1:  Chemical structures of specific aromatic molecules identified in
the ISM.  Figure 2:  UV absorption cross sections (per C atom) of
cyanonaphthalenes “synthesized” from the experimental data of 1-CNN over
$0.21\,{\rm\mu m}<\lambda<0.33\,{\rm\mu m}$ (red line; Perkampus 1992), and of
naphthalene over $0.14\,{\rm\mu m}<\lambda<0.21\,{\rm\mu m}$ (black line;
Robin 1975) and over $0.04\,{\rm\mu m}<\lambda<0.14\,{\rm\mu m}$ (blue line;
Koch et al. 1972).  Figure 3:  IR absorption cross sections (per C atom) of
cyanonaphthalenes (red line: 1-CNN; blue line: 2-CNN) calculated from the
vibrational frequences and intensities computed by Bauschlicher (1998) using
B3LYP/4-31G. Each vibrational line is assigned a width of 30$\,{\rm cm}^{-1}$,
consistent with the natural line width expected from a vibrationally excited
PAH molecule (see Allamandola et al. 1999).

Figure 4:  Left panel (a): IR emission spectra of 1-CNN (red line) and 2-CNN
(blue line) illuminated by the MMP83 ($U$ = 1) interstellar radiation field.
Right panel (b): IR emission spectra of 1-CNN illuminated by the MMP83
radiation field with different intensities (red solid line: $U$ = 1; blue
shot-dashed line: $U$ = 1,000; black long-dashed line: $U$ = 105).

Figure 5:  IR emission spectra of 1-CNN illuminated by radiation fields of
different starlight spectra: $T_{\rm eff}$ = 40,000$\,{\rm K}$ and $U$ = 104
(top left panel, a), $T_{\rm eff}$ = 22,000$\,{\rm K}$ and $U$ = 104 (top
right panel, b), $T_{\rm eff}$ = 8,000$\,{\rm K}$ and $U$ = 105 (bottom left
panel, c), and MMP83 radiation field attenuated by a visual extinction of
$A_{V}=1.8\,{\rm mag}$ (bottom right panel, d).  Figure 6:  Vibrational energy
probability distribution functions for the excited vibrational states of
1-CNN, exposed to radiation fields of different intensities ($U$) and
different spectral shapes: MMP83 (red), attenuated MMP83 (black), and stars of
different effective temperatures (blue: $T_{\rm eff}=4\times 10^{4}\,{\rm K}$
and $U=10^{4}$, cyan: $T_{\rm eff}=2.2\times 10^{4}\,{\rm K}$ and $U=10^{3}$,
and magenta: $T_{\rm eff}=8\times 10^{3}\,{\rm K}$ and $U=10^{5}$).

Figure 7:  Left panel (a): Comparison of the IR emission spectra of 1-CNN of a
natural line width of $\gamma=30\,{\rm cm}^{-1}$ (red line) for each
vibrational transition with that of $\gamma=10\,{\rm cm}^{-1}$ (blue line).
The molecule is illuminated by the MMP83 ($U$ = 1) interstellar radiation
field. Right panel (b): Same as (a) but for a comparison of $\gamma=30\,{\rm
cm}^{-1}$ (red line) with $\gamma=50\,{\rm cm}^{-1}$ (blue line).


#  Bouncing universe for deformed non-minimally coupled inflation model

S. Upadhyay<EMAIL_ADDRESS>Department of Physics, K.L.S. College,
Magadh University, Nawada-805110, Bihar, India

###### Abstract

In this paper, we consider a non-minimally coupled gravity model to study the
bouncing universe. The holographic principle has various effects on the
bouncing universe. We choose some suitable new variables and achieve the new
Hamiltonian and Lagrangian which have harmonic oscillator form. The
corresponding Lagrangian is deformed by non-commutative geometry. In order to
have a solution for the bouncing universe we specify the potential in the
equation state. In that case, we draw the equation of state in terms of time
and show that the equation of state crosses $-1$. Such bouncing behavior leads
us to apply some conditions on $\theta$ and $\beta$ from non-commutative
geometry. Here, we can also check the system’s stability due to the
deformation of the non-minimally coupled gravity model. In order to examine
the stability of the system we obtain the variation of pressure with respect
to density energy. Also, we draw the variation of pressure with respect to
energy density and show the stability condition.

PACS: 04.50.Kd; 04.20.Fy; 95.36.+X.

Keywords: Non-commutative geometry; Non-minimally coupled gravity; Bouncing
universe; Dark energy.

##  I Introduction

The moving from the accelerated collapse era to the expanding era without a
singularity lead us to have a bouncing universe. The bounce realization due to
the holographic principle in the early universe studied recently [1]. The
holographic dark energy plays a crucial role in describing the late-time
universe, and found in agreement with observations [2, 3, 4, 5, 6]. In that
case, many cosmological solutions for the bouncing universe tested in the last
decades [7]. The subject of bouncing investigated in several model as brane-
world scenarios of cosmology. Here, if the bulk space is taken to be a charged
AdS black hole, then we have the bouncing universe [8]. In this case, brane
makes a contracting to expanding phase transition. As we know from the
singularity theorem in four-dimensions, such a bounce cannot occur because of
the certain energy condition. Hence, in order to produce the bouncing
universe, it is necessary that the bulk geometry contribute a negative energy
density to the effective stress-energy tensor of the brane configuration [9].
As we know, several authors give further developments for these bouncing
brane-world scenario [10, 11, 12, 13, 14, 15, 16]. Now, it is believed that
the universe is in accelerated expansion phase. For this cosmic acceleration,
we have seen some observations from type Ia supernovae [17, 18, 19] in
associated with Large Scale Structure (LSS) [20, 21] and Cosmic Microwave
Background (CMB) anisotropy [22, 23]. On the other hand the dark energy with
negative pressure lead us to this cosmic acceleration. In order to have such
phenomena, the theories trying to modify Einstein equations. On the other
hand, in modern cosmology, the dark energy can be explained by equation of
state parameter as $\omega=\frac{p}{\rho}$, where $p$ and $\rho$ are pressure
and energy density, respectively. In that case, for the accelerated expansion
the equation of state dark energy must be satisfied by $\omega<-\frac{1}{3}$.
As we know, there are two puzzles from cosmological constant as dark energy.
One of them is that, cosmological constant is about 120 orders of magnitude
smaller than its natural expectation. The other is cosmological coincidence.
In order to solve these problems, various model of dark energy have been
introduced such as quintessence [24, 25, 26] or k-essence [27], scaler-field
dark energy model including tachyon [28, 29], ghost condensate [30, 31] and
quintom [32, 33, 34, 35, 36, 37, 38, 39]. Also there is a unified dark energy-
dark matter model based on Chaplygin gas equation of state and its extensions
[40, 41, 42, 43, 44, 45, 46, 47, 48, 49]. We note here the analysis of the
dark energy properties from recent observation mildly favour model with
$\omega$ crossing $-1$ in the near past. In that case there are other
proposals to describe accelerating expansion of the universe include
interacting dark energy models [50, 51], brane-world models [52, 53], and
holographic dark energy models [54, 55, 56, 57]. The difficulty of realizing
$\omega$ crossing over $-1$ in quintessence and phantom-like models for the
first time have been shown by Refs. [33, 36, 37]. Because when the equation of
state $\omega$ approaches to $-1$, the dark energy perturbation would be
divergent [58, 59, 60, 61]. The most important point here is that the quintom
scenario of dark energy very well understand the nature of dark energy with
$\omega$ across $-1$ [62, 63].
All above physical information give us motivation to work non-minimally
coupled gravity with deformed phase space. As we know, such non-minimal
coupling is very successful [64, 65]. The simple example of such model is
Higgs boson which act the inflation field in the form of $\xi\varphi^{2}R$. In
that case $\varphi$ and $R$ are Higgs field and Ricci scaler respectively.
Here, we are going to study the dark energy and bouncing universe $\omega$
with usual metric formalism in non-minimally coupling model. We deformed the
non-minimally coupling to gravity model. Then, we compare the results of
deformed and non-deformed corresponding model to each other. Here, we note
that the bouncing universe and dark energy in non-minimally coupled gravity
with deformed phase space lead us to have some new results. Here, another
deformed approach will be Finsler geometry. In future one can realize the
relation between the above deformed method, non-commutative geometry and
Finsler geometry.
So, in section II, we introduce the non-minimal coupling to the gravity model.
Also, we take Friedmann-Roberson-Walker (FRW) background and write the
corresponding Lagrangian. In section III, first of all we assume
$\xi=\frac{-1}{6}$, $\kappa=0$ and write the corresponding Lagrangian. In
order to deform such a model we chose some changes of variables. Also, again
we obtain the Lagrangian and Hamiltonian of the non-minimal model by new
variables. The corresponding Hamiltonian will be form of harmonic oscillator,
and it is useful to investigate the gravitational theories. On the other hand,
we review some non-commutative (NC) geometries [66]. We transform the
classical phase space variables of $x_{i}$, $y_{i}$, $P_{x_{i}}$ and
$P_{y_{i}}$. Also, we achieve some deformed algebra for the new variables
which is satisfied by anti-commutative relations. Here, the NC geometry helps
us to write the Hamiltonian and Lagrangian for the non-minimal coupling to
gravity model with the new variable. In section IV, we have some cosmology
solutions for the deformed non-minimally model and investigate the dark energy
and bouncing universe. The equation of state plays important role in the
bouncing universe. Here also, we have some figures for the equation of state
cross $-1$. In the last section, we have some results and suggestions for the
deformed model in the bouncing universe.

##  II Lagrangian of non-minimal coupling to gravity

Now, we consider the general model of inflation with non-minimal coupling
which has a following action,

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{1}{2}(M^{2}+\xi\varphi^{2})g^{\mu\nu}R_{\mu\nu}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi-V(\varphi)\right],$
(1)

where $M$, $R_{\mu\nu}$, $\xi$ and $\varphi$ are mass scales, Ricci tensor,
dimensionless coupling constant and the inflation field respectively. If we
back to the Einstein frame, we need the following conformal transformation,

$g_{\mu\nu}\rightarrow\Omega(\varphi)^{-1}g_{\mu\nu},$ (2)

where

$\Omega(\varphi)=\frac{M^{2}+\xi\varphi^{2}}{M_{p}^{2}}.$ (3)

We do not need such transformation because we take general form of action (1),
and employ the FRW metric background with $\kappa=0$. Here, we try first to
study the bouncing universe and dark energy to usual non-minimal coupling
model (1). In order to have some cosmology solutions, one can write the
corresponding Lagrangian with the following equation,

$\mathcal{L}=3M^{2}(\kappa
a-a\dot{a}^{2})+3\xi(-a\dot{a}^{2}\varphi^{2}-2a^{2}\dot{a}\varphi\dot{\varphi}+\kappa
a\varphi^{2})+\frac{1}{2}a^{3}\dot{\varphi}^{2}-V(\varphi)a^{3},$ (4)

where we used,

$\displaystyle\sqrt{-g}$ $\displaystyle=$ $\displaystyle a^{3},$
$\displaystyle R$ $\displaystyle=$ $\displaystyle
6(\frac{\ddot{a}}{a}+\frac{\kappa+\dot{a}^{2}}{a^{2}}).$ (5)

The equation (4) is the Lagrangian of non-minimal coupling to gravity. In the
next section, we are going to take the above Lagrangian and apply some
transformations to the phase space variables.

##  III Deformation of non-minimal coupling to gravity model

Now, we take the equation (4) and assume $\xi=\frac{-1}{6}$ and $\kappa=0$, to
rewrite the Lagrangian as following,

$\mathcal{L}=3M^{2}(-a\dot{a}^{2})+\frac{1}{2}(a\dot{a}^{2}\varphi^{2}+2a^{2}\dot{a}\varphi\dot{\varphi}+a^{3}\dot{\varphi}^{2})-V(\varphi)a^{3}.$
(6)

Here, we choose the following change of variables,

$\begin{array}[]{ccc}x_{1}=a\sinh(\frac{\varphi
a}{\sqrt{2}}),&x_{2}=a\cosh\sqrt{3a}M,&x_{3}=\sinh(\sqrt{V(\varphi)}a^{\frac{3}{2}}),\\\
y_{1}=a\cosh(\frac{\varphi
a}{\sqrt{2}}),&y_{2}=a\sinh\sqrt{3a}M,&y_{3}=\cosh(\sqrt{V(\varphi)}a^{\frac{3}{2}}),\end{array}$
(7)

and

$\begin{array}[]{c}x_{4}+y_{4}=1+\sqrt{V(\varphi)}a^{\frac{3}{2}},\\\
x_{4}-y_{4}=1-\sqrt{V(\varphi)}a^{\frac{3}{2}}.\end{array}$ (8)

We put the equations (7) and (8) into the equation (6) and obtain the
following Lagrangian,

$\mathcal{L}=\sum\limits_{i=1}^{4}((\dot{x}_{i}^{2}-\dot{y}_{i}^{2})+(x_{i}^{2}-y_{i}^{2})).$
(9)

By using the equation
$\mathcal{H}=P_{x_{i}}\dot{x}_{i}+P_{y_{i}}\dot{y}_{i}-\mathcal{L}$ one can
write the corresponding Hamiltonian as,

$\mathcal{H}=\frac{1}{4}\sum\limits_{i=1}^{4}((P_{x_{i}}^{2}-P_{y_{i}}^{2})+\omega_{i}^{2}(x_{i}^{2}-y_{i}^{2})),$
(10)

where

$\begin{array}[]{c}P_{x_{i}}=\frac{\partial\mathcal{L}}{\partial\dot{x_{i}}}=2\dot{x}_{i},\\\
P_{y_{i}}=\frac{\partial\mathcal{L}}{\partial\dot{y_{i}}}=-2\dot{y}_{i}.\end{array}$
(11)

It is obvious that the Hamiltonian (10) has a harmonic oscillator form with
the following frequency,

$\displaystyle\omega_{i}^{2}=-4.$ (12)

Here, we note that the Hamiltonian (10) is very useful to obtain the solution
of any cosmological system. Now, we are going to deform the non-minimal
coupling to the gravity system. In order to do such a process, we need to give
some transformations to the classical phase space. So, the first, we review
such transformation as an NC geometry. Therefore, before applying NC geometry
to the new Hamiltonian (10), we give some explanation of such geometry. As we
know, in commutative space we have the following Poisson brackets,

$\\{x_{i},x_{j}\\}=0,\quad\\{P_{x_{i}},P_{x_{j}}\\}=0,\quad\\{x_{i},P_{x_{j}}\\}=\delta_{ij},$
(13)

where $i=1,2$. But in NC geometry the Moyal product plays an important role in
the deformation phase space variables. Also, the quantum effects can be
dissolved by the following Moyal brackets,

$\\{f,g\\}_{\alpha}=f\star_{\alpha}g-g\star_{\alpha}f.$ (14)

The brackets (14) are based on the Moyal product which is given by,

$(f\star_{\alpha}g)(x)=exp[\frac{1}{2}\alpha^{ab}\partial_{a}^{(1)}\partial_{b}^{(2)}]f(x_{1})g(x_{2})|_{x_{1}=x_{2}=x}.$
(15)

By using the equations (14) and (15), one can find the following relation,

$\\{x_{i},x_{j}\\}_{\alpha}=\theta_{ij},\quad\\{x_{i},P_{j}\\}_{\alpha}=\delta_{ij}+\sigma_{ij},\quad\\{P_{i},P_{j}\\}=\beta_{ij}.$
(16)

So, in that case the transformations of the classical phase space variables
are given by,

$\hat{x}_{i}=x_{i}+\frac{\theta}{2}P_{y_{i}},\quad\hat{y}_{i}=y_{i}-\frac{\theta}{2}P_{x_{i}},\quad\hat{P}_{x_{i}}=P_{x_{i}}-\frac{\beta}{2}y_{i},\quad\hat{P}_{y_{i}}=P_{y_{i}}+\frac{\beta}{2}x_{i},\\\
$ (17)

and

$\\{\hat{y},\hat{x}\\}=\theta,\quad\\{\hat{x},\hat{P}_{x}\\}=\\{\hat{y},\hat{P}_{y}\\}=1+\sigma,\quad\\{\hat{P}_{y},\hat{P}_{x}\\}=\beta,$
(18)

where $\sigma=\frac{\beta\theta}{2}$. In order to construct the deformed
Hamiltonian for the non-minimal coupling to gravity model, we apply the new
variables (17) on the Hamiltonian (10). Hence, the Hamiltonian (10) in
deformed form will be as,

$\hat{\mathcal{H}}=\frac{1}{4}\sum_{i=1}^{4}((P_{x_{i}}^{2}-P_{y_{i}}^{2})-\gamma_{i}^{2}(y_{i}P_{x_{i}}+x_{i}P_{y_{i}})+{\tilde{\omega}_{i}^{2}}(x_{i}^{2}-y_{i}^{2})),$
(19)

where $\omega_{i}^{2}$ can be changed by the following equation,

${\tilde{\omega}_{i}}^{2}=\frac{{\omega_{i}}^{2}-\frac{\beta^{2}}{4}}{1-{\omega_{i}}^{2}\frac{\theta^{2}}{4}},\quad{\gamma_{i}}^{2}=\frac{\beta-{\omega_{i}}^{2}\theta}{1-{\omega_{i}}^{2}\frac{\theta^{2}}{4}}.$
(20)

Here, two parameters of NC geometry as $\theta$ and $\beta$ play important
role for the deformation of theory. When we investigate the bouncing universe,
we will explain the role of two parameters $\theta$ and $\beta$ with some
figures. The above deformed Hamiltonian lead us to have a new Lagrangian which
is deformed form of original Lagrangian of non-minimal coupling to gravity
model.
Now, for the corresponding model, we arrange the deformed Lagrangian as,

$\hat{\mathcal{L}}=\frac{1}{4}\sum_{i=1}^{4}((P_{x_{i}}^{2}-P_{y_{i}}^{2})-{\hat{\gamma}_{i}^{2}}(y_{i}P_{x_{i}}+x_{i}P_{y_{i}})-{\hat{\omega}_{i}^{2}}(x_{i}^{2}-y_{i}^{2})),$
(21)

where

${\hat{\omega}_{i}}^{2}=\frac{{\omega_{i}}^{2}+\frac{\beta^{2}}{4}}{1+{\omega_{i}}^{2}\frac{\theta^{2}}{4}},\quad{\hat{\gamma}_{i}}^{2}=\frac{\beta+{\omega_{i}}^{2}\theta}{1+{\omega_{i}}^{2}\frac{\theta^{2}}{4}}.$
(22)

We use the equations (11) and (21) to obtain the deformed Lagrangian with
respect to $a$, $\dot{a}$, $\varphi$, $\dot{\varphi}$ and $V(\varphi)$ which
is given by,

$\displaystyle\hat{\mathcal{L}}=-3M^{2}a\dot{a}^{2}+\frac{1}{2}(a\dot{a}^{2}\varphi^{2}+2a^{2}\dot{a}\varphi\dot{\varphi}+a^{3}\dot{\varphi}^{2})-\lambda_{1}^{2}V(\varphi)a^{3}+\lambda_{2}^{2}(\frac{\dot{\varphi}a^{3}}{\sqrt{2}}+\frac{a^{2}\dot{a}\varphi}{\sqrt{2}}-\sqrt{3}M\dot{a}a^{\frac{3}{2}}),$
(23)

where

$\lambda_{1}^{2}=\frac{1-\frac{\beta^{2}}{16}}{1+\theta^{2}},\quad\lambda_{2}^{2}=\frac{2\theta-\frac{1}{2}\beta}{1-\theta^{2}}.$
(24)

The above Lagrangian is deformed shape of original Lagrangian given by the
equation (6). So the above deformed Lagrangian shows us how the NC parameters
effect to the corresponding theory. We note here, for the non-minimal model
there were several papers for investigating the bouncing universe, and they
have shown that the equation of state cross $-1$. But, we use the deformed
non-minimal coupling to gravity model and obtain the equation of state. The
corresponding equation of state for the deformed case lead us to investigate
the bouncing universe.

##  IV Bouncing universe of deformed non-minimal coupling to gravity model

According to the conservation of energy,

$E_{\hat{\mathcal{L}}}=P_{a}\dot{a}+P_{\varphi}\dot{\varphi}-\hat{\mathcal{L}}$
(25)

where

$P_{a}=\frac{\partial\hat{\mathcal{L}}}{\partial\dot{a}},\quad
P_{\varphi}=\frac{\partial\hat{\mathcal{L}}}{\partial\dot{\varphi}}$ (26)

Here, by using definition of Hubble parameter $H=\frac{\dot{a}}{a}$, and
assumption $E_{\hat{\mathcal{L}}}\equiv 0$, we will arrive at,

$H^{2}=\frac{1}{6M^{2}}(H\varphi+\dot{\varphi})^{2}+\frac{\lambda_{1}^{2}}{3M^{2}}V(\varphi).$
(27)

The Euler-Lagrangian equation help us to obtain the following equations,

$\ddot{\varphi}+3H\dot{\varphi}+(2H^{2}+\dot{H})\varphi+\sqrt{2}\lambda_{2}^{2}H+\lambda_{1}^{2}V^{{}^{\prime}}(\varphi)=0,$
(28)

and

$\dot{H}=-\frac{2(H\varphi+\dot{\varphi})^{2}+\sqrt{2}\lambda_{2}^{2}(H\varphi+\dot{\varphi})+\lambda_{1}^{2}V^{{}^{\prime}}(\varphi)\varphi}{6M^{2}}.$
(29)

And then, by using $\dot{H}=-\frac{1}{2}(\rho+p)$, one can rewrite the
equation (29) as,

$p+\rho=\frac{\frac{2}{3}(H\varphi+\dot{\varphi})^{2}+\frac{\sqrt{2}}{3}\lambda_{2}^{2}(H\varphi+\dot{\varphi})+\frac{\lambda_{1}^{2}}{3}V^{{}^{\prime}}(\varphi)\varphi}{M^{2}}.$
(30)

Then, by using $p+\rho=(1+\omega)\rho$ and assuming $\dot{\varphi}=0$ and
$\omega=-1$, from the equation (24), and $\theta^{2}=\beta^{2}=0$, one can
obtain the following equation,

$\frac{2}{3}(H\varphi)^{2}+\frac{\sqrt{2}}{3}(2\theta-\frac{1}{2}\beta)(H\varphi)+\frac{1}{3}V^{{}^{\prime}}(\varphi)\varphi=0.$
(31)

The bouncing behavior of deformed of non-minimally coupled to gravity lead us
to arrange the suitable potential of system. For this reason, we need to check
the necessary condition for the bouncing universe, such condition help us to
obtain the corresponding potential. As we know, during the contracting phase,
the scale factor $a(t)$ is decreasing, it means $\dot{a}<0$, For the expanding
phase, we have $\dot{a}>0$. Also, we note here for a period of time at
bouncing point and around this point we have $\dot{a}=0$ and $\ddot{a}>0$
respectively. At the bouncing point, we have a transition from $H<0$ to $H\geq
0$ [49].

Figure 1: $H$ with respect to $t$, where $M=1$, $\varphi=0.5$, $\theta=2$.

In the Fig. 1 we show numerical solution of the equation (29) where $H<0$ to
$H>0$ transition illustrated. We can see modified parameters are necessary to
have bouncing universe. Green dotted line of the Fig. 1 show that at the
special time the Hubble expansion parameter change the sign from plus to
minus. In absence of modification, the Hubble expansion parameter is a
constant.
We can also find analytical expression, in such case the suitable bouncing
universe can obtained by,

$H=\frac{-\frac{\sqrt{2}}{2}\lambda_{2}^{2}\pm\sqrt{\frac{1}{2}\lambda_{2}^{4}-2\lambda_{1}^{2}V^{{}^{\prime}}(\varphi)\varphi}}{2\varphi}.$
(32)

In order to have bouncing universe we have to consider three cases, which help
us to obtain the exactly potential for the non-minimally coupled to gravity
model.
The first case: if we assume $H=0$, one can obtain
$V^{{}^{\prime}}(\varphi)=0$, therefore $V(\varphi)=constant$.
The second case: if the Hubble parameter $H$ is positive, then
$V^{{}^{\prime}}(\varphi)<0$ as a result $V(\varphi)\propto\varphi^{-n}$.
The third case: if $H$ is negative then there is no condition for potential,
therefore we suppose $V(\varphi)\propto\varphi^{-n}$. The pointed potential
gives us opportunity to study the cosmology evolution of equation of state for
the proposed model.
In order to explore the possibility of $\omega$ across -1, we have to check
out $\frac{d}{dt}(p+\rho)\neq 0$ when $\omega\neq-1$,

$\frac{d}{dt}(p+\rho)=\frac{4}{3}(H\varphi+\dot{\varphi})(\dot{H}\varphi+H\dot{\varphi}+\ddot{\varphi})+\frac{\sqrt{2}}{3}\lambda_{2}^{2}(\dot{H}\varphi+H\dot{\varphi}+\ddot{\varphi})+\frac{\lambda_{1}^{2}}{3}(V^{{}^{\prime\prime}}(\varphi)\dot{\varphi}\varphi+V^{{}^{\prime}}(\varphi)\dot{\varphi}),$
(33)

by assuming $\dot{\varphi}=\ddot{\varphi}=0$, one can rewrite following
equation,

$\frac{d}{dt}(p+\rho)=\frac{4}{3}(H\varphi)(\dot{H}\varphi)+\frac{\sqrt{2}}{3}\lambda_{2}^{2}(\dot{H}\varphi),$
(34)

therefore $\frac{d}{dt}(p+\rho)\neq 0$, otherwise
$H\varphi=-\frac{\sqrt{2}}{4}\lambda_{2}^{2}$ (it means that $H\varphi$ is
constant). This is not true because $H\varphi$ is function of time. So, the
above result proves the equation of state cross -1, therefore we have bouncing
universe. The equation of state given by,

$\omega=\frac{\frac{(H\varphi)^{2}}{6}+\frac{\sqrt{2}}{3}\lambda_{2}^{2}H\varphi-\frac{\lambda_{1}^{2}}{\varphi^{2}}-\lambda_{1}^{2}V(\varphi)}{\frac{(H\varphi)^{2}}{6}+\lambda_{1}^{2}V(\varphi)}.$
(35)

In the plots of the Fig. 1 we draw equation of state parameter in terms of
time to see effect of deformation parameters separately. In the left plot we
can see the effect of $\beta$ while in the right plot we can see the effect of
$\beta$.

Figure 2: The variation $\omega$ with respect to $t$, where
$V(\varphi)=\frac{V_{0}}{2\varphi^{2}}$, $\varphi=0.5$, $V_{0}=0.5$. (a)
$\theta=2$; (b) $\beta=0.1$

It is illustrated that the variation with $\beta$ is infinitesimal while
increasing $\theta$ increases value of the $\omega$ (see Fig. 2 (b)). We can
see from Fig. 2 that the equation of state parameter yields to a constant at
the late time. It is shown that in presence of the deformation parameters, the
equation of state parameter across -1 at earlier time. The values of
parameters fixed according to the stability condition which discussed below.
Here, also we can check the stability of system due to deformation of the non-
minimally coupled to gravity model. In order to examine the stability of
system we need to obtain the variation of pressure with respect to density
energy of system. Hence, we investigate the stability of non-minimal coupling
to gravity with deformation phase space variables. In that case we want to
consider the stability by useful functions as $C_{s}^{2}=\frac{dp}{d\rho}$. In
order to have stability with non-minimal deformation model such function must
become more than zero. So, if we want to investigate the stability of system
the above mentioned function can be expressed by sound speed in prefect liquid
[67]. By using the equations (27), (30) and $H^{2}=\frac{\rho}{3}$, one can
obtain the sound speed,

$\displaystyle\frac{dp}{d\rho}$ $\displaystyle=$
$\displaystyle\frac{\frac{1}{3}(H\varphi+\dot{\varphi})(\dot{H}\varphi+H\dot{\varphi}+\ddot{\varphi})+\frac{\sqrt{2}}{3}\lambda_{2}^{2}(\dot{H}\varphi+H\dot{\varphi}+\ddot{\varphi})}{(H\varphi+\dot{\varphi})(\dot{H}\varphi+H\dot{\varphi}+\ddot{\varphi})+\lambda_{1}^{2}V^{{}^{\prime}}(\varphi)\dot{\varphi}}$
(36) $\displaystyle+$
$\displaystyle\frac{\frac{\lambda_{1}^{2}}{3}(V^{{}^{\prime\prime}}(\varphi)\dot{\varphi}\varphi+V^{{}^{\prime}}(\varphi)\dot{\varphi})-\lambda_{1}^{2}V^{{}^{\prime}}(\varphi)\dot{\varphi}}{(H\varphi+\dot{\varphi})(\dot{H}\varphi+H\dot{\varphi}+\ddot{\varphi})+\lambda_{1}^{2}V^{{}^{\prime}}(\varphi)\dot{\varphi}}.$

If we assume $\dot{\varphi}=\ddot{\varphi}=0$, then $c_{s}^{2}\equiv C_{s}$
lead us to have following condition,

$\frac{1}{3}+\frac{\sqrt{2}}{3}\frac{\lambda_{2}^{2}}{H\varphi}>0.$ (37)

So, the above condition for the stability of non-minimally coupled to gravity
in case of deformed phase space lead us to have $4\theta>\beta$. Here, we draw
the variation of pressure with respect to the variation of energy density in
terms of time ($C_{s}=\frac{dp}{d\rho}$) and show that how the stability
correspond to two parameters $\theta$ and $\beta$ from NC geometry. It is
illustrated by the Fig. 3 that non-deformed case is completely stable as well
as the deformed case with condition $\beta=4\theta$ (see solid red line of the
Fig . 3).

Figure 3: The variation $\frac{dp}{d\rho}$ with respect to $t$, where
$V(\varphi)=\frac{V_{0}}{2\varphi^{2}}$, $\varphi=0.5$, $V_{0}=0.5$,
$\theta=0.5$.

##  V Conclusion

In this paper, we introduced non-minimally coupled gravity and wrote
Lagrangian. Also, by using some suitable variables, we obtained new
Hamiltonian and Lagrangian. In order to deform the above mentioned system we
arranged the harmonic oscillator form of Hamiltonian. We used the NC geometry
and deformed the new Hamiltonian and Lagrangian. Now, we take advantage from
dark energy solution and investigated bouncing universe in deformed non-
minimal coupling to gravity model. The condition of bouncing universe in
deformed model lead us to arrange the potential in equation of state. In order
to have bouncing solution we drew the equation of state in terms of time. It
shown that the equation of state crossing $-1$. Also, we obtained some
suitable relation between the $\theta$ and $\beta$ parameters of NC geometry
to satisfy by bouncing universe. We found special condition ($\beta=4\theta$)
where the stability of deformed model behaves like non-deformed case. In the
case of $\beta\leq 4\theta$ the deformed model is completely stable. Otherwise
the model is initially stable which yields to unstable phase at the late time.
Hence it is possible to have stable/unstable phase transition. In order to
have solution of dark energy, we arranged the potential in the equation of
state of deformed non-minimally coupled to gravity. In that case we draw the
equation of state in terms of time and shown that the equation of state cross
$-1$. Finally, we employed the pressure and energy density and investigated
the stability of deformed non-minimally coupled to gravity. Here, we obtained
the relation between $\theta$ and $\beta$ in NC for the approving the
stability of system. Also here, we had some figure for the describing of
stability with some special values of parameters of NC geometry. The relation
between the above deformation of corresponding model and Finsler geometry will
be important in future for the describing some parameter in physics. Also,
such relation give us more degree of freedom to compare the obtained results
of theory with data. As we know the results from Finsler geometry have a
relation with non-commutative geometry. So, we hope to continue this research
for the different cosmological model in future. For instance, recently
bouncing universe considered in the contexts of generalized cosmic Chaplygin
gas and variable modified Chaplygin gas, now it is interesting to consider
bouncing universe in the context of extended Chaplygin gas [68, 69, 70, 71,
72] and use method of this paper. It is also possible to use the same method
of this paper to study the holographic Barrow dark energy model [73, 74].

## References

  * [1] S. Nojiri, S. D. Odintsov and E. N. Saridakis, Holographic bounce, Nucl. Phys. B 949, 114790 (2019)
  * [2] X. Zhang and F. Q. Wu, Constraints on holographic dark energy from type Ia supernova observations, Phys. Rev. D 72, 043524 (2005).
  * [3] M. Li, X. D. Li, S. Wang and X. Zhang, Constraints on holographic dark energy from type Ia supernova observations, JCAP 0906, 036 (2009)
  * [4] J. Lu, E. N. Saridakis, M. R. Setare and L. Xu, Observational constraints on holographic dark energy with varying gravitational constant, JCAP 1003, 031 (2010)
  * [5] Q. G. Huang and Y. G. Gong, Supernova Constraints on a holographic dark energy model, JCAP 0408 (2004) 006 [astro-ph/0403590].
  * [6] B. Wang, E. Abdalla and R. K. Su, Constraints on the dark energy from holography,  Phys. Lett. B 611, 21 (2005)
  * [7] M. Novello and S. E. Perez Bergliaffa, Bouncing Cosmologies, Phys. Rept. 463, 127 (2008)
  * [8] S. Mukherji and M. Peloso, Bouncing and cyclic universes from brane models, Phys. Lett. B 547, 297 (2002)
  * [9] D. N. Vollick, Negative Energies on the Brane, Gen. Relativ. Gravit. 34, 1 (2002)
  * [10] A. Biswas, S. Mukherji and S. S. Pal, Nonsingluar cosmologies from branes, Int. J. Mod. Phys. A 19, 557 (2004)
  * [11] D. H. Coule, Does brane cosmology have realistic principles?, Class. Quantum Grav. 18, 4265 (2001)
  * [12] S. Foffa, Bouncing pre-big bang on the brane, Phys. Rev. D 68, 043511 (2003)
  * [13] J. Hovdebo and R. Myers, Bouncing braneworlds go crunch!, JCAP 011, 012 (2003)
  * [14] P. Kanti and K. Tamvakis, Challenges and obstacles for a bouncing universe in brane models, Phys. Rev. D 68, 024014 (2003)
  * [15] A. J. Medved, Of Bounces, Branes and Bounds JHEP 0305, 008 (2003)
  * [16] Y. S. Myung, Bouncing and cyclic universes in the charged AdS bulk background, Class. Quantum Grav. 20, 935 (2003)
  * [17] P. Astier and et al., The Supernova Legacy Survey: measurement of $\Omega_{M}$, $Omega_{\Lambda}$ and $w$ from the first year data set,  Astron. Astrophys. 447, 31 (2006)
  * [18] S. Perlmutter and et al., Measurements of $\Omega$ and $\Lambda$ from 42 High-Redshift Supernovae (Supernova Cosmology Project Collaboration), Astrophys. J. 517, 565 (1999)
  * [19] A. G. Riess and et al., Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant (Supernova Search Team Collaboration), Astron. J. 116, 1009 (1998)
  * [20] K. Abazajian and et al., The Second Data Release of the Sloan Digital Sky Survey (SDSS Collaboration), Astron. J. 128, 502 (2004)
  * [21] K. Abazajian and et al., The Third Data Release of the Sloan Digital Sky Survey (SDSS Collaboration), Astron. J. 129, 1755 (2005)
  * [22] D. N. Spergel and et al., First-Year Wilkinson Microwave Anisotropy Probe (WMAP)* Observations: Determination of Cosmological Parameters (WMAP Collaboration), Astrophys. J. Suppl. 148, 175 (2003)
  * [23] D. N. Spergel and et al., Three-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Implications for Cosmology, Astrophys. J. Suppl. 170, 377 (2007)
  * [24] R. R. Caldwell, R. Dave and P. J. Steinhardt, Cosmological Imprint of an Energy Component with General Equation of State, Phys. Rev. Lett. 80, 1582 (1998)
  * [25] P. J. Steinhardt, L. M. Wang and I. Zlatev, Cosmological tracking solutions, Phys. Rev. D  59, 123504 (1999)
  * [26] C. Wetterich, Cosmology and the fate of dilatation symmetry, Nucl. Phys. B 302, 668 (1988)
  * [27] C. Armendariz-Picon, V. F. Mukhanov and P. J. Steinhardt, Essentials of $k$-essence, Phys. Rev. D 63 103510 (2001)
  * [28] T. Padmanabhan, Accelerated expansion of the universe driven by tachyonic matter, Phys. Rev. D 66, 021301 (2002)
  * [29] A. Sen, Tachyon Matter, JHEP 0207, 065 (2002)
  * [30] N. Arkani-Hamed, H. C. Cheng, M. A. Luty and S. Mukohyama, Ghost Condensation and a Consistent Infrared Modification of Gravity, JHEP 0405, 074 (2004)
  * [31] F. Piazza and S. Tsujikawa, Dilatonic ghost condensate as dark energy, JCAP 0407, 004 (2004)
  * [32] P. S. Apostolopoulos and N. Tetradis, Late acceleration and $w=-1$ crossing in induced gravity, Phys. Rev. D 74, 064021 (2006)
  * [33] B. Feng, X. Wang and X. Zhang, Dark Energy Constraints from the Cosmic Age and Supernova, Phys. Lett. B 607, 35 (2005)
  * [34] B. Feng, M. Li, Y. Piao and X. Zhang, Oscillating Quintom and the Recurrent Universe, Phys. Lett. B 634, 101 (2006)
  * [35] Z. Guo, Y. Piao, X. Zhang and Y. Z. Zhang, Cosmological Evolution of a Quintom Model of Dark Energy, Phys. Lett. B 608, 177 (2005)
  * [36] M. Li, B. Feng and X. Zhang, A Single Scalar Field Model of Dark Energy with Equation of State Crossing -1, JCAP 0512, 002 (2005)
  * [37] X. F. Zhang and T. T. Qiu, Avoiding the Big-Rip Jeopardy in a Quintom Dark Energy Model with Higher Derivatives, Phys. Lett. B 642, 187 (2006)
  * [38] H. S. Zhang and Z. H. Zhu, Crossing $w=-1$ by a single scalar on a Dvali-Gabadadze-Porrati brane, Phys. Rev. D 75, 023510 (2007)
  * [39] X. F. Zhang, H. Li, Y. S. Piao and X. M. Zhang, Two-field Models of Dark Energy with Equation of State Across -1, Mod. Phys. Lett. A 21, 231 (2006)
  * [40] A. Y. Kamenshchik, U. Moschella and V. Pasquier, Chaplygin-like gas and branes in black hole bulks, Phys. Lett. B 487, 7 (2000)
  * [41] A. Yu. Kamenshchik, U. Moschella and V. Pasquier, An alternative to quintessence, Phys. Lett. B 511, 265 (2001)
  * [42] H. Saadat and B. Pourhassan, FRW Bulk Viscous Cosmology with Modified Chaplygin Gas in Flat Space, Astro. Space Science 343, 783 (2013)
  * [43] H. Saadat and B. Pourhassan, FRW bulk viscous cosmology with modified cosmic Chaplygin gas, Astro. Space Science 344, 237 (2013)
  * [44] N. Bilic, G. B. Tupper and R. D. Viollier, Unification of Dark Matter and Dark Energy: the Inhomogeneous Chaplygin Gas, Phys. Lett. B 535 17 (2002)
  * [45] M. C. Bento, O. Bertolami and A. A. Sen, Generalized Chaplygin gas, accelerated expansion, and dark-energy-matter unification, Phys. Rev. D 66, 043507 (2002)
  * [46] J. C. Fabris, S. V. B. Goncalves and P. E. de Souza, Density perturbations in an Universe dominated by the Chaplygin gas, Gen. Relativ. Gravit. 34, 53 (2002)
  * [47] B. Pourhassan, Viscous Modified Cosmic Chaplygin Gas Cosmology, Int. J. Mod. Phys. D 22, 1350061 (2013)
  * [48] E. O. Kahya, B. Pourhassan and S. Uraz, Constructing an inflaton potential by mimicking modified Chaplygin gas, Phys. Rev. D  92, 103511 (2015)
  * [49] B. Pourhassan, Unified universe history through phantom extended Chaplygin gas, Can. J. Phys. 94, 659 (2016)
  * [50] L. Amendola, Coupled quintessence, Phys. Rev. D 62, 043511 (2000)
  * [51] M. Szydlowski, A. Kurek and A. Krawiec, Top ten accelerating cosmological models, Phys. Lett. B 642, 171 (2006)
  * [52] C. Deffayet, G. R. Dvali and G. Gabadadze, Accelerated universe from gravity leaking to extra dimensions, Phys. Rev. D 65, 044023 (2002)
  * [53] M. R. Setare, Bulk–brane interaction and holographic dark energy, Phys. Lett. B 642, 421 (2006)
  * [54] A. G. Cohen, D. B. Kaplan and A. E. Nelson, Effective Field Theory, Black Holes, and the Cosmological Constant, Phys. Rev. Lett. 82, 4971 (1999)
  * [55] P. Horava and D. Minic, Probable Values of the Cosmological Constant in a Holographic Theory, Phys. Rev. Lett. 85, 1610 (2000)
  * [56] E. Elizalde, S. Nojiri, S. D. Odintsov and P. Wang, Dark Energy: Vacuum Fluctuations, the Effective Phantom Phase, and Holography, Phys. Rev. D 71, 103504 (2005)
  * [57] M. R. Setare, J. Zhang and X. Zhang, Statefinder diagnosis in a non-flat universe and the holographic model of dark energy, JCAP 0703, 007 (2007)
  * [58] R. R. Caldwell and M. Doran, Dark-energy evolution across the cosmological-constant boundary, Phys. Rev. D 72, 043527 (2005)
  * [59] W. Hu, Crossing the phantom divide: Dark energy internal degrees of freedom, Phys. Rev. D 71, 047301 (2005)
  * [60] A. Vikman, Can dark energy evolve to the phantom?, Phys. Rev. D 71, 023515 (2005)
  * [61] G. B. Zhao, J. Q. Xia, M. Li, B. Feng and X. Zhang, Perturbations of the Quintom Models of Dark Energy and the Effects on Observations, Phys. Rev. D 72, 123515 (2005)
  * [62] J. Sadeghi, M. R. Setare and A. Banijamali, String inspired quintom model with non-minimally coupled modified gravity, Phys. Lett. B 678, 164 (2009)
  * [63] M. R. Setare, J. Sadeghi and A. Banijamali, Bouncing Universe with Non-minimally Coupled Quintom Matter,  Phys. Lett. B 669, 9 (2008)
  * [64] P. A. R. Ade and et al., Planck 2015 results - XX. Constraints on inflation, Astron. Astrophys. 594, A20 (2016)
  * [65] T. Markkanen, T. Tenkanen, V. Vaskonen and H. Veermäe, Quantum corrections to quartic inflation with a non-minimal coupling: metric vs. Palatini, JCAP 03, 029 (2018)
  * [66] J. Sadeghi, B. Pourhassan, Z. Nekouee and M. Shokri, Deformation of the quintom cosmological model and its consequences, Int. J. Mod. Phys. D  27, 1850025 (2018)
  * [67] J. Sadeghi, F. Milani and A. R. Amani, Bouncing universe with the non-minimally coupled scaler field and its reconstructing, Mod. Phys. Lett. A 24, 2363 (2009)
  * [68] E. O Kahya, B. Pourhassan, E. O Kahya and B. Pourhassan, Modern Physics Letters A 30, 1550070 (2015), Mod. Phys. Lett. A 30, 1550070 (2015)
  * [69] B. Pourhassan, Extended Chaplygin Gas in Horava-Lifshitz Gravity, Physics of the Dark Universe 13, 132 (2016)
  * [70] E.O. Kahya, M. Khurshudyan, B. Pourhassan, R. Myrzakulov and A. Pasqua, Higher order corrections of the extended Chaplygin gas cosmology with varying $G$ and $\Lambda$, Eur. Phys. J. C 75, 43 (2015)
  * [71] J. Sadeghi, H. Farahani and B. Pourhassan, Interacting Holographic Extended Chaplygin Gas and Phantom Cosmology in the Light of BICEP2, Eur. Phys. J. Plus 130, 84 (2015)
  * [72] E.O. Kahya and B. Pourhassan, Observational constraints on the extended Chaplygin gas inflation, Astro. Space Science 353, 677 (2014)
  * [73] E. N. Saridakis, Barrow holographic dark energy, Phys. Rev. D 102, 123525 (2020).
  * [74] S. Maity and P. Rudra, Inflation driven by Barrow holographic dark energy, Journal of Holography Applications in Physics, 2(1), 1 (2022)