text:
\section{Introduction}

\textit{\textbf{Aspects of Nuclear Matter:-}}
Nuclear matter is a hypothetical, infinitely large system of nucleons ($N=Z$) with coulomb interaction turned off. The matter inside the compact stars, the early universe, a nucleon gas,
or quark-gluon plasma ({\bf QGP}) are some of the types of nuclear matter. The only physical parameters that remains to characterize such a system is the bulk binding energy and the saturation density, letting the total number of nucleons go to infinity.

\begin{figure}[ht]
\vskip -0.2in
\begin{center}
\includegraphics[width=7cm,angle=0]{phase1.eps}
\vskip -0.2in
\caption{Theoretical Phase diagram for nuclear matter.}
\end{center}
\end{figure}

Fig. 1 shows the conventional picture of different aspects and phase of nuclear matter with the variation of baryon density ($\rho_B$) and temperature (T). At low temperature and high density, roughly of the order of ten times normal nuclear density, we approach the matter configuration which may be similar to that present in the core of neutron stars. The discovery of neutron stars in the form of pulsars has been a major stimulus to dense matter studies, which serves as an ideal laboratory to study matter at extreme conditions. They are the most compact objects observed till date, confined to a radius of nearly 10 km with the density ranging from $(3~-~10)~\rho_0$ ($\rho_0$ is the normal nuclear matter density) in the core of these stars. Bound by gravity, the evolution and constitution of these stars represents a beautiful amalgamation of all the known forces in nature.

Theoretically, the framework of Quantum Hadrodynamics \cite{wal74,serot86} present an elegant and consistent treatment of finite nuclei as well as infinite nuclear matter, which seem to provide solution to the so called ``{\it the Coester band}" problem \cite{co01}. A realistic nuclear equation of state must satisfy certain minimum criteria quantified as the {\it ``nuclear saturation properties''}, which are the physical constants of nature. In this context, the two most important quantities which play vital role and are known to have substantial impact on the EOS are the nucleon effective mass and the nuclear incompressibility \cite{compact}, which are not very well determined and posses large uncertainty. Nucleon effective mass or the medium mass modification of nucleon in nuclear medium is a consequence of the Dirac field and forms an essential element for the success of relativistic phenomenology. On the other hand, the nuclear incompressibility derived from nuclear measurements and astrophysical observations exhibit a broad range of values $K = (180 - 800)$ MeV \cite{glen88}. Further the non-relativistic and the relativistic models fails to agree to a common consensus. The non-relativistic calculations predict the compression moduli in the range $K=(210-240) MeV$ \cite{k2}, whereas, relativistic calculations predicts it in the range $(200-300)~MeV$ \cite{nl3,k6}. Apart from that we are inevitably marred by the uncertainty in the determination of mass of the scalar meson ($\sigma$-meson). The attractive force resulting from the scalar sector is responsible for the intermediate range attraction which, along with the repulsive vector forces provides the saturation mechanism for nuclear matter \cite{serot86}. The estimate from the Particle Data Group quotes the mass of this scalar meson `$f_0 (600)$' or $\sigma-$meson in the range $(400 - 1200)$ MeV \cite{pdg}. A recent estimate however, for sigma meson mass is found to be $513 \pm 32$ MeV \cite{mura02}. Therefore in order to address these issues, one needs a model that has the desired attributes of the relativistic framework and which can be successfully applied to various nuclear force problem both in the vicinity of $\rho_0$ as well as at higher densities.

\vskip 0.1in
\textit{\textbf{Neutron Star Structure:-}} Matter at extremes opens up numerous possibilities, which may be present in neutron star core as shown in Fig. 2. Possible scenarios such as the transition from nuclear and hadronic matter, to exotic states involving pionic and  Kaonic Bose-Einstein condensation, to bulk quark matter and quark matter in droplets, including super-conducting states, as well as strange quark matter, have been proposed \cite{compact,weber,QCD}. Hyperons can form in the neutron star cores, when the nucleon chemical potential are large enough to compensate the mass difference between the nucleon and the hyperon. Also, in dense matter, the energy of the pion ($\pi$ meson) and the kaon ($K$) is modified by inter-particle interactions and if it becomes sufficiently low, they will form the condensate, i.e., being bosons, they will accumulate the same ground state \cite{condense}. The presence of the condensates is also known to enhance the neutrino cooling \cite{cool}, in addition to the nucleon direct URCA (Unrecordable Cooling Agents) process \cite{urca}. It is also widely believed that the attractive interactions between the quarks will lead to pairing and color superconductivity. For three massless flavors, the condensation pattern that minimizes the free energy is known as the color-flavor-locked (CFL) scheme \cite{cfl}. A detailed description of these novel phases can be found in Ref. \cite{sanjay, strange1}.

\begin{figure}\centering
\includegraphics[width=5cm,angle=0]{ns-phase1.eps}
\caption{Possible phases structure in neutron star core.}
\end{figure}

\vskip 0.1in
\textit{\textbf{Observations from Neutron star:-}}
Observations from neutron stars can lead to an understanding of the state of their interiors and the key unknowns such as the stars maximum mass and radius. Apart from that, the other key observables are the gravitational redshift, the central density of the star, the moments of inertia and the pulsar timings. These informations can be inferred from the photons, ranging from radio waves to X-rays, and also those involving neutrinos and gravity waves. The largest precisely known neutron star mass is only $1.44~ M_{\odot}$, but in case of binaries with white dwarf companions, the mass might be somewhere in the vicinity of 2 $M_{\odot}$ or even larger. The radius of the star is found to be sensitive to the properties of the nuclear matter near the saturation density ($\rho_0$), such as the density dependence of the nuclear symmetry energy \cite{lat01}. However, the radius measurements are not that precisely known unlike the mass observations, but the upper limits of the radius values are inferred from the thermal emission of cooling neutron stars, the gravitational redshift measurements, or from the crustal properties such as the pulsar glitches, star quakes and cooling timescales. The vital information from these observations can be used to constrain the nuclear equation of state at high densities.

\begin{figure}
\begin{center}
\includegraphics[width=7cm,angle=0]{ns-mass1.eps}
\caption {Observed mass of neutron stars in binaries. The figure is taken from Ref \cite{sanjay}. ($PSR+NS$) is the binary of a pulsar and a neutron star and ($PSR+WD$) is the binary of a pulsar with a white dwarf companion. HMXBs' are the high-mass X-ray binaries and `LMXBs' are the low-mass X-ray binaries. The two objects plotted as crosses, the HMXB 4U $1700-37$ and the LMXB 2S $0921-630$ are blackhole candidates, but may be neutron stars. Error bars are $1\sigma$ errors.}
\end{center}
\end{figure}

The measured mass of 40 neutron stars in binary systems are summarized in Fig. 3. In case of `$PSR J0737-3039$', for the A and B system, both neutron stars are detected as radio pulsars, whereas in the three systems with the companion marked by `?' may be a white dwarf. In a high-mass X-ray binary (HMXB), the companion is a massive star $M > 10M_{\odot}$, whereas in a low-mass X-ray binary (LMXB), the companion's mass is below $1M_{\odot}$. The most accurately measured masses are from timing observations of the radio binary pulsars such as in case of binary pulsar PSR 1913+16, in which the masses are $M=1.3867 \pm 0.0002 M_{\odot}$ and $M=1.4414 \pm 0.0002 M_{\odot}$, respectively. It is found that the mass determinations in binaries with white dwarf companions shows a broader range of neutron star masses than binary neutron star pulsars. The neutron star, PSR J0751+1807 with mass 2.1 $M_{\odot}$, is about 4$\sigma$ from the canonical value of 1.4 $M_{\odot}$. The system of $Vela X-1$ however, has lower mass limit (1.6 to 1.7$M_{\odot}$).

{\it Observed radii:} From the study of the measured flux of the nearest neutron star $RX J1856.5-3754$, with the Hubble space telescope, Kaplan et.al., estimated the stars radius
\cite{ob5}. The analysis puts the radiation radius of the star $R_{\infty}$ = $\frac{R}{\sqrt{1-2GM/R}} \simeq 15$km corresponding to $R\simeq12$ km for $M=1.4M_{\odot}$. Recently, a few surveys from the non-accreting X-ray binaries, the so called the Quiescent stars
$CXOU132619.7-472910.8$ in $NGC 5139$ gives $R_{\infty}=14.3 \pm 2.1$km (90\% CL) \cite{ob6}. In case of the neutron star $4U1636-53$, the compactness parameter `$M/R$' $<0.16$ implies a radius $R>12-13$km for the star with $M=1.4M_{\odot}$.

Some other observables from compact stars are the gravitational redshift measurement of the star, the moment of inertia or the estimated central density of the star. Overall, there are large uncertainties in these estimates. Further stringent conditions can be imposed on the EOS of dense matter from the laboratory estimates of the heavy-ion data, such as constraints on nuclear incompressibility, symmetry energy or nucleon effective mass.

\section{Effective chiral Model}

Chiral models \cite{ch01,ch02} have been developed and were applied to nuclear matter studies. We start with an effective chiral Lagrangian and calculate the density effects in nuclear matter in the mean field approach. To have a realistic description of the dense neutron
star matter, we consider the effective Lagrangian of the chiral model generalized
to include the lowest lying octet of baryons ($n,p,\Lambda^{0},\Sigma^{-,0,+},\Xi^{-,0}$)
interacting through the exchange of the pseudo-scalar meson $\pi$, the scalar
meson $\sigma$, the vector meson $\omega$ and the iso-vector $\rho-$meson,
and is given by \cite{tkj06}:

\begin{eqnarray}
\label{lag}
{\cal L}&=& \bar\psi_B~\left[ \big(i\gamma_\mu\partial^\mu
         - g_{\omega B}\gamma_\mu\omega^\mu
         - \frac{1}{2}g_{\rho B}{\vec \rho}_\mu\cdot{\vec \tau}
            \gamma^\mu\big ) \right]~ \psi_B \nonumber \\
         &-& \bar\psi_B~\left[g_{\sigma B~}~\big(\sigma + i\gamma_5
             \vec \tau\cdot\vec \pi \big)\right]~ \psi_B
        + \frac{1}{2}\big(\partial_\mu\vec \pi\cdot\partial^\mu\vec\pi
\nonumber \\
&&
        + \partial_{\mu} \sigma \partial^{\mu} \sigma\big)
        - \frac{\lambda}{4}\big(x^2 - x^2_0\big)^2
        - \frac{\lambda B}{6}\big(x^2 - x^2_0\big)^3 \nonumber \\
        &-& \frac{\lambda C}{8}\big(x^2 - x^2_0\big)^4
      - \frac{1}{4} F_{\mu\nu} F_{\mu\nu}
        + \frac{1}{2}{g_{\omega B}}^{2}x^2 \omega_{\mu}\omega^{\mu} \nonumber \\
        &-& \frac {1}{4}{\vec R}_{\mu\nu}\cdot{\vec R}^{\mu\nu}
        + \frac{1}{2}m^2_{\rho}{\vec \rho}_{\mu}\cdot{\vec \rho}^{\mu}\ .
\end{eqnarray}

The above Lagrangian represents the interaction of baryons $\Psi_B$ with the aforesaid mesons. We have the kinetic and the non-linear terms in the pseudoscalar-isovector pion field `$\vec \pi$', the scalar field `$\sigma$', and with $x^2= {\vec \pi}^2+\sigma^{2}$. Finally in the last two lines, we have the field strength and the mass term for the vector field `$\omega$' and the iso-vector field `$\vec \rho$' meson. The terms in eqn. (1) with the subscript $'B'$ should be interpreted as sum over the states of all baryonic octets. In this paper we shall be concerned only with the normal non-pion condensed state of matter, so we take $<\vec \pi>=0$.
The interaction of the scalar and the pseudoscalar mesons with the vector boson generates a dynamical mass for the vector bosons through spontaneous breaking of the chiral symmetry with scalar field getting the vacuum expectation value $x_0$. Then the masses of the baryons, the scalar and the vector mesons, are respectively given by

\begin{eqnarray}
m_B = g_{\sigma B} x_0,~~ m_{\sigma} = \sqrt{2\lambda} x_0,~~
m_{\omega} = g_{\omega B} x_0\ .
\end{eqnarray}

Using the meson field equations for $\omega$ and $\sigma$-meson along with the corresponding energy density and pressure expression for symmetric nuclear matter, the parameters of the present model are then evaluated. Here we fix the saturation density ($\rho_0 = 0.153~fm^{-3}$) and for a desired value of nucleon effective mass ($Y~=~m^{\star}/m~=~(0.75 - 0.90)$), we obtain the nuclear matter parameters of the model enlisted in Table I. For details of the procedure involved in obtaining the parameters, one can refer to \cite{tkj08}.

We now go directly to the total energy density `$\varepsilon$' and pressure `$P$'
for a given baryon density in terms of the dimensionless variable $Y=x/x_0$
which is given as:

\begin{eqnarray}
\label{ep0}
\varepsilon
&=&
 \frac{2}{\pi^2}\int^{k_B}_0 k^{2}dk{\sqrt{k^2+m_B^{\star 2}}}
         +  \frac{m_B^2(1-Y^2)^2}{8c_{\sigma B}} \nonumber \\
        &-& \frac{m_B^2 B}{12c_{\omega B}c_{\sigma B}}(1-Y^2)^3
        + \frac{m_B^2 C}{16c_{\omega B}^2c_{\sigma B}}(1-Y^2)^4 \nonumber \\
        &+& \frac{1}{2Y^2}{c_{\omega_{B}} \rho_B^2}
        +\frac{1}{2}m_{\rho}^{2}\rho_{03}^2 \nonumber \\
        &+& \frac{1}{\pi^{2}}\sum_{\lambda=e,\mu^{-}}\int^
           {k_\lambda}_0 k^{2}dk{\sqrt{k^2+m^2_{\lambda}}}\ ,
\end{eqnarray}
\begin{eqnarray}
P &=&
          \frac{2}{3\pi^2}\int^{k_B}_0\frac{k^{4}dk}
          {{\sqrt{k^2+m_B^{\star 2}}}}
        - \frac{m_B^2(1-Y^2)^2}{8c_{\sigma B}} \nonumber \\
        &+& \frac{m_B^2 B}{12c_{\omega B}c_{\sigma B}}(1-Y^2)^3
        - \frac{m_B^2 C}{16c_{\omega B}^2c_{\sigma B}}(1-Y^2)^4 \nonumber \\
        &+& \frac{1}{2Y^2}{c_{\omega_{B}} \rho_B^2}
        + \frac{1}{2}m_{\rho}^{2}\rho_{03}^2\ \nonumber \\
        &+& \frac{1}{3\pi^2}\sum_{\lambda=e,\mu^{-}}\int^{k_\lambda}_0
          \frac{k^{4}dk}{{\sqrt{k^2+m^2_{\lambda}}}}
\end{eqnarray}

The terms in eqns. (3) and (4) with the subscript $`B'$ should be interpreted
as sum over all the states of the baryonic octets.
The meson field equations for the $\sigma$, $\omega$ and $\rho-$mesons are then solved
self-consistently at a fixed baryon density to obtain the respective field strengths.
The EOS for the $\beta-$equilibrated for the hyperon rich matter is obtained
with the requirements of conservation of total baryon number and charge
neutrality condition \cite{tkj06}. Using the computed EOS for the neutron star sequences, we calculate the structural properties of neutron stars with and without hyperon core.

\section{Results}

Nuclear matter saturation is a consequence of the interplay between the attractive (scalar) and the repulsive (vector) forces and hence the variation in the coupling strength effects other related properties as well. Fig. \ref{effm}(A) reflects the same, where we have plotted the nuclear incompressibility for the evaluated parameter sets of the present model as a function of the nucleon effective mass. For better correlation between them, the corresponding ratio of the scalar and vector coupling is also indicated. On comparison with the incompressibility bound inferred from heavy ion collision experiment (HIC)\cite{data02}, we find that the EOS with lower nucleon effective mass is ruled out. The model favors EOS for which the nucleon effective mass $m^{\star}/m > 0.82$. It can also be seen that the EOS becomes much softer with increasing ratio of $C_{\sigma}/C_{\omega}$.

\begin{table*}
\caption{Parameter sets of the effective chiral model that satisfies the nuclear matter saturation properties such as binding energy per nucleon $B/A - m = -16.3 ~MeV$, nucleon effective mass $Y = m^{\star}/m = (0.75 - 0.90)$ and the asymmetry energy coefficient is $J \approx 32$ MeV at saturation density $\rho_0$ $=0.153 fm^{-3}$. The nucleon, the vector mesons ($\omega$ \& $\rho_0^3$) masses are taken to be 939 MeV, 783 MeV and 770 MeV respectively and $c_{\sigma} = (g_{\sigma}/ m_{\sigma})^2$, $c_{\omega} = (g_{\omega}/ m_{\omega})^2$ and $c_{\rho} = (g_{\rho}/ m_{\rho})^2$ are the corresponding coupling constants. $B = b/m^2$ and $C = c/ m^4$ are the higher order constants in the scalar field. Also given is the scalar meson mass `$m_{\sigma}$', the pion decay constant `$f_{\pi}$' and the nuclear matter incompressibility ($K$) at $\rho_0$. The maximum mass and radius of neutron star composed of nucleon only matter ($M_N, ~~R_N$) and hyperon rich matter $(M_H,~~~R_H)$ for static case are also tabulated. The last column shows the central density of the star composed of nucleon ($\varepsilon_N$) and hyperon ($\varepsilon_H$).}
\begin{center}
\begin{tabular}{cccccccccccccccc}
\hline
\hline
\multicolumn{1}{c}{set}&
\multicolumn{1}{c}{$c_{\sigma}$}&
\multicolumn{1}{c}{$c_{\omega}$} &
\multicolumn{1}{c}{$c_{\rho}$} &
\multicolumn{1}{c}{$B$} &
\multicolumn{1}{c}{$C$} &
\multicolumn{1}{c}{$m_{\sigma}$} &
\multicolumn{1}{c}{$Y$} &
\multicolumn{1}{c}{$f_{\pi}$} &
\multicolumn{1}{c}{$K$} &
\multicolumn{1}{c}{$M_N$($M_H$)} &
\multicolumn{1}{c}{$R_N~ (R_H)$} &
\multicolumn{1}{c}{$\varepsilon_{N}~(\varepsilon_{H})$} \\
\multicolumn{1}{c}{ } &
\multicolumn{1}{c}{($fm^2$)} &
\multicolumn{1}{c}{($fm^2$)} &
\multicolumn{1}{c}{($fm^2$)} &
\multicolumn{1}{c}{($fm^2$)} &
\multicolumn{1}{c}{($fm^4$)}&
\multicolumn{1}{c}{(MeV)} &
\multicolumn{1}{c}{ } &
\multicolumn{1}{c}{(MeV)} &
\multicolumn{1}{c}{($MeV$)} &
\multicolumn{1}{c}{$(M/ M_{\odot})$} &
\multicolumn{1}{c}{$(km)$} &
\multicolumn{1}{c}{$(10^{15} g cm^{3})$} \\
\hline
1   &5.916  &3.207  &5.060  &1.411   &1.328    &691.379  &0.75  &110.185  &1098 &2.49 (2.04)& 14.2 (20.8) & 1.26 (0.59) \\
2   &6.047  &3.126  &5.087  &0.822   &0.022    &675.166  &0.76  &111.601  &916 &2.47 (2.02)& 14.2 (20.7) & 1.26 (0.62) \\
3   &6.086  &3.031  &5.107  &0.485   &0.174    &662.642  &0.77  &113.346  &809 &2.44 (1.98)& 14.0 (20.5)& 1.30 (0.62) \\
4   &6.005  &2.933  &5.131  &0.582   &2.650    &656.183  &0.78  &115.238  &737 &2.41 (1.95)& 13.9 (20.4) & 1.35 (0.62) \\
5   &6.172  &2.825  &5.155 &-0.261   &0.606    &635.287  &0.79  &117.403  &638 &2.38 (1.92)& 13.7 (20.2) & 1.40 (0.65) \\
6   &6.223  &2.709  &5.178 &-0.711   &0.748    &619.585  &0.80  &119.890  &560 &2.34 (1.89)& 13.6 (20.0) & 1.44 (0.65)\\
7   &6.325  &2.585  &5.200 &-1.381   &0.089    &600.270  &0.81  &122.740  &491 &2.30 (1.85)& 13.4 (19.8) & 1.49 (0.68)\\
8   &6.405  &2.451  &5.222 &-1.990   &0.030    &580.876  &0.82  &126.039  &440 &2.26 (1.81)& 13.2 (19.5) & 1.54 (0.68)\\
{\bf 9}   & {\bf 6.474}  & {\bf 2.323}  & {\bf 5.242} &{\bf -2.533}   &{\bf 0.300}  &{\bf 562.500}  &{\bf 0.83}  &{\bf 129.465}  &{\bf 391} &{\bf 2.22 (1.77)} & {\bf 13.0 (19.3)} & {\bf 1.60 (0.71)}\\
10  &6.598  &2.159  &5.265 &-3.340   &0.445    &536.838  &0.84  &134.378  &344 &2.16 (1.72)& 12.8 (19.0) & 1.65 (0.71)\\
{\bf 11}  &{\bf 6.772}  &{\bf 1.995}  &{\bf 5.285} &{\bf -4.274}   &{\bf 0.292}  &{\bf 509.644}  &{\bf 0.85}  &{\bf 139.710}  &{\bf 303} &{\bf 2.10 (1.66)}& {\bf 12.6 (18.7)} & {\bf 1.76 (0.75)} \\
12  &7.022  &1.823  &5.305 &-5.414   &0.039    &478.498  &0.86  &146.131  &265 &2.04 (1.62)& 12.3 (18.3) & 1.82 (0.78)\\
{\bf 13}  &{\bf 7.325}  &{\bf 1.642}  &{\bf 5.324} &{\bf -6.586}   &{\bf 0.571} &{\bf 444.614}  &{\bf 0.87}  &{\bf 153.984}  &{\bf 231} &{\bf 1.97 (1.56)} & {\bf 12.0 (17.9)} & {\bf 1.94 (0.82)} \\
14  &7.865  &1.451  &5.343 &-8.315   &0.502  &403.303  &0.88  &163.824  &199 &1.89 (1.51)& 11.8 (17.5) & 2.08 (0.86) \\
15  &8.792  &1.249  &5.362 &-10.766  &0.354    &353.960  &0.89  &176.552  &168 &1.80 (1.44)& 11.4 (16.9) & 2.30 (0.94)\\
16  &7.942  &1.041  &5.388 &-6.908   &15.197   &339.910  &0.90  &193.437  &163 &1.71 (1.37)& 11.0 (16.3) & 2.54 (1.03)\\
\hline
\hline
\end{tabular}
\end{center}
\end{table*}

\begin{figure}[ht]
\vskip 0.3in
\begin{center}
\includegraphics[width=4cm,height=4cm,angle=0]{effm-k.eps}
\hskip 0.1in
\includegraphics[width=4cm,height=4cm,angle=0]{sigma-k.eps}
\end{center}
\caption{(A)- Nuclear matter incompressibility as a function of the nucleon effective mass for the parameters of the present model at $\rho_0 = 0.153 fm^{-3}$. Also plotted is the corresponding ratio of scalar to vector coupling on the opposite x-axis. (B)- Incompressibility as a function of obtained sigma meson mass for various parameter sets. The upper and the lower limit for incompressibility inferred from Heavy Ion Collision data \cite{data02} $K = (167-380)$ MeV is shown with horizontal lines. Recent experimental scalar meson mass limit ($m_{\sigma} = 513 \pm 32$) MeV \cite{mura02} is depicted with vertical lines}
\label{effm}
\end{figure}

Figure \ref{effm}(B) shows the variation of incompressibility as a function of scalar meson mass obtained for various parameter sets. Recent experimental estimate for scalar meson mass  $m_{\sigma} = 513 \pm 32$ MeV \cite{mura02} is compared with the present calculation. Here we find that the EOS with $Y = (0.84 - 0.86)$ seems to agree with the combined constraint from the HIC flow data and the experimental meson mass range.

\begin{figure}
\begin{center}
\includegraphics[width=4cm,height=4cm,angle=0]{sigma-fpi.eps}
\includegraphics[width=4cm,height=4cm,angle=0]{fpi-y.eps}
\end{center}
\caption{(A)- Sigma meson mass as a function of the pion decay constant at a fixed saturation density ($\rho_0 = 0.153 fm^{-3}$). (B)- Nucleon effective mass as a function of the pion decay constant. The experimental limit of the scalar meson mass and the pion decay constant are indicated with the horizontal and the vertical lines.}
\label{fpi}
\end{figure}

Figure \ref{fpi}(A) shows the obtained sigma mass as a function of the vacuum value of the pion decay constant. The experimental bound of the pion decay constant seems to agree with slightly higher value of $m_{\sigma}$, which agree with the upper bound of the experimental bound on $m_{\sigma}$ \cite{mura02}. Figure \ref{fpi}(B) shows the nucleon effective mass $Y = m^{\star}/m$ as a function of the pion decay constant. The constraint of $f_{\pi}$ agree with EOS with $Y = (0.82 -0.84)$, however the corresponding incompressibility lies on the higher side of presently acceptable bounds (\cite{k2} - \cite{k6}).

\begin{figure}
\vskip 0.3in
\begin{center}
\includegraphics[width=4cm,height=4cm,angle=0]{snm-p-expt.eps}
\includegraphics[width=4cm,height=4cm,angle=0]{expt-pnm.eps}
\end{center}
\caption{Comparison of the selected EOS with the heavy ion collision data at high densities \cite{data02}. (A)- Comparison with the Symmetric Nuclear Matter data (SNM). (B)- Comparison with the Pure Neutron Matter data (PNM).}
\label{snm}
\end{figure}

Fig. \ref{snm}(A) display the pressure as a function of baryon density up to nearly 6$\rho_0$ for the selected parameters of the model for symmetric nuclear matter. The shaded region corresponds to the experimental HIC data \cite{data02} for symmetric nuclear matter (SNM). Among the three theoretical calculations shown, the EOS with Y = 0.85 \& 0.87 agree very well with the collision data. Precisely, the third set (K = 231 MeV) completely agree with the flow data in the entire density span of $2 < \rho_B/\rho_0 < 4.6$. In Fig. \ref{snm}(B), the case of pure neutron matter (PNM) is compared with the experimental flow data. The experimental flow data is categorized in terms of stiff or soft based on whether the density dependence of the symmetry energy term is strong or week \cite{pra88}. The EOS predicted by the present model seems to rather lie on the softer regime. However, the EOS with $Y = 0.87, K = 231$ MeV though satisfy the combined constraint rather well, is not consistent with the vacuum value of the pion decay constant.

\begin{figure}
\begin{center}
\includegraphics[width=6cm,height=6cm,angle=0]{npe-hyp-stat.eps}
\end{center}
\caption{({\it Upper panel:-}) The mass of the star obtained as a function of the central density of the star. ({\it Lower panel:-}) The mass of the star as a function of the radius of the star. The two cases considered here are the star composed of nucleons only ($n, p, e, \mu$) shown with the filled circles and the one composed with hyperons, shown with the triangles.}
\label{star}
\end{figure}

In Fig \ref{star}, we plot the results for the global properties of the static star, obtained from the tabulated parameters of the model for two cases. In one, the star is composed purely of $\beta-$equilibrated nucleon only matter ($n,p,e,\mu$) and in the other case, the star is composed of matter which is Hyperon rich. To study hyperon rich matter, we fixed the hyperon couplings for all the parameter sets, so as to yield the binding energy of $\Lambda^0$ at $\rho_0 = -30 MeV$. This is obtained by fixing the scalar coupling $x_{\sigma}=x_{\sigma N}/ x_{\sigma H} = 0.70$ and varying the strength of the vector counterpart to yield the binding of $\Lambda$ as suggested by hypernuclei experiments. In the upper panel of the figure, the maximum mass of the star is plotted as a function of the central density of the star. Here it can be seen that the star with nucleon matter is much massive than those composed of hyperons. However in case of star composed of nucleons, the central density ranges from ($4.5 - 9.5$)$\times ~\rho_0$, but with hyperons, the central density of the star falls in the range $\approx$~($2 - 4$)$\times ~\rho_0$ and the resulting radius increases for the later case. For all the cases given in Table I, we find nearly ($46 -48$) \% increase in the star radius, when we move from nucleonic star to hyperon rich star matter. The corresponding decrease in maximum mass is $\approx (18 - 20)$\%. This feature is very much evident from the lower panel of Fig. \ref{star}, where we have plotted the maximum mass of the star obtained as a function of the star radius. In the transition of nucleon only matter to hyperon rich neutron star matter, the compactness ratio (M/R) of the star falls from ($0.15 - 0.17$) to ($0.09 - 0.10$).

\section{Summary}

The effective chiral model provides a natural framework to interlink the standard state properties of nuclear matter with the vacuum correlations, such as the pion decay constant. We find that the pion decay constant is experimentally well known quantity in comparison to other derived quantities such as the nuclear incompressibility and $\sigma-$meson mass, that can put stringent constraint on the model parameters. Experimentally determined effective mass from scattering of neutron over $Pb$ nuclei \cite{nuclei} seems to go well with the present model. Both of them favor higher value for nucleon effective mass. On a comparative analysis of the resulting EOS with that of the HIC data for symmetric nuclear matter as well as pure neutron matter, parameter set with $Y = 0.85; K \approx 300$ MeV seems to be the ideal parameterization of the present model. The resulting scalar meson mass $m_{\sigma} \approx 510 MeV$, is also consistent with the experimentally observed masses \cite{mura02,aitala01}. We then applied the model to study the global properties of neutron stars with and without hyperons. We find that the star composed of hyperons has a lower value of central density and results in larger radius than the star without hyperons. This in turn results in the decrease in the compactness parameter of the star. Although the resulting higher central density of the nucleonic star seems to be unphysical on account that at those higher densities, it is very unlikely that nucleons would retain its identity. Further, from the astrophysical point of view (observed neutron star masses), none of the parameters can be ruled out, however it shall be interesting to study the rotational attributes of the star with hyperons/ mixed/ quark phase in order to constrain the parameters of the model. Work is in progress in this direction \cite{tkj09}.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Automatic music generation using neural networks has attracted much attention.
There are two classes of music generation approaches, symbolic music generation~\cite{hadjeres2016deepbach}\cite{magenta2016}\cite{yang2017midinet} and audio music generation~\cite{oord2016wavenet}\cite{mehri2016samplernn}.
In this study, we focus on symbolic melody generation, which requires learning from sheet music.

Many music genres such as pop music consist of melody and harmony.
Since usually beautiful harmonies can be ensured by using legitimate chord progressions which have been summarized by musicians, we only focus on melody generation, similar to some recent studies~\cite{magenta2016}\cite{yang2017midinet}\cite{colombo2017deep}\cite{pmlr-v80-roberts18a}.
This greatly simplifies the melody generation problem.

Melody is a linear succession of musical notes along time.
It has both short time scale such as notes and long time scale such as phrases and movements,
which makes the melody generation a challenging task.
Existing methods generate pitches and rhythm simultaneously \cite{magenta2016} or sequentially \cite{chu2016song} using Recurrent Neural Networks~(RNNs), but they usually work on the note scale without explicitly modeling the larger time-scale components such as rhythmic patterns.
It is difficult for them to learn long-term dependency or structure in melody.

Theoretically, an RNN can learn the temporal structure of any length in the input sequence, but in reality, as the sequence gets longer it is very hard to learn long-term structure.
Different RNNs have different learning capability, e.g., LSTM~\cite{hochreiter1997long} performs much better than the simple Elman network.
But any model has a limit for the length of learnable structure, and this limit depends on the complexity of the sequence to be learned.
To enhance the learning capability of an RNN, one approach is to invent a new structure.
In this work we take another approach: increase the granularity of the input.
Since each symbol in the sequence corresponds to longer segment than the original representation, the same model would learn longer temporal structure.

To implement this idea, we propose a Hierarchical Recurrent Neural Network (HRNN) for learning melody. It consists of three LSTM-based sequence generators --- Bar Layer, Beat Layer and Note Layer.
The Bar Layer and Beat Layer are trained to generate bar profiles and beat profiles, which are designed to represent the high-level temporal features of melody.
The Note Layer is trained to generate melody conditioned on the bar profile sequence and beat profile sequence output by the Bar Layer and Beat Layer.
By learning on different time scales, the HRNN can grasp the general regular patterns of human composed melodies in different granularities, and generate melody with realistic long-term structures. This method follows the general idea of granular computing~\cite{bargiela2012granular}, in which different resolutions of knowledge or information is extracted and represented for problem solving.
With the shorter profile sequences to guide the generation of note sequence, the difficulty of generating note sequence with well-organized structure is alleviated.

\section{Related Work}

\subsection{Melody Generation with Neural Networks}
There is a long history of generating melody with RNNs.
A recurrent autopredictive connectionist network called CONCERT is used to compose music~\cite{mozer1994neural}.
With a set of composition rules as constraints to evaluate melodies, an evolving neural network is employed to create melodies~\cite{chen2001creating}.
As an important form of RNN, LSTM~\cite{hochreiter1997long} is used to capture the global music structure and improve the quality of the generated music~\cite{eck2002first}.
Boulanger-Lewandowski, Bengio, and Vincent explore complex polyphonic music generation with an RNN-RBM model~\cite{boulanger2012modeling}.
Lookback RNN and Attention RNN are proposed to tackle the problem of creating melody's long-term structure~\cite{magenta2016}.
The Lookback RNN introduces a handcrafted lookback feature that makes the model repeat sequences easier while the Attention RNN leverages an attention mechanism to learn longer-term structures.
Inspired by convolution, two variants of RNN are employed to attain transposition invariance~\cite{johnson2017generating}.
To model the relation between rhythm and melody flow, a melody is divided into pitch sequence and duration sequence and these two sequences are processed in parallel~\cite{colombo2016algorithmic}.
This approach is further extended in \cite{colombo2017deep}.
A hierarchical VAE is employed to learn the distribution of melody pieces in~\cite{pmlr-v80-roberts18a}, the decoder of which is similar to our model.
The major difference is that the higher layer of its decoder uses the automatically learned representation of bars, while our higher layers use predefined representation of bars and beats which makes the learning problem easier.
Generative Adversarial Networks~(GANs) have also been used to generate melodies. For example,
RNN-based GAN~\cite{mogren2016c} and CNN-based GAN~\cite{yang2017midinet} are employed to generate melodies, respectively.
However, the generated melodies also lack realistic long-term structures.

Some models are proposed to generate multi-track music.
A 4-layer LSTM is employed to produce the key, press, chord and drum of pop music seperately~\cite{chu2016song}.
With pseudo-Gibbs sampling, a model can generate highly convincing chorales in the style of Bach~\cite{colombo2017deep}.
Three GANs for symbolic-domain multi-track music generation were proposed~\cite{DBLP:conf/aaai/DongHYY18}.
An end-to-end melody and arrangement generation framework XiaoIce Band was proposed to generate a melody track with accompany tracks with RNN~\cite{Zhu:2018:XBM:3219819.3220105}.

\subsection{Hierarchical and Multiple Time Scales Networks}
The idea of hierarchical or multiple time scales has been used in neural network design, especially in the area of natural language processing.
The Multiple Timescale Recurrent Neural Network~(MTRNN) realizes the self-organization of a functional hierarchy with two types of neurons ``fast'' unit and ``slow'' unit~\cite{yamashita2008emergence}.
Then it is shown that the MTRNN can acquire the capabilities to recognize, generate, and correct sentences in a hierarchical way: characters grouped into words, and words into sentences~\cite{hinoshita2011emergence}.
An LSTM auto-encoder is trained to preserve and reconstruct paragraphs by hierarchically building embeddings of words, sentences and paragraphs~\cite{li2015hierarchical}.
To process inputs at multiple time scales, the Clockwork RNN is proposed, which partitions the hidden layers of RNN into separate modules~\cite{koutnik2014clockwork}.
Different from the Clockwork RNN, we integrate the prior knowledge of music in constructing the hierarchical model and feed multiple time scales of features to different layers.

\section{Music Concepts and Representation}
\label{Sec:BMT}
We first briefly introduce some basic music concepts and their properties to familiarize the readers who do not have a music background, then explain how the concepts are represented in the model.

\subsection{Basic Music Concepts}

\begin{figure}
\centering
\includegraphics[width=0.48\textwidth]{Figs/melody.pdf}
\caption{A typical form of melody. The time signature of this musical piece is 4/4. The numerator means a bar contains 4 beats, and the denominator means the time length of 1 beat is a quarter note.}
\label{Fig:Melody}
\end{figure}

As shown in Fig.~\ref{Fig:Melody}, melody, often known as tune, voice, or line, is a linear succession of musical notes, and each note represents the pitch and duration of a sound. Several combined notes form a beat that determines the rhythm based on which listeners would tap their fingers when listening to music.
A bar contains a certain number of beats in each musical piece.
Time signature~(e.g., 3/4) specifies which note value is to be given in each beat by the denominator and the number of beats in each bar by the numerator.
Each musical piece has a key chosen from 12 notes in an octave.
Key signature, such as C$\sharp$ or B$\flat$, designates which key the current musical piece is.
The musical piece can be transformed to different keys while maintaining the general tone structure.
Therefore we can transpose all of the musical pieces to key C, while maintaining the relative relationship between notes.
Shifting all musical pieces to the same key makes it easier for the model to learn the relative relationship between notes.
The generated pieces can be transposed to any key.

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{Figs/profiles_aug_v2.pdf}
\caption{
Samples of beat profiles and bar profiles.
Here we use notes with same pitch to illustrate rhythm in beat and bar.
The rhythm represented by beat profile 5 and 6 are related with the rhythm of the previous beat so they are shown with two beats where the first beats are all quarter notes.
}
\label{Fig:Profiles}
\end{figure*}

\begin{figure}
\centering
\includegraphics[width=0.8\linewidth]{Figs/Representation.pdf}
\caption{An example of melody representation. \textbf{Top}: A melody with the length of one bar. \textbf{Bottom}: Representation of the melody, in which the N means a no-event and the O means a note-off event. Since the fourth note is not immediately followed by any note, a note-off event is necessary here.}
\label{Fig:Representation}
\end{figure}

\begin{figure*}[htpb]
\centering
\includegraphics[width=\textwidth]{Figs/Arch_v3.pdf}
\caption{Architecture of HRNN. From top to bottom are Bar Layer, Beat Layer and Note Layer respectively. Inner layer connections along time are shown with black lines. Connections between layers are shown with green lines, blue lines and red lines. }
\label{Fig:Arch}
\end{figure*}

\subsection{Melody Representation}
\label{Sec:MelodyRepresentation}
To simplify the problem, we only chose musical pieces with the time signature 4/4. This is a widely-used time signature.
According to the statistics on the Wikifonia dataset described in Section~\ref{Sec:Exp}, about 99.83\% of notes have pitches between C2 and C5.
Thus, all notes are octave-shifted to this range. Then there are 36 options for a pitch of a note (3 octaves and each octave has 12 notes).
To represent duration, we use event messages in the Midi standard.
When a note is pressed, a note-on event with the corresponding pitch happens; and when the note is released, a note-off event happens.
For a monophonic melody, if two notes are adjacent, the note-on event of the latter indicates the note-off event of the former, and the note-off event of the former is therefore not needed.
In this study, every bar was discretized into 16 time steps.
At every time step, there are 38 kinds of events (36 note-on events, one note-off event and one no-event), which are exclusive.
One example is shown in Fig.~\ref{Fig:Representation}.
In this way, note-on events mainly determine the pitches in the melody and no-events mainly determine the rhythm as they determine the duration of the notes.
So a 38-dimensional one-hot vector is used to represent the melody at every time step.

\subsection{Rhythmic Patterns and Profiles}

Rhythmic patterns are successions of durations of notes which occur periodically in a musical piece.
It is a concept on a larger time scale than the note scale and is important for melodies' long-term structure.
Notice that in this model we do not encode the melody flow because it is hard to find an appropriate high-level representation of it.

Two features named \textit{beat profile} and \textit{bar profile} are designed, which are high-level representations of a whole bar and beat, respectively.
Compared with individual notes, the two profiles provide coarser representations of the melody.
To construct the beat profile set, all melodies are cut into melody clips with a width of one beat and binarized at each time step with $1$ for an event~(note-on events and note-off event) and $0$ for no-event at this step.
Then we cluster all these melody clips into several clusters via the K-Means algorithm and use the cluster centers as our beat profiles.
Given a one beat melody piece, we can binarize it in the same manner and choose the closest beat profile as its representation.
The computation of bar profile is similar, except that the width of melody clip is changed to one bar.
Based on the well-known elbow method, the numbers of clusters for beat profiles and bar profiles are set to be 8 and 16 respectively.
In Fig.~\ref{Fig:Profiles}, some frequently appeared profiles are shown with notes.

\section{Hierarchical RNN for Melody Generation}
\label{Sec:Methods}

\subsection{Model Architecture}\label{Method:HMG}
HRNN consists of three event sequence generators: Bar Layer, Beat Layer and Note Layer, as illustrated in Fig.~\ref{Fig:Arch}.
These layers are used to generate bar profile sequence, beat profile sequence and note sequence, respectively.

The lower-level generators generate sequences conditioned on the sequence output by the higher-level generators.
So to generate a melody, one needs to first generate a bar profile sequence and a beat profile sequence in turn.
Suppose that we want to generate a melody piece with the length of one bar, which is represented as $n_t, ..., n_{t+15}$(see Fig.~\ref{Fig:Arch}).
First, the Bar Layer generates a bar profile $B_{t}$ with the last bar profile $B_{t-16}$ as input.
Then the Beat Layer generates 4 beat profiles $b_t, ..., b_{t+12}$ with  $b_{t-4}$ as input conditioned on the bar profile $B_t$.
To generate the notes $n_t, n_{t+1},...n_{t+3}$, the Note Layer is conditioned on both $B_t$ and $b_t$;
to generate the notes $n_{t+4}, ..., n_{t+7}$, the Note Layer is conditioned on both $B_t$ and $b_{t+4}$; and so on.
In this way, each bar profile is a condition for the 16 generated notes and each beat profile is a condition for the 4 generated notes.

All of the three layers use LSTM but the time scales of the input are different.
Theoretically, the Beat Layer and Bar Layer can learn 4 and 16 times longer temporal structure than the Note Layer, respectively.
Note that it is difficult to quantify the length of temporal structure learned in a model, since ``temporal structure'' is an abstract concept and its characterization is still an open problem.
We could only probe the difference in length produced by different models indirectly by measuring the quality of the generated sequences using behavior experiments (see Section 5).

To explicitly help RNN memorize recent events and potentially repeat them, a Lookback feature was proposed for \textit{the Lookback RNN}~\cite{magenta2016}.
A user study suggested that the RNN with Lookback feature outperforms basic RNN~\cite{yang2017midinet} so we also use it in our model\footnote{For fair comparison in experiments, all models were equipped with this feature.}.
The lookback distance is 2 and 4 for the Bar Layer, 4 and 8 for the Beat Layer, 4 and 8 for the Note Layer.
Therefore, the Note Layer without the condition of the Beat layer and Bar layer is equivalent to \textit{the Lookback RNN}.

\subsection{LSTM-Based Event Sequence Generator}
Bar profiles, beat profiles and notes can be abstracted as events, which can be generated by RNN.
It might be better to use different models for generating different types of events, but for simplicity we use the same LSTM-based event sequence generator for the Bar Layer, Beat Layer and Note Layer.

The event sequence generator $G_\theta$ is trained by solving the following optimization problem:
\begin{equation}\label{equ:condition}
\max_\theta \sum_{y \in \mathcal{Y}} \sum_{t=1}^{\text{len}(y)} \log p(y_t|y_0,...,y_{t-1}, c_t)
\end{equation}
where $\theta$ are the parameters of the generator, $y$ is a sequence sampled in the event sequences dataset $\mathcal{Y}$.
And $y_t$ denotes the $t$-th event in $y$, $c_t$ denotes the condition for $y_t$.

LSTM is used to predict the conditional probability in Eq.~\eqref{equ:condition}, which is characterized by input gates $i_t$, output gates $o_t$ and forgetting gates $f_t$~\cite{hochreiter1997long}:
\begin{equation}
\begin{aligned}
i_t & = \sigma(W_{ix}x_t + W_{im}m_{t-1}) \\
f_t & = \sigma(W_{fx}x_t + W_{fm}m_{t-1}) \\
o_t & = \sigma(W_{ox}x_t + W_{om}m_{t-1}) \\
c_t & = f_t \odot c_{t-1} + i_t \odot \tanh(W_{cx} x_t + W_{cm} m_{t-1}) \\
m_t & = o_t \odot c_t \\
p_t & = \text{Softmax}(m_t)
\end{aligned}
\end{equation}
where $W_{ix}$, $W_{im}$, $W_{fm}$, $W_{fm}$, $W_{ox}$ and $W_{om}$ are trainable parameters, $\odot$ denotes the element-wise multiplication and $\sigma(\cdot)$ denotes the sigmoid function.
The $y_{t-1}$ is used as input $x_t$.

The lookback feature is added to the Bar Layer, Beat Layer and Note Layer, to help the model memorize recent events and potentially repeat them. The lookback distance is 2 and 4 for the Bar Layer, 4 and 8 for the Beat Layer, 4 and 8 for the Note Layer.

During generation, given a primer sequence as an initial input sequence, the LSTM network generates the distribution $p_0$ over all candidate events.
The next event was chosen by sampling over $p_0$.
The successive events are generated according to $p(y_t|y_0,...,y_{t-1}, c_t)$.

\section{Experiments}
\label{Sec:Exp}

Evaluating the performance of the models for melody generation is difficult.
The main reason is that measuring the quality of the generated melodies is subjective and it is hard to find an objective metric.

We evaluated three generative models, HRNN-1L, HRNN-2L and HRNN-3L mainly based on behavioral experiments.
HRNN-3L is the model we described in the previous section.
HRNN-2L is the HRNN-3L without the Bar Layer while HRNN-1L is the HRNN-3L without the Bar Layer and the Beat Layer.
Note that HRNN-1L is actually the \textit{Lookback RNN} developed by Google Magenta \cite{magenta2016}.
The music pieces generated by the models were not post-processed.

All melodies used in experiments were publicly available \footnote{\url{https://www.dropbox.com/s/vnd6hoq9olrpb5g/SM.zip?dl=0}}.

\subsection{Implementation Details}
All LSTM networks used in experiments had two hidden layers and each hidden layer had 256 hidden neurons.
They were trained with Adam algorithm~\cite{kingma2014adam} and the initial learning rate was 0.001.
The minibatch size was 64.
The $\beta_1$, $\beta_2$ and $\epsilon$ of Adam optimizer were set to 0.9, 0.999, 1e-8.
To avoid over-fitting, dropout with ratio 0.5 was adopted for every hidden layer of LSTM and validation-based early stopping~(see Fig.~\ref{Fig:Accuracy}) was employed so that the training was stopped as soon as the loss on the validation set increased for 5 times in a row~(the model is evaluated on validation set every 20 training iterations).

In each generation trial, primer sequences~(both profiles and melodies) were randomly picked from the validation dataset.
For the Bar Layer and Beat Layer, one profile is given as the primer.
For the Note Layer, the length of the primer sequence is 1 beat.
Beam search with a beam of size 3 was used in all experiments.

\subsection{Dataset}\label{SubSec:Dataset}
We collected 3,859 lead sheets with the time signature of 4/4 in MusicXML format from \url{http://www.wikifonia.org}. We have made these lead sheets publicly available\footnote{\url{https://www.dropbox.com/s/x5yc5cwjcx2zuvf/Wikifonia.zip?dl=0}}
90\% of the lead sheets were used as training set and the other 10\% were used as validation set.
The speed of most music pieces in the dataset is 120 beats per minute.
To guarantee the correct segmentation of melodies, all melodies started with weak beats were removed so that we can take bar as a basic unit.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{Figs/Figure_5.pdf}
\caption{The accuracy curves of Note Layer for training and validation dataset. \textbf{Left}: no-event accuracy curves. \textbf{Right}:  event~(both note-on event and note-off event) accuracy curves. Arrows indicate iterations at which training stopped to prevent over-fitting.}
\label{Fig:Accuracy}
\end{figure}

\subsection{Guiding Effect of Profiles}

To verify whether the beat and bar profiles can guide the generation of melody, we plotted the Note Layer's accuracy curves of in Fig.~\ref{Fig:Accuracy}.
Here both event accuracy~(accuracy of the note-on event and note-off event; chance level is $1/ 37$) and no-event accuracy~(accuracy of the no-event; chance level is $1/2$) are plotted.

With beat and bar profiles, the Note Layer learned the pattern of no-event quickly and easily.
For models with profiles, the accuracy of no-event increased to nearly 100\% at about 200 iterations while the model without profile converged slowly and over-fitting started after about 2000 iterations.
Since rhythm is encoded by no-event~(see Section~\ref{Sec:MelodyRepresentation}), this showed that the Note Layer successfully utilized the rhythm provided by the beat and bar profiles.
The accuracy of note-on and note-off events also improved,
which means models with profiles not only did a good job in predicting rhythm, but also in predicting pitch.

\begin{figure}
\centering
Beat profiles: "1, 1, 1, 2, ..., 1, 1, 1, 2"{
\includegraphics[width=0.48\textwidth]{Figs/44488.pdf}
\label{Fig:beat_profile_44488}
} \\
Beat profiles: "2, 1, 2, 1, ..., 2, 1, 2, 1"{
\includegraphics[width=0.48\textwidth]{Figs/884884.pdf}
\label{Fig:beat_profile_2121}
} \\
Beat profiles: "2, 5, 2, 5, ..., 2, 5, 2, 5"{
\includegraphics[width=0.48\textwidth]{Figs/2525.pdf}
\label{Fig:beat_profile_2525}
}
\caption{Melodies generated with given beat profile sequences.}
\label{Fig:beat_profile}
\end{figure}

\begin{figure}
\centering
Bar profiles: "1, 2, 1, 2, 1, 2, 1, 2"{
\includegraphics[width=0.48\textwidth]{Figs/bar-12.pdf}
} \\
Bar profiles: "2, 3, 2, 3, 2, 3, 2, 3"{
\includegraphics[width=0.48\textwidth]{Figs/bar-23.pdf}
} \\
Bar profiles: "3, 1, 3, 1, 3, 1, 3, 1"{
\includegraphics[width=0.48\textwidth]{Figs/bar-31.pdf}
}
\caption{Melodies generated with given bar profile sequences.}
\label{Fig:bar_profile}
\end{figure}

With given profile sequences, the Note Layer will generate melodies with rhythm represented by profile sequence.
To show this, we used handcrafted profile sequences to guide the generation of the Note Layer.
Fig.~\ref{Fig:beat_profile} and Fig.~\ref{Fig:bar_profile} show generated melodies given beat profile sequences to HRNN-2L and bar profile sequences to HRNN-3L~(profile index in Fig.~\ref{Fig:Profiles}).
The results verified that the generated melodies are strongly constrained by the given profile sequence patterns.
The same conclusion can be obtained using fixed beat profiles and bar profiles extracted from existing melodies.
We extracted beat profiles and bar profiles of children's rhymes \textit{Twinkle, Twinkle Little Star} and generated melody conditioned on these profiles.
The result is shown in Fig.~\ref{Fig:Twinkle}. The audio files can be found in \textbf{Supplementary Materials}. The rhythm of the generated melody is unison with the original melody, which suggests that the beat profiles and bar profiles effectively guided the generation of melody.

\begin{figure}
\centering
Original melody{
\includegraphics[width=0.48\textwidth]{Figs/twinkle_origin.pdf}
} \\
Generated melody{
\includegraphics[width=0.48\textwidth]{Figs/twinkle_hrnn.pdf}
}
\caption{Original melody of Twinkle, Twinkle Little Star and  generated melody by Note Layer of HRNN-3L, given profiles of the original melody.}
\label{Fig:Twinkle}
\end{figure}

\subsection{Qualitative  Comparison}

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{Figs/123l_v3.pdf}
\caption{Melodies generated by HRNN-1L, HRNN-2L and HRNN-3L.
}
\label{Fig:123l}
\end{figure*}

The strong guiding effect of profiles implies that the Note Layer could output good melodies if higher layers could generate good profile sequences.
Since note sequences are much longer than their profile sequences, learning the latter should be easier than learning the former using the same type of model.
Thus, compared to HRNN-1L, melodies generated by HRNN-2L and 3L model should be more well-organized and keep better long-term structures.
The qualitative comparison verified this point.
Three typical music pieces generated by HRNN-1L, HRNN-2L, HRNN-3L with the same primer note were shown in Fig.~\ref{Fig:123l}.
The melody generated by HRNN-1L has basic rhythm, but also irregular rhythmic patterns.
And the melodies generated by HRNN-2L and HRNN-3L contain less irregular rhythmic patterns.

\subsection{Comparison of Different Number of Layers}
Three human behavior experiments were conducted to evaluate melodies generated by models.
For this propose, we built an on-line website where people could listen to melodies and give their feedback.
To model real piano playing scenario, sustain pedal effect was added to all model generated and human composed musical pieces evaluated in these experiments. This was achieved by extending all notes' duration so that they ended at the end of the corresponding bars.

\subsubsection{Two-Alternative Forced Choice Experiment}
\label{subsec:AFC}
We randomly provided subjects pairs of melodies with the length of 16 bars (about 32 seconds) and asked them to vote  (press one of two buttons in the experiment interface) which melody sounded better in every pair.
This is the two-alternative forced choice~(2AFC) setting. Subjects had infinite time for pressing the buttons after they heard the melodies. Pressing the button started a new trial.

Three types of pairs were compared:  HRNN-3L versus HRNN-1L, HRNN-2L versus  HRNN-1L and HRNN-3L versus HRNN-2L.
Each model generated a set of 15 melodies and in every trial two melodies were randomly sampled from the two corresponding sets.
Different types of pairs were mixed and randomized in the experiment.

Call for participants advertisement was spread in a social media.
1637 trials were collected from 103 IP addresses~(Note that one IP address may not necessarily correspond to one subject).
The results are shown in Fig.~\ref{Fig:Exp}.
In nearly two-thirds of trials, melodies generated by hierarchical models were favored (Pearson's chi-squared test, $p=3.96\times10^{-10}$ for HRNN-3L versus HRNN-1L and $p=2.70\times10^{-8}$ for HRNN-2L versus HRNN-1L).
In addition, subjects voted more for melodies generated by HRNN-3L than by HRNN-2L~($p=3.38\times10^{-6}$)

\begin{figure}
\centering
\includegraphics[width=\linewidth]{Figs/SCORE_v7.pdf}
\caption{Results of the 2AFC experiment~(\textbf{left}) and the melody score experiment~(\textbf{right}).}
\label{Fig:Exp}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.7\linewidth]{Figs/Figure_11.pdf}
\caption{Results of the music turing test.}
\label{Fig:turing}
\end{figure}

\subsubsection{Melody Score Experiment}
To quantitatively measure the quality of melodies generated by different models and verify the conclusion obtained in the online experiment,  we invited 18 subjects between ages of 18 and 25 to score these melodies.

Every subject was asked to score every melody used in 2AFC experiment with 5 levels: 5 the best and 1 the worst.  It took each subject about 24 minutes to finish the experiment.

We calculated the average score of every melody in Fig.~\ref{Fig:Exp}.
The results verified that the two additional layers improved the quality of melodies generated by the single-layer model~(two-tailed test, $p = 0.00576$ for HRNN-3L versus HRNN-1L).

\subsubsection{Control the Number of Parameters}
In the above experiment, the number of parameters of HRNN-3L (2.84M) was three times that of HRNN-1L (0.94M).
That might be the reason why HRNN-3L performed better.
So we trained a HRNN-1L model with 450 hidden neurons~(2.79M) and conducted another 2AFC experiment~(as described in Section 5.5.1) to compare their performances.
A total of 203 trials were collected from 21 IP addresses. In 127 trials (62.6\%, $p = 3.44\times 10^{-4}$), melodies generated by HRNN-3L were favored, which is similar to the result (63.3\%) in comparison with HRNN-1L with fewer parameters~(Fig.~\ref{Fig:Exp} left).
The results indicate that the better performance of the hierarchical structure was not mainly due to the increased number of parameters.

\subsubsection{Music Turing Test}
To compare the quality of the melodies generated by the models and the quality of melodies composed by human composers, a music ``Turing test" was carried out.
Only two models, HRNN-3L and HRNN-1L, were tested.
We found that without chord, it was difficult for the models to generate melodies that could fool human.
So chords were added as a condition of the Beat Layer and Note Layer in training and generation.
Chords and primer sequences used in the generation of a melody were extracted from the same music piece in the validation set.

A total of 50 musical pieces containing 8 bars were randomly chosen from the validation set as human composed music.
Then the HRNN-1L and the HRNN-3L both generated 25 melodies with the length of 8 bars. Sample
We provided subjects music pieces from these 100 examples and asked them to judge if they were composed by human by pressing one of two buttons in the experiment interface. Subjects had infinite time for pressing the buttons after they heard the melodies. Pressing the button started a new trial.
Feedback about the correctness of the choice was provided immediately after the subjects made the choice in every trial.
Then the subjects had a chance to learn to distinguish human-composed and machine-composed melodies, which made it hard for the models to pass the Turing Test.

In this experiment, we collected 4185 trials from 659 IP addresses, among which 1018 music pieces were generated by HRNN-1L, 1003 by HRNN-3L and 2164 by human. As shown in Fig.~\ref{Fig:turing}, 33.69\% of music pieces generated by HRNN-3L were thought to be human composed (or real), which is higher than the result of HRNN-1L, 28.68\%.

It is seen that not all music pieces sampled from the original dataset were thought to be composed by humans (only 74.95\% were correctly classified).
This implies that some music pieces generated by the models sounded better than human composed pieces,
and that the quality of the dataset is not very high.

\subsection{Comparison with Other Models}
\label{subsec:com_midinet}

Though many models have been proposed for melody generation, to the best of our knowledge, only the lookback RNN [21], attention RNN [21], MidiNet [23] and MusicVAE [19] have public available source codes.
These models represent the state-of-the-art in the area of melody generation.
It was reported in~\cite{yang2017midinet} that the attention RNN had similar performance to the lookback RNN.
Our previous experiments have shown that HRNN-3L performed better than the lookback RNN, i.e. HRNN-1L.
We then compared MidiNet and MusicVAE with HRNN-3L based on human evaluations.

\subsubsection{MusicVAE}

MusicVAE is a variatiaonal autoencoder that can generate melodies with the length of 16 bars.
We compared the HRNN-3L model with a MusicVAE trained on our dataset~(with the same training setting in the original paper) and their pretrained MusicVAE using the 2AFC setting separately.
Each model generated 15 melodies with the length of 16 bars.
In each experiment, we randomly provided subjects 20 pairs of melodies and subjects were asked to vote for the better sounded melody.

In the comparison between HRNN-3L and MusicVAE that was trained on our dataset, 435 trials were collected from 17 IP addresses.
In  317 trials~(72.6\%,  $p =1.41 \times 10^{-21}$), melodies generated by HRNN-3L were favored.
We found the rhythm of melody generated by this MusicVAE is chaotic.
One reason might be that the size of our dataset is too small compared with the dataset used in~\cite{pmlr-v80-roberts18a}.

We then compared HRNN-3L and Pretrained-MusicVAE.
 461 trials were collected from  21 IP addresses. In 293 trials~(63.5\%, $p =5.82 \times 10^{-9}$), melodies generated by HRNN-3L were favored.
We generated 200 melodies with the Pretrained-MusicVAE and the statistics on 200 melodies show that about 40\% of notes generated by Pretrained-MusicVAE had pitches lower than C2, which made some melodies sound strange.

\subsubsection{MidiNet}
Another 2AFC experiment was used to compare HRNN-3L with MidiNet~\cite{yang2017midinet}.
The MidiNet was trained on our dataset with the same training setting in the original paper.
Since MidiNet required chords as an input, chords were used as a condition for MidiNet and HRNN-3L.
MidiNet generated 25 melodies with the length of 8 bars conditioned on chords.
The 25 melodies of HRNN-3L used in the Music Turing Test were used here for comparison..
In this 2AFC experiment, we randomly provided subjects pairs of melodies (HRNN-3L versus MidiNet) and asked them to vote for the better sounded melody in every pair.
290 trials were collected from 28 IP addresses.
In 226 trials (77.9\%, $p =1.85 \times 10^{-21}$), melodies generated by HRNN-3L were favored.

\section{Discussions}
\label{Sec:Con}
In this paper, we present a hierarchical RNN model to generate melodies.
Two high-level rhythmic features, beat profile and bar profile, are designed to represent rhythm at two different time scales respectively.
The human behavior experiment results show that the proposed HRNN can generate more well-organized melodies than the single-layer model.
In addition, the proposed HRNN, though very simple, can generate better melodies than the well-known models MusicVAE and MidiNet.

In the Music Turing Test, only 33.69\% pieces generated by the proposed model were thought to be composed by human.
This proportion is still far from our expectation and there is still a long way to go for developing a perfect automatic melody generator.
However, under current technology, HRNN has achieved good enough results.
On one hand, one should notice that automatic generation of other forms of data is in the similar stage. For example, many state-of-the-art machine learning models trained on natural images~\cite{zhu2017unpaired}\cite{isola2016image} generate no more than 30\% images that can fool human.
On the other hand, the dataset used in this study is not good enough (only 74.95\% human composed pieces were thought to be composed by human) which has hurt the performance of the model.
If the model is trained on a dataset in which nearly all human composed pieces can be correctly classified, one may expect that about  44.9\% (=33.69/74.95) pieces generated by the model would fool human subjects.

The proposed approach of course have many limitations which should be considered in future. First, since we quantized a bar into 16 time step, the encoding could not
represent triplet or other types of rhythm. Second, we only selected musical pieces with 4/4 time signature from the original dataset for training. More time signatures should be taken into consideration
to improve the capability of the model. Third, we only considered beats and bars as larger units than notes, and did not consider phrases which are often present in pop music, since they are not labeled in the dataset. With larger time scale units, the model may output pieces with longer temporal structure.

\section*{Acknowledgements}
This work was supported in part by the National Natural Science Foundation
of China under Grant Nos. 61332007, 61621136008 and 61620106010.


text:
\subsection*{Acknowledgments}}
The authors would like to thank Xiaoming Xie and Mianheng Jiang for enlightening discussions, especially to Chunli Luo for her experimental assist. This work has been supported by the National Fundamental Research Program (under Grant No. 2011CB921300, 2013CB336800 and 2011CBA00300), the National Natural Science Foundation of China, the Chinese Academy of Science, and the Shandong Institute of Quantum Science \& Technology Co., Ltd.


text:
\section{Introduction}

Topological superfluid  (or superconductor) (TS) has a full pairing
gap in the bulk and is labeled by a non-zero integer topological
invariant \cite{hasan, qi}. From the famous bulk-boundary
correspondence such a topological integer ensures the existence of
gapless excitations on the boundary of the system, in other words
Majorana fermions (MF) \cite{majorana} in vortex core of pairing
order parameter. Roughly, MFs are neither fermions nor bosons but
non-Abelian anyons \cite{ia} and play an important role for the
realization of fault-tolerant topological quantum computation (TQC)
\cite{sarma}. The application prospect of MFs makes TS become one of
the hottest frontiers.

In the condensed matter physics some practical two-dimensional
systems have been theoretically proposed to realize TS \cite{read,
tsut, sato1, lee, fu, linder, jau, anj1}. In terms of these systems
the entrance into TS requires subtle adjustment of Hamiltonian and
it is very difficult in condensed matter physics, although MFs have
been detected in InSb nanowires contacted with one normal (Au) and
one superconducting electrode (NbTiN) \cite{rea}. In the light of
the disadvantage for condensed matter, TS has been also suggested in
cold Fermi gases owing to their many controllable advantages and
operabilities. Following the successful observation p-wave Feshbach
resonance (FR), Gurarie et al. \cite{gurarie} show that degenerate
Fermi gases near a p-wave FR naturally give a concrete realization
of TS. Zhang et al. \cite{zhang} propose to create TS directly from
an s-wave interaction making use of an artificially generated
spin-orbit coupling (SOC). In fact, SOC have been realized in a
neutral atomic Bose-Einstein condensate (BEC) by dressing two atomic
spin states with a pair of lasers and the same technique is also
feasible for cold Fermi gases \cite{e1, e2}. Realizing that in a
dual transformation SOC is formally equivalent to a p-wave
superfluid gap, Sato et al. \cite{sato} suggest to artificially
generate the vortices of SOC by using lasers carrying orbital
angular momentum. In terms of the latter two ways, SOC and a large
magnetic field are crucial in order to enter into TS.

In this paper we suggest to create TS from spinless Fermi gases in
p-band honeycomb optical lattices with so-called on-site rotation
(OSR), that rotates every lattice site around its own center but
keeps the whole lattice intact and has been realized for triangular
optical lattices \cite{gem}. As a matter of fact p-band Fermi gases
in honeycomb optical lattices in absence of OSR have shown many
interesting characteristics, such as ferromagnetism \cite{ferr} and
Wigner crystallization \cite{wigner, wigner1} associating with flat
bands, f-wave superfluidity with conventional pairing interaction
\cite{fwave}. The motivation for this paper comes from the Wu's work
on quantum anomalous Hall effect in the same system \cite{wu}. Under
single particle picture, Wu found that an arbitrary non-zero OSR not
only breaks time-reversal symmetry and changes the topological
properties of the system, but also drives a topological phase
transition when OSR is beyond a critical value. Here we add on-site
attraction interaction between p-band Fermi atoms into Hamiltonian
and ask a question whether OSR can drive a phase transition into TS.
The results are positive and OSR brings phase transitions not only
from normal superfluid (NS) to TS, but also among different TS. From
another perspective our work also can be considered as an extension
to \cite{fwave}, where f-wave superfluidity without OSR is
discussed. Thus we also investigate the effects of OSR on f-wave
superfluidity.

Experimentally the route to realize TS suggested here is also
feasible. On the one hand by placing two electro-optic modulators at
two of three laser beams which coherently superpose to form a
honeycomb lattice, OSR is available as illustrated in \cite{machi}.
On the other hand due to Pauli exclusion principle the occupation of
p-band is very convenient as long as the lowest s-band is fulfilled.
In addition, on-site attraction interaction can be enhanced by using
atoms with large magnetic moments, such as ${}^{167}$Er with $m =
7\mu_B$ on which laser cooling has been performed \cite{35}. In
contrast to \cite{zhang, sato}, where a pair of extra lasers and a
large magnetic field are needed to produce an effective SOC and
split two SOC bands respectively, our system is much simpler.

The organization of this paper is as follows. In section 2, we give
the model and at the mean-field level investigate the ground state
of the system by numerically minimizing the thermodynamic potential.
In section 3 by calculating TKNN number $I_{TKNN}$ \cite{TKNN} of
occupied bands addressing the topological properties of the model,
the topological phase diagram is obtained. In addition we also
investigate the properties of edge states to prove our results. A
brief conclusion is given in section 4.

\section{Model and Mean-Field Ground State}

The honeycomb optical lattice was realized experimentally by using
three laser beams with co-planar propagating wavevectors quite some
time age \cite{salomon}. It is well known that a honeycomb lattice
is not a Bravais lattice and there are two inequivalent sites in a
unit cell, denoted by A and B respectively. Fulfilling the lowest
s-band and defining three unit vectors
$\vec{e}_1=\frac{\sqrt{3}}{2}\vec{e}_x+\frac{1}{2}\vec{e}_y$,
$\vec{e}_2=-\frac{\sqrt{3}}{2}\vec{e}_x+\frac{1}{2}\vec{e}_y$ and
$\vec{e}_3=-\vec{e}_y$, the Hamiltonian of p-band honeycomb optical
lattices with OSR is
\begin{eqnarray}
H=t_{\|}\sum_{\vec{r}\in A, i}\left[p_{\vec{r},
i}^{\dag}p_{\vec{r}+\vec{e}_i, i}+H.C.\right]- U\sum_{\vec{r}\in
A\oplus
B}p_{\vec{r}x}^{\dag}p_{\vec{r}y}^{\dag}p_{\vec{r}y}p_{\vec{r}x}-\Omega
\sum_{\vec{r}\in A\oplus B}\hat{l}_{\vec{r}, z}-\mu\sum_{\vec{r}\in
A\oplus B}\hat{n}_{\vec{r}}, \label{1}
\end{eqnarray}
where $p_{\vec{r}, i}=(p_{\vec{r}, x} \vec{e}_x+p_{\vec{r}, y}
\vec{e}_y)\cdot\vec{e}_i$ and $p_{\vec{r}, x}$ ($p_{\vec{r}, y}$) is
the annihilation operator for $p_x$ ($p_y$) band at the lattice site
$\vec{r}$. $\hat{n}_{\vec{r}}=p_{\vec{r}, x}^{\dag} p_{\vec{r}, x}+
p_{\vec{r}, y}^{\dag} p_{\vec{r}, y}$ and $\hat{l}_{\vec{r},
z}=-i(p_{\vec{r}, x}^{\dag} p_{\vec{r}, y}- p_{\vec{r}, y}^{\dag}
p_{\vec{r}, x})$ represent particle number and orbital angular
moment operators. $t_{\|}$ is the nearest-neighbor hopping matrix
element of atoms in $\sigma$ bonds and positive due to the odd
parity of the p-orbital. $U$ ($>0$), $\Omega$ ($>0$) and $\mu$ are
the on-site interaction strength, on-site rotation angular velocity
and chemical potential, respectively. Note that we have neglected
the nearest-neighbor atom hopping of $\pi$ bonds and supposed the
nearest neighbor distance in the lattice to be unit.

When $U=0$, introducing operator $\phi(k)=[p_{Ax}(k), p_{Ay}(k),
p_{Bx}(k), p_{By}(k)]^T$ and making a unitary transformation
$\phi_n(k)=U_{nm}(k)\Psi_m(k)$, Hamiltonian can be diagonalized
exactly. Meanwhile four energy bands can be obtained. Wu found two
of four bands always are topological for any nonzero OSR and the
others can be topological only if OSR is beyond a critical value
\cite{wu}. On the basis of this findings, Wu proposed an orbital
analogue of the quantum anomalous Hall effect, arising from orbital
angular momentum polarization due to OSR. With $\Omega=0$, Lee et
al. discussed f-wave superfluidity and charge density wave (CDW) in
this system at the mean-field level \cite{fwave}. Their results show
that away from the half filling the system is f-wave superfluidity,
while around the half filling superfluidity and CDW coexist and the
system is a supersolid. Although superfluidity exists all the time,
it is not topological as stated below.

Following the same spirit in \cite{fwave} we decouple interaction
term into CDW channel
\begin{eqnarray}
H_{int}^{CDW} =\sum_{\tau=x,y}\left[\sum_{\vec{r}\in
A}(-\frac{n}{2}U-\frac{\Delta_{CDW}}{2})p_{\vec{r}, \tau}^{\dag}
p_{\vec{r}, \tau} + \sum_{\vec{r}\in
B}(-\frac{n}{2}U+\frac{\Delta_{CDW}}{2})p_{\vec{r}, \tau}^{\dag}
p_{\vec{r}, \tau}\right] \label{2}
\end{eqnarray}
and pairing channel
\begin{eqnarray}
H_{int}^{pairing}&=&-\sum_k\left[\Delta_A p_{Ax}^{\dag}(k)
p_{Ay}^{\dag}(-k) + \Delta_B p_{Bx}^{\dag}(k) p_{By}^{\dag}(-k) +
H.C. \right] \nonumber\\
&=& -\sum_{k'}\left[\Delta_{nm}(k')\Psi_{n}^{\dag}(k')
\Psi_{m}^{\dag}(-k') + H.C.\right] \label{3}
\end{eqnarray}
where $n=<\hat{n}_{\vec{r}_A}+\hat{n}_{\vec{r}_B}>/2$ is filling
factor of every site,
$\Delta_{CDW}=U<\hat{n}_{\vec{r}_A}-\hat{n}_{\vec{r}_B}>/2$,
$\Delta_{A}=U\sum_k <p_{Ay}(-k) p_{Ax}(k)>$, $\Delta_{B}=U\sum_k
<p_{By}(-k) p_{Bx}(k)>$ are order parameters for CDW and
superfluidity. In (\ref{3}) we also express the pairing channel
using quasiparticle $\Psi(k)$. In this representation
\begin{eqnarray}
\Delta_{nm}(k)&=&\Delta_A\left[ U_{1n}^{\ast}(k)U_{2m}^{\ast}(-k) -
U_{2n}^{\ast}(k)U_{1m}^{\ast}(-k)\right]\nonumber\\ &+&
\Delta_B\left[ U_{3n}^{\ast}(k)U_{4m}^{\ast}(-k) -
U_{4n}^{\ast}(k)U_{3m}^{\ast}(-k)\right]. \label{4}
\end{eqnarray}

After the mean-field approximation, the Hamiltonian (\ref{1})
becomes a BdG Hamiltonian $H=[\phi^{\dag}(k), \phi(-k)] H_k
[\phi(k), \phi^{\dag}(-k)]^T$ and the properties of system are
completely decided by the $8\times 8$ matrix $H_k$. Diagonalizing
$H_k$, we attain spectrum $\epsilon_i(k)$ and correspondingly
eigenvectors $\varphi_i(k)$ ($i=1,2,\cdot\cdot\cdot, 8$). Due to
particle-hole symmetry inherent in this BdG Hamiltonian, the
spectrum are symmetric about zero energy and we assume
$\epsilon_1(k)=-\epsilon_8(k)>0$, $\epsilon_2(k)=-\epsilon_7(k)>0$,
$\epsilon_3(k)=-\epsilon_6(k)>0$, $\epsilon_4(k)=-\epsilon_5(k)>0$.
Then the thermodynamical potential at zero temperature is
\begin{eqnarray}
F=\frac{1}{2}\sum_k
\left[-4\mu-\epsilon_1(k)-\epsilon_2(k)-\epsilon_3(k)-\epsilon_4(k)\right]+\frac{N}{U}
|\Delta_A|^2+ \frac{N}{U} |\Delta_B|^2+\frac{N}{2U} \Delta_{CDW}^2,
\label{5}
\end{eqnarray}
where $N$ is the number of the unit cell. Below we numerically
minimize thermodynamic potential $F$ about $\Delta_A$, $\Delta_B$
and $\Delta_{CDW}$ for fixed interaction strength $U$. Without loss
of generality we choose $\Delta_A$ to be real,
$\Delta_B=|\Delta_B|e^{i\theta}$ and $U/t_{\|}=3.0$.

Fig.1 shows the solutions of the Hamiltonian (\ref{1}) at the
mean-field level for changing chemical potential $\mu$ and OSR
$\Omega$. Due to the particle-hole symmetry we only concentrate on
negative chemical potential. Fig.1(a) describes the variation of
$\Delta_{CDW}$. For $\Omega=0$ CDW is robust, but when $\Omega$ is
beyond a critical value $\Omega_c$ it vanishes suddenly. This is due
to the fact that the appearance of OSR changes the band structures
of single particle and breaks the nesting condition for CDW.
Numerically we find $\Omega_c/t_{\|}\approx 0.4\sim 0.6$ and is
monotonically increasing as the function of chemical potential.
Fig.1(b) shows the effect of OSR on particle density and further
exemplifies that the variations of band structures driven by OSR
cause nonmonotonic behavior of particle density. In contrast,
superfluid order parameters $\Delta_{A}$, $\Delta_{B}$ are more
interesting and shown in (c) and (d). On the one hand for
$\Omega>\Omega_c$, $\Delta_{A}=|\Delta_{B}|$ and with the increase
of OSR superfluid order smoothly decreases until disappearance. This
suppression mechanism of superfluidity consists in time-reversal
symmetry broken caused by OSR. While on the other hand for
$\Omega<\Omega_c$ $\Delta_{A}$ is still decreasing but $\Delta_{B}$
is increasing with $\Omega$. In fact the increase of $\Delta_{B}$
originates from the redistribution of particle density between sites
A and B, in other words the decrease of $\Delta_{CDW}$ as seen in
(a). Thus at the mean-field level our calculation suggests that (1)
OSR weakens stabilities of CDW and superfluidity and (2) for
$\Omega<\Omega_c$, superfluidity and CDW coexist and the system is a
supersolid.

The optimization of $\theta$ leads to $\theta=\pi$ for all
parameters we choose. Below we discuss the effects of OSR on pairing
symmetry for $\Omega>\Omega_c$. From \cite{fwave} without OSR and
away from the half filling the intraband pairings in (\ref{4}) have
f-wave symmetry with three nodal lines of $k_x = 0$, $k_y = \pm
k_x/\sqrt{3}$ and $\pi/3$ rotation symmetry [Fig.2(d)]. On
introducing OSR, in terms of the pairing magnitude, nodal lines
degenerate into some disconnected regions where intraband gap
disappears, and $\pi/3$ rotation symmetry retains. However real and
imaginary parts of pairing break $\pi/3$ into $\pi$ rotation
symmetry. In Fig.2(a) (b) and (c) as an example we show the
magnitude, real and imaginary parts of $\Delta_{11}$.

\section{Topological Phase Diagram and Majorana Fermion Modes}

In this section we discuss topological properties of Hamiltonian
(\ref{1}). In terms of our system, it explicitly breaks the
time-reversal symmetry due to OSR. Thus TKNN number $I_{TKNN}$ plays
a central role in deciding topological nature of the system
\cite{TKNN}. TKNN number is defined, by eigenvectors $\varphi_i(k)$
($i=5,6,7,8$) corresponding to negative energy spectrm of the matrix
$H_k$, into $ I_{TKNN}=\frac{1}{2\pi i}\int d^2k\, Tr\, dA$, where
$A$ is a matrix one-form $A_{ij}=A_{ij}^{\nu}(k)dk_{\nu}$ with
$A_{ij}^{\nu}(k)=\varphi_i^{\dag}(k) \nabla_{k_{\nu}}\varphi_j(k)$.
By numerically calculating TKNN number \cite{fukui}, we show the
topological phase diagram of the system in Fig.3. For parameter
region we choose, there are four different subregions labeled by
$I_{TKNN}=1,0,-1,2$ respectively. Moreover by comparison with
Fig.1(a) it is easily found that the boundary between $I_{TKNN}=0$
and other TKNN numbers in the direction of $\Omega$ coincides with
that of CDW disappearance. This finding is very important and
ensures that topological order of our system is not topological CDW
\cite{satom}. According to the criteria for TS \cite{satom}
$I_{TKNN}=2$ corresponds to Abelian TS while $I_{TKNN}=1,-1$ are
non-Abelian TS. Thus Fig.3 tells us that OSR drives topological
phase transition not only from NS to TS, but also between different
TS. Here we mention a fact the energy gap of the bulk spectrum
closes when topological phase transitions between topologically
distinct phases occur.

From the bulk-edge correspondence, a non-trivial bulk topological
number implies the existence of gapless edge states localized on
open edges of the system. Cold Fermi gases with sharp edges may be
realized along the lines proposed in \cite{goldman}. In order to
understand the relation between $I_{TKNN}$ and the number of edge
states, we study the Hamiltonian (\ref{1}) with the open boundary
condition along the zigzag edge of the honeycomb lattice. The
resulting excitation spectrum are depicted in Fig.4 for
representative parameter choices. Very explicitly the number of
gapless states for every edge is one-to-one correspondence with the
TKNN number. For $I_{TKNN}=\pm 1$ ($I_{TKNN}= 2$) there are one
(two) pair(s) of gapless states, while for $I_{TKNN}=0$, gapless
state does not exist. Due to particle-hole symmetry, in terms of
gapless states, they are Majorana fermion modes. It should also be
remembered that the core of a vortex is topologically equivalent to
an edge which has been closed on itself. The edge modes we describe
are therefore equivalent to the Majorana fermions known to exist in
the core of vortices of p-wave superfluids \cite{lewenstein}.

\section{Conclusions}

In conclusion at the mean-field level we have investigated the
effects of OSR on CDW and superfluidity for p-band spinless Fermi
gases in honeycomb optical lattices. We found that OSR weakens the
stabilities of CDW and superfluidity simultaneously, although
superfluidity can survives a larger OSR. This conclusion leads to
another important result that once CDW drops out the system enters
into topological superfluidity. By numerically calculating the TKNN
number we obtained topological phase diagram of the system. In
addition edge states, i.e. bulk-boundary correspondence are also
investigated.

\section*{Acknowledgement}

The work was supported by National Natural Science Foundation of
China under Grant No. 10675108 and Foundation of Yancheng Institute
of Technology under Grant No. XKR2010007.