text:
\section{Introduction}
Ionized (H\,{\sc{ii}}) regions are known to trigger the formation of stars by means of various physical mechanisms (see Elmegreen~\cite{elm98} and Deharveng et al.~\cite{deh05} for a review). Several H\,{\sc{ii}} regions have been studied individually in the context of
triggered star formation, focusing on the associated neutral material and the young stellar population (Zavagno et al.~\cite{zav06}; Deharveng et al.~\cite{deh09}; Pomar\`es et al.~\cite{pom09}; Bieging et al.~\cite{bie09}). These studies have shown that the expansion of H\,{\sc{ii}} regions can trigger the formation of new stars of all masses. The $\it{Spitzer}$-GLIMPSE survey of the Galactic plane (Benjamin et al.~\cite{ben03}) detected nearly 600 bubbles (Churchwell et al. \cite{chu06}). Deharveng et al. (\cite{deh10}) selected a series of 102 ionized bubbles and studied the star formation in
their surroundings using $\it{Spitzer}$-GLIMPSE and MIPSGAL (Carey et al. \cite{car09}), radio (MAGPIS; Helfand et al. \cite{hel06} and VGPS; Stil et al. \cite{sti06}), and ATLASGAL (Schuller et al. \cite{sch09}) data.
They show that 86\% of these bubbles enclose H\,{\sc{ii}} regions, and that more than 20\% of 64 bubbles (for which the ATLASGAL angular resolution is sufficient to resolve the spatial distribution of cold dust) show massive star formation on their borders. This indicates that triggering is important in the creation of massive stars and that hot photodissociation regions (PDRs) are a good place to study the earliest phases of massive-star formation.

Our long-term aim is to use Hi-GAL (Molinari et al.~\cite{mol10}), combined with existing infrared, submillimeter, and radio surveys, to study the influence of H\,{\sc{ii}} regions on triggering the formation of stars in our Galaxy. Hi-GAL's extended wavelength coverage towards the far-infrared and its unprecedented sensitivity offer a unique opportunity to detect an embedded population of young sources that are not detected at shorter wavelengths. This allows us to observe intermediate and high-mass YSOs over the complete range of evolutionary stages. The unprecedented resolution of Hi-GAL also offers the opportunity to accurately characterize the physical nature of the sources by means of a detailed fit to their spectral energy distribution (SED).

In this Letter we study the bubble-shaped ionized region, N49, from the Churchwell et al.~(\cite{chu06}) catalogue to illustrate the purpose of our project. The N49 bubble was studied by Watson et al.~(\cite{wat08}) and by Deharveng et al.~(\cite{deh10}). However, these studies had no information in the 70--500\,$\mu$m range. This information, obtained with the PACS and SPIRE data presented here, allows us to discuss the star formation of this region in detail.

\section{Hi-GAL results on N49}
We use PACS and SPIRE images obtained in parallel mode at 70, 160, 250, 350 and 500\,$\mu$m (at resolutions of 5\arcsec, 11.4\arcsec, 17.9\arcsec, 25\arcsec and 35.7\arcsec, respectively) to determine the influence of H\,{\sc{ii}}\ regions on their surroundings using
the N49 H\,{\sc{ii}}\ bubble as an example. Our general goal is to quantify the efficiency of H\,{\sc{ii}}\ regions on triggering star formation around the perimeter of their radiation-driven expanding bubble.
Figure~\ref{N49} shows images of this region at different wavelengths, including the new PACS and SPIRE images.
ATLASGAL 870\,$\mu$m contours are superimposed on all images.
   \begin{figure*}[ht]
   \centering
   \includegraphics[angle=0,width=180mm]{14587Fig1.jpg}
   \caption{N49: a) Three-color composite image with 8\,$\mu$m {\it Spitzer}-GLIMPSE (blue), 24~$\mu$m MIPSGAL (green), 20~cm MAGPIS (red). ATLASGAL 870~$\mu$m contours are superimposed on all images (levels are 0.01, 0.1, 0.22, 0.31, 0.5, and 2~Jy/beam). ATLASGAL condensations 1 to 4 are identified. We have identified the two massive stage I YSOs (YSO\#1 and YSO\#3) discussed by Watson et al.~(\cite{wat08}). The ultracompact H\,{\sc{ii}}\ region is also identified; b) PACS 70\,$\mu$m image; c) PACS 160\,$\mu$m image; d) SPIRE 250\,$\mu$m image; e) SPIRE 350\,$\mu$m image; f) Two-color composite image with PACS 70\,$\mu$m (blue) and SPIRE 350\,$\mu$m image (red). The field is 5$\arcmin\,\times$ 5$\arcmin$, centered on $l$=28.83$\degr$, $b$=$-$0.23$\degr$
  }
              \label{N49}%
    \end{figure*}
This H\,{\sc{ii}}\ region is centered on $l$=28.83$\degr$, $b$=$-$0.23$\degr$ and located at 5.5~kpc (Anderson et al. \cite{and09}). It has a diameter of 4~pc. A compact H\,{\sc{ii}}\ region lies on its border (see Fig.~\ref{N49}a).  ATLASGAL 870\,$\mu$m emission shows a fragmented shell of collected material surrounding the ionized region (green contours on Fig.~\ref{N49}). Three condensations are seen surrounding the ionized region and a secondary peak is detected farther away. These are referred to as condensations 1 through 4, respectively, and are indicated in Fig.~\ref{N49}a.
Using ATLASGAL data, Deharveng et al.~(\cite{deh10}) derived a mass for the entire shell of 4250~M$_{\odot}$, and the mass of the brightest condensation in this shell (condensation 1 in Fig.~\ref{N49}a) of 2350~M$_{\odot}$. The NH$_3$ velocity measurements show that the condensations 1, 2, and 3 are associated with the ionized region (Wienen et al., $\it{in\,preparation}$). No measurement exists for condensation 4 but $^{13}$CO (1-0) data from the Galactic Ring survey (Jackson et al.~\cite{jac06}) show that it is also associated with N49. Two young stellar objects (identified YSO\#1 and YSO\#3 in Fig.~\ref{N49}a) are seen as bright 24\,$\mu$m sources towards condensation 1, the brightest at 870\,$\mu$m. These two YSOs are massive and in an early evolutionary stage (Watson et al. \cite{wat08}). Using high resolution 6.7~GHz methanol maser measurements, Cyganowski et al.~(\cite{cyg09}) suggest that a rotating disk is associated with YSO\#3.
An ultracompact H\,{\sc{ii}}\ region nearby YSO\#3 is shown in Fig.~\ref{N49}a as a bright, but more diffuse, source at 24\,$\mu$m. According to Deharveng et al.~(\cite{deh10}), this region may have been ionized by a B0V star.

Figure~\ref{N49}b shows the PACS 70\,$\mu$m image with  superimposed ATLASGAL 870\,$\mu$m contours. We see that YSO\#3 dominates the emission at longer wavelengths. The ultracompact H\,{\sc{ii}}\ region, bright at 24\,$\mu$m on the MIPSGAL image, is not bright in the 70-500\,$\mu$m range. A 70\,$\mu$m source is detected towards condensation 2 and 3. These sources are not detected at 24\,$\mu$m.

Figure~\ref{N49}d shows the SPIRE 250\,$\mu$m image with ATLASGAL 870\,$\mu$m contours superimposed. The sources previously detected at 70\,$\mu$m  towards condensations 2 and 3 are seen. Weak emission is observed towards condensation 4. This emission was not detected at 70\,$\mu$m and 160\,$\mu$m.

The PACS images clearly show the hot PDR that surrounds the H\,{\sc{ii}}\ region. The 70\,$\mu$m and the 8\,$\mu$m emission of the PDR are very similar. The only difference comes from the PACS 70\,$\mu$m source seen towards condensation 2 and observed as a dark spot on the 8\,$\mu$m image. This indicates that condensation 2 and the PACS IR source are probably located in front of the PDR. The situation is different for condensation 3, for which no clear absorption is seen at 8\,$\mu$m suggesting that the condensation is located behind the PDR. The two color-composite image shown in Fig.~\ref{N49}f, PACS 70\,$\mu$m (blue) and SPIRE 350\,$\mu$m (red), indicate that parts of the region where star formation is occurring are denser or colder or both, as indicated by the red emission.

The brightest IR source observed around N49, probably corresponding to YSO\#3, is located along an IR dark cloud seen on the 8\,$\mu$m GLIMPSE image. Since blobs are observed along this filament on larger scales, the origin of the formation of this bright IR source is questionable. Was it produced by the fragmentation and collapse of the shell surrounding N49 by means of the collect and collapse process (Elmegreen \& Lada~\cite{elm77}), or by the compression of a pre-existing blob by the N49 (and the ultracompact) H\,{\sc{ii}}\ region(s)? Morphological considerations favor the first option as a pre-existing clump would have distorted the ionization front, which is not what we observe here. We discuss the star formation observed around N49 in the framework of the collect and collapse model in Sect.~3.

\subsection{The spectral energy distribution fitting}
We used the on-line SED fitting tool of Robitaille et al.~(\cite{rob07}) to derive the physical properties of the sources observed towards ATLASGAL condensations bordering N49. Fluxes were measured at 70, 160, 250, 350, 500, and 870\,$\mu$m by rebinning all the images to the PACS 70\,$\mu$m resolution, using a square aperture of 100$\arcsec$, 50$\arcsec$, 70$\arcsec$, and 40$\arcsec$ for condensations 1, 2, 3, and 4, respectively (2.67, 1.33, 1.87, and 1~pc for the N49 distance of 5.5~kpc), and subtracting emission from a local background.

Watson et al.~(\cite{wat08}) identified seven YSOs towards N49. Their YSO\#5 and \#6, observed near condensation 4, are probably not related because they are not observed at longer wavelengths on the $\it{Herschel}$ images. Watson et al.~(\cite{wat08}) did not detect the red sources observed with PACS and SPIRE towards the ATLASGAL condensations 2 and 3. As seen in their Fig.~17, their SED fits are not constrained at wavelengths beyond 24\,$\mu$m due to a lack of data. The Hi-GAL coverage between 70 and 500\,$\mu$m, together with the ATLASGAL 870\,$\mu$m emission, strongly constrain the SED of the observed sources, allowing us to derive their physical properties.

In Table~1, we present our results of the SED fitting for all sources, taking into account fits with $\chi^2-\chi^2_{best}$ per datapoint $ < 3 $.
\begin{table}[h]
\caption{Results of the SED fitting}
\begin{tabular}{l r c c c}
  \hline
               &      &      &    &    \\
  & $M_{\rm star}$  &
   $\dot{M}_{\rm envelope}$ & $M_{\rm envelope}$ & $L$ \\
     & ($M_{\sun}$)  & (10$^{-3}$ $M_{\sun}~yr^{-1}$) & ($M_{\sun}$)&  (10$^3$ $L_{\sun}$) \\
  \hline
Cond.~1 &  20--30 & 8--10 & 5000  & 50\\
Cond.~2 & 8--10  & 4 & 8 & 1  \\
Cond.~3 & 8--15 & 2--5 & 2000--5000 & 2--5  \\
\hline
\end{tabular}
\end{table}

Towards condensation 1, a high luminosity, high mass source is found that has a high accretion rate. The values given in Table~1 are for the entire condensation 1, including YSO\#1 and YSO\#3. The 70\,$\mu$m emission is dominated by YSO\#3.  This is an extended green object and is associated with a methanol maser (Cyganowski et al.~\cite{cyg08}, \cite{cyg09}), which implies a massive YSO with an age younger than 3.5$\times$10$^{4}$ years (Breen et al.~\cite{bre10}). However, no emission is seen in its direction at 20~cm on the MAGPIS image. The high accretion rate derived for condensation~1 (see Table~1) may prevent the development of the ionized region (Churchwell~\cite{chu02}). The emission may also be optically thick at 20~cm.
The results of the fit for condensation~1 are shown in Fig.~\ref{fit}.
   \begin{figure}
   \centering
   \includegraphics[angle=0,width=90mm]{14587Fig2.jpg}
   \caption{Results of the fit to the fluxes measured towards condensation 1 in N49 using the on-line SED fitting tool of Robitaille et al.~(2007). The best-fit (dark solid line) is strongly constrained by the Hi-GAL wide wavelength coverage. The grey lines show all the range of models. The 24\,$\mu$m flux is taken from Watson et al.~(2008).}
              \label{fit}%
    \end{figure}

The sources observed towards condensations 2 and 3 are also massive objects with high accretion rates.
The emission that is barely detected at 250\,$\mu$m towards condensation 4 cannot be fit using the Robitaille et al. SED fitting tool because no data at shorter wavelengths exist to constrain the fit. This object is probably a cold core with no internal source of heating. However, it is interesting to detect this emission farther away from the ionization front as it may represent an earlier evolutionary stage in the star formation process. The 870\,$\mu$m ATLASGAL peak flux of 0.26~Jy/beam indicates a column density of N(H$_{2}$)=  6.8$\times$10$^{21}$~cm$^{-2}$ (2.2$\times$10$^{22}$~cm$^{-2}$) for a temperature of 20~K (10~K). This represents a low density condensation compared to the values derived by Deharveng et al.~(\cite{deh10}) around Galactic bubbles where star formation is clearly observed.

In all cases, the age of the source is not well constrained by the fitting procedure. Other indicators, such as chemical clocks, outflows, and masers, are needed to more tighly constrain the ages.
\section{Discussion}
Following the analytical approach of Whitworth et al.~(\cite{whi94}, see their Sect.~5), we discuss the observed star formation around N49 in the context of the collect and collapse process. The present radius of N49 (2~pc), its estimated age (0.5--1~Myr), and the spectral type of its ionizing star (O5V, Watson et al.~\cite{wat08}, Everett \& Churchwell~\cite{eve10}) allow us to derive a mean value for the initial density of the surrounding medium, n$_0 \simeq$~10$^4$~cm$^{-3}$.
Assuming this density and following the formula given by Whitworth et al.~(\cite{whi94}), we can estimate the time at which the fragmentation begins, $t_{\rm{fragment}} \simeq$0.5~Myr and the radius and column density of the shell at this time, $R_{\rm{fragment}} \simeq$1.55~pc and $N_{\rm{fragment}} \simeq$1.6$\times$10$^{22}$~cm$^{-2}$, which are similar to the values we derived towards condensation 4. All these values strongly depend on the sound speed in the shocked gas a$_s$ (which is between 0.2 and 0.6 km/s) and n$_0$, the initial density of the medium. The values we give here are derived using a$_s$=0.2 km/s and n$_0$=10$^4$ cm$^{-3}$. If n$_0$ is lower, we obtain a $t_{\rm{fragment}}$ higher than the estimated age for N49, which is not possible because the fragments are already seen. The same happens if we increase a$_s$. However, the first set of values that provide acceptable results for the start time of the fragmentation and radius (and column density) of the shell at that time does not provide acceptable results for the mass of the fragments, which are too low (around 7.5~M$_{\odot}$) compared to the masses derived by Deharveng et al.~(\cite{deh10}) using the ATLASGAL measurements (of about 2500~M$_{\odot}$, 240~M$_{\odot}$, and 355~M$_{\odot}$ for condensations 1, 2 and 3, respectively). As discussed below, the wind dynamics that exists towards N49 can modify this simplified approach. A finer age determination for the sources observed towards the ATLASGAL condensations may help in more tightly constraining this model.

The Hi-GAL images and results of the SED fits show that five massive YSOs (YSO\#1 and \#3, the UC H\,{\sc{ii}}\ region, and IR sources towards condensations 2 and 3) are formed (or are in the process of forming) on the border of N49. This region is ionized by an O5V star and stellar winds are believed to play a crucial role in the formation of this bubble. Everett \& Churchwell~(\cite{eve10}) studied in detail the case of N49 as a wind-blown bubble. They conclude that the survival of dust within the ionized gas in the presence of winds requires a high density in the surrounding medium. This high density together with greater dynamical motion caused by stellar winds might have increased the efficiency of triggered star formation, in terms of number and mass of the second-generation stars. Other dusty wind-blown bubbles will be observed in Hi-GAL and the unbiased access to the YSOs population will help us to understand whether winds play a crucial role in increasing the efficiency of triggered star formation.

\section{Conclusions}
We have presented the first PACS and SPIRE images from Hi-GAL of the bubble-shaped Galactic H\,{\sc{ii}}\ region N49. This region is used as an illustration of the study dedicated to the star formation triggered by Galactic H\,{\sc{ii}}\ regions that we plan to lead for the whole survey, combining the Hi-GAL results with other infrared and radio surveys of the Galactic plane. We have shown that:
   \begin{itemize}
      \item{The Hi-GAL SPIRE and PACS images allow us to study the distribution of young sources towards N49. The far-IR fluxes have been measured and strongly constrain the spectral energy distribution of these sources. This allows us to characterize their properties. }
     \item{The PACS images reveal the existence of red young stellar objects towards two ATLASGAL condensations, sources that had not been previously detected at shorter wavelengths. The ultracompact H\,{\sc{ii}}\ region is not clearly seen in the $\it{Herschel}$ range. The bright YSO\#3 Watson et al. observed to be coincident with condensation 1 dominates the emission at longer wavelengths.}
      \item{SED fits for the 3 sources detected by $\it{Herschel}$ towards millimeter condensations using the Robitaille et al.~(\cite{rob07}) model show that these sources are young and massive. However, their age has not been constrained and other indicators are needed to refine the discussion of star formation history in this region.}
	\item{Five massive stars are forming in the N49 PDR. The high star formation efficiency in N49 may be due to the presence of winds from the first generation massive star.}
   \end{itemize}
The study of N49 using a multiwavelength approach shows that Hi-GAL enables measurements for the crucial far-IR range to be made that are essential to infer the properties of the YSOs. Hi-GAL clearly provides important insight into star formation triggered by expanding H\,{\sc{ii}}\ regions. Seventy-six H\,{\sc{ii}}\ regions of all shapes are detected in the $l$=30$\degr$ field and a higher density of Hi-GAL sources is clearly observed towards these regions, indicating that star formation triggered by H\,{\sc{ii}}\ regions may be an important process.

\begin{acknowledgements}
    Part of this work was supported by the ANR (\emph{Agence Nationale pour la Recherche}) project
      "PROBeS", number ANR-08-BLAN-0241 (LA).
\end{acknowledgements}


text:
\section{Introduction}

The apparent detection of B-mode polarization of the Cosmic Microwave Background (CMB) on $l \sim 100$ by the BICEP2 observation has stirred significant interests in primordial gravitational waves  \cite{Ade:2014xna}. If confirmed to be primordial \cite{Flauger:2014qra, Mortonson:2014bja},  this detection implies the existence of primordial gravitational waves with the scalar-to-tensor ratio  $r= 0.2^{\, +0.07}_{\, -0.05}$.
In the near future, the BICEP2 detection of $r$ will be cross-checked by many on going and forthcoming experiments. For example, Planck \cite{Planck16}, SPTPol \cite{Austermann12}, ACTPol \cite{Niemack10}, PolarBear \cite{PolarBear} and CLASS \cite{Eimer12}.

These forthcoming experiments can  probe more detailed properties for the tensor mode.  For example, when the statistical features of the primordial perturbations are not isotropic, a large number of new observables arises, including the correlation functions with different multipole moment $l$, and TB and EB cross-correlations, which have been forbidden by the isotropic statistics of the primordial fluctuations \cite{Watanabe:2009ct, Chen:2013eaa, new paper}.

Motivated by the future data, in this paper we study inflationary dynamics producing statistical
 anisotropies in details. We investigate the anisotropic inflation scenario with a charged scalar inflaton coupled to the gauge field. This model is studied originally in \cite{Emami:2010rm} at the background level and its predictions for statistical anisotropies in curvature perturbation power spectrum were studied in \cite{Emami:2013bk}. However,  the perturbations in tensor sectors and their correlation with the scalar perturbations  have  not been explored so far.  We show that the gauge coupling induces large statistical anisotropies in tensor perturbations as compared to model of anisotropic inflation with no charge coupling \cite{Ohashi:2013qba}. We show that while statistical anisotropies in temperature power spectrum can be small as required by the Planck data \cite{Kim:2013gka}, the tensor mode can develop significant statistical anisotropies.  Therefore it is important for the forthcoming experiments to look for the statistical anisotropies in the B-mode polarization even though the statistical anisotropies in temperature map is well-constrained. For this purpose, in this work we calculate the primordial correlations of the curvature perturbation $\zeta$, and the two tensor modes $h_+$ and $h_\times$. The CMB temperature and polarization correlations TT, TE, EE, TB, EB and BB are then calculated from these primordial perturbations, for the same and different multipole moments.

A novel result in our analysis is that  the anisotropies in the tensor sector also contribute to the TT correlation on the CMB. However, the transfer function from primordial tensor to CMB temperature  decays towards large $l$. Thus the TT anisotropy coming from the primordial tensor has a decaying amplitude and is highly suppressed after $l \gtrsim 100$. As a result, we naturally obtain anisotropies at low multipoles of TT without modifying the high multipoles. This scale-dependent
anisotropies will have a better fit to the CMB anomalies.

Another motivation of the present work is that, the Planck experiment puts a limit of tensor-to-scalar ratio $r<0.12$. This result is in tension with the BICEP2 detection. The tension is reported to be as unlikely as $0.1\%$. Recently, it is proposed \cite{Contaldi:2014zua} that with anti-correlated curvature and tensor perturbations, the tension between Planck and BICEP2 may be reconciled. However, our detailed study show that such a mechanism does not work. As we shall show, there is indeed an anti-correlation between the curvature and tensor perturbations. However the contribution to the TT power spectrum is cancelled when summing over the angular mode $m$.

The rest of the paper is organized as follows. In Section \ref{model} we present our model of anisotropic inflation and review the background dynamics and the perturbations. In Section
\ref{correlations} we calculate the anisotropic scalar and the tensor power spectra and the scalar-tensor
cross-correlations. In Section \ref{numerical} we present the predictions of our model for CMB.
We relegate the technical details into appendices.

\section{Anisotropic Inflation}
\label{model}

Here we review anisotropic inflation.   For earlier works on various aspects of anisotropic inflation see
\cite{Watanabe:2009ct, Ohashi:2013pca, Ohashi:2013qba, Kim:2013gka,
Ohashi:2013mka, Emami:2010rm, Thorsrud:2012mu,
Dulaney:2010sq, Gumrukcuoglu:2010yc, Watanabe:2010fh, Bartolo:2012sd,
Funakoshi:2012ym, Yamamoto:2012sq, Emami:2013bk,
Shiraishi:2013vja, Abolhasani:2013zya,
Abolhasani:2013bpa, Ramazanov:2013wea, Nurmi:2013gpa, Urban:2013spa, Fujita:2013qxa, Shiraishi:2013oqa, Jain:2012vm, Kanno:2010nr, Murata:2011wv, Bhowmick:2011em, Hervik:2011xm,  Boehmer:2007ut, Koivisto:2008xf, Yokoyama:2008xw, Baghram:2013lxa, Watanabe:2010bu, Yamamoto:2012tq, Emami:2011yi, Emami:2014tpa, BeltranAlmeida:2011db, Rodriguez:2013cj, Lyth:2013sha, TuanQ}, for a review of anisotropic inflation see \cite{Soda:2012zm, Maleknejad:2012fw}. Also see \cite{Pereira:2007yy, Pitrou:2008gk, Pitrou:2012ge} for related works on primordial anisotropies.

As mentioned before, the model contains  a $U(1)$ gauge field  $A_\mu$, as in Maxwell theory, which is turned on at the background level.  However, it is well-known that the Maxwell theory suffers from the conformal invariance on expanding backgrounds.
Therefore, the background gauge field $A_\mu(t)$ is diluted exponentially during inflation. Furthermore,  the quantum excitations of the gauge field $\delta A_\mu (t, \bfx)$ are not scale-invariant. One interesting mechanism to break the conformal invariance is to couple the gauge field to the inflaton field non-trivially.

With these general discussions in mind we present our model of anisotropic inflation.  The model contains  a complex inflaton field $\phi$ which is charged under the $U(1)$ gauge field $A_\mu$ with the electric charge (coupling) $\e$. The action is given by
\ba
 \label{action} S= \int
d^4 x  \sqrt{-g} \left [ \frac{M_P^2}{2} R - \frac{1}{2} D_\mu \varphi
\,  \overline {D^\mu \varphi} -   \frac{f^{2}(\varphi)}{4} F_{\mu \nu} F^{\mu
\nu}  - V(\varphi, \overline \varphi) \right] \, ,
\ea
where $M_P$ is the reduced Planck mass.

As usual, $F_{\mu \nu}$ is  the gauge field strength given by
\ba F_{\mu \nu} = \nabla_\mu A_\nu
- \nabla_\nu A_\mu  = \partial_\mu A_\nu - \partial_\nu A_\mu \, .
\ea
In addition,  the covariant derivative associated with the $U(1)$ gauge field is given by
\ba
D_\mu \varphi = \partial_\mu  \varphi + i \e \,  \varphi  \, A_\mu \, .
\ea
where $\e$ represents  the electric gauge coupling. The model of anisotropic inflation based on the above action was studied at the background level
 in \cite{Emami:2010rm} and its  analysis for curvature perturbations were performed in
 \cite{Emami:2013bk}.

Note that in Maxwell theory $f(\varphi)=1$. However, as explained above,  we need a time-dependent gauge kinetic coupling $f(\varphi)$ in order to break the conformal invariance such that the background gauge field survives the exponential expansion and the gauge field excitations acquire a nearly scale-invariant power spectrum.
As we shall verify below one requires $f \propto a^{-2}$  in order to obtain a scale-invariant power spectrum for the  gauge field quantum excitations.  This choice of the gauge kinetic coupling corresponds to a constant background electric field energy density during inflation.

The complex scalar field can be decomposed into a radial part and an axial part via $\varphi(x) = \phi(x) \,  e^{i \theta(x)}$. As usual, we assume the model is axially symmetric in field space so the potential
$V$ and $f(\varphi)$ are only functions of $\varphi \bar \varphi=  \phi ^2$. It is convenient to go to the unitary gauge in which $\theta=0$. In this gauge, $\varphi$ becomes real-valued and in the analysis below we  take $\varphi=\varphi^* = \phi$.

We consider the coordinate system in which the background gauge field is turned on along the $x$-direction so  $A_\mu = (0, A_x(t), 0,0)$. With the background gauge field in the $x$-direction the
background space-time becomes anisotropic, taking the form of  type I Bianchi Universe, with the metric
\ba
\label{bian0}
\label{metric} ds^2 = - dt^2 + e^{2\alpha(t)}\left( e^{-4\sigma(t)}d x^2
+e^{2\sigma(t)}(d y^2 +d z^2) \right) \, .
\ea
Note that the metric (\ref{bian0}) still has a residue two-dimensional rotational symmetry on the $y-z$ plane. In this convention $\dot \alpha$ measures the averaged Hubble expansion while $\dot \sigma(t)$ measures the level of anisotropic expansion.  However, on the observational grounds, as we shall see in next Section, the level of anisotropy in curvature perturbations is not more than few percent which subsequently is translated  into the conclusion that
 $\dot \sigma /\dot \alpha \ll 1$. As a result, we can treat the analysis prturbatively in term of background  anisotropy.

The background fields equations are  \cite{Emami:2010rm}
\ba
\label{back-A-eq}
\partial_t{\left(  f^2(\phi) e^{\alpha + 4 \sigma} \dot A_x        \right)}& =& - \e^2 \phi^2 e^{\alpha + 4 \sigma}  A_x \\
\label{back-rho-eq}
\ddot\phi+3\dot \alpha\dot \phi+ V_{,\phi}+ \left(
-f(\phi)f_{,\phi}(\phi)\dot A_x^2 +\e^2 \phi A_x^2   \right) e^{-2\alpha+4\sigma}&=&0  \\
\label{Ein1-eq}
\frac{1}{2}\dot
\phi^2+V(\phi)+ \left(   \frac{1}{2}f^2(\phi)\dot
A_x^2 +\frac{\e^2\phi^2}{2}A_x^2 \right) e^{-2\alpha+4\sigma}
&=&
3 M_P^2 \left(   \dot \alpha^2-\dot \sigma^2 \right)  \\
\label{Ein2-eq}
V(\phi)+  \left(  \frac{1}{6}f^2(\phi)\dot
A_x^2+\frac{\e^2\phi^2}{3}A_x^2  \right)e^{-2\alpha+4\sigma}
&=& M_P^2 \left( \ddot \alpha    + 3 \dot \alpha^2 \right)  \\
\label{anisotropy-eq}
\left(\frac{1}{3}f^2(\phi)\dot A_x^2  -\frac{\e^2\phi^2}{3}A_x^2    \right) e^{-2\alpha+4\sigma}
&=& M_P^2\left( 3\dot \alpha \dot \sigma+ \ddot \sigma      \right)\, ,
\ea
in which a dot indicates derivative with respect to $t$ and $f_{, \phi}= \partial_\phi f$ and so on.

Here we summarize the main features of the above equations. Eq. (\ref{back-A-eq}) is the Maxwell equation in the inflationary background with the gauge coupling $\e$ appearing in the source term.
Eq. (\ref{back-rho-eq}) is the modified Klein-Gordon equation in the presence of the gauge field.
As we see the scalar field dynamics are affected by the gauge field via the last two terms in the bracket. The first term in the bracket in Eq. (\ref{back-rho-eq}) comes from the gauge kinetic term while the last term  comes from the charge coupling $\e$.
As we shall see, these two terms play crucial roles both at the background and the perturbation levels. The remaining equations are the Einstein equations with Eq. (\ref{anisotropy-eq}) controlling the dynamics of the anisotropy $\dot \sigma$.

In general it is not easy to solve the above set of equations analytically, even in the slow-roll limit. However, as we mentioned before,  the level of anisotropy in curvature perturbations is small
which also yields $\dot \sigma / \dot{\alpha} \ll 1$.  This means that the background expansion is nearly
isotropic.  Therefore,  as in in conventional models of inflation, the background expansion is controlled  by
the potential term $V$. However, we expect the gauge field also to contribute in the background expansion in the form of electric field energy density $E_x^2$ where $E_x = F_{0 x}$.
In order for the background to be nearly isotropic  we require that the electric field energy density to be very small compared to $V$. It is convenient to parameterize the gauge field (electric) energy density
by the parameter $\Omega_A$ via
\ba
\label{Omega-def}
\Omega_A \equiv \frac{\dot A_x^2 f(\phi)^2 e^{-2 \alpha}}{2 V} \, .
\ea
In order for the anisotropy to be small we assume $\Omega_A \ll 1$.

Alternatively, for the perturbation analysis, it is more convenient to express the background metric (\ref{bian0}) in the following form
\ba
\label{Bianchi-metric}
ds^2 = a(\eta)^2 (d \eta^2 + d x^2) + b(\eta)^2 ( d y^2 + d z^2)
\ea
where $a= e^{\alpha - 2 \sigma}$ and $b=e^{\alpha + \sigma}$ and the conformal time  $\eta$ is defined via via $dt = a(\eta) d \eta$.

We will work in the slow roll limit where the change in expansion rates in all directions are small. Defining the average Hubble expansion rate via $\dot \alpha \equiv H$, the
slow-roll parameters are given by
\ba
\epsilon_H \equiv - \frac{\dot H}{H^2}  \quad , \quad
\eta_H \equiv \epsilon_H - \frac{\ddot H}{2 H \dot H}

\ea
We work in  the slow-roll limit where $\epsilon_H , \eta_H \ll 1$ and  to leading order in slow-roll parameters and anisotropy $a \simeq b \simeq -1/H \eta$.

Before we study the attractor solution in next subsection, let us discuss the effects of the gauge coupling $\e$. The main effect of the gauge coupling is captured by the third term in bracket in Eq. (\ref{back-rho-eq}). From this term we see that the interaction $\e^2 \phi^2 A_\mu A^\mu$ induces a
time-dependent mass for the inflaton. Since this induced mass is exponentially time-dependent, we expect it becomes important only towards the end of inflation where the exponential growth
of the gauge field becomes significant.  As studied in details in
\cite{Emami:2010rm} inflation ends when the induced mass from the  back-reaction $\e^2 \phi^2 A_\mu A^\mu$ becomes comparable to the bare inflaton mass $m$. Therefore, in order for inflation to sustain  long enough, the back-reaction $\e^2 \phi^2 A_\mu A^\mu$ is negligible during much of the period of inflation and it only controls the mechanism of end of inflation.  As studied  in \cite{ Emami:2010rm} the end of inflation depends logarithmically on $\e$ where $
\alpha_e \sim -\frac{\ln \e}{2} + ...
$
where dots indicate the dependence on other parameters such as the inflaton value and its mass.
Therefore, the larger is the gauge coupling $\e$, the shorter is the period of inflation.
As explained above this is easily understood from the fact that the induced mass  for the inflaton field from  the Higgs mechanism, $\e^2 A_\mu A^\mu $, becomes comparable to the bare inflaton mass.

In the remaining analysis of the background dynamics we neglect the effects of $\e$ during
inflation. In particular, one can neglect the source term in Maxwell equation  (\ref{back-A-eq})
during much of period of inflation.   In this approximation one can easily solve the Maxwell equation (\ref{back-A-eq}) to get
\ba
\label{gaugefield}
\dot{A_{x}}= f(\phi)^{-2}e^{-\alpha(t)-4\sigma(t)}p_{A} \, ,
\ea
where  $p_{A}$ is a constant of integration.  Of course, this approximation breaks down near the end of inflation where the induced mass term from the interaction $\e^2 A_\mu A^\mu \phi^2$ becomes important.

\subsection{The Attractor Solution}
\label{background formulas}

As we discussed above  the anisotropy is small and  the average Hubble expansion rate in Eq. (\ref{Ein1-eq}) mainly comes from the isotropic potential term. However, the back-reactions of the gauge field on the inflaton field induce an effective mass term for the inflaton  as captured by the last term in Eq. (\ref{back-rho-eq}). Therefore, the dynamics of the  inflaton field is affected by the  back-reactions of the gauge field.  A key  observation was made in  \cite{Watanabe:2009ct} where it is shown that
for a general form of slow-roll potential $V$ and  with an appropriate choice of $f(\phi)$ the system reaches an attractor regime where $ \Omega_A \sim \eH$.  As a result,  the modified Klein-Gordon equation still admits a slow-roll solution but now with  a modified effective mass for the inflaton which is  induced from the back-reactions of the gauge field.

Let us see under what condition $\Omega_A$ reaches a near constant value during the attractor regime. Combining Eqs. (\ref{gaugefield})  and (\ref{Omega-def}) we obtain
\ba
\label{R-def2}
\Omega_A = \frac{p_A^2}{2 V} f(\phi)^{-2} e^{-4 \alpha - 8 \sigma} \, .
\ea
In order for $\Omega_A$ to reach a nearly constant value, and neglecting the contribution of $\sigma$
in the small anisotropy limit,   one has to choose $f(\phi)$ such that $f(\phi) \propto e^{-2 \alpha(t)} = a(t) ^{-2}$. Now to find $f$ as a function of $\phi$, our job is to solve the background expansion equation and  express $a(t)$ in terms of $\phi$. In the small anisotropy limit, and for a given potential $V(\phi)$,  the background isotropic expansion is given by
\ba
\label{a-scale}
a(t) \propto \exp \left[ - \int d \phi \frac{V}{ V_{,\phi}} \right] \, .
\ea
This indicates that  if one chooses
\ba
\label{f-scale}
f \propto \exp \left[ -n  \int d \phi \frac{V}{V_\phi} \right] \, ,
\ea
then one obtains $f \propto a^{n}$.

The exact form of $f(\phi)$  depends on $V(\phi)$.  To be specific, in this work we consider the  simple  chaotic potential
\ba
\label{Chaoticpotential}
V=  \frac{1}{2} m^2 | \varphi|^2 =
\frac{1}{2} m^2 \phi^2 \, .
\ea
Plugging this form of potential into Eq. (\ref{f-scale}) yields   \cite{Watanabe:2009ct}
\ba
\label{f-form0}
f(\phi) = \exp {\left( \frac{c\phi^2}{2 M_P^2}  \right)} \, ,
\ea
where $c$ is a constant.

Now we examine under what conditions the system allows for an attractor solution where the anisotropy reaches a sub-dominant but constant value.  In the small anisotropy limit, the background equation is given by the potential term  and we have
\ba
\label{Friedmann}
3 M_P^2 \dot \alpha^2 \simeq V = \frac{m^2 \phi^2}{2} \, .
\ea
On the other hand, in the slow-roll regime, the scalar field equation, Eq. (\ref{back-rho-eq}), is given by
\ba
\label{KG-eq1}
3 \dot \alpha \dot \phi = -m^2 \phi + \frac{c\,  p_A^2 \,  \phi }{M_P^2} f(\phi)^{-2} e^{-4\alpha} \, .
\ea
Note that in order to get the second term in the right hand side  the solution for $\dot A_x$ and
 the form of $f(\phi)$, given respectively in  Eq. (\ref{gaugefield}) and  Eq. (\ref{f-form0}), have been used.  Eliminating $\dot \alpha$ from  Eqs. (\ref{Friedmann})  and (\ref{KG-eq1}) we obtain
an equation for $\phi$ in terms of the number of e-folds $\alpha$ as
\ba
\phi \frac{d \phi}{d \alpha} = - 2 M_P^2 + \frac{2 c p_A^2}{m^2} e^{- c \phi^2/M_P^2} e^{- 4 \alpha} \, .
\ea
One can easily solve  this differential equation to get
\ba
\label{sol-general}
e^{-4 \alpha} e^{-c \phi^2/M_P^2} = \frac{m^2 (c-1) M_P^2}{c^2 p_A^2} \left[ 1+ C\,  e^{-4 (c-1) \alpha} \right]^{-1} \, ,
\ea
where $C $ is a constant of integration. Now the important point to note  is that for
$c\ge 1$, the second term in the square bracket above decays exponentially  during inflation
and we quickly reaches the attractor regime
\ba
\label{att1}
e^{-4 \alpha} e^{-c \phi^2/M_P^2} \simeq \frac{m^2 (c-1) M_P^2}{c^2 p_A^2} \, .
\ea
This indicates that $\Omega_A$,  defined in Eq. (\ref{R-def2}), reaches a constant value during the attractor regime as claimed.  More specifically, from Eq. (\ref{R-def2}) we obtain
\ba
\Omega_A = \frac{m^2 M_P^2}{2 V} \frac{c-1}{c} \, .
\ea
On the other hand, using  the remaining Einstein equation,  one can also verify that the slow-roll parameter is
\ba
\epsilon_H \equiv -\frac{\dot H}{H^2} = -\frac{\ddot \alpha}{\dot \alpha^2} = \frac{2 M_P^2}{c \phi^2} \, .
\ea
As a result, we conclude
\ba
\label{R-app}
\Omega_A = \frac{c-1}{2c}\epsilon_{H} = \frac{1}{2}I\epsilon_{H} \, ,
\ea
where we have introduced the anisotropy parameter $I$ via
\ba
 I \equiv\frac{c-1}{c} \, .
 \ea
 Note that in order for the gauge field to survive the exponential expansion and the system reaches the attractor regime we require $I \ge 0$.

It is also instructive to look at the scalar field equation in the attractor regime.  Plugging  Eq. (\ref{att1}) back into the scalar field equation  (\ref{KG-eq1})  we get
\ba
\label{klinGordon}
M_P^{-2}\frac{d \phi}{d \alpha} \simeq -\frac{V_\phi}{ V} +
\frac{c-1}{c}\frac{V_\phi}{ V}  = \frac{-2 c}{\phi}
\ea
This equation implies  that the back-reaction of the gauge field has reduced the effective mass of the inflaton such that    $m \rightarrow m - \frac{c-1}{c}m = \frac{m}{c}$.

The above equation can be easily solved to find $\phi$ as a function of the number of e-foldings $\alpha$ as
\ba
\label{klinN}
\phi_{e}^2 - \phi^2 = 4 M_P^2 \alpha (1-I) \, ,
\ea
where $\phi_e$ represents  the value of $\phi$ at the end of inflation. Note that the above solution is valid until end of inflation when the effects of the gauge coupling $\e$ starts to dominate. As we discussed before, there is additional induced mass from the term $\e^2 A_\mu A^\mu \phi^2$ which
terminate inflation quickly when it becomes comparable to the bare mass $m$.

\subsection{Perturbations}

Here we present cosmological perturbations in anisotropic inflation. The perturbation analysis
for various models of inflation were studied in
\cite{Dulaney:2010sq, Gumrukcuoglu:2010yc, Watanabe:2010fh, Bartolo:2012sd,
Funakoshi:2012ym, Yamamoto:2012sq, Emami:2013bk,
Shiraishi:2013vja, Abolhasani:2013zya,
Abolhasani:2013bpa, Ohashi:2013qba}.
The general form of the metric and matter perturbations have been studied in \cite{Emami:2013bk}, see also \cite{Watanabe:2010fh}.  For the perturbations in  the metric we note that there are some non-dynamical degrees of freedom, namely $\delta g_{0 \mu}$, which should be integrated out in order to calculate the dynamical action. Integrating out these non-dynamical metric degrees of freedom  lead us to extra terms in the action. Therefore, one important question is what the leading contributions in the final action are. Are they from the metric sector perturbations or from the matter sector perturbations?
Fortunately, as it has been verified  in \cite{Emami:2013bk}, it turns out that the metric sector perturbations are  either slow-roll suppressed or would cancel with each other. As a result,  we conclude that  the leading terms in the interactions come from the matter sector.  In other word,  we do not need to consider the perturbations from the non-dynamical degrees of freedom  in the metric. To simplify the situation further, we go to flat gauge where the curvature perturbations is given by the
inflaton perturbations $\zeta = - \frac{H}{\dot \phi} \delta \phi$. As a result, the metric has no scalar perturbations and we are left with the simple form of metric perturbations
\ba
\label{dynamical metric}
ds^2 = a(\eta)^2 \left(-d\eta^2 + \left[\delta_{ij}+ h_{ij} \right]dx^i dx^j \right ) \, .
\ea
Note that since we work in small anisotropy limit, we can set $a=b$ to leading order.
The corrections in our results below from using this assumption
will be suppressed by additional factors of $I$ or $\epsilon_H$
which are small.  The perturbations $h_{ij} $ represents the tensor modes subject to the transverse and traceless conditions  $\partial_i h_{ij} = 0$ and $ h_{ii}=0$ where the repeated indices are summed.  We denote the two independent polarizations of the metric by $h_\times$ and  $h_+$.

As for the perturbations in gauge field  sector  there is one non-dynamical degree of freedom,  $\delta A_{0}$, which must be integrated out from the action. However, similar to the case of   non-dynamical degrees of freedom from the metric perturbations, it turns out that the new terms from integrating out $\delta A_{0}$ are also sub-leading. As a result, the leading interaction terms in the total Lagrangian come from  the dynamical degrees of freedom $\delta A_i$.

In order to simplify the analysis further, we can use the remaining  two-dimensional rotational symmetry on the $y-z$ plane to set $k_{z} \equiv 0$ so $\overrightarrow{k} = (k_x , k_{y} , 0) =
k \, (\cos \theta, \sin \theta, 0 )$ in Fourier space. In addition, since the gauge field  has three polarizations in the unitary gauge,
two transverse and one longitudinal polarizations, we  can choose the following ansatz for the gauge field perturbations, i.e. $\delta A_i$,
\ba
\delta A_\mu^{(S)} = (\delta A_0, \delta A_1, \partial_y M, 0) \quad \quad , \quad \quad
\delta A_\mu^{(V)} =(0, 0, 0, D) \, .
\ea
Where we have defined, $\delta A_1, \delta A_2$ and $\delta A_3$ to be $\delta A_x, \partial_y M$ and $D$ respectively. Here $A_\mu^{(V)}$ refers to one transverse polarizations in the vector sector  while $A_\mu^{(S)}$
represents the two polarizations in the scalar sector.
Furthermore, we can decompose the two polarizations in $A_\mu^{(S)}$ into  one transverse and one longitudinal polarizations
as follows \cite{Emami:2013bk}
\begin{align}
\label{D12}
D_{1}&\equiv \delta A_{1} -ik \cos{\theta}M \\
D_{2}&\equiv \cos{\theta} \delta A_{1} + ik \sin^2{\theta}M \, .
\end{align}
In this decomposition $D_{1}$ represents  the transverse polarization while $D_{2}$ refers to the longitudinal polarization of the gauge field.  However, as it has been demonstrated in \cite{Emami:2013bk}, the interactions containing the longitudinal mode are  exponentially suppressed during inflation and can be neglected from the analysis. Physically, this is understandable since the interactions containing the
longitudinal mode $D_2$ originate from the ``Higgs mechanism'' via the interaction $\e^2 A_\mu A^\mu \phi^2$ which are exponentially suppressed during much of the period of inflation as discussed before.

We can quantize the curvature perturbation and the gauge field perturbations as usual.
For the curvature perturbation, note that we work in the flat gauge so
\ba
\label{zeta-phi}
\zeta = - \frac{H}{\dot \phi} \delta \phi  = \frac{\delta \phi}{M_P \sqrt{2 \epsilon_H}}.
\ea
Expanding the quantum operator $\widehat \zeta$ in terms of the annihilation and the creation operator
$a(\bfk)$ and $a^{\dagger}(\bfk)$ we have
\ba
\widehat \zeta (\bfx , \eta) =  \int \frac{d^3k}{(2\pi)^{3/2}} e^{i \mathbf{k}.\mathbf{x}}
\widehat\zeta(\mathbf{k},\eta) \quad , \quad
\widehat\zeta(\mathbf{k},\eta) = \zeta(k, \eta) a(\bfk) +   \zeta^*(k, \eta)  a^{\dagger}(- \bfk)
\ea
where the creation and the annihilation operators satisfy the usual commutation relation
$ [ a(\bfk),    a^\dagger  (\bfk')] = \delta^{(3)} (\bfk-\bfk' )$.

The wave function  of the curvature perturbation has the standard form of  the excitations of
a  massless scalar field on a dS background
\ba
\label{zeta-wave}
\zeta_k(\eta) = \frac{i H \eta }{M_{P}\sqrt{2\epsilon_{H}k} }  \left( 1 - \frac{i}{k\eta} \right) e^{-ik\eta} \, .
\ea
The power spectrum of the curvature perturbations is given by
\ba
\langle   \widehat{\zeta}(\mathbf{k_1})  \widehat{\zeta} (\mathbf{k_2}) \rangle
= (2 \pi)^3 \delta^{(3)} (\bfk_1 + \bfk_2) P_\zeta (k_1) \quad , \quad
\calP_\zeta \equiv   \frac{ k_1^3  }{ 2 \pi^2}  P_\zeta(k_1)
\ea
In particular, the power spectrum for the free isotropic theory is
\ba
\label{calP0}
\calP_\zeta^{(0)} = \frac{H^2}{8 \pi^2 \epsilon_H M_P^2}\, .
\ea

Similarly, the quantum excitations of the gauge field perturbations $\widehat D_{1\bf k}(\eta)$ and $\widehat D_\bfk(\eta)$ can be expanded in terms of their annihilation and creation operators with the wave functions
\ba
\label{D-wave}
\sin{\theta} D_{1 k}(\eta)  = D_k(\eta) =
 \frac{i}{f \sqrt{2k}} \left( 1 - \frac{i}{k\eta} \right) e^{-ik\eta} \, .
\ea

Now we present our decomposition of the tensor perturbations  $h_{ij}$ into $h_\times$ and $h_+$ polarizations following the method of  \cite{Ohashi:2013qba}.  Decomposing
$h_{ij}$ into $e_{ij}^{(s)}(\bfk)$ in Fourier space, the traceless  and transverse conditions, $h_{ii} = h_{ij,j} =0$,
yields
\ba
e_{i i}^{(s)}(\bfk) = 0 \quad , \quad  k_j e_{i j}^{(s)}(\bfk) = 0  \, ,
\ea
with $s= \times, +$ representing the two polarizations. In addition, we choose the following normalization
\ba
e^{(s)}_{ij}(\mathbf{k}) e^{*(s')}_{ij}(\mathbf{k}) = \delta_{ss'} \, ,
\ea
where $*$ represents the complex-conjugation. Note that we also have $e^{(s)}_{ij}(\mathbf{k}) = e^{*(s)}_{ij}(\mathbf{-k})$.

The  quantum operators $\widehat{h}_{ij}(\mathbf{k},\eta) $ in Fourier space
are represented in terms  of the annihilation and creation operators by
\ba
\label{tensor}
\widehat{h}_{ij}(\mathbf{k},\eta) = \sum _{s=+,\times}  \widehat{h}_{s}(\mathbf{k},\eta)
e_{ij}^{(s)}(\bfk)
 \quad , \quad \widehat{h}_{s}(\mathbf{k},\eta)=
h_{s}(k, \eta)a_{s}(\mathbf{k})+ h^{*}_{s}(k, \eta)a^{\dag}_{s}(-\mathbf{k}) \, ,
\ea
with  the commutation relations  $ [ a_{s}(\bfk),    a_{s}^\dagger  (\bfk')] = \delta_{s s'} \delta^{(3)} (\bfk-\bfk' )$.

As we mentioned before, we chose the convention that
$\mathbf{k} = k(\cos{\theta}, \sin{\theta}, 0)$. With this choice, the polarizations
$e^{+}_{ij}(\mathbf{k})$ and $e^{\times}_{ij}(\mathbf{k})$ become
\begin{align}
\label{polarization2}
e^{+}_{ij}(\mathbf{k}) =\frac{1}{\sqrt{2}} \left( \begin{array}{ccc}
 \sin^2{\theta} & -\sin{\theta}\cos{\theta} & 0 \\
  -\sin{\theta}\cos{\theta} & \cos^2{\theta} & 0 \\
  0 & 0 & -1 \\
\end{array} \right) ~~~,~~~
e^{\times}_{ij}(\mathbf{k}) = \frac{i}{\sqrt{2}}\left( \begin{array}{ccc}
  0 & 0 & -\sin{\theta} \\
  0 & 0 & \cos{\theta} \\
  -\sin{\theta} & \cos{\theta} & 0 \\
\end{array} \right) ~.
\end{align}
Using Eq. (\ref{tensor}) and Eq. (\ref{polarization2}), we find the following expression for the Fourier mode of the tensor field
\begin{align}
\label{polarization3}
\widehat{h}_{ij}(\mathbf{k}) =\frac{1}{\sqrt{2}} \left( \begin{array}{ccc}
 \widehat{h}_{+}\sin^2{\theta} & -\widehat{h}_{+}\sin{\theta}\cos{\theta} & -i\widehat{h}_{\times}\sin{\theta} \\
  -\widehat{h}_{+}\sin{\theta}\cos{\theta} & \widehat{h}_{+}\cos^2{\theta} & i\widehat{h}_{\times}\cos{\theta} \\
  -i\widehat{h}_{\times}\sin{\theta} &  i\widehat{h}_{\times}\cos{\theta}& -\widehat{h}_{+} \\
\end{array} \right) \, .
\end{align}
We will use this expression later on when  calculating the cross-correlation between the tensor mode and the curvature as well as the gauge field.

The profile of the tensor excitations has the standard form
\ba
\label{hs-wave}
{h}_{s}(k,\eta) = \frac{2 i H\eta }{M_{P}\sqrt{2 k}}\left(1-\frac{i}{ k\eta} \right)e^{-ik\eta} ~~~,~~~ (s= +, \times) \, .
\ea
The power spectrum of the tensor perturbations is given by
\ba
\langle   \widehat{h}_{ij}(\mathbf{k_1})  \widehat{h}_{ij}(\mathbf{k_2}) \rangle
= (2 \pi)^3 \delta^{(3)} (\bfk_1 + \bfk_2) P_h(k_1) \quad , \quad
\calP_h \equiv   \frac{ k_1^3  }{ 2 \pi^2}  P_h(k_1)
\ea
In the absence of anisotropy the power spectrum has the standard form
\ba
\label{calPh0}
\calP_h^{(0)} = \frac{2 H^2}{\pi^2 M_P^2}  = 16 \epsilon_H  \calP_\zeta^{(0)} \, .
\ea
Therefore, defining the tensor-to scalar ratio $r \equiv \calP_h/\calP_\zeta  $ we have $r = 16 \epsilon_H$
for the isotropic theory.

\subsection{The Interaction Lagrangian}
\label{interactions}

Having presented the background in some details, here we separate the Lagrangian into the free field part and interaction part. Here and below we call the latter the interaction Lagrangian.
The starting Lagrangian from the action (\ref{action}) is
\ba
\label{interaction Lagrangian1}
L = -\frac{a^4}{4}f(\phi)^2 F_{\mu \nu} F^{\mu \nu} - \frac{a^4}{2}  \e^2 \phi^2 A_{\mu} A^{\mu} \, .
\ea
Expanding the above action around the background values,  neglecting the contributions of the non-dynamical field $\delta A_0$ which are sub-leading as discussed before, and using the relation
$\left(\frac{\partial f^2}{\partial \phi}\right)\delta \phi = 4 f^2 \zeta$,
the interaction Lagrangians in the Fourier space is calculated as (see Apendix for further details)
\ba
\label{Lzetah}
L_{\zeta h_+} &=&
 - \frac{3\sqrt{2}}{2} I \epsilon_{H} M_{P}^2\sin^2{\theta}a^2\left(-\eta\right)^{-2}\left( \zeta^{*} {h}_{+} + c.c.  \right)  + \frac{\e^2 \sqrt{2}}{6}  I \epsilon_{H}M_{P}^4 \sin^2{\theta} \left(\frac{a ^4}{f^2}\right)\left( \zeta^{*} {h}_{+} + c.c. \right)  \\
\label{LzetaA}
 L_{\zeta D_1} &=& -2 M_{P}\sqrt{3I \epsilon_{H}}\sin^2{\theta} \left(\frac{a f}{\eta}\right) \left( \zeta ^{*} D'_{1} + c.c.\right) - 2 \e^2 M_{P}^3 \sqrt{\frac{I \epsilon_{H}}{3}}\sin^2{\theta} \left(\frac{a^3}{f}\right) \left( \zeta ^{*} D_{1} + c.c.\right) \\
 \label{Lh+D1}
 L_{h_+  D_1} &=& \frac{M_{P}}{2} \sqrt{\frac{3I\epsilon_{H}}{2}}\sin^2{\theta}\left( \frac{fa}{\eta}\right) \left( D^{'*}_{1} {h}_{+} +  c.c.\right) + \sqrt{\frac{I}{6\epsilon_{H}}} \e^2 M_{P}^3 \sin^2{\theta}\left( \frac{a^3}{f}\right)\left( D^{*}_{1} {h}_{+} + c.c.\right) \\
\label{LhxzD}
 L_{h_\times  D} &=&
 \frac{M_{P}}{2} \sqrt{\frac{3I \epsilon_{H}}{2}}\sin{\theta} \left(\frac{fa}{\eta}\right)\left( i D^{'} {h}^{*}_{\times} + c.c.\right) + \sqrt{\frac{I}{6\epsilon_{H}}}\e^2 M_{P}^3 \sin{\theta} \left(\frac{a^3}{f}\right)\left( i D{h}^{*}_{\times} + c.c. \right)
\ea
where c.c stands for complex conjugation.

The above interaction Lagrangians are needed in order to calculate the anisotropy corrections in
$\langle \zeta \zeta \rangle,  \langle h_s  h_{s'} \rangle$ and the cross-correlations $\langle \zeta h_s \rangle$.  Note that in the free (isotropic) theory with $I=\e =0$ there is no anisotropy corrections in power spectra and $\langle \zeta h_s \rangle$ as expected.

\section{Anisotropic Correlations  }
\label{correlations}

Having calculated the interaction Lagrangians as given in Eqs. (\ref{Lzetah})-(\ref{LhxzD})
now we are ready to calculate the anisotropic correlation functions by using the in-in formalism.
For this purpose, we need to obtain the interaction Hamiltonian from the interaction Lagrangian.
One should notice that $H_{int} = -L_{int}$ is not necessary true with kinetically coupled interactions.
So it is worth to check it in this model before proceeding with the in-in calculation of the correlation functions. We have calculated it in the Appendix \ref{int-hamilton}. It turns out that the above formula is true for the whole of the interactions except $H_{\zeta h_{+}}$. So in the following we use $H_{int} = -L_{int}$ everywhere except that in $H_{\zeta h_{+}}$, special care is taken of.

We are interested in anisotropic contributions in $\langle \zeta_\bfk \zeta_\bfk^* \rangle  $, $\langle h_{s \, \bfk} h_{s' \, \bfk}^* \rangle  $  and $ \langle \zeta_{ \bfk} h_{s' \, \bfk}^* \rangle$. We calculate each term in turn.
Note that the wave function of the free theory for $\zeta_k, D_k, D_{1 k}$ and $h_{s k}$ are given in Eqs. (\ref{zeta-wave}),  (\ref{D-wave}) and (\ref{hs-wave}).

\subsection{Anisotropies in curvature power spectrum}

\label{zeta-power}

Here we calculate the anisotropic contributions in curvature perturbation power spectrum $\langle \zeta_\bfk \zeta_\bfk^* \rangle  $.  We denote the change in curvature perturbation power spectrum from the anisotropic sources by $\delta \langle \zeta_\bfk \zeta_\bfk^* \rangle  $.  This analysis were performed in  \cite{Emami:2013bk} and here we outline the analysis briefly.

We use the in-in formalism to take care of the corrections from the bilinear coupling terms \cite{Chen:2009zp, Weinberg:2005vy, Chen:2010xka, Wang:2013zva}. The leading order corrections in curvature perturbation power spectrum are given by
\ba
\label{delta-P-zeta}
\delta \langle \zeta_\bfk \zeta_\bfk^* \rangle  = - \int_{\eta_{0}}^{\eta_{e}} d\eta_{1} \int_{\eta_{0}}^{\eta_{1}}d\eta_{2} \left \langle
\bigg{[} L_{I}(\eta_{2}) , \bigg{[} L_{I}(\eta_{1}) ,
 \zeta_\bfk (\eta_e) \zeta_\bfk^*(\eta_e)  \bigg{]}\bigg{]} \right \rangle  \, ,
\ea
where $L_I$ represents the interaction Lagrangian.  The lower limit of the integral should be set
$k \eta_0 \rightarrow -\infty$  corresponding to initial  modes being deep inside the horizon. However, as studied in \cite{Bartolo:2012sd, Emami:2013bk}, the interactions responsible for anisotropies operate on super-horizon scales so to a good approximation one can safely take $k \eta_0 =-1$ corresponding to the time when the mode leaves the horizon.
The upper limit of the above integral as usual corresponds to $k \eta_e \simeq 0$.

The interaction Lagrangians relevant to
$\delta \langle \zeta_\bfk \zeta_\bfk^* \rangle$ are $L_{\zeta h_+}$ and  $L_{\zeta D_1}$  given
in Eqs. (\ref{Lzetah}) and (\ref{LzetaA}).  A look at these two equations show that $L_{\zeta h_+}$
is suppressed compared to $L_{\zeta D_1}$ by the factor $\sqrt{I \eH} \ll 1$. Therefore, the leading order anisotropic corrections in  curvature perturbation power spectrum comes from  $L_{\zeta D_1}$.
In addition, $L_{\zeta D_1}$ has two independent terms denoted by $L^{(1)}_{\zeta D_1}$ and
$L^{(2)}_{\zeta D_1}$:
\ba
L^{(1)}_{\zeta D_1} \equiv
 -2 M_{P}\sqrt{3I \epsilon_{H}}\sin^2{\theta} \left(\frac{a f}{\eta}\right) \left( \zeta ^{*} D'_{1} + c.c.\right)
  \quad , \quad
L^{(2)}_{\zeta D_1} \equiv
- 2 \e^2 M_{P}^3 \sqrt{\frac{I \epsilon_{H}}{3}}\sin^2{\theta} \left(\frac{a^3}{f}\right) \left( \zeta ^{*} D_{1} + c.c.\right)
\ea
Depending on whether one chooses either $L^{(1)}_{\zeta D_1}$ or $L^{(2)}_{\zeta D_1}$ in place of  $L_I(\eta_1)$ and $L_I(\eta_2)$ in the   integral in Eq. (\ref{delta-P-zeta}), there are four possible contributions in  $\delta \langle \zeta_\bfk \zeta_\bfk^* \rangle$ denoted by $\delta \langle \zeta_\bfk \zeta_\bfk^* \rangle_{ij}$ where $i, j=1, 2$ with the  assumption  that
$L_I(\eta_1) = L^{(i)}_{\zeta D_1}$  and  $L_I(\eta_2) = L^{(j)}_{\zeta D_1}$.
For example,    $\delta \langle \zeta_\bfk \zeta_\bfk^* \rangle_{12}$ means
$L_I(\eta_1) = L^{(1)}_{\zeta D_1}$  and  $L_I(\eta_2) = L^{(2)}_{\zeta D_1}$. With this identification we have
\ba
\label{correction zeta}
\delta  \bigg{\langle} \zeta_{\bfk}(\eta_{e}) \zeta_{\bfk}^*(\eta_{e})\bigg{\rangle} =
\delta \bigg{\langle} \zeta_{\bfk}(\eta_{e}) \zeta_{\bfk}^*(\eta_{e})\bigg{\rangle} _{11}
+\delta \bigg{\langle} \zeta_{\bfk}(\eta_{e}) \zeta_{\bfk}^*(\eta_{e})\bigg{\rangle} _{12}
+ \delta \bigg{\langle} \zeta_{\bfk}(\eta_{e}) \zeta_{\bfk}^*(\eta_{e})\bigg{\rangle} _{21}
+ \delta \bigg{\langle} \zeta_{\bfk}(\eta_{e}) \zeta_{\bfk}^*(\eta_{e})\bigg{\rangle} _{22}
\ea
The details of the in-in analysis are presented in Appendix \ref{in-in}. As a sample analysis, here we present the integral form of
$\delta  \bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) \zeta_{\mathbf{k}_{1}}(\eta_{e})^*\bigg{\rangle}_{11} $
which is  (here and hence after, the momentum conservation $\delta$-function is omitted to save writing)
\ba
\delta  \bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) \zeta_{\mathbf{k}_{1}}^*(\eta_{e})\bigg{\rangle}_{11}  &&=  384 I \eH M_P^2  \sin^4{\theta} \
\int_{\eta_{0}}^{\eta_{e}} d\eta_{1} \left(\frac{a f}{\eta} \right)_{\eta_1}
\im \left[ \zeta_k (\eta_{1})\zeta_k^{*}(\eta_{e})\right]  \\
&& ~~~~~~~~~~~~~~~~~~~~~~~\times \int_{\eta_{0}}^{\eta_{1}}d\eta_{2}
 \left(\frac{a f}{\eta} \right)_{\eta_2}
\im \left[ \zeta_k(\eta_{2})\zeta_k^{*}(\eta_{e})  D_{1k}^{' *} (\eta_{1}) D_{1k}^{'}(\eta_{2})
\right] \, .
\ea
Expanding the integrand for small $k \eta$ arguments and assuming $k \eta_0 = -1$ and
$k \eta_e =0$ as explained above, the above integral yields
\ba
\delta  \bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) \zeta_{\mathbf{k}_{1}}(\eta_{e})\bigg{\rangle}_{11}&=& \frac{6 I N^2}{k_{1}^3 \epsilon_{H}} \left(\frac{H}{M_{P}} \right)^2 \sin^2{\theta}
\ea
in which $N= - \ln (-k \eta_e)$ represents the number of e-folds when the mode
$k$ has left the horizon. Taking $k$ to be the CMB scales we need $N \sim 60$ in order to solve
the flatness and the horizon problem.

Performing the same procedure for other integrals, we obtain  \cite{Emami:2013bk}
\ba
\delta  \bigg{\langle} \zeta_{\mathbf{k}}(\eta_{e}) \zeta_{\mathbf{k}}^*(\eta_{e})\bigg{\rangle}_{12}&=& -\frac{31 }{490}  \frac{\e^2 I}{k^3 {\epsilon_{H}}} \sin^2{\theta} \\
\delta  \bigg{\langle} \zeta_{\mathbf{k}}(\eta_{e}) \zeta_{\mathbf{k}}^*(\eta_{e})\bigg{\rangle}_{21}&=& -\frac{ I \e^2 N}{7 k^3 \epsilon_H}   \sin^2{\theta} \\
\delta  \bigg{\langle} \zeta_{\mathbf{k}}(\eta_{e}) \zeta_{\mathbf{k}}^*(\eta_{e})\bigg{\rangle}_{22}&=& \frac{9}{2156}  \frac{\e^4 I}{k_{1}^3}\left(\frac{M_{P}}{m} \right)^2 \sin^2{\theta}
\ea
So combining these four contributions, and assuming $N \gg 1$,  we have
\ba
\label{correction zeta final}
\delta  \bigg{\langle} \zeta_{\mathbf{k}}(\eta_{e}) \zeta_{\mathbf{k}}^*(\eta_{e})\bigg{\rangle}
&\simeq&
\left( \frac{6 I N^2}{ \epsilon_{H}} \frac{H^2}{M_{P}^2}
-\frac{ I \e^2 N}{7  \epsilon_H}
 +\frac{9 \, \e^4 I}{2156} \frac{M_{P}^2}{m^2}  \right) \left(\frac{\sin^2{\theta}}{k^3}\right)\\
 &=&   \frac{6 I N^2}{ \epsilon_{H}} \frac{H^2}{M_{P}^2}    \frac{\sin^2{\theta}}{k^3} F (\beta)
\ea
where we have defined
\ba
\label{beta-def}
\beta \equiv \frac{\e^2 }{42 N}  \left(\frac{M_P}{H} \right)^2  \quad , \quad
F(\beta) \equiv 1- \beta + \frac{9}{22} \beta^2 \, .
\ea
We have also written $m$ in terms of $H$ and $\epsilon_H$ by using $m^2 = \left(3 \epsilon_H  H^2 \right)$.\\
Note that $\beta$ is a measure of the gauge field coupling $\e^2$. In particular, in the model of \cite{Watanabe:2010fh, Bartolo:2012sd, Ohashi:2013qba} with $\e=0$ we have  $F(\beta)=1$.
With $M_P/H \sim 10^{5}$, and with $\e \gtrsim 10^{-4}$ we obtain $\beta \gtrsim 1$. For larger value of $\e$ we see that  $F(\beta)$ grows  like $\beta^2$.

The anisotropic power spectrum $\delta \calP_\zeta= \frac{k^3}{2 \pi^2} \delta {\langle} \zeta_{\mathbf{k}}(\eta_{e}) \zeta_{\mathbf{k}}^*(\eta_{e}) {\rangle}$ therefore is
\ba
\label{delta-calP-zeta}
\delta \calP_\zeta= \frac{3 I N^2 H^2 }{\eH \pi^2 M_P^2} F(\beta) \sin^2 \theta \, .
\ea
Correspondingly, the total anisotropic power spectrum  $\calP_\zeta$ is
\ba
\calP_\zeta =  \calP_\zeta^{(0)} \left[ 1+ 24 I N^2  F(\beta)  \sin^2 \theta
\right] \, ,
\ea
where $\calP_\zeta^{(0)}$ represents the isotropic power spectrum for the free theory. Note that in the limit  when $\e=\beta=0$ so $F(\beta)=1$, our result for  $\delta \calP_\zeta$ agrees with the result
in  \cite{Watanabe:2010fh, Bartolo:2012sd, Ohashi:2013qba}.

Now defining the anisotropy estimator $g_*$ via
\ba
\label{g*}
\calP_\zeta =  \calP_\zeta^{(0)} \left[ 1+  g_* \left( \widehat\bfk \cdot \widehat\bfp \right)^2 \right]
\ea
where $\bfp$ is the preferred anisotropic direction in the sky (the $x$-direction in our example), we
obtain
\ba
\label{g*2}
g_* = -24 I N^2 F(\beta) \,
\ea
It is important to note that the form of $g_*$ we have defined here is with respect to the primordial curvature power spectrum. In the TT and other correlation functions, the anisotropy not only comes from the $g_*$ here, but also comes from an ``effective $g_*$'' contribution from the tensor sector. Such an ``effective $g_*$'' has a scale dependence on the TT and other correlations because the tensor mode is decaying after it returns to the horizon. We will return to this issue later.
Also we see that the gauge coupling $\e$ appears  in $g_*$ via the parameter $\beta$.
Taking $| g_*| \lesssim 10^{-2}$ from the Planck data  constraint  \cite{Kim:2013gka} we require
$I F(\beta) \lesssim 10^{-6}$.

One may ask what the theoretical limits on the value of $\e$ or the parameter $\beta$ are. First, we have to make sure that we get enough number of e-folds of inflation at the background level. As we mentioned before, the number of e-folds depends logarithmically on $\e$
so as studied in \cite{Emami:2010rm} one can take say $\e <0. 1$ to get  a long enough
period of inflation. In addition, our assumption in parametrizing the anisotropy was that
the anisotropic  power spectrum  is smaller than the isotropic power spectrum, i.e.
$|g_*| < 1$ so our perturbative approach using the leading order in-in formalism is valid.
Therefore, demanding $|g_*| < 1$ we need  $I F(\beta) < 10^{-4}$. We have presented the contour plot of the allowed range of $I$ and $\e$ in Fig. \ref{Ie-plot}. As can be seen, we need $\e \lesssim 10^{-3}$ in order not to produce too much anisotropy in tensor perturbations (to be discussed in next subsection). Therefore, with $\e \lesssim 10^{-3}$, and $M_P/H \sim 10^5$ we have $\beta \lesssim 3$.

Checking the behavior of the function $F(\beta)$ for the approximate allowed range
$\beta \lesssim 3$ indicates that  $g_*$ has  a weak dependence on $\e$.
One can check that  for $0\le  \beta \le 3$, $F(\beta)$
takes the value in the range $0.4 \lesssim    F(\beta) \lesssim 1.7 $.  As we shall see this conclusion plays important roles for the predictions of our model for various cross-correlations.  While the anisotropy in curvature perturbation is under control for the above range  of $\beta$, the tensor perturbations become highly anisotropic when $\beta \gtrsim 1$.

Here we pause to mention one conceptual problem associated with anisotropic inflation. As seen from
Eq. (\ref{g*2}) the amplitude of anisotropy in scalar power spectrum scales like $N^2$. If inflation is prolonged in the past, this yields a large value of $g_*$ and the system becomes highly anisotropic. Therefore, our treatment of taking the anisotropies as small corrections to the isotropic
FRW background will be invalid. As pointed out in  \cite{Bartolo:2012sd} this corresponds to IR gauge field fluctuations which have left the horizon in the past inflationary history and contributed to the classical background trajectory. Therefore, a prolonged period of inflation will bring more contributions from these anisotropic IR modes which can destroy the near isotropy of the background.
In order to prevent this to happen, we demand that the total period of anisotropic inflation is under control, say less than few hundreds of e-folds.

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.6\textwidth]{fig/param_range}
  \caption{The  allowed range of $I$ and $\e$ for different values of $g_*$.  The shaded region
  is beyond our current scope because a large value of $\e$ induces too much anisotropy in tensor power spectrum
  which undermines our perturbative approach, i.e. $\delta \calP_h $ in Eq. (\ref{delta-Ph})
  becomes comparable to the   isotropic power spectrum $\calP_h^{(0)}$. In addition, note that we require $\e < 0.1$ in order to get large enough number of e-folds. As a result the range  $\e < 0.1$ is consistent both at the background and perturbation levels.}
  \label{Ie-plot}
\end{figure}

\subsection{Correction to the tensor power spectrum}

Now we calculate the anisotropy in tensor power spectra
 $\langle {h}_{\bfk \times } {h}_{\bfk \times }^*\rangle$ and $\langle {h}_{\bfk + } {h}_{\bfk + }^*\rangle$.

Let us  start with  $\langle {h}_{\times } {h}_{\times }^*\rangle$.  The relevant  interaction Lagrangian is
$L_{D {h}_{\times } } =  L_{D {h}_{\times } }^{(1)} +  L_{D {h}_{\times } }^{(2)} $ where
$ L_{D {h}_{\times } }^{(1)}$ and  $ L_{D {h}_{\times } }^{(2)}$ respectively are the first and the second terms in Eq. (\ref{LhxzD}).   Following the same convention as  in our analysis for
anisotropies in curvature perturbation power spectrum we have
\ba
\delta \bigg{\langle} {h}_{\times } {h}_{\times }^*\bigg{\rangle} &=& - \int_{\eta_{0}}^{\eta_{e}} d\eta_{1}\int_{\eta_{0}}^{\eta_{1}} d\eta_{2}\bigg{[}L_{D {h}_{\times }  } ,\bigg{[}L_{D {h}_{\times}}, {h}_{\times\mathbf{k}} {h}_{\times\mathbf{k}}\bigg{]}\bigg{]} \nonumber\\
&=& \delta \bigg{\langle} \widehat{h}_{\times }\widehat{h}_{\times }\bigg{\rangle}_{11} +  \delta \bigg{\langle} \widehat{h}_{\times }\widehat{h}_{\times }\bigg{\rangle}_{12} + \delta \bigg{\langle} \widehat{h}_{\times }\widehat{h}_{\times }\bigg{\rangle}_{21} + \delta \bigg{\langle} \widehat{h}_{\times }\widehat{h}_{\times }\bigg{\rangle}_{22}
\ea
The details of the in-in analysis are presented in Appendix \ref{in-in}. The results for each contribution are
\ba
\label{tensor times 1}
\delta \bigg{\langle} {h}_{\times } {h}_{\times }^*\bigg{\rangle}_{11}
&=& \left(\frac{12}{k^3}\right)\left(\frac{H}{M_{P}}\right)^2 I \epsilon_{H} N^2\sin^2{\theta}
\\
\label{tensor times 2}
\delta \bigg{\langle} {h}_{\times}  {h}_{\times}^*\bigg{\rangle}_{12}&=&  -\left(\frac{4}{7 k^3}\right)N I\e^2\sin^2{\theta}\\
\label{tensor times 3}
\delta \bigg{\langle} {h}_{\times}  {h}_{\times}^*\bigg{\rangle}_{21}&=&  -\left(\frac{62}{245 k^3}\right) I \e^2 \sin^2{\theta}
\\
\label{tensor times 4}
\delta \bigg{\langle} {h}_{\times}  {h}_{\times}^*\bigg{\rangle}_{22}
&=& \left(\frac{6 }{539 k^3 }\right) \left(\frac{M_{P}}{H}\right)^2 \left(\frac{I \e^4}{\epsilon_{H}}\right) \sin^2{\theta}
\ea

Summing up those four terms we get
\ba
\label{correction hh times final}
\delta \bigg{\langle} \widehat{h}_{\times\mathbf{k}_{1}}  \widehat{h}_{\times\mathbf{k}_{2}}\bigg{\rangle}
&\simeq& \left( 12I \epsilon_{H}N^2 \frac{H^2}{M_{P}^2} -\frac{4}{7}N I\e^2 + \frac{6 I \e^4}{539 \eH} \frac{M_{P}^2 }{H^2 } \right) \frac{\sin^2{\theta}}{k^3} \nonumber\\
&=&  12I \epsilon_{H}N^2 \frac{H^2}{M_{P}^2} F(\hat \beta)  \frac{\sin^2{\theta}}{k^3}
\ea
where the function $F(x)$ is defined in Eq. (\ref{beta-def}) and $\hat \beta$ is related to
$\beta$ via
\ba
\hat \beta \equiv \frac{2 \beta}{\epsilon_H} \, .
\ea
Note the crucial point that $\hat \beta$ is enhanced compared to $\beta$ by the factor $1/\epsilon_H$.
As we have discussed before, from the observational and the theoretical constraints on $\beta$ we have
$\beta \lesssim 1$. Now, from the above relation between $\hat \beta$ and $\beta$, we see that
$\hat \beta$ can be as large as $100$ for $\epsilon_H \sim 0.01$. As we shall see shortly, the anisotropy in tensor power spectra becomes very strong for large $\e$ so there will be  upper bound on $\e$ and
$\hat \beta$.

Now we calculate   $\langle {h}_{+ }{h}_{+}^*\rangle$. In this case the relevant interaction Lagrangians are  $L_{D_{1} {h}_{+ } }$ and $ L_{\zeta {h}_{+ }  }$ so we have
\ba
\label{tensor plus}
\delta \bigg{\langle} {h}_{+} {h}_{+}\bigg{\rangle} &=& - \int_{\eta_{0}}^{\eta_{e}} d\eta_{1}\int_{\eta_{0}}^{\eta_{1}} d\eta_{2}\bigg{[}L_{D_{1} {h}_{+ }  } ,\bigg{[}L_{D_{1} {h}_{+}}, {h}_{+} {h}_{+}\bigg{]}\bigg{]} - \int_{\eta_{0}}^{\eta_{e}} d\eta_{1}\int_{\eta_{0}}^{\eta_{1}} d\eta_{2}\bigg{[}L_{\zeta {h}_{+ }  } ,\bigg{[}L_{\zeta {h}_{+}}, {h}_{+} {h}_{+}\bigg{]}\bigg{]}
\ea
Comparing $L_{D_{1} {h}_{+ } }$ and $ L_{\zeta {h}_{+ }  }$ we see that $ L_{\zeta {h}_{+ }  }$ is suppressed compared to $L_{D_{1} {h}_{+ } }$ by a factor $\sqrt{I} \ll 1$
so to leading order in anisotropy we can neglect the contribution from  $ L_{\zeta {h}_{+ }  }$ in
$\langle {h}_{+ }{h}_{+}^*\rangle$. Now the analysis is exactly the same as what we performed
in  $\langle {h}_{\times }{h}_{\times}^*\rangle$ and therefore
\ba
\label{tensor plus final}
\delta \bigg{\langle} {h}_{+} {h}_{+}\bigg{\rangle} =
\delta \bigg{\langle} {h}_{\times} {h}_{\times}\bigg{\rangle} \, .
\ea

To summarize,  the anisotropy in total tensor power spectrum is
\ba
\label{delta-Ph}
\delta \calP_h &=& 2 \left( \frac{k^3}{2 \pi^2}   \right)
  \delta \bigg{\langle} {h}_{\times} {h}_{\times}^*\bigg{\rangle} \nonumber\\
&=&   24 I \epsilon_{H}N^2 \frac{H^2}{M_{P}^2} F(\hat \beta)  \sin^2{\theta} ~.
\ea
so the total tensor power spectrum is
\ba
\calP_h = \calP_h^{(0)} \left( 1 +  6   I \epsilon_{H}N^2 F(\hat \beta) \sin^2 \theta ~.
\right)
\ea
This is an interesting formula indicating that the effects of the gauge coupling is very strong
in tensor power spectrum anisotropy. This is because $\hat \beta = 2 \beta/\epsilon_H$ so with
$\beta \sim 1$ we gate  $\hat \beta \sim 100$  and therefore
$\delta \calP_h /\calP_h^{(0)} \simeq 24 I N^2 (\beta^2/\epsilon_H) = |g_*| \beta^2/\epsilon_H$.
With $| g_*| \sim O(\epsilon_H)$, which is consistent with the observational constraints and with
$\beta \sim 1$, one easily gets to the regime in which  $\delta \calP_h /\calP_h^{(0)} \simeq1 $.
Note that the epsilon enhancement for the charged interaction is a very special feature of this model, from the specific form of the potential. Explicitly, since our interaction includes $\e^2 \phi^2 A^2$, the charged contribution in $\langle \zeta \zeta\rangle$ comes from $(\e^2 \phi A \delta\phi \delta A )^2$, while the charged contribution in $\langle hh\rangle$ comes from $(\e^2 \phi^2 A h_{ij} \delta A)^2$. So we see that the ratio between these two effects controls with $\phi^2$.
In the chaotic inflation, $\phi^2$ is proportional to $1/\epsilon$. \footnote{On the other hand, for the symmetry breaking potential, the ratio would be $\epsilon$, so no hope to see any enhancement for this case.}
This signals the strong dependence of the tensor anisotropies to the gauge coupling. Of course, we can not trust our analysis when we approach the limit $\delta \calP_h /\calP_h^{(0)} \simeq1 $. This is because, we have followed a perturbative approach and only kept the leading interaction terms in
our in-in analysis. Our situation  is in contrast to models of anisotropic inflation with a real inflaton field, as studied in  \cite{Ohashi:2013qba}, where $\e= \beta=0$,  so $F(\beta) =1$ and
$\delta \calP_h /\calP_h^{(0)} = - g_* \eH /4$ which is highly suppressed.

Demanding that $\delta \calP_h < \calP_h^{(0)}$ so our theoretical analysis is under perturbative control, we obtain the following upper bound
on the parameter $\beta$
\ba
\label{beta-cons}
\beta \lesssim  \sqrt{\frac{\eH}{ |g_* | }} \, .
\ea
With $\epsilon_H \sim 10^{-2}$  and $| g_*| \lesssim 10^{-2}$ we conclude that $\beta \lesssim 1$
in order for the anisotropic contribution in tensor power spectrum to be under control, corresponding to  $\e \lesssim 10^{-3}$.  The contour plot of  $I $ versus $\e$ is shown in Fig. \ref{Ie-plot}. The strong constraints on the allowed range of $\e$ comes from the tensor power spectrum.

\subsection{Cross-Correlation between $\zeta$ and $h_{ij}$}

Here we calculate the cross-correlation $\zeta$ and $h_{ij}$. We should calculate the following terms,
\ba
\bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) {h}_{+\mathbf{k}_{2}}(\eta_{e}) \bigg{\rangle} &=& i \int_{\eta_{0}}^{\eta_{e}} d\eta_{1} \bigg{\langle}   \bigg{[} H_{\zeta {h}_{+}} , \zeta_{\mathbf{k}_{1}}{h}_{+\mathbf{k}_{2}}\bigg{]}   \bigg{\rangle}
- \int_{\eta_{0}}^{\eta_{e}} d\eta_{1}\int_{\eta_{0}}^{\eta_{1}} d\eta_{2}
\bigg{\langle} \bigg{[}L_{\zeta D_{1}} ,\bigg{[}L_{D_{1} {h}_{+}}, \zeta_{\mathbf{k}_{1}}{h}_{+\mathbf{k}_{2}}\bigg{]}\bigg{]}  \bigg{\rangle} \nonumber\\
&& -\int_{\eta_{0}}^{\eta_{e}} d\eta_{1}\int_{\eta_{0}}^{\eta_{1}} d\eta_{2}
\bigg{\langle} \bigg{[} L_{D_{1} {h}_{+}},\bigg{[} L_{\zeta D_{1}}, \zeta_{\mathbf{k}_{1}}{h}_{+\mathbf{k}_{2}}\bigg{]}\bigg{]}  \bigg{\rangle}
\nonumber\\
&\equiv& \bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) {h}_{+\mathbf{k}_{2}}(\eta_{e})\bigg{\rangle}_{1} + \bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) {h}_{+\mathbf{k}_{2}}(\eta_{e})\bigg{\rangle}_{2} + \bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) {h}_{+\mathbf{k}_{2}}(\eta_{e})\bigg{\rangle}_{3}
\ea
where the indices 1, 2 and 3 indicate the above three integrals respectively. The nested integrals 2 and 3 each have four different contributions as in previous analysis so in total we have nine contributions in the above cross correlation.  The details of the analysis are given in Appendix \ref{in-in} and here we present the final result:
\ba
\label{final cross correlation}
\bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) {h}_{+\mathbf{k}_{2}}(\eta_{e})\bigg{\rangle} &\simeq&
I \left(- 6\sqrt{2} \frac{ N^2 H^2 }{M_{P}^2}  + \frac{  \sqrt{2} \e^2 N}{7 \epsilon_{H} }  - \frac{3\sqrt{2} \e^4 }{1078 \epsilon_{H} } \frac{M_{P}^2}{H^2}  \right) \frac{\sin^2{\theta}}{k^3}
\ea
Eq. (\ref{final cross correlation}) is the final result for the cross-correlation.

Finally one can easily check that
\ba
\bigg{\langle} \zeta_{\mathbf{k}_{1}}(\eta_{e}) {h}_{\times\mathbf{k}_{2}}(\eta_{e})\bigg{\rangle}=0 \, .
\ea
This is because at the second order level $\zeta $ does not see ${h}_{\times}$.

The power spectrum of $\langle \zeta h \rangle$ cross-correlation is therefore
\ba
\label{calP-zh}
\calP_{\zeta h} &=& \frac{k^3}{2 \pi^2 } \bigg{\langle} \zeta_{\mathbf{k}}(\eta_{e}) {h}_{+\mathbf{k}^*}(\eta_{e})\bigg{\rangle}  \nonumber\\
&=& -24 \sqrt{2}  I N^2 \eH \calP_\zeta^{(0)} G(\beta) \sin^2 \theta
\ea
where the function $G(\beta)$ is defined via
\ba
\label{G-beta}
G(\beta) \equiv 1- \frac{\beta }{\epsilon_H} + \frac{9}{11} \frac{\beta^2 }{\epsilon_H} \, .
\ea
For typical value of $\epsilon_H \ll 1$, the function  $G(\beta)$ has two positive roots $\beta_1$ and $\beta_2$ where  $\beta_1  \ll 1$ and $\beta _2 \gtrsim  1$. Thus the function $G(\beta)$ is negative in the range $\beta_1 < \beta < \beta_2$ while it is positive beyond this region. Therefore, with appropriate
choice of $\e$ or $\beta$ the cross-correlation $\calP_{\zeta h} $ can have both signs, i.e. $\zeta$ and $h$
can be either correlated or anti-correlated.

Alternatively, one can also write $\calP_{\zeta h} $ in terms of $g_*$ as (note that $g_* <0$)
\ba
\calP_{\zeta h} =  \sqrt{2} g_*  \eH \calP_\zeta^{(0)} \frac{G(\beta)}{F(\beta)} \, .
\ea
Note that in the limit where $\e= \beta=0$ so $F(\beta) = G(\beta) =1$, the above expression coincides with the result obtained in  \cite{Ohashi:2013qba}. With $\beta \sim1$ the ratio $\calP_{\zeta h}/\calP_{\zeta}^{(0)} $ in our model is about one or two orders of magnitude bigger than the result in
 \cite{Ohashi:2013qba} in which $\beta =0$.

Before closing this Section and presenting our numerical results for various correlations, let us summarize the main results of our model.
The anisotropy in curvature perturbation power spectrum is given by Eq. (\ref{delta-calP-zeta})
with $g_*$ given in Eq. (\ref{g*}). As discussed below  Eq. (\ref{g*}), $\delta \calP_\zeta$ and $g_*$ depend weakly on $\beta$ so we do not get strong constraints on the value of $\e$ from the constraints on curvature perturbations anisotropies. On the other hand, the anisotropic tensor power spectrum is given in Eq. (\ref{delta-Ph}). The crucial point is that $\delta \calP_h$ scales
with $F (\hat \beta)$ where $\hat \beta = 2 \beta/\epsilon_H$. With $\beta \sim 1$  we get
$\hat \beta \sim 100$ which yields  an enhancement $\sim 10^4$ from the function  $F (\hat \beta)$.
This is a novel effect indicating that while the scalar perturbations are well-constrained to be statistically symmetric, the tensor perturbations show strong directional dependence. This is the motivation for careful scrutiny of B-mode polarizations for the TB, EB and BB correlations in the upcoming Planck polarization maps.
Finally the cross-correlation of scalar-tensor,  $\calP_{\zeta h}$,  is given in Eq. (\ref{calP-zh}). The situation here is a hybrid of the above two limits of $\delta \calP_h$ and $\delta \calP_\zeta$.
For large enough value of $\e$, i.e. with $  \beta \sim 1$, we get an enhancement of order 10-100 from
the function $G(\beta)$.

\section{From the primordial fluctuations to the CMB}
\label{numerical}
In this section, we shall relate the calculation of the above sections to CMB anisotropies, using the same method as in \cite{Watanabe:2010bu}. We use the spin weighted spherical harmonics in our analysis, \cite{Kamionkowski:1996zd, Hu:1997hp}.

To calculate the CMB anisotropies, one has to project the three-momentum onto two dimensional spherical harmonics. The $\int d^3k$ integral breaks up into two pieces with the presence of statistical anisotropy -- the radius part and the angular part. The correlators of CMB observables takes the form
\begin{align}\label{eq:aXaX}
\langle a^{X_1}_{l_1, m_1} a^{X_2}_{l_2, m_2}\rangle = 4 \pi
\int \frac{dk}{k} \Delta_{l_1}^{i_1 X_1}(k) \Delta_{l_2}^{i_2 X_2}(k)
\int  d\Omega ~  [{}_{i_1}Y^*_{l_1m_1}(\theta, \phi)] [{}_{i_2}Y_{l_2m_2}(\theta, \phi)] P^{i_1, i_2}(k, \theta, \phi)~,
\end{align}
where $X^i$ takes value (T, E, B), which are the temperature anisotropy, the E-mode and B-mode respectively. Here although the dimensionless power spectrum in principle depends on $k$, we approximate it as scale invariant. The effects on the plots are tinny\footnote{On the other hand, the scale dependence is needed once data analysis is to be performed for this model and the modification is straightforward if the scale dependence can be written in a factorizable form $P^{i_1, i_2}(k, \theta, \phi) = \sum_n f_n(k)g_n(\theta,\phi)$.}. Also, in our case there is no $\phi$-dependence for the power spectrum and the $\theta$-dependence takes the shape $P^{i_1, i_2} = P^{i_1, i_2}(\sin^2 \theta)$. With the absence of $\phi$-dependence, the momentum along the $\phi$-rotation is conserved and thus $\langle a^{X_1}_{l_1, m_1} a^{X_2}_{l_1, m_1}\rangle$ is non-vanishing only when $m_1 = m_2$.

However, the rotational symmetry on the $\theta$ direction is broken. As a result, the correlation $\langle a^{X_1}_{l_1, m_1} a^{X_2}_{l_1, m_1}\rangle$ is not restricted to $l_1=l_2$. Instead, other than those diagonal correlations, we also have $l_1=l_2\pm1$ for TB and EB, and $l_1=l_2\pm 2$ for TT, TE, EE and BB respectively.

In \eqref{eq:aXaX}, the $i_1$ and $i_2$ are indices indicating the spin of the component. And the ${}_{i}Y^*_{lm}(\theta, \phi)$ is the spin-$i$-weighted spherical harmonics. On the $P^{i_1, i_2}(\sin^2 \theta)$ side, the correlation functions on the spin bases takes the form
\begin{align}
  P^{0,0} = P^{\zeta\zeta}~, \quad P^{0,\pm2} = P^{\pm2,0} = \frac{1}{\sqrt2} P^{0,+}~, \quad
  P^{\pm2, \pm2} = \frac{1}{2} \left( P^{++} + P^{\times\times}\right)~, \quad P^{\pm2, \mp2} = \frac{1}{2} \left( P^{++} - P^{\times\times}\right)~.
\end{align}
Note that  $P^{\zeta\zeta} \equiv P_{\zeta}$, $P^{0,+} \equiv P_{\zeta h_+}$, $P^{++}\equiv P_{h_+}$ and $P^{\times\times} \equiv P_{h_\times}$ components are given by  equations \eqref{correction zeta final}, \eqref{final cross correlation}, \eqref{correction hh times final} and \eqref{tensor plus final} respectively.

Among those angular integrals $\int  d\Omega ~  [{}_{i_1}Y^*_{l_1m_1}(\theta, \phi)] [{}_{i_2}Y^*_{l_2m_2}(\theta, \phi)] P^{i_1, i_2}(\sin^2\theta)$, we would like to emphasize a particularly interesting integral:
\begin{align}
 \int  d\Omega ~  [{}_{0}Y^*_{lm}(\theta, \phi)] [{}_{2}Y_{lm}(\theta, \phi)] \sin^2\theta ~.
\end{align}
Note that this term picks up cross correlation of $\langle\zeta h_+\rangle$, and maps it onto TT or other CMB anisotropies. It is claimed in \cite{Contaldi:2014zua} that such a term, with anti-correlated $\langle\zeta h_+\rangle$ (which is indeed possible in our case), suppresses the TT power spectrum and thus reconciles the current tension between BICEP2 and Planck.

However, we shall show that this is unfortunately not true. By expanding ${}_{2}Y^*_{lm}(\theta, \phi)$ into the spin 0 spherical harmonics, the $\sin^2\theta$ factor is cancelled and the integral ends up to be
\begin{align}
  \langle T_l^\zeta T_l^h\rangle \propto \sum_{m=-l}^{l} \int  d\Omega ~  [{}_{0}Y^*_{lm}(\theta, \phi)] [{}_{2}Y_{lm}(\theta, \phi)] \sin^2\theta \propto \sum_{m=-l}^{l} [3m^2-l(l+1)] = 0~,
\end{align}
where $T_l^\zeta $ and $T_l^h$ denotes the contribution from $\zeta$ and $h$ to the temperature fluctuations, respectively. Thus no contribution from $\zeta h_+$ to $\langle TT\rangle$ could be observed once summing over $m$. On the other hand, before summing over $m$, the correlation does exist. But this does not help reconciling the tension between Planck and BICEP2 because the data analysis by the Planck team (for the purpose of tensor to scalar ratio) indeed sums over $m$. The above conclusion can actually be generalized to any primordial power spectrum, not restricted to isotropic inflation. The details of the general no-go result is presented in \cite{Emami:2014xga}.

The transfer function part encodes very complicated late time physics. Fortunately this part does not have angular dependence thus one can use the standard Boltzmann code for the calculation. Here ``the Cosmic Linear Anisotropy Solving System'' (CLASS) \cite{Blas:2011rf} is used for solving those transfer functions, where the cosmological parameters are chosen to be the same as those used by the BICEP2 group.

Now for the purpose of studying statistical anisotropies, we are not to sum over $m_1$ and $m_2$ (unlike the case above) because the summation would average away some signals of the anisotropy. Rather, we leave $m_1=m_2=m$ free and look into $m=0$ and $m=l$ cases respectively for illustration.

Three sets of parameters are examined numerically:
\begin{itemize}
\item Case I: Real Inflaton field with no electric gauge coupling:  In this class of models we plot $I=10^{-7}$ and $\e=0$. This case is the same as previous studies on the chargeless scalar field \cite{Watanabe:2010bu}. On the plots, Case I is shown in blue color.
\item Case II: Balanced: In this class of models we plot $I=10^{-7}$ and $\e=10^{-3}$. Such choice of parameters does not significantly modify $g_*$ of the scalar sector. However, the gravitational sector is largely modified because the two point correlation functions with tensors are more sensitive to the charge of the complex inflaton, as we have seen in the previous sections. On the plots, Case II is shown in black color.
\item Case III: Charge coupling dominated: In this class of models we plot $I=10^{-11}$ and $\e=0.025$. In this case  $g_*$, which is a measure of the scalar power spectrum anisotropy,
is considerably smaller than Cases I and II since $g_*$
is mostly sensitive to $I$ than $\e$. However, the charge contribution is large in the tensor sector and is marginally under control to calculate the statistical anisotropy perturbatively (where the anisotropic term is about 34\% of the isotropic term). On the plots, Case III is shown in red color.
\end{itemize}

By choosing the above parameters, we have taken into consideration that the anisotropies in the temperature correlations cannot be large, making use of Planck data in range $2 \leq l \leq 2000$ \cite{Kim:2013gka}. On the other hand, if scale dependence of the anisotropy is allowed, the constraint may not be so strong. In that case, we would have larger effects from anisotropies and our plots are shifted upwards.

The figures are plotted in logarithm scales because sometimes different cases differ by order-of-magnitude. However, we note that in various cases the correlations could go negative. To represent as detailed information as possible, we shall use solid lines to denote the logarithm of a quantity, where the quantity is positive (for example, $\log C_{l, l+2}^{TT}$) and dashed lines for multiplying $-1$ before taking logarithm (for example, $\log (- C_{l, l+2}^{TT})$), where the quantity is negative.

\subsection{The TT and BB modes}

With the above choice of parameters and conventions, here we plot the TT and BB power spectra in Figs.~\ref{fig:TT0}, \ref{fig:BB0}, \ref{fig:TT2} and \ref{fig:BB2}. In Figs.~\ref{fig:TT0} and \ref{fig:BB0} the $l_1=l_2$ part of the correlation function is plotted and in Figs.~\ref{fig:TT2} and \ref{fig:BB2} with $l_2 = l_1 +2$. The left and right panels corresponds to $m=0$ and $m=l_1$ respectively. From Figs.~\ref{fig:TT0}, \ref{fig:BB0}, one observes that the anisotropic modification to the standard TT and BB power spectra are slight but still visible.

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_TT_dL0}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_TT_dL0}
  \caption{\label{fig:TT0} The TT correlation at $l_2=l_1$. Here and hence after, the blue curve denotes Case I, with $I=10^{-7}$ and $\e=0$; the black curve denotes Case II, with $I=10^{-7}$ and $\e=10^{-3}$. The red curve denotes Case III, with $I=10^{-11}$ and $\e=0.025$. The left panel is for $m=0$ and the right panel is for $m=l_1$. It is important to note that the anisotropic contribution in Case III is considerably larger than the anisotropic contribution in case I and II at low $l$. But the contribution decays at high $l$. Here the black and blue curves are not visibly distinguishable from each other, which also coincides with the isotropic power spectrum, because of the tight constraint on the $g_*$ of the scalar sector.}
\end{figure}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_BB_dL0}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_BB_dL0}
  \caption{\label{fig:BB0} The $m=0$ (left) and $m=l_1$ (right) plots for BB correlation with $l_2=l_1$. }
\end{figure}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_TT_dL2}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_TT_dL2}
  \caption{\label{fig:TT2} The $m=0$ (left) and $m=l_1$ (right) plots for TT correlation with $l_2=l_1+2$. Here and hence after, the dashed lines denote the plotted quantity (here $C_{l, l+2}^{TT}$) is negative along this line segment, and thus we plot $-C_{l, l+2}^{TT}$ on the logarithm scales.}
\end{figure}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_BB_dL2}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_BB_dL2}
  \caption{\label{fig:BB2} The $m=0$ (left) and $m=l_1$ (right) plots for BB correlation with $l_2=l_1+2$.}
\end{figure}

\subsection{TB and EB correlations}

In the anisotropic case, the TB and EB correlations are opened up, with $l_2 = l_1 \pm 1$. In Fig.~\ref{fig:TBEB1}, those cross correlations are plotted.

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/mL_TB_dL1}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_EB_dL1}
  \caption{\label{fig:TBEB1} The $m=l_1$ plots for TB (left) and EB (right) correlation with $l_2=l_1+1$.}
\end{figure}

\subsection{TE and EE correlations}

In Figs. \ref{fig:TE0}, \ref{fig:EE0}, \ref{fig:TE2} and \ref{fig:EE2}, the TE and EE correlations for $l_2=l_1$ and $l_2=l_1+2$ are plotted respectively.

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_TE_dL0}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_TE_dL0}
  \caption{\label{fig:TE0} The $m=0$ (left) and $m=l_1$ (right) plots for TE correlation with $l_2=l_1$.}
\end{figure}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_EE_dL0}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_EE_dL0}
  \caption{\label{fig:EE0} The $m=0$ (left) and $m=l_1$ (right) plots for EE correlation with $l_2=l_1$.}
\end{figure}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_TE_dL2}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_TE_dL2}
  \caption{\label{fig:TE2} The $m=0$ (left) and $m=l_1$ (right) plots for TE correlation with $l_2=l_1+2$.}
\end{figure}

\begin{figure}[!h]
  \centering
  \includegraphics[width=0.45\textwidth]{fig/m0_EE_dL2}
  \hspace{0.05\textwidth}
  \includegraphics[width=0.45\textwidth]{fig/mL_EE_dL2}
  \caption{\label{fig:EE2} The $m=0$ (left) and $m=l_1$ (right) plots for EE correlation with $l_2=l_1+2$.}
\end{figure}

\section{Conclusion and discussion}

In this work we have studied  statistical anisotropies in model of anisotropic inflation from a charged inflaton field. More specifically, we work in large field inflation model with the chaotic potential  $V(\varphi) = \frac{1}{2}m^2 |\varphi|^2$ with the conformal coupling $f $ given
in Eq. (\ref{f-form0}).  It is worth to emphasis that our results, especially the enhancement of the tensor mode, is because we work in large field model. We have calculated the anisotropies generated in $\langle \zeta \zeta \rangle $ and $\langle h_{s}  h_{s} \rangle$ with the tensor polarizations $s=+, \times$. In addition, we also calculated the cross-correlation $\langle \zeta  h_{s} \rangle$ which is sourced in the  anisotropic inflationary background. Our general conclusion is that
the anisotropy in curvature perturbation power spectrum $\delta \calP_\zeta$ can be small so one can easily satisfy the observational bound  \cite{Kim:2013gka} $|g_*| \lesssim 10^{-2} $.  However, the effects of gauge coupling $\e$ appears strongly in tensor perturbations. In particular, we have shown that while $g_*$ is small enough to be within the observational bound, $\delta \calP_h$ induced from anisotropies  can be comparable to the background tensor power spectrum $\calP_h$.

The induced anisotropy of TT spectrum from the primordial tensor fluctuations can obtain an ``effective $g_*$''  as large as a few percents at low $l$. This effective $g_*$ contribution is decaying towards high $l$, characterized by the tensor to temperature transfer function. It remains interesting to explore CMB anisotropies with this new profile of TT anisotropy.

We would also like to mention that the mechanism we have used to generate small $l$ CMB anomalies from the tensor sector is quite general. It may also be applied to other kinds of CMB anomalies such as non-Gaussianities, hemispherical asymmetry, etc. In particular, the possibility that
hemispherical asymmetry in tensor modes which are generated from the long mode modulations  \cite{Dai:2013kfa, Abolhasani:2013vaa} may be behind the tension between the BICEP2 and
Planck has been emphasized in \cite{Chluba:2014uba}.
It is also another interesting way to get out of this discrepancy as it has been shown by \cite{Kamionkowski:2014faa}.
We hope we can address these issues in a future work.

Note that we have assumed that the anisotropies are small and have expanded those anisotropies at linear order. This is a very good assumption for models with no electric charge since both the scalar sector and the tensor sector remains nearly isotropic as studied in \cite{Ohashi:2013qba}.
However, in the presence of the  gauge coupling, things become different. As described above
with large enough value of $\e$ the tensor sector becomes completely anisotropic. Thus our assumption of small anisotropy is now a purely technical simplification. It is interesting to generalize our perturbative technique in analyzing the tensor perturbations   into a non-perturbative calculation, in which case larger cross-correlations with more significant observational signatures can be calculated.

 Based on the correlations $\langle \zeta \zeta \rangle $, $\langle h_{s}  h_{s} \rangle$
 and $\langle \zeta h_s \rangle $ obtained in this model, we have calculated the observational signatures on the CMB temperature and polarization maps. Before the BICEP2 detection of B-modes polarization, the probe of statistical anisotropies in the cosmological perturbations were confined to the TT, TE and EE maps, which are dominated by the anisotropy of the scalar sector. However, this statement is now changed after  the recent detection of the B-modes. The TB, EB and BB correlations are now available for testing the anisotropies. Especially, we have shown that when the gauge interaction dominates over  the interaction between the inflaton and the gauge field, the tensor sector can be more anisotropic than the scalar sector, resulting in enhanced correlations involving B-modes. For example, in Fig~\ref{fig:TBEB1} one can observe that the TB correlation in the model with the
 charged coupling can be order-of-magnitude larger compared to model with no electric charge coupling,  which is also considerably larger than the size of the BB correlation near the reionization bump at $l\lesssim 10$.

Having this said note that the sensitivities on different CMB correlations we studied here also crucially depend on the design of experiments and the proprieties of the detectors. It would be very interesting to perform experiment-specific investigations and dig into data to hunt for the anisotropic correlations.

\section*{Acknowledgment}
It is a pleasure to thank A. A. Abolhasani, S. Baghram,   N. Bartolo,  J. Chluba,  D. Jeong,  E. Komatsu,  M. Peloso  and  S.-C.~Su for helpful discussion. YW is supported by a Starting Grant of the European Research Council (ERC STG grant 279617), and the Stephen Hawking Advanced Fellowship.

redpajama_set_name: RedPajamaArXiv