# Probing Dark Matter Isocurvature with Primordial Non-Gaussianity

Michael Geller  Soubhik Kumar  Lian-Tao Wang

###### Abstract

Multiple fields can become dynamical during the inflationary epoch. We
consider an example where a light field acquires isocurvature fluctuations
during inflation and contributes to the dark matter abundance at late times.
Interactions between the light field and the adiabatic sector contribute to
mixed adiabatic-isocurvature non-Gaussianity (NG). We show the resulting form
of NG has a different kinematic dependence than the ‘local shape’ commonly
considered, and highlight the parameter space where a dedicated search is
expected to significantly improve the current Planck sensitivity. We interpret
our results in the context of the QCD axion and illustrate how the proposed NG
searches can improve upon the existing searches for isocurvature power
spectrum and bispectrum.

##  1 Introduction

Non-Gaussianities (NG) in primordial density perturbations can reveal the
nature of their dynamical origin and related physics (for reviews see [1, 2]).
The search for, and the non-observation of NG in the correlations of both
adiabatic and isocurvature perturbations [3] have already helped inform us in
early universe model building. Among the various measures of NG, the three-
point correlation function, i.e., the bispectrum, has been probed extensively,
especially using the Cosmic Microwave Background (CMB) and also Large-Scale
Structure (LSS) (e.g., [4, 5, 6]). Assuming translational and rotational
invariance, the primordial bispectrum depends only on the magnitudes
$k_{1},k_{2},k_{3}$ of the three fluctuation modes which form a triangle,
$\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3}=0$. Different inflationary
scenarios exhibit bispectra peaks for different ‘shapes’ of the triangles.

Among these shapes, a primary focus has been on the bispectrum with the so-
called ‘local shape’ [7]. This is particularly relevant when multiple fields
during inflation are present. In such scenarios initial isocurvature
fluctuations can convert into curvature fluctuations on superhorizon scales,
generating the local shape NG in the process. In particular, the curvaton
scenario [8, 9, 10] predict the local NG parameter $f_{\rm NL}^{\rm loc}\sim
1$ [11, 12, 13], a target for various upcoming surveys. Since the isocurvature
perturbations are constrained to be subdominant compared to the adiabatic
perturbations, the standard expectation is that the associated bispectrum,
involving only adiabatic perturbations, would be the leading discovery channel
for interesting dynamics.

In this work, we point out that bispectrum involving both adiabatic and
isocurvature perturbation can be indispensable in uncovering the fundamental
dynamics leading to NG. In particular, we consider an example in which the
adiabatic density perturbation is generated by a curvaton [8, 9, 10]. At the
same time, there is another light field, with an approximate shift symmetry,
that contributes to the dark matter (DM) abundance. The interaction between
the curvaton and the additional light field can come from a kinetic mixing,
which preserves the shift symmetry and is straightforward to UV complete. In
this case, we show that correlations involving two adiabatic and one DM
isocurvature mode provide the leading channel in observing the existence of
the additional light field and proving such a coupling.111This model also
predict sizable NG with more than one isocurvature perturbation. However, due
to the limited sensitivities in the observations, these are less useful.
Furthermore, the bispectrum involving only adiabatic fluctuations can also be
larger than the minimal expectation $f_{\rm NL}^{\rm loc}\sim 1$.

In addition, the predicted shape of NG in this class of models is not of the
local type as considered in [14, 15, 16, 17, 18, 19, 20, 21, 22], but rather
of the type considered in, e.g., [23, 24, 25] (in the context of adiabatic
NG), mediated by derivative interactions. We quantify the difference in shape
following the approach of [26]. We find the ‘cosine’ [26] between the shapes
employed in [3] and the ones studied here is $\lesssim 0.5$, indicating a
factor of few stronger sensitivity could be obtained by performing a dedicated
analysis just by using the existing Planck data. In addition, for certain
values of local NG parameters, the ‘cosine’ [26] is $\sim 0.1$, implying an
even better improvement would be feasible. This highlights the importance of
optimizing the template used in interpreting the relevant observations.

A natural application of the setup considered here is to a scenario in which
the light field is the QCD axion. We demonstrate that the bispectrum involving
two adiabatic and one isocurvature mode is indeed the leading observable for
the axion-curvaton interaction. In particular, it is more powerful than both
the imprint of isocurvature power spectrum and the isocurvature-only
bispectrum in revealing the existence of the axion sector. We find a dedicated
search would be able to probe novel regions of axion parameter space,
unconstrained by the current stellar cooling and CMB bounds.

The rest of the paper is organized as follows. In Sec. 2 we summarize the
shapes and the existing limits derived in the previous literature. This also
sets up the notation we use in Sec. 3 to describe the model. We compute the
predicted shapes and their overlap with the local shape to highlight the
importance of a dedicated analysis. We also interpret these results in the
context of the QCD axion. In Sec. 4 we comment on other model-building
possibilities. We also show how observable adiabatic-isocurvature NG can arise
in the more model-agnostic effective field theory of inflation framework. We
conclude in Sec. 5.

##  2 Notation and Existing Limits

In this section, we briefly summarize the existing searches for isocurvature
bispectrum carried out by the Planck collaboration [3]. This will also
establish the convention we use for the following discussion. The templates
used to search for isocurvature bispectrum is a generalized version of the
template for local NG:

$\begin{split}B^{IJK}(k_{1},k_{2},k_{3})&\equiv\langle\Phi^{I}(\mathbf{k}_{1})\Phi^{J}(\mathbf{k}_{2})\Phi^{K}(\mathbf{k}_{3})\rangle^{\prime}\\\
&=2f_{\rm NL}^{I,JK}P_{\Phi}(k_{2})P_{\Phi}(k_{3})+2f_{\rm
NL}^{J,KI}P_{\Phi}(k_{1})P_{\Phi}(k_{3})+2f_{\rm
NL}^{K,IJ}P_{\Phi}(k_{1})P_{\Phi}(k_{2}),\end{split}$ (2.1)

where $\Phi(\mathbf{k})$ is the Bardeen potential [27]. The symbol ′ denotes
that momentum-conserving delta functions have been taken out, i.e.,
$\langle\cdots\rangle\equiv\langle\cdots\rangle^{\prime}(2\pi)^{3}\delta(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})$.
The labels $I,J,K$ take values $a$ or $i$, and denote whether we have an
adiabatic mode or an isocurvature mode, respectively. The power spectrum is
denoted as
$P_{\Phi}(k)\equiv\langle\Phi(\mathbf{k})\Phi(-\mathbf{k})\rangle^{\prime}$
with $k=|\mathbf{k}|$. There is an exchange symmetry over the two indices
after the comma in $f_{\rm NL}^{I,JK}$. Therefore, for the parametrization in
(2.1), six possible quantities determine the strength of NG: $f_{\rm
NL}^{a,aa}$, $f_{\rm NL}^{a,ai}$, $f_{\rm NL}^{a,ii}$, $f_{\rm NL}^{i,aa}$,
$f_{\rm NL}^{i,ai}$, $f_{\rm NL}^{i,ii}$. The current observational constraint
on these for DM isocurvature, from Planck 2018 temperature and polarization
data, is given by [3]

$\begin{split}&f_{\rm NL}^{a,aa}=-2.5\pm 5.0,~{}~{}f_{\rm NL}^{a,ai}=-10\pm
10,~{}~{}f_{\rm NL}^{a,ii}=-450\pm 520,\\\ &f_{\rm NL}^{i,aa}=20\pm
28,~{}~{}f_{\rm NL}^{i,ai}=-32\pm 46,~{}~{}f_{\rm NL}^{i,ii}=-290\pm
210.\end{split}$ (2.2)

Here each of the six parameters is analyzed independently, i.e., it is assumed
that only one of these parameters is non-zero at a time. However, for the
scenario we describe later, all these parameters are expected to be non-zero.
A joint analysis, considering all the six parameters, is more appropriate and
the constraints are given by [3]

$\begin{split}&f_{\rm NL}^{a,aa}=4\pm 10,~{}~{}f_{\rm NL}^{a,ai}=-14\pm
21,~{}~{}f_{\rm NL}^{a,ii}=-3100\pm 1500,\\\ &f_{\rm NL}^{i,aa}=96\pm
52,~{}~{}f_{\rm NL}^{i,ai}=190\pm 180,~{}~{}f_{\rm NL}^{i,ii}=-640\pm
400.\end{split}$ (2.3)

To compare with theory predictions, it is convenient to relate the above
parameters in terms NG defined in terms of primordial curvature $\zeta$ and
isocurvature $S$ perturbations [19]:

$\begin{split}\tilde{B}^{ABC}(k_{1},k_{2},k_{3})&\equiv\langle
X^{A}(\mathbf{k}_{1})X^{B}(\mathbf{k}_{2})X^{C}(\mathbf{k}_{3})\rangle^{\prime}\\\
&=\tilde{f}_{\rm NL}^{A,BC}P_{\zeta}(k_{2})P_{\zeta}(k_{3})+\tilde{f}_{\rm
NL}^{B,CA}P_{\zeta}(k_{1})P_{\zeta}(k_{3})+\tilde{f}_{\rm
NL}^{C,AB}P_{\zeta}(k_{1})P_{\zeta}(k_{2}).\end{split}$ (2.4)

The indices $A,B,C$ can each be either $\zeta$ or $S$, and the fluctuation $X$
is accordingly determined, i.e., $X^{\zeta}\equiv\zeta$ and $X^{S}\equiv S$.
To relate $f_{\rm NL}$ with $\tilde{f}_{\rm NL}$, we use the relations
$\Phi=3\zeta/5$ and $\Phi=S/5$, valid during matter domination, for
superhorizon adiabatic and isocurvature modes, respectively. Furthermore, in
the Planck search, it is assumed that the power spectrum of $S$ has the same
spectral index as that of $\zeta$. Thus for the purpose of constraining
isocurvature NG, any difference in the magnitude of $P_{S}$ and $P_{\zeta}$ is
absorbed into the parameters $\tilde{f}_{\rm NL}$, and $P_{S}$ is taken to be
equal to $P_{\zeta}$. Using this, the numerical factors relating $f_{\rm NL}$
with $\tilde{f}_{\rm NL}$ can be derived [3]

$\begin{split}\tilde{f}_{\rm NL}^{\zeta,\zeta\zeta}=(6/5)f_{\rm
NL}^{a,aa},~{}\tilde{f}_{\rm NL}^{\zeta,\zeta S}=(2/5)f_{\rm
NL}^{a,ai},~{}\tilde{f}_{\rm NL}^{\zeta,SS}=(2/15)f_{\rm NL}^{i,aa},\\\
\tilde{f}_{\rm NL}^{S,\zeta\zeta}=(18/5)f_{\rm NL}^{i,aa},~{}\tilde{f}_{\rm
NL}^{S,\zeta S}=(6/5)f_{\rm NL}^{i,ai},~{}\tilde{f}_{\rm
NL}^{S,SS}=(2/5)f_{\rm NL}^{i,ii}.\end{split}$ (2.5)

In terms of these, the Planck constraints in (2.3) read:

$\begin{split}&\tilde{f}_{\rm NL}^{\zeta,\zeta\zeta}=4.8\pm
12,~{}~{}\tilde{f}_{\rm NL}^{\zeta,\zeta S}=-5.6\pm 8.4,~{}~{}\tilde{f}_{\rm
NL}^{\zeta,SS}=-413\pm 200,\\\ &\tilde{f}_{\rm NL}^{S,\zeta\zeta}=346\pm
187,~{}~{}\tilde{f}_{\rm NL}^{S,\zeta S}=228\pm 216,~{}~{}\tilde{f}_{\rm
NL}^{S,SS}=-256\pm 160.\end{split}$ (2.6)

Therefore, the current constraints on $\tilde{f}_{\rm NL}$ are ${\cal
O}(100)$, except for $\tilde{f}_{\rm NL}^{\zeta,\zeta\zeta}$ and
$\tilde{f}_{\rm NL}^{\zeta,\zeta S}$, for which the constraints are an order
of magnitude stronger. For the forecasted sensitivity of future surveys, such
as LiteBIRD and CMB-S4, see Ref. [22].

While certain multifield models predict the bispectrum shape dictated by
(2.4), it is not the only possibility. We will show in the next section that
in scenarios involving light, (pseudo)-Goldstone bosons, DM isocurvature NG
can naturally arise with a shape different from the local shape (2.4). Such
interactions are mediated by derivative couplings, and therefore the shape is
peaked at the equilateral triangle configurations $k_{1}\simeq k_{2}\simeq
k_{3}$. Correspondingly, the overlap between the equilateral shape and the
local shape in (2.4) is small. This implies that a search carried out using
(2.4) will not be optimal in detecting the actual equilateral shape. To have a
simple measure of this, the ‘cosine’ between two shapes can be computed [26].
First, a dot product is defined:

$\displaystyle B_{1}\cdot
B_{2}\equiv\sum_{\mathbf{k}_{i}}{\frac{\tilde{B}^{ABC}_{1}(k_{1},k_{2},k_{3})\tilde{B}^{ABC}_{2}(k_{1},k_{2},k_{3})}{\sigma_{k_{1}}^{2}\sigma_{k_{2}}^{2}\sigma_{k_{3}}^{2}}},$
(2.7)

with $\sigma_{k}^{2}$ being the variance of a given $k$-mode and the sum
running over all momentum triangle configurations. The cosine is then defined
as,

$\begin{split}\cos(B_{1},B_{2})\equiv\frac{B_{1}\cdot B_{2}}{(B_{1}\cdot
B_{1})^{1/2}(B_{2}\cdot B_{2})^{1/2}}.\end{split}$ (2.8)

A cosine value close to unity implies a template designed to detect shape
$B_{1}$ would also be optimal in detecting NG which has shape $B_{2}$, while a
cosine value closer to zero would imply the opposite. We will compute the
associated cosines in Sec. 3.

##  3 Model of Dark Matter Isocurvature

The cosmological evolution we consider is similar to the curvaton scenario [8,
9, 10]. While this is not the most minimal scenario, we argue in Sec. 4 why
the standard slow-roll inflationary scenario, in the absence of any curvaton
dynamics, gives a suppressed adiabatic-isocurvature NG. On the other hand, in
Sec. 4 we also show that the EFT of inflation framework [28] can give
observable NG, while being agnostic about the full homogeneous dynamics. From
that perspective, the model in this section can be considered a ‘UV
completion’ of the EFT description and the EFT cutoff scale can be explained
in terms of the fundamental UV parameters.

###  3.1 Setup

Consider an inflaton field $\phi$ that gives rise to the homogeneous spacetime
expansion during inflation. The fluctuations of $\phi$ are much smaller than
$\sim$$10^{-5}$, the typical size of CMB anisotropies. Instead, these
anisotropies originate from a ‘curvaton’ field $\sigma$ which behaves as a
spectator field during inflation, having a subdominant energy density compared
to $\phi$. In the post-inflationary universe, $\sigma$ eventually dominates
the total energy density, and fluctuations of $\sigma$ determine the
primordial density fluctuations. Eventually, $\sigma$ decays and gives rise to
the standard radiation-dominated universe with the fluctuations of $\sigma$
determining the CMB anisotropies. The dimensionless primordial scalar power
spectrum, originating from $\sigma$, is given by

$\begin{split}{\cal P}_{\zeta}(k)\equiv{k^{3}\over
2\pi^{2}}\langle\zeta(\mathbf{k})\zeta(\mathbf{-k})\rangle^{\prime}\equiv{k^{3}\over
2\pi^{2}}P_{\zeta}(k)={4\over 9}{H^{2}\over
4\pi^{2}\sigma_{0}^{2}},\end{split}$ (3.1)

where $H$ is the Hubble scale and $\sigma_{0}$ is the ‘misaligned’ field value
of $\sigma$ when the $k$-mode exits the horizon during inflation. The tilt is
given by,

$\begin{split}n_{s}-1=-2\epsilon+{2\over 3}\eta_{\sigma},\end{split}$ (3.2)

where $\epsilon=-\dot{H}/H^{2}$ and $\eta_{\sigma}=m_{\sigma}^{2}/H^{2}$ is
given by the mass of the curvaton $m_{\sigma}$. The tensor-to-scalar ratio
${\rm r}$ is determined by

$\begin{split}{\rm r}={8H^{2}/M_{\rm pl}^{2}\over
4H^{2}/(9\sigma_{0}^{2})}=18\sigma_{0}^{2}/M_{\rm pl}^{2}.\end{split}$ (3.3)

In particular, different from the single field scenario, ${\rm r}$ is not
determined by $\epsilon$, since the scalar fluctuations originate from a field
that is different from the one that drives the homogeneous expansion.

We now consider the effect of another scalar field $\chi$, that gives (a
subcomponent of) the DM abundance, on this cosmology. In this subsection, we
remain agnostic about the particle nature of $\chi$ and only impose that it is
a light field during inflation, obeying an approximate shift symmetry:
$\chi\rightarrow\chi+{\rm constant}$. Later in this section, we consider the
case where $\chi$ is the QCD axion [29, 30, 31, 32] that solves the strong CP
problem and also contributes to the DM abundance [33, 34, 35]. In both these
scenarios, we treat $\chi$ as effectively massless during inflation and it
only starts to oscillate when the Hubble scale falls below the mass of $\chi$,
after which $\chi$ abundance dilutes as matter. In our scenario, this happens
after $\sigma$ decay such that $\chi$ never dominates the energy density of
the universe before the matter-radiation inequality. The total curvature
perturbation $\zeta$ after inflation on superhorizon scales is given by,

$\begin{split}\zeta=f_{r}\zeta_{r}+f_{\sigma}\zeta_{\sigma}+f_{\chi}\zeta_{\chi},\end{split}$
(3.4)

such that $\sum_{i=r,\sigma,\chi}f_{i}=1$ with
$f_{i}=\dot{\rho}_{i}/\dot{\rho}_{\rm tot}$ and $\rho_{\rm
tot}=\sum_{i=r,\sigma,\chi}\rho_{i}$ is the total energy density. $\zeta_{i}$
denotes the individual gauge-invariant perturbation in each species, for a
review see [36]. The label $r$ denotes radiation decay products of the
inflaton. As in the curvaton paradigm, we assume
$\zeta_{r}\ll\zeta_{\sigma}\sim 10^{-5}$. Furthermore, the fractional energy
density in inflaton decay products, $f_{r}$, decreases with time since both
$\sigma$ and $\chi$ dilute slower than radiation. Therefore, after $\sigma$
comes to dominate the energy density, we can approximate

$\begin{split}\zeta\approx(1-f_{\chi})\zeta_{\sigma}+f_{\chi}\zeta_{\chi}=\zeta_{\sigma}+{1\over
3}f_{\chi}S_{\chi},\end{split}$ (3.5)

where $f_{\chi}\approx 1-f_{\sigma}$ and
$S_{\chi}=3(\zeta_{\chi}-\zeta_{\sigma})$ denotes the isocurvature fluctuation
of $\chi$ relative to $\sigma$. $S_{\chi}$ is given by
$(\delta\rho_{\chi}/\rho_{\chi}-\delta\rho_{\sigma}/\rho_{\sigma})\approx\delta\rho_{\chi}/\rho_{\chi}$
where we have used the isocurvature initial condition
$\delta\rho_{\chi}+\delta\rho_{\sigma}=0$ along with the fact that
$\rho_{\sigma}\gg\rho_{\chi}$ much before matter-radiation equality. We can
write $S_{\chi}\approx 2\delta\chi/\chi_{0}$, where $\chi_{0}$ is the
homogeneous field value of the $\chi$ field. The power spectrum of $S_{\chi}$
is then given by

$\begin{split}{\cal
P}_{S_{\chi}}=\left({H\over\pi\chi_{0}}\right)^{2}.\end{split}$ (3.6)

If $\chi$ is the QCD axion then $\chi_{0}=\theta_{i}F_{a}$ where $\theta_{i}$
is the initial misalignment angle and $F_{a}$ is the axion decay constant.

###  3.2 Interactions

To describe the interaction between $\chi$ and $\sigma$, we impose a (softly-
broken) shift symmetry on the system. In the UV, such a shift symmetry can
originate from gauge invariance, as we argue later. We also demand a
consistent EFT expansion such that higher order terms that we do not consider
contribute subdominantly, and the derivative expansion is under control. The
leading interactions are then given by

$\begin{split}{\cal L}=-{1\over 2}(\partial\sigma)^{2}-{1\over
2}(\partial\chi)^{2}-{1\over 2}(\partial
s)^{2}+\varepsilon\partial\sigma\partial\chi-{1\over
2}m^{2}_{s}s^{2}-(\partial\chi)^{2}{s\over v_{s}}+\cdots.\end{split}$ (3.7)

Here and below, we use the shorthand notation
$(\partial\sigma)^{2}\equiv\partial_{\mu}\sigma\partial^{\mu}\sigma$,
$\partial\sigma\partial\chi\equiv\partial_{\mu}\sigma\partial^{\mu}\chi$, and
similarly for other quantities. We have included a mass term for $\sigma$, as
needed to explain the tilt of the scalar power spectrum, while $\chi$ is
massless. We also have a kinetic mixing $\varepsilon$ between $\chi$ and
$\sigma$, and have kept the radial mode $s$ associated with $\chi$, defined as
$\Sigma=(s+v_{s})\exp(i\chi/v_{s})/\sqrt{2}$. One can remove the kinetic
mixing by performing a field redefinition,

$\begin{split}\chi\rightarrow\chi+\varepsilon\sigma.\end{split}$ (3.8)

This field shift gives rise to higher dimensional derivative mixing terms,

$\begin{split}(\partial\chi)^{2}{s\over
v_{s}}\rightarrow(\partial\chi)^{2}{s\over
v_{s}}+2\varepsilon(\partial\chi)(\partial\sigma){s\over
v_{s}}+\varepsilon^{2}(\partial\sigma)^{2}{s\over v_{s}}.\end{split}$ (3.9)

The last term above gives a subdominant correction to the curvaton kinetic
term when $s$ is set to some homogeneous value $\lesssim v_{s}$. Assuming the
radial mode mass $m_{s}>H$, we can integrate it out,222See [37] for a
complementary study where the effects of on-shell dynamics of $s$ were
studied. Such effects give oscillatory adiabatic-isocurvature NG as in the
context of cosmological collider physics [38, 39]. and that gives rise to
various dimension-8 local interactions such as,

$\begin{split}{1\over
2m_{s}^{2}v_{s}^{2}}\left((\partial\chi)^{4}+4\varepsilon(\partial\chi)^{2}(\partial\chi\partial\sigma)+2\varepsilon^{2}(\partial\chi)^{2}(\partial\sigma)^{2}\right.\\\
\left.+4\varepsilon^{2}(\partial\chi\partial\sigma)^{2}+4\varepsilon^{3}(\partial\sigma)^{2}(\partial\chi\partial\sigma)+\varepsilon^{4}(\partial\sigma)^{4}\right).\end{split}$
(3.10)

The homogeneous part of the curvaton $\sigma_{0}$, which follows the EOM
$\dot{\sigma}_{0}\approx-m_{\sigma}^{2}\sigma_{0}/(3H)$, gives correction to
the leading order Lagrangian considered in (3.7).333Owing to the negligible
mass of $\chi$, similar corrections coming from $\dot{\chi}_{0}$ can be
neglected. For consistency, we demand such corrections to be parametrically
small.

  * •

Correction to curvaton kinetic term:

$\begin{split}{\varepsilon^{4}\dot{\sigma}_{0}^{2}\over
m_{s}^{2}v_{s}^{2}}<1.\end{split}$ (3.11)

  * •

Contribution to the speed of curvaton fluctuation:

$\begin{split}c_{\sigma}^{2}=\left(1+{4\dot{\sigma}_{0}^{2}\varepsilon^{4}\over
m_{s}^{2}v_{s}^{2}}\right)^{-1}\approx 1.\end{split}$ (3.12)

  * •

Correction to $\chi$ kinetic term:

$\begin{split}{\varepsilon^{2}\dot{\sigma}_{0}^{2}\over
m_{s}^{2}v_{s}^{2}}<1.\end{split}$ (3.13)

  * •

Contribution to the speed of $\chi$ fluctuation:

$\begin{split}c_{\chi}^{2}=\left(1+{4\dot{\sigma}_{0}^{2}\varepsilon^{2}\over
m_{s}^{2}v_{s}^{2}}\right)^{-1}\approx 1.\end{split}$ (3.14)

  * •

Correction to curvaton-$\chi$ mixing:

$\begin{split}{2\varepsilon^{3}\dot{\sigma}_{0}^{2}\over
m_{s}^{2}v_{s}^{2}}<\varepsilon.\end{split}$ (3.15)

Therefore, the most stringent among these requirements is given by (3.13),
(3.14), and (3.15) with $\varepsilon<1$. We can summarize these restriction by
defining
$\kappa\equiv\varepsilon^{2}\dot{\sigma}_{0}^{2}/(m_{s}^{2}v_{s}^{2})$ and
imposing $\kappa<1$.

#### Origin of kinetic mixing.

Eq. (3.7) describes kinetic mixing $\varepsilon$ between two light fields
$\sigma$ and $\chi$. From a UV perspective, $\varepsilon$ can arise in several
ways, including a scenario involving one flat extra spatial dimension.
Consider two $U(1)$ gauge fields, $A_{M}$, and $B_{M}$, with gauge couplings
$g_{5,A}$ and $g_{5,B}$, and field strengths $F_{MN}$ and $G_{MN}$,
respectively. The 5D Lagrangian is given by

$\begin{split}{\cal L}_{5\rm D}=-{1\over 4g_{5,A}^{2}}\int{\rm
d}^{5}x\sqrt{-g_{5D}}F^{MN}F_{MN}-{1\over 4g_{5,B}^{2}}\int{\rm
d}^{5}x\sqrt{-g_{5D}}G^{MN}G_{MN}\\\ -{\varepsilon_{5}\over
g_{5,A}g_{5,B}}\int{\rm d}^{5}x\sqrt{-g_{5D}}F^{MN}G_{MN}+\cdots.\end{split}$
(3.16)

The gauge kinetic mixing $\varepsilon_{5}$ can arise by integrating out
massive 5D fermions that are charged under both the $U(1)$’s. To reduce to a
4D theory, we can choose orbifold boundary conditions such that only $A_{5}$
and $B_{5}$ have zero modes, while $A_{\mu}$ and $B_{\mu}$ do not. Then the
low-energy 4D Lagrangian is given by

$\begin{split}{\cal L}_{4\rm D,~{}eff}=-{2\pi R\over 2g_{5,A}^{2}}\int{\rm
d}^{4}x\sqrt{-g_{4D}}(\partial_{\mu}A_{5})^{2}-{2\pi R\over
2g_{5,B}^{2}}\int{\rm d}^{4}x\sqrt{-g_{4D}}(\partial_{\mu}B_{5})^{2}\\\ -{4\pi
R\varepsilon_{5}\over g_{5,A}g_{5,B}}\int{\rm
d}^{4}x\sqrt{-g_{4D}}(\partial^{\mu}A_{5})(\partial_{\mu}B_{5})+\cdots.\end{split}$
(3.17)

Canonically normalizing the $A_{5}$ and $B_{5}$ kinetic terms, we can rewrite
the above as the first three terms of (3.7) by identifying $A_{5}$ and $B_{5}$
with $\sigma$ and $\chi$, respectively, along with
$\varepsilon=2\varepsilon_{5}$. The shift symmetries in the 4D effective field
theory come from the 5D gauge transformations of $A_{5}$ and $B_{5}$.

###  3.3 Magnitude of Bispectrum

Eq. (3.10) determines the cubic interaction vertices between $\sigma$ and
$\chi$. This can be obtained by expressing the fields in terms of a
homogeneous background and a fluctuating piece:
$\sigma(t,\mathbf{x})=\sigma_{0}(t)+\delta\sigma(t,\mathbf{x})$ and
$\chi(t,\mathbf{x})=\chi_{0}(t)+\delta\chi(t,\mathbf{x})$. Practically,
$\dot{\chi}_{0}\approx 0$ since $\chi$ is approximately massless during
inflation, and therefore we will use
$\partial_{\mu}\delta\chi(t,\mathbf{x})\approx\partial_{\mu}\chi(t,\mathbf{x})$.
With that the cubic interactions are given by,

$\begin{split}{\cal L}=-{\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}}\left(2\varepsilon(\partial\chi)^{2}\dot{\chi}+2\varepsilon^{2}\dot{\delta\sigma}(\partial\chi)^{2}+4\varepsilon^{2}(\partial\chi\partial\delta\sigma)\dot{\chi}\right.\\\
\left.+2\varepsilon^{3}(\partial\delta\sigma)^{2}\dot{\chi}+4\varepsilon^{3}\dot{\delta\sigma}(\partial\chi\partial\delta\sigma)+2\varepsilon^{4}\dot{\delta\sigma}(\partial\delta\sigma)^{2}\right)+\cdots.\end{split}$
(3.18)

These determine the various types of mixed adiabatic-isocurvature bispectra.
The associated shapes have been computed in [40] which we can express in terms
of the parameters from (3.18). To that end, we define a kinematic function

$\begin{split}Q(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})={1\over
k_{1}^{3}k_{2}^{3}k_{3}^{3}}\left(\frac{2k_{1}^{2}k_{2}^{2}k_{3}^{2}}{k_{123}^{3}}-\frac{(\mathbf{k}_{2}\cdot\mathbf{k}_{3})k_{1}^{2}}{k_{123}}\left(1+\frac{k_{23}}{k_{123}}+\frac{2k_{2}k_{3}}{k_{123}^{2}}\right)\right),\end{split}$
(3.19)

where $k_{i}=|\mathbf{k}_{i}|$ and $k_{ij\cdots
m}=|\mathbf{k}_{i}|+|\mathbf{k}_{j}|+\cdots+|\mathbf{k}_{m}|$. We also define
$Q_{\rm
full}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})=Q(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})+Q(\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{1})+Q(\mathbf{k}_{3},\mathbf{k}_{1},\mathbf{k}_{2})$.
The various three-point functions can then be expressed as

$\begin{split}\langle\sigma(\mathbf{k}_{1})\sigma(\mathbf{k}_{2})\sigma(\mathbf{k}_{3})\rangle^{\prime}&={\varepsilon^{4}\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}}Q_{\rm
full}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}),\\\
\langle\chi(\mathbf{k}_{1})\sigma(\mathbf{k}_{2})\sigma(\mathbf{k}_{3})\rangle^{\prime}&={\varepsilon^{3}\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}}Q_{\rm
full}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}),\\\
\langle\chi(\mathbf{k}_{1})\chi(\mathbf{k}_{2})\sigma(\mathbf{k}_{3})\rangle^{\prime}&={\varepsilon^{2}\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}}Q_{\rm
full}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}),\\\
\langle\chi(\mathbf{k}_{1})\chi(\mathbf{k}_{2})\chi(\mathbf{k}_{3})\rangle^{\prime}&={\varepsilon\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}}Q_{\rm
full}(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3}),\end{split}$ (3.20)

where $\sigma(\mathbf{k})$ and $\chi(\mathbf{k})$ are Fourier transforms of
$\delta\sigma(t,\mathbf{x})$ and $\delta\chi(t,\mathbf{x})$, respectively. To
convert the above into correlation function between $\zeta$ and $S$, we can
use $\zeta=(2/3)(\delta\sigma/\sigma_{0})$ and $S=f_{\rm DM}S_{\chi}=2f_{\rm
DM}(\delta\chi/\chi_{0})$, where $f_{\rm DM}$ is the fraction of DM abundance
determined by $\chi$.444The DM isocurvature is defined as $S=3(\zeta_{\rm
DM}-\zeta_{\gamma})$. Writing $\zeta_{\rm DM}=f_{\rm DM}\zeta_{\chi}+(1-f_{\rm
DM})\zeta_{c}$, where $\zeta_{c}$ is the fluctuation in the other (adiabatic)
component of DM, and setting $\zeta_{c}=\zeta_{\gamma}$, we get $S=f_{\rm
DM}S_{\chi}$. While the shapes in (3.20) are different from the ones used in
the Planck search (2.4), we can make a rough comparison by defining the
following quantities [40]:

$\begin{split}\tilde{f}^{\zeta\zeta\zeta}_{\rm
NL,equi}\equiv&{\langle\zeta(\mathbf{k})\zeta(\mathbf{k})\zeta(\mathbf{k})\rangle^{\prime}\over
P_{\zeta}(k)^{2}}={7\over 6}{81\sigma_{0}^{4}\over 4H^{4}}\left({2H\over
3\sigma_{0}}\right)^{3}{\varepsilon^{4}H^{2}\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}},\\\ \tilde{f}^{\zeta\zeta S}_{\rm
NL,equi}\equiv&{\langle\zeta(\mathbf{k})\zeta(\mathbf{k})S(\mathbf{k})\rangle^{\prime}\over
P_{\zeta}(k)^{2}}={7\over 6}{81\sigma_{0}^{4}\over 4H^{4}}\left({2H\over
3\sigma_{0}}\right)^{2}\left(2Hf_{\rm
DM}\over\chi_{0}\right){\varepsilon^{3}H^{2}\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}},\\\ \tilde{f}^{\zeta SS}_{\rm
NL,equi}\equiv&{\langle\zeta(\mathbf{k})S(\mathbf{k})S(\mathbf{k})\rangle^{\prime}\over
P_{\zeta}(k)^{2}}={7\over 6}{81\sigma_{0}^{4}\over 4H^{4}}\left({2H\over
3\sigma_{0}}\right)\left(2Hf_{\rm
DM}\over\chi_{0}\right)^{2}{\varepsilon^{2}H^{2}\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}},\\\ \tilde{f}^{SSS}_{\rm NL,equi}\equiv&{\langle
S(\mathbf{k})S(\mathbf{k})S(\mathbf{k})\rangle^{\prime}\over
P_{\zeta}(k)^{2}}={7\over 6}{81\sigma_{0}^{4}\over 4H^{4}}\left(2Hf_{\rm
DM}\over\chi_{0}\right)^{3}{\varepsilon H^{2}\dot{\sigma}_{0}\over
m_{s}^{2}v_{s}^{2}}.\end{split}$ (3.21)

To compare with observational bounds on isocurvature power spectrum, we define

$\begin{split}\beta_{\rm iso}\equiv{P_{S}\over
P_{\zeta}+P_{S}}\Rightarrow{\beta_{\rm iso}\over 1-\beta_{\rm
iso}}={P_{S}\over P_{\zeta}}={9\sigma_{0}^{2}f_{\rm
DM}^{2}\over\chi_{0}^{2}}.\end{split}$ (3.22)

In terms of $\beta_{\rm iso}$, (3.21) reads:

$\begin{split}\tilde{f}^{\zeta\zeta\zeta}_{\rm
NL,equi}&={7H\dot{\sigma}_{0}\varepsilon^{4}\sigma_{0}\over
m_{s}^{2}v_{s}^{2}}=-{21H^{2}\kappa\varepsilon^{2}\over m_{\sigma}^{2}},\\\
\tilde{f}^{\zeta\zeta S}_{\rm
NL,equi}&={7H\dot{\sigma}_{0}\varepsilon^{3}\sigma_{0}\over
m_{s}^{2}v_{s}^{2}}\left({\beta_{\rm iso}\over 1-\beta_{\rm
iso}}\right)^{1/2}=-{21H^{2}\kappa\varepsilon\over
m_{\sigma}^{2}}\left({\beta_{\rm iso}\over 1-\beta_{\rm iso}}\right)^{1/2},\\\
\tilde{f}^{\zeta SS}_{\rm
NL,equi}&={7H\dot{\sigma}_{0}\varepsilon^{2}\sigma_{0}\over
m_{s}^{2}v_{s}^{2}}\left({\beta_{\rm iso}\over 1-\beta_{\rm
iso}}\right)=-{21H^{2}\kappa\over m_{\sigma}^{2}}\left({\beta_{\rm iso}\over
1-\beta_{\rm iso}}\right),\\\ \tilde{f}^{SSS}_{\rm
NL,equi}&={7H\dot{\sigma}_{0}\varepsilon\sigma_{0}\over
m_{s}^{2}v_{s}^{2}}\left({\beta_{\rm iso}\over 1-\beta_{\rm
iso}}\right)^{3/2}=-{21H^{2}\kappa\over\varepsilon
m_{\sigma}^{2}}\left({\beta_{\rm iso}\over 1-\beta_{\rm
iso}}\right)^{3/2},\end{split}$ (3.23)

where $\kappa=\varepsilon^{2}\dot{\sigma}_{0}^{2}/(m_{s}^{2}v_{s}^{2})<1$ to
obey the constraints (3.13)–(3.15). We have used the homogeneous EOM for
$\sigma$, $3H\dot{\sigma}_{0}\approx-m_{\sigma}^{2}\sigma_{0}$. To obtain an
estimate of the magnitude of NG, we choose the following benchmark parameter
values: $m_{\sigma}^{2}=0.03H^{2}$, $\varepsilon=0.3$, $v_{s}=10H$, and
$m_{s}=2H$. We also use $H/\sigma_{0}=4.3\times 10^{-4}$, determined by the
scalar power spectrum, and $\beta_{\rm iso}=0.038$, current upper limit on the
magnitude of the DM isocurvature power spectrum [41]. These parameter choices
imply $\kappa\approx 0.1$ and

$\begin{split}\tilde{f}^{\zeta\zeta\zeta}_{\rm
NL,equi}\approx-7.7,\tilde{f}^{\zeta\zeta S}_{\rm
NL,equi}\approx-5.1,\tilde{f}^{\zeta SS}_{\rm
NL,equi}\approx-3.4,\tilde{f}^{SSS}_{\rm NL,equi}\approx-2.2.\end{split}$
(3.24)

Although the bounds in (2.6) were derived using a different shape, we can make
a rough comparison using (3.24). For example, the Planck constraint (2.6)
gives $\tilde{f}_{\rm NL}^{\zeta,\zeta S}=-5.6\pm 8.4$, comparable to the
values obtained in (3.24). We note that the NG involving two adiabatic and one
isocurvature perturbation, $\tilde{f}^{\zeta\zeta S}_{\rm NL,equi}$, is the
leading probe of DM isocurvature NG. While $\tilde{f}^{\zeta SS}_{\rm
NL,equi}$ and $\tilde{f}^{SSS}_{\rm NL,equi}$ are not too different in
strengths from $\tilde{f}^{\zeta\zeta S}_{\rm NL,equi}$, the observational
bounds on NG involving more than one isocurvature perturbation is an order of
magnitude weaker (2.6). Therefore, in the following, we focus on the
correlation function $\langle\zeta\zeta S\rangle$ to analyze the difference in
shapes between (3.20) and the ones considered in the Planck search (2.4).
Similar steps can be repeated for models where $\langle\zeta SS\rangle$ and
$\langle SSS\rangle$ are relevant. In Fig. 1 we show the absolute size of
$\tilde{f}_{\rm NL,equi}^{\zeta\zeta S}$ as a function of $v_{s}/H$ and
$\varepsilon$.

Figure 1: Absolute size of $\tilde{f}_{\rm NL,equi}^{\zeta\zeta S}$. We fix
$m_{s}=2H$, $m_{\sigma}^{2}=0.03H^{2}$, and saturate the current isocurvature
power spectrum bound $\beta_{\rm iso}=0.038$. In the region above the orange
line $\kappa>1$ and the corrections discussed in (3.13)–(3.15) become
important. We can roughly recast the current constraint on this parameter
space and that rules out the region above the upper white line, labeled
‘$f>125$’, where $\tilde{f}_{\rm NL,equi}^{\zeta\zeta S}>125$. A dedicated
template would improve this bound by roughly a factor of 3, and thereby be
sensitive to the region above the lower white line labeled ‘$f>42$’ (see Sec.
3.4 for more details).

###  3.4 Shape of Bispectrum

In Fig. 2, we compare the two classes of shapes given by (2.4) and (3.20). For
convenience, we rewrite the shapes by ignoring the overall prefactors which
are not relevant for examining the shapes,

$\begin{split}\langle\zeta(\mathbf{k}_{1})\zeta(\mathbf{k}_{2})S(\mathbf{k}_{3})\rangle^{\prime}_{\rm
local}&\propto\tilde{f}_{\rm NL}^{\zeta,\zeta S}\left({1\over
k_{2}^{3}}{1\over k_{3}^{3}}+{1\over k_{1}^{3}}{1\over
k_{3}^{3}}\right)+\tilde{f}_{\rm NL}^{S,\zeta\zeta}{1\over k_{1}^{3}}{1\over
k_{2}^{3}},\\\
\langle\zeta(\mathbf{k}_{1})\zeta(\mathbf{k}_{2})S(\mathbf{k}_{3})\rangle^{\prime}_{\rm
derivative}&\propto
Q(\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3})+Q(\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{1})+Q(\mathbf{k}_{3},\mathbf{k}_{1},\mathbf{k}_{2}).\end{split}$
(3.25)

The local shape is symmetric under
$\mathbf{k_{1}}\leftrightarrow\mathbf{k_{2}}$ and is only dependent on the
ratio $\tilde{f}_{\rm NL}^{\zeta,\zeta S}/\tilde{f}_{\rm NL}^{S,\zeta\zeta}$.
On the other hand, the derivative shape is symmetric under the cyclic
permutation of $\mathbf{k_{1}},\mathbf{k_{2}},\mathbf{k_{3}}$. To compare the
shapes we choose $k_{1}$ to be the largest side of the triangle and define the
ratios $x_{2}=k_{2}/k_{1}$ and $x_{3}=k_{3}/k_{1}$. Since at the moment we are
interested in the shape, not the overall magnitude, we normalize both the
shapes to be unity at $x_{2}=x_{3}=1$. We also restrict the shapes in the
region $x_{2}+x_{3}\geq 1$ (triangle inequality). From Fig. 2 we note that the
local shape is not symmetric under the exchange $x_{2}\leftrightarrow x_{3}$,
while the derivative shape is, as expected from (3.25). For the chosen
parameter points, the two shapes are significantly different.

To quantify the difference between the two shapes, we can compute the cosine
between them as per (2.8) for a range of values of $\tilde{f}_{\rm
NL}^{\zeta,\zeta S}$ and $\tilde{f}_{\rm NL}^{S,\zeta\zeta}$. The result is
shown in Fig. 3 for the absolute value of the cosine. The cosine depends only
on the ratio $\tilde{f}_{\rm NL}^{S,\zeta\zeta}/\tilde{f}_{\rm
NL}^{\zeta,\zeta S}$, as evident from (3.25). In particular, the result is
symmetric for $(\tilde{f}_{\rm NL}^{\zeta,\zeta S},\tilde{f}_{\rm
NL}^{S,\zeta\zeta})\rightarrow-(\tilde{f}_{\rm NL}^{\zeta,\zeta
S},\tilde{f}_{\rm NL}^{S,\zeta\zeta})$. We utilize this feature to show the
cosine as a function of this ratio in the right panel of Fig. 3. The minimum
value of the cosine is reached for $\tilde{f}_{\rm
NL}^{S,\zeta\zeta}=-2\tilde{f}_{\rm NL}^{\zeta,\zeta S}$. For this value the
local shape (3.25) vanishes in the equilateral limit. We see for the entire
range of the parameter space, the cosine is smaller than $0.5$ in magnitude.
In fact, for regions of parameter space, the cosine can be $\sim 0.1$. This
implies, even without including new data, performing a reanalysis with the
dedicated template can improve the sensitivity by an order of magnitude
(inverse of the cosine).

Figure 2: Left. Local shape (3.25) for $\tilde{f}_{\rm NL}^{\zeta,\zeta
S}=-5.6$ and $\tilde{f}_{\rm NL}^{S,\zeta\zeta}=346$, as per the central
values in (2.6). Right. Derivative shape (3.25). We have imposed the triangle
inequality $x_{2}+x_{3}\geq 1$. Here $x_{2}=k_{2}/k_{1}$ and
$x_{3}=k_{3}/k_{1}$ and both the bispectra are normalized to unity at
$x_{2}=x_{3}=1$. The two shapes are significantly different from each other,
as further quantified via the cosine in Fig. 3.

Figure 3: Left. Cosine between the local shape and the derivative shape (3.25)
for a range of $\tilde{f}_{\rm NL}^{\zeta,\zeta S}$ and $\tilde{f}_{\rm
NL}^{S,\zeta\zeta}$ motivated by (2.6). Changing the signs of both
$\tilde{f}_{\rm NL}^{\zeta,\zeta S}$ and $\tilde{f}_{\rm NL}^{S,\zeta\zeta}$
gives the same absolute value of the cosine. Rescaling both $\tilde{f}_{\rm
NL}^{\zeta,\zeta S}$ and $\tilde{f}_{\rm NL}^{S,\zeta\zeta}$ by the same
factor also gives the same value of the cosine, also seen in the above plot.
Right. These features are utilized in the right panel where we show the cosine
as a function of the ratio $\tilde{f}_{\rm NL}^{S,\zeta\zeta}/\tilde{f}_{\rm
NL}^{\zeta,\zeta S}$. We see that by using the correct shape, we can improve
the sensitivity by at least a factor of two in comparison with a search
treating $\tilde{f}_{\rm NL}^{\zeta,\zeta S}$ and $\tilde{f}_{\rm
NL}^{S,\zeta\zeta}$ as free parameters. At the same time, there could be one
order of magnitude difference in sensitivity for certain choices of
$\tilde{f}_{\rm NL}^{\zeta,\zeta S}$ and $\tilde{f}_{\rm NL}^{S,\zeta\zeta}$.

We can also do a very rough recast of the existing bound by calculating the
(3D) ‘fudge factor’ as in [26] which accounts for both the cosine and the norm
of the signal.555Note, a better recast can be obtained by computing the 2D
fudge factors [26], as well as, including the correlation between
$\tilde{f}_{\rm NL}^{\zeta,\zeta S}$ and $\tilde{f}_{\rm NL}^{S,\zeta\zeta}$.
For $\tilde{f}_{\rm NL}^{\zeta,\zeta S}=-5.6$ and $\tilde{f}_{\rm
NL}^{S,\zeta\zeta}=346$, the fudge factor is $\approx 0.067$ and the cosine is
$\approx 0.36$. This means the current constraint on $\tilde{f}_{\rm
NL,equi}^{\zeta\zeta S}$, based on (2.6), is roughly $\sim 8.4/0.067\approx
125$. Including the derivative shape template would tighten this bound by a
factor of the inverse cosine, $\sim 3$. We illustrate this improvement in Fig.
1.

###  3.5 Interpretation for the QCD Axion

Figure 4: QCD Axion parameter space as a function of the decay constant
$f_{a}$ and the inflationary Hubble scale $H$. We choose $\varepsilon=0.3$,
$m_{s}=H$, $m_{\sigma}^{2}=0.03H^{2}$, and $\theta_{i}=1$. In the region above
the dashed line $|\tilde{f}_{\rm NL}^{S,SS}|>100$ (cf. (3.31)) and is ruled
out by the current isocurvature NG searches. We also show the cooling bound on
the QCD axion [42] and bound from isocurvature power spectrum (Iso PS) [41].
In the region labeled $\kappa>1$, the restrictions (3.13), (3.14), and (3.15)
are not obeyed. We show the projected sensitivity of a search for the
templates proposed in this work for two benchmark choices $\tilde{f}_{\rm
NL,equi}^{\zeta\zeta S}=1,10$ which would be able to cover the regions above
the teal lines.

While the above discussion has been general, we now interpret the constraints
in terms of the QCD axion, which plays the role of $\chi$ in the model
discussed above. The cosmological abundance when the Peccei-Quinn symmetry is
broken before inflation is given by (ignoring anharmonic effects) [43],

$\begin{split}\Omega_{a}h^{2}\approx 0.12\left({F_{a}\over 9\times
10^{11}~{}{\rm GeV}}\right)^{1.165}\theta_{i}^{2}.\end{split}$ (3.26)

The isocurvature fraction is given by,

$\begin{split}{\beta_{\rm iso}\over 1-\beta_{\rm iso}}={P_{S}\over
P_{\zeta}}=\left(\Omega_{a}\over\Omega_{\rm
DM}\right)^{2}{H^{2}\over\pi^{2}F_{a}^{2}\theta_{i}^{2}{\cal
P}_{\zeta}}.\end{split}$ (3.27)

Using Eq. (3.27) to trade initial misalignment $\theta_{i}$ with other
parameters, and saturating the current constraint $\beta_{\rm iso}<0.038$
[41], we can derive

$\begin{split}H<2.6\times 10^{7}~{}{\rm
GeV}\left({0.12\over\Omega_{a}h^{2}}\right)^{1/2}\left({F_{a}\over 9\times
10^{11}~{}{\rm GeV}}\right)^{0.42}.\end{split}$ (3.28)

Requiring $|\tilde{f}^{\zeta\zeta S}_{\rm NL,equi}|>10$ implies,

$\begin{split}H^{3}>10\pi F_{a}\theta_{i}{\cal
P}_{\zeta}^{1/2}\left(\Omega_{\rm
DM}\over\Omega_{a}\right)\left(m_{\sigma}^{2}\over
21\kappa\varepsilon\right).\end{split}$ (3.29)

The parameter $\kappa$ is related to the other parameters as,
$\kappa=\varepsilon^{2}(m_{\sigma}^{2}\sigma_{0})^{2}/(9H^{2}m_{s}^{2}F_{a}^{2})$.
There is also a model-independent contribution to isocurvature NG. This arises
because the isocurvature fluctuation contains a term quadratic in the axion
fluctuation,

$\begin{split}S={\Omega_{a}\over\Omega_{\rm DM}}\left(2{\left(\delta a\over
a\right)}+\left(\delta a\over a\right)^{2}\right).\end{split}$ (3.30)

This implies an isocurvature NG [15] that is independent of the kinetic mixing
parameter $\varepsilon$,

$\begin{split}\tilde{f}_{\rm NL}^{S,SS}={1\over 2}{\Omega_{\rm
DM}\over\Omega_{a}}{P_{S}^{2}\over P_{\zeta}^{2}}.\end{split}$ (3.31)

Using (3.27) we can project the bound on $\tilde{f}_{\rm NL}^{S,SS}$ in terms
of $H$ and $f_{a}$. The parameter space ruled out by the current Planck
constraint on $\tilde{f}_{\rm NL}^{S,SS}$ is shown via the dashed line in Fig.
4. We summarize the other constraints in Fig. 4. The regions above the teal
solid and dot-dashed lines would be covered by a future search of
$\langle\zeta\zeta S\rangle$ NG with the shapes considered in this work. Such
a search would improve upon the $\langle SSS\rangle$ NG constraints by an
order of magnitude in parts of the parameter space, and go beyond the standard
isocurvature power spectrum searches for the QCD axion.

##  4 Generalities

In the above, we have given an example in the context of the curvaton scenario
in which observable isocurvature NG can arise. It is natural to ask how
general our result is, and clarify our non-minimal choice involving curvaton
dynamics. To that end, we explore two possibilities: first, a model-agnostic
EFT description based on [28]; second, a slow-roll inflationary scenario
without the curvaton dynamics. As we show below, we find that while the EFT
description can allow a stronger interaction between the adiabatic and
isocurvature fluctuation with an observable NG signal, its UV completion in a
slow-roll inflaton model results in a significantly suppressed signal. It
remains to be seen whether an observable signal can arise in a scenario where
the scalar perturbations emerge from the inflaton.

###  4.1 Goldstone EFT

In the Goldstone EFT of inflation [28] the inflaton is modeled as a Goldstone
boson arising from spontaneous breaking of time-translation symmetry. We can
then write an EFT for the fluctuation of the adiabatic mode, along with any
other spectator field such as $\chi$, considered above. Such an EFT including
various operators up to quartic terms has been described in [40], building up
on [44]. Here we just select a particular cubic term to illustrate that
isocurvature NG can be large and observable while being in theoretical
control. In particular, we focus on a term

$\begin{split}-{1\over\Lambda_{\rm
EFT}^{2}}(\partial\pi_{c})^{2}\dot{\chi},\end{split}$ (4.1)

allowed by the symmetries of the theory. The field $\pi_{c}$ is the
canonically normalized inflaton fluctuation field and the EFT cutoff scale is
denoted by $\Lambda_{\rm EFT}$, which can be below the symmetry breaking scale
$f_{\pi}=(2M_{\rm pl}^{2}|\dot{H}|c_{s})^{1/4}\approx 60H$. However, as long
as $\Lambda_{\rm EFT}\gg H$, we can have a controlled theory of fluctuations
between $H$ and $\Lambda_{\rm EFT}$.

The operator in (4.1) has the same structure as the one in (3.18) with the
replacements, $\delta\sigma\rightarrow\pi_{c}$,
$2/(3\sigma_{0})\rightarrow-H/f_{\pi}^{2}$ and
$2\varepsilon^{3}\dot{\sigma}_{0}/(m_{s}^{2}v_{s}^{2})\rightarrow
1/\Lambda_{\rm EFT}^{2}$. Thus the associated strength of the NG is given by,

$\begin{split}\tilde{f}_{\rm NL,equi}^{\zeta\zeta S}={14\over
3}\left(f_{\pi}\over H\right)^{4}\left({2f_{\rm
DM}H\over\chi_{0}}\right)\left(H^{2}\over 2\Lambda_{\rm
EFT}^{2}\right).\end{split}$ (4.2)

For example choosing $\Lambda_{\rm EFT}=5H$ and saturating the isocurvature
power spectrum constraint $\beta_{\rm iso}<0.038$, gives $\tilde{f}_{\rm
NL,equi}^{\zeta\zeta S}\approx 67$ from the above. The benchmark value
considered in (3.24) corresponds to $\Lambda_{\rm EFT}\approx 18H$ in the
context of this Goldstone EFT and correspondingly, the value in (3.24) is an
order of magnitude smaller.

###  4.2 Slow-roll EFT without the Curvaton Field

We now consider similar NG but within the slow-roll EFT where we assume a
slowly rolling inflaton field $\phi$ provides both the homogeneous expansion
and the adiabatic density fluctuations. In particular, we assume there is no
curvaton field in this scenario, while there is still a spectator field $\chi$
that sources DM isocurvature. The leading terms governing the interactions
between $\phi$ and $\chi$ are given by a set of terms similar to (3.10),

$\begin{split}{c_{1}\over\Lambda^{4}}(\partial\phi)^{4}+{c_{2}\over\Lambda^{4}}(\partial\phi)^{2}(\partial\phi\partial\chi)+{c_{3}\over\Lambda^{4}}(\partial\phi)^{2}(\partial\chi)^{2}+{c_{4}\over\Lambda^{4}}(\partial\phi\partial\chi)^{2}+{c_{5}\over\Lambda^{4}}(\partial\phi\partial\chi)(\partial\chi)^{2}+{c_{6}\over\Lambda^{4}}(\partial\chi)^{4}.\end{split}$
(4.3)

Requiring a controlled derivative expansion then forces

$\begin{split}\Lambda>\sqrt{\dot{\phi}_{0}},\end{split}$ (4.4)

from the term with the coefficient $c_{1}\sim 1$. However, adopting a simple
power counting scheme where all $c_{i}\sim 1$, implies a small value for
isocurvature NG. To see this we can focus on the term with the coefficient
$c_{2}$ which has the same structure as one term in (3.10) with the
replacements $2\varepsilon^{3}/(m_{s}^{2}v_{s}^{2})\rightarrow
c_{2}/\Lambda^{4}$, and $\sigma\rightarrow\phi$. This implies,

$\begin{split}\tilde{f}_{\rm NL,equi}^{\zeta\zeta S}={14\over
3}{\dot{\phi}_{0}^{2}\over H^{4}}\left({2f_{\rm
DM}H\over\chi_{0}}\right){c_{2}H^{2}\dot{\phi}_{0}\over
2\Lambda^{4}}.\end{split}$ (4.5)

Choosing $\Lambda=2\sqrt{\dot{\phi}_{0}}$ to respect the above derivative
expansion bound, and saturating the isocurvature power spectrum constraint
$\beta_{\rm iso}<0.038$, gives $\tilde{f}_{\rm NL,equi}^{\zeta\zeta S}\approx
0.03$ from the above. This is significantly smaller than the values considered
previously. This suppression originates from requiring (1) the same field that
drives the inflationary expansion also generates the adiabatic fluctuations;
and (2) the same EFT cutoff for the adiabatic source and the isocurvature
source. The examples in Sec. 4.1 and Sec. 3 evade these restrictions by being
agnostic about the origin of the homogeneous inflationary expansion and having
another adiabatic source field (curvaton), respectively.

###  4.3 Effects of Shift Symmetry Breaking

One may wonder whether the suppression in the previous subsection originates
from requiring shift symmetric couplings between the inflaton and $\chi$. To
explore this possibility by taking a ‘bottom-up’ approach, we consider a non-
shift symmetric (NS) operator

$\begin{split}{\phi\over\Lambda_{\rm
NS}}(\partial\phi\partial\chi).\end{split}$ (4.6)

To have a controlled description of this interaction, the scale $\Lambda_{\rm
NS}$ has to be at least bigger than the field space traversed by the inflaton
in one Hubble time, i.e., $\Lambda_{\rm NS}>\dot{\phi}_{0}/H$. The induced NG
is given by,

$\begin{split}\tilde{f}_{\rm NL,equi}^{\zeta\zeta
S}\sim{\dot{\phi}_{0}^{2}\over H^{4}}\left(f_{\rm
DM}H\over\chi_{0}\right)\left(H\over\Lambda_{\rm NS}\right).\end{split}$ (4.7)

Saturating $\Lambda_{\rm NS}\sim\dot{\phi}_{0}/H$, this gives the same level
of suppression as (4.5). Hence even allowing for non-shift symmetric
couplings, the isocurvature NG is still suppressed.

##  5 Conclusion

In this work, we highlighted the importance of NG in the bispectrum involving
two adiabatic and one isocurvature perturbations. We argue that it can be
naturally sizable. Due to the good sensitivities to this NG observable in CMB
and LSS, it is a promising place to look for new dynamics in the early
universe. As an example, we present a simple model in which the adiabatic
perturbation is generated by a curvaton. There is also an additional light
scalar field responsible for generating the isocurvature perturbation, which
can also be part of dark matter. Both the curvaton and the new light field
have an approximate shift symmetry. We can consistently add a kinetic mixing
between them which preserves all the shift symmetries. We demonstrate such a
mixing can naturally arise in UV models. In the context of this model, we show
that the bispectrum
$\langle\zeta(\mathbf{k}_{1})\zeta(\mathbf{k}_{2})S(\mathbf{k}_{3})\rangle^{\prime}$
is the best observable to probe the presence of the light field and its
coupling to the curvaton. This model also predicts sizable NG in bispectra
involving more than one isocurvature perturbation. However, those are less
useful due to limited sensitivities in CMB and LSS observations. We interpret
our result in the scenario that the additional light scalar is the QCD axion,
and demonstrate that the bispectrum with one isocurvature perturbation is
indeed the most promising probe.

The accurate prediction of the shape of the bispectrum has an important impact
on its search. The bispectrum involving one or more isocurvature perturbations
has been searched for assuming the so-called local shape. In this work, we
emphasize that it can naturally be of a very different shape, and demonstrate
this fact by our example. We quantify the difference between the predicted
shape and the local shape and find that it could be significant. Hence, using
the correct shape can lead to improvement in sensitivity by at least a factor
of two, and in some cases almost an order of magnitude.

There will be a large influx of data from CMB and LSS observations in the
coming decades, see e.g. [45] for a summary. Developing effective and
comprehensive strategies will be crucial for extracting the maximal amount of
information. Given the multitude of possibilities, it is essential to explore
more scenarios and identify potential new signatures. Our work is a step in
this direction, and much more needs to be done. For example, as discussed in
Section 4, it is quite generic to expect a sizable NG signal of the kind
discussed in this paper from the point of view of the effective field theory
of inflation. At the same time, it remains to be seen whether this case can be
UV completed in a satisfying way. Regarding complete models, we have only
considered a special class, mostly motivated by enforcing shift symmetries.
There are good motivations to think this is an important ingredient. At the
same time, it would be interesting to explore more model building directions,
which may lead to new observables and new shapes of NG. We leave all of these
promising directions for future work.

## Acknowledgment

We thank J. Colin Hill and Ben Safdi for the helpful discussions. SK is
supported partially by the National Science Foundation (NSF) grant PHY-2210498
and the Simons Foundation. LTW is supported by the Department of Energy grant
DE-SC0013642. MG is supported by the Israel Science Foundation under Grant No.
1302/19 and 1424/23. MG is also supported by the US-Israeli BSF grant 2018236
and the NSF-BSF grant 2021779. LTW, MG, and SK are grateful to the Aspen
Center for Physics, supported by NSF grant PHY-2210452, for hospitality while
this work was in progress. We are also grateful to the Maryland Center for
Fundamental Physics, University of Maryland, for hosting the Cosmological
Probes of New Physics workshop where this project was initiated.

## References

  * [1] N. Bartolo, E. Komatsu, S. Matarrese, and A. Riotto, “Non-Gaussianity from inflation: Theory and observations,” Phys. Rept. 402 (2004) 103–266, arXiv:astro-ph/0406398.
  * [2] X. Chen, “Primordial Non-Gaussianities from Inflation Models,” Adv. Astron. 2010 (2010) 638979, arXiv:1002.1416 [astro-ph.CO].
  * [3] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. IX. Constraints on primordial non-Gaussianity,” Astron. Astrophys. 641 (2020) A9, arXiv:1905.05697 [astro-ph.CO].
  * [4] eBOSS Collaboration, E.-M. Mueller et al., “Primordial non-Gaussianity from the completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey II: measurements in Fourier space with optimal weights,” Mon. Not. Roy. Astron. Soc. 514 no. 3, (2022) 3396–3409, arXiv:2106.13725 [astro-ph.CO].
  * [5] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, “Constraints on Single-Field Inflation from the BOSS Galaxy Survey,” Phys. Rev. Lett. 129 no. 2, (2022) 021301, arXiv:2201.07238 [astro-ph.CO].
  * [6] G. D’Amico, M. Lewandowski, L. Senatore, and P. Zhang, “Limits on primordial non-Gaussianities from BOSS galaxy-clustering data,” arXiv:2201.11518 [astro-ph.CO].
  * [7] E. Komatsu and D. N. Spergel, “Acoustic signatures in the primary microwave background bispectrum,” Phys. Rev. D 63 (2001) 063002, arXiv:astro-ph/0005036.
  * [8] K. Enqvist and M. S. Sloth, “Adiabatic CMB perturbations in pre - big bang string cosmology,” Nucl. Phys. B 626 (2002) 395–409, arXiv:hep-ph/0109214.
  * [9] D. H. Lyth and D. Wands, “Generating the curvature perturbation without an inflaton,” Phys. Lett. B 524 (2002) 5–14, arXiv:hep-ph/0110002.
  * [10] T. Moroi and T. Takahashi, “Effects of cosmological moduli fields on cosmic microwave background,” Phys. Lett. B 522 (2001) 215–221, arXiv:hep-ph/0110096. [Erratum: Phys.Lett.B 539, 303–303 (2002)].
  * [11] D. H. Lyth, C. Ungarelli, and D. Wands, “The Primordial density perturbation in the curvaton scenario,” Phys. Rev. D 67 (2003) 023503, arXiv:astro-ph/0208055.
  * [12] N. Bartolo, S. Matarrese, and A. Riotto, “On nonGaussianity in the curvaton scenario,” Phys. Rev. D 69 (2004) 043503, arXiv:hep-ph/0309033.
  * [13] M. Sasaki, J. Valiviita, and D. Wands, “Non-Gaussianity of the primordial perturbation in the curvaton model,” Phys. Rev. D 74 (2006) 103003, arXiv:astro-ph/0607627.
  * [14] N. Bartolo, S. Matarrese, and A. Riotto, “Nongaussianity from inflation,” Phys. Rev. D 65 (2002) 103505, arXiv:hep-ph/0112261.
  * [15] M. Kawasaki, K. Nakayama, T. Sekiguchi, T. Suyama, and F. Takahashi, “Non-Gaussianity from isocurvature perturbations,” JCAP 11 (2008) 019, arXiv:0808.0009 [astro-ph].
  * [16] D. Langlois, F. Vernizzi, and D. Wands, “Non-linear isocurvature perturbations and non-Gaussianities,” JCAP 12 (2008) 004, arXiv:0809.4646 [astro-ph].
  * [17] C. Hikage, K. Koyama, T. Matsubara, T. Takahashi, and M. Yamaguchi, “Limits on Isocurvature Perturbations from Non-Gaussianity in WMAP Temperature Anisotropy,” Mon. Not. Roy. Astron. Soc. 398 (2009) 2188–2198, arXiv:0812.3500 [astro-ph].
  * [18] M. Kawasaki, K. Nakayama, T. Sekiguchi, T. Suyama, and F. Takahashi, “A General Analysis of Non-Gaussianity from Isocurvature Perturbations,” JCAP 01 (2009) 042, arXiv:0810.0208 [astro-ph].
  * [19] D. Langlois and B. van Tent, “Hunting for Isocurvature Modes in the CMB non-Gaussianities,” Class. Quant. Grav. 28 (2011) 222001, arXiv:1104.2567 [astro-ph.CO].
  * [20] C. Hikage, M. Kawasaki, T. Sekiguchi, and T. Takahashi, “Extended analysis of CMB constraints on non-Gaussianity in isocurvature perturbations,” JCAP 03 (2013) 020, arXiv:1212.6001 [astro-ph.CO].
  * [21] D. Langlois and B. van Tent, “Isocurvature modes in the CMB bispectrum,” JCAP 07 (2012) 040, arXiv:1204.5042 [astro-ph.CO].
  * [22] T. Montandon, G. Patanchon, and B. van Tent, “Isocurvature modes: joint analysis of the CMB power spectrum and bispectrum,” JCAP 01 (2021) 004, arXiv:2007.05457 [astro-ph.CO].
  * [23] P. Creminelli, “On non-Gaussianities in single-field inflation,” JCAP 10 (2003) 003, arXiv:astro-ph/0306122.
  * [24] P. Creminelli, A. Nicolis, L. Senatore, M. Tegmark, and M. Zaldarriaga, “Limits on non-gaussianities from wmap data,” JCAP 05 (2006) 004, arXiv:astro-ph/0509029.
  * [25] L. Senatore, K. M. Smith, and M. Zaldarriaga, “Non-Gaussianities in Single Field Inflation and their Optimal Limits from the WMAP 5-year Data,” JCAP 01 (2010) 028, arXiv:0905.3746 [astro-ph.CO].
  * [26] D. Babich, P. Creminelli, and M. Zaldarriaga, “The Shape of non-Gaussianities,” JCAP 08 (2004) 009, arXiv:astro-ph/0405356.
  * [27] J. M. Bardeen, “Gauge Invariant Cosmological Perturbations,” Phys. Rev. D 22 (1980) 1882–1905.
  * [28] C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, “The Effective Field Theory of Inflation,” JHEP 03 (2008) 014, arXiv:0709.0293 [hep-th].
  * [29] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [30] R. D. Peccei and H. R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D16 (1977) 1791–1797.
  * [31] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * [32] F. Wilczek, “Problem of Strong p and t Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [33] J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. 120B (1983) 127–132.
  * [34] L. F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. 120B (1983) 133–136.
  * [35] M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. 120B (1983) 137–141.
  * [36] K. A. Malik and D. Wands, “Cosmological perturbations,” Phys. Rept. 475 (2009) 1–51, arXiv:0809.4944 [astro-ph].
  * [37] X. Chen, J. Fan, and L. Li, “New inflationary probes of axion dark matter,” JHEP 12 (2023) 197, arXiv:2303.03406 [hep-ph].
  * [38] X. Chen and Y. Wang, “Quasi-Single Field Inflation and Non-Gaussianities,” JCAP 04 (2010) 027, arXiv:0911.3380 [hep-th].
  * [39] N. Arkani-Hamed and J. Maldacena, “Cosmological Collider Physics,” arXiv:1503.08043 [hep-th].
  * [40] M. Freytsis, S. Kumar, G. N. Remmen, and N. L. Rodd, “Multifield positivity bounds for inflation,” JHEP 09 (2023) 041, arXiv:2210.10791 [hep-th].
  * [41] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020) A10, arXiv:1807.06211 [astro-ph.CO].
  * [42] M. Buschmann, C. Dessert, J. W. Foster, A. J. Long, and B. R. Safdi, “Upper Limit on the QCD Axion Mass from Isolated Neutron Star Cooling,” Phys. Rev. Lett. 128 no. 9, (2022) 091102, arXiv:2111.09892 [hep-ph].
  * [43] Particle Data Group Collaboration, R. L. Workman and Others, “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [44] L. Senatore and M. Zaldarriaga, “The Effective Field Theory of Multifield Inflation,” JHEP 04 (2012) 024, arXiv:1009.2093 [hep-th].
  * [45] A. S. Chou et al., “Snowmass Cosmic Frontier Report,” in Snowmass 2021.  11, 2022.  arXiv:2211.09978 [hep-ex].


# Shorter Pulses, Smaller Errors

Quantum Circuit Optimization via Parameterized Pulses

Finn Voichick  0000-0002-1913-4178<EMAIL_ADDRESS>University of MarylandCollege
ParkUnited States ,  Leonidas Lampropoulos  0000-0003-0269-9815
<EMAIL_ADDRESS>University of MarylandCollege ParkUnited States and  Robert
Rand  0000-0001-6842-5505 University of ChicagoUnited States<EMAIL_ADDRESS>
(2023)

###### Abstract.

We propose a technique for optimizing quantum programs by temporally
stretching pre-calibrated pulses. As an example, we modify a three-qubit
Toffoli gate implementation by using an off-the-shelf numerical optimization
algorithm to shorten the cross-resonance pulses in the sequence. Preliminary
quantum process tomography results suggest that our strategy sometimes halves
a Toffoli gate’s error in practice, increasing process fidelity from around
60% to around 80%. Unlike existing quantum control techniques, ours takes
seconds to converge, demonstrating its potential utility when incorporated
into a general-purpose compiler pass that improves both the time and the
accuracy of quantum programs.

quantum process tomography, analog quantum computing, cross-resonance pulse

††copyright: none††conference: 4th International Workshop on Programming
Languages for Quantum Computing; 2024; London††price: ††ccs: Hardware Quantum
computation††ccs: Software and its engineering Compilers††ccs: Theory of
computation Timed and hybrid models††ccs: Theory of computation Quantum
information theory

##  1\. Introduction

Quantum programs rely on multiple kinds of optimizations. At one end of the
spectrum, computer scientists design optimizing compilers for high-level
quantum circuits, treating a program as a long sequence of discrete gates and
incrementally replacing subsequences with smaller ones. At the other end,
physicists use experiments, computational models, and numerical algorithms to
engineer and optimize the low-level pulses that make up a quantum computer’s
base gate set. Both strategies can improve a quantum program’s performance,
but the two work on different scales and use different techniques. This
compartmentalization creates a gap: high-level optimizers view individual
gates as atomic entities and overlook the impact of noise, while low-level
quantum control tends to be slow and specialized for a small set of gates.

This work draws attention to the opportunities in the middle ground of this
optimization spectrum. Though precisely modeling a quantum system’s dynamics
can be a computationally expensive task, it is possible to reasonably
approximate the impact of decoherence and to take this noise into account in a
general-purpose optimizing compiler without sacrificing performance.

As the title suggests, our main observation is that a hardware pulse of
shorter duration tends to have less error, as there is less opportunity for
cross-talk and decoherence. While there is a tendency to view a quantum
computer’s operations as atomic “logic gates,” the underlying microwave pulses
that implement these gates inhabit a continuous parameter space. Shortening a
pulse’s duration changes its behavior in a predictable way, and a circuit
optimizer can leverage this fact to produce higher-fidelity programs. We rely
on a quasi-Newton method to find parameters that shorten pulses while still
implementing the target unitary. Such numerical optimization techniques are
uncommon in compilers, but this is an unusual use case where they apply –
changing a program’s real gate parameters changes the output, and we use a
continuous and differentiable function to model these effects.

Though this work is still preliminary and further experiments are needed, we
believe that our techniques show promise for a general-purpose compiler
optimization pass that is noise-aware and avoids the computational expense of
full quantum optimal control. The rest of this technical report makes the
following contributions:

  * •

A description and Python implementation of a new technique for optimizing
quantum circuits, consisting of two main components:

    * –

An extension of ibm’s existing gate set. We temporally stretch the pre-
calibrated cross-resonance pulses, creating an additional variable-length
pulse.

    * –

Software that optimizes a quantum circuit to take advantage of these
parameterized pulses, generating shorter pulse sequences without significantly
altering circuit semantics.

  * •

A case study demonstrating that this technique can quickly converge to a pulse
sequence that improves a quantum Toffoli gate implementation. Experimental
quantum process tomography data collected from ibm’s superconducting machines
put the process fidelity of the original implementation below 60%, while the
fidelity of our modified version exceeds 80%.

  * •

A discussion of integrating of our technique into a general-purpose compiler
pass, which could improve performance on large-scale quantum circuits.

##  2\. Technique

Maintainers of quantum hardware provide a set of instructions available on
their machines, typically a gate set universal for quantum computation. For
example, some of Ibm’s machines use the gates in Fig 1. On quantum devices,
multiqubit gates require more time and thus produce more error than their
single-qubit counterparts (Sheldon et al., 2016; Kirchhoff et al., 2018).
Using Ibm’s superconducting qubits, single-qubit gate errors are on the order
of $10^{-4}$, while two-qubit gate errors are on the order of $10^{-2}$.

$\displaystyle R_{z}(\theta)$
$\displaystyle\defeq\begin{bmatrix}e^{-i\theta/2}&0\\\
0&e^{i\theta/2}\end{bmatrix}$ $\displaystyle X$
$\displaystyle\defeq\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}$
$\displaystyle\sqrt{X}$
$\displaystyle\defeq\frac{1}{\sqrt{2}}\begin{bmatrix}1&-i\\\
-i&1\end{bmatrix}$ cnot $\displaystyle\defeq\begin{bmatrix}1&0&0&0\\\
0&1&0&0\\\ 0&0&0&1\\\ 0&0&1&0\end{bmatrix}$ Figure 1. Basis gates for some of
Ibm’s superconducting quantum computers

Our main insight is that shortening the two-qubit pulses can reduce error in
quantum circuits. By extending the gate set beyond the fixed-duration cnot
gate, execution involves less time performing two-qubit operations. This
technique requires an alternative compilation method that takes advantage of
parameterized gates. Our implementation is available in a public git
repository (Voichick et al., 2023).

$\displaystyle\operatorname{eCR}(\theta)$
$\displaystyle\defeq\begin{bmatrix}0&0&\cos\theta&i\sin\theta\\\
0&0&i\sin\theta&\cos\theta\\\ \cos\theta&-i\sin\theta&0&0\\\
-i\sin\theta&\cos\theta&0&0\end{bmatrix}$
$\displaystyle\operatorname{eCR}\left(\frac{\pi}{4}\right)$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}0&0&1&i\\\ 0&0&i&1\\\
1&-i&0&0\\\ -i&1&0&0\end{bmatrix}$ Figure 2. The echo-cross-resonance gate
matrix, maximally entangling when $\theta=\frac{\pi}{4}$ Figure 3. A pre-
calibrated $\operatorname{eCR}(\frac{\pi}{4})$ gate on device ibm_lagos,
qubits 0 and 1 Figure 4. A custom $\operatorname{eCR}(\frac{\pi}{8})$ gate on
ibm_lagos, qubits 0 and 1

###  2.1. Pulse Stretching

Ibm recalibrates all of the pulse sequences on their quantum machines daily.
To implement multi-qubit (entangling) operations, they use an “echo-cross-
resonance” (eCR) pulse sequence, based on the $Z\otimes X$ Hamiltonian
interaction (Córcoles et al., 2013). The relationship between pulse amplitude
and Hamiltonian is nonlinear, and the eCR sequence mitigates the effect of
extraneous terms such as $Z\otimes I$ and $I\otimes X$ (Malekakhlagh et al.,
2020). Consisting of two opposing cross-resonance (CR) pulses sandwiching a
Pauli-X gate on the control qubit, this sequence can approximate the family
$\operatorname{eCR}(\theta)$ of unitaries in Fig. 2, and ibm regularly
recalibrates the $\operatorname{eCR}\left(\frac{\pi}{4}\right)$ pulse, which
can implement a cnot gate when composed with single-qubit gates (Krantz et
al., 2019). By customizing the calibrated CR pulses to be shorter, we obtain
the full family of $\operatorname{eCR}$ gates. This strategy is similar in
principle to that of Gokhale et al. 2020: By bootstrapping from existing
calibrations, we avoid a potentially costly closed-loop tune-up process
(Kirchhoff et al., 2018).

As an example, Fig. 3 depicts the pulse sequence for an
$\operatorname{eCR}(\frac{\pi}{4})$ gate. In arbitrary units, it plots the
amplitudes of the microwave pulses applied to three channels over the course
of around five hundred nanoseconds. Channel u0 is the coupler between qubits 0
and 1, the d1 pulse is an active cancellation pulse on qubit 1 (Sheldon et
al., 2016), and the d0 pulse is an intermediate Pauli-X gate applied to qubit
0. Fig. 4 depicts a smaller rotation with higher fidelity. To generate this
pulse envelope, we leave the Pauli-X gate unchanged but shorten cross-
resonance segments, treating the angle of rotation as roughly proportional to
the area under the curve. We maintain the same shape for the rising and
falling of the pulse amplitude because speeding up this increase and decrease
can lead to leakage into higher energy levels.

###  2.2. Optimization

To synthesize a quantum circuit implementing a desired unitary, we employ an
off-the-shelf quasi-Newton optimization algorithm, the Bfgs method (Zhu et
al., 1997). Our model and cost function are considerably simpler than those
typically seen in quantum optimal control (Koch et al., 2022), which makes it
less precise but scalable to larger circuits. We rely only on the error
metrics that ibm already computes as part of its daily calibration routine.

We assume depolarizing noise that depends on pulse length. For simplicity, we
restrict attention to the dominant Kraus operator in the resulting
superoperator, treating a gate as a norm-decreasing operator rather than a
(norm-preserving) unitary one. Consequently, if the
$\operatorname{eCR}(\theta)$ matrix defined in the previous section describes
the ideal unitary, our noise model treats it instead as
$(1-\varepsilon)^{\sin^{2}\theta}\operatorname{eCR}(\theta)$, where
$\varepsilon$ is the estimated cnot error rate. We use the $\sin^{2}\theta$
exponent not based on rigorous experimental data, but simply because it is
periodic and differentiable everywhere, and it works in practice. A maximally
entangling
$\theta\in\\{\ldots,-\frac{\pi}{4},\frac{\pi}{4},\frac{3\pi}{4},\ldots\\}$
incurs a maximum error of $\varepsilon$, while a zero-duration
$\theta\in\\{\ldots,0,\frac{\pi}{2},\ldots\\}$ requires only single-qubit
gates with comparably negligible error. The $\varepsilon$ parameter can vary
across different qubit pairs, enabling the optimizer to base its pulse
durations on the noisiness of the underlying hardware.

[column sep=2mm,row sep=2mm] & U(θ_1, ϕ_2, λ_3) [2]eCR(θ_7) U(θ_8, ϕ_9, λ_10)
[2]eCR(θ_14) U(θ_15, ϕ_16, λ_17)
U(θ_4, ϕ_5, λ_6) U(θ_11, ϕ_12, λ_13) U(θ_18, ϕ_19, λ_20)

Figure 5. A two-qubit circuit skeleton with twenty parameters

Using Qiskit’s built-in basis gate conversion procedures, we begin by
translating an input circuit into the gate set
$\\{R_{z},R_{x},\operatorname{eCR}(\pi/4)\\}$. We fix the resulting circuit
skeleton so that the circuit semantics is defined only by a vector of real
parameters, such as the twenty parameters in Fig. 5. If the (norm-decreasing)
semantics matrix of the resulting circuit is $V$ and the target unitary is
$U$, we use a cost function of $1-|\operatorname{trace}(V^{\dagger}U)|^{2}$.
Because this cost function applies only sinusoidal and polynomial functions to
the input parameters $\theta$, we can quickly compute its Jacobian, making the
optimization algorithm more efficient.

##  3\. Evaluation

We applied this technique to optimize a Toffoli gate. Under the assumption of
all-to-all connectivity, a Toffoli gate decomposition requires six cnot gates
(Shende and Markov, 2009), but superconducting qubits tend to have limited
connectivity (Hazra et al., 2021), so we assume that the three qubits are
arranged linearly with each connected only to its nearest neighbors. Using
Qiskit’s libraries (Wille et al., 2019), we arrived at an implementation with
ten entangling gates, where the final three cnot gates implement a swap gate
used to correct for Qiskit’s discrepency between physical and virtual qubits.
This is not the optimal circuit for this gate set and coupling map, but we
were unable to improve upon this circuit using several existing optimizers (Xu
et al., 2023; Quetschlich et al., 2023; Sivarajah et al., 2020), so it seems a
reasonable point of comparison.

[column sep=2mm,row sep=2mm] & 1 U 1 U 1 U 1
U U U U U U U U -1 -1
-1 -1 U -1 U -1 U

Figure 6. A skeleton of the ten-cnot Toffoli gate used as input for our
optimizer

Running on a low-end consumer laptop, our optimization routine takes a few
seconds to converge to a locally optimal solution. The resulting program uses
ten corresponding eCR pulse sequences of varying duration, some with the full
$\frac{\pi}{4}$ rotation, but most markedly shorter.

To compare the original ten-cnot Toffoli implementation to our version with
shorter pulses, we use quantum process tomography (qpt), a standard procedure
for characterizing a quantum process by preparing and measuring the qubits in
a variety of configurations (O’Brien et al., 2004). Thus far, we have
performed full qpt on three of Ibm’s superconducting machines. On device
ibm_lagos, we used the physical qubits indexed 0, 1, and 3, with qubit 1 as
the target. On ibm_nairobi, we used qubits 0, 1, and 2, with qubit 1 as the
target. On ibm_perth, we used qubits 1, 2, and 3, with qubit 1 as the target.
Table 1 shows the qpt results.

Table 1. Toffoli gate fidelity comparison | Process fidelity
---|---
Machine |  Cnot-based | Shortened eCR
ibm_lagos | 49% | 82%
ibm_nairobi | 78% | 80%
ibm_perth | 64% | 78%

##  4\. Prospects

This project is a work in progress, but we believe that our preliminary data
shows promise for a general-purpose compiler pass. Though it would not be
feasible to optimize a large-scale circuit in its entirety using this method,
it would be feasible as a peephole optimization for subcircuits below a
tractable size. Our current strategy uses the input circuit’s existing
structure of multiqubit gates, but future work could explore modifying this
skeleton by pruning short-duration gates or adding more entangling gates. We
could also fine-tune our error model further, perhaps changing the
$\sin^{2}(\theta)$ exponent to something more physically accurate.

Measuring the performance of this optimizer on large-scale programs could be a
challenge. Researchers often compare optimizers based on gate count (Xu et
al., 2023), but our strategy improves circuit fidelity without affecting gate
count at all. Analytical approximations of noise often differ from
experimental results, and full qpt requires a number of circuits exponential
in the number of qubits, which prohibits any large-scale use. One option is to
restrict attention to a class of circuits that is classically simulable, but
this is a bit limited for testing a general-purpose optimizer; another is to
use a novel analytical metric such as nanoseconds spent performing cross-
resonance pulses, but this would require more experiments to establish its
reliability.

The strategies that we present here apply to programs written in the “digital”
quantum circuit model. While it draws inspiration from analog quantum
computing for the pulse-level control, optimized circuits performing
Hamiltonian simulation will still underperform relative to fully analog
simulation, but it is our hope that these optimizations can make near-term
circuit-based algorithms more feasible.

More generally, we argue that quantum compiler optimizations should not view
quantum gates as atomic and perfect. While there is a tendency to treat
quantum circuit optimizations like any other compiler optimization, the
continuously-parameterized nature of quantum operations means that there are
optimizations that traditional methods cannot discover. This is an opportunity
for programming language researchers to draw on well-developed algorithms from
other kinds of computational optimization, such as machine learning and
numerical analysis.

## References

  * (1)
  * Córcoles et al. (2013) A. D. Córcoles, Jay M. Gambetta, Jerry M. Chow, John A. Smolin, Matthew Ware, Joel Strand, B. L. T. Plourde, and M. Steffen. 2013.  Process Verification of Two-Qubit Quantum Gates by Randomized Benchmarking.  _Phys. Rev. A_ 87 (March 2013), 030301.  Issue 3. https://doi.org/10.1103/PhysRevA.87.030301 arXiv:1210.7011 [quant-ph]
  * Gokhale et al. (2020) Pranav Gokhale, Ali Javadi-Abhari, Nathan Earnest, Yunong Shi, and Frederic T. Chong. 2020.  Optimized Quantum Compilation for Near-Term Algorithms with OpenPulse. In _53rd Annual IEEE/ACM International Symposium on Microarchitecture_ (Athens). IEEE, New York, 186–200.  https://doi.org/10.1109/MICRO50266.2020.00027 arXiv:2004.11205 [quant-ph]
  * Hazra et al. (2021) Sumeru Hazra, Anirban Bhattacharjee, Madhavi Chand, Kishor V. Salunkhe, Sriram Gopalakrishnan, Meghan P. Patankar, and R. Vijay. 2021.  Ring-Resonator-Based Coupling Architecture for Enhanced Connectivity in a Superconducting Multiqubit Network.  _Phys. Rev. Appl._ 16 (Aug. 2021), 024018.  Issue 2. https://doi.org/10.1103/PhysRevApplied.16.024018
  * Kirchhoff et al. (2018) Susanna Kirchhoff, Torsten Keßler, Per J. Liebermann, Elie Assémat, Shai Machnes, Felix Motzoi, and Frank K. Wilhelm. 2018.  Optimized Cross-Resonance Gate for Coupled Transmon Systems.  _Phys. Rev. A_ 97 (April 2018), 042348.  Issue 4. https://doi.org/10.1103/PhysRevA.97.042348 arXiv:1701.01841 [quant-ph]
  * Koch et al. (2022) Christiane P. Koch, Ugo Boscain, Tommaso Calarco, Gunther Dirr, Stefan Filipp, Steffen J. Glaser, Ronnie Kosloff, Simone Montangero, Thomas Schulte-Herbrüggen, Dominique Sugny, and Frank K. Wilhelm. 2022.  Quantum Optimal Control in Quantum Technologies: Strategic Report on Current Status, Visions and Goals for Research in Europe.  _EPJ Quantum Technology_ 9, 1 (July 2022), 60 pages.  https://doi.org/10.1140/epjqt/s40507-022-00138-x arXiv:2205.12110 [quant-ph]
  * Krantz et al. (2019) P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver. 2019.  A Quantum Engineer’s Guide to Superconducting Qubits.  _Applied Physics Reviews_ 6, 2 (June 2019), 021318.  https://doi.org/10.1063/1.5089550 arXiv:1904.06560 [quant-ph]
  * Malekakhlagh et al. (2020) Moein Malekakhlagh, Easwar Magesan, and David C. McKay. 2020.  First-Principles Analysis of Cross-Resonance Gate Operation.  _Phys. Rev. A_ 102 (Oct. 2020), 042605.  Issue 4. https://doi.org/10.1103/PhysRevA.102.042605 arXiv:2005.00133 [quant-ph]
  * O’Brien et al. (2004) J. L. O’Brien, G. J. Pryde, A. Gilchrist, D. F. V. James, N. K. Langford, T. C. Ralph, and A. G. White. 2004.  Quantum Process Tomography of a Controlled-NOT Gate.  _Phys. Rev. Lett._ 93 (Aug. 2004), 080502.  Issue 8. https://doi.org/10.1103/PhysRevLett.93.080502 arXiv:quant-ph/0402166
  * Quetschlich et al. (2023) Nils Quetschlich, Lukas Burgholzer, and Robert Wille. 2023.  Compiler Optimization for Quantum Computing Using Reinforcement Learning. In _2023 60th ACM/IEEE Design Automation Conference (DAC)_. IEEE, New York, 1–6.  https://doi.org/10.1109/DAC56929.2023.10248002 arXiv:2212.04508 [quant-ph]
  * Sheldon et al. (2016) Sarah Sheldon, Easwar Magesan, Jerry M. Chow, and Jay M. Gambetta. 2016.  Procedure for Systematically Tuning up Cross-Talk in the Cross-Resonance Gate.  _Phys. Rev. A_ 93 (June 2016), 060302.  Issue 6. https://doi.org/10.1103/PhysRevA.93.060302 arXiv:1603.04821 [quant-ph]
  * Shende and Markov (2009) Vivek V. Shende and Igor L. Markov. 2009.  On the CNOT-Cost of Toffoli Gates.  _Quantum Info. Comput._ 9, 5 (May 2009), 461–486.  arXiv:0803.2316 [quant-ph]
  * Sivarajah et al. (2020) Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington, and Ross Duncan. 2020.  t—ket⟩: a retargetable compiler for NISQ devices.  _Quantum Science and Technology_ 6, 1 (Nov. 2020), 014003.  https://doi.org/10.1088/2058-9565/ab8e92 arXiv:2003.10611 [quant-ph]
  * Voichick et al. (2023) Finn Voichick, Leonidas Lampropoulos, and Robert Rand. 2023.  Shorter Pulses code repository.  https://gitlab.umiacs.umd.edu/finn/shorter-pulses
  * Wille et al. (2019) Robert Wille, Rod Van Meter, and Yehuda Naveh. 2019.  IBM’s Qiskit Tool Chain: Working with and Developing for Real Quantum Computers. In _2019 Design, Automation and Test in Europe Conference and Exhibition_. IEEE, New York, 1234–1240.  https://doi.org/10.23919/DATE.2019.8715261
  * Xu et al. (2023) Amanda Xu, Abtin Molavi, Lauren Pick, Swamit Tannu, and Aws Albarghouthi. 2023.  Synthesizing Quantum-Circuit Optimizers.  _Proc. ACM Program. Lang._ 7, PLDI, Article 140 (June 2023), 25 pages.  https://doi.org/10.1145/3591254 arXiv:2211.09691 [cs.PL]
  * Zhu et al. (1997) Ciyou Zhu, Richard H. Byrd, Peihuang Lu, and Jorge Nocedal. 1997.  Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-Scale Bound-Constrained Optimization.  _ACM Trans. Math. Softw._ 23, 4 (Dec. 1997), 550–560.  https://doi.org/10.1145/279232.279236