\section{Introduction}

The key to new discoveries and precision measurements at the Large Hadron Collider (LHC) is the ability of excellent muon identification, triggering and precise momentum reconstruction with high efficiency.
The muon system of the Compact Muon Solenoid (CMS) experiment~\cite{CMS_experiment} and the muon reconstruction algorithms have been designed to achieve these goals.
The focus of this paper is the contribution of the RPC sub-system to the muon identification and reconstruction capabilities, alongside the Drift Tubes (DT) and Cathode Strip Chambers (CSC) independent sub-systems.

Muon reconstruction is performed using the central tracker and with all three muon detector systems.
Due to the redundancy provided by the combination of DT (CSC) and RPC in the barrel (endcap) region,
it is expected that the inclusion of the RPC measurements in the muon reconstruction
can improve the reconstruction of low momentum muons, especially in the case where the muons escape through the gaps between the wheels in which the chambers are divided along the axis parallel to the beam,
leaving hits in only one DT station.
Thus, the analysis presented here shows the impact of the RPC hits on the offline reconstruction
of the muons from the resonance decays of $J/\psi \rightarrow \mu^+\mu^-$ events in the region of $p_T$ below 20 GeV/\emph{c} and $Z^{0} \rightarrow \mu^+\mu^-$ events in the region of $p_T$ above.
Furthermore, the number of reconstructed hits in the RPC chambers which could be associated to a muon track has been compared to the Monte Carlo (MC) expectations.

\section{Detector Layout}

A schematic view of the CMS detector is shown in Figure~\ref{fig:MuonSys}.
The muon system covers the pseudo-rapidity region $|\eta|$ < 2.4~\cite{MUO-10-004}. DTs and CSCs detect muons in the regions of $|\eta|$ < 1.2 and 0.9 < $|\eta|$~<~2.4, respectively, and are complemented by the RPC system covering the range~of~$|\eta| $~<~1.6~\cite{RPC_chambers}.

For the study of the RPC contribution to muon reconstruction, we define three~regions in the detector, referred to as barrel ($|\eta|<0.8$), transition ($0.8<|\eta|<1.2$), and endcap ($1.2<|\eta|<1.6$).
The current design of the RPC endcap is not expected to be suitable for high particle rates at high~$\eta$ in the scenario of an LHC luminosity going up to $10^{34-35}~{\mathrm{cm}}^{-2} \cdot {\mathrm{s}}^{-1}$. There is ongoing work in identifying suitable technologies to instrument the high $\eta$ region
up to 2.4 matching the CSC system.
The RPCs are used at trigger level and also in the standard offline muon reconstruction.

In the barrel region the muon chambers are organised in four coaxial stations,
interleaved with iron return yokes.
A station is an assembly of chambers around a fixed value of $r$ in the barrel (or $z$ in the endcap).
Each of the two inner stations (MB1, MB2) contains layers of DTs sandwiched between two layers of RPCs, whereas each of the two outer stations (MB3, MB4) consists of one layer of RPCs and layers of DTs.
Along $z$, these barrel stations are grouped into five wheels, which are in turn divided into twelve $\phi$ sectors.
Similarly in the $r$ direction in the endcaps there are rings of CSCs, increasing with the radial distance from the beam line.
Each endcap ring is composed of 36 chambers covering the full azimuthal range.
The endcap region is composed of three iron disks holding a total of three RPC planes and four CSC planes on each side (plus/minus) along $z$.

\begin{figure}
  \begin{center}
    \texttt{\includegraphics[width=0.66\textwidth]{MuonSys-mod3}}
    \caption{
    Layout of one quadrant of CMS. The four DT stations in the barrel (MB1-MB4, green), the four CSC stations in the endcap (ME1-ME4, blue), and the RPC stations (red) are shown.
    }
    \label{fig:MuonSys}
    \vspace{-0.5cm}
  \end{center}
\end{figure}

\section{RPC Hits in Muon Reconstruction}
Muon tracking in CMS can be performed with the silicon tracker, and with either three or four stations of muon chambers installed outside the solenoid, sandwiched between steel layers serving both as hadron absorbers and as a return yoke for the magnetic field~\cite{PLB_chargeRatio}.

Three types of muon-track reconstruction were designed for muons originating from the LHC proton-proton collisions.
Muon tracks can be reconstructed by using hits in the muon detectors alone (standalone muons). Alternatively, the reconstruction can combine hits in the muons detectors with those in the central tracker (global muons). The muon system can also be used purely to tag extrapolated tracks from the central tracker; such tracks qualify as tracker muons if at least one muon segment (i.e. a short track stub made of DT or CSC hits) matches the extrapolated track.

For the reconstruction of the standalone- or global-muon track, at least two measurements, one of which must be of the DT or CSC segment, must be matched within the muon system.
The RPC hits are used by default in the standard muon reconstruction.
Figure~\ref{fig:hitDistr_STA} shows the number of RPC hits associated to global muons for the barrel, transition and endcap regions.
Figure~\ref{fig:fraction_STA} shows the average number of RPC hits of global muons as a function of $\eta$ and $\phi$ for given $\eta$ regions.
The dips in the barrel $\eta$ around $\pm$0.25 and $\pm$0.8 are due to the gaps between the 5 wheels of the yoke.
Similarly, the oscillations in the $\phi$ distribution for RPC hits with $|\eta|$ < 1.2 are due to the cracks between the adjacent sectors.
There is a small bias from the $Z^{0}$ event selection with at least one muon triggered.
The error bars indicate the statistical uncertainty and are smaller than the points.

\begin{figure}
  \begin{center}
    \texttt{\includegraphics[width=0.8\textwidth]{hitDistr_STA}}

    \caption{
    Number of RPC hits per global muon in $Z^{0}$ events for the barrel, endcap, transition minus and plus regions (clock-wise, starting from top-left). Results are shown for both data (dots) and MC (solid line). The error bars indicate the statistical uncertainty and are smaller than the points.
    }
    \label{fig:hitDistr_STA}
    \vspace{-0.4cm}
  \end{center}
\end{figure}

\section{Muon Efficiencies including or excluding RPC hits}
In this study we measure the efficiency of the selections
with given $p_T$ thresholds for low (3.5 < $p_T$ < 7 GeV/\emph{c} and 7 < $p_T$ < 20 GeV/\emph{c})
and high $p_T$ muons ($p_T$ > 20 GeV/\emph{c}) based on the following muon identification variables. $\chi^{2}$ divided by degrees of freedom in the global-muon track fit is less than 10, and its corresponding tracker track must use more than 10 silicon tracker hits
and have a transverse impact parameter $|d_{xy}| <$ 2 mm with respect to the primary vertex.
The minimal number of measurement points in the tracker is also required for the muons from the $J/\psi$ decays.

\begin{figure}
  \begin{center}
    \texttt{\includegraphics[width=0.44\textwidth]{AveHitVsEta_STA}}\hspace{0.15cm}\texttt{\includegraphics[width=0.55\textwidth]{AveHitVsPhi_STA}}
    \caption{
    Average number of RPC hits of global muons in $Z^{0}$ events vs. $\eta$ (left) and $\phi$ for the barrel, transition, and endcap minus and plus regions (right four plots, clock-wise, starting from top-left).  The error bars indicate the statistical uncertainty and are smaller than the points.
    }
    \label{fig:fraction_STA}

  \end{center}
\end{figure}

\begin{figure}
  \begin{center}
    \texttt{\includegraphics[width=0.45\textwidth]{AbsEff_pt_Data}}\hspace{0.5cm}\texttt{\includegraphics[width=0.45\textwidth]{AbsEff_pt_MC}}
    \texttt{\includegraphics[width=0.45\textwidth]{AbsEff_eta_Data}}\hspace{0.5cm}\texttt{\includegraphics[width=0.45\textwidth]{AbsEff_eta_MC}}
    \caption{

    For $Z^{0}$ events, efficiencies measured on data (left) and MC (right) as a function of $p_T$ (top) and $\eta$ (bottom).
    The black points are the global muon in the standard reconstruction ($\epsilon_{\mathrm{withRPC}}$) and the red open triangles are the global muon reconstruction by using hits in the tracker and DT/CSC only ($\epsilon_{\mathrm{withoutRPC}}$).
    }
    \label{fig:muEffAbs}
    \vspace{-0.4cm}
  \end{center}
\end{figure}

A global-muon track can be fitted combining hits from the tracker track and standalone-muon track, using the Kalman-filter technique~\cite{Kalman}.
The efficiency for global muons in the standard CMS reconstruction that uses RPC hits in the global fit ($\epsilon_{\mathrm{withRPC}}$) was measured and compared with the efficiency for muon reconstruction when RPC hits are excluded from the global fit ($\epsilon_{\mathrm{withoutRPC}}$).
In this way the RPC contribution to the global muon reconstruction efficiency is evidenced.

For these efficiency measurements, one can use the events with well known resonance decays into a pair of muons
and consider one of the muons as the tag to probe the efficiency using the other muon.
The tag is a well isolated muon, while the probe has less stringent selection requirements without isolation. We used the samples of $Z^{0}$ and $J/\psi$ events corresponding to an integrated luminosity of about $1~\rm{fb^{-1}}$ and $4~\rm{fb^{-1}}$, respectively: the simulated processes are $Z^{0} \rightarrow \mu^+\mu^-$ and $J/\psi \rightarrow \mu^+\mu^-$ only.
Events are selected from samples collected using the high-$p_T$ single-muon trigger
and the dedicated $J/\psi$ triggers to enrich the purity of the low-$p_T$ single-muon trigger by selecting events in which the muon track can be paired to a tracker track of opposite charge yielding an invariant mass close to that of the $J/\psi$ peak.
The tag is a global muon with the following additional requirements: associated with at least two track segments in different muon stations and at least one muon chamber hit included in the global-muon track fit. To further purify the sample the isolation cuts based on the tracker and calorimeters have been applied only for high-$p_T$ muons.

To measure efficiency for high-$p_T$ muons, the probe is a tracker track fulfilling the following selection criteria,
$p_{T}$ > 20 GeV/\emph{c}, $|\eta| <$ 1.8 and at least 10 hits in the tracker, so there is no bias from the muon sub-detectors.
The tag and probe tracks are also required to share the same primary vertex.
The invariant mass of the tag and probe tracks is fitted to extract the signal event yields for both the denominator and the numerator on the efficiency ratio, requiring all probes in the numerator match a global muon candidate reconstructed with and without
the RPC, respectively.

\begin{figure}
  \begin{center}
    \texttt{\includegraphics[width=0.45\textwidth]{effRatio_eta}}\hspace{0.5cm}\texttt{\includegraphics[width=0.45\textwidth]{effRatio_pt}}
    \caption{

    For $Z^{0}$ events, the ratio between the efficiencies with and without the inclusion of the RPC hits in the global muon reconstruction ($\epsilon_{\mathrm{withoutRPC}}/\epsilon_{\mathrm{withRPC}}$).
Results are shown for both data (dots) and MC (red open squares) as a function of $\eta$ (left) and $p_T$ (right), respectively.
    }
    \label{fig:muEffRatio}
    \vspace{-0.3cm}
  \end{center}
\end{figure}

Figure~\ref{fig:muEffAbs} shows
the efficiency of the muons coming from $Z^{0}$ decays, for both data and MC samples,
by applying a tag-and-probe technique to obtain almost unbiased estimates.
The standard global muon reconstruction with RPC is more efficient than the global muon reconstruction without RPC over the full $p_T$ and $\eta$ ranges, because it requires at least one muon segment if a segment matches at least one RPC hit within the muon system.
The efficiency drops with excluding the RPC in the global muon reconstruction,
which varies from a few percent to 3\% depending on $\eta$.

The ratio between the efficiencies for global muons reconstructed with and without the RPC ($\epsilon_{\mathrm{withoutRPC}}/\epsilon_{\mathrm{withRPC}}$) is shown in Figure~\ref{fig:muEffRatio}
as a function of $\eta$ and $p_T$, respectively.
The results on data are compared with the ones extracted applying the same procedure on MC samples.
Overall, there is good agreement between data and simulation in both the efficiencies and shapes of the distributions. Some discrepancies result from imperfect simulation of local detector conditions, affecting for example the muon identification efficiency~\cite{MUO-10-004}.

\begin{figure}
  \begin{center}

    \texttt{\includegraphics[width=0.45\textwidth]{RelEff_Medium_pt_PLOT}}\\
    \vspace{0.3cm}
    \texttt{\includegraphics[width=0.45\textwidth]{RelEff_Medium_eta_PLOT_pt_bin0}}\hspace{0.5cm}\texttt{\includegraphics[width=0.45\textwidth]{RelEff_Medium_eta_PLOT_pt_bin1}}
    \caption{

    For $J/\psi$ events, efficiency for muons that have been reconstructed with the standard algorithm (including RPC hits) to also be reconstructed once RPC hits are excluded.
Results are shown as a function of $p_T$ (top) and $\eta$ (bottom) for both data (dots) and MC (red open squares). The bottom plots show efficiencies for muons with 3.5 < $p_T$ < 7 GeV/\emph{c} (left) and 7 < $p_T$ < 20~GeV/\emph{c} (right), respectively, in the events selected by $J/\psi$ triggers with different thresholds on $p_T$.
    }
    \label{fig:muEffRel}
    \vspace{-0.3cm}
  \end{center}
\end{figure}

In the case of the $J/\psi$ events, combinatorial backgrounds from other tracks in the event are high.
Hence, for the muon efficiency at low $p_T$, the probe is a global muon with at least 10 hits in the tracker.
Figure~\ref{fig:muEffRel} shows the efficiency for muons that have been reconstructed with the standard algorithm (including RPC hits) to also be reconstructed once RPC hits are excluded.
The top plot shows the $p_T$ dependence of efficiency and the only significant discrepancy is around the turn-on of the efficiency curves, where the efficiency in data is systematically higher than that in the simulation (in the region $|\eta|$ < 1.2 as shown in the bottom-left plot). This discrepancy arises from a difference in the widths of the track-to-segment pulls between data and simulation.
The bottom plots show efficiencies for muons
in the events selected by $J/\psi$ triggers with different thresholds~on~$p_T$.

Some of tag-probe pairs may come from the resonance decays, while  other tag-probe pairs will be combinatorial background, where the probe is usually a charged hadron.
This combinatorial background can be removed in a simultaneous fit to the invariant mass spectra for passing and failing probes with identical signal shape and appropriate background shapes. The normalizations of the signal shapes in the two spectra are used to compute the efficiency.

The uncertainty on the fitted efficiency is determined from the likelihood function.
As normalizations of signal and background, efficiency of the background, and parameters controlling the shapes of the signal and background are all parameters of the fit, the uncertainty includes the contributions from the background subtraction procedure.
When the background is large, as in the lower mass resonances,
the uncertainty on the efficiency obtained by the fit is purely statistical~\cite{MUO-10-004}.
The typical statistical uncertainty on the $J/\psi$ results presented here
was estimated to be about $1\%$.

As a result, exclusion of the RPC hits in the global muon reconstruction degrades the muon efficiency by 1$\%$ on average, and to about 3$\%$ in the gaps between the wheels corresponding to certain eta region,
with tracker tracks from $Z^{0}$ decays for $p_T$ > 20 GeV/\emph{c} identified as possible muon candidates by applying the tag-and-probe method.
We see an overall good agreement between data and MC simulation within the
uncertainties.
Similar results are obtained with the muons coming from $J/\psi$ decays with 7 < $p_T$ < 20 GeV/\emph{c}.
Preliminary results on collision data at higher luminosities were checked and differences were found to be negligible.

\section{Conclusions}

The standard CMS muon reconstruction is performed using the central tracker and with all three muon detector systems: DT, CSC, and RPC.
The RPC contribution to muon reconstruction has been studied using
data collected in proton-proton collisions at $\sqrt s$ = 7 TeV during the 2011 LHC run. Tag-and-probe techniques were used to measure offline reconstruction efficiency for the muons in the region of $p_T$ above 3.5 GeV/\emph{c} using the well known resonance decays.
Using RPC hits in the muon track reconstruction improves efficiency of reconstructing $J/\psi$ and $Z^{0}$ events with muons in the gaps between the 5 wheels of the yoke (at $|\eta|$ = 0.25 and 0.8) and the transition region between the barrel outer wheels and the endcap disks ($0.8<|\eta|<1.2$).
The efficiency for reconstructing muons in the region of $p_T$ above 7 GeV/\emph{c} improves up to about 3$\%$ by the RPC chambers, in good agreement with simulation.

\acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC machine. We thank the technical and administrative staff at CERN and other CMS institutes, and acknowledge support from BMWF and FWF (Austria); FNRS and FWO (Belgium); CNPq, CAPES, FAPERJ, and FAPESP (Brazil); MEYS (Bulgaria); CERN; CAS, MoST, and NSFC (China); COLCIENCIAS (Colombia); MSES (Croatia); RPF (Cyprus); MoER, SF0690030s09 and ERDF (Estonia); Academy of Finland, MEC, and HIP (Finland); CEA and CNRS/IN2P3 (France); BMBF, DFG, and HGF (Germany); GSRT (Greece); OTKA and NKTH (Hungary); DAE and DST (India); IPM (Iran); SFI (Ireland); INFN (Italy); NRF and WCU (Korea); LAS (Lithuania); CINVESTAV, CONACYT, SEP, and UASLP-FAI (Mexico); MSI (New Zealand); PAEC (Pakistan); MSHE and NSC (Poland); FCT (Portugal); JINR (Armenia, Belarus, Georgia, Ukraine, Uzbekistan); MON, RosAtom, RAS and RFBR (Russia); MSTD (Serbia); SEIDI and CPAN (Spain); Swiss Funding Agencies (Switzerland); NSC (Taipei); ThEP, IPST and NECTEC (Thailand); TUBITAK and TAEK (Turkey); NASU (Ukraine); STFC (United Kingdom); DOE and NSF (USA).


\section{The CMS experiment at the LHC}
\label{section:intro}

The Large Hadron Collider (LHC) has become operational in 2009. High-energy
physics runs took place in 2010 and 2011
with proton-proton
collisions at a center-of-mass energy of 7 TeV,
and in 2012,
with proton-proton
collisions at a center-of-mass energy of 8 TeV.
The maximal instantaneous luminosity reached
$3.5 \cdot 10^{33}$ cm$^{-2}$ s$^{-1}$ in 2011 and almost twice
this value in 2012 at the time this article was written.

The Compact Muon Solenoid (CMS)~\cite{cms1, cms2} Collaboration,
one of the six experiments
currently operating at the LHC, consists of over 3000 scientists, engineers and
graduate students from 173 institutes in 40 Countries.

The central feature of the
CMS detector
is a superconducting solenoid,
of 6 m internal diameter, providing a field of 3.8 T.
Within the field volume are the silicon pixel and strip tracker,
the lead-tungstate crystal electromagnetic calorimeter, and the
brass-scintillator hadron calorimeter. Muons are measured in
gas-ionization detectors embedded in the steel return yoke.
In addition to the barrel and endcap detectors, CMS has
extensive forward calorimetry, assuring very good hermeticity
with pseudorapidity coverage up to high values ($|\eta| < 5 $).

Muons with pseudorapidity in the range $| \eta | < 2.4 $ are measured
with detection planes made of three technologies:
Drift Tube chambers (DT), Cathode Strip Chambers (CSC)
and Resistive Plate Chambers (RPC).
Matching the muons
to the tracks measured in the silicon tracker gives
a transverse momentum (\pt) resolution between 1\% and 5\%, for $\pt$ values up
to 1 TeV.

In 2011, CMS has recorded 5.20~\fbi of data out
of 5.72~\fbi delivered by the LHC, for an efficiency of 91\%. Roughly 93\%
of the recorded data has been certified as ``golden'' for all physics analyses.
An average of 98\% of the subdetector channels are operational and in the
readout.
In 2012, at the time this article was submitted, 9.98~\fbi of data out
of 10.66~\fbi delivered had already been recorded.

\section{The CMS RPCs}

\begin{figure}[hbp]
\includegraphics[width=.49\textwidth]{view}
\includegraphics[width=.49\textwidth]{rpc_gap}
\caption{Left: Schematic view, in the R-z plane, of one quadrant of the CMS detector.
The position of the RPC chambers is shown (not to scale).
Right: Schematic layout of a double-gap barrel chamber composed by two sub-units, called rolls.
The readout strip plane is also shown.
}
\label{fig:detector}
\end{figure}

The RPC~\cite{cms_rpc1, cms_rpc2} detectors are implemented in CMS
as a dedicated trigger system,
both in the barrel and in the endcap regions.
Figure~\ref{fig:detector} shows a schematic view of
one quarter of the CMS detector
in the R-z plane (Fig.~\ref{fig:detector}, left)
and the layout of a double-gap barrel chamber (Fig.~\ref{fig:detector}, right).
Two gas gaps, of 2 mm each, are formed by two parallel bakelite electrodes, with
one single plane of copper read-out strips  in-between. The two gaps feature
2 mm thickness and have a bulk resistivity of the order of $10^{10} \Omega$~cm.

The barrel RPC system consists of five wheels,
installed at $| \eta |< $ 0.8 and $| \mathrm{z} |< 7 $~m,
subdivided into 12 azimuthal sectors, each one equipped
with six radial layers of RPCs. Six endcap disks, three on the positive and
three on the negative endcap side, are divided into 36 azimuthal sectors,
with two radial rings in each one. They assure a full coverage
up to $| \eta |< $ 1.6.
In total, 480 barrel chambers and
432 endcap chambers are installed, adding up to 68136 barrel strips and
41472 endcap strips, respectively, covering a total surface of about 3000 m$^2$.
The CMS RPCs work in saturated avalanche mode and use a three-component,
non-flammable gas mixture composed of 95.2\% C$_2$H$_2$F$_4$ (R134a),
4.5\% iC$_4$H$_{10}$ (isobutane) and 0.3\% SF$_6$.
Water vapor is added in order to maintain
the relative humidity at constant values and allow for
constant bakelite resistivity.

\section{Resolution and efficiency studies: the method}
\label{section:method}

The presence, next to each RPC chamber,
of either a DT (in the barrel) or a CSC chamber (in the endcap), allows
to profit by the redundancy of the CMS muon system in order
to define the RPC hit efficiency independently of final physics ``objects''.
Muon track segments, reconstructed in the multi-layer detectors DT and CSC,
are linearly extrapolated to the RPC strip plane
and used to predict RPC hits in a fiducial region around the extrapolated
impact points, as illustrated in Figure~\ref{fig:extrapolation}.
Each extrapolated hit is matched to the closest RPC cluster, which
is formed by contiguous fired strips.
The typical cluster size corresponds to
about two strips, as shown in Section~\ref{section:stability}.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.45\textwidth]{extrap1.eps}
\caption{Schematic drawing of the extrapolation method. A DT or CSC
track segment is linearly extrapolated to the RPC strip plane. RPC hits
are sought in the acceptance region around the extrapolated impact point.
}
\label{fig:extrapolation}
\end{center}
\end{figure}

An additional requirement is applied to ensure that only extrapolated segments
associated to muon tracks with hits in the central tracker
are taken into account, while DT/CSC segments
probably originated by neutral background particles are discarded.
The three muon detectors employ different technologies and materials, thus
they have different sensitivities to the various backgrounds.

The extrapolation method allows to measure both the RPC spatial resolution, through
the residuals, and the hit efficiency.
The measured spatial resolution increases with increasing strip widths
and it ranges from 0.8 cm (for inner detector layers, characterized by smaller
strip pitches) to 1.3 cm (for outer
layers), both in the Barrel and in the Endcap.
The strip pitch dimensions are between 2.28 and 4.10 cm in the barrel and between
1.95 and 3.63 cm in the endcap.

\section{Detector uniformity. The High Voltage Scan}
\label{section:scan}

High Voltage (HV) scans were performed at the beginning of the 2011 and of the 2012
proton-proton LHC running,
aiming at determining the optimal operating HV for each individual chamber.
 The variation of
the environmental pressure $\mathrm P$ and the temperature $\mathrm T$ inside the CMS cavern
was taken into account as described in Section~\ref{section:stability}.

Collision data were taken at eleven different values of the effective HV ($\mathrm{HV_{eff}}$,
defined in Eq.~(\ref{eqn:corr}) below),
from 8.5~kV to 9.7~kV. It is worth mentioning that  only a
negligible amount of data (about 3~\pbi \,  out of 5.72 \fbi in 2011,
and about 6~\pbi in 2012)
was discarded, {\it i.e.} not included in the ``golden'' sample
for physics analyses, because of the RPC HV calibration.
A dedicated data stream was used, containing
information from the muon detectors and the first level trigger.

The efficiency ($\epsilon$) dependence on $\mathrm{HV_{eff}}$ is parametrized
by a sigmoidal
response function that can be written as:

\begin{equation}
\nonumber
\mathrm{
\epsilon = \frac{\epsilon_{max}} {1 + e^{-\lambda (HV_{eff} - HV_{50\%} )}}
}
\end{equation}
\noindent
where $\epsilon_{max}$ is the asymptotic efficiency for HV $\rightarrow
\infty$, the
$\lambda$ coefficient is proportional to the sigmoid slope at the inflection point, and
the High Voltage value HV$_{50\%}$ is the inflection point of the
sigmoid, for which 50\% of  $\epsilon_{max}$ is reached.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.75\textwidth]{P4}
\caption{
Efficiency measurement (in \%), as a function of the effective
High Voltage
$\mathrm{HV_{eff}}$ (in kV) for a typical chamber. The ``knee'' and the
optimal working point, shown in the picture, are defined in the text.
}
\label{fig:knee}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.65\textwidth]{HV50_new_nolegend}
\caption{
HV$_{50\%}$ (in kV) distribution for the barrel (left curves, red and blue)
and the endcap (right curves, green and black) rolls. The red and green curves
show 2012 data.
HV$_{50\%}$ is the inflection point
of the sigmoid, for which 50\% of the asymptotic efficiency is obtained.
}
\label{fig:hv50}
\end{center}
\end{figure}

Figure~\ref{fig:knee}
shows the efficiency as a function of the
effective High Voltage for a typical endcap chamber. The ``knee'' is the
$\mathrm{HV_{eff}}$ value for which 95\% of the asymptotic efficiency is reached.
The optimal working point (WP) is chosen beyond the knee, to ensure high efficiency,
and in the plateau region, to minimize the dependence on the environmental
parameters.
The WP is then defined for each individual chamber as the
knee value plus 100 V (for barrel chambers) or 150 V (for endcap chambers).
The difference between barrel and endcap reflects the different trigger algorithms.
The endcap trigger algorithm, requiring three coincidences out of three planes,
is more sensitive to efficiency variations.
The resulting WP values are averaged for chambers fed by the same HV supplier.

The results of the sigmoidal fit are highlighted in
Figure~\ref{fig:hv50},
presenting the HV$_{50\%}$ distribution
for barrel and endcap. The different HV$_{50\%}$ average values (around
8.9~kV for the barrel and 9.2~kV for the endcap) might be due to different
construction techniques. As shown in Figure~\ref{fig:hv50},
a high level of uniformity is obtained both for the barrel and for the
endcap chambers, with RMS values of the
HV$_{50\%}$ distributions of the order of 60~V and 80~V, respectively.

\section{Detector stability}
\label{section:stability}

 The variation of
the environmental pressure $\mathrm P$ and the temperature $\mathrm T$ inside the CMS cavern
was taken into account
by rescaling~\cite{scaling} to the
chosen reference values ($\mathrm P_0 = 965 $~mbar and
$\mathrm T_0 = 293$~K):

\begin{equation}
\nonumber
\mathrm{
HV_{eff}(P,T) = HV \cdot  \frac{P_0}{P} \cdot \frac{T}{T_0},
}
\label{eqn:corr}
\end{equation}
\noindent
where $\mathrm{HV_{eff}}$ is the resulting effective HV.
The dominant effect is due to the pressure variation: a 1\% variation, of the order of 10~mbar,
produces a sizeable $\mathrm{HV_{eff}(P,T)}$  difference of about 100 V.
Starting from July 2011, the $\mathrm{HV_{eff}(P,T)}$  correction of
Eq.~(\ref{eqn:corr}) was automatically implemented.

Figure~\ref{fig:eff} and Figure~\ref{fig:cls} show preliminary results
obtained with 2011 data. The efficiency and the cluster size stability
as a function of time have improved after the introduction of the
automatic $\mathrm{HV_{eff}(P,T)}$  correction of Eq.~(\ref{eqn:corr}).
The barrel
efficiency fluctuations (Fig~\ref{fig:eff}, left), mainly due to pressure
variations in the CMS cavern,
are reduced from about $ \pm 1\% $ to about $ \pm 0.5\% $.
The higher average efficiency ($\sim 97\%$ compared to $\sim 96\%$) in the first
part of 2011 is due to the choice of 965 mbar as a reference value for
rescaling. In the first part of the year, (Fig~\ref{fig:eff}) the
atmospheric pressure was on average lower than 965 mbar, giving rise
to higher $\mathrm{HV_{eff}}$ values and therefore to higher efficiency
values with respect to the second part of the year.

An increased stability, with reduced fluctuations, is also observed in
the endcap cluster size as a function of time (Fig~\ref{fig:eff}, right).
As mentioned above for the efficiency, the same considerations apply to the average values
before and after the $\mathrm{HV_{eff}(P,T)}$  correction.

Fig~\ref{fig:cls} presents the average cluster size in the barrel (left) and in
the endcap (right) as a function of the atmospheric pressure, before and
after the automatic correction. Both in the barrel and in the endcap, a
clear anti-correlation is shown in the first part of 2011, when the correction
was not applied.
The slight positive correlation in the second part of the year might hint
at an over-correction in Eq.~(\ref{eqn:corr}), currently under study.

\begin{figure}[hbtp]
\includegraphics[width=.5\textwidth]{eff_barrel}
\includegraphics[width=.5\textwidth]{CLSE}
\caption{Left: Barrel efficiency (in \%) as a function of the run number,
{\it i.e.} as a function of time, for runs taken between April 2011 and October
2011. The two regions, before and after the automatic $\mathrm{HV_{eff}(P,T)}$  correction,
are shown in the plot.
Right: Endcap cluster size (in number of strips), as a function of the run number,
for the same run range between April and October 2011.
}
\label{fig:eff}
\end{figure}

\begin{figure}[hbtp]
\includegraphics[width=.5\textwidth]{corr_barrel}
\includegraphics[width=.5\textwidth]{corr_endcap}
\caption{Average cluster size for the barrel (left) and for the endcap (right),
in number of strips, as a function of the atmospheric pressure P (in mbar),
for runs before (left, blue curve) or after (right, red curve) the
automatic $\mathrm{HV_{eff}(P,T)}$  correction. Preliminary results with
2011 data are shown.
}
\label{fig:cls}
\end{figure}

\section{Conclusions}

This note summarize new results highlighting the CMS RPC stability
and uniformity over the 2011 and 2012 data taking periods.
HV scans have been performed at the beginning of 2011 and again at
the beginning of 2012. They have been extremely effective,
allowing to select the optimal
operating HV values for each individual RPC chamber and to obtain a high level
of uniformity both in the barrel and in the endcap.
A new method for determining the
RPC hit efficiency is used by the CMS RPC Collaboration and automatic
$\mathrm{HV_{eff}(P,T)}$  corrections
are in place since July 2011.
Those efforts result in increased efficiency stability and increased
cluster size stability as a function of pressure and time.

\acknowledgments

We congratulate our colleagues in the CERN accelerator departments for
the excellent performance of the LHC machine. We thank the technical and
administrative staff at CERN and other CMS institutes, and acknowledge
support from BMWF and FWF (Austria); FNRS and FWO (Belgium);
CNPq, CAPES, FAPERJ, and FAPESP (Brazil); MEYS (Bulgaria);
CERN; CAS, MoST, and NSFC (China); COLCIENCIAS (Colombia); MSES (Croatia);
RPF (Cyprus); MoER, SF0690030s09 and ERDF (Estonia); Academy of Finland, MEC,
and HIP (Finland); CEA and CNRS/IN2P3 (France); BMBF, DFG, and HGF (Germany);
GSRT (Greece); OTKA and NKTH (Hungary); DAE and DST (India); IPM (Iran);
SFI (Ireland); INFN (Italy); NRF and WCU (Korea); LAS (Lithuania); CINVESTAV,
CONACYT, SEP, and UASLP-FAI (Mexico); MSI (New Zealand); PAEC (Pakistan);
MSHE and NSC (Poland); FCT (Portugal); JINR (Armenia, Belarus, Georgia,
Ukraine, Uzbekistan); MON, RosAtom, RAS and RFBR (Russia); MSTD (Serbia);
SEIDI and CPAN (Spain); Swiss Funding Agencies (Switzerland); NSC (Taipei);
ThEP, IPST and NECTEC (Thailand); TUBITAK and TAEK (Turkey); NASU (Ukraine);
STFC (United Kingdom); DOE and NSF (USA).


\section{Background} \label{sec:introduction}

Cosmological inflation \cite{Starobinsky:1980te,Guth:1980zm,Albrecht:1982wi,Hawking:1981fz,Linde:1981mu,Linde:1983gd}, typically assumed to be driven by one or more scalar fields, offers a compelling explanation for the origin of structure in our universe (see \cite{Lidsey:1995np,Lyth:1998xn, Lyth:2009zz,Liddle:2000cg} for reviews). But if inflation is indeed the source of cosmic structure, the field or fields driving this inflationary period have so far not been identified, and a great many candidates have been suggested, including scalar fields with non-canonical kinetic terms \cite{ArmendarizPicon:1999rj, Garriga:1999vw, Silverstein:2003hf}.

During inflation, quantum fluctuations are promoted to classical perturbations, and the statistics of these perturbations are probed by observations (for example \cite{Peiris:2003ff,Komatsu:2010fb}). In recent years, interest in the statistics beyond the two point function, which is parametrised by the power spectrum, has intensified. In particular, numerous investigations of the three-point function, parametrised by the bispectrum, have been undertaken \cite{Maldacena:2002vr, Seery:2005wm, Chen:2006nt}. The bispectrum is potentially a powerful diagnostic, for example it may allow non-canonical single field inflation, which is capable of producing an observable bispectrum \cite{Chen:2006nt}, to be distinguished from canonical inflation, for which the bispectrum is typically too small ever to be observed \cite{Maldacena:2002vr}.

\subsection{Three-form inflation}
In addition to scalar fields, other ways of driving inflation have been proposed. These include inflation driven by vector fields (for example \cite{Ford:1989me,Koivisto:2008xf,Golovnev:2008cf}) and higher form fields \cite{Germani:2009iq,KN1,KN2,Germani:2009gg, Koivisto:2009sd,DeFelice:2012jt}. For such models, the techniques needed to make observational predictions are less well explored than for scalar field models. In this paper, our interest lies in the recently proposed inflationary scenarios involving a massive three-form field. Induced potentials for three-forms were studied in the context of string theory in Refs.\cite{Gubser:2000vg,Frey:2002qc}. Inflationary scenarios with three-forms were first investigated by Germani \& Kehagias \cite{Germani:2009iq,Germani:2009gg}, who focused on a non-minimally coupled three-form, with the couplings chosen to make the theory exactly equivalent to a canonical and minimally coupled scalar field theory.
The behaviour of gravitational waves in this models was studied in Ref.\cite{Kobayashi:2009hj}.
Later Koivisto \& Nunes \cite{KN1, KN2} considered a minimally coupled massive three-form, which has a novel and much richer associated phenomenology. They found that if
particular conditions on the potential were satisfied, the model was ghost-free, and moreover could support inflation. More recently, Koivisto \& Urban \cite{Urban:2012ib} studied the non-Gaussian signatures resulting from the magnetic fields generated from the coupling of a three-form to electromagnetism.

Here we follow Koivisto \& Nunes and consider  the action for a minimally coupled, canonical three-form
\beta \label{3formS}
S = -\int d^4 x \sqrt{-g} \left(\frac{1}{2} R - \frac{1}{48}F^2 - V(A^2) \right)\,,
\end{equation}
where $8\pi G = 1$, $F= \nabla  A$ is a four-form, $A$ a three-form, and we have adopted  a compact notation in which indices are suppressed once the valance of an object is specified,  squaring a quantity denotes contraction of all indices (for example $F^2 \equiv F_{\mu\nu\rho\sigma} F^{\mu\nu\rho\sigma} \equiv \nabla_\mu A_{\nu\rho\sigma} \nabla^\mu A^{\nu\rho\sigma}$), and a dot ($\cdot$) denotes contraction on the first index (for example $\nabla \cdot A=\nabla_\mu A^\mu_{ \rho \sigma}$). Here and throughout Greek indices label space-time dimensions and Roman indices spatial dimensions.

In the context of inflationary cosmology,
the unperturbed three-form must be consistent with the universal symmetries of homogeneity and isotropy. The background three-form can then be written in the form
\beta
A_{ijk}=a(t)^3\epsilon_{ijk}\chi(t)\,
\end{equation}
where $a(t)$ is the scale factor of the universe as a function of cosmic time, and all other elements of $A$ are zero. Hence the unperturbed three-form satisfies $A^2 =6 \chi^2$ \footnote{Note that in Ref. \cite{KN2}, the three form quantity was denoted by $X$ not $\chi$, but here we reserve $X$ for use when we discuss non-canonical scalar field models to match notation commonly used in the literature.}. The dynamics of the universe is then governed by the behaviour of the scalar quantity $\chi(t)$ which is directly related to the three-form. As emphasised elsewhere \cite{KN1,KN2}, when written in terms of this scalar quantity the equations of motion which govern the behaviour of the universe, and the role of the three-form potential, are straightforward to interpret. In particular, adopting units where $8\pi G = 1$, the Friedmann and Raychaudhuri equations are given by
\begin{eqnarray}
\label{background}
3 H^2 &=&  \frac{1}{2}\left(\dot{\chi} + 3H\chi  \right)^2 + V(\chi^2) \nonumber \\
\dot{H} &=& -\frac{1}{2}V_{,\chi}\chi\,,
\end{eqnarray}
respectively.

On the other hand, as we will review in \S \ref{sec:dualities}, a series of dualities exist which relate the action \eqref{3formS} to the action for a massive vector field, and in turn to the action for a massive, and in general non-canonical, scalar field of the form
\beta \label{eq:Paction}
S = -\frac{1}{2}\int d^4 x \sqrt{-g} \left ( R - 2 P(X, \phi) \right)\,,
\end{equation}
where $\phi$ is a scalar field and $X=-\nabla^\mu \phi \nabla_\mu \phi$. In principle, therefore, if one wished to understand the dynamics and observational predictions of a particular three-form inflationary model, one could move to working with the dual non-canonical scalar field model, for which suitable techniques are already readily available \cite{Garriga:1999vw,Maldacena:2002vr,Seery:2005wm,Chen:2006nt}. Unfortunately the mapping between theories is rather complicated, and the analytical form of the dual non-canonical scalar field theory is often rather unwieldy or may not even exist, even for fairly simple three-form potentials.

With this in mind, Koivisto \& Nunes \cite{KN1,KN2} advocated working directly with the three-form. Utilising the Einstein field equations, they derived the background equations of motion, and perturbing both the three-form and the metric, the first order perturbation equations. Finally, they formed the Sasaki-Mukhanov equation for the curvature perturbation on comoving hypersurfaces, $\zeta$. This equation follows from the second order perturbed action written purely in terms of $\zeta$, which they thus implied. In addition they were able to address questions such as the stability of the theory, and the  power spectrum produced by three-form inflationary models. Later De Felice, Karwan \& Wongjun \cite{DeFelice:2012jt}  directly perturbed the action, and, following the standard approach, rewrote the result purely in terms of the curvature perturbation. In this way they obtained similar results,  though curiously their perturbed action did not match exactly that of Koivisto \& Nunes.

\subsection{Overview}

In this paper,  we initially follow Koivisto \& Nunes by working with the three-form theory itself rather than any dual. In contrast to their study, however, we manipulate the action directly, perturbing it to second order and writing the result purely in terms of $\zeta$. Doing so we find agreement with their result. Then we extend existing work by deriving the third order action, which is required to calculate the bispectrum of the curvature perturbation produced during inflation. We also explain why the two previously derived second order actions, by Koivisto \& Nunes and De Felice, Karwan \& Wongjun \cite{DeFelice:2012jt}  respectively, differed slightly. These calculations are covered in \S \ref{sec:actions}. In \S \ref{sec:dualities}, we review the dualities which relate the three-form theory to a non-canonical scalar-field theory. We then show that despite the complexity of the mapping between the theories, we may use this duality to provide an alternative derivation of our second and third order action, providing confirmation of our results. Finally, we employ our results to calculate the bispectrum in two explicit examples of three-form inflation, one of which produces a large non-Gaussianity, in \S \ref{sec:NG}. We conclude in \S \ref{sec:conclusions}.

\section{Perturbations and observables} \label{sec:actions}

In order to make observationally testable predictions for inflation driven by a  three-form field, we must understand the statistics of perturbations produced during the inflationary era. In particular, it is usual to focus on $\zeta$ \cite{Bardeen:1983qw,Lyth:1984gv}, the co-moving curvature perturbation, which is conserved on super-horizon scales if the dynamics is adiabatic \cite{Lyth:2004gb,Rigopoulos:2003ak}.  This property helps  to connect the statistics of the Fourier modes of $\zeta$ soon after they cross the horizon with observations of later epochs.

To this end, we wish to  perturb the action for the three-form theory (see \cite{Mukhanov:1990me,Malik:2008yp,Malik:2008im} for reviews of cosmological perturbation theory). The aim is to determine the form of the autonomous action for $\zeta$ in the setting of three form-inflation. Then by following the standard procedure of quantising $\zeta$ and fixing initial conditions far inside the cosmological horizon, the same-time correlation functions of $\zeta$, produced during inflation, can be calculated using the In-In formalism.   The second order action is required to calculate the two-point correlation function, parametrised by the power spectrum (see for example \cite{Lidsey:1995np} for a review), while the three-point correlation function, parametrised by the bispectrum, requires knowledge of the third order interactions in the theory, and hence the third order action \cite{Maldacena:2002vr, Seery:2005wm, Chen:2006nt}.

The procedure to calculate the action is to perturb the metric and the three-form, and insert the result into Eq.~\eqref{3formS}. Using various constraint equations, it is then possible to write the result purely in terms of only the perturbed quantity $\zeta$, and background quantities. The action can then be separated out order by order. We now present the details of this calculation.

\subsection{The perturbed action for three form inflation}

Following the approach of Maldacena \cite{Maldacena:2002vr}, we consider the Arnowitt-Deser-Misner form of the metric
\beta
{\rm d} s^2 = -N^2 {\rm d} t^2 + h_{ij} \left ( {\rm d} x^i + N^i {\rm d} t \right) \left( {\rm d}x^j + N^j {\rm d} t \right),
\end{equation}
where $N$ is the lapse function, $N_i$ the shift vector, and $h_{ij}$ the spatial metric. In this notation the action for the three-form takes the form
\beta
\label{ADMaction}
S=\frac{1}{2}\int {\rm d} t {\rm d} x^3 \sqrt{h} N \left( R^{(3)} - \frac{1}{24} F^2 - 2 V(A^2) \right )+
\frac{1}{2} \int {\rm d} t {\rm d} x^3 \sqrt{h} N^{-1} \left( E_{ij} E^{ij} -E^2 \right )\,,
\end{equation}
where $h={\rm det} h_{ij}$, $R^{(3)}$ is the Ricci scalar associated with the spatial metric,
and $E_{ij} =1/2 \dot{h}_{ij} - N_{(i|j)}$ is the extrinsic curvature. Considering only scalar perturbations,
when perturbed about an FRW metric the spatial metric takes the form
\beta
h_{ij} = a^2 e^{2 \zeta} \delta_{ij}\,,
\end{equation}
where we have used the spatial gauge freedom available to make the spatial metric diagonal, and where $\zeta$ denotes the curvature
perturbation.
The lapse function and shift vector are also perturbed to give
\beta
N_i= \psi_{,i} + \tilde{N}_i, ~~~~ N = 1+ \tilde{\alpha}\,,
\end{equation}
where $\tilde{N}_{i,i}=0$ and $\psi$ and $\tilde \alpha$ can be written as an expansion in powers of $\zeta$ as $\tilde \alpha= \tilde \alpha_1 + \tilde \alpha_2$ and $\psi = \psi_1 +\psi_2$.

We must also account for perturbations in the matter sourcing the Universe's evolution, which in this case is the
three-form field. The most general form of the perturbed $3$-form is
\beta
A_{0ij} = a(t)\epsilon_{ijk} (\alpha_{,k}+\alpha_k), ~~~~ A_{ijk} = a(t)^3\epsilon_{ijk}(\chi(t) + \alpha_0)\,.
\end{equation}
The vector perturbation $\alpha_k$ together with the vector perturbation in the metric (which we did not write down) decouple from scalar perturbations in
the second order action and also decay as usual \cite{KN2}. At higher
order they can also be neglected since our interest here is in the interaction terms which contribute to the scalar-scalar-scalar three-point function.
The metric tensor perturbation also decouples in the second order action and can be neglected at higher order for the same reason.
There are two scalar degrees of freedom, $\alpha$ and $\alpha_0$, which just like the metric perturbations can be expanded in powers of $\zeta$.

All the scalar perturbed quantities for the both the metric and the matter are then substituted into Eq.~\eqref{ADMaction}, which is expanded to third order.

\subsection{Fixing the temporal gauge}
\label{gaugechoice}
Thus far in our calculation we have only fixed the spatial part of the gauge freedom. At first order in perturbation theory,
under a general gauge transformation, $x^\mu \to x^{u} + \xi^{\mu}$ where $\xi^\mu =  ( \xi^0,\xi_{,i} )$ the scalar three-form
perturbations transform as
\beta
\alpha_0 \to \alpha_0-\dot{K} \xi^0+X \nabla^2 \xi\,, ~~~~ \alpha \to \alpha - a X \xi\,.
\end{equation}
The freedom associated with $\xi$ is already fixed, but not that associated with $\xi^0$. Considering the transformations above, it is clear that we could fix the gauge by setting $\alpha_0$ to zero, which was the choice made by De Felice, Karwan \& Wongjun \cite{DeFelice:2012jt}. This can also be done at second order.

However, it is important to note that this does not correspond to the choice which is commonly made,
for example by Maldacena \cite{Maldacena:2002vr}, Seery \& Lidsey \cite{Seery:2005wm} and Chen et al. \cite{Chen:2006nt},
who select the comoving gauge, defined at linear order by the condition
\beta
\delta T^0_i =0\,.
\end{equation}
In models sourced by a scalar field, this is equivalent to
setting $\delta \phi$ to zero, while in the present case it implies the following condition
\beta
\label{gauge}
\alpha = a^2 \chi \psi \,.
\end{equation}
Enforcing this condition at all orders fixes our temporal gauge freedom, and we can use this relation to replace $\alpha$ with $\psi$, or visa versa.
We will later  confirm explicitly that this choice is consistent with that made in work on scalar field inflationary models.

Before moving on with our calculation, we note that in their study Koivisto \& Nunes \cite{KN1,KN2} worked in the longitudinal gauge,
but wrote down the action for the curvature perturbation in the comoving gauge, $\zeta$. The choice of
De Felice, Karwan \& Wongjun \cite{DeFelice:2012jt} selects a subtly different gauge, and
explains why the second order actions in these two works do not coincide\footnote{It is important to note, that there is no problem with picking a gauge other than the comoving gauge and then forming $\zeta$ out of quantities in the chosen gauge at the end of the calculation. For example  in Refs \cite{Maldacena:2002vr,Seery:2005gb} the consistency of perturbing the action in the flat gauge and then calculating the statistics of $\zeta$ was shown.}.

\subsection{The second order action}

Moving forward with our calculation, we use the condition \eqref{gauge} to substitute for $\alpha$ in the action. We find a rather cumbersome expression containing terms in $\psi$, $\alpha_0$, $\tilde \alpha$ and $\zeta$. The first order part of the action is zero once the background equations \eqref{background} of motion are utilised (or equivalently varying the first order part of the action leads to the background equations of motion).  Our aim now is to find substitutions for $\psi$, $\alpha_0$ and $\tilde \alpha$ in terms of $\zeta$, and hence to obtain an action only in terms of $\zeta$. As long as these substitutions follow from a constraint equation, and do not alter the order of the action, we do no harm by using them to write the action in terms of $\zeta$ alone. First we focus on the second order part of the action.

In the present case the simplest way to proceed is to utilise four constraint equations, which are not independent of each other, for the three quantities we wish to eliminate. The first is the equation of motion for $N$, or equivalently $\tilde \alpha$
\begin{eqnarray}
\label{eqmN}
&~&3 H (\dot\zeta-H\tilde\alpha) -\frac{1}{2} (\dot \chi + 3H\chi) (\dot \alpha_0 + 3H\alpha_0) +\frac{1}{2} (\dot\chi+3H\chi)^2 (\tilde \alpha + 3\zeta) \nonumber \\
&~&-\frac{1}{2} V_{,\chi}\alpha_0 + \frac{3}{2} \chi V_{,\chi} \zeta
+\frac{1}{2} a^2 \chi (\dot \chi + 3 H \chi) \frac{\partial^2}{a^2} \psi - \frac{\partial^2}{a^2} (H\psi+\zeta) = 0.
\end{eqnarray}
The next is the equation of motion for $N_i$, or equivalently $\psi$
\beta
\tilde \alpha = \frac{\dot \zeta}{H} .
\label{eqpsi}
\end{equation}
Then we can also use the equation
which follows from the anti-symmetry of the four-form, and an equation which is most easily seen
from the conservation of the energy momentum tensor. These are
\begin{eqnarray}
\label{antisymm}
\alpha_0 &=& 3 \chi \zeta - \frac{V_{,\chi}}{V_{,\chi\chi}} \tilde \alpha, \\
\label{consTmn}
\dot \alpha_0 &=& a^2 \chi \frac{\partial^2}{a^2}\psi + (\dot\chi + 3 H \chi) (\tilde \alpha + 3 \zeta) - 3 H \alpha_0 ,
\end{eqnarray}
respectively. One can verify that all these equations are self-consistent both with each other and
with the variation of the action with respect to $\tilde \alpha$, $\psi$ and $\alpha_0$.

We find that the substitutions which consistently satisfy the equations above are
\begin{eqnarray}
\label{subalpha}
\tilde\alpha &=& \frac{\dot \zeta}{H}, \\
\label{subpsi}
\psi &=& -\frac{\zeta}{H} + \sigma, \\
\label{subsigma}
\partial^2 \sigma &=& a^2 \frac{\epsilon}{c_s^2} \dot\zeta,
\end{eqnarray}
with $\epsilon = -\dot{H}/H^2 = \chi V_{,\chi}/2 H^2$ and the speed of sound is  given by \cite{KN2}
\begin{equation}
\label{csound}
c_s^2 = \frac{V_{,\chi \chi} \chi}{V_{,\chi}}.
\end{equation}

We note that these substitutions are only first order, while we might expect the second order action to include terms which come from $\tilde \alpha_2$ etc..  As highlighted by Maldacena \cite{Maldacena:2002vr}, however, that these relations come from a constraint equation is advantageous, since it means that only the first order term in the substitutions for the perturbed quantities in terms of $\zeta$ is required. The second order part of each term (for example $\tilde \alpha_2$) will multiply the zeroth order constraint in the second order action. Looking ahead to the third order action, we note that even there we will only need the first order substitution. This is because the second order parts  multiply the first order constraint in the third order action, and so disappear. While the third order part will multiply the zeroth order constraint in the third order action. That this occurs can be verified explicitly.

The next step is to make all these substitutions in the second order part of the action so as to write it purely in terms of $\zeta$. After proceeding in this way, utilising the background equations \eqref{background} and  integrating by parts in a suitable manner we arrive at the action
\beta
\label{action2a}
S_2= \int {\rm d}t {\rm d}^3x \left [a^3 \frac{\Sigma}{H^2}\dot\zeta^2 - a \epsilon (\partial \zeta)^2 \right ]\,,
\end{equation}
with
\beta
\label{sigma}
\Sigma=\frac{H^2 \epsilon}{c_s^2}.
\end{equation}
We note that in this calculation the $\psi$ substitution is not actually used, as $\psi$ naturally cancels out of the second order action after the other substitutions. This second order action matches precisely that derived by Koivisto \& Nunes \cite{KN2}.

\subsection{The third order action}

Now we wish to isolate the third order part of the action in terms of $\zeta$ alone. Once again we can make the substitutions above in the third order part. After numerous integrations by parts we arrive at
\begin{eqnarray}
\label{action3a}
S_3&=&\int \left \{ dt d^3x [-\epsilon a \zeta (\partial \zeta)^2 - a^3 (\Sigma + 2 \lambda) \frac{\dot{\zeta}^3}{H^3}\right.+ \frac{3 a^3 \epsilon}{c_s^2} \zeta \dot{\zeta}^2 \nonumber\\
	&+&\left. \frac{1}{2 a}
\left ( 3 \zeta - \frac{\dot{\zeta}}{H}\right ) \left ( \partial_i \partial_j \psi \partial_i \partial_j\psi - \partial^2 \psi \partial^2 \psi\right) -2 a^{-1} \partial_i\psi\partial_i\zeta \partial^2 \psi\right \}
\end{eqnarray}
where
\beta
\label{lambda}
\lambda = -\frac{1}{12} \frac{V_{,\chi}^3 V_{,\chi\chi\chi}}{V_{,\chi\chi}^3}\,.
\end{equation}
and where for simplicity we have not explicitly substituted for $\psi$, but could do so using Eq.~\eqref{subpsi}.

At this juncture we pause to make an observation. The form of both the second and third order action, is precisely that found for a $P(X, \phi)$ action of the form \eqref{eq:Paction}, except that in that case, $c_s^2$ and $\lambda$ are written in terms of derivatives of $P(X,\phi)$ with respect to $X$ \cite{Seery:2005wm}. One might have expected this, given the formal equivalence of the theories. However, as we have noted previously, the mapping between them is rather complex, and working directly with the scalar dual for a three-form model with a given potential is often cumbersome or indeed impossible, if the analytic form of the inverse potential is not known\footnote{This can be the case, for example, if the potential is a polynomial of degree five or above.}. We will see below, however, that using the duality in a more formal manner allows us to use the $P(X,\phi)$ description to recover Eqs. \eqref{csound} and \eqref{lambda}, effectively providing an alternative derivation for the perturbed three-form action to third order, and providing a powerful confirmation of our results thus far. We now proceed to that calculation.

\section{Dual theories} \label{sec:dualities}

In this section we review a number of dualities which inter-relate $p$-form theories in four space-time dimensions. These have been discussed elsewhere \cite{Germani:2009iq,Germani:2009gg, Koivisto:2009sd,KN2}, and our aim here is to employ them as a practical tool. To this end we show how one of these dualities can be used to confirm our earlier results for the perturbed action of a three-form theory, utilising previous results derived for a non-canonical scalar field theory by other authors \cite{Seery:2005wm,Chen:2006nt}.

\subsection{Equations of motion for the three-form}

We begin by considering a Palatini-type action, which treats the three-form $A$ and the four-form $F$ as independent variables, and which  is equivalent to \eqref{3formS} up to boundary terms. The Lagrangian for the matter part of the action is given by
\beta \label{parent}
{\cal L}_{1} = \frac{1}{48}F^2 - \frac{1}{6}A \nabla \cdot F - V(A^2)\,.
\end{equation}
Variations of the associated action relate $A$ and $F$.  The resulting equations of motion are
\begin{eqnarray}
F &=& -4 \left[ \nabla A \right] \nonumber \\
\nabla \cdot F &=& -12 A V'\left(A^2\right), \label{eom1}
\end{eqnarray}
where a prime denotes differentiation with respect to the argument in brackets. Integrating the middle term of ${\cal L}_1$ by parts inside the action, one finds
\beta \label{ibp}
{\cal L}_{2} = \frac{1}{48}F^2 + \frac{1}{6}F \left[ \nabla A \right] - V(A^2),
\end{equation}
which consequently shares the same equations of motion. The first equation of  \eqref{eom1} now appears as a constraint equation and may be substituted back into ${\cal L}_{2}$ to confirm that we do indeed recover~\eqref{3formS}. We note that we are always free to perform such an integration by parts and this procedure leaves the dynamics of the theory invariant. Moreover, we are free to substitute constraint equations  which do not change the order back into the action.

\subsection{Dual actions}

Our primary aim is to rewrite Lagrangians~\eqref{parent} and~\eqref{ibp} in terms of the Hodge ($*$) dual fields to $A$ and $F$. In appendix \ref{appendix1}, for clarity, we give some pedagogical detail regarding such duals. We recall that any $p$-form has a dual ($d-p$)-form, where $d$ is the number of space-time dimensions, four in our context. In particular, the three-form $A$ and four-form $F$ that make up~\eqref{ibp} can be expressed in terms of their duals as
\begin{align}
\label{hodge}
&(\star F) = \frac{1}{4!}\epsilon_{\alpha\beta\gamma\delta}F^{\alpha\beta\gamma\delta} \equiv \Phi   &F_{\alpha\beta\gamma\delta} = - \epsilon_{\alpha\beta\gamma\delta} \Phi \nonumber\\
&(\star A)_{\alpha} = \frac{1}{3!} \epsilon_{\alpha\beta\gamma\delta} A^{\beta\gamma\delta} \equiv B_\alpha &A_{\beta\gamma\delta} = - \epsilon_{\alpha\beta\gamma\delta} B^\alpha
\end{align}

The Hodge duals to $F$ and $A$ enable us to recast the original theory \eqref{3formS} into a scalar-vector description, with Lagrangian
\begin{eqnarray} \label{scalarvector}
{\cal L}_{3} &=& -\frac{1}{2}\Phi^2 - \Phi \nabla \cdot B - V\left(- 6 B^2\right),
\end{eqnarray}
which follows from Lagrangian~\eqref{ibp}.  The equations of motion for $\Phi$ and $B_\mu$ can now be obtained either by varying~\eqref{scalarvector} or equivalently by substituting Hodge duals into~\eqref{eom1}. They are
\begin{eqnarray}
\Phi &=& - \nabla \cdot B , \nonumber \\
\nabla \Phi &=& - 12 B V'(-6B^2).
\label{eq:eomsBPhi}
\end{eqnarray}
The first equation of motion now appears as a constraint equation with respect to ${\cal L}_{3}$ and can be substituted  back to express our original theory~\eqref{parent} as a pure vector theory. Integrating the middle term of ${\cal L}_{3}$ by parts (or equivalently substituting $A$ and $F$ for their duals in ${\cal L}_{2} $), however, one finds that the converse is true, and the second equation of motion in \eqref{eq:eomsBPhi} appears as a constraint and may be substituted into the action to remove the vector field in favour of the scalar (or equivalently its dual four form). With an eye on calculating inflationary observables, and in particular the 3-point correlation function, the scalar picture is particularly intriguing, since, as we have discussed, it opens up the possibility of using existing machinery for  dealing with scalar field models of inflation. The final set of equivalent actions for a four-form, three-form, vector and scalar respectively are
\begin{eqnarray} \label{dualtheories}
{\cal L}_{IV}(F,\nabla \cdot F) &=& -\frac{1}{48}F^2 + 2 A^2(\nabla \cdot F) V' \left(A^2(\nabla \cdot F)\right) - V(A^2(\nabla \cdot F)), \\
{\cal L}_{III}(A, \nabla A)) &=& -\frac{1}{3}\left[ \nabla A \right]^2 - V(A^2), \\
{\cal L}_{I}(B,\nabla \cdot B) &=& \frac{1}{2}(\nabla \cdot B)^2  - V\left(-6 B^2\right), \\
{\cal L}_{0}(\Phi, \nabla \Phi) &=& -\frac{1}{2}\Phi^2 -12 B^2(\nabla \Phi) V' \left(-6B^2(\nabla \Phi) \right) - V\left(-6 B^2(\nabla \Phi)\right).
\end{eqnarray}
Here $V\left(-6 B^2(\nabla \Phi)\right)$ and $V(A^2(\nabla \cdot F))$ indicate that the second equation of motion in \eqref{eq:eomsBPhi} has to be used in order to express $A^2$ in terms of $\nabla \cdot F$ and $B^2$ in terms of $\nabla \Phi$. Note that we cannot simultaneously use the first equation of motion ($\Phi = - \nabla \cdot B $) to substitute for $-\Phi \nabla \cdot B$ in the action in order to replace it with $\Phi^2$, since this would change the order of the action.

It is interesting to examine the form of these dual theories. The potential for $A/B$ essentially gets mapped into a non-canonical kinetic term for the $F/\Phi$  (4-form/scalar) theory respectively. The canonical kinetic terms in the $A/B$ picture, on the other hand, give rise to simple quadratic potential terms in the corresponding 4-form/scalar theories. This is important for several reasons. First, in this way an effective non-canonical scalar theory arises from a very simple three-form theory.\footnote{The reader might have noticed that there are differences in factors and signs when comparing~\eqref{ibp} and~\eqref{dualtheories} with corresponding equations presented in~\cite{KN2}. The version presented here corrects a small number of typographical mistakes in that work.}
Second, this immediately tells us that
all scalar models dual to the three-form share
the same simple quadratic potential. This is particularly important for standard
slow-roll inflation\footnote{By this we here mean $\epsilon,\eta... \ll 1$ as well as requiring no rapidly varying speed of sound.}, where in the dual scalar picture the potential dominates over the kinetic terms. The
fact that all
models share the same dual scalar potential then implies that the form of the original
three-form potential (which turns into a non-canonical kinetic term) is not important
when computing, for example, the spectral
index $n_s$. We will return to this point later.

Finally, one may wonder how one can come up with an effective single scalar theory dual to a three-form theory, which in principle possesses more physical degrees of freedom. However, starting with the most general canonical and minimally coupled 3-form action in 4$d$, as we do here, guarantees that such a dual single scalar field description always exists. This is the case, because 1) the canonical kinetic term for the three-form dualises to a $\Phi^2$ potential 2) the three-form potential is a function of $A^2$ only, because in 4d this is the only covariant scalar combination that can be built from a 3-form - thus the potential only depends on one effective degree of freedom: $A^2$ and 3) these two degrees of freedom, $\Phi$ and $A^2$, are related via an equation of motion, leaving only one effective independent degree of freedom, thus explaining the existence of a dual single scalar description.

\subsection{The perturbed action from scalar duality}

Now we wish to show that one can derive Eqs. \eqref{csound} and \eqref{lambda} using the results above. In particular, we employ the expressions
\beta
B^2 = -\chi^2, \hspace{1cm} X \equiv -g^{\mu\nu} \nabla_\mu \phi \nabla_\nu \phi = -12^2 B^2 \left(V'(-6B^2)\right)^2 \,,
\end{equation}
which follow from the definition of the Hodge dual to $A$, Eq.~\eqref{hodge}, and from Eq.~\eqref{eq:eomsBPhi} respectively, to find
\beta
\label{derv}
X = V_{,\chi}^2, \hspace{1cm} {\rm and} \hspace{1cm}
\frac{\partial \chi}{\partial X} = \frac{1}{2 V_{,\chi} V_{,\chi \chi}}\,.
\end{equation}
From the discussion of dualities, and in particular ${\cal L}_{0}$, we see that the $P(X,\phi)$ theory dual to our minimally coupled three-form theory was parametrised by
\begin{eqnarray}
P(X,\phi) &=& - \frac{1}{2} \phi^2 -12B^2 V'(-6B^2) -V(-6B^2)\nonumber \\
&=&  -\frac{1}{2} \phi^2 + \chi V_{,\chi} - V(6\chi^2).
\end{eqnarray}
Differentiating this expression with respect to $X$ and using \eqref{derv} we find $P_{,X}$, $P_{,XX}$ and $P_{,XXX}$ which upon substitution into
\beta
\label{csTF}
c_s^2 = \frac{P_{, X}}{P_{, X}+2  X P_{,X X}}, \hspace{1cm} \lambda = X^2 P_{,XX} +\frac{2}{3} X^3 P_{,XXX}\,,
\end{equation}
which are the expressions for $c_s^2$ and $\lambda$ for a $P(X,\phi)$ theory \cite{Seery:2005wm}, confirms the expression for the speed of
sound \eqref{csound} and for $\lambda$ \eqref{lambda} we found for the three-form theory.

This constitutes the result we are looking for. If we had only known about the perturbed action in terms of $\zeta$ from a $P(X,\phi)$ theory, we could have used the argument above to write that perturbed action   solely in terms of three-form quantities, hence providing an alternative route to our earlier result.

The power of this approach is that we can also readily probe higher order statistics, such as the trispectrum, since the quartic action has already been calculated for $P(X,\phi)$ theories \cite{Huang:2006eha,Arroja:2008ga}. This involves the quantity $\Pi$, defined in Ref.~\cite{Arroja:2008ga}, which for completeness, we calculate to be
\beta
\Pi = \frac{1}{40} \frac{V_{,\chi}^3 V_{,\chi\chi\chi}}{V_{,\chi\chi}^3}
+ \frac{3}{40} \frac{V_{,\chi}^4 V_{,\chi\chi\chi}^2}{V_{,\chi\chi}^5}
- \frac{1}{40} \frac{V_{,\chi}^4 V_{,\chi\chi\chi\chi}}{V_{,\chi\chi}^4},
\end{equation}
in terms of the three-form quantity $\chi$.

\subsection{Confirmation of our gauge choice}
Finally we would like to confirm that the gauge we have chosen in \S \ref{gaugechoice} is equivalent to the gauge chosen when the perturbed
action for the $P(X,\phi)$ theory is calculated \cite{Seery:2005wm, Chen:2006nt}. There the choice is usually expressed as $\delta \phi = 0$. Utilising the expression
\beta
X=12 \chi^2 \left(V'(6\chi^2)\right)^2\,,
\end{equation}
and perturbing both sides, we find this leads to the expression
\beta
\alpha_0 = 3 \chi \zeta - \frac{V_{,\chi}}{V_{,\chi\chi}} \tilde\alpha,
\end{equation}
which is the same as the anti-symmetry condition \eqref{antisymm} after we fixed the gauge choice $\alpha = a^2 \chi \psi$. It is then possible to verify that the symmetry constraint on the three-form gives Eq. \eqref{gauge}, and the other substitutions remain unchanged. In this way the equivalence of our gauge choice with the usual choice made in scalar field theories is confirmed, in contrast with the choice made in other works \cite{DeFelice:2012jt}.

\section{Non-Gaussianities}
\label{sec:NG}
\subsection{Correlation functions}

We will now proceed to utilise our results thus far to compute  observables for the three-form inflationary theory under consideration. Because our action is identical to that for a non-canonical scalar field (except with $c_s^2$ and $\lambda$ being expressed in terms of background three-form quantities), the calculation follows precisely that of Garriga \& Mukhanov \cite{Garriga:1999vw} to calculate the power spectrum, and  Seery \& Lidsey \cite{Seery:2005wm} and Chen et al. \cite{Chen:2006nt} to calculate the three-point function. Here we simply provide an overview of the main steps and the important results.

The 2-point correlation function for the curvature perturbation $\zeta$ is defined as
\begin{equation} \label{2point}
\langle \zeta(\textbf{k}_1)\zeta(\textbf{k}_2)\rangle = (2\pi)^5
\delta^3({\bf k}_1+{\bf k}_2)  \frac{P_\zeta}{2 k_1^3}\,,
\end{equation}
and is calculated using the second order action Eq.~\eqref{action2a}. As found by \cite{Garriga:1999vw} the end result is
\begin{equation}
P_\zeta \equiv \frac{1}{2\pi^2}k^3\left\vert\zeta_k\right\vert^2 =\left. \frac{1}{2 (2 \pi)^2\epsilon}\frac{{ H}^2}{{ c}_sM_{\rm Pl}^2}\right|_*\,,
\label{zetaPk}
\end{equation}
where $*$ indicates that the expression is evaluated at horizon crossing $c_s k = a H$. The spectral index $n_s$ is then given by
\beta \label{ns}
1-n_s = 2\epsilon +\frac{\dot \epsilon}{\epsilon H} +\frac{\dot c_s}{c_s H}\,.
\end{equation}

The slow-roll approximation (i.e. neglecting $\ddot\chi$ in the equation of motion for $\chi$) implies that
\beta
\label{chiN}
\frac{\chi_{,N}}{\chi} \approx -\frac{2}{3}\frac{1}{\chi^2}\left(1-\frac{3}{2}\chi^2\right)\epsilon .
\end{equation}
Substituting this into the Friedmann equation we obtain \cite{KN1}
\beta
\label{epsiloneq}
\epsilon \equiv - \frac{\dot H}{H^2} \approx \frac{3}{2} \chi \frac{V_{,\chi}}{V} \left(1- \frac{3}{2}\chi^2\right) .
\end{equation}
Eliminating the term in brackets from the last two equations we can write
\beta
\frac{\chi_{,N}}{\chi} = -\frac{4}{9}\frac{1}{\chi^2} \frac{V}{\chi V_{,\chi}} \epsilon^2 .
\end{equation}
This expression allows us to easily compute $\dot \epsilon$ from the definition of $\epsilon$ to find $\dot\epsilon/\epsilon H = 2\epsilon + {\cal O}(\epsilon^2)$. Similarly, we find that $\dot c_s/c_s H \approx
{\cal O}(\epsilon^2)$
and also $\dot \lambda/\lambda H \approx {\cal O}(\epsilon^2)$.
This means that to first order in the slow-roll parameters, the scalar spectral index is simply $n_s = 1 - 4\epsilon$.

The ratio of tensor to scalar perturbations was computed in Ref.\cite{KN2} and shown to be related to the slow roll parameter $\epsilon$ as\footnote{A more exact answer, fully taking into account that tensor and scalar modes freeze out at different times for models with $c_s \ne 1$, can be found in \cite{Agarwal:2008ah}. Here we just note that this extra effect means that the actual tensor-to-scalar ratio $r$ is in fact smaller than naively expected from $r = 16 c_s \epsilon$ for $c_s < 1$.}
\beta
\label{eqr}
r = 16 c_s \epsilon,
\end{equation}
mirroring the analogous expression for non-canonical scalar field models. One can also compute the spectral index for tensor perturbations $n_t$, finding $n_t = -2\epsilon$ to first order in slow-roll parameters as usual.

The third order action \eqref{action3a} is needed to calculate the three-point correlation function.
First, however, it must be written in a suitable form, which follows from the use of the $\psi$ substitution, and is \cite{Seery:2005wm,Chen:2006nt}
\begin{eqnarray} \label{action3b}
S_3&=&\int {\rm d}t {\rm d}^3x \left\{
-a^3 \left[\Sigma\left(1-\frac{1}{c_s^2}\right)+2\lambda\right] \frac{\dot{\zeta}^3}{H^3}
+\frac{a^3\epsilon}{c_s^4}(\epsilon-3+3c_s^2)\zeta\dot{\zeta}^2 \right.
\nonumber \\ &+&
\frac{a\epsilon}{c_s^2}(\epsilon-2\epsilon_s+1-c_s^2)\zeta(\partial\zeta)^2-
2a \frac{\epsilon}{c_s^2}\dot{\zeta}(\partial
\zeta)(\partial \sigma) \nonumber \\ &+& \left.
\frac{a^3\epsilon}{2c_s^2}\frac{d}{dt}\left(\frac{\eta}{c_s^2}\right)\zeta^2\dot{\zeta}
+\frac{\epsilon}{2a}(\partial\zeta)(\partial
\sigma) \partial^2 \sigma+\frac{\epsilon}{4a}(\partial^2\zeta)(\partial
\sigma)^2+ 2 f(\zeta)\left.\frac{\delta L}{\delta \zeta}\right\vert_1 \right\},
\end{eqnarray}
where $f(\zeta) \delta L/\delta \zeta|_1$ indicates terms proportional to the functional derivative of the Lagrangian evaluated at first order in $\zeta$.
This would be zero if $\zeta$ was Gaussian. Such terms can be removed by a field redefinition (though one
must recall the redefinition when forming correlation functions so as to form corrections of $\zeta$ itself).
This form of the action thus identifies the relevant interaction vertices which contribute towards the
three-point function.

At tree level in quantum field theory, and in the interaction picture, the In-In (equal time)
three-point correlation function is given by the expression
\begin{equation} \label{interaction}
\langle
\zeta(t,\textbf{k}_1)\zeta(t,\textbf{k}_2)\zeta(t,\textbf{k}_3)\rangle=
-i\int_{t_0}^{t}{\rm d}t^{\prime}\langle[
\zeta(t,\textbf{k}_1)\zeta(t,\textbf{k}_2)\zeta(t,\textbf{k}_3),H_{\rm int}(t^{\prime})]\rangle ~,
\end{equation}
where $H_{\rm int}$ is the Hamiltonian evaluated at third order in the perturbations and follows directly
from~\eqref{action3b}. Vacuum expectation values are evaluated with respect to the interacting vacuum $|\Omega \rangle$.
By convention the 3-point correlation function is parametrised by the amplitude ${\cal A}$.
\begin{equation} \label{amplidef}
\langle \zeta(\textbf{k}_1)\zeta(\textbf{k}_2)\zeta(\textbf{k}_3)\rangle = (2\pi)^7
\delta^3({\bf k}_1+{\bf k}_2+{\bf k}_3) P_\zeta^{\;2} \frac{1}{\Pi_j k_j^3}{\cal A}\,,
\end{equation}
where, again by convention the power spectrum $P_\zeta$ in the above formula is calculated for the mode
$K = k_1 + k_2 + k_3$. Evaluating \eqref{interaction}, one can determine ${\cal A}$

In principle ${\cal A}$ is a general function of the three Fourier modes, which are related by the
condition $\textbf{k}_1 + \textbf{k}_2 + \textbf{k}_3=0$ (within the slow-roll approximation the full form of ${\cal A}$ is given in Ref.~\cite{Chen:2006nt}).
For a given shape of non-Gaussianity, however, (see figure~\ref{fig1}) the size of non-Gaussianity can
be adequately characterised by a single-value measure $f_{\rm NL}$. For an equilateral shape (i.e. one peaking in the limit $k_1 \sim k_2 \sim k_3$), this can be defined as~\cite{Khoury:2008wj}
\beta
f^{\rm equil}_{\rm NL} = 30\frac{{\cal A}_{k_1=k_2=k_3}}{K^3}\,,\label{fnl}
\end{equation}
where amplitudes are matched at $k_1=k_2=k_3=K/3$. Note that we follow the WMAP sign convention here, where positive $f_{\rm NL}$ physically corresponds to negative-skewness for the temperature fluctuations. The parameters controlling the overall size of $f_{\rm NL}^{\rm equil}$ are $c_s,\Sigma,\lambda$, which we now summarise for convenience:
\begin{align} \label{3formCSL}
&c_s^2 = \frac{V_{,\chi\chi} \chi}{V_{,\chi}},&
&\Sigma = \frac{V_{,\chi}^2}{2 V_{,\chi\chi}},&
&\lambda = - \frac{V_{,\chi}^3 V_{,\chi\chi\chi}}{12 V_{,\chi\chi}^3}.&
\end{align}
Following~\cite{Chen:2006nt} we may now compute $f^{\rm equil}_{\rm NL}$ in the slow-roll regime, finding the result at leading order to be\footnote{Note that \eqref{eqfnl} differs from the expression found by Chen et al.~\cite{Chen:2006nt} by an overall sign. This is because Chen et al. use a sign convention opposite to that used by WMAP and throughout this paper.}
\beta
\label{eqfnl}
f_{\mathrm{NL}}^{\rm equil} \approx \frac{5}{81}\left(\frac{1}{c_s^2}-1- 2 \frac{\lambda}{\Sigma}\right)
- \frac{35}{108} \left(\frac{1}{c_s^2}-1\right) +
 {\cal O}\left(\epsilon,\frac{\epsilon}{c_s^2}, \epsilon \frac{\lambda}{\Sigma},\frac{\dot \lambda}{\lambda H}\right).
\end{equation}

This expression clearly shows that models with small speed of sound can lead to a large non-Gaussian signal. In the context of three-forms, a hyperbolic potential was discussed in Ref. \cite{DeFelice:2012jt}, which yields a small speed of sound during inflation, though close to unity in the oscillatory phase. Here we explore two other simple potentials.

\subsection{Example I: A power-law potential} \label{subsec-power}

We now consider some concrete examples, and begin by considering a three-form model self-interacting through
a simple power-law potential
\beta \label{powerlawL}
{\cal L} = -\frac{1}{48}F^2  - V_0 A^{2p}\,,
\end{equation}
that is, $V(A^2) = V_0 A^{2p} = V_0 \left(6 \chi^2\right)^p$, where $p$ is a constant.

This is a special example, where the equivalent $P(X,\phi)$ theory is relatively simple, and is given by
\beta
{\cal L}_{\phi} = (2p - 1) \left(\frac{1}{V_0}\right)^{1/(2p-1)}\left( \frac{X}{24 p^2}\right)^{\frac{p}{2p-1}} - \frac{1}{2}\phi^2,
\end{equation}
where one now sees explicitly that the 3-form potential has been mapped into a non-canonical kinetic term for the effective scalar $\phi$.
As long as $2p - 1 \ne 0$, one can check that equations~\eqref{csound}, \eqref{sigma} and~\eqref{csTF} equivalently yield
\begin{eqnarray}
c_s^2 &=& 2p - 1 \,,\nonumber \\
\frac{\lambda}{\Sigma} &=&  \frac{1}{3}\frac{1 - p}{2p-1}\,, \label{csls}
\end{eqnarray}
and we can immediately make observational predictions for the theory, finding several interesting results.

Here both the speed of sound and $\lambda/\Sigma$ are constant, and we have $n_s -1 = -4\epsilon$. The constancy of the speed of sound has the interesting consequence that the observational requirement of obtaining a (near)-scale-invariant spectral index $n_s$ forces $\epsilon$ to be close to zero and slowly-varying here (compare \eqref{ns} and also \cite{Khoury:2008wj}). For this power-law model, slow-roll therefore becomes an observational requirement in contrast to models with varying speed of sound.

We now need to calculate this quantity $N$ $e$-folds before the end of inflation.
It was shown in \cite{KN1} that for a power law potential, the value of the field at this time can be estimated to be
\beta
\chi_N^2 = \frac{2}{3} - \frac{4}{18 p} \frac{1}{1+2N},
\end{equation}
which upon substitution in \eqref{epsiloneq} gives
\beta
\epsilon_N \approx \frac{1}{1+2N}.
\end{equation}
Assuming that $N\approx 60$ is required to solve the horizon problem, we predict that the spectral index on observationally
relevant scales for a three-form with power law potential is
\beta
n_s \approx 0.97,
\end{equation}
and independent of the value of the power $p$. Naturally the spectral index will be closer to scale-invariant the longer inflation lasts.
Substituting for $c_s^2$ and $\lambda/\Sigma$ into \eqref{eqfnl} we obtain the dependence of $f_{NL}^{\rm equil}$ on $c_s^2$
illustrated in Fig.~\ref{tejo1fnl}. It clearly  shows that a small speed of sound leads to a large $f_{NL}^{\rm equil}$ as expected.
\begin{figure}[tp]
\begin{center}$
\begin{array}{c}
\includegraphics[width=0.7\linewidth]{fnlPower.png}
\end{array}$
\end{center}
\caption{
Dependence of $f_{\rm NL}^{equil}$ on ${c}_s^2$ for the power law potential $V = V_0 (6\chi^2)^p$ and $N = 60$. A large and generically negative non-Gaussian amplitude is found. The shaded region is disallowed by the WMAP 2$\sigma$ bound $-214 < f_{\rm NL}^{equil} < 266$.\cite{Komatsu:2010fb} We recall that $c_s^2 = 2p - 1$ for this model.
\label{tejo1fnl}}
\end{figure}
Substituting for the speed of sound for a power law potential in Eqs.\eqref{eqr} and \eqref{eqfnl} we can relate $f_{NL}^{\rm equil}$ and $r$. Figure~\ref{fnlr} illustrates how these two quantities are related, and shows the region of the parameter space ($r,f_{NL}$) that is allowed given current bounds.
\begin{figure}[tp]
\begin{center}
\includegraphics[width=0.7\linewidth]{fnlvsrPower.png}
\caption{\label{fnlr} The solid lines show how the parameters $f_{NL}^{\rm equil}$ and $r$ are related to each other for the power law potential, $V = V_0 (6\chi^2)^{p}$ and for $N = 50, 60, 80$ from top to bottom. The shaded region is disallowed by the WMAP 2$\sigma$ bound $-214 < f_{\rm NL}^{equil} < 266$.\cite{Komatsu:2010fb}}
\end{center}
\end{figure}

Phenomenologically we find that no large non-Gaussian amplitude of the enfolded or orthogonal types (which have peaks in the folded limit $k_1 \sim k_2 \sim 2 k_3$) can be generated here in contrast to generic single field inflation models. To see why, it is useful to notice that one may express $\lambda/\Sigma$ as \cite{Seery:2005wm,Khoury:2008wj}
\beta \label{fX}
\frac{\lambda}{\Sigma} = \frac{1}{6} \left( \frac{2 f_X + 1}{c_s^2} - 1 \right)
\end{equation}
where
\beta
f_X = \frac{\epsilon \epsilon_s}{3 \epsilon_X}, \quad\quad\quad \epsilon_s = \frac{\dot{c_s}}{H c_s}, \quad\quad\quad \epsilon_X = -\frac{\dot{X}}{H^2}\frac{\partial H}{\partial X}.
\end{equation}
Now, since the speed of sound is constant for the power-law model, $\epsilon_s, f_X$ are identically zero here. A large and predominantly orthogonal (or enfolded) amplitude, however, requires a negative and non-zero $f_X$ (assuming positive $\epsilon$ and $0 \leq c_s \leq 1$)\cite{Senatore:2009gt}\footnote{In general $P(X,\phi)$ inflation models a predominantly orthogonal or folded shape can be generated by finely balancing the contributions from the $\dot{\zeta}^3$ interaction vertex (which depends on $\lambda/\Sigma$) against the other vertices such that the generically predominant equilateral shape contributions cancel out. In this way the otherwise subdominant orthogonal or folded configurations are brought out. As stated above, this delicate cancellation relies on $f_X$-dependent contributions though, whereas $f_X = 0$ identically for the power-law example given.}.
As such the non-Gaussian shape found here is always predominantly equilateral. Furthermore, since parameters $\lambda$ and $\Sigma$ are simply related to $c_s$ via \eqref{fX}, the non-Gaussian amplitude is completely controlled by $c_s^2$ for a power-law potential. This means that $f_{\rm NL}^{\rm equil} \sim {\cal O}(c_s^{-2})$ for subluminal and positive speed of sound $c_s$, i.e. $p \; \epsilon \left[1/2,1\right]$. Figure~\ref{fig1} illustrates these points. Finally note that, analogous to e.g. DBI inflation, any sizable level of non-Gaussianity has a negative $f_{\rm NL}^{\rm equil}$ associated with it here.
\begin{figure}[tp]
\begin{center}$
\begin{array}{c}
\includegraphics[width=0.7\linewidth]{shapePower.png} \\
\end{array}$
\end{center}
\caption{
The dimensionless bispectrum ${\cal A}(k_1,k_2,k_3)/k_1 k_2 k_3$ for the power law potential, $V = V_0 (6\chi^2)^{p}$ with $p = 0.505$, corresponding to $c_s = 0.1$, and with $N = 60$. Triangular shapes denote the equilateral, squeezed/local and enfolded limit clockwise from top right. A predominantly equilateral shape is found.
\label{fig1}}
\end{figure}

\subsection{Example II:  Exponential potential}

Let us now consider an alternative potential. Even though it is fairly simple in its analytical form, the exponential potential, $V = V_0 \exp(\beta A^2)  = V_0 \exp(6 \beta \chi^2)$, becomes very complex when written as a $P(X,\phi)$ theory,
\beta
{\cal L}_0 = \left(W(x) -1\right) V_0 \exp\left(\frac{1}{2} W(x)\right) - \frac{1}{2}\phi^2,
\end{equation}
where $W(x)$ is the Lambert-W function and $x = X/12\beta V_0^2$. Dealing with this model in the $P(X,\phi)$ description is therefore very difficult. This is a good example showing that performing the calculations in the original three-form theory is far simpler than going to the $P(X,\phi)$ description.

Inflation ends when $\epsilon \approx 1$ which for this potential takes place when
\beta
\chi_e^2 = \frac{1}{3} + \frac{1}{3} \left(1-\frac{1}{3\beta}\right)^{1/2} .
\end{equation}
This expression implies that inflation only ends if $\beta >1/3$. In particular, if $\beta \gg 1/3$ we can approximate
\beta
\chi_e^2 \approx \frac{2}{3} - \frac{1}{18\beta}.
\end{equation}
In this case, the solution of \eqref{chiN} and \eqref{epsiloneq} is non-analytical and we must resort to an approximation. Defining $\chi^2 = 2/3 - y$, since we know that inflation occurs very close to $2/3$, we obtain
\beta
\chi_N^2 = \frac{2}{3} - \frac{1}{18\beta} \, \frac{1}{1+\sqrt{6} N}.
\end{equation}
The slow-roll parameter $N$ $e$-folds before the end of inflation is then $\epsilon_N \approx 1/(1+\sqrt{6}N)$. Here, $c_s^2 = 1+ 12 \beta \chi^2$ and it can be verified that $\dot{c}_s/c_sH  = -\epsilon^2/\chi^2 \sim {\cal O}(\epsilon^2)$, hence, for $N = 60$, we obtain for this model
\beta
n_s \approx 0.97,
\end{equation}
independent of the value of $\beta$. The fact that the value of the scalar spectral index is independent of parameter $\beta$, or $p$ in the case of the power law potential, and we obtain the same value in both examples, should not come as a surprise. The reason becomes clear when considering the
dual $P(X,\phi)$ theory, which for any three-form potential, has the same quadratic scalar potential. Since we are
in a slow-roll regime, the functional form of the kinetic term is not important with respect
to the potential which leads to identical results for the spectral index.

It was found above that the choice of the exponential three-form potential offers an exit from inflation only if $\beta > 1/3$.
Taking $\chi_N^2 \approx 2/3$, this puts the lower bound, $c_s^2 \gtrsim 11/3$, which means that the speed of sound is superluminous.

The amplitude of the three-point function is controlled by $c_s^2$ and $\lambda/\Sigma$, which are given by
\beta
c_s^2 = 1 + 12 \beta \chi^2, \quad\quad\quad \frac{\lambda}{\Sigma} = -\frac{6 \beta^2 \chi^2 (1 + 4 \beta \chi^2)}{c_s^4}.
\end{equation}
Varying $\beta$ and substituting for $c_s^2$ and $\lambda/\Sigma$ $N$ $e$-folds before the end of inflation
into \eqref{eqfnl}, we obtain the dependence of $f_{\rm NL}^{equil}$ on $c_s^2$ which we show in Fig.~\ref{tejo2fnl} for $N = 60$.
\begin{figure}[tp]
\begin{center}$
\begin{array}{c}
\includegraphics[width=0.7\linewidth]{fnlExpo.png}
\end{array}$
\end{center}
\caption{
Dependence of $f_{\rm NL}^{equil}$ on ${c}_s^2$ for the exponential potential, $V = V_0 \exp(6\beta\chi^2)$ and $N = 60$. A small and generically positive non-Gaussian amplitude is found. Notice the lower bound $c_s^2 \gtrsim 11/3$ for this model.
\label{tejo2fnl}}
\end{figure}
In Fig.~\ref{tejo2rfnl} we see how $r$ and $f_{\rm NL}^{equil}$ relate to each other. In particular, we observe that the large values of $c_s^2$ render this model  disfavoured by current bounds on the ratio of tensor to scalar perturbations, $r \lesssim 0.2$ for $N = 60$.\footnote{Note that this bound significantly depends on the assumption of no running spectral index $n_s$. If a running $n_s$ is allowed the bound weakens to $r < 0.49$.} Of course, if we allow for more e-folds of inflation this reduces the relevant $\epsilon$ and hence $r$. A crude estimate for the minimal amount of inflation to bring this model into agreement with current $2\sigma$ bounds on $r$ may be obtained by assuming a speed of sound right at the lower bound of $c_s \gtrsim \sqrt{11/3}$ and $r \sim 0.2$, yielding $N \gtrsim 62$ for the exponential model considered here.
\begin{figure}[tp]
\begin{center}
\includegraphics[width=0.7\linewidth]{fnlvsrExpo.png}
\caption{\label{tejo2rfnl} The solid lines show how the parameters $f_{NL}^{\rm equil}$ and $r$ are related to each other for the exponential potential, $V = V_0 \exp(6 \beta \chi^2)$ and for $N = 50, 60, 80$ from top to bottom. The shaded region is disallowed by the WMAP 2$\sigma$ bound $r \lesssim 0.2$.\cite{Komatsu:2010fb}. End points for the solid lines at $f_{\rm NL}^{equil} \sim 0.215$ correspond to the lower bound $c_s^2 \gtrsim 11/3$ for this model.}
\end{center}
\end{figure}
Finally we see in Fig.~\ref{tejo2ampl} that the sign of $f_{\rm NL}^{equil} \sim 0.2$ is positive, albeit with a rather small amplitude. The non-Gaussian shape here is predominantly equilateral, but also picks up contributions in the enfolded/orthogonal limit.\footnote{Note that considering departure from slow-roll could modify this statement, e.g. potentially rendering the shape predominantly enfolded or orthogonal\cite{Noller:2011hd,Noller:2012ed}. For further details on fast-roll corrections also see\cite{Burrage:2011hd,Ribeiro:2012ar}.}
\begin{figure}[tp]
\begin{center}$
\begin{array}{c}
\includegraphics[width=0.7\linewidth]{shapeExpo.png} \\
\end{array}$
\end{center}
\caption{
The dimensionless bispectrum ${\cal A}(k_1,k_2,k_3)/k_1 k_2 k_3$ in the slow-roll limit for the exponential potential, $V = V_0 \exp(6\beta\chi^2)$ with $\beta = 1$, corresponding to $c_s \approx 3$ and for $N = 60$. Triangular shapes denote the equilateral, squeezed/local and enfolded limit clockwise from top right. $f_{\rm NL}^{equil} \sim 0.2$ and a predominantly equilateral shape is found. \label{tejo2ampl}}
\end{figure}

\section{Conclusions}
\label{sec:conclusions}

In this paper, we have considered three-form theories of inflation and their observational footprints. The salient results may be summarised as follows.
\begin{itemize}
\item We have calculated the action for a three-form theory perturbed about a FRW background, isolating the perturbed quadratic and cubic action expressed
purely in terms of the comoving curvature perturbation $\zeta$, as well as background and three-form quantities.

\item We then reviewed dualities between various $p$-form models in four dimensional space-times. In particular, these allow the
linking of the three-form theory under consideration to an effective scalar description with a $P(X,\phi)$ action, due to the fact that the most
general canonical and minimally coupled 3-form action in 4d has only one independent degree of freedom.
Since the three-form potential is mapped to the scalar field's kinetic term in a non-trivial way, such setups generically result in highly
non-canonical $P(X,\phi)$ theories, which do, however, all share a simple $\Phi^2$ potential.

\item Despite the complexity of the mapping, we were, however, able to use this duality
as a formal tool to give an alternative derivation
for the perturbed quadratic and cubic action. To do so we used existing ``non-Gaussian machinery'' for $P(X,\phi)$ theories, but
translated the end-result back to three-form quantities (essentially re-expressing $c_s, \Sigma$ and $\lambda$ in terms of the three-form variable $\chi$).
We showed that the result is identical to that obtained via the direct three-form computation.

\item Finally, using the tools developed, we  explored the observational and particularly non-Gaussian features of two example three-form models:
\begin{itemize}
\item First, we investigated
three-form inflation with a power-law potential. This simple
setup produces a constant speed of sound and generically
equilateral non-Gaussian shapes with negative $f_{\mathrm{NL}}^{equil}$, similar to
DBI-type inflation. Observational limits on this form of non-Gaussianity constrain this model.

\item Secondly, we explored the phenomenology for a three-form inflationary model with an exponential potential. We found
the non-Gaussianity was once again predominantly of equilateral shape, but in this case was unobservably small. We found, however,
that the model  can be constrained from limits on the tensor-to-scalar ratio $r$.
\end{itemize}
\item An interesting general result is that the spectral index for all three-form inflationary models is $n_s \approx 0.97$
to leading order in slow-roll when $60$ e-folds of inflation occur.  This was shown explicitly for the two models at hand, but will hold much more generally. The value of $n_s$ is therefore uniquely predicted once $N$, the number of e-folds of inflation, is specified, independent of the exact form of the three-form potential and as long as $\epsilon \ll 1$.
\end{itemize}

In conclusion, three-form inflation has a number of interesting features. Simple three-form potentials lead to
interesting models of inflation with potentially large non-Gaussian signatures, but with a spectral index which is
within the preferred WMAP bounds. It would be interesting to extend our studies of explicit three-form potentials to more complex
potentials which have no analytic form for the dual scalar theory; to use our framework to compute the trispectrum; to consider a multi three-form action.  We defer these questions to future studies.

\section*{Acknowledgments}
DJM acknowledges support from the Science and Technology Facilities Council
[grant number ST/J001546/1]. JN acknowledges support from the STFC.
NJN is supported by a Ci\^encia 2008 research contract funded by FCT and  through the projects PEst-OE/FIS/U
I2751/2011 and  grant CERN/FP/116398/2010. The authors would like to thank Antonio de Felice, Tomi Koivisto and Karim Malik for helpful discussions.


\section*{Methods}
\textbf{Setup.}
The master laser is a frequency stabilized Ti:Sapphire laser (MBR-110, Coherent) with frequency $\nu_0$ of 340.69656 THz. The frequency of the master laser is stabilized to a low-drift Fabry-P\'{e}rot Interferometer (FPI) with PDH technique. The FPI is placed in vacuum housing with temperature controlled to a stability of 5 mK. The linewidth of the master laser is below 50 kHz.

Part of the laser output is fed into a cavity enhanced frequency doubler to generate blue light. The blue light is controlled with an -200-MHz AOM1. The photon pairs are filtered with the optical gratings and the etalon, combined with SMF collection. The signal photons are collected with SMF and ready for use. Another part of the master laser output is used for AFC preparation. The AOM2 with center frequency of -250 MHz is used in double-pass configuration to down shift the laser frequency. The AOM3 and AOM4 have a center frequency of 200 MHz and are also used in double-pass configuration. Thus the center frequency of the AFC preparation matches with the photon source at 340.69646 THz.

The AFC preparation light and single-photon pulses are combined with 10:90 BS and collected with SMF. The light focus to a diameter of 100 $\mu m$ with a lens (f=250 mm). The sample is placed in a cryostat (Oxford Instruments, SpectromagPT) at a temperature of 1.5 K and with a superconducting magnetic field of 0.3 T in the horizontal direction. The two parallel Nd$^{3+}$:YVO$_4$ crystals' c-axes are placed in the horizontal direction. Because each crystal only strongly absorbs $H$-polarized light, the $H$-polarized components of input photons are stored in the first crystal, and the $V$-polarized components are stored in the second crystal after polarization rotation by the HWP. The detection efficiencies for both photons are approximately 35\% with 50-Hz dark counts.

The initial absorption depth of the sample is $d\simeq8.0$. With optimized spectral tailoring, we can prepare AFC with extremely small background absorption $d_0\leq0.05$. The Zeeman population relaxation lifetime is measured by spectral-hole burning, $T_Z=42.9$ ms in the current crystals, which is 7 times longer than that of 10 ppm doped crystals \cite{SHB}. With two-pulse photon echoes, we measured an optical coherence time $T_2=28.4$ us (homogeneous linewidth $\Gamma_h=11$ kHz). Because of the better coherence properties and the more efficient optical pumping achieved in the current sample, the AFC efficiency is significantly higher than that obtained in 10 ppm doped crystals \cite{polarization}.

The preparation and measurement timing is controlled by three arbitrary function generators (AFG, Tektronix, AFG3252). It takes 12.5 ms for AFC preparation. To avoid the fluorescence noise caused by the classical pump light, the measurement cycle begins after waiting for time $T_w=1.5$ ms after the preparation cycle is completed. $N$ trials of single photon pulses are stored in the sample in the 10-ms measurement phase. The complete preparation and measurement cycles are repeated at a frequency of 40 Hz. A schematic drawing of the timing sequence of the experiment is given in the Supplemental Information.

\textbf{Details of the photon source.}
The phase matching temperature to generate degenerate photon pairs is 44.50 $^{\circ}$C.
The power of the 440 nm pump light is 20 mW for all the measurements. With continuous pump, the coincidence rate of the unfiltered photon pairs is 1.4 MHz. The heralding efficiency is 0.20. The diffraction efficiency of the gratings (1200 lines per millimeter) is 0.88. The transmission efficiency of the etalon is 0.95. The etalon has a bandwidth of 700 MHz and a free spectral range of 50 GHz. The temperature drift of the etalon is smaller than 10 mK. After spectral filtering, the coincidence rate of the photon pairs is 6.1 kHz and the heralding efficiency is 0.10. The photon source is switched into pulsed operations for interface with the AFC. The temporal width of blue pump light is 32 ns. The transmission efficiency from the photon source to the final SPD for the signal photons is approximately 0.23.

A figure of merit for the nonclassical nature of the photon correlations is the normalized cross-correlation function $g^{(2)}_{s i}=p_{s i}/p_sp_i$, where $p_s$ ($p_i$) is the probability to detect a signal (idler) photon and $p_{s i}$ is the probability to detect a coincidence in a specified time window. In practice, $p_{s i}$ and $p_sp_i$ are determined by the number of coincidences in the time window centered on and away from the coincidence peak, respectively. Assuming that the second-order auto-correlations of signal and idler $g^{(2)}_{x}$ (where x=`s' for signal or `i' for idler) satisfying $1\leq g^{(2)}_{x} \leq2$, then the non-classicality is proved by measuring $g^{(2)}_{s i}$ greater than 2 \cite{entangle,ref1,g2}.
We note that, compared with the unfiltered source, $g^{(2)}_{s i}$ rises significantly after spectral filtering. We measured $g^{(2)}_{s i} \approx 452$ between the transmitted signal photon and the idler photon. The temporal width of the coincidence window is 2 ns and the detection window for noise is 400 ns away from the transmitted signal peak. The $g^{(2)}_{s i}$ drops after storage because of low recall efficiency, the temporal broadening and the rising contributions of dark counts. The temporal width of the coincidence window for retrieved signal is 10 ns. The detection window for related noise is placed at 400 ns away from the retrieved signal peak. An example of the measurement on the $g^{(2)}_{s i}$ is given in the Supplemental Information.

We are most grateful to A. J. Leggett for his comments on the preliminary version of this
manuscript and to C. Brukner and M. Paternostro for their feedback on the revised version. This
work was supported by the National Basic Research Program (2011CB921200), the CAS, National Natural
Science Foundation of China (Grant Nos. 11274289, 11325419, 61327901.), the Fundamental Research Funds for
the Central Universities (wk2030380004 and wk2470000011), and the EU STREP project PAPETS.