On the primordial correlation of gravitons with
gauge fields

Rajeev Kumar<EMAIL_ADDRESS>P. Jishnu<EMAIL_ADDRESS>Martin S<EMAIL_ADDRESS>
aDepartment of Physics, Indian Institute of Science,
Bangalore 560012, India
bCP3-Origins, Centre for Cosmology and Particle Physics Phenomenology,
University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark

###### Abstract

We calculate the primordial correlation of gravitons with an abelian gauge
field non-minimally coupled through a dynamical dilaton field or a volume
moduli during inflation in the early universe. In particular, we compute the
cross-correlation of a tensor mode with two gauge field modes and the
corresponding correlation functions for the associated magnetic and electric
fields using the in-in formalism. Moreover, using semi-classical methods, we
show that the three-point cross-correlation functions satisfy new consistency
relations (soft theorems) in the squeezed limit. Our findings exhibit a
complete agreement of the full in-in results with the new consistency
relations. An interesting consequence of our scenario is the possibility of a
novel correlation of the primordial tensor mode with the primordial curvature
perturbation. Finally, we discuss how these correlation functions are
imprinted on cosmological observables today and the applications to scenarios
of inflationary magnetogenesis.

###### Contents

  1. 1 Introduction
  2. 2 Dynamics of primordial tensors and gauge fields during inflation
  3. 3 Cross-correlation of inflationary tensor perturbation with primordial gauge fields
    1. 3.1 A complete calculation using in-in formalism
      1. 3.1.1 $\langle\gamma A_{\mu}A^{\mu}\rangle$
      2. 3.1.2 $\langle\gamma B_{\mu}B^{\mu}\rangle$
      3. 3.1.3 $\langle\gamma E_{\mu}E^{\mu}\rangle$
    2. 3.2 The magnetic and electric non-linearity parameters
    3. 3.3 The squeezed limit
  4. 4 Semi-classical derivation of the consistency relations
  5. 5 A novel correlation of tensor and curvature perturbations
  6. 6 Conclusions and discussions
  7. A Useful integrals
  8. B Validity of the semi-classical approach
  9. C The in-in formalism and the dynamical correction term

##  1 Introduction

In effective field theories with extra dimensions, as derived from string
theory, it is natural to expect that gauge fields are non-minimally coupled in
the early universe through a dynamical dilaton field or the moduli of the
internal dimensions. This breaks the conformal invariance of the gauge field
action and enables the amplification of the gauge fields during inflation. If
gauge fields are enhanced in the early universe, during the inflationary
epoch, it opens the possibility that such a primordial gauge field could have
left a detectable imprint on the observable universe today. It is well known
that observable primordial gravitational waves may also be created during
inflation by the same causal mechanism, and we therefore expect a non-
vanishing non-Gaussian cross-correlation between primordial gravitational
waves and primordial gauge fields, which we will calculate here. Focussing for
simplicity on abelian gauge fields, we will consider a model where conformal
invariance is broken by a non-minimal coupling between an evolving scalar
field and the kinetic term of the gauge field [1, 2], which takes the form

$\mathcal{L}_{\textrm{gauge}}=\lambda(\phi)F_{\mu\nu}F^{\mu\nu}~{}.$ (1)

Such models are often referred to as the kinetic coupling models. There could
also be a non-minimal coupling with the parity violating term of the gauge
field. However, we shall not consider those models in this paper. In all these
scenarios of broken conformal symmetry, one also assumes that the conformal
invariance is restored by the end of inflation (or by the end of the reheating
epoch). Moreover, such a model only breaks the conformal invariance, and not
the gauge invariance.

A further motivation for studying this type of models stems from the fact that
astrophysical observations indicate the presence of coherent magnetic fields
of femto-Gauss strength on large cosmological scales (Mpc scales or larger)
exceeding the scales of galaxies and galaxy clusters [3, 4, 5]. It is unclear
if any astrophysical process can produce magnetic fields with coherence
lengths of cosmological sizes, and it has therefore been speculated that such
fields must have an inflationary origin. After all, inflation is known to be
the perfect mechanism for producing the very large, even super-horizon,
correlations of primordial fluctuations, needed to explain the observed large
scale correlations in the cosmic microwave background (CMB) anisotropies.

For inflationary magnetogenesis to work, it requires a breaking of the
conformal invariance of electromagnetism during inflation, and therefore imply
a departure from standard electromagnetism as in eq. (1)[2, 1, 6]. But even
then, back-reaction arising due to the generated electromagnetic fields and
issues with entering regimes of strong coupling severely restricts the
possibility for significant magnetic fields to be generated during inflation
[7]. Only in some very special cases, it is possible for the observed magnetic
fields to have been generated during inflation [8, 9, 10, 11, 12, 13]. It
would therefore be very important to understand if magnetic fields observed on
cosmological scales really have an inflationary origin or not.

One possibility is to look at the correlation of the cosmological magnetic
fields with the inflaton perturbation, or equivalently the gauge-invariant
curvature perturbation. Such correlations have earlier been discussed in [14,
15, 16, 17, 18], and also in [19, 20] in a different set-up. The idea being
that magnetic fields, if they are indeed generated during inflation, will be
correlated with the primordial curvature perturbation. On the other hand, if
they are generated after inflation, no such correlation will exist. Assuming
an exact model for generating the magnetic fields during inflation, one can
compute the expected correlation. Since one of the few models, which can
indeed generate significant primordial magnetic fields during inflation, takes
the general form in eq. (1), much attention was focussed on this model in the
previous work, and in [16, 17] the non-Gaussian correlation function of the
curvature perturbation with the magnetic field

$\left<\zeta(\bf{k_{1}})B(\bf{k_{2}})\cdot B(\bf{k_{3}})\right>~{},$ (2)

was calculated. It was shown that such a non-Gaussian correlation function
satisfies a new consistency relation in the squeezed limit and its magnitude
becomes quite large in the flattened configuration in the Fourier space.

Here we will consider the natural extension of that work, and look at the non-
trivial cross correlations of primordial gauge fields and primordial
gravitational waves. Our main objective is therefore to compute the primordial
non-Gaussian correlations functions

$\left<\gamma(\bf{k_{1}})A(\bf{k_{2}})\cdot
A(\bf{k_{3}})\right>\,,\qquad\left<\gamma(\bf{k_{1}})B(\bf{k_{2}})\cdot
B(\bf{k_{3}})\right>\,,\qquad\left<\gamma(\bf{k_{1}})E(\bf{k_{2}})\cdot
E(\bf{k_{3}})\right>~{},$ (3)

where $\gamma$ is the tensor mode, ${\bf A}$ is our abelian gauge boson
(vector potential) and ${\bf B}$, ${\bf E}$ are the associated magnetic and
electric fields, respectively. Since these are equal time correlation
functions, we will calculate them using the full in-in formalism [21].

As a way of checking our final results, we propose new semi-classical
consistency relations (soft theorems444Since the semi-classical relations
relate a higher order correlation function involving an infrared (soft mode)
to a lower order correlation function, the semi-classical relations are also
sometimes called the soft theorems with a reference to the Weinberg’s soft
photon and graviton theorems [22], which can be related to the asymptotic
symmetries and the gravitational memory of flat space [23, 24]. For a
discussion of the relation between semi-classical relations/soft theorems,
asymptotic symmetries and gravitational memory in inflation, see [25, 26, 27]
and also [28, 28, 29, 30, 31, 32].) in the limit wherein the wavelength of the
graviton goes to infinity (the squeezed limit). The semi-classical consistency
relation for correlation functions only involving the curvature perturbation,
$\zeta$, or the tensor mode, $\gamma$, was first proposed by Maldacena in [33,
34] (and also subsequently by others [35, 36]), hence they are also sometimes
called the Maldacena consistency relations. A related consistency relation for
gauge fields was first discussed in [16], wherein a magnetic consistency
relation for the correlation function $\left<\zeta\,{\bf B}\cdot{\bf
B}\right>$ was proposed. Here we will consider the natural extension to the
correlation functions above in eq. (3), and show that the full in-in
correlation function in fact satisfies the semi-classical consistency relation
in the appropriate limits555The semi-classical relations serve an important
role as a cross check of the correctness of the full in-in correlation
functions (not only Maldacena’s three-point function). As was first shown in
[37], it can also be used to calculate the contribution to the four-point
function in the limit where the momentum of the exchanged graviton goes to
zero, the so called counter collinear limit, by simply taking the square of
the three-point function in the squeezed limit. In [38], it was further
realized that one can also extract the infrared divergent part of higher order
loop diagrams using such semi-classical relations..

One of the problems with generating large magnetic fields during inflation, is
that the universe becomes an almost perfect conductor after reheating, which
very rapidly erases any magnetic field which would have been generated during
inflation. Models of successful magnetogenesis therefore typically require a
very low scale of inflation [11, 12], or at least a very low reheating
temperature [13]. This problem is absent if the magnetic field belongs to a
dark gauge group. Therefore we can more easily imagine that large dark
magnetic fields could have been generated during inflation. While such dark
magnetic fields are not directly observable in the same way as the standard
magnetic fields of electromagnetism, they do leave indirect imprints in the
curvature perturbation of the universe. The total curvature perturbations is
then the sum of the ”usual” contribution, $\zeta_{0}$ generated independently
of the magnetic field, and the contribution, $\zeta_{B}$, imprinted by
magnetic fields

$\zeta(\tau)\simeq\zeta_{0}(\tau)+\zeta_{B}(\tau)~{},$ (4)

where the curvature perturbation induced by any magnetic field is [39]

$\zeta_{B}(\tau)=\int_{\tau_{0}}^{\tau}d\ln\tilde{\tau}\lambda(\tilde{\tau})\frac{B_{i}B^{i}}{3H^{2}\epsilon}~{}.$
(5)

This expression is valid for both dark and ordinary magnetic fields. In the
remainder of this paper, it can refer to both cases. The magnetic fields
generated during inflation follow a Gaussian statistics to the leading order
in perturbations which implies that the induced curvature perturbation
$\zeta_{B}$ is generally a non-Gaussian field. Various imprints of this
induced component such as non-Gaussianities and anisotropies in the power
spectrum and bispectrum have also been greatly studied [40, 41, 42, 43, 44,
45, 46, 47].

One very interesting consequence of such an induced contribution is that this
enables the possibility of a direct non-vanishing correlation of the
primordial tensor mode with the primordial curvature perturbation of the
symbolic form (we will be more precise later)

$\left<\gamma\zeta\right>\simeq\left<\gamma\zeta_{B}\right>\propto\left<\gamma{\bf
B\cdot B}\right>\neq 0~{}.$ (6)

A similar contribution will arise from the corresponding electric fields as
well. We will therefore propose the two-point correlation function of
primordial gravitational waves with the primordial curvature perturbation, as
a novel observable, which can possibly be used to probe the existence of
primordial (dark) gauge fields666Some scaling arguments of such a mixed
correlator based on the special conformal transformations have been discussed
in [48] when there is an explicit parity violating term in the Lagrangian..

This paper is organized as follows. In the next section we will quickly review
the dynamics of tensor perturbations during inflation and the kinetic coupling
model as a mechanism for the production of large scale primordial gauge fields
during inflation. In section 3, We calculate the full bispectrum associated
with the cross correlations of the inflationary tensor perturbation with gauge
fields using the very general in-in formalism and discuss the extent of non-
linearities using the Fourier space shape functions. In section 4, we derive
the consistency relations for these non-Gaussian cross correlations using a
semi-classical approach. In section 5, we calculate a novel non-vanishing
correlation of the primordial tensor mode and the primordial curvature
perturbation induced by the gauge fields. Finally, in section 6, we summarise
our results and conclude with a discussion. In appendix A, we have listed some
useful integrals appearing in the in-in results of various cross correlations.
In appendix B, we discuss the validity of the semi-classical approach used to
derive our consistency relations and in appendix C, we discuss the presence of
the dynamical correction terms arising in the cross correlation involving
electric fields in the in-in formalism.

Throughout this paper, we work in natural units with $\hbar=c=1$, and the
Planck mass $M_{\rm Pl}^{2}=1/8\pi G$ is set to unity. Our metric convention
is $(-,+,+,+)$.

##  2 Dynamics of primordial tensors and gauge fields during inflation

In this section, we shall briefly discuss the dynamical evolution of
primordial tensor modes and abelian gauge fields during inflation. The
homogeneous and isotropic background during the inflationary expansion is
described by the spatially flat FLRW metric

$\displaystyle ds^{2}=-dt^{2}+a^{2}(t)\,d{\bf
x}^{2}=a^{2}(\tau)\left(-d\tau^{2}+d{\bf x}^{2}\right),$ (7)

where $\tau$ is the conformal time, defined by $d\tau=dt/a$ and $a(\tau)$ is
the scale factor. The perturbed FLRW metric in the presence of tensor modes is
given by

$\displaystyle
ds^{2}=-dt^{2}+a^{2}(t)\bigl{[}e^{\gamma}\bigr{]}_{ij}dx^{i}dx^{j}\approx-
dt^{2}+a^{2}(t)\bigl{[}\delta_{ij}+\gamma_{ij}\bigr{]}dx^{i}dx^{j}~{},$ (8)

where $\gamma_{ij}$ is the metric tensor perturbation which is transverse and
traceless i.e. $\gamma_{{ij},j}=0$ and $\gamma_{ii}=0$. These conditions lead
to only two independent radiative degrees of freedom of tensor perturbations
which correspond to the two polarisations of gravitational waves (GW). To the
linear order in metric fluctuations, there is typically no active source of
GWs during inflation. However, quantum fluctuations of $\gamma_{ij}$ are
parametrically amplified by the quasi-exponential expansion of the universe.
In order to describe this phenomenon, one needs to quantise the canonical
degrees of freedom associated with the tensor perturbations. Following the
usual quantisation formalism, the corresponding mode expansion of tensor
perturbation is defined as

$\displaystyle\gamma_{ij}({\bf x},\tau)=\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\sum_{s=\pm 2}\left[\gamma_{k}(\tau)\,e^{i{\bf k}\cdot{\bf
x}}\,\epsilon^{s}_{ij}(\hat{\bf k})\,b_{\bf k}^{s}+h.c.\right],$ (9)

where $\epsilon_{ij}^{s}$ is the polarisation tensor corresponding to the
helicity $s$, with the normalisation
$\epsilon_{ij}^{s}\epsilon_{ij}^{*s^{\prime}}=2\delta_{ss^{\prime}}$ and the
creation and annihilation operators do satisfy the condition, $[b_{\bf
k}^{s},b^{s^{\prime}\dagger}_{\bf k^{\prime}}]=(2\pi)^{3}\delta^{(3)}({\bf
k}-{\bf k^{\prime}})\delta_{ss^{\prime}}$. The tensor helicity mode
$\gamma^{s}({\bf k},\tau)$ can be defined through the following equation,

$\displaystyle\gamma_{ij}({\bf x},\tau)=\sum_{s=\pm 2}\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\left[\epsilon^{s}_{ij}(\hat{\bf k})\gamma^{s}({\bf
k},\tau)\,e^{i{\bf k}\cdot{\bf x}}\right]\,.$ (10)

We can now define the two point function in Fourier space as

$\displaystyle\langle\gamma^{s}({\bf k},\tau)\gamma^{s^{\prime}}({\bf
k^{\prime}},\tau)\rangle=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}+\bm{\mathrm{{k^{\prime}}}})P_{\gamma}(k,\tau)\delta_{ss^{\prime}}$
(11)

where the power spectrum of tensor modes is given by
$P_{\gamma}(k,\tau)=|\gamma_{k}(\tau)|^{2}$. It is often useful to define a
dimensionless tensor power spectrum ${\cal P}_{\gamma}$ as

${\cal
P}_{\gamma}(k,\tau)=\frac{k^{3}}{2\pi^{2}}P_{\gamma}(k,\tau)=\frac{k^{3}}{2\pi^{2}}|\gamma_{k}(\tau)|^{2},$
(12)

which is usually evaluated in the super-horizon limit $|k\tau|\to 0$. For a
nearly de-Sitter background $\tau\simeq-1/aH$ where $H$ is the Hubble
parameter, the Bunch-Davies normalised solution for the mode function
$\gamma_{k}(\tau)$ is given by

$\gamma_{k}(\tau)=\frac{H}{\sqrt{k^{3}}}(1+ik\tau)\,e^{-ik\tau}~{},$ (13)

which, in the super-horizon limit, leads to ${\cal
P}_{\gamma}\sim(H/2\pi)^{2}$. The tensor power spectrum thus measures the
energy scale during inflation and is exactly scale invariant in a de-Sitter
background (and nearly scale invariant in a slow roll inflation).

As mentioned already in the introduction, it is well known that gauge fields
described by the canonical Maxwell’s action are not amplified in the FLRW
background due to the conformal symmetry and thus, a necessary condition for
amplification of gauge fields during inflation is to break this conformal
invariance. One of the simplest possibility which has been explored to a large
extent in the literature is the so called kinetic coupling model which breaks
the conformal invariance (while preserving gauge invariance) by explicitly
coupling gauge fields to the inflaton as [2]

$S_{A}=-\frac{1}{4}\int
d^{4}x\sqrt{-g}\,\lambda(\phi)F_{\mu\nu}F^{\mu\nu}+\int
d^{4}x\sqrt{-g}\,j^{\mu}A_{\mu},$ (14)

where the gauge field tensor is defined as
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and the second term
represents the interaction where $j^{\mu}$ is the four current density. For
simplicity, we shall neglect the interaction term and assume that there are no
free charges. In the Coulomb gauge wherein $A_{0}=0$ and
$\partial_{i}A^{i}=0$, the quadratic action for the gauge field $A_{i}$
becomes

$S_{A}=\frac{1}{2}\int
d^{3}x\,d\tau\lambda(\phi)\left({A_{i}^{\prime}}^{2}-\frac{1}{2}(\partial_{i}A_{j}-\partial_{j}A_{i})^{2}\right).$
(15)

Since we are interested in the dynamical amplification of the quantum gauge
field originating from the vacuum, we promote the classical gauge field
$A_{i}$ to an operator and impose the canonical quantization conditions.
Following this, we define the usual mode expansion as

$A_{i}({\bf x},\tau)=\int\frac{d^{3}{\bf k}}{(2\pi)^{3}}\sum_{\lambda=\pm
1}\left[A_{k}(\tau)\,e^{i{\bf k}\cdot{\bf
x}}\,\epsilon^{(\lambda)}_{i}(\hat{\bf{k}})\,a_{\bf k}^{\lambda}+h.c.\right],$
(16)

and impose the standard commutation relations

$[a_{\bf k}^{\lambda},a^{\lambda^{\prime}\dagger}_{\bf
k^{\prime}}]=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}-\bm{\mathrm{{k}}}^{\prime})\delta_{\lambda\lambda^{\prime}}~{}.$
(17)

Here $A_{k}(\tau)$ is the mode function in Fourier space and the polarization
vectors will satisfy the usual relations:
$\bm{\mathrm{{k}}}\cdot\epsilon^{\sigma}(\hat{k})=0$,
$\epsilon^{\sigma}(\hat{\bf k})\cdot{\epsilon^{\sigma^{\prime}}}^{*}(\hat{\bf
k})=\delta_{\sigma\sigma^{\prime}}$, and
$\sum_{\sigma=\pm}\epsilon_{i}^{\sigma}(\hat{\bf
k}){\epsilon_{j}^{\sigma}}^{*}(\hat{\bf k})=\delta_{ij}-k_{i}k_{j}/k^{2}$.

One can now define a canonically normalized gauge field $v_{i}=S(\tau)A_{i}$
where $S^{2}(\tau)=\lambda(\phi(\tau))$. In this case, the quadratic action
for the gauge field takes the canonical form as

$S_{A}=\frac{1}{2}\int d^{3}x\,d\tau\left[v_{i}^{\prime
2}-(\partial_{j}v_{i})^{2}+\frac{S^{\prime\prime}}{S}v_{i}^{2}\right]~{},$
(18)

and the equation of motion for the rescaled mode function,
$v_{k}=S(\tau)A_{k}$, takes the form of a harmonic oscillator with a time
dependent mass term as

$v_{k}^{\prime\prime}+\left(k^{2}-\frac{S^{\prime\prime}}{S}\right)v_{k}=0~{}.$
(19)

In order to compute the two point and higher order correlators, one needs to
define the coupling function $\lambda(\phi)$. For instance, the dilatonic
coupling corresponds to $\lambda(\phi)\sim{\rm exp}(2\phi)$ [2, 49]. Since
$\phi$ must be dynamical, it turns out to be more convenient to parametrize
the time dependence of the coupling function as
$\lambda(\phi(\tau))=\lambda_{I}(\tau/\tau_{I})^{-2n}$ where $\tau_{I}$ is the
conformal time at the end of inflation, and $\lambda_{I}$ is the coupling
evaluated at the end of inflation $\tau=\tau_{I}$.

For a power law coupling $\lambda(\tau)=\lambda_{I}(\tau/\tau_{I})^{-2n}$, one
finds that $S^{\prime\prime}/S=(\nu^{2}-1/4)/\tau^{2}$ with $\nu=(n+1/2)$ and
the mode solution can be written in terms of the Hankel functions. In the sub-
horizon limit, the mode solution must be normalized to the Bunch-Davies vacuum

$v_{k}(\tau)=\frac{1}{\sqrt{2k}}e^{-ik\tau}~{},$ (20)

which leads to the full solution for $v_{k}$ as

$\displaystyle
v_{k}(\tau)=\frac{\sqrt{\pi}}{2}e^{i\pi(n+1)/2}\sqrt{-\tau}H^{(1)}_{n+\frac{1}{2}}(-k\tau)~{}.$
(21)

With this solution, the mode function $A_{k}(\tau)$ is expressed as

$\displaystyle
A_{k}(\tau)=\frac{1}{\sqrt{\lambda_{I}}}\frac{\sqrt{\pi}}{2}e^{i\pi(n+1)/2}\sqrt{-\tau}\left(\frac{\tau}{\tau_{I}}\right)^{n}H^{(1)}_{n+\frac{1}{2}}(-k\tau)~{},$
(22)

where $H_{n}^{(1)}(x)$ is a Hankel function of the first kind. The two point
correlation function of gauge fields $A_{i}$ can be obtained as follows,

$\left<A_{i}(\tau,\bm{\mathrm{{k}}})A_{j}(\tau,\bm{\mathrm{{k}}}^{\prime})\right>=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}+\bm{\mathrm{{k}}}^{\prime})\left(\delta_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right)\left|A_{k}(\tau)\right|^{2}.$
(23)

With respect to an observer characterised by the 4-velocity vector $u^{\mu}$,
one can covariantly define the electric field $E_{\mu}$ and magnetic field
$B_{\mu}$ as [50].

$E_{\mu}=F_{\mu\nu}u^{\nu},\qquad B_{\mu}={{}^{*}F}_{\mu\nu}u^{\nu}\,,$ (24)

where ${{}^{*}F}_{\mu\nu}$ is the dual of the electromagnetic field tensor
$F_{\mu\nu}$ which is defined by the relation

${{}^{*}F}_{\mu\nu}=\frac{1}{2}\eta_{\mu\nu\alpha\beta}F^{\alpha\beta},$ (25)

with $\eta_{\mu\nu\alpha\beta}$ being the totally antisymmetric permutation
tensor of space-time with $\eta_{0123}=\sqrt{-g}$. In this work, the covariant
magnetic and electric fields are defined with respect to a comoving observer
having 4-velocity $u^{\mu}=(\frac{1}{a},0,0,0)$. The power spectra of the
magnetic field and electric field are defined as

$\displaystyle\left<B_{\mu}(\tau,\bm{\mathrm{{k}}})B^{\mu}(\tau,\bm{\mathrm{{k}}}^{\prime})\right>$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}+\bm{\mathrm{{k}}}^{\prime})P_{B}(k,\tau),$
(26)
$\displaystyle\left<E_{\mu}(\tau,\bm{\mathrm{{k}}})E^{\mu}(\tau,\bm{\mathrm{{k}}}^{\prime})\right>$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}+\bm{\mathrm{{k}}}^{\prime})P_{E}(k,\tau).$
(27)

Note that, both $P_{B}$ and $P_{E}$ are dimension-full. We can now easily
express $P_{B}$ and $P_{E}$ in terms of the mode function $A_{k}(\tau)$ as

$\displaystyle P_{B}(k,\tau)$ $\displaystyle=$ $\displaystyle
2\frac{k^{2}}{a^{4}}|A_{k}(\tau)|^{2}$ (28) $\displaystyle P_{E}(k,\tau)$
$\displaystyle=$ $\displaystyle\frac{2}{a^{4}}|A^{\prime}_{k}(\tau)|^{2}$ (29)

which can be readily calculated using the obtained mode solution (22) for the
case of a power law coupling [49, 50, 14, 17]. For the cases of $n=2$ and
$n=-3$, the magnetic field spectrum turns out to be scale invariant i.e.
$P_{B}\propto k^{-3}$ while the electric field spectrum behave as
$P_{E}\propto k^{-1}$ and $P_{E}\propto k^{-5}$, respectively. In both these
cases, there arise various issues such as the strong coupling and the back
reaction problem which have been discussed to a great extent in the literature
[49, 7, 50, 11, 12].

##  3 Cross-correlation of inflationary tensor perturbation with primordial
gauge fields

In this section, we shall discuss our calculations of the cosmological cross-
correlations of primordial tensor perturbation with gauge fields using the in-
in formalism while in the next section, we shall discuss the semiclassical
approach to arrive at these correlators. As expected, both these approaches
lead to the same results for these cross-correlations in the squeezed limit.

###  3.1 A complete calculation using in-in formalism

In order to compute the correlation function during inflation, we adopt a very
powerful tool of the in-in formalism which can be used to calculate various
correlation functions at equal times. A detailed discusses of the in-in
formalism can be found in appendix C. In this formalism, the expectation value
of an operator $\mathcal{O}$ at time $\tau_{I}$ is given by

$\left\langle\mathcal{O}(\tau_{I})\right>=\left\langle
0\right|{\bar{\mathrm{T}}}\left(e^{i\int_{-\infty}^{\tau_{I}}d\tau H_{\rm
int}}\right)\mathcal{O}(\tau_{I})\mathrm{T}\left(e^{-i\int_{-\infty}^{\tau_{I}}d\tau
H_{\rm int}}\right)\left|0\right>$ (30)

where $H_{\rm int}$ is the interaction Hamiltonian and $\mathrm{T}$,
$\bar{\mathrm{T}}$ are the time ordering and anti-time ordering operators,
respectively. The quantum state $\left|0\right>$ refers to the physical vacuum
state of the theory which we take to be the standard Bunch-Davies vacuum
state. We shall first apply the in-in formalism to calculate the cross-
correlations of the tensor mode with gauge fields.

####  3.1.1 $\langle\gamma A_{\mu}A^{\mu}\rangle$

Let us first compute the cross-correlation of a tensor mode with two gauge
field modes, i.e., a correlator of the form $\langle\gamma
A_{\mu}A^{\mu}\rangle$. We first find that

$\displaystyle\gamma A_{\mu}A^{\mu}=\frac{1}{a^{2}}\gamma
A_{i}A_{i}-\frac{1}{a^{2}}\gamma\gamma_{ij}A_{i}A_{j}+{\cal O}(\gamma^{3})$
(31)

where all the repeated indices in the above and following expressions are
summed over (irrespective of up or down). By using the master formula of the
in-in formalism (30), we can calculate $\langle\gamma A_{\mu}A^{\mu}\rangle$
perturbatively which, to the leading order, is given by

$\displaystyle\langle\gamma({\bf k}_{1},\tau_{I})A_{\mu}({\bf
k}_{2},\tau_{I})A^{\mu}({\bf k}_{3},\tau_{I})\rangle$ $\displaystyle=$
$\displaystyle\left\langle 0\right|\gamma({\bf k}_{1},\tau_{I})A_{\mu}({\bf
k}_{2},\tau_{I})A^{\mu}({\bf k}_{3},\tau_{I})\left|0\right>$ (32)
$\displaystyle+$ $\displaystyle i\int_{-\infty}^{\tau_{I}}d\tau\left\langle
0\right|[H_{\rm int}(\tau),\gamma({\bf k}_{1},\tau_{I})A_{\mu}({\bf
k}_{2},\tau_{I})A^{\mu}({\bf k}_{3},\tau_{I})]\left|0\right>,$

where $H_{\rm int}$ is the cubic order interaction Hamiltonian which can be
calculated from equation (14),

$H_{\rm int}(\tau)=\frac{1}{2}\int
d^{3}x\,\lambda(\tau)\left(\gamma^{ij}A^{\prime}_{i}A^{\prime}_{j}-\gamma^{ij}\delta^{kl}(\partial_{i}A_{k}\partial_{j}A_{l}+\partial_{k}A_{i}\partial_{l}A_{j})+2\gamma^{ij}\delta^{kl}\partial_{i}A_{k}\partial_{l}A_{j}\right).$
(33)

Figure 1: The leading order cosmological diagram contributing to all equal
time $\tau_{I}$ cross-correlations of our interest. In addition to the
contribution from this diagram, there will also be additional corrections for
each correlator, as discussed in the text. For $\left\langle\gamma
A_{\mu}A^{\mu}\right>$ and $\left\langle\gamma B_{\mu}B^{\mu}\right>$, there
is a kinematical correction term. But for the case of $\left\langle\gamma
E_{\mu}E^{\mu}\right>$, there is not only a kinematical correction term but
also a dynamical correction term.

By substituting (31) in (32), we observe that the first term on the right hand
side (RHS) of equation (32) gives a non-zero contribution. It can be easily
calculated from the Fourier mode functions (13) and (22) as,

$\displaystyle\left\langle 0\right|\gamma({\bf k}_{1},\tau_{I})A_{\mu}({\bf
k}_{2},\tau_{I})A^{\mu}({\bf
k}_{3},\tau_{I})\left|0\right>=\frac{1}{2}(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\epsilon_{ij}\frac{k_{2i}k_{2j}}{\tilde{k}_{2}^{2}}P_{\gamma}(k_{1},\tau_{I})P_{A}(\tilde{k}_{2},\tau_{I}).$
(34)

Here, ${\bf\tilde{k}}_{2}={\bf k_{2}}+\frac{1}{2}{\bf k_{1}}$,
$P_{\gamma}(k,\tau)$ is the power spectrum of the tensor modes, as defined in
(11) and $P_{A}(k,\tau)$ is the power spectrum of the gauge fields defined as,

$\displaystyle\langle A_{\mu}({\bf k},\tau)A^{\mu}({\bf
k^{\prime}},\tau)\rangle=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}+\bm{\mathrm{{k}}}^{\prime})P_{A}(k,\tau).$
(35)

Thus, the vacuum expectation value of $\gamma A_{\mu}A^{\mu}$ does contribute
to this non-Gaussian correlator at the leading order. Since this contribution
is originating purely from the kinematical considerations without using the
dynamics of gauge fields, we refer to such contributions as the kinematical
correction terms. As we shall see later, similar kinematical terms will also
appear for the cases of $\langle\gamma B_{\mu}B^{\mu}\rangle$ and
$\langle\gamma E_{\mu}E^{\mu}\rangle$, respectively.

Now, the second term in the RHS of equation (32) is calculated using the
cosmological diagrammatic rules [21]. The leading order diagram contributing
to this correlator is shown in figure 1. Substituting the interaction
Hamiltonian $H_{\rm int}$ from equation (33) in equation (32) and after some
simplifications, we obtain the total contribution as,

$\displaystyle\langle\gamma({\bf k}_{1})A_{\mu}({\bf k}_{2})A^{\mu}({\bf
k}_{3})\rangle$ $\displaystyle=$ $\displaystyle(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\Biggl{[}\frac{1}{2}\frac{k_{2i}k_{2j}}{\tilde{k}_{2}^{2}}P_{\gamma}(k_{1})P_{A}(\tilde{k}_{2})$
(36) $\displaystyle+$
$\displaystyle\left(\delta_{mi}-\frac{k_{2m}k_{2i}}{k_{2}^{2}}\right)\left(\delta_{mj}-\frac{k_{3m}k_{3j}}{k_{3}^{2}}\right)\left(\frac{\mathcal{I}_{1}}{a^{2}}+({\bf
k_{2}\cdot k_{3}})\frac{\mathcal{I}_{2}}{a^{2}}\right)$ $\displaystyle+$
$\displaystyle k_{2i}k_{3j}\left(1+\frac{({\bf k_{2}\cdot
k_{3}})^{2}}{k_{2}^{2}k_{3}^{2}}\right)\frac{\mathcal{I}_{2}}{a^{2}}-k_{3i}k_{2n}\left(\delta_{mn}-\frac{k_{3m}k_{3n}}{k_{3}^{2}}\right)\left(\delta_{mj}-\frac{k_{2m}k_{2j}}{k_{2}^{2}}\right)\frac{\mathcal{I}_{2}}{a^{2}}$
$\displaystyle-$ $\displaystyle
k_{2j}k_{3n}\left(\delta_{mn}-\frac{k_{2m}k_{2n}}{k_{2}^{2}}\right)\left(\delta_{mi}-\frac{k_{3m}k_{3i}}{k_{3}^{2}}\right)\frac{\mathcal{I}_{2}}{a^{2}}\Biggr{]},$

where the integrals $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ can be expressed
in terms of the mode functions of tensor mode and gauge field and their time
derivatives as

$\displaystyle\mathcal{I}_{1}$ $\displaystyle=$ $\displaystyle 2\,{\rm
Im}\left[\gamma_{k_{1}}(\tau_{I})A_{k_{2}}(\tau_{I})A_{k_{3}}(\tau_{I})\int
d\tau\lambda(\tau)\gamma_{k_{1}}^{*}(\tau)A_{k_{2}}^{{}^{\prime}*}(\tau)A_{k_{3}}^{{}^{\prime}*}(\tau)\right],$
(37) $\displaystyle\mathcal{I}_{2}$ $\displaystyle=$ $\displaystyle 2\,{\rm
Im}\left[\gamma_{k_{1}}(\tau_{I})A_{k_{2}}(\tau_{I})A_{k_{3}}(\tau_{I})\int
d\tau\lambda(\tau)\gamma_{k_{1}}^{*}(\tau)A_{k_{2}}^{*}(\tau)A_{k_{3}}^{*}(\tau)\right],$
(38)

and the explicit time dependence, $\tau_{I}$, of the operators and the power
spectra in the equation (36) is omitted for notational convenience. It is
useful to rewrite these integrals in (37) and (38) in terms of another set of
auxiliary integrals, denoted by $\tilde{\mathcal{I}}_{n}^{(1)}$ and
$\tilde{\mathcal{I}}_{n}^{(2)}$, as follows

$\displaystyle\mathcal{I}_{1}$ $\displaystyle=$
$\displaystyle-|\gamma^{(0)}_{k_{1}}(\tau_{I})|^{2}|A^{(0)}_{k_{2}}(\tau_{I})||A^{(0)}_{k_{3}}(\tau_{I})|~{}k_{2}k_{3}~{}\tilde{\mathcal{I}}_{n}^{(1)},$
(39) $\displaystyle\mathcal{I}_{2}$ $\displaystyle=$
$\displaystyle-|\gamma^{(0)}_{k_{1}}(\tau_{I})|^{2}|A^{(0)}_{k_{2}}(\tau_{I})||A^{(0)}_{k_{3}}(\tau_{I})|~{}\tilde{\mathcal{I}}_{n}^{(2)},$
(40)

where $\gamma^{(0)}_{k}(\tau_{I})$ and $A^{(0)}_{k}(\tau_{I})$ are the
asymptotic super horizon expansion of the mode functions (13) and (22), given
by

$\displaystyle|\gamma^{(0)}_{k}(\tau_{I})|$ $\displaystyle=$
$\displaystyle\frac{H}{\sqrt{k^{3}}}$ (41)
$\displaystyle|A_{k}^{(0)}(\tau_{I})|$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{1}{\sqrt{\lambda_{I}}}\;\frac{2^{n}}{\sqrt{2\pi}}\;\frac{\Gamma(n+\frac{1}{2})}{\sqrt{k}}(-k\tau_{I})^{-n},\;\text{for}\;n>-\frac{1}{2}\\\
\\\
\frac{1}{\sqrt{\lambda_{I}}}\;\frac{\sqrt{2\pi}}{2^{n+1}\Gamma(n+\frac{3}{2})}\;\frac{1}{\sqrt{k}}\;(-k\tau_{I})^{n+1},\;\text{for}\;n<-\frac{1}{2}\end{cases}$
(42)

and the explicit forms of the integrals $\tilde{\mathcal{I}}_{n}^{(1)}$ and
$\tilde{\mathcal{I}}_{n}^{(2)}$ corresponding to both $n>-\frac{1}{2}$ and
$n<-\frac{1}{2}$ are given in Appendix A.

One can straightforwardly substitute (37) and (38) into (36) and after
considerable simplification, the final result for the correlator
$\langle\gamma A_{\mu}A^{\mu}\rangle$ can be written as

$\displaystyle\langle\gamma({\bf k}_{1})A_{\mu}({\bf k}_{2})A^{\mu}({\bf
k}_{3})\rangle$ $\displaystyle=$ $\displaystyle(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf k_{3}})\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}$
(43) $\displaystyle\times$ $\displaystyle\left[{\cal
M}_{A}(k_{1},k_{2},k_{3},\tau_{I}){\cal
W}_{A}(k_{1},k_{2},k_{3})+\frac{1}{2}\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}P_{\gamma}(k_{1})P_{A}(\tilde{k}_{2})\right],$

with

$\displaystyle{\cal M}_{A}(k_{1},k_{2},k_{3},\tau_{I})$ $\displaystyle=$
$\displaystyle\frac{|\gamma^{(0)}_{k_{1}}(\tau_{I})|^{2}|A^{(0)}_{k_{2}}(\tau_{I})||A^{(0)}_{k_{3}}(\tau_{I})|}{a^{2}}\;,$
(44) $\displaystyle{\cal W}_{A}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\left[({\bf k_{2}}+{\bf k_{3}})^{2}-({\bf k_{2}}\cdot{\bf
k_{3}})\right]\frac{k_{2}\tilde{\mathcal{I}}_{n}^{(1)}}{k_{3}}+k_{2}^{2}\tilde{\mathcal{I}}_{n}^{(2)}.$
(45)

Here the subscript ‘$A$’ in ${\cal M}_{A}$ and $W_{A}$ refers to the gauge
field $A_{\mu}$. Note that, in the final result (43), the first contribution
arises from the interaction term corresponding to the leading order
cosmological diagram in figure 1 while the second term is the pure kinematical
correction term which is independent of the coupling function $\lambda(\tau)$.

####  3.1.2 $\langle\gamma B_{\mu}B^{\mu}\rangle$

In this section, we will calculate the correlator $\langle\gamma
B_{\mu}B^{\mu}\rangle$. The covariant magnetic field $B_{\mu}$ is defined with
respect to a comoving observer as in (24). Then $\gamma B_{\mu}B^{\mu}$ can be
expressed in terms of gauge field tensor as follows

$\displaystyle\gamma B_{\mu}B^{\mu}=\frac{1}{2a^{4}}\gamma
F_{ij}F_{ij}-\frac{1}{a^{4}}\gamma\gamma_{ij}F_{il}F_{jl}+{\cal
O}(\gamma^{3})\;.$ (46)

One can evaluate the correlator by substituting (46) and (33) into (30).

$\displaystyle\langle\gamma({\bf k}_{1},\tau_{I})B_{\mu}({\bf
k}_{2},\tau_{I})B^{\mu}({\bf k}_{3},\tau_{I})\rangle$ $\displaystyle=$
$\displaystyle\left\langle 0\right|\gamma({\bf k}_{1},\tau_{I})B_{\mu}({\bf
k}_{2},\tau_{I})B^{\mu}({\bf k}_{3},\tau_{I})\left|0\right>$ (47)
$\displaystyle+$ $\displaystyle i\int_{-\infty}^{\tau_{I}}d\tau\left\langle
0\right|[H_{\rm int}(\tau),\gamma({\bf k}_{1},\tau_{I})B_{\mu}({\bf
k}_{2},\tau_{I})B^{\mu}({\bf k}_{3},\tau_{I})]\left|0\right>$

As discussed in the case of the $\gamma A_{\mu}A^{\mu}$ correlator, the first
term in the RHS of (47) leads to a non zero contribution to the vacuum
expectation value of $\gamma B_{\mu}B^{\mu}$. This kinematical correction term
is calculated using (46) to yield the following contribution

$\displaystyle\left\langle 0\right|\gamma({\bf k}_{1},\tau_{I})B_{\mu}({\bf
k}_{2},\tau_{I})B^{\mu}({\bf
k}_{3},\tau_{I})\left|0\right>=-\frac{1}{2}(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\epsilon_{ij}\frac{k_{2i}k_{2j}}{\tilde{k}_{2}^{2}}P_{\gamma}(k_{1},\tau_{I})P_{B}(\tilde{k}_{2},\tau_{I})$
(48)

where ${\bf\tilde{k}}_{2}={\bf k}_{2}+\frac{1}{2}{\bf k}_{1}$. Once again, we
notice that the vacuum expectation value of $\gamma B_{\mu}B^{\mu}$
contributes at the leading order to the non-Gaussian correlator. However, we
would like to emphasise that these kinematical correction terms were not
explored in earlier works [17, 51].

The leading order contribution of the second term in the RHS of (47) is
calculated using the cosmological diagrammatic rules [21]. Then one would
obtain the total contribution as,

$\displaystyle\langle\gamma({\bf k}_{1},\tau_{I})B_{\mu}({\bf
k}_{2},\tau_{I})B^{\mu}({\bf k}_{3},\tau_{I})\rangle$ $\displaystyle=$
$\displaystyle-(2\pi)^{3}\delta^{(3)}({\bf k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}$ (49) $\displaystyle\times$
$\displaystyle\Biggl{[}\frac{k_{2}^{2}\mathcal{I}_{1}}{a^{4}}+\left\\{({\bf
k_{2}}+{\bf k_{3}})^{2}-({\bf k_{2}}\cdot{\bf
k_{3}})\right\\}\frac{k_{2}^{2}\mathcal{I}_{2}}{a^{4}}\Biggr{.}$
$\displaystyle+$
$\displaystyle\Biggl{.}\frac{1}{2}\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}P_{\gamma}(k_{1},\tau_{I})P_{B}(\tilde{k}_{2},\tau_{I})\Biggr{]},$

where $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ are the same integrals, as given
in (37) and (38). By inserting (39) and (40) into (49), the final result for
the correlation function becomes

$\displaystyle\left\langle\gamma({\bf k}_{1})B_{\mu}({\bf k}_{2})B^{\mu}({\bf
k}_{3})\right>$ $\displaystyle=$ $\displaystyle(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf k_{3}})\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}$
(50) $\displaystyle\times$ $\displaystyle\Biggl{[}{\cal
M}_{B}(k_{1},k_{2},k_{3},\tau_{I}){\cal
W}_{B}(k_{1},k_{2},k_{3})-\frac{1}{2}\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}P_{\gamma}(k_{1})P_{B}(\tilde{k}_{2})\Biggr{]},$

with

$\displaystyle{\cal M}_{B}(k_{1},k_{2},k_{3},\tau_{I})$ $\displaystyle=$
$\displaystyle\frac{|\gamma^{(0)}_{k_{1}}(\tau_{I})|^{2}|A^{(0)}_{k_{2}}(\tau_{I})||A^{(0)}_{k_{3}}(\tau_{I})|}{a^{4}}\
=\ \frac{1}{a^{2}}{\cal M}_{A}(k_{1},k_{2},k_{3},\tau_{I}),$ (51)
$\displaystyle{\cal W}_{B}(k_{1},k_{2},k_{3})$ $\displaystyle=$ $\displaystyle
k_{2}^{3}k_{3}~{}\tilde{\mathcal{I}}_{n}^{(1)}+\left[({\bf k_{2}}+{\bf
k_{3}})^{2}-({\bf k_{2}}\cdot{\bf
k_{3}})\right]k_{2}^{2}\,\tilde{\mathcal{I}}_{n}^{(2)}.$ (52)

These contributions borrow similar structure as derived earlier in the case of
the correlator $\left<\zeta\,{\bf B}\cdot{\bf B}\right>$ in [16, 17, 51]. Note
the presence of the two non-trivial integrals which can be exactly calculated
for various values of $n$. However, the second term in (50) is the pure
kinematical correction term which was not pointed out in earlier works [17,
51].

####  3.1.3 $\langle\gamma E_{\mu}E^{\mu}\rangle$

We now turn to the calculation of the correlator $\langle\gamma
E_{\mu}E^{\mu}\rangle$, i.e., non-Gaussian cross correlation of tensor mode
with the electric fields which is a bit more subtle than the other
correlators. Using equation (24), one can express the covariant electric field
in terms of gauge fields as

$\displaystyle E_{\mu}({\bf x},\tau)=\frac{d}{d\tau}A_{\mu}({\bf
x},\tau)+i[H_{\rm int},A_{\mu}],$ (53)

where the extra commutator $i[H_{\rm int},A_{\mu}]$ arises due to the fact
that the operators are all implicitly understood to be in the interaction
picture, as discussed in Appendix C. With this definition, the observable
$\gamma E_{\mu}E^{\mu}$ is expressed as follows

$\displaystyle\gamma
E_{\mu}E^{\mu}=\frac{1}{a^{4}}\left[\gamma\frac{dA_{i}}{d\tau}\frac{dA_{i}}{d\tau}+i\gamma\left(\frac{dA_{i}}{d\tau}[H_{\rm
int},A_{i}]+[H_{\rm
int},A_{i}]\frac{dA_{i}}{d\tau}\right)-\gamma\gamma_{ij}\frac{dA_{i}}{d\tau}\;\frac{dA_{j}}{d\tau}\right]+{\cal
O}(\gamma^{3})~{}.$ (54)

Again, using the master formula of the in-in formalism (30), we find that

$\displaystyle\langle\gamma({\bf k}_{1},\tau_{I})E_{\mu}({\bf
k}_{2},\tau_{I})E^{\mu}({\bf k}_{3},\tau_{I})\rangle$ $\displaystyle=$
$\displaystyle\left\langle 0\right|\gamma({\bf k}_{1},\tau_{I})E_{\mu}({\bf
k}_{2},\tau_{I})E^{\mu}({\bf k}_{3},\tau_{I})\left|0\right>$ (55)
$\displaystyle+$ $\displaystyle i\int_{-\infty}^{\tau_{I}}d\tau\left\langle
0\right|[H_{\rm int}(\tau),\gamma({\bf k}_{1},\tau_{I})E_{\mu}({\bf
k}_{2},\tau_{I})E^{\mu}({\bf k}_{3},\tau_{I})]\left|0\right>$

As discussed earlier, the first term in the RHS of the above equation produces
a non zero contribution to the vacuum expectation value of $\gamma
E_{\mu}E^{\mu}$ due to the second and the third term in the RHS of equation
(54). The contributions arising from the second term can be referred as the
dynamical correction term and the third term as the kinematical correction
term, respectively. In the super horizon limit, they can be calculated as
follows,

$\displaystyle\left\langle 0\right|\gamma({\bf k}_{1},\tau_{I})E_{\mu}({\bf
k}_{2},\tau_{I})E^{\mu}({\bf k}_{3},\tau_{I})\left|0\right>$ $\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}$
(56) $\displaystyle\times$
$\displaystyle\Bigg{[}-\left(\frac{k_{2}^{2}+k_{3}^{2}+{\bf k_{2}}\cdot{\bf
k_{3}}}{k_{3}^{2}}\right)\left(|A^{\prime}_{k_{3}}(\tau_{I})|^{2}+|A^{\prime}_{k_{2}}(\tau_{I})|^{2}\right)|\gamma_{k_{1}}(\tau_{I})|^{2}$
$\displaystyle+$
$\displaystyle\frac{1}{2}\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}P_{\gamma}(k_{1},\tau_{I})P_{E}(\tilde{k}_{2},\tau_{I})\Bigg{]}.$

Unlike previous sections, the above vacuum expectation value is not only
contributed by a kinematical correction term but also a dynamical correction
term. Finally, the leading order contribution arising from the second term in
the RHS of (55) is calculated using the cosmological diagrammatic rules [21]
and the final result becomes

$\displaystyle\left\langle\gamma({\bf k}_{1})E_{\mu}({\bf k}_{2})E^{\mu}({\bf
k}_{3})\right>$ $\displaystyle=$ $\displaystyle(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}\Bigg{[}-\left(\frac{({\bf
k_{2}}+{\bf k_{3}})^{2}-({\bf k_{2}}\cdot{\bf
k_{3}})}{k_{3}^{2}}\right)\frac{\mathcal{J}_{1}}{a^{4}}-k_{2}^{2}\,\frac{\mathcal{J}_{2}}{a^{4}}$
(57) $\displaystyle-$ $\displaystyle\left(\frac{k_{2}^{2}+k_{3}^{2}+{\bf
k_{2}}\cdot{\bf
k_{3}}}{k_{3}^{2}}\right)\left(\frac{|A^{\prime}_{k_{3}}(\tau_{I})|^{2}+|A^{\prime}_{k_{2}}(\tau_{I})|^{2}}{a^{4}}\right)|\gamma_{k_{1}}(\tau_{I})|^{2}$
$\displaystyle+$
$\displaystyle\frac{1}{2}\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}P_{\gamma}(k_{1})P_{E}(\tilde{k}_{2})\Bigg{]},$

wherein the two integrals $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$ are defined
in terms of the mode functions of tensor mode and gauge field and their time
derivatives as

$\displaystyle\mathcal{J}_{1}$ $\displaystyle=$ $\displaystyle 2\,{\rm
Im}\left[\gamma_{k_{1}}(\tau_{I})A^{\prime}_{k_{2}}(\tau_{I})A^{\prime}_{k_{3}}(\tau_{I})\int
d\tau\lambda(\tau)\gamma_{k_{1}}^{*}(\tau)A_{k_{2}}^{{}^{\prime}*}(\tau)A_{k_{3}}^{{}^{\prime}*}(\tau)\right],$
(58) $\displaystyle\mathcal{J}_{2}$ $\displaystyle=$ $\displaystyle 2\,{\rm
Im}\left[\gamma_{k_{1}}(\tau_{I})A^{\prime}_{k_{2}}(\tau_{I})A^{\prime}_{k_{3}}(\tau_{I})\int
d\tau\lambda(\tau)\gamma_{k_{1}}^{*}(\tau)A_{k_{2}}^{*}(\tau)A_{k_{3}}^{*}(\tau)\right].$
(59)

Again, it becomes useful to define the two auxiliary integrals
$\tilde{\mathcal{J}}_{n}^{(1)}$ and $\tilde{\mathcal{J}}_{n}^{(2)}$ associated
with $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$, in a similar manner as earlier,
by

$\displaystyle\mathcal{J}_{1}$ $\displaystyle=$
$\displaystyle-|\gamma^{(0)}_{k_{1}}(\tau_{I})|^{2}|A^{\prime(0)}_{k_{2}}(\tau_{I})||A^{\prime(0)}_{k_{3}}(\tau_{I})|~{}k_{2}k_{3}~{}\tilde{\mathcal{J}}_{n}^{(1)}$
(60) $\displaystyle\mathcal{J}_{2}$ $\displaystyle=$
$\displaystyle-|\gamma^{(0)}_{k_{1}}(\tau_{I})|^{2}|A^{\prime(0)}_{k_{2}}(\tau_{I})||A^{\prime(0)}_{k_{3}}(\tau_{I})|~{}\tilde{\mathcal{J}}_{n}^{(2)}$
(61)

where $\gamma^{(0)}_{k}(\tau_{I})$ and $A^{\prime(0)}_{k}(\tau_{I})$ are the
asymptotic super horizon values of the mode functions of the tensor mode and
the time derivative of mode functions of the gauge field respectively. In the
limit corresponding to the end of inflation $|k\tau_{I}|\to 0$, we find that
$\tilde{\mathcal{J}}_{n}^{(1)}\simeq\tilde{\mathcal{I}}_{n}^{(1)}$ and
$\tilde{\mathcal{J}}_{n}^{(2)}\simeq\tilde{\mathcal{I}}_{n}^{(2)}$, as shown
in appendix A. Note that, a few more terms appear in the result for the
correlator $\langle\gamma E_{\mu}E^{\mu}\rangle$ as compared to $\langle\gamma
B_{\mu}B^{\mu}\rangle$ since the covariant electric field is defined in terms
of the time derivative of the gauge field while the covariant magnetic field
only contains spatial derivatives of the gauge field. The final result can now
be written in a compact form as

$\displaystyle\left\langle\gamma({\bf k}_{1})E_{\mu}({\bf k}_{2})E^{\mu}({\bf
k}_{3})\right>$ $\displaystyle=$ $\displaystyle(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf k_{3}})\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}$
(62) $\displaystyle\times$ $\displaystyle\Bigg{[}\sum_{\alpha=1}^{2}{\cal
M}_{\alpha E}{\cal W}_{\alpha
E}+\frac{1}{2}\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}P_{\gamma}(k_{1})P_{E}(\tilde{k}_{2})\Bigg{]},$

where

$\displaystyle{\cal M}_{1E}(k_{1},k_{2},k_{3},\tau_{I})$ $\displaystyle=$
$\displaystyle\frac{|\gamma^{(0)}_{k_{1}}(\tau_{I})|^{2}|A^{\prime(0)}_{k_{2}}(\tau_{I})||A^{\prime(0)}_{k_{3}}(\tau_{I})|}{a^{4}},$
(63) $\displaystyle{\cal M}_{2E}(k_{1},k_{2},k_{3},\tau_{I})$ $\displaystyle=$
$\displaystyle-|\gamma_{k_{1}}(\tau_{I})|^{2}\left(\frac{|A^{\prime}_{k_{2}}(\tau_{I})|^{2}+|A^{\prime}_{k_{3}}(\tau_{I})|^{2}}{a^{4}}\right),$
(64) $\displaystyle{\cal W}_{1E}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\left[({\bf k_{2}}+{\bf k_{3}})^{2}-({\bf k_{2}}\cdot{\bf
k_{3}})\right]\frac{k_{2}\tilde{\mathcal{I}}_{n}^{(1)}}{k_{3}}+k_{2}^{2}\tilde{\mathcal{I}}_{n}^{(2)},$
(65) $\displaystyle{\cal W}_{2E}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\frac{k_{2}^{2}+k_{3}^{2}+{\bf k_{2}}\cdot{\bf
k_{3}}}{k_{3}^{2}}.$ (66)

Note that, the above result has been obtained in the super horizon limit. This
correlator has a similar structure as $\left<\gamma B_{\mu}B^{\mu}\right>$,
albeit an extra contribution which arises precisely due to the definition of
the covariant electric field in equation (53). We refer to this contribution
as the dynamical correction term to the correlator. Moreover, the last term in
equation (62) is similar to the kinematical correction term as appeared
earlier for $\left<\gamma B_{\mu}B^{\mu}\right>$ in equation (50). As a
result, the associated non-linearity parameter becomes constant in the
squeezed limit, as we shall discuss in the following section.

###  3.2 The magnetic and electric non-linearity parameters

Following the approach of our previous papers [16, 17] for the case of a
cross-correlation of the comoving curvature perturbation with primordial gauge
fields, we can similarly define convenient and useful non-linearity parameters
to characterise the cross-correlations of the inflationary tensor perturbation
with primordial gauge fields. Let us first define the bispectra associated
with $\left<\gamma B_{\mu}B^{\mu}\right>$ and $\left<\gamma
E_{\mu}E^{\mu}\right>$ as,

$\displaystyle\left<\gamma(\bm{\mathrm{{k}}}_{1})B_{\mu}(\bm{\mathrm{{k}}}_{2})B^{\mu}(\bm{\mathrm{{k}}}_{3})\right>$
$\displaystyle\equiv$
$\displaystyle(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\,{\cal
B}_{\gamma BB}(k_{1},k_{2},k_{3})~{},$ (67)
$\displaystyle\left<\gamma(\bm{\mathrm{{k}}}_{1})E_{\mu}(\bm{\mathrm{{k}}}_{2})E^{\mu}(\bm{\mathrm{{k}}}_{3})\right>$
$\displaystyle\equiv$
$\displaystyle(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\,{\cal
B}_{\gamma EE}(k_{1},k_{2},k_{3})~{},$ (68)

where the bispectra only depend on the magnitude of the wavevectors due to
spatial isotropy of the background expansion. The strength of the magnetic and
electric field bispectra can be characterised by defining the non-linearity
parameters $b_{NL}^{\gamma}$ and $e_{NL}^{\gamma}$ as follows

$\displaystyle{\cal B}_{\gamma BB}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\frac{1}{2}b_{NL}^{\gamma}P_{\gamma}(k_{1})\big{[}P_{B}(k_{2})+P_{B}(k_{3})\big{]},$
(69) $\displaystyle{\cal B}_{\gamma EE}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\frac{1}{2}e_{NL}^{\gamma}P_{\gamma}(k_{1})\big{[}P_{E}(k_{2})+P_{E}(k_{3})\big{]},$
(70)

where $P_{\gamma}$, $P_{B}$ and $P_{E}$ are the power spectra of the tensor
perturbations (11), magnetic field (28) and electric field (29), respectively
and the RHS of these definitions are written in a symmetric manner. If the two
non-linearity parameters $b_{NL}^{\gamma}$ and $e_{NL}^{\gamma}$ are momentum
independent, they correspond to a local shape of the bispectra in the Fourier
space which can be obtained from the relations [16, 17]

$\displaystyle{\bf B}$ $\displaystyle=$ $\displaystyle{\bf
B}^{G}+\frac{1}{2}b_{NL}^{\gamma}\gamma^{G}{\bf B}^{G}~{},$ (71)
$\displaystyle{\bf E}$ $\displaystyle=$ $\displaystyle{\bf
E}^{G}+\frac{1}{2}e_{NL}^{\gamma}\gamma^{G}{\bf E}^{G}~{},$ (72)

where ${\bf B}^{G}$, ${\bf E}^{G}$ and $\gamma^{G}$ are the Gaussian random
fields. In the absence of any parity violating interactions, as is the case in
our scenario, the Gaussian tensor mode $\gamma^{G}$ represents both the
polarisation degrees of freedom associated with the graviton. Similar
definitions also allow one to compute the higher order cross correlations of
curvature perturbations with gauge fields [39].

For the case of the $\langle\gamma B_{\mu}B^{\mu}\rangle$ correlator, we can
now simply compare the equations (67) and (69) with (50) and read off the
bispectrum which leads to the following result for the non-linearity parameter
$b_{NL}^{\gamma}$ for $n+1/2>0$ as

$\displaystyle
b_{NL}^{\gamma}=\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}\Biggl{[}\left(\frac{(k_{2}k_{3})^{-(n+\frac{1}{2})}}{k_{2}^{1-2n}+k_{3}^{1-2n}}\right){\cal
W}_{B}(k_{1},k_{2},k_{3})-\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}\left(\frac{\tilde{k}_{2}^{{1-2n}}}{k_{2}^{1-2n}+k_{3}^{1-2n}}\right)\Bigg{]},$
(73)

where, once again, the first term in the above result is due to the
interaction Hamiltonian and the second term is the pure kinematical correction
term. Similarly, on comparing (68) and (70) with (57), yields the result for
the non-linearity parameter $e_{NL}^{\gamma}$ for $n-1/2>0$ as

$\displaystyle
e_{NL}^{\gamma}=\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}\Biggl{[}\left(\frac{(k_{2}k_{3})^{(\frac{3}{2}-n)}}{k_{2}^{3-2n}+k_{3}^{3-2n}}\right){\cal
W}_{1E}(k_{1},k_{2},k_{3})-{\cal
W}_{2E}(k_{1},k_{2},k_{3})+\left(\frac{k_{2}}{\tilde{k}_{2}}\right)^{2}\\!\\!\left(\frac{\tilde{k}_{2}^{{3-2n}}}{k_{2}^{3-2n}+k_{3}^{3-2n}}\right)\Bigg{]}.$
(74)

These general expressions for $b_{NL}^{\gamma}$ and $e_{NL}^{\gamma}$ are one
of the main results of this paper and they depend only on the wave numbers
$k_{1}$, $k_{2}$ and $k_{3}$. The presence of
$\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})$
in various cross-correlations constrains the wave vectors to form a closed
triangle in the Fourier space. In order to understand the extent of these non-
linearities, it is useful to study $b_{NL}^{\gamma}$ and $e_{NL}^{\gamma}$ for
different shapes of the triangle configuration. Conventionally, such shape
functions are studied by defining the dimensionless variables
$x_{1}=\frac{k_{1}}{k_{2}}$ and $x_{3}=\frac{k_{3}}{k_{2}}$ while $k_{2}$ is
set at an arbitrary scale. For an isotropic background, the shape function
only depends on any of the two momentum ratios. This allows us to first
calculate the following two expressions as

$\displaystyle\left(\frac{\tilde{k}_{2}}{k_{2}}\right)^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1+x_{3}^{2}-\frac{x_{1}^{2}}{2}\right),$ (75)
$\displaystyle\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\frac{(x_{3}^{2}-x_{1}^{2}-1)^{2}}{4x_{1}^{2}}-1\right).$
(76)

In order to calculate $b_{NL}^{\gamma}$ and $e_{NL}^{\gamma}$, we also need to
evaluate the integrals $\tilde{\mathcal{I}}_{n}^{(1)}$ and
$\tilde{\mathcal{I}}_{n}^{(2)}$ for different values of $n$ which has already
been done in [17]. For the most interesting case of $n=2$ which leads to a
scale invariant spectrum of magnetic fields, we present the results for
$\tilde{\mathcal{I}}_{n}^{(1)}$ and $\tilde{\mathcal{I}}_{n}^{(2)}$ in
Appendix A and use them to plot the extent of these non-linearity parameters
as contour plots in the momentum space in figure 2. There are a few
interesting points to note about the plots. First, both $b_{NL}^{\gamma}$ and
$e_{NL}^{\gamma}$ remain very small in the squeezed limit which is
characterised by the top left corner, i.e., $x_{1}\to 0,x_{3}\to 1$. Close to
the flattened limit defined by $x_{1}\to 2,x_{3}\to 1$, $b_{NL}^{\gamma}$
becomes large but $e_{NL}^{\gamma}$ becomes very small in the exact flattened
limit. However, $e_{NL}^{\gamma}$ becomes large in a somewhat different limit
$x_{1}\to 1,x_{3}\to 0$ in which $b_{NL}^{\gamma}$ remains small. Note that,
the allowed region in the space of $x_{1}-x_{3}$ plane is constrained by the
triangle inequality in the Fourier space.

Figure 2: The extent of the non-linearity parameters $b_{NL}^{\gamma}$ (left)
and $e_{NL}^{\gamma}$ (right) corresponding to different triangular
configuration are plotted for the case of $n=2$. The color legends
representing the magnitude are also shown below each plot. While the top left
corner of these plots corresponds to the squeezed limit, i.e.,
$(x_{1},x_{3})=(0,1)$, the top right corner represents the flattened shape,
i.e., $(x_{1},x_{3})=(2,1)$. In the squeezed limit, $b_{NL}^{\gamma}$ is
relatively small while it is large in the nearly flattened shape. However,
$e_{NL}^{\gamma}$ remains small both in the exact squeezed and flattened limit
but it becomes large as $x_{3}$ becomes smaller.

###  3.3 The squeezed limit

Let us now consider the squeezed limit wherein the wavelength of the graviton
mode is much larger than the corresponding wavelengths of the electromagnetic
fields. Thus, a squeezed shape non-Gaussianity implies a correlation of very
large scale with very small scales. Using our general results of the in-in
calculations from the previous sections, we can easily write down the results
for all the correlators in the squeezed limit. However, as we shall discuss
later, one can also obtain the squeezed limit non-Gaussianity by using semi-
classical methods. In this limit, we have $\bm{\mathrm{{k}}}_{1}\to 0$ and
$\bm{\mathrm{{k}}}_{2}\to-\bm{\mathrm{{k}}}_{3}\equiv\bm{\mathrm{{k}}}$. The
two integrals of interest $\tilde{\mathcal{I}}_{n}^{(1)}$ and
$\tilde{\mathcal{I}}_{n}^{(2)}$ simplify as

$\displaystyle\tilde{\mathcal{I}}_{n}^{(1)}$ $\displaystyle=$
$\displaystyle\begin{cases}\pi\int^{\tau_{I}}d\tau\tau
J_{n-1/2}(-k\tau)Y_{n-1/2}(-k\tau),\quad n>-\frac{1}{2}\\\
-\pi\int^{\tau_{I}}d\tau\tau J_{n-1/2}(-k\tau)Y_{n-1/2}(-k\tau),\quad
n<-\frac{1}{2}\end{cases}$ (77) $\displaystyle\tilde{\mathcal{I}}_{n}^{(2)}$
$\displaystyle=$ $\displaystyle\begin{cases}\pi\int^{\tau_{I}}d\tau\tau
J_{n+1/2}(-k\tau)Y_{n+1/2}(-k\tau),\quad n>-\frac{1}{2}\\\
-\pi\int^{\tau_{I}}d\tau\tau J_{n+1/2}(-k\tau)Y_{n+1/2}(-k\tau),\quad
n<-\frac{1}{2}\end{cases}$ (78)

where $J_{\alpha}$ and $Y_{\alpha}$ are the Bessel functions of the first and
second kind, respectively. These integrals are straightforward to evaluate and
for integer values of $n$, one can show the following analytical result
as777Note that, this result also holds true for real non-integer values of $n$
which can be verified numerically. [17]

$\pi\int^{\tau_{I}}d\tau\tau J_{n}(-k\tau)Y_{n}(-k\tau)=\frac{n}{k^{2}},$ (79)

and evidently,
$\tilde{\mathcal{I}}_{n}^{(2)}=\tilde{\mathcal{I}}_{n+1}^{(1)}$. In the
squeezed limit, $\bm{\mathrm{{k}}}_{1}\to 0$, so using this result in all the
correlators that we have obtained earlier, and after some simplification, we
find the following results

$\displaystyle\langle\gamma(\bm{\mathrm{{k}}}_{1})A_{\mu}(\bm{\mathrm{{k}}}_{2})A^{\mu}(\bm{\mathrm{{k}}}_{3})\rangle$
$\displaystyle=$ $\displaystyle\begin{cases}(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{i}k_{j}}{k^{2}}\left(n+\frac{1}{2}\right)P_{\gamma}(k_{1})P_{A}(k),\;\text{if}\;n>-\frac{1}{2}\\\
-(2\pi)^{3}\delta^{(3)}({\bf k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{i}k_{j}}{k^{2}}\left(n-\frac{1}{2}\right)P_{\gamma}(k_{1})P_{A}(k),\;\text{if}\;n<-\frac{1}{2}\end{cases}$
(80)
$\displaystyle\langle\gamma(\bm{\mathrm{{k}}}_{1})B_{\mu}(\bm{\mathrm{{k}}}_{2})B^{\mu}(\bm{\mathrm{{k}}}_{3})\rangle$
$\displaystyle=$ $\displaystyle\begin{cases}(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{i}k_{j}}{k^{2}}\left(n-\frac{1}{2}\right)P_{\gamma}(k_{1})P_{B}(k),\;\text{if}\;n>-\frac{1}{2}\\\
-(2\pi)^{3}\delta^{(3)}({\bf k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{i}k_{j}}{k^{2}}\left(n+\frac{1}{2}\right)P_{\gamma}(k_{1})P_{B}(k),\;\text{if}\;n<-\frac{1}{2}\end{cases}$
(81)
$\displaystyle\langle\gamma(\bm{\mathrm{{k}}}_{1})E_{\mu}(\bm{\mathrm{{k}}}_{2})E^{\mu}(\bm{\mathrm{{k}}}_{3})\rangle$
$\displaystyle=$ $\displaystyle\begin{cases}(2\pi)^{3}\delta^{(3)}({\bf
k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{i}k_{j}}{k^{2}}\left(n-\frac{1}{2}\right)P_{\gamma}(k_{1})P_{E}(k),\;\text{if}\;n>\frac{1}{2}\\\
-(2\pi)^{3}\delta^{(3)}({\bf k_{1}}+{\bf k_{2}}+{\bf
k_{3}})\epsilon_{ij}\frac{k_{i}k_{j}}{k^{2}}\left(n+\frac{1}{2}\right)P_{\gamma}(k_{1})P_{E}(k),\;\text{if}\;n<\frac{1}{2}\end{cases}$
(82)

The above results comprise of a set of novel consistency relations for the
cross correlations involving the tensor perturbation and gauge fields. Note
that, in the absence of other source terms, there exists an electromagnetic
duality which states that the equations of motion and the power spectra of the
electric and magnetic fields can be obtained from each other under a
simultaneous exchange of ${\bf B}\to{\bf E}$, ${\bf E}\to-{\bf B}$ and
$n\to-n$ [52]. To our knowledge, it has not been known if this duality is also
maintained at the level of higher order correlation functions beyond the power
spectra. In our case, we observe that this duality is indeed preserved in the
squeezed limit of the cross correlations involving the electric and magnetic
fields. In particular, the results in equation (81) can be obtained from
equation (82) using the electromagnetic duality and vice versa. The results in
both the regimes of the $\langle\gamma E_{\mu}E^{\mu}\rangle$ correlator can
be correctly obtained from the corresponding results for the $\langle\gamma
B_{\mu}B^{\mu}\rangle$ correlator under the electromagnetic duality
transformations. In any case, it will be interesting to examine whether the
electromagnetic duality is maintained for the higher order correlation
functions away from the squeezed limit.

##  4 Semi-classical derivation of the consistency relations

In this section, we want to apply semi-classical techniques to derive the
three-point correlation function of the form
$\left<\gamma(\tau_{I},\bm{\mathrm{{k}}}_{1})Y_{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{2})Y^{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{3})\right>$
for $Y_{\mu}=\left\\{A_{\mu},B_{\mu},E_{\mu}\right\\}$ in the squeezed limit
$k_{1}\ll k_{2}\approx k_{3}$. We would use the similar approach developed
earlier to compute loop corrections to the primordial power spectrum during a
de-Sitter/slow roll background [38, 53]. To start with, the correlator
$\langle\gamma Y_{\mu}Y^{\mu}\rangle$ in real space can be written as

$\displaystyle\left<\gamma(\bm{\mathrm{{x}}})Y_{\mu}({\bm{\mathrm{{x}}}})Y^{\mu}({\bm{\mathrm{{x}}}})\right>=\int\frac{d^{3}\bm{\mathrm{{k}}}_{1}}{(2\pi)^{3}}\int\frac{d^{3}\bm{\mathrm{{k}}}_{2}}{(2\pi)^{3}}\int\frac{d^{3}\bm{\mathrm{{k}}}_{3}}{(2\pi)^{3}}e^{i\bm{\mathrm{{x}}}\cdot(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})}\left<\gamma(\bm{\mathrm{{k}}}_{1})Y_{\mu}(\bm{\mathrm{{k}}}_{2})Y^{\mu}(\bm{\mathrm{{k}}}_{3})\right>.$
(83)

In the limit where one of the momenta is much smaller than the other two i.e.
$k_{1}\ll k_{2}\approx k_{3}$, the wavelength corresponding to the graviton
$\gamma$ is very large and is frozen to a constant amplitude $\gamma^{B}$
outside the horizon. As a result, the effect due to a long wavelength graviton
mode $\gamma^{B}$ on the two-point correlator of gauge fields can be
understood by rescaling the background as
$ds^{2}=-dt^{2}+a^{2}(t)d\tilde{x}^{2}$ with $d\tilde{x}^{2}\to
dx^{2}+\gamma^{B}_{ij}dx^{i}dx^{j}$. Thus we can compute the two point
correlation function $\langle Y_{\mu}Y^{\mu}\rangle$ in the modified
background using the similar arguments of [38], by Taylor expanding around the
unperturbed background in real space as

$\displaystyle\left<Y_{\mu}(\bm{\mathrm{{x}}})Y^{\mu}(\bm{\mathrm{{x}}})\right>_{B}=\left<Y_{\mu}(\bm{\mathrm{{x}}})Y^{\mu}(\bm{\mathrm{{x}}})\right>_{0}+\gamma_{ij}^{B}\frac{\partial}{\partial\gamma^{B}_{ij}}\left.\left<Y_{\mu}({\tilde{\bf
x}})Y^{\mu}({\tilde{\bf x}})\right>\right|_{\gamma^{B}=0}+\dots\,,$ (84)

and in the Fourier space as

$\displaystyle\left<Y_{\mu}(\bm{\mathrm{{x}}})Y^{\mu}(\bm{\mathrm{{x}}})\right>_{B}$
$\displaystyle=$
$\displaystyle\int\int\frac{d^{3}\bm{\mathrm{{k}}}_{2}}{(2\pi)^{3}}\frac{d^{3}\bm{\mathrm{{k}}}_{3}}{(2\pi)^{3}}e^{i\bm{\mathrm{{x}}}\cdot\bm{\mathrm{{k}}}_{2}}e^{i\bm{\mathrm{{x}}}\cdot\bm{\mathrm{{k}}}_{3}}\left<Y_{\mu}({\bm{\mathrm{{k}}}_{2}})Y^{\mu}({\bm{\mathrm{{k}}}_{3}})\right>_{B}$
(85) $\displaystyle=$
$\displaystyle\left<Y_{\mu}(\bm{\mathrm{{x}}})Y^{\mu}(\bm{\mathrm{{x}}})\right>_{0}$
$\displaystyle+\left(\gamma^{B}_{ij}(\bm{\mathrm{{x}}}_{0})\frac{\partial}{\partial\gamma^{B}_{ij}}+\frac{1}{2}\gamma^{B}_{ij}(\bm{\mathrm{{x}}}_{0})\gamma^{B}_{lk}(\bm{\mathrm{{x}}}_{0})\frac{\partial}{\partial\gamma^{B}_{ij}}\frac{\partial}{\partial\gamma^{B}_{kl}}\right)$
$\displaystyle\times\left.\left[\int\int\frac{d^{3}\tilde{\bm{\mathrm{{q}}}}_{1}}{(2\pi)^{3}}\frac{d^{3}\tilde{\bm{\mathrm{{q}}}}_{2}}{(2\pi)^{3}}e^{i\bm{\mathrm{{\tilde{x}}}}\cdot\bm{\mathrm{{\tilde{q}}}}_{1}}e^{i\tilde{\bm{\mathrm{{x}}}}\cdot\bm{\mathrm{{\tilde{q}}}}_{2}}\left<Y_{\mu}({\tilde{\bm{\mathrm{{q}}}}}_{1})Y^{\mu}({\tilde{\bm{\mathrm{{q}}}}}_{2})\right>_{B}\right]\right|_{\gamma^{B}=0}\
.$

Since the long wavelength mode $\gamma_{B}$ is almost constant over the Hubble
horizon, and in particular, over the scales of variation of the short
wavelength modes, we can choose $\bm{\mathrm{{x}}}_{0}=\bm{\mathrm{{x}}}$. We
can now express the modified coordinates due to the long wavelength graviton
as

${\tilde{x}}^{2}=\left[e^{\gamma^{B}}\right]_{ij}x^{i}x^{j}=x^{2}+\gamma^{B}_{ij}x^{i}x^{j}+\dots~{},$
(86)

and using that momentum transforms inversely, define

${\tilde{q}}^{2}=\left[e^{-\gamma^{B}}\right]^{ij}q_{i}q_{j}=q^{2}-\gamma^{B}_{ij}q_{i}q_{j}+\dots~{},$
(87)

such that
$\bm{\mathrm{{\tilde{x}}}}\cdot\bm{\mathrm{{\tilde{q}}}}=\bm{\mathrm{{x}}}\cdot\bm{\mathrm{{q}}}$,
and

$\frac{\partial}{\partial\gamma_{ij}^{B}}=\frac{\partial\tilde{q}^{2}}{\partial\gamma_{ij}^{B}}\frac{\partial}{\partial\tilde{q}^{2}}=-q_{i}q_{j}\frac{\partial}{\partial\tilde{q}^{2}}~{}.$
(88)

It is useful to note that $\det[\exp(\gamma^{B})]=1$, so that the Jacobian of
the transformation $\tilde{q}\to q$ is trivial, and
$\delta^{(3)}(\tilde{\bm{\mathrm{{q}}}}_{1}+\tilde{\bm{\mathrm{{q}}}}_{2})=\delta^{(3)}(\bm{\mathrm{{q}}}_{1}+\bm{\mathrm{{q}}}_{2})$,
so we can write

$\displaystyle\left<Y_{\mu}(\bm{\mathrm{{\tilde{q}}}}_{1})Y^{\mu}(\bm{\mathrm{{\tilde{q}}}}_{2})\right>_{B}$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{q}}}_{1}+\bm{\mathrm{{q}}}_{2})\left[1+\left(-\gamma^{B}_{ij}q_{i}q_{j}+\frac{1}{2}\gamma^{B}_{il}\gamma^{B}_{lj}q_{i}q_{j}+\dots\right)\frac{\partial}{\partial
q^{2}}\right.$ (89)
$\displaystyle\left.+\frac{1}{2}\left((\gamma^{B}_{ij}q_{i}q_{j})^{2}+\dots\right)\frac{\partial}{\partial
q^{2}}\frac{\partial}{\partial q^{2}}+\dots\right]P_{Y}(q).$

In the squeezed limit, due to the rescaled background by the long wavelength
graviton mode, one can write a three point correlation function in terms of
the modified two point function as

$\displaystyle\lim_{k_{1}\to
0}\left<\gamma(\tau_{I},\bm{\mathrm{{k}}}_{1})Y_{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{2})Y^{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{3})\right>=\left<\gamma(\tau_{I},\bm{\mathrm{{k}}}_{1})\left<Y_{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{2})Y^{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{3})\right>_{B}\right>.$
(90)

Let us now justify that the arguments used to obtain the consistency relations
might only be trusted for specific regions of $n$. To illustrate this, let’s
write the evolution equation for the canonically normalised gauge field, i.e.,
equation (19) again as

$v_{k}^{\prime\prime}+\left(k^{2}-\frac{S^{\prime\prime}}{S}\right)v_{k}=0~{}.$
(91)

This equation is identical to the evolution equation of a massless scalar
field in the FLRW spacetime. However, the only difference is that the gauge
field experiences an effective scale factor $S$ which we refer to as the pump
field. For the massless scalar field in de-Sitter with $a\sim 1/\tau$, there
exists an effective event horizon and therefore, the long and the short
wavelength modes decouple on the superhorizon scales. This allows one to
capture all the physical effects of the long wavelength mode by rescaling the
background for the short wavelength modes to calculate the three point cross-
correlations involving long and short wavelength modes of the scalar field
fluctuations or the curvature perturbations. Similar arguments can be applied
for any other correlator in the squeezed limit as far as the super horizon
modes are fixed by the values of the background quantities at horizon
crossing.

If we inspect the mode function equation (91) above in the long wavelength
limit $k^{2}\ll S^{\prime\prime}/S$, the solution can be found by simple
integration to be given by

$v_{k}(\tau)\sim C_{1}S(\tau)+C_{2}S(\tau)\int\frac{d\tau}{S^{2}(\tau)}$ (92)

where $C_{1}$ and $C_{2}$ are $k$-dependent integration constants. It is clear
that if $S={\sqrt{\lambda}}\sim\tau^{-n}$, then for $\tau\to 0$, the first
term will dominate for $n>-1/2$, which means that the solution for
$A_{k}=v_{k}/S$ will freeze to a constant in the long wavelength limit. A
trivial rescaling of the field is of course not important, what is important
is that the solution is in a super horizon attractor regime. If the solution
is dominated by the attractor, then the rescaling of the momentum by a long
mode will enter trivially as a rescaling of the initial normalisation through
$C_{1}$. However, if the non-attractor is important, the rescaling of momentum
will also enter through the lower limit of the time integral multiplying
$C_{2}$ and thus modify the dynamics non-trivially. The existence of an
attractor solution of the long wavelength modes leads to the constraint
$n>-1/2$. Therefore, the semi-classical derivation for $\langle\gamma
A_{\mu}A^{\mu}\rangle$ can only be trusted for $n>-1/2$. A similar condition
can also be obtained for the corresponding correlator involving the graviton
and magnetic fields. However, we find that the properly rescaled electric
field experiences an effective scale factor, $\frac{1}{\sqrt{\lambda}}$, which
reverses the condition on $n$ for electric field. As a result, we conclude
that the semi-classical derivation for the graviton electric fields cross-
correlator can only be trusted for $n<1/2$.

By substituting equation (89) into equation (90) and after further
simplifications, we will obtain the desired correlators in the squeezed limit.
For $n>-1/2$, we get the following consistency relations

$\lim_{k_{1}\to
0}\left<\gamma(\tau_{I},\bm{\mathrm{{k}}}_{1})A_{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{2})A^{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{3})\right>=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\left(n+\frac{1}{2}\right)\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}P_{\gamma}(k_{1})P_{A}(k_{2}),$
(93) $\lim_{k_{1}\to
0}\left<\gamma(\tau_{I},\bm{\mathrm{{k}}}_{1})B_{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{2})B^{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{3})\right>=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\left(n-\frac{1}{2}\right)\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}P_{\gamma}(k_{1})P_{B}(k_{2}),$
(94)

Similarly, for $n<1/2$, we get

$\lim_{k_{1}\to
0}\left<\gamma(\tau_{I},\bm{\mathrm{{k}}}_{1})E_{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{2})E^{\mu}(\tau_{I},\bm{\mathrm{{k}}}_{3})\right>=-(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2}+\bm{\mathrm{{k}}}_{3})\left(n+\frac{1}{2}\right)\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}}P_{\gamma}(k_{1})P_{E}(k_{2}).$
(95)

In the squeezed limit, these consistency relations are quite general as they
have been derived only by using the semi-classical techniques. These relations
do exactly match with the squeezed limit results of the full in-in correlates,
as obtained in (80), (81) and (82). From these correlators, one can now read
off the non-linearity parameters $b^{\gamma}_{NL}$ and $e^{\gamma}_{NL}$ as

$\displaystyle b^{\gamma}_{NL}$ $\displaystyle=$
$\displaystyle\left(n-\frac{1}{2}\right)\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}},\quad
n>-1/2$ (96) $\displaystyle e^{\gamma}_{NL}$ $\displaystyle=$
$\displaystyle-\left(n+\frac{1}{2}\right)\epsilon_{ij}\frac{k_{2i}k_{2j}}{k_{2}^{2}},\quad
n<1/2.$ (97)

It is important to mention here that these are indeed the local non-linearity
parameters as the factor $\epsilon_{ij}k_{2i}k_{2j}/{k_{2}^{2}}$ is momentum
independent in the squeezed limit, as evident from equation (76). Moreover,
both the parameters in the squeezed limit are proportional to $n$ and are of
order unity as expected, following the discussion of our earlier paper [16].
This behaviour is also quite evident from the contour plots in figure 2
wherein the top left corner indicates the extent of these non-linearity
parameters in the squeezed limit.

##  5 A novel correlation of tensor and curvature perturbations

As mentioned in the introduction, an interesting consequence of these three
point correlators is the possibility of a direct non-vanishing two-point
correlation of the primordial tensor mode with the primordial curvature
perturbation. We refer to this two-point function as a novel observable as it
usually vanishes in the single field models of inflation and is exclusively
induced by the excitation of gauge fields during inflation.

The presence of primordial gauge fields can leave indirect imprints on the
spectrum of curvature perturbation during inflation. This can be illustrated
by splitting the total curvature perturbations into the “usual” contribution,
$\zeta_{0}$ generated independently of the gauge field, and the contributions
imprinted by magnetic and electric fields,

$\displaystyle\zeta(\tau)=\zeta_{0}(\tau)+\zeta_{B}(\tau)+\zeta_{E}(\tau),$
(98)

where $\zeta_{B}$ and $\zeta_{E}$ are induced by the corresponding magnetic
and electric fields derived from the gauge fields only. Now, we can express
the two-point correlation function of the tensor mode with the curvature
perturbation as follows,

$\left<\gamma({\bf k_{1}})\zeta({\bf k_{2}})\right>=\left<\gamma({\bf
k_{1}})\zeta_{B}({\bf k_{2}})\right>+\left<\gamma({\bf k_{1}})\zeta_{E}({\bf
k_{2}})\right>.$ (99)

To this end, we are going to generically denote $\zeta_{B}$ and $\zeta_{E}$ by
$\zeta_{Y}$. Following [39], they can be written in the Fourier space as

$\displaystyle\zeta_{Y}({\bf
k},\tau_{I})=\int_{\tau_{0}}^{\tau_{I}}d\ln\tau\lambda(\tau)\int\frac{d^{3}\bm{\mathrm{{q}}}}{(2\pi)^{3}}\frac{Y_{i}({\bf
q})Y^{i}({\bf k-q})}{3H^{2}\epsilon}~{},$ (100)

where $\epsilon$ is the first slow roll parameter, $H$ is the Hubble
parameter, and $Y=\\{E,B\\}$. We can now express the two-point correlator
$\left<\gamma\zeta_{Y}\right>$ in terms of the three point function
$\langle\gamma Y_{\mu}Y^{\mu}\rangle$ as

$\displaystyle\left<\gamma({\bf k_{1}},\tau_{I})\zeta_{Y}({\bf
k_{2}},\tau_{I})\right>=\int_{\tau_{0}}^{\tau_{I}}d\ln\tau\lambda(\tau)\int\frac{d^{3}\bm{\mathrm{{q}}}}{(2\pi)^{3}}\frac{\left<\gamma({\bf
k}_{1},\tau_{I})Y_{i}({\bf q},\tau)Y^{i}({\bf
k_{2}-q},\tau)\right>}{3H^{2}\epsilon}~{}.$ (101)

To compute $\left<\gamma\zeta\right>$, we need to evaluate this three
dimensional integral in Fourier space. However, note that, these three point
correlators appearing in the integrand of above equation are not equal time
correlation functions [39]. Thus, we first have to convert them in terms of an
equal time correlation function as

$\displaystyle\left<\gamma({\bf k_{1}},\tau_{I})\zeta_{Y}({\bf
k_{2}},\tau_{I})\right>=\int_{\tau_{0}}^{\tau_{I}}d\tau\frac{\lambda(\tau_{I})}{\tau}\left(\frac{\tau}{\tau_{I}}\right)^{\alpha_{Y}}\int\frac{d^{3}\bm{\mathrm{{q}}}}{(2\pi)^{3}}\frac{\left<\gamma({\bf
k}_{1},\tau_{I})Y_{i}({\bf q},\tau_{I})Y^{i}({\bf
k_{2}-q},\tau_{I})\right>}{3H^{2}\epsilon}$ $\displaystyle\times\;\theta(|{\bf
k_{2}}-{\bf
q}|-k_{0})\theta(q-k_{0})\theta\left(-\frac{1}{\tau}-q\right)\theta\left(-\frac{1}{\tau}-|{\bf
k_{2}}-{\bf q}|\right),$ (102)

where

$\displaystyle\alpha_{B}$ $\displaystyle=$
$\displaystyle\begin{cases}4-2n,\;n>-1/2\\\ 6+2n,\;n<-1/2\end{cases}$ (103)
$\displaystyle\alpha_{E}$ $\displaystyle=$
$\displaystyle\begin{cases}6-2n,\;n>1/2\\\ 4+2n,\;n<1/2\end{cases}$ (104)

Here, we have assumed that only super horizon modes of gauge fields are
contributing non trivially to these integrals with $\theta$ functions
stressing this fact and $k_{0}\equiv a_{0}H_{0}$ is the horizon crossing scale
at the time $\tau_{0}$ which we assume to be the onset of generation of
magnetic fields.

In general, the momentum integral appearing in equation (5) is somewhat non-
trivial to evaluate. However, one can analytically calculate these correlators
by assuming that the magnetic fields are produced after the horizon exit of
the graviton mode corresponding to the momentum $k_{1}$ and thus, the primary
contribution to above integral comes from the local shape of the bispectrum.
Thus in the limit $k_{1}\ll q$, we can use the local Ansatzes (71) and (72) to
get an estimate of these correlators. We begin with the expansions given in
equations (69) and (70), and take $b^{\gamma}_{NL}$ and $e^{\gamma}_{NL}$ to
be momentum independent. Let’s first focus on the
$\left\langle\gamma\zeta_{B}\right>$ which can be evaluated

$\displaystyle\left<\gamma({\bf k_{1}})\zeta_{B}({\bf
k_{2}})\right>=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2})\;b^{\gamma}_{NL}\left(\frac{\lambda(\tau_{I}){\cal
P}_{B}(k_{2})}{3H^{2}\epsilon}\right)P_{\gamma}(k_{1})\;F(N_{0},N_{2},\alpha_{B})~{},$
(105)

with

$\displaystyle
F(N_{0},N_{2},\alpha_{Y})=\frac{e^{N_{2}\alpha_{Y}}}{\alpha_{Y}}\Biggl{[}\frac{1-e^{-N_{0}\alpha_{Y}}}{\alpha_{Y}}-N_{0}\Biggr{]}~{}.$
(106)

Here, $N_{0}$ denotes the number of e-folds from the onset of magnetic field
generation to the end of inflation, $N_{2}$ denotes the number of e-folds from
the horizon exit of the mode $k_{2}$ to the end of inflation and ${\cal
P}_{B}$ is the dimensionless power spectrum of magnetic field which is defined
similar to equation (12). Similarly, one can also define ${\cal P}_{E}$ and
obtain $\left\langle\gamma\zeta_{E}\right>$ in terms of it as

$\displaystyle\left<\gamma({\bf k_{1}})\zeta_{E}({\bf
k_{2}})\right>=(2\pi)^{3}\delta^{(3)}(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2})\;e^{\gamma}_{NL}\left(\frac{\lambda(\tau_{I}){\cal
P}_{E}(k_{2})}{3H^{2}\epsilon}\right)P_{\gamma}(k_{1})\;F(N_{0},N_{2},\alpha_{E})~{}.$
(107)

In the regime $n>0$, say, for $n=2$, the magnetic field spectrum is scale
invariant, i.e., $P_{B}\sim k^{-3}$ while the electric field spectrum decays
so the correlator $\left<\gamma\zeta_{B}\right>$ would be dominant over
$\left<\gamma\zeta_{E}\right>$. However, in the negative $n$ regime, say, for
$n=-3$, $P_{B}$ is again scale invariant but $P_{E}$ grows and as a result,
$\left<\gamma\zeta_{E}\right>$ would dominate over
$\left<\gamma\zeta_{B}\right>$. In any case, $\left<\gamma\zeta\right>$ will
always receive a contribution from one of these correlators and thus, will be
non-vanishing.

In order to corroborate these results with observations, one has to rewrite
them in terms of the energy densities of the magnetic and electric fields,
defined by $\rho_{B}=3H^{2}\Omega_{B}$ and $\rho_{E}=3H^{2}\Omega_{E}$,
respectively. With this, we can now express $\Omega_{Y}$ as

$\displaystyle\Omega_{Y}(\tau_{I})=\left(\frac{\lambda(\tau_{I}){\cal
P}_{Y}(k_{2},\tau_{I})}{6H^{2}}\right)\left(\frac{1-e^{-N_{0}\alpha_{Y}}}{\alpha_{Y}}\right)e^{N_{2}\alpha_{Y}}~{},$
(108)

which allows us to write our results as

$\displaystyle\frac{\left\langle\gamma({\bf k_{1}})\zeta_{B}({\bf
k_{2}})\right>}{\left\langle\gamma({\bf k_{1}})\gamma({\bf k_{2}})\right>}$
$\displaystyle=$ $\displaystyle
b^{\gamma}_{NL}\left(\frac{\Omega_{B}(\tau_{I})}{\epsilon}\right)\;\left(\frac{2}{\alpha_{B}}-\frac{2N_{0}}{1-e^{-N_{0}\alpha_{B}}}\right)~{},$
(109) $\displaystyle\frac{\left\langle\gamma({\bf k_{1}})\zeta_{E}({\bf
k_{2}})\right>}{\left\langle\gamma({\bf k_{1}})\gamma({\bf k_{2}})\right>}$
$\displaystyle=$ $\displaystyle
e^{\gamma}_{NL}\left(\frac{\Omega_{E}(\tau_{I})}{\epsilon}\right)\;\left(\frac{2}{\alpha_{E}}-\frac{2N_{0}}{1-e^{-N_{0}\alpha_{E}}}\right)~{}.$
(110)

From these results, it is interesting to note that
$\left\langle\gamma\zeta_{B}\right>$ vanishes for a scale invariant magnetic
field spectrum defined by $n=2,-3$ such that $\alpha_{B}=0$. Similarly, the
correlator $\left\langle\gamma\zeta_{E}\right>$ vanishes for a scale invariant
electric field spectrum defined by $n=3,-2$ such that $\alpha_{E}=0$. The
energy density of magnetic fields at the end of inflation can be related to
the amplitude of magnetic fields today, as discussed in [39], by

$\displaystyle\Omega_{B}(\tau_{I})\sim
10^{-7}\left(\frac{B_{\text{today}}}{10^{-9}\text{G}}\right)^{2},$ (111)

which is obtained by using the fact the energy density of magnetic fields can
only be a very small fraction of the CMB energy density. Phenomenologically,
it is motivated to look at the regime in which the magnetic fields dominate.
For the relevant choice $-1/2<n\ll 2$, we obtain

$\displaystyle\bigg{|}\frac{\left\langle\gamma\zeta_{B}\right>}{\left\langle\gamma\gamma\right>}\bigg{|}\sim
2\times
10^{-5}\;b^{\gamma}_{NL}\left(\frac{B_{\text{today}}}{10^{-9}\text{G}}\right)^{2}\left(\frac{0.01}{\epsilon}\right)N_{0}.$
(112)

Similarly, for the same range of $n$, we can also express
$\left\langle\gamma\zeta_{E}\right>$ by defining a parameter
$\delta=\rho_{E}/\rho_{B}$, as

$\displaystyle\bigg{|}\frac{\left\langle\gamma\zeta_{E}\right>}{\left\langle\gamma\gamma\right>}\bigg{|}\sim
2\times
10^{-5}\;\delta\;e^{\gamma}_{NL}\left(\frac{B_{\text{today}}}{10^{-9}\text{G}}\right)^{2}\left(\frac{0.01}{\epsilon}\right)N_{0}~{}.$
(113)

At later epochs when the conductivity of the universe becomes large, the
electric fields are supposed to vanish and magnetic fields do evolve
adiabatically until today. In such a case, the ratio $\delta$ should be an
extremely small quantity and thus,
$\left\langle\gamma\zeta\right>\sim\left\langle\gamma\zeta_{B}\right>$. The
above result can be used to estimate the strength of
$\left\langle\gamma\zeta_{B}\right>$ if we assume that the long wavelength
mode associated with the graviton is at the CMB scales. Finally, for a nG
magnetic field strength today and with $\epsilon\sim 0.01$, $N_{0}\sim
10^{2}$, $b^{\gamma}_{NL}\sim{\cal O}(1)$, we find that
$\left\langle\gamma\zeta\right>\sim 10^{-3}\left\langle\gamma\gamma\right>$.
This indicates that the contribution of this novel correlator to the
cosmological observables remains sufficiently small for the case wherein the
magnetic fields are generated after the horizon exit of the graviton mode.
However, if the magnetic fields were generated before the horizon exit of the
graviton mode, one needs to use the full bispectrum obtained from the in-in
formalism to estimate the extent of this novel correlator. In addition, from
fig. (2), we get an indication that the result can be a few orders of
magnitude larger away from the squeezed limit. It would be interesting to
examine this result for a general case more carefully in the future.

##  6 Conclusions and discussions

Higher order cosmological correlations generally provide new insights into the
non-trivial interactions in the early universe. This paper has calculated a
new set of such cosmological correlation functions and their corresponding
semi-classical consistency relations. More precisely, we have computed the
non-Gaussian correlation functions of primordial gravitons produced during
inflation, with gauge fields that are non-minimally coupled in the early
universe. Our correlation functions are derived, both in terms of the vector
potential, corresponding to the correlation of gravitons with a vector boson,
and in terms of the associated electric and magnetic fields for the $U(1)$
gauge boson. We have obtained the full general results for these cross-
correlations for the first time in the literature by taking into account the
correct time evolution operator in the interaction picture of the in-in
formalism. Moreover, we have also shown the presence of the leading order
correction terms in all the correlators that had not been noticed earlier. In
the squeezed limit, we found that the non-linearity parameters remain small
but can be large for the other triangular shapes of the bispectra in Fourier
space.

In order to check our results, we have derived new semi-classical relations
analogous to known semi-classical relations in the literature (also sometimes
called the consistency relations or soft theorems, etc.). Some well-known
examples are the Maldacena consistency relations [33] for the three-point
correlation functions, the SSV relations [37] for the four-point functions,
and the GS relations for the infrared part of loop diagrams [38]. These
relations are generally valid for curvature perturbation modes and graviton
modes, but they are not directly translatable to the present case. In the
minimally coupled case, the gauge field sector is invariant under the
conformal transformations and therefore does only feel the expansion of
spacetime through trivial conformal rescaling of the flat spacetime results.
Thus, it is evident that the non-trivial contribution to the correlation
functions arises due to the time dependence of the non-minimal coupling that
breaks the conformal invariance. Therefore all the correlation functions
computed here are proportional to the time derivative of the non-minimal
coupling function. However, a suitable redefinition of the pump field in terms
of the non-minimal coupling function and a gauge field redefinition allows us
to write the action for the polarization modes of the gauge fields in the
identical form to a scalar field in a fictitious expanding universe, with an
expansion given by the time dependence of the non-minimal coupling. For a non-
minimal coupling corresponding to a fictitious inflationary universe, we show
that our new semi-classical relations correctly recover the squeezed limit of
the full in-in calculations. We checked the agreement with our new consistency
relations both for the vector field itself, and also in terms of the
corresponding electric and magnetic fields.

Finally, we have also discussed an interesting possible observational
consequence of our scenario. In particular, it is worthwhile keeping in mind,
that even if the gauge field is not directly observable today, as in the case
of dark sector gauge fields, it may be indirectly observable through its
imprints on the curvature perturbations. We then calculated the novel
correlation between one graviton mode and one curvature perturbation mode,
induced by the three-point correlation function of one graviton mode and two
gauge field modes. We stress that such a two-point correlator vanishes in
single field inflationary models but is non-trivial in scenarios in which
gauge fields are simultaneously excited during inflation due to their non-
conformal coupling. It would be interesting to study the observational
consequences of such correlators in terms of non-zero cross-correlations
involving temperature and polarisation anisotropies in the CMB and understand
the prospects of their detectability with both ground and space based upcoming
CMB experiments.

## Acknowledgments

RKJ would like to acknowledge financial support from the new faculty seed
start-up grant of IISc, the Core Research Grant CRG/2018/002200 from the
Science and Engineering Research Board, Department of Science and Technology,
Government of India and the Infosys Foundation, Bangalore through the Infosys
Young Investigator award. MSS is supported by Villum Fonden grant 13384 and
Independent Research Fund Denmark grant 0135-00378B.

##  Appendix A Useful integrals

For $n>-\frac{1}{2}$, the integrals $\tilde{\mathcal{I}}_{n}^{(1)}$ and
$\tilde{\mathcal{I}}_{n}^{(2)}$ can be explicitly written as,

$\displaystyle\tilde{\mathcal{I}}_{n}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{\pi^{3}}{2}\frac{2^{-2n-1}}{\Gamma^{2}(n+1/2)}(-k_{2}\tau_{I})^{n+1/2}(-k_{3}\tau_{I})^{n+1/2}$
(114) $\displaystyle\times$ $\displaystyle{\rm
Im}\biggl{[}(1+ik_{1}\tau_{I})e^{-ik_{1}\tau_{I}}H_{n+1/2}^{(1)}(-k_{2}\tau_{I})H_{n+1/2}^{(1)}(-k_{3}\tau_{I})\Biggr{.}$
$\displaystyle\times$
$\displaystyle\biggl{.}\int^{\tau_{I}}d\tau\tau(1-ik_{1}\tau)e^{ik_{1}\tau}H_{n-1/2}^{(2)}(-k_{2}\tau)H_{n-1/2}^{(2)}(-k_{3}\tau)\biggr{]}~{},$
$\displaystyle\tilde{\mathcal{I}}_{n}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{\pi^{3}}{2}\frac{2^{-2n-1}}{\Gamma^{2}(n+1/2)}(-k_{2}\tau_{I})^{n+1/2}(-k_{3}\tau_{I})^{n+1/2}$
(115) $\displaystyle\times$ $\displaystyle{\rm
Im}\biggl{[}(1+ik_{1}\tau_{I})e^{-ik_{1}\tau_{I}}H_{n+1/2}^{(1)}(-k_{2}\tau_{I})H_{n+1/2}^{(1)}(-k_{3}\tau_{I})\biggr{.}$
$\displaystyle\times$
$\displaystyle\biggl{.}\int^{\tau_{I}}d\tau\tau(1-ik_{1}\tau)e^{ik_{1}\tau}H_{n+1/2}^{(2)}(-k_{2}\tau)H_{n+1/2}^{(2)}(-k_{3}\tau)\biggr{]}~{}.$

Similarly for $n<-\frac{1}{2}$,

$\displaystyle\tilde{\mathcal{I}}_{n}^{(1)}$ $\displaystyle=$
$\displaystyle\pi\;\frac{2^{2n}\;\Gamma^{2}(n+3/2)}{(-k_{2}\tau_{I})^{n+1/2}(-k_{3}\tau_{I})^{n+1/2}}$
(116) $\displaystyle\times$ $\displaystyle{\rm
Im}\biggl{[}(1+ik_{1}\tau_{I})e^{-ik_{1}\tau_{I}}H_{n+1/2}^{(1)}(-k_{2}\tau_{I})H_{n+1/2}^{(1)}(-k_{3}\tau_{I})\Biggr{.}$
$\displaystyle\times$
$\displaystyle\biggl{.}\int^{\tau_{I}}d\tau\tau(1-ik_{1}\tau)e^{ik_{1}\tau}H_{n-1/2}^{(2)}(-k_{2}\tau)H_{n-1/2}^{(2)}(-k_{3}\tau)\biggr{]}~{},$
$\displaystyle\tilde{\mathcal{I}}_{n}^{(2)}$ $\displaystyle=$
$\displaystyle\pi\;\frac{2^{2n}\;\Gamma^{2}(n+3/2)}{(-k_{2}\tau_{I})^{n+1/2}(-k_{3}\tau_{I})^{n+1/2}}$
(117) $\displaystyle\times$ $\displaystyle{\rm
Im}\biggl{[}(1+ik_{1}\tau_{I})e^{-ik_{1}\tau_{I}}H_{n+1/2}^{(1)}(-k_{2}\tau_{I})H_{n+1/2}^{(1)}(-k_{3}\tau_{I})\biggr{.}$
$\displaystyle\times$
$\displaystyle\biggl{.}\int^{\tau_{I}}d\tau\tau(1-ik_{1}\tau)e^{ik_{1}\tau}H_{n+1/2}^{(2)}(-k_{2}\tau)H_{n+1/2}^{(2)}(-k_{3}\tau)\biggr{]}~{}.$

In the limit corresponding to the end of inflation $|k\tau_{I}|\to 0$, one
could show that
$\tilde{\mathcal{J}}_{n}^{(1)}\sim\tilde{\mathcal{I}}_{n}^{(1)}$ and
$\tilde{\mathcal{J}}_{n}^{(2)}\sim\tilde{\mathcal{I}}_{n}^{(2)}$. To achieve
this result, let’s rewrite the integrals (58) and (59) as follows,

$\displaystyle\mathcal{J}_{1}$ $\displaystyle=$ $\displaystyle 2\,{\rm
Im}\Biggl{[}\gamma_{k_{1}}(\tau_{I})A_{k_{2}}(\tau_{I})A_{k_{3}}(\tau_{I})\left(\frac{A^{\prime}_{k_{2}}(\tau_{I})}{A_{k_{2}}(\tau_{I})}\right)\left(\frac{A^{\prime}_{k_{3}}(\tau_{I})}{A_{k_{3}}(\tau_{I})}\right)$
(118) $\displaystyle\times$ $\displaystyle\int
d\tau\lambda(\tau)\gamma_{k_{1}}^{*}(\tau)A_{k_{2}}^{{}^{\prime}*}(\tau)A_{k_{3}}^{{}^{\prime}*}(\tau)\Biggr{]},$
$\displaystyle\mathcal{J}_{2}$ $\displaystyle=$ $\displaystyle 2\,{\rm
Im}\Biggl{[}\gamma_{k_{1}}(\tau_{I})A_{k_{2}}(\tau_{I})A_{k_{3}}(\tau_{I})\left(\frac{A^{\prime}_{k_{2}}(\tau_{I})}{A_{k_{2}}(\tau_{I})}\right)\left(\frac{A^{\prime}_{k_{3}}(\tau_{I})}{A_{k_{3}}(\tau_{I})}\right)$
(119) $\displaystyle\times$ $\displaystyle\int
d\tau\lambda(\tau)\gamma_{k_{1}}^{*}(\tau)A_{k_{2}}^{*}(\tau)A_{k_{3}}^{*}(\tau)\Biggr{]}.$

Using equation (22), the time derivative of mode function can be obtained,

$\displaystyle
A^{\prime}_{k}(\tau)=-kA_{k}(\tau)\;\left[\frac{H^{(1)}_{n-\frac{1}{2}}(-k\tau)}{H^{(1)}_{n+\frac{1}{2}}(-k\tau)}\right]~{}.$
(120)

By definition, i.e. from (39), (40), (60) and (61), we can see that,

$\displaystyle\frac{{\mathcal{J}}_{1}}{{\mathcal{I}}_{1}}$ $\displaystyle=$
$\displaystyle\frac{|A^{\prime(0)}_{k_{2}}(\tau_{I})|}{|A^{(0)}_{k_{2}}(\tau_{I})|}\frac{|A^{\prime(0)}_{k_{3}}(\tau_{I})|}{|A^{(0)}_{k_{3}}(\tau_{I})|}\;\left(\frac{\tilde{\mathcal{J}}_{n}^{(1)}}{\tilde{\mathcal{I}}_{n}^{(1)}}\right)$
(121) $\displaystyle\frac{{\mathcal{J}}_{2}}{{\mathcal{I}}_{2}}$
$\displaystyle=$
$\displaystyle\frac{|A^{\prime(0)}_{k_{2}}(\tau_{I})|}{|A^{(0)}_{k_{2}}(\tau_{I})|}\frac{|A^{\prime(0)}_{k_{3}}(\tau_{I})|}{|A^{(0)}_{k_{3}}(\tau_{I})|}\;\left(\frac{\tilde{\mathcal{J}}_{n}^{(2)}}{\tilde{\mathcal{I}}_{n}^{(2)}}\right).$
(122)

Now, let us take the limit $|k\tau_{I}|\to 0$ and consider the equations
(116), (120) and (121), we will obtain,

$\tilde{\mathcal{J}}_{n}^{(1)}\sim\tilde{\mathcal{I}}_{n}^{(1)}.$ (123)

Similarly by considering the equations (117), (120) and (122),

$\tilde{\mathcal{J}}_{n}^{(2)}\sim\tilde{\mathcal{I}}_{n}^{(2)}.$ (124)

One can evaluate these integrals for different values of $n$ which are listed
in the appendix of [17]. For completeness, we present the result below for the
most interesting case of $n=2$ which has been used to plot the non-linearity
parameters in Sec. 3.2, as

$\displaystyle\tilde{\mathcal{I}}_{2}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{-1}{(k_{2}k_{3})^{3/2}k_{t}^{2}}$ (125)
$\displaystyle\times$
$\displaystyle\left[-k_{1}^{3}-2k_{1}^{2}(k_{2}+k_{3})-2k_{1}(k_{2}^{2}+k_{2}k_{3}+k_{3}^{2})-(k_{2}+k_{3})(k_{2}^{2}+k_{2}k_{3}+k_{3}^{2})\right]$

and

$\displaystyle\tilde{\mathcal{I}}_{2}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{-1}{(k_{2}k_{3})^{5/2}k_{t}^{2}}$ (126)
$\displaystyle\times$
$\displaystyle\left[(k_{1}+k_{2})^{2}(-3k_{1}^{3}-3k_{1}^{2}k_{2}-k_{2}^{3})+(k_{1}+k_{2})(-9k_{1}^{3}-6k_{1}^{2}k_{2}-2k_{2}^{3})k_{3}\right.$
$\displaystyle+$
$\displaystyle(-9k_{1}^{3}-6k_{1}^{2}k_{2}-2k_{1}k_{2}^{2}-2k_{2}^{3})k_{3}^{2}$
$\displaystyle-$
$\displaystyle\left.2(2k_{1}^{2}+k_{1}k_{2}+k_{2}^{2})k_{3}^{3}-2(k_{1}+k_{2})k_{3}^{4}-k_{3}^{5}+3k_{1}^{3}k_{t}^{2}(\gamma+\ln(-k_{t}\tau_{I}))\right]$

where $k_{t}=k_{1}+k_{2}+k_{3}$ and $\gamma$ is the Euler gamma constant.

##  Appendix B Validity of the semi-classical approach

Here, we shall quickly review the validity of the semi-classical approach for
the three-point correlators involving tensor perturbation. For the semi-
classical derivation of single field consistency relations, one argues that
the effect of the long wavelength metric perturbation can be absorbed into a
local coordinate transformation. Then the two point functions in a modified
background can be calculated by rescaling the coordinates. And use them to
calculate the squeezed limit three-point functions. We can see that the same
arguments will be valid also in our scenario. So let’s try to establish this
fact at the action level. To start with, the gauge field action of our
interest is given by

$\displaystyle S_{EM}=-\frac{1}{4}\int
d^{4}x\sqrt{-g}\,\lambda(\phi)F_{\mu\nu}F^{\mu\nu}.$ (127)

Now we can calculate the leading order action by introducing tensor
perturbation,

$\displaystyle S=S^{(0)}-\frac{1}{4}\int
d^{4}x\lambda(\phi)\left[2F_{0i}F_{0j}\gamma^{ij}-2F_{ij}F_{im}\gamma^{mj}\right],$
(128)

where $\gamma_{ij}$ is the tensor perturbation and $S^{(0)}$ denotes the free
part of the action. So the second term in RHS of equation (128) corresponds to
interacting part of the action. One can express this interacting part in the
Fourier space and study its squeezed limit. In the squeezed limit, wherein the
momentum corresponding to $\gamma_{ij}$ is much smaller than the gauge fields,
we can write,

$\displaystyle S=S^{(0)}-\frac{1}{2}\int d\tau\lambda(\tau)\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\Big{[}A^{\prime}_{i}({\bf k})A^{\prime}_{j}(-{\bf
k})\gamma_{ij}-k^{2}A_{i}({\bf k})A_{j}(-{\bf
k})\gamma_{ij}-k_{i}k_{j}\gamma_{ij}A^{2}\Big{]}~{}.$ (129)

Since $\gamma_{ij}$ is time independent in the squeezed limit, we can use the
on-shell conditions to simplify equation (129) as follows,

$\displaystyle\int d\tau A^{\prime}_{i}({\bf k})A^{\prime}_{j}(-{\bf k})$
$\displaystyle=$ $\displaystyle\int
d\tau\left(\frac{\partial}{\partial\tau}(A_{i}\lambda(\tau)A^{\prime}_{j})-A_{i}\lambda(\tau)A^{\prime\prime}_{j}-A_{i}\lambda^{\prime}(\tau)A^{\prime}_{j}\right)$
(130) $\displaystyle=$ $\displaystyle-\int
d\tau\lambda(\tau)A_{i}\left(A^{\prime\prime}_{j}+\frac{\lambda^{\prime}}{\lambda}A^{\prime}_{j}\right)$
$\displaystyle=$ $\displaystyle\int d\tau\lambda(\tau)k^{2}A_{i}({\bf
k})A_{j}(-{\bf k})~{},$

where we have used the equation of motion of $A_{i}$ in the last line of the
above equation. Substitute equation (130) into (129) and simplify,

$\displaystyle S=S^{(0)}+\frac{1}{2}\int d\tau\lambda(\tau)\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}k_{i}k_{j}\gamma_{ij}A^{2}~{}.$ (131)

$S^{(0)}$ can be rewritten as,

$\displaystyle S^{(0)}=\frac{1}{2}\int d\tau\lambda(\tau)\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\left((A^{\prime}_{i})^{2}-k^{2}A^{2}\right)~{}.$ (132)

Then the total action becomes,

$\displaystyle S$ $\displaystyle=$ $\displaystyle\frac{1}{2}\int
d\tau\lambda(\tau)\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\left((A^{\prime}_{i})^{2}-(k^{2}-k_{i}k_{j}\gamma_{ij})A^{2}\right)$
(133) $\displaystyle=$ $\displaystyle\frac{1}{2}\int
d\tau\lambda(\tau)\int\frac{d^{3}{\bf
k}}{(2\pi)^{3}}\left((A^{\prime}_{i})^{2}-\tilde{k}^{2}A^{2}\right)~{},$

where $\tilde{k}^{2}=k^{2}-k_{i}k_{j}\gamma_{ij}$. One can identify this as
the rescaled momentum in the rescaled background, i.e.,
$ds^{2}=-dt^{2}+a^{2}(t)d\tilde{x}^{2}$ with $d\tilde{x}^{2}\to
dx^{2}+\gamma^{B}_{ij}dx^{i}dx^{j}$. So this establishes the fact that the
effect of the long wavelength graviton mode can be absorbed into a local
coordinate transformation.

##  Appendix C The in-in formalism and the dynamical correction term

The electric field is defined through the time-derivative of the vector
potential. When promoting these to operators and taking expectation values, we
need to be careful and distinguish between time-derivatives of operators in
the Heisenberg picture and the interaction picture. The relation between the
electric field and the time-derivative of the vector potential is simplest in
the Heisenberg picture where only the operators depend on time. We will
therefore start by considering the time derivative of an operator in the
Heisenberg picture, and then translate that into the corresponding operator in
the interaction picture where also the state depends on time. In this appendix
we will carefully distinguish between interaction picture and Heisenberg
picture operators, but in the main text all discussion is in the interaction
picture and all operators are implicitly understood to be in the interaction
picture unless explicitly defined otherwise.

Let $\mathcal{O}(\tau)$ be an observable of our interest which is constructed
out of field operators and their derivatives. The time evolution of this
observable is given by the Heisenberg equation of motion

$\displaystyle\frac{d}{d\tau}\mathcal{O}^{H}(\tau)=i[H_{\rm
tot}(\tau),\mathcal{O}^{H}(\tau)].$ (134)

One can use the standard procedure in quantum field theory and express the
solution in terms of a unitary operator $U$ as

$\displaystyle\mathcal{O}^{H}(\tau)=U^{-1}(\tau,\tau_{i})\mathcal{O}^{H}(\tau_{i})U(\tau,\tau_{i}),$
(135)

where $\tau_{i}$ denotes the initial time and $U$ obeys the following equation

$\displaystyle i\frac{\partial}{\partial\tau}U(\tau,\tau_{i})=H_{\rm
tot}(\tau)U(\tau,\tau_{i}),$ (136)

and satisfies the following properties

$\displaystyle U(\tau,\tau)=1,\qquad
U^{-1}(\tau,\tau_{i})U(\tau,\tau_{i})=1,\qquad
U(\tau_{3},\tau_{2})U(\tau_{2},\tau_{1})=U(\tau_{3},\tau_{1}).$ (137)

Thus, $U$ can be interpreted as a time evolution operator which evolves the
operator $\mathcal{O}$ from an initial time $\tau_{i}$ to a later time $\tau$.

In this work, we are interested in the expectation value of
$\mathcal{O}(\tau)$, i.e.,

$\displaystyle\left\langle\mathcal{O}(\tau)\right>\equiv\left\langle
in\right|\mathcal{O}^{H}(\tau)\left|in\right>=\left\langle
in\right|U^{-1}(\tau,\tau_{i})\mathcal{O}^{H}(\tau_{i})U(\tau,\tau_{i})\left|in\right>,$
(138)

where $|in\rangle$ denotes the initial quantum state of the system. Since we
are working in the Heisenberg picture, the states are independent of time.
Now, our main goal is to calculate equation (138). In order to do so, one has
to split the total Hamiltonian as follows

$\displaystyle H_{\rm tot}(\tau)=H_{0}(\tau)+H_{\text{int}}(\tau),$ (139)

where $H_{0}(\tau)$ and $H_{\text{int}}(\tau)$ denote the free and the
interacting parts of the Hamiltonian, respectively. Since we know the
evolution of the system under free Hamiltonian, i.e. without any interaction,
we can make use of it to understand the time evolution under full Hamiltonian.
This can be easily illustrated by introducing the interaction picture. To this
end, let us define $U_{0}(\tau,\tau_{i})$ as

$\displaystyle
i\frac{\partial}{\partial\tau}U_{0}(\tau,\tau_{i})=H_{0}(\tau)U_{0}(\tau,\tau_{i}),$
(140)

so one can rewrite the expectation value of $\mathcal{O}(\tau)$ as in equation
(138) as follows

$\displaystyle\left\langle
in\right|U^{-1}(\tau,\tau_{i})\mathcal{O}^{H}(\tau_{i})U(\tau,\tau_{i})\left|in\right>=\left\langle
in\right|F^{-1}(\tau,\tau_{i})U_{0}^{-1}(\tau,\tau_{i})\mathcal{O}^{H}(\tau_{i})U_{0}(\tau,\tau_{i})F(\tau,\tau_{i})\left|in\right>,$
(141)

with

$\displaystyle F(\tau,\tau_{i})=U_{0}^{-1}(\tau,\tau_{i})U(\tau,\tau_{i}).$
(142)

Then, we can define the operator $\mathcal{O}$ in the interaction picture in
the following manner as

$\displaystyle\mathcal{O}^{I}(\tau)=U_{0}^{-1}(\tau,\tau_{i})\mathcal{O}^{H}(\tau_{i})U_{0}(\tau,\tau_{i}).$
(143)

Using this definition, we can now express the expectation value of
$\mathcal{O}(\tau)$ as

$\displaystyle\left\langle\mathcal{O}(\tau)\right>=\left\langle
in\right|F^{-1}(\tau,\tau_{i})\mathcal{O}^{I}(\tau)F(\tau,\tau_{i})\left|in\right>,$
(144)

with

$\displaystyle
i\frac{\partial}{\partial\tau}F(\tau,\tau_{i})=H^{I}_{\text{int}}(\tau)F(\tau,\tau_{i}),$
(145)

and the superscript $I$ refers to the interaction picture. The solution to the
above equation can be written in terms of $H^{I}_{\text{int}}$ as

$\displaystyle
F(\tau,\tau_{i})=T\left(\text{exp}\left[-i\int_{\tau_{i}}^{\tau}d\tau^{\prime}H^{I}_{\text{int}}(\tau^{\prime})\right]\right),$
(146)

where $T$ is the time ordering operator. Using this equation, one can
calculate the expectation value of an operator $\mathcal{O}$ at time $\tau$,

$\displaystyle\left\langle\mathcal{O}(\tau)\right>=\left\langle
in\right|\bar{T}\left(e^{i\int_{\tau_{i}}^{\tau}d\tau^{\prime}H^{I}_{\text{int}}(\tau^{\prime})}\right)\mathcal{O}^{I}(\tau)T\left(e^{-i\int_{\tau_{i}}^{\tau}d\tau^{\prime}H^{I}_{\text{int}}(\tau^{\prime})}\right)\left|in\right>.$
(147)

This result is the master formula of the in-in formalism for calculating equal
time correlation functions. It was first developed by Schwinger [54] and
others [55, 56, 57, 58] and then later applied to cosmology [59, 60].

In this paper, we are interested in the correlators $\gamma A_{\mu}A^{\mu}$,
$\gamma B_{\mu}B^{\mu}$ and $\gamma E_{\mu}E^{\mu}$. The construction of
interaction picture of $\gamma A_{\mu}A^{\mu}$ and $\gamma B_{\mu}B^{\mu}$ is
trivial but $\gamma E_{\mu}E^{\mu}$ requires more attention, as it involves
the time derivative of the vector potential. Let us therefore compare the
time-derivative of the operator in the Heisenberg picture and in the
interaction picture.

The time derivative of the field has a clear meaning as only the operators
evolve in time in the Heisenberg picture,

$\frac{d}{d\tau}\left\langle\mathcal{O}(\tau)\right>=\left\langle
in\right|\frac{d}{d\tau}\mathcal{O}^{H}(\tau)\left|in\right>~{}.$ (148)

Now, identifying $\mathcal{Q}^{H}(\tau)\equiv(d/d\tau)\mathcal{O}^{H}(\tau)$
and using (143), we have
$\mathcal{Q}^{I}(\tau)=U_{0}^{-1}(\tau,\tau_{i})\mathcal{Q}^{H}(\tau_{i})U_{0}(\tau,\tau_{i})$.
Using equation (134), we get

$\displaystyle\mathcal{Q}^{I}(\tau)=U_{0}^{-1}(\tau,\tau_{i})\left(\frac{d}{d\tau}\mathcal{O}^{H}(\tau)\right)\bigg{|}_{\tau=\tau_{i}}U_{0}(\tau,\tau_{i})$
$\displaystyle=$ $\displaystyle iU_{0}^{-1}(\tau,\tau_{i})[H_{\rm
tot}(\tau_{i}),\mathcal{O}^{H}(\tau_{i})]U_{0}(\tau,\tau_{i})$ (149)
$\displaystyle=$ $\displaystyle
iU_{0}^{-1}(\tau,\tau_{i})[H_{0}(\tau_{i})+H_{\rm
int}(\tau_{i}),\mathcal{O}^{H}(\tau_{i})]U_{0}(\tau,\tau_{i})$
$\displaystyle=$ $\displaystyle i[H_{0}(\tau),\mathcal{O}^{I}(\tau)]+i[H_{\rm
int}^{I}(\tau),\mathcal{O}^{I}(\tau)]$ $\displaystyle=$
$\displaystyle\frac{d}{d\tau}\mathcal{O}^{I}(\tau)+i[H_{\rm
int}^{I}(\tau),\mathcal{O}^{I}(\tau)].$

We conclude, that if we have an operator, which in the Heisenberg picture
would take the form of a time derivative of an operator, i.e.
$\mathcal{Q}^{H}(\tau)\equiv(d/d\tau)\mathcal{O}^{H}(\tau)$, then the
interaction picture version of that operator is the time derivative of the
interaction picture operator plus and extra piece, i.e.
$\mathcal{Q}^{I}(\tau)=(d/d\tau)\mathcal{O}^{I}(\tau)+i[H_{\rm
int}^{I}(\tau),\mathcal{O}^{I}(\tau)]$.

Applying this observation to the electric field in the interaction picture
$E_{\mu}^{I}$, we find

$\displaystyle
E^{I}_{i}(\tau)=\frac{1}{a(\tau)}\left[\left(\frac{d}{d\tau}A_{i}^{I}(\tau)\right)+i[H_{\rm
int}^{I}(\tau),A_{i}^{I}(\tau)]\right].$ (150)

Since all the operators in the above equation are in the interaction picture,
let’s suppress the corresponding superscript $I$ for brevity. Then one can
express the observable $\gamma E_{\mu}E^{\mu}$, to the leading order, in the
interaction picture as follows,

$\displaystyle\gamma
E_{\mu}E^{\mu}=\frac{1}{a^{4}}\left[\gamma\frac{dA_{i}}{d\tau}\frac{dA_{i}}{d\tau}+i\gamma\left(\frac{dA_{i}}{d\tau}[H_{\rm
int},A_{i}]+[H_{\rm
int},A_{i}]\frac{dA_{i}}{d\tau}\right)-\gamma\gamma_{ij}\frac{dA_{i}}{d\tau}\;\frac{dA_{j}}{d\tau}\right].$
(151)

We refer to the second term as the dynamical correction term and the third as
the kinematical correction term.

## References

  * [1] M. S. Turner and L. M. Widrow, Inflation Produced, Large Scale Magnetic Fields, Phys. Rev. D37 (1988) 2743.
  * [2] B. Ratra, Cosmological ’seed’ magnetic field from inflation, Astrophys. J. 391 (1992) L1–L4.
  * [3] F. Tavecchio et al., The intergalactic magnetic field constrained by Fermi/LAT observations of the TeV blazar 1ES 0229+200, Mon. Not. Roy. Astron. Soc. 406 (2010) L70–L74, [arXiv:1004.1329].
  * [4] A. Neronov and I. Vovk, Evidence for strong extragalactic magnetic fields from Fermi observations of TeV blazars, Science 328 (2010) 73–75, [arXiv:1006.3504].
  * [5] W. Chen, J. H. Buckley, and F. Ferrer, Search for GeV $\gamma$ -Ray Pair Halos Around Low Redshift Blazars, Phys. Rev. Lett. 115 (2015) 211103, [arXiv:1410.7717].
  * [6] L. M. Widrow, Origin of Galactic and Extragalactic Magnetic Fields, Rev. Mod. Phys. 74 (2002) 775–823, [astro-ph/0207240].
  * [7] V. Demozzi, V. Mukhanov, and H. Rubinstein, Magnetic fields from inflation?, JCAP 0908 (2009) 025, [arXiv:0907.1030].
  * [8] R. Durrer, L. Hollenstein, and R. K. Jain, Can slow roll inflation induce relevant helical magnetic fields?, JCAP 1103 (2011) 037, [arXiv:1005.5322].
  * [9] C. T. Byrnes, L. Hollenstein, R. K. Jain, and F. R. Urban, Resonant magnetic fields from inflation, JCAP 1203 (2012) 009, [arXiv:1111.2030].
  * [10] R. K. Jain, R. Durrer, and L. Hollenstein, Generation of helical magnetic fields from inflation, arXiv:1204.2409.
  * [11] R. J. Ferreira, R. K. Jain, and M. S. Sloth, Inflationary magnetogenesis without the strong coupling problem, JCAP 10 (2013) 004, [arXiv:1305.7151].
  * [12] R. J. Ferreira, R. K. Jain, and M. S. Sloth, Inflationary Magnetogenesis without the Strong Coupling Problem II: Constraints from CMB anisotropies and B-modes, JCAP 06 (2014) 053, [arXiv:1403.5516].
  * [13] T. Kobayashi and M. S. Sloth, Early Cosmological Evolution of Primordial Electromagnetic Fields, Phys. Rev. D 100 (2019), no. 2 023524, [arXiv:1903.02561].
  * [14] R. R. Caldwell, L. Motta, and M. Kamionkowski, Correlation of inflation-produced magnetic fields with scalar fluctuations, Phys. Rev. D84 (2011) 123525, [arXiv:1109.4415].
  * [15] L. Motta and R. R. Caldwell, Non-Gaussian features of primordial magnetic fields in power-law inflation, Phys. Rev. D85 (2012) 103532, [arXiv:1203.1033].
  * [16] R. K. Jain and M. S. Sloth, Consistency relation for cosmic magnetic fields, Phys. Rev. D 86 (2012) 123528, [arXiv:1207.4187].
  * [17] R. K. Jain and M. S. Sloth, On the non-Gaussian correlation of the primordial curvature perturbation with vector fields, JCAP 02 (2013) 003, [arXiv:1210.3461].
  * [18] M. Biagetti, A. Kehagias, E. Morgante, H. Perrier, and A. Riotto, Symmetries of Vector Perturbations during the de Sitter Epoch, JCAP 07 (2013) 030, [arXiv:1304.7785].
  * [19] D. Chowdhury, L. Sriramkumar, and M. Kamionkowski, Cross-correlations between scalar perturbations and magnetic fields in bouncing universes, JCAP 01 (2019) 048, [arXiv:1807.05530].
  * [20] D. Chowdhury, L. Sriramkumar, and M. Kamionkowski, Enhancing the cross-correlations between magnetic fields and scalar perturbations through parity violation, JCAP 10 (2018) 031, [arXiv:1807.07477].
  * [21] S. B. Giddings and M. S. Sloth, Cosmological diagrammatic rules, JCAP 1007 (2010) 015, [arXiv:1005.3287].
  * [22] S. Weinberg, Infrared photons and gravitons, Phys. Rev. 140 (1965) B516–B524.
  * [23] A. Strominger and A. Zhiboedov, Gravitational Memory, BMS Supertranslations and Soft Theorems, JHEP 01 (2016) 086, [arXiv:1411.5745].
  * [24] A. Strominger, Lectures on the Infrared Structure of Gravity and Gauge Theory, arXiv:1703.05448.
  * [25] R. Z. Ferreira, M. Sandora, and M. S. Sloth, Asymptotic Symmetries in de Sitter and Inflationary Spacetimes, JCAP 04 (2017) 033, [arXiv:1609.06318].
  * [26] R. Z. Ferreira, M. Sandora, and M. S. Sloth, Patient Observers and Non-perturbative Infrared Dynamics in Inflation, JCAP 02 (2018) 055, [arXiv:1703.10162].
  * [27] R. Z. Ferreira, M. Sandora, and M. S. Sloth, Is Patience a Virtue? Cosmic Censorship of Infrared Effects in de Sitter, Int. J. Mod. Phys. D 26 (2017), no. 12 1743019, [arXiv:1705.06380].
  * [28] L. Hui, A. Joyce, and S. S. C. Wong, Inflationary soft theorems revisited: A generalized consistency relation, JCAP 02 (2019) 060, [arXiv:1811.05951].
  * [29] A. Kehagias and A. Riotto, Inflation and Conformal Invariance: The Perspective from Radial Quantization, Fortsch. Phys. 65 (2017), no. 5 1700023, [arXiv:1701.05462].
  * [30] K. Hinterbichler, A. Joyce, and J. Khoury, Inflation in Flatland, JCAP 01 (2017) 044, [arXiv:1609.09497].
  * [31] K. Hinterbichler, L. Hui, and J. Khoury, Conformal Symmetries of Adiabatic Modes in Cosmology, JCAP 08 (2012) 017, [arXiv:1203.6351].
  * [32] K. Hinterbichler, L. Hui, and J. Khoury, An Infinite Set of Ward Identities for Adiabatic Modes in Cosmology, JCAP 01 (2014) 039, [arXiv:1304.5527].
  * [33] J. M. Maldacena, Non-Gaussian features of primordial fluctuations in single field inflationary models, JHEP 05 (2003) 013, [astro-ph/0210603].
  * [34] J. M. Maldacena and G. L. Pimentel, On graviton non-Gaussianities during inflation, JHEP 09 (2011) 045, [arXiv:1104.2846].
  * [35] P. Creminelli and M. Zaldarriaga, Single field consistency relation for the 3-point function, JCAP 0410 (2004) 006, [astro-ph/0407059].
  * [36] C. Cheung, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, On the consistency relation of the 3-point function in single field inflation, JCAP 02 (2008) 021, [arXiv:0709.0295].
  * [37] D. Seery, M. S. Sloth, and F. Vernizzi, Inflationary trispectrum from graviton exchange, JCAP 03 (2009) 018, [arXiv:0811.3934].
  * [38] S. B. Giddings and M. S. Sloth, Semiclassical relations and IR effects in de Sitter and slow-roll space-times, JCAP 1101 (2011) 023, [arXiv:1005.1056].
  * [39] S. Nurmi and M. S. Sloth, Constraints on Gauge Field Production during Inflation, JCAP 07 (2014) 012, [arXiv:1312.4946].
  * [40] S. Yokoyama and J. Soda, Primordial statistical anisotropy generated at the end of inflation, JCAP 0808 (2008) 005, [arXiv:0805.4265].
  * [41] N. Barnaby, R. Namba, and M. Peloso, Observable non-gaussianity from gauge field production in slow roll inflation, and a challenging connection with magnetogenesis, Phys. Rev. D85 (2012) 123523, [arXiv:1202.1469].
  * [42] J. Soda, Statistical Anisotropy from Anisotropic Inflation, Class. Quant. Grav. 29 (2012) 083001, [arXiv:1201.6434].
  * [43] N. Bartolo, S. Matarrese, M. Peloso, and A. Ricciardone, The anisotropic power spectrum and bispectrum in the $f(\phi)F^{2}$ mechanism, Phys.Rev. D87 (2013) 023504, [arXiv:1210.3257].
  * [44] D. H. Lyth and M. Karciauskas, The statistically anisotropic curvature perturbation generated by $f(\phi)^{2}F^{2}$, JCAP 05 (2013) 011, [arXiv:1302.7304].
  * [45] A. A. Abolhasani, R. Emami, J. T. Firouzjaee, and H. Firouzjahi, $\delta N$ formalism in anisotropic inflation and large anisotropic bispectrum and trispectrum, JCAP 08 (2013) 016, [arXiv:1302.6986].
  * [46] M. Shiraishi, E. Komatsu, M. Peloso, and N. Barnaby, Signatures of anisotropic sources in the squeezed-limit bispectrum of the cosmic microwave background, JCAP 05 (2013) 002, [arXiv:1302.3056].
  * [47] T. Fujita and S. Yokoyama, Higher order statistics of curvature perturbations in IFF model and its Planck constraints, JCAP 09 (2013) 009, [arXiv:1306.2992].
  * [48] J. P. Beltrán Almeida, J. Motoa-Manzano, and C. A. Valenzuela-Toledo, Correlation functions of sourced gravitational waves in inflationary scalar vector models. A symmetry based approach, JHEP 09 (2019) 118, [arXiv:1905.00900].
  * [49] J. Martin and J. Yokoyama, Generation of Large-Scale Magnetic Fields in Single-Field Inflation, JCAP 0801 (2008) 025, [arXiv:0711.4307].
  * [50] K. Subramanian, Magnetic fields in the early universe, Astron. Nachr. 331 (2010) 110–120, [arXiv:0911.4771].
  * [51] M. Shiraishi, S. Saga, and S. Yokoyama, CMB power spectra induced by primordial cross-bispectra between metric perturbations and vector fields, JCAP 1211 (2012) 046, [arXiv:1209.3384].
  * [52] M. Giovannini, Electric-magnetic duality and the conditions of inflationary magnetogenesis, JCAP 04 (2010) 003, [arXiv:0911.0896].
  * [53] S. B. Giddings and M. S. Sloth, Cosmological observables, IR growth of fluctuations, and scale-dependent anisotropies, Phys. Rev. D84 (2011) 063528, [arXiv:1104.0002].
  * [54] J. S. Schwinger, Brownian motion of a quantum oscillator, J. Math. Phys. 2 (1961) 407–432.
  * [55] K. T. Mahanthappa, Multiple production of photons in quantum electrodynamics, Phys. Rev. 126 (1962) 329–340.
  * [56] P. M. Bakshi and K. T. Mahanthappa, Expectation value formalism in quantum field theory. 2., J. Math. Phys. 4 (1963) 12–16.
  * [57] L. V. Keldysh, Diagram technique for nonequilibrium processes, Zh. Eksp. Teor. Fiz. 47 (1964) 1515–1527.
  * [58] K.-c. Chou, Z.-b. Su, B.-l. Hao, and L. Yu, Equilibrium and Nonequilibrium Formalisms Made Unified, Phys. Rept. 118 (1985) 1–131.
  * [59] R. D. Jordan, Effective Field Equations for Expectation Values, Phys. Rev. D 33 (1986) 444–454.
  * [60] E. Calzetta and B. L. Hu, Closed-time-path functional formalism in curved spacetime: Application to cosmological back-reaction problems, Phys. Rev. D 35 (Jan, 1987) 495–509.