\section{PC in single-field inflation}
We first review the effects of a PC field $\sigma$ on adiabatic density perturbations seeded by the quantum fluctuations $\delta \phi$ of the inflaton field, as proposed in \cite{Chen:2011zf}.
The initial condition for the classical motion is that it is restricted to slow-roll in the $\phi$ direction of field space, while the massive $\sigma$ field is stuck at its origin. It is assumed that there is a feature in the scalar potential that leads to a sudden transfer of energy to excite $\sigma$ at time $t_0$, after which $\phi$, however, continues to slow-roll.
The classical dynamics of $\sigma$ after $t_0$ is governed by the homogeneous equation of motion,
\begin{align}
    \partial^{2}_{t}\sigma + 3H\, \partial_{t} \sigma+ M^2 \sigma =0 \, ,
\end{align}
where $M$ is the mass of the $\sigma$ field.
In the limit $M \gg H$, the solution has an oscillatory form,
\begin{align}\label{eq:sig_oscil}
    \sigma_{\rm osc}(t)  \simeq  \theta(t-t_0) \, \sigma_{0}
     \left(\frac{a(t)}{a(t_0)}\right)^{-3/2} \cos[M (t-t_0)  ] \,.
\end{align}
The oscillating field energy in an inflationary background imparts a small oscillatory correction to the Hubble constant,
\begin{align}\label{eq:Hosc}
        \Delta H_{\rm osc} \approx \theta(t-t_0) \frac{\sigma_{0}^2 M}{8 M_{\rm Pl}^2} \left(\frac{a(t)}{a(t_0)}\right)^{-3} \cos[2M (t-t_0)] \, ,
\end{align}
where $M_{\rm Pl}$ is the Planck scale.
This in turn imparts an oscillatory component to the slow-roll parameters, $\Delta \epsilon_{\rm osc} = -\partial_{t} (\Delta H_{\rm osc})/H^{2}$ and $\Delta \eta_{\rm osc} \approx \partial_{t} (\Delta  \epsilon_{\rm osc})/(H \epsilon)$. Here $H$ and $\epsilon$ are the unperturbed parameters in the absence of $\sigma$ oscillations.

These oscillatory components then introduce scale-invariance-breaking features in the power spectrum of adiabatic perturbations, ${\cal P}_{\zeta_{\phi}}$, defined as
\begin{align}
    \< \zeta_{\phi}(\vec{k}) \zeta_{\phi}(\vec{k}') \> \equiv (2\pi)^3 \delta(\vec{k}+\vec{k}') {\cal P}_{\zeta_{\phi}}(k) \, ,
\end{align}
where $\zeta_{\phi}$ are the gauge-invariant perturbations capturing $\delta \phi$, as measured in the primordial fluctuations of the CMB and LSS.
This \textit{irreducible} contribution is predominantly given by
\begin{align}\label{eq:cmbSIBgrav}
    \left(\frac{\Delta {\cal P}_{\zeta_{\phi}}}{{\cal P}^{(0)}_{\zeta_{\phi}}}\right)_{\Delta H_{\rm osc}}  \simeq & \,\, \theta(k-k_g) \frac{\sqrt{\pi}}{2} \frac{( M^{2} \sigma_{0}^{2})}{ \,\dot{\phi}_{0}^2} \left(\frac{M}{H}\right)^{1/2} \nonumber \\
    & \quad \times \left( \frac{k}{k_g}\right)^{-3} \cos\left[ \frac{2M}{H} \log \frac{k}{k_g}\right] \,.
\end{align}
where ${\cal P}^{(0)} $ denotes the unperturbed scale-invariant spectrum $\propto 1/k^3$.
We can understand the key features of this  calculation intuitively, by noting that ordinarily ${\cal P}_{\zeta_{\phi}} \propto 1/\epsilon$, so that we expect roughly
$\Delta {\cal P}_{\zeta_{\phi}}/{\cal P}^{(0)}_{\zeta_{\phi}} \propto \Delta \epsilon_{\rm osc}/\epsilon$.
Now, $\Delta \epsilon_{\rm osc}(t)$ represents a high-frequency ($2M$ by Eq.~(\ref{eq:Hosc})) gravitational source for production of inflaton pairs, so their physical momentum on production, $k e^{-Ht} \approx M$, by the local Minkowski energy conservation.
The softest co-moving momentum $k_g$ produced from this gravitational source is given by $k_g  e^{-H t_0} \approx M$. Solving for $t, t_0$ in terms of $k, k_g$ results in the dependence on the latter seen in Eq.~(\ref{eq:cmbSIBgrav}).

Similarly, $\Delta \eta_{\rm osc} $ gives the dominant irreducible PC contribution to the 3-point correlators, which can be written in a following normalized form:
\begin{align}
 F  \equiv \frac{5}{6} \frac{\<\zeta_{\phi}(\vec{k}_1) \zeta_{\phi}(\vec{k}_2) \zeta_{\phi}(\vec{k}_3)\>}{{\cal P}_{\zeta_{\phi}}(k_1) {\cal P}_{\zeta_{\phi}}(k_2) + (k_1 \rightarrow k_3) + (k_2 \rightarrow k_3)}\,.
\end{align}
In the equilateral triangle configuration $k_{1} \sim k_2 \sim k_3$,
\begin{align}\label{eq:cmbSIB3pnt}
F_{\Delta H_{\rm osc}} &\simeq  \theta(K-2k_g)  \, f_{\rm NL}^{\rm PC} \nonumber \\
&\qquad \times \left( \frac{K}{2 k_g}\right)^{-3} \cos\left[2 \frac{M}{H} \log \frac{K}{2 k_g}\right] \,,
\end{align}
where $K=\sum_{i=1}^{3} |\vec{k}_i|$, and the amplitude of NG can be evaluated to be $f_{\rm NL}^{\rm PC} \sim \frac{10}{9} \frac{\sqrt{\pi}}{4} \frac{(M^{2} \sigma_{0}^{2})}{\dot{\phi}_{0}^2} \left(\frac{M}{H}\right)^{5/2}$.

Apart from these irreducible contributions, the coherent $\sigma$ field may also directly decay to inflaton fluctuations.
Consider the lowest dimensional interaction  consistent with an (approximate) inflaton shift-symmetry, ${\cal L}_{\rm int} \sim (\sigma \, \partial_{\mu} \phi \partial^{\mu} \phi)/\Lambda_{\phi} $.
The PC contribution to ${\cal P}_{\zeta_{\phi}}$ from the decay of $\sigma_{\rm osc}$ through this interaction can be evaluated using the in-in formalism \cite{Chen:2014cwa, chen2022classical} (see \cite{Weinberg:2005vy} for a review of the in-in formalism),
\begin{align}\label{eq:cmbSIBdec}
   \left( \frac{\Delta {\cal P}_{\zeta_{\phi}}}{{\cal P}^{(0)}_{\zeta_{\phi}}}\right)_{\rm dec} \simeq & \,\, \theta(k - k_d)  \sqrt{2 \pi} \frac{ \sigma_0}{\Lambda_{\phi}}  \left(\frac{M}{H}\right)^{1/2} \nonumber\\
   & \quad \times \left( \frac{k}{k_d}\right)^{-3/2} \cos\left[\frac{M}{H} \log \frac{k}{k_d}\right] \,.
\end{align}
Here, $\sigma_{\rm osc}$ directly acts as a high-frequency ($M$ from Eq.~(\ref{eq:sig_oscil})) source for inflaton production instead of the irreducible gravitational source $\Delta H_{\rm osc}$.
Therefore, the $k$-dependence can be obtained using a similar argument given for Eq.~(\ref{eq:cmbSIBgrav}).
With this source, a co-moving mode $k$ is produced when $k e^{-Ht} \approx M/2$. The softest such mode produced through direct decay is then $k_d e^{-Ht_0}=M/2 $.
Since the interaction is linear in $\sigma$,
the contribution in Eq.~(\ref{eq:cmbSIBdec}) matter-dilutes
 $\propto k^{-3/2}$, as compared to the $\propto k^{-3}$ dependence in \Eq{eq:cmbSIBgrav}.

Such PCs have been searched for in the Planck CMB data (see e.g. analysis in \cite{Chen:2014cwa,Braglia:2021rej}).
While there are potential candidates, none are statistically significant.
In the low to medium $k$-range, the constraints on the amplitudes of PC features in adiabatic perturbations are close to the ultimate cosmic variance limit.
While future probes like CMB polarization, LSS, and 21-cm tomography will improve sensitivity, they are typically  effective for PCs at high $k$ (small scales)  \cite{Braglia:2021rej, Chen:2016vvw, Chen:2016zuu}.
The PCs could be undetectable if $\sigma_0$ is  small, or, for the direct decay contribution, if $\Lambda_{\phi}$ is  large.
Considering the valuable information encoded in the PC features,
we ask if they can be probed by other cosmological observables.\footnote{
Refs.~\cite{Braglia:2020taf, Fumagalli:2020nvq} have also studied PC-like scale-invariance-breaking features in scalar-induced SGWB.
There are however important qualitative differences: 1) They consider features that appear at very small scales, which can not be probed by CMB as they are highly Silk-damped. In our work, PC occurs on scales at which the CMB is undamped, but are too small there for detection.
2) Their features appear in the {\it frequency} spectrum of SGWB, unlike in the anisotropies in our case.
}

\section{PC in multi-field inflation}
In the following, we study one such possibility where the excited $\sigma$ field decays dominantly to another light ``spectator'' field, $\chi$, present during inflation.
If  $m_{\chi}^{2} \ll H^2$, it will develop an (approximately) scale-invariant spectrum, ${\cal P}^{(0)}_{{\cal S}_{\chi}} \propto 1/k^3$, similar to that of the inflaton.
Here, ${\cal S}_{\chi} \equiv \zeta_{\chi}-\zeta_{\phi}$ is the isocurvature perturbation of $\chi$, independent of the fluctuations of $\zeta_{\phi}$ when the energy density in isocurvature is subdominant.
To naturally maintain the lightness of $\chi$, we again consider a shift-symmetric effective interaction of the form $\mathcal{L}_{\rm int} \sim  (\sigma \, \partial_{\mu}\chi \partial^{\mu} \chi)/\Lambda_{\chi}$, $\Lambda_{\chi} \ll \Lambda_{\phi}$. This gives a PC feature in $P_{{\cal S}_{\chi}}$ completely analogous to Eq.~(\ref{eq:cmbSIBdec}),\footnote{The contribution from $\Delta H_{\rm osc}$ is small, $(\Delta {\cal P}_{{\cal S}_{\chi}}/{\cal P}^{(0)}_{{\cal S}_{\chi}}) \propto (M^2 \sigma_{0}^2)/(H^2 M_{\rm Pl}^{2})$, and hence we only focus on the direct decay of $\sigma$ to $\delta \chi$.}
\begin{align}\label{eq:chiSIB}
   \frac{\Delta {\cal P}_{{\cal S}_{\chi}}}{{\cal P}^{(0)}_{{\cal S}_{\chi}}} &\simeq  \theta(k-k_d)  \alpha
    \left( \frac{k}{k_d}\right)^{-3/2} \cos\left[\frac{M}{H} \log \frac{k}{k_d}\right] \nonumber \\
    \text{where}\quad \alpha &= \sqrt{\frac{2 \pi M}{ H}} \frac{  \sigma_0}{\Lambda_{\chi}}\,.
\end{align}
There are some constraints on  $\Lambda_{\chi}$ to ensure effective field theory (EFT)  control. We require $ (\sigma_0/\Lambda_{\chi}) < 1$ to keep the correction to the kinetic energy of $\chi$ small, and $(\dot{\chi}_{0}/\Lambda_{\chi}^{2}),\, M/\Lambda_{\chi} < 1$ to ensure that no energy scales exceed the scale of non-renormalizability at which the EFT must break down.
The perturbative control of the in-in diagrams requires the PC feature to be small compared to the scale-invariant spectrum, giving a constraint $\alpha <1$.
We see qualitatively that it is possible to have small $\sigma_0$ and large $\Lambda_{\phi}$ so that PC feature is small in ${\cal P}_{\zeta_{\phi}}$, while $\sigma_0/\Lambda_{\chi}$ gives a sizeable PC contribution in ${\cal P}_{{\cal S}_{\chi}}$.

\section{Non-minimal reheating scenario}
We now illustrate how the primordial fluctuations of $\chi$, including sizeable PC, may be probed by a potentially
observable SGWB.
\begin{figure}[t]
	\centering
		\includegraphics[width=0.9\columnwidth]{ghost.pdf}
		\caption{A schematic of the model with modified post-inflationary history showing energy densities in various sectors as a function of the scale factor (time).}
		\label{fig:ghost}
\end{figure}
We will employ the cosmological model developed in Ref.~\cite{Geller:2018mwu} which showed how  the SGWB can have a large isocurvature component, making it uncorrelated with the CMB. We will add the new features of the $\sigma$ oscillations, back-reaction on the inflationary expansion and decays into $\chi$, perturbatively.
The model is shown schematically in \Fig{fig:ghost}.
After the end of inflation, $\chi$ reheats the Visible Sector (VS), which is presumed to undergo a strong first order PT, and whose remnants are  just the SM and dark matter (DM). If the PT occurs as the universe is passing through temperatures of order (multi-)TeV, this would result in a SGWB at frequencies observable at proposed GW detectors.
The inflaton on the other hand is taken to reheat a Hidden Sector of particle physics (HS) that interacts very weakly with the VS.
Some time after the PT, the HS particles decay to the VS, thermalizing with it before  Big Bang Nucleosynthesis. This leaves the SGWB as the only source of isocurvature.
The constraints from isocurvature and $N_{\rm eff}$ can be cast in terms of the ratio of energy density to critical density today as $\delta \Omega_{\rm GW} h^2 \lesssim 2\times 10^{-11}$ and  $\Omega_{\rm GW} h^2 \lesssim 2\times 10^{-6} $, respectively.
We will later show that our benchmarks easily satisfy these constraints.

The GW are sourced by  sub-horizon physics at the time of PT, and hence the fluctuations $\delta \rho_{\rm GW}$ on super-horizon scales are directly related to the initial fluctuations
in the VS after it is reheated by $\chi$ (but before the HS decays into it),
\begin{align}
    \delta_{\rm GW} \sim (\delta_{\rm VS})_{\rm initial} \sim \frac{ \delta \chi}{\chi_{0}} \,,
\end{align}
where $\delta \equiv \delta \rho / \rho$ and $\chi_{0}$ is the original homogeneous VEV of the $\chi$ field.
More precisely, the GW anisotropy can be written as \cite{Kumar:2021ffi, Malik:2008im}
\begin{align}\label{eq:delGW}
    \delta_{\rm GW} = -\frac{4}{3}\zeta_{\phi}+4 {\cal S}_{\chi} \left(1-\frac{4}{3}f_{\rm VS}\right) \approx 4 {\cal S}_{\chi} \,,
\end{align}
where we consider $f_{\rm VS} \ll 1$ and  ${\cal S}_{\chi} \gtrsim \zeta_{\phi}$.

The CMB anisotropies receive two independent contributions since the final fluctuations in the VS are a combination of the decay products of the HS and the initial VS,
\begin{align}
    \delta_{\gamma} \sim \left( \delta_{\rm VS} \right)_{\rm final} \sim  \left(\delta_{\rm HS}\right)_{\rm initial} + f_{\rm VS} \left(\delta_{\rm VS}\right)_{\rm initial} \,.
\end{align}
Here, $f_{\rm VS} \equiv (\rho_{\rm VS})_{\rm initial}/\rho_{\rm total}$.
In terms of the independent gauge-invariant fluctuations,
\begin{align}\label{eq:ghostCMBpert}
    \delta_{\gamma} \approx -\frac{4}{5} \zeta_{\phi} -\frac{4}{5}f_{\rm VS} {\cal S}_{\chi}.
\end{align}
From Eq.~(\ref{eq:ghostCMBpert}), we see that $f_{\rm VS} \ll 1$ can ensure a subdominant contribution of ${\cal S}_{\chi}$ to CMB anisotropies. The GW map in this case will be (mostly) uncorrelated with the CMB since the ${\cal S}_{\chi}$ is independent of $\zeta_{\phi}$.
An important subtlety of the multi-field reheating is that
the early Sachs-Wolfe (SW) effect suppresses CMB anisotropies but not the GW anisotropies, $(\delta_{\rm GW}/\delta_{\gamma}) \approx 5 ({\cal S}_{\chi}/\zeta_{\phi})$, as can be seen from Eqs.~(\ref{eq:delGW}) and (\ref{eq:ghostCMBpert}).
This is because the SW effect is {\it anti-correlated} with the dominant adiabatic density fluctuations from $\phi$, leading to reduced anisotropy in the observed CMB \cite{sachs1967perturbations}. On the other hand, the SW contribution, which is predominantly sourced by $\rho_{\rm HS}$, is {\it uncorrelated} with $\chi$ fluctuations, and hence there is no analogous cancellation in GW anisotropies.

Considering the high precision of CMB measurements, we must consider the PC feature induced by ${\cal S}_{\chi}$ into the CMB after the HS decay into the VS, Eq.~(\ref{eq:ghostCMBpert}), as well as
the irreducible PC feature,  Eqs.~(\ref{eq:cmbSIBgrav}) and (\ref{eq:cmbSIB3pnt}), despite their subdominance.
From Eq.~(\ref{eq:ghostCMBpert}), ${\cal S}_{\chi}$ induces
\begin{align}\label{eq:cmbSIBind}
    \left(\frac{\Delta {\cal P}_{\gamma}}{{\cal P}_{\gamma}} \right)_{\rm ind} \approx f^{2}_{\rm VS} \left(\frac{A_{\chi}}{A_{\phi}} \right) \left(\frac{\Delta {\cal P}_{{\cal S}_{\chi}}}{{\cal P}_{{\cal S}_{\chi}}} \right) \,,
\end{align}
where
${\cal P}^{(0)}_{j}(k) = 2 \pi^2 A_j k^{-3}$.
The amplitude $A_{\phi} \approx 2 \times 10^{-9}$ is fixed by the CMB.
An analogous relation holds for PC in LSS.

The anisotropies for a 2D map (CMB or PT-SGWB) are expressed in the multipole basis,
\begin{align}
    \langle \delta(\hat{n}) \delta(\hat{n}') \rangle
    =  \sum_{\ell} \frac{(2 \ell +1)}{4 \pi} C_{\ell} P_{\ell}(\cos \theta)\, ,
\end{align}
where  $\hat{n}\cdot\hat{n}' = \cos \theta$ and $P_{\ell}$ are the Legendre polynomials.
For higher multipoles $\ell \gtrsim O(100)$ (which will be of interest to us later), the expression for $C_{\ell}$ can be simplified using flat-sky approximation (see for example \cite{Munoz:2015eqa}):
\begin{align}\label{eq:flatskyCl}
    C_{\ell} \approx \frac{1}{r^{2}} \int \frac{d k_{\parallel}}{2\pi} {\cal P}_{\delta}\left(k_{\parallel}, k_{\perp} \sim \frac{\ell}{r} \right) \,,
\end{align}
where ${\cal P}_{\delta}$ is the power spectrum of the density perturbations, $k_{\parallel (\perp)}$ is the component of $\vec{k}$ parallel (perpendicular) to the line-of-sight, and  $r$ is the co-moving radius from us to the surface of the map that we see today.
The correlation can be separated into a nearly scale-invariant part and a small PC correction, $ C_{\ell} = C^{(0) }_{\ell} + \Delta C_{\ell}$, where $C^{(0)}_{\ell} [\ell (\ell+1)] \equiv A$ is approximately constant. From eqs.~(\ref{eq:delGW}) and (\ref{eq:flatskyCl}), $A_{\rm GW} \simeq 32 \pi A_{\chi}$.
The PC feature in SGWB can be computed using \Eqs{eq:chiSIB}{eq:flatskyCl},
\begin{align}\label{eq:chiSIBl}
    \frac{\Delta C^{\rm GW}_{\ell}}{C^{(0) {\rm GW}}_{\ell}} \approx & \,\, \theta(\ell-\ell_d)\,  \tilde{\alpha}  \left(\frac{\ell}{\ell_d}\right)^{-3/2} \cos\left[ \frac{M}{H} \log \frac{\ell}{\ell_d}\right]  \,,
\end{align}
where the amplitude of the PC feature in the multipole basis, $\tilde{\alpha} = \sqrt{2\pi H/M} \, \alpha$, is smaller than $\alpha$ due to the 2-dimensional nature the map.
$\ell_d \approx k_d r_{\rm GW}$ corresponds to the multipole where the PC feature appears in the GW map.
We can similarly evaluate the amplitude of PCs in CMB using eqs.~(\ref{eq:cmbSIBgrav}), (\ref{eq:cmbSIBdec}), and (\ref{eq:cmbSIBind}) instead.

With sufficient experimental sensitivity, the detection of the PC feature is ultimately limited by  cosmic variance.
If we search in the parameter regime $H \ell_d/M \gg 1$, then our signal oscillates slowly as a function of $\ell$, in which case we can mitigate cosmic variance by averaging the data over a bin of size $1 < \Delta \ell < \pi H \ell_d/M$ around each $\ell$.
Thus the signal adds coherently while the binned cosmic variance is
\begin{align}\label{eq:cvbin}
    \frac{\delta C_{\ell}^{bin}}{C_{\ell}} \approx  \sqrt{\frac{2}{(2\ell+1) \Delta \ell}} \,.
\end{align}
We determine the observability of the PC feature by comparing its largest amplitude (at $\ell_d$) with the smallest compatible  binned cosmic variance,
\begin{align}\label{eq:alphaObs}
     \tilde{\alpha} > \sqrt{\frac{2 (M/H)}{\pi (2\ell_d +1) \ell_d}} \, .
\end{align}
PCs with smaller amplitudes can be observable if they appears at larger $\ell_d$, as the cosmic variance drops. It is however experimentally challenging to probe high $\ell$ anisotropies. Keeping these issues in mind, we will consider $\ell_d \sim 100$ in this work.

We also want to estimate the strength of GW signal in our model.
During PT, processes like bubble collisions, sounds waves, and turbulence give rise to GW \cite{Caprini:2015zlo}.
After production, GW decouple and redshift as radiation until today.
An estimate for the peak energy density today as a fraction of critical energy density \cite{Hindmarsh:2013xza} is \footnote{We have only considered bubble wall collisions here, but the contribution from sound waves is generally comparable and even dominates at larger $\beta$ as it scales linearly with ($H_{\rm PT}/\beta$).}
\begin{align}\label{eq:gwMonopole}
    \Omega_{\rm GW} h^2 \sim 1.3 \times 10^{-6} \left(\frac{H_{\rm PT}}{\beta}\right)^2 f_{\rm VS}^{2} \,.
\end{align}
Here $\beta$ is the rate of bubble nucleation, $H_{\rm PT}$ is the Hubble expansion rate during PT, and we are considering a strong PT where most of the latent heat goes into accelerating the bubble walls.
The factor of $f_{\rm VS}$ appears since the PT is occurring in the VS sector that is subdominant in the energy density during PT.

The energy density fluctuation in the $\ell^{th}$ multipole is given by
\begin{align}\label{eq:gwAnisotropy}
    \delta \Omega_{\rm GW} h^2 (\ell) &\sim \Omega_{\rm GW} h^2 \times \frac{\sqrt{A_{\rm GW}}}{\ell}

\end{align}
The factor of $1/\ell$ appears because $C^{(0)}_{\ell} \propto [\ell (\ell +1)]^{-1}$ by scale-invariance.
The quantity of interest to us is the deviation from this scale-invariant spectrum due to the PC feature, $\Delta C_{\ell}^{\rm GW}$.
From Eq.~(\ref{eq:chiSIBl}),
\begin{align}\label{eq:gwSIB}
     \delta \Omega_{\rm GW}^{\rm PC} h^2 (\ell_d) \sim
    & 1.3 \times 10^{-6} \left(\frac{H_{\rm PT}}{\beta}\right)^2 f_{\rm VS}^{2} \frac{\sqrt{A_{\rm GW}}}{\ell_d} \left(\frac{\tilde{\alpha}}{2}\right).

\end{align}

\section{Benchmarks}
We will now illustrate our mechanism with two benchmarks.
Consider the following set of parameters,
\begin{align}\label{eq:benchmark-common}
    \sigma_0 \simeq 2H ,\quad \Lambda_{\chi} \simeq 40 H,\quad   M = 10H , \quad  \ell_d = 100

\end{align}
Let us first consider the PC signal in SGWB.
From eqs.~(\ref{eq:chiSIB}) and (\ref{eq:chiSIBl}), the amplitude of the PC feature is $\tilde{\alpha} = 0.31$.
This easily satisfies the condition of observability over cosmic variance, $\tilde{\alpha} > 0.018$, given in Eq.~(\ref{eq:alphaObs}).
A realization of SGWB map containing this PC feature and a random sample of cosmic variance is shown in Fig.~(\ref{fig:benchmark_plots}) (a) without binning and (b) with binning using $\Delta \ell = 6$. We see that even with a modest bin size, the observability of signal is significantly improved, with the signal-to-noise ratio $\sim 7.6$.

There are also irreducible PC corrections to the maps of adiabatic perturbations as discussed before.
Currently, the strongest constraints on these corrections come from the CMB \cite{Braglia:2021rej}.
Therefore, in our first benchmark, we consider a PC whose contributions in the 2D map of CMB are completely hidden under  cosmic variance, and yet are  observable in the SGWB.
In the future, 3D maps such as LSS and 21-cm are expected to achieve higher sensitivity as they can access more modes than a 2D map, and hence have lower cosmic variance.
It is then interesting to ask whether there are PCs which can not be detected even in a perfect 3D map of adiabatic perturbations, making isocurvature SGWB possibly the only channel of discovery.
We demonstrate this possibility in the second benchmark.

Consider the first benchmark along with Eq.~(\ref{eq:benchmark-common}),
\begin{align}
    (1) \quad f_{\rm VS} = 0.1, \, \frac{A_{\chi}}{A_{\phi}}  = 3 , \, \frac{\beta}{H_{\rm PT}} = 5 \,.
\end{align}
The PC corrections in CMB also appear around $\ell_{d}$ and $\ell_g$, as the co-moving radii to the surface of recombination and PT are approximately the same.
We take $\Lambda_{\phi} \gg H$ such that the contribution in Eq.~(\ref{eq:cmbSIBdec}) is negligible.
Eq.~(\ref{eq:cmbSIBgrav}) along with $\dot{\phi}_{0} \simeq (60H)^2$ \cite{Akrami:2018odb} gives a correction $\left(\Delta C^{\gamma}_{\ell_g} / C^{(0) \gamma}_{\ell_g}\right)_{\Delta H_{\rm osc}} \simeq 6 \times 10^{-5}$, and
Eq.~(\ref{eq:cmbSIBind}) gives $\left(\Delta C_{\ell_d}^{\gamma} / C_{\ell_d}^{(0)\gamma}\right)_{\rm ind} \simeq 0.01$.
This is below the cosmic variance as can be checked with Eq.~(\ref{eq:alphaObs}) (even after accounting for a smaller cosmic variance at $\ell_{g} \approx 200$).
This benchmark is actually conservative in that the current constraints on the PC features from Planck data are even weaker \cite{Braglia:2021rej}, only constraining the amplitudes of PC features to be $\lesssim 0.08$.
Access to higher multipoles $\ell \sim {\cal O} (1000)$ in the CMB is not really an advantage as the $\ell^{-3/2}$ (or $\ell^{-3}$) dependence in the PC features makes them fairly localised near $\ell_d$ (or $\ell_g$).

We now estimate the experimental sensitivity required to probe this PC in SGWB.
The strength of SGWB depends on the details of the PT.
Here, we consider a relatively strong PT
(see \cite{Caprini:2015zlo} for models that could undergo strong PTs).
From Eqs.~(\ref{eq:gwMonopole}) and (\ref{eq:gwSIB}), we get the isotropic component, $\Omega_{\rm GW} h^2 \simeq 5  \times 10^{-10} $, and the small PC variation in anisotropies, $\delta \Omega_{\rm GW}^{\rm PC} h^2 (\ell_d) \simeq 4 \times  10^{-15}$.
Note that the isocurvature and $N_{\rm eff}$ bounds mentioned earlier are easily satisfied.
While the isotropic SGWB could be within LISA sensitivity, $\Omega_{\rm GW}^{\text{\tiny LISA}} h^2 \gtrsim 10^{-14}$ \cite{caprini2016science},
the anisotropies require more advanced future space-based detectors.
For reference, the projected sensitivity of the proposed BBO mission for large angular scales is $\Omega_{\rm GW}^{\text{\tiny BBO}} h^2 \gtrsim 10^{-17}$ \cite{Corbin:2005ny, Schmitz:2020syl} and that of ultimate-DECIGO is $\Omega_{\rm GW}^{\text{\tiny DECI}} h^2 \gtrsim 10^{-19}$ \cite{Kudoh:2005as}.
Detection of the PC will require a detector that maintains such a high sensitivity out to $\ell \sim O(100)$.

Now consider a second benchmark with smaller $f_{\rm VS}$,
\begin{align}
       (2) \quad f_{\rm VS} = 0.01, \, \frac{A_{\chi}}{A_{\phi}}  = 3 , \, \frac{\beta}{H_{\rm PT}} = 5 \,.
\end{align}
The PC corrections in this case from eqs.~(\ref{eq:cmbSIBgrav}) and (\ref{eq:cmbSIBind}) are respectively $\left(\Delta {\cal P}_{k_g} / {\cal P}_{k_g}^{(0)}\right)_{\Delta H_{\rm osc}} \simeq 8 \times 10^{-5}$ and
$\left(\Delta {\cal P}_{k_d} / {\cal P}_{k_d}^{(0)}\right)_{\rm ind} \simeq 10^{-4}$.
3D maps of adiabatic perturbations like LSS and 21-cm are expected to have lower cosmic variance as they can access more modes than the 2D CMB.
Cosmic variance for an idealised 3D map is $\left(\delta {\cal P}/{\cal P}^{(0)}\right) \sim \sqrt{(2 k_{\rm min}^{3})/(4 \pi k^{2} \Delta k)}$, where $k_{\rm min}$ in the largest co-moving mode that can be probes today.
Then, using the largest possible bin size, $\Delta k \sim \pi k H/M$, the cosmic variance at $k_d \approx 100 k_{\rm min}$ is $\left(\delta {\cal P}/{\cal P}^{(0)}\right) \sim 7\times 10^{-4}$, and at $k_g \approx 200 k_{\rm min}$ is $\left(\delta {\cal P}/{\cal P}^{(0)}\right) \sim 2\times 10^{-4}$.
We see that the PC correction in the second case can not be detected even with an ideal 3D map.

The strength of SGWB in this benchmark is smaller, with monopole, $\Omega_{\rm GW} h^2 \simeq 5  \times 10^{-12} $, and the PC feature, $\delta \Omega_{\rm GW}^{\rm PC} h^2 (\ell_d) \simeq 4 \times  10^{-17}$.
While this case is more challenging experimentally, the impossibility of detection even in a perfect adiabatic map makes SGWB a powerful (and possibly the only) probe of such a PC.

The final check is the NG contribution to adiabatic perturbations as shown in \Eq{eq:cmbSIB3pnt}. In both benchmarks, the imparted NG is $f_{\rm NL}^{\rm PC} \sim 5 \times 10^{-3}$, which is well below the CMB constraint $f_{\rm NL} < {\cal O}(5-50)$ (the exact constraint depends on the shape) \cite{Akrami:2019izv}.
A conservative estimate of the cosmic variance for a 3D map gives $\delta f_{\rm NL}^{\rm PC} \sim 1/(\sqrt{N_{k_g}} \zeta_{\phi}) \sim 10^{-2}$, where $N_{k_g}$ is the number of 3-point configurations satisfying $k_1+k_2+k_3 = 2 k_g$ and $\vec{k_1}+\vec{k_2}+\vec{k_3}=0$, and the anisotropy $\zeta_{\phi} \sim 10^{-5}$.
Again, the imparted NG in our benchmarks is safely below  cosmic variance.

\begin{figure*}[ht]
    \subfloat[\footnotesize{Unbinned}]{\includegraphics[width=0.45\linewidth]{benchmark_GW_unbin_2.pdf}}
    \hspace{3em}

    \subfloat[\footnotesize{Binned, $\Delta \ell = 6$}]{\includegraphics[width=0.45\linewidth]{benchmark_GW_bin_2.pdf}}

	\caption{Primordial spectrum $C^{\rm GW}_{\ell}$ with the PC feature using benchmark parameters from Eq.~(\ref{eq:benchmark-common}). The yellow line is the expected signal where the onset of PC feature is happening at $\ell_d=100$ and the green dots incorporate a realization of cosmic variance. Figures (a) and (b) correspond to GW map without and with binning ($\Delta \ell = 6$) respectively.
    }
	\label{fig:benchmark_plots}
\end{figure*}

\begin{figure*}[ht]
    \subfloat[\footnotesize{$\ell_d=100$}]{\includegraphics[width=0.45\linewidth]{parameter_plot_alt.pdf}}
   \hspace{3em}
    \subfloat[\footnotesize{$\ell_d=1000$}]{\includegraphics[width=0.45\linewidth]{parameter_plot_alt2.pdf}}
	\caption{
	Analysis of the $\alpha-M$ parameter space showing excluded regions from following considerations: binned cosmic variance (blue), CMB constraints from contribution to $\Delta {\cal P}_{\gamma}$ (yellow) and $f_{\rm NL}$ (green).
	The dashed lines show contours of PC signal strength $\delta \Omega_{\rm GW}^{PC} h^2 (\ell_d)$ as given in Eq.~(\ref{eq:gwSIB})
	such that the parameter regions above these lines can be probed with corresponding experimental sensitivities.
	The analysis is carried out for (a) $\ell_d =100$, as in the benchmark cases, and (b) $\ell_d =1000$, where the cosmic variance uncertainty is smaller.
	See the text for more details.
    }
	\label{fig:parameter_plot}
\end{figure*}

\section{broader analysis of the parameter space}
Let us now broaden  the exploration of the  parameter space of PC in SGWBs beyond the specific benchmarks above.
We continue to fix the start of the PC feature at $\ell_d=100$, as in the benchmark cases, and analyze observability of PC in $\alpha-M$ space, taking both of these as our independent variables.
The requirement of $\alpha <1$ and $M>H$ ensures $(\sigma_0/\Lambda_{\chi}) < 1$, satisfying the EFT perturbativity constraint.

Requiring that the (envelope) of the PC feature is larger than the binned cosmic variance  excludes the blue region shown in figure~(\ref{fig:parameter_plot}).
Requiring Eqs.~(\ref{eq:cmbSIBgrav}) and (\ref{eq:cmbSIB3pnt}) to lie below current CMB sensitivity, $(\Delta {\cal P}/{\cal P}^{(0)})_{\gamma} < 0.01$ and $f_{\rm NL} < 10$, along with
the restriction
$M < \Lambda_{\chi}$ needed for EFT control,
exclude the yellow and green regions respectively in  \Fig{fig:parameter_plot}.
The contribution from HS decay in Eq.~(\ref{eq:cmbSIBind}) can be kept below CMB exclusion with a conservative choice of $f_{\rm VS}^{2} = 3 \times 10^{-3}$ (for $A_{\chi}/A_{\phi} =3$), while $f_{\rm VS}$ can in principle be larger for smaller $\alpha$.
The contours of different signal strengths,  $\Omega^{\rm PC}_{\rm GW} h^2$, at given $\ell_d$ are also shown in \Fig{fig:parameter_plot} with dashed lines.

We have assumed that the decay rate $\Gamma(\sigma \rightarrow \delta \chi \, \delta \chi) < 3H$ such that the dilution of the classical $\sigma$ field is primarily due to the Hubble expansion, $\propto k^{-3/2}$ dependence in the PC signal. For the interaction considered, $\Gamma(\sigma \rightarrow \delta \chi \delta \chi) \approx \frac{M^3}{32 \pi \Lambda_{\chi}^{2}}$ \cite{Kofman:1997yn}.
For the largest mass $M \simeq \Lambda_{\chi}$, we see that $\Gamma < 3H$ for $M\lesssim 300 H$.
Also, the PC oscillation would be unobservable if it is faster than the integer discretization of $\ell$.
Since the oscillation frequency is fastest near $\ell_d$, we require $(M/H) < \pi \ell_d$.
For $\ell_d = 100$, we again get a similar upper limit $M \lesssim 300H$, although CMB constraints become dominant well before reaching this limit as can be seen in \Fig{fig:parameter_plot}~(a).

At least in principle, PC features with smaller amplitudes can be observable if they appear at larger $\ell_d$, as the cosmic variance becomes smaller at higher multipoles. An example is given in \Fig{fig:parameter_plot}~(b) for $\ell_d=1000$. This highlights the gains that can be achieved with higher angular resolution and sensitivity in GW experiments.

\section{Oscillating heavy field from a polynomial potential}
The above analysis considered an unspecified ``sharp feature'' in the potential \cite{Chen:2011zf,Braglia:2021rej} to excite the heavy field, leading to the simplest form for the PC features. Here we construct a  smoother, more familiar polynomial potential for $\sigma, \phi$, inspired by Hybrid Inflation \cite{Linde:1993cn, Wang:2018tbf}.
In this model, the rolling of the inflaton field triggers a ``waterfall'' for the $\sigma$ field. The waterfall however does not end inflation, but leads to a second phase of inflation.
\begin{align}\label{eq:hybridpot}
    V(\phi, \sigma) = V_{\rm inf}(\phi) - \kappa \sigma+\frac{1}{2}(g\phi^2-\mu^2)\sigma^2+ \frac{\lambda}{4} \sigma^4 \, .
\end{align}
We will take $V_{\rm inf}(\phi)$ to dominate the energy density and hence drive inflation.  $\sigma$ has a time-dependent mass as $\phi$ rolls towards 0,
\begin{align} \label{eq:Msigma}
    M^{2}_{\sigma}(t) = g \phi^{2}(t)-\mu^2 + 3 \lambda \<\sigma\>^2 \,.
\end{align}
Initially, $g \phi^2 \gg \mu^2$ and $\< \sigma \> \approx 0$, $V \approx V_{\rm inf}$. As $\phi$ rolls,  $\sigma$ becomes tachyonic as $g \phi^2 < \mu^2$. The transition occurs at $t_0$ where  $\phi(t_0) = \phi_c$ and
\begin{align}
    g \phi_{c}^{2} = \mu^2 \,.
\end{align}
From here, $\sigma$ starts rolling rapidly towards the new minimum in about $\Delta t \sim H^{-1}$ when the kinetic energy is large enough to overcome Hubble friction.
When $\sigma$ becomes large enough, the last term of Eq.~(\ref{eq:Msigma}) ensures $\sigma$ is non-tachyonic again and
oscillates around the adiabatic minimum $ M_{\sigma}(t)/\sqrt{2 \lambda}$, which approaches a maximum value $\mu/\sqrt{\lambda}$ at late time.
Here we are taking $\kappa$ to be small such that its role in the $\sigma$ dynamics after $t_0$ is subdominant.
Earlier, it breaks the degeneracy of the two $\sigma$ minima when $\sigma$ becomes tachyonic.
This ensures selection of a unique vacuum ($\<\sigma\> >0$ for $\kappa >0$) in the entire observable universe, avoiding the domain wall problem \cite{Linde:1993cn}\footnote{In \cite{Wang:2018tbf}, this is done by considering $\sigma$ to be the radial direction in the quasi-single field model of inflation. Then the centrifugal force due to the turning trajectory acts as the source of explicit symmetry breaking.}.

Consider the following model parameters:
$\mu = 50H$, $g=1.9\times10^{-8}$, $\phi_c = 3.6 \times 10^{5} H$, $\lambda = 0.1$, and $\kappa= 0.1 H^3$ (with $H/M_{\rm pl} = 10^{-5}$ and $\dot{\phi}_{0} \simeq 3.6 \times 10^{3} H^2$ \cite{Akrami:2018odb}).
Let us first check that the back-reaction of $\sigma$ dynamics on $\phi$ is small.
The maximum energy density in $\sigma$ is subdominant to the inflationary potential energy, $\mu^4/\lambda \sim 10^{7} H^4 \ll V_{\rm inf} \simeq 3H^2 M_{\rm pl}^{2} \sim 10^{10} H^4$.
The kinetic energy of $\phi$ is given by $3H \dot{\phi}_{0} \approx -\partial_{\phi} V(\phi,\sigma) = -\partial_{\phi}V_{\rm inf} - g \<\sigma\>^2 \phi$.
The correction $g \<\sigma\>^2 \phi \lesssim O(10) H^3  \ll \partial_{\phi}V_{\rm inf} \sim 10^4 H^3$. This ensures slow roll of the inflaton even after the tachyonic transition in $\sigma$.

\begin{figure}[t]
   \includegraphics[width=0.85\linewidth]{hybrid_sigma.pdf}
	\caption{The classical dynamics of $\sigma$ field after the tachyonic transition at $t_0$ (blue) and the minimum $ M_{\sigma}(t)/\sqrt{2\lambda}$ (yellow). The full dynamics shows oscillations about the minimum.}
	\label{fig:hybrid_sigma}
\end{figure}
\begin{figure}[t]
	\centering
		\includegraphics[width=0.95\linewidth]{hybrid_deltaPchi_corr.pdf}
		\caption{The PC feature in ${\cal P}_{{\cal S}_{\chi}}$ in the model with polynomial potential. The blue dots correspond to the numerical evaluation of $\Delta {\cal P}_{{\cal S}_{\chi}}$, where the overall shift is the result of $\<\sigma\> \neq 0$. For reference, the yellow dashed line shows a matching signal from a simplified scenario where $\sigma$ is instead excited suddenly, resulting in a PC with constant amplitude $\alpha \approx 0.125$ and mass $M \approx 14.8 H$.}
		\label{fig:hybrid_sig_k}
\end{figure}

In \Fig{fig:hybrid_sigma}, the full $\< \sigma (t)\>$ dynamics is shown in blue along with the expected minimum at $ M_{\sigma}(t)/\sqrt{2\lambda}$ shown in yellow.
We can see that the dynamics is given by an oscillatory component superimposed on a slowly shifting minimum.
We consider the same direct coupling to $\chi$ as before, $\mathcal{L}_{\rm int} \sim (\sigma \, \partial_{\mu} \chi \partial^{\mu} \chi)/\Lambda_{\chi}$. The slowly shifting minimum induces a tilt in the otherwise scale-invariant spectrum, while the oscillatory component adds a PC feature.
We take $\Lambda_{\chi} \sim 300 H$ so as to ensure $(\< \sigma \> / \Lambda_{\chi}) < 0.5 $.

Figure~\ref{fig:hybrid_sig_k} shows the fractional correction to ${\cal P}_{{\cal S}_{\chi}}$ evaluated numerically (blue) and
a matching signal (yellow) of the type considered in the previous sections arising from a more sharply featured potential with $\alpha \approx 0.125$ and $M \approx 14.8 H$.
The overall shift in ${\cal P}_{{\cal S}_{\chi}}$ is the result of $\<\sigma(t)\> \neq 0$.
The largest amplitude of PC feature is around $k_d \sim 100 k_{\rm min}$, corresponding to $\ell_d \sim 100$ as expected.
The PC is within the observable parameter space, as can be seen from \Fig{fig:parameter_plot}.
In this model,
the irreducible PC contributions to the CMB from eqs.~(\ref{eq:cmbSIBgrav}) and (\ref{eq:cmbSIB3pnt}) are $(\Delta {\cal P}/{\cal P}^{(0)})_{\gamma} \simeq 0.001$ and $f_{\rm NL} \simeq 0.1$ respectively, which are below CMB sensitivity.
Comparing to the yellow signal due to sudden excitation of the $\sigma$ field, we see that our waterfall model gives a PC that
dilutes slightly slower and oscillates with higher frequency as $M_{\sigma}$ increases with time.

\section{Conclusion}
We have shown that the power spectrum of SGWB superhorizon anisotropies could reveal striking oscillations (see figure~(\ref{fig:benchmark_plots})), which could only arise from extremely heavy fields during inflation.
However, any such new physics introduces irreducible (gravitational) corrections to the spectrum of adiabatic perturbations.
Remarkably,  it is possible for the SGWB oscillations to be clearly visible while the  CMB and LSS remain insensitive to the irreducible corrections, despite their high precision and statistics.
While we have considered specific PC features, the mechanism can be extended to other localised features and to non-Gaussian correlators.
We illustrated a model in which the oscillations in anisotropy appeared at $\ell \sim {\cal O}(100)$ in order to maintain better theoretical control and minimize cosmic variance. Similar features could appear at $\ell \sim {\cal O}(10)$, but would require stronger couplings and a more rigorous data analysis. However, lower experimental angular resolution would then suffice.  In our smoother model, the primordial clock ``beats'' irregularly, but its non-trivial
frequency $M(k)$ seen in \Fig{fig:hybrid_sig_k}
sheds light on the smooth potential and mechanism exciting $\sigma$.

\textit{Acknowledgements:}
We would like to thank Soubhik Kumar, Yuhsin Tsai, Naren Manjunath, and Peter Shawhan for useful discussions. This work was supported by NSF grant PHY-1914731 and the Maryland Center for Fundamental Physics.


\section{Introduction}
Correlation functions are the backbone of quantum field theory.
In the realm of gauge/string correspondence, correlation functions of $n$~single-trace gauge-invariant operators in the planar limit describe the scattering of $n$~closed strings, whose world-sheet is a punctured sphere. Under special circumstances this sphere simplifies further and factorizes into a product of two disks. (Reminiscent of KLT type relations connecting open and closed strings.) This can happen in at least two ways as depicted in figure \ref{octagonExample}.

One option uses supersymmetry. In $\mathcal{N}=4$ SYM say, we can take single-trace half-BPS operators made out of \emph{many} fields and use $R$-symmetry polarizations in order to form a thick frame with many propagators sequentially connected drawing a closed polygon. The frame is SUSY protected from quantum corrections and thus the correlator splits into decoupled dynamics happening inside and outside the polygon \cite{Frank,FrankBootstrap}.

An alternative way to decouple the inside and outside of the polygon is to make the polygon's edges approach the light-cone. If we send to zero the coupling~$g^2 $ as well as the distances between consecutive cusps~$x_{i,i+1}^2 $ with the so-called \textit{cusp times}
\begin{equation}
\label{DS}
t_k^2  \equiv \,
g^2 \log(x_{k-1,k}^2)\log(x_{k,k+1}^2)\,,
\end{equation}
held fixed, we obtain the decoupling without any requirement on the frame's thickness \cite{stampedes1}.
This is the regime we will focus on. This decoupling is dynamical: in this limit a fast classical particle travels around the null perimeter decoupling the quantum dynamics that are confined inside and outside the polygon \cite{Alday:2010zy,AldayBissi,BGV,BGHV,stampedes1}. The \emph{null polygon} limit (\ref{DS}) is a double-scaling limit re-summing
the maximal power of logarithms at each perturbative order. Since it only depends
on very universal one-loop features of the conformal gauge theory, most of our results (if not all) should apply to any $4d$ conformal gauge theory with a planar limit and at least one adjoint scalar \cite{QCD}.

\begin{figure}[t]
\includegraphics[scale=0.28]{talk_4pt_fact}
\caption{Correlators with a large BPS frame (left) and doubly-scaled correlators with operators approaching the cusps of a null polygon (right) both factorize into products of disk topology correlators (middle).}
\label{octagonExample}
\end{figure}

The key objects of this paper are these null polygons of disk topology. We think of them as a tree level skeleton which are then dressed by quantum corrections at loop order. The skeleton is described by a set of internal \textit{bridges} which are bundles of propagators connecting non-consecutive cusps. For squares and pentagons we have a single choice of such bridges (see figures~\ref{octagonExample} and~\ref{5pt_coord}) while for hexagons and higher there are multiple different topologies (see for instance figure \ref{hexagonBridges}).

Belitsky and Korchemsky found a remarkable result for the null square as  \cite{BK0,BK}  $\texttt{square}_{\ell} = e^{-s^2} \tau_{\ell}(s)$ where $\ell$ is the number of bridges, $s^2=\sum t_i^2$ and $\tau_{\ell}$ is a solution of the Toda lattice equation
\begin{equation}
\label{toda}
(s\partial_s)^2 \log \tau_{\ell}(s) = s^2  {\tau_{\ell+1}(s)\tau_{\ell-1}(s)}/{\tau_{\ell}^2(s)}\,.
\end{equation}
The square with no bridges is simply given by the Sudakov exponential $e^{-s^2}$ so that $\tau_0=1$; the exponential can be thought of as a square null Wilson loop \cite{Alday:2010zy,AldayBissi,BK0}. The square with one bridge $\tau_1$, \textit{the seed}, is a non-trivial function of $s$. Together with $\tau_0=1$ and the Toda equation it gives us recursively all other squares. Amusingly the seed itself can be fixed by the simple requirement
that loop corrections in presence of $\ell$ bridges kick-in at $(\ell+1)$ loops: $\texttt{square}_{\ell} =1 + O(s^{2(\ell+1)})$ uniquely fixes $\tau_1$ to be a modified Bessel function $I_0(2s)$. All other squares then evaluate to simple determinants of Bessel functions.

In this letter we conjecture an extension of this result and of the Toda system to any null polygon with any bridge configuration. We will encounter a novel hierarchy of Toda equations, determinantal solutions and generalized Gross-Witten-Wadia matrix integrals \cite{GWW,Brezin}.

\section{From Stampedes to Data}

We first consider pentagons.
If we take two non-consecutive edges to approach a null limit there is a tool we can use to capture the corresponding leading logarithms: \textit{the stampedes} \cite{stampedes1}. There are five different pairs of edges we can take to be null as depicted in figure \ref{pentagonChannels}.
\begin{figure}[t]
\includegraphics[scale=0.385]{letter_channel}
\caption{Three of the five channels for the stampede computation of $\texttt{pentagon}$. (Reflections around the vertical axis yield the two missing channels.) Thick edges are null and expanded in $\texttt{top}$ and $\texttt{bot}$ states, then acted on by the $SL(2)$ hamiltonian $\mathbf{H}$ which captures the associated leading $\log$.  A red dot in $j$ corresponds to sub-leading cusp time $t^2_j \to 0$. } \label{pentagonChannels}
\end{figure}

Consider the choice $x_{12}$ and $x_{45}$ for concreteness; the other choices are treated similarly. A convenient choice of coordinates for this limit is depicted in figure \ref{5pt_coord}. In terms of these coordinates the relevant \textit{cusp times} simply read
\begin{equation}
\begin{array}{ll}
t_1^2={\color{blue} g^2 \log(\tfrac{1}{\bar z})}{\color{magenta} \log(z)} \, ,
& t_2^2={\color{blue} g^2 \log(\tfrac{1}{\bar z})} {\color{green(pigment)} \log(\tfrac{1}{v})} \,, \\
 t_4^2 ={\color{blue}g^2 \log(\bar{w})} {\color{magenta} \log(\tfrac{1}{1-w})}  \, , & t_5^2= {\color{blue}g^2 \log(\bar{w})}  {\color{magenta} \log(z)} \,.
 \end{array}
  \label{timesZ}
\end{equation}
The colouring is meant to highlight the different origin of these logarithms. The magenta ones arise once we Taylor expand the top and bottom edges in terms of top and bottom states
\begin{equation}
\label{states}
\mathcal{O}_1\mathcal{O}_2 \to \sum_{J=0}^{\infty}\frac{z^J}{J!} \,  | \texttt{top}_J \rangle\,,\,\,
\mathcal{O}_4\mathcal{O}_5 \to \sum_{J=0}^{\infty}\frac{w^{-J}}{J!} \,  \langle \texttt{bot}_J |\,,
\end{equation}
where
\begin{align*}
\begin{aligned}
| \texttt{top}_J \rangle = {  | (\partial_z^J X) \underbrace{X ... X}_{h+k+1}\rangle } \,,\,\,
 \langle \texttt{bot}_J |= { \langle  \underbrace{X... X}_{1+h} \partial_w^J(  \underbrace{X ... X}_{k+1})|}\,.
\end{aligned} \notag
\end{align*}
Here we performed already three major simplifications which we learned from \cite{stampedes1}. First, the factorization of the polygon in the limit \eqref{DS} restricts the analysis only to its interior (or, alternatively, exterior) so that $\texttt{top}$ and $\texttt{bot}$ states are mapped to vectors of an \textit{open} spin chain where each site is one field that belongs to the inside of the polygon. Second, the one-loop dilation operator that generates quantum corrections is replaced by the $SL(2)$ spin chain hamiltonian $\mathbf{H}$. That is, the leading action of the dilation on the null polygon is generated only by the \emph{hopping} of light-cone derivatives and we replaced every scalar by $X$. Finally, the thickness $K$ of the frame is irrelevant in this limit where boundary-dependant terms are sub-leading, hence we set $K=1$ \cite{generalLC}.
\begin{figure}[t]
\includegraphics[scale=0.49]{5pt_fact}
\caption{To parametrize the pentagon we take four points on the plane as indicated in the figure and one point $
x_3$ outside the plane.
The five parameters match the five independent cross-ratios. In the null limit ${\bar z}, 1/z, 1/\bar{w}, 1-{w}, v$ go to zero.
}
\label{5pt_coord}
\end{figure}

The blue terms in (\ref{timesZ}) appear through the action of the Hamiltonians that generate the stampedes
\begin{eqnarray}
\label{stampede}
&&\texttt{stampede}^{h,k}_{J,J'}=  \\
&&\qquad =\langle \texttt{bot}_{J'} |   e^{g^2 \log(\bar{w})  \mathbf{H} } \big(\mathbb{I} \otimes |v\rangle \langle v |  \big) e^{g^2 \log(1/\bar{z})  \mathbf{H}} | \texttt{top}_{J} \rangle\,.  \nonumber
\end{eqnarray}
Dependence on $h,k$ is implicit on the rhs. The stampede is given by bottom and top time evolution together with the overlap with the impurity at the right-most site which is the sole effect of the last operator, the off the plane one. For that we use the simple scalar product
\begin{equation*}
\langle J |v\rangle \langle v | J' \rangle  = \frac{1}{v^{J}}  \times 1 \,,
\end{equation*}
see appendix B for more details. In the end we have
\begin{equation}
\label{F_zwv}
\mathcal{F}_{h,k}\equiv \sum_{J,J'=0} \frac{z^J}{w^{J'}} \texttt{stampede}^{h,k}_{J,J'}\,.
\end{equation}
For example for two equal bridges of length $1$ the pentagon starts at two loops and evaluating the smallest stampedes for the lowest spins $J,J'\leq 1$ yields
\begin{eqnarray}
\mathcal{F}_{1,1} &=& 1+ \frac{z}{v}(-g^6 \log^3 \bar z+\dots)+ \nonumber \\
&&\!\!\! + \frac{z}{w}(g^4 \log^2 \bar{z}/\bar w +4 g^6 \log^3 \bar{z}/\bar w +g^6 \log^3 \bar{z} +\dots) \nonumber \\
&&\!\!\!+\frac{1}{w} (-g^4 \log^2 \bar{w}+ 4 g^6 \log^3\bar w+\dots)+\dots \,.
 \label{expansion}
\end{eqnarray}
Such expansions are very easy to generate to large orders in $z$ and $1/w$ and $1/v$ (i.e. large spins) and to relatively large order in $g^2$ (we have gone to five loops, i.e.~$g^{10}$). They are what we call our \textit{data}.

\begin{figure}[t]
\includegraphics[scale=0.52]{5pt_hop}
\caption{The Taylor expansion of null edges (dashed) at order $J/J'$ around $0/\infty$ creates $\texttt{top}_{J}/\texttt{bot}_{J'}$ states, depicted here for $h=3$, $k=2$. Each chain site is a circle: filled ones are those where derivatives are injected. The states $\texttt{bot}$ and $\texttt{top}$ evolve by the $SL(2)$ action at order $(g^2 \log 1/\bar z)^m$ and $(g^2 \log\bar w)^n$. The resulting states can have derivatives distributed  among all sites. Overlap site-by-site is illustrated by dotted lines.}
\label{fig:np_fact}
\end{figure}

\section{From Data to Symbols to Pentagons}
We now need to re-sum (\ref{expansion}) which is an expansion with small $z, 1/v, 1/w$ while we are after the null limit where these three variables are either finite or very large.
This is easy. We found that at $L$ loops we can express the full functions of these three variables in terms of
iterated Goncharov polylogarithms of trancendentality $L$ which can be nicely coded in terms of symbols \cite{symbols}. The symbol ansatz is generated as a linear combination of \emph{words} $a_1 \otimes a_2 \otimes \cdots \otimes a_L$ of \emph{letters} $a_i$ picked from the set
\begin{equation}
\label{symbolSet}
\{z,w,v,1-z,1-w,1-v,w-z,v-z,v- w\}\,.
\end{equation}
(more in appendix C.)
At two loops, for instance, the expansion (\ref{expansion}) is nothing but the Taylor expansion of
\begin{equation}
\notag
 \log^2 \!\bar w\left[\text{Li}_2({1}/{w}) +\tfrac{1}{2} \log ^2\!\left(1\!-\!{1}/{w}\right)\right]\!-\! \log^2\! ({\bar w}/{\bar z} )\text{Li}_2(z/w)\,,
\end{equation}
whose symbol is
\begin{equation}
\la{symbol11}
\!\!  \log^2 \!\bar w\left[  \frac{1}{1-w} \otimes \frac{1-w}{w}\right] \! + \log^2\! ({\bar w}/{\bar z} )\left[ \frac{z-w}{w} \otimes \frac{z}{w} \right]\,.
\end{equation}
Stampedes count movements of bosons (the derivatives) in a spin chain leading to nice rational numbers with combinatorial interpretations. There are no $\pi$'s or $\zeta$'s in expansions like (\ref{F_zwv}) so the symbol misses no information in this case! We simply write an ansatz as in (\ref{symbol11}) with unknown rational coefficients multiplying each word which we then fix by matching with expansions such as \eqref{F_zwv}.

We can then expand the resulting functions in the null limits of interest. For example
\begin{equation*}
\lim_{z\to \infty} \lim_{w \to 1} \lim_{v\to 0}\mathcal{F}_{h,k} = \texttt{pentagon}_{h,k}(t_1,t_2,0,t_4,t_5)\,.
\end{equation*}
Applied to the example $h=k=1$ we get
\begin{align*}
\begin{aligned}
&\left.\texttt{pentagon}_{1,1}\right|_{t_3=0} =1- \tfrac{1}{4} \left(2\, t_1^2
   t_5^2+t_1^4+t_4^4+t_5^4\right)\\&\qquad\qquad\qquad + \tfrac{1}{36} (4\, t_4^6 - 3\, t_1^2 t_2^4 - t_2^6) +\tfrac{1}{9} (t_1^2 + t_5^2)^3+\dots
\end{aligned}
\end{align*}
Considering different channels each time, we can generate a huge amount of null \texttt{pentagons} with a single cusp time set to zero. This allowed us to fix all \texttt{pentagons} up to a \textit{single} constant at five loops when the first completely mixed term $t_1^2t_2^2t_3^2t_4^2t_5^2$ appears. This term can never be detected if we send any of the five cusp times to zero. The $\texttt{pentagon}_{1,1}$ for example is given in appendix D up to five loops.

\section{A gift from the stampedes}
How can we use all these \texttt{pentagons} computed to high perturbative order
to unveil an underlying structure leading to their exact expressions?
We looked for a set of differential equations for these objects akin to the Toda equation (\ref{toda}) for the square. We made an educated ansatz for such equations and used the stampede data to fix their form. This leads to the main formulae of this letter:
\begin{align}
\begin{aligned}
\label{Todas}
&\left(\mathcal D_1+\mathcal D_2\right)\left(\mathcal D_1+\mathcal D_5\right) \log \mathbb{P}_{h,k}= 4t_1^2 \, \frac{\mathbb{P}_{h+1,k} \mathbb{P}_{h-1,k}}{\mathbb{P}_{h,k}^2}\,,\\
&\left(\mathcal D_3+\mathcal D_2\right)\left(\mathcal D_3+ \mathcal D_4\right) \log \mathbb{P}_{h,k}= 4t_3^2\,  \frac{\mathbb{P}_{h,k+1} \mathbb{P}_{h,k-1}}{\mathbb{P}_{h,k}^2}\,,
\end{aligned}
\end{align}
where $\mathcal{D}_j=t_j \partial/\partial t_j$ and
\begin{equation}
\texttt{pentagon}_{h,k}=e^{-t_1^2-\dots-t_5^2}\times \mathbb{P}_{h,k}(t_j)\,,\,\,\,\, \mathbb{P}_{0,0}=1\,.
\end{equation}

 A second crucial piece of information from the stampede are the reductions of $\texttt{pentagon}_{h,k}$ when only specific cusp times are switched on. In particular, when the times $t_1,t_3$ are set to zero equations \eqref{Todas} imply factorization. Setting most t's to zero degenerates the pentagon into squares so we expect these functions to be given by squares. We find precisely such square reductions as
 \begin{equation}
 \label{pentared}
 \mathbb{P}_{h,k}(0,t_2,0,t_4,t_5)= {\tau}_{h}(t_5) {\tau}_{h+k}(t_2) {\tau}_{k}(t_4)\,,
\end{equation}
where each $\tau$-function has an index given by the total number of bridges emitted by the $i$-th cusp.

Henceforth we take the stampedes gift (\ref{Todas},\ref{pentared}) as our starting point and explain how these equations fix \emph{all} polygons at \emph{any} loop order.
\section{Toda Bootstrap}
Under the assumption that $\mathbb{P}_{h,k}$ has a regular Taylor expansion around $t_j ^2=0$
\begin{equation}
\mathbb{P}_{h,k} =\sum_{n_1\dots n_5} \mathbbm{c}_{h,k}(n_1,\dots,n_5) \times  t_1^{2n_1}\cdots t_5^{2n_5}\,, \la{startBootstrap}
\end{equation}
the first of \eqref{Todas} can be used to express any coefficient $\mathbbm{c}(n_1,\dots,n_5)$ as a function of coefficients $\mathbbm{c}_{h,k},\mathbbm{c}_{h\pm 1,k}$ with lower value of $n_1$, due to the factor $4t_1^2$ in the r.h.s. of the equation. Likewise, the second equation in \eqref{Todas} can be used to express each $\mathbbm{c}_{h,k}$ via coefficients with lower $n_3$. The iterative application of the two equations eventually reduces the problem to the determination of
\begin{equation}
\notag
\mathbbm{c}_{h,k}(0,n_2,0,n_4,n_5)\,,\,\,\, \forall \,h,k \in \mathbb{N}\,.
\end{equation}
But these coefficients are precisely those determined by the reduction (\ref{pentared}) so we are done!

Pertinent for higher polygon generalizations, an important observation
 is that $\mathbb{P}_{h,0}$ only only depends on times which are either end-points of the bridges $(t_2,t_5)$ or sandwiched between two such points $(t_1)$; the two times $t_3$ and $t_4$ do not show up. For the first such pentagon we managed to identify its exact re-summation as $\mathbb{P}_{1,0} = \phi(t_1,t_2,t_5)$ and $\mathbb{P}_{0,1} = \phi(t_3,t_2,t_4)$, where
\begin{equation}
\phi(x,y,z) = \sum^{\infty}_{n=0}\sum_{m}^n \sum_{r}^{n-m} \frac{x^{2(n-m-r)} y^{2 m}z^{2 r}}{(n-m)!(n-r)!m!r!}\,. \label{phi}
\end{equation}
Any solution $\mathbb{P}_{h,k}$ with $h,k\geq 1$ can be cast in the compact form of a nested determinant
\begin{align}
\begin{aligned}
\label{determinant}
&\mathbb{P}_{h,k}\! =\! \frac{\begin{vmatrix}\left(\mathcal D_1\!+\!\mathcal D_2\right)^{i-1}\left(\mathcal D_1\!+\!\mathcal D_5\right)^{j-1} \mathbb{P}_{1,k}
\end{vmatrix}_{i,j\leq h}}{(2t_1)^{h-h^2} \,\mathbb{P}_{0,k}^{h-1}} \,,\\
&\mathbb{P}_{1,k}\! =\! \frac{ \begin{vmatrix} \left(\mathcal D_3\!+\!\mathcal D_2\right)^{i-1}\left(\mathcal D_3\!+\!\mathcal D_4\right)^{j-1} \mathbb{P}_{1,1}
\end{vmatrix}_{i,j\leq k}}{
(2t_3)^{k-k^2} \mathbb{P}_{1,0}^{k-1}   }
\,.
\end{aligned}
\end{align}

\begin{figure}[t]
\includegraphics[scale=0.235]{higher_pt_letter_v1}
\caption{A decagon example with blue bridges, stretching between next-to-neighbouring points, and red bridges. Dashed lines cut the null polygon into factorized contributions, each subject to a set of coupled Toda equations. On the right we depicted the boundary data at $t_j = 0$ which fixes the solution.}
\label{fig:factor_n}
\end{figure}
This representation follows iteratively from \eqref{Todas} and from the fact that $\mathbb{P}_{0,k}$ and $\mathbb{P}_{h,0}$ are functions of three cusp times only. Formulae \eqref{determinant} make manifest that everything can be reduced to \emph{seed} functions $\mathbb{P}_{0,1}$ and $\mathbb{P}_{1,1}$. The former is given by (\ref{phi}) while the latter can be bootstrapped by the self-consistency argument presented above; we present it in appendix D up to five loops.
\section{Higher polygons}
Based on the experience collected on \texttt{pentagon} functions, we can readily conjecture the form of any null polygon. We express a polygon with $n$ edges as the product of an exponential of cusp times -- reminiscent of the Sudakov term in a null polygonal Wilson loop -- and a function $\mathbb H$ that depends on the configuration of bridges
\begin{equation}
\texttt{polygon} = e^{-\sum_i t_i^2} \times \mathbb{H}(t_1,\dots,t_n)\,.
\end{equation}
When there are no bridges, as for instance for correlators of short operators, we known from \cite{BGV} that $\mathbb{H}$ is one.

Let $\ell$ be the number of bridges stretching between $t_{i-1},t_{i+1}$, we denote $\hat t_i$ the cusp time comprised between the end-points. We focus on this bridge and write $\mathbb{H}_{\ell}$ leaving other bridge indices implicit. We can state that:
\begin{itemize}
\item there is a Toda equation of type \eqref{Todas} in $t_{i-1},\hat t_i,t_{i+1}$ relating $\mathbb{H}_{\ell}$ to $\mathbb{H}_{\ell \pm 1}$. When $\hat t_i=0$ it becomes a factorization condition $\mathcal{D}_{i-1} \mathcal{D}_{i+1} \mathbb{H}_{\ell}=0$.
\item the function $\mathbb{H}$ depends only on cusp times that emit bridges (e.g. $t_{i-1},t_{i+1}$) or are sandwiched between the two endpoints of a bridge (e.g. $\hat t_{i}$).
\item when $\hat t_{j}=0$, the function $\mathbb{H}$ factorizes into a product of squares $\tau_{n}(t_k)$, one for each cusp $k$ emitting a total of $n$ bridges.
\end{itemize}
This set of statements traces back to $\texttt{pentagon}$ observations. For higher polygons the novel feature is the presence of both blue bridges in figure \ref{fig:factor_n} -- connecting next-to-nearest neighbouring cusps $t_{i\pm1}$ and governed by a $2d$ Toda equation -- and the red bridges which connect further apart cusps. To complete the conjecture we add:
\begin{itemize}
\item The polygon $\mathbb{H}$ is factorized into the product of functions $\mathbb{X}$ of those cusps $t_i$ and $\hat t_i$ comprised between the end-points of a continuous path of blue bridges. (In figure \ref{fig:factor_n}, for instance, we have three such blue paths of lengths $2$, $1$ and $0$ and thus the correlator would factor into three factors.) Paths of length $0$ correspond to $\texttt{square}_n(t_j)$ where n is the total number of bridges emitted by the cusp $t_j$.
\end{itemize}
Let us apply these conjectures to the general \texttt{hexagon}. A planar hexagon can have three different bridge configurations as depicted in figure \ref{hexagonBridges}. Cases $(a)$ and $(b)$ in the figure feature two blue bridges, connected and disconnected respectively. Hence they are subject to two $2d$ Toda equations as does \texttt{pentagon} to which they are closely related.
Case $(a)$ is factorized as
\begin{equation}
\mathbb{H}^{(a)}_{h,l,k} =  \mathbb{X}_{h,l,k}(t_1,\hat t_2,t_3, \hat t_4,t_5) \tau_{l}(t_6)\,.\notag
\end{equation}
The function $\mathbb{X}$ solves two coupled equations in $(t_1,\hat t_2,t_3)$ and $(t_3,\hat t_4,t_5)$ moving $h$ and $k$. (Notice that for $l=0$ the same holds for $\texttt{pentagon}$, which means $\mathbb{X}_{h,0,k}=\mathbb{P}_{h,k}$.) Case $(b)$ features two disconnected blue bridges, i.e. it solves decoupled equations and is thus
factorized into two functions of three cusps and expressed via $\texttt{pentagon}$ functions already determined in this letter,
\begin{equation}
\mathbb{H}^{(b)}_{h,l,k}=\mathbb{P}_{l,h}(0,t_6,\hat t_1,t_2,0) \mathbb{P}_{l,k}(0,t_3,\hat t_4,t_5,0)\,.\notag
\end{equation}
Case $(c)$ features three connected blue bridges
and thus is solution to three pair-wise coupled $2d$ Toda equations in $(t_1,\hat t_2,t_3)$, $(t_3,\hat t_4,t_5)$ and $(t_5,\hat t_6,t_1)$ moving $h$, $k$ and $l$ respectively. The boundary data for the three cases are
\begin{eqnarray*}
\mathbb{H}^{(a)}_{h,k,l}(t_1,0,t_3,0,t_5,t_6)&=& \tau_{h}(t_1) \tau_{h+k+l}(t_3) \tau_{k}(t_5) \tau_l(t_6) \,, \nonumber\\
\mathbb{H}^{(b)}_{h,k,l}(0,t_2,t_3,0,t_5,t_6)&=&\tau_{h}(t_2) \tau_{h+l}(t_6) \tau_{k+l}(t_3)  \tau_{k}(t_5)\,,\nonumber \\
\mathbb{H}^{(c)}_{h,k,l}(t_1,0,t_3,0,t_5,0)&=& \tau_{h+l}(t_1) \tau_{h+k}(t_3) \tau_{k+l}(t_5)\,.
\end{eqnarray*}
We tested these various statements against stampede data at three loops for a few six-point correlators \cite{hexa_hop}.

\begin{figure}[t]
\includegraphics[scale=0.36]{hexa_topo}
\caption{The three possible configurations of bridges for the $\texttt{hexagon}$ in the planar limit. Blue bridges are associated to $2d$ Toda equations of type \eqref{Todas}.}
\label{hexagonBridges}
\end{figure}

\section{Discussion}

In this paper we used the stampede technology developed in \cite{stampedes1} to unveil new integrable structures in the null limits of conformal gauge theories. The final picture is quite compact. In the double-scaling limit~\eqref{DS}, correlators are characterized by internal bridges which can be blue (if they connect next-to-neighbouring cusps) or red (if they connect farther separated vertices). With zero blue bridges the correlators factorize into products of square functions (the boundary data); finite blue lengths can then be constructed recursively by means of Toda equations \eqref{Todas} \cite{GrishaHint}.

Some of the solutions to these PDEs are rather novel. If we take a five point correlator with two non-zero cusp times, for instance, we find
\begin{align*}
\begin{aligned}
&\mathbb{P}_{h,k}(t_1,t_2,0,0,0) = \det_{1\leq i,j\leq {h+k}} \begin{vmatrix} z_i^{i-j} I_{i-j}\left(2 z_i \right)\end{vmatrix}\,,
\end{aligned}
\end{align*}
where $I_j$ are modified Bessel functions of the first kind and $z_i=\sqrt{t_1^2+t_2^2}$ if $i\leq h$ and $z_i=t_2$ for $i>h$. It admits a
representation as an exotic matrix integral
\begin{equation}
\notag
 \int \! \frac{d \theta_1}{2\pi}\cdots\!\int\! \frac{d\theta_{\ell}}{2\pi} \Delta(Z) e^{\text{Tr}(Z+Z^{\dagger})} \prod_{k} (z_k e^{i \theta_k})^{1-k},
 \end{equation}
where $Z$ is a matrix with eigenvalues $\{z_j e^{i\theta_j}\}$ and $\Delta$ is the Vandermonde determinant. (Details in appendix E.) When $x_i=y$ this integral reduces to the Gross-Witten-Wadia matrix model \cite{GWW}. What is the physical origin of these matrix integrals?

It would be fascinating to connect the structures found here to other exact approaches to quantum gauge theories. We can envisage three natural connections.

A first connection would be with conformal bootstrap explorations. In \cite{AldayBissi,BGV,BGHV,Antunes} the conformal bootstrap was applied to the study of null correlation functions. There, no simplifying double-scaling limit was taken. On the other hand the analysis was restricted to correlation functions of the lightest fields of the theory which translates into polygons with no internal bridges.
Indeed, if we take the Sudakov and recoil factors bootstrapped there, in the limit (\ref{DS}) we would obtain the universal exponentials
\begin{equation}
e^{-t_1^2-\dots - t_n^2}\,
\end{equation}
 showing up everywhere in our letter. When the bridges are non-zero, dressing these exponentials we have the very nontrivial solutions to the Toda hierarchy. Would be remarkable to fit them inside a bootstrap analysis. Combining the technology developed here with a bootstrap approach might lead to a more general null limit for any polygon. The physical picture should be one where in the null limit correlators become Wilson loops dressed by further cusp insertions of adjoint fields. The latter then propagate from cusp to cusp interacting with the planar chromodynamic flux tube of the theory.

 We computed the hexagon for any bridge topology as depicted in figure \ref{hexagonBridges}; from an OPE perspective the last one resembles a snowflake while the other two mimic the structure of comb channel decompositions. Despite some encouraging first steps \cite{Goncalves:2019znr, Rosenhaus:2018zqn}, there is no analytic result for comb channels for $n > 4$ points. Can one use our hexagon results to improve this state of affairs?

A second connection would be with integrability. In planar $\mathcal{N} = 4$ SYM four-point disk correlators were computed for any kinematics in \cite{Frank,FrankBootstrap} by means of bootstrap plus hexagonalization \cite{hexagons,Hexagonalization,morehexagons2,Determinant} and preliminary higher polygons explorations were started in \cite{decagons,Hexagonalization2}. Would be very interesting to  connect our results with these. In the hexagonalization formalism, for instance, the limit~\eqref{DS} is governed by an interesting region where all rapidities and bound-state numbers of the mirror particles are (analytically continued to) zero. The square was identified with an amplitude in the Coulomb branch \cite{FrankSimon}; are the doubly-scaled higher polygons found here higher-point scattering amplitudes in some interesting kinematics?

A last connection could be with supersymmetric localization. Determinants of Bessel functions and Toda equations  constantly pop up in such supersymmetric studies, see e.g. \cite{GiombiShota,KolyaBremss,JaumeZohar} --  is this a coincidence?

Finally, can we fix the so called pentagon at all loops following a bootstrap \`{a}-la \cite{Frank,stampedes1} for the square (called octagon there)? A strategy would be to construct a symbol ansatz for five point functions of large operators. (We already know some of its letters from (\ref{symbolSet}) arising from the stampede analysis here.) Full null polygon limits will not be enough to fix the ansatz but single light-cone limits – which are still fully captured by the stampedes – might suffice once combined with simple additional physical conditions such as Steinmann relations and single-valuedness in the physical sheet.

\begin{acknowledgments}
\section*{Acknowledgments}
\label{sec:acknowledgments}

We thank B.~Basso, A.V.~Belitski,  C.~Bercini, V.~Goncalves, A.~Homrich, S.~Komatsu and especially G.P.~Korchemsky for discussions and valuable comments on our draft.
Research at the Perimeter Institute is supported in part by the Government of Canada
through NSERC and by the Province of Ontario through MRI. This work was additionally
supported by a grant from the Simons Foundation (Simons Collaboration on the Nonperturbative Bootstrap \#488661) and ICTP-SAIFR FAPESP grant 2016/01343-7 and FAPESP grant 2017/03303-1. \end{acknowledgments}


\section{Introduction}

In recent years papers appeared which discuss magic and mana properties for different models, such as CFT, knot theory and others \cite{ConfMagic,KnotMagic,Chaos}. These quantities characterize how far is the certain quantum mechanical state defined for such a model from the element of the Clifford group \cite{Cliff}. According to the Gottesmann-Knill theorem \cite{GKnill}, Clifford group elements can be effectively modeled on a classical computer. Thus it is claimed that ``magic'' is in effect a non-classicality of a certain state, and mana measures this non-classicality. These properties can be important if one discusses these properties in relation to the quantum computations.

The Gottesman-Knill theorem is based on a fact that Clifford group is a finite subgroup of studied group $G$, which is a tensor product of several $SU(N)$'s. However, it is not the only finite subgroup. One can define infinitely many such subgroups for the same group $G$.  Among these the defining property of the Clifford group is its connection to the sigma matrices. From the point of view of quantum computing there is no need to demand this. Thus depending on the set of problems one wants to present to quantum computer, mana can be defined differently. Our claim is that mana is in fact a relative rather than absolute property.

In the present paper we will present how the Clifford group is usually defined and how it can be modified to get other finite subgroups.
We will apply this new mana definition to studying knot states.

Knot theory is a widely studied subject with lots of relations to other theories. Among others there are connections between knot theory and quantum computations which provide both approaches to calculate knot polynomials using quantum algorithms, as well as describing quantum algorithms as some knot configuration in effective topological field theory \cite{tqc}-\cite{tqcKol2}. This involves calculating knots using unitary matrices through the Reshetikhin-Turaev algorithms \cite{RT}-\cite{RTfin}. Specifically for some particular series of knots any quantum algorithm can be described as consecutive approximations by a series of knots \cite{tqcKol,tqcKol2}.

However, in the present paper we discuss different approach to the knot theory. Mana and magic are properties of the quantum states (density matrices) rather than unitary operations. There is a way to define quantum states, corresponding to a knot \cite{KnotMagic}, using ideas of topological field theories \cite{Ati,Witt}. Matrix elements of this density matrix are made from the knot polynomials in special points. Thus classicality of such states provides us with some information on how can these knot invariants be calculable on classical computer.

Paper is organised as follows. In Chapter 2 we define Clifford group, which is a finite subgroup of the $SU(N)$ group. In Chapter 3 we provide a definition of mana as it is given in other papers on the subject, such as \cite{ConfMagic,KnotMagic,Chaos}. In Chapter 4 we discuss ambiguity in mana definition and show how can the definition be modified to give mana connected to a differen finite subgroup of $SU(N)$. In Chapter 4 we define quantum mechanical states, describing different knots, according to \cite{Ati,Witt,KnotMagic}. In Chapter 5 we study how mana looks like for the knot states and how it can be changed by defining Mana differently.

\section{Clifford group}

Clifford group was first defined by D. Gottesman \cite{Cliff}. Let us take a system of $d$ orthormal states $|k>$, $k=0\ldots d-1$. We take a pair of operators $z$ and $x$:

\begin{equation}
\begin{array}{l}
Z=\sum\limits_{k=0}^{d-1}\omega^k|k><k|,\ \ \ \ \omega=e^{\frac{2\pi i}{d}},
\\
\\
X=\sum\limits_{k=0}^{d-1}|(k+1) \text{mod}\ d><k|.
\end{array}
\label{eq:ZX}
\end{equation}
Using these operators one could define generalized Pauli operators
\begin{equation}
T_{aa^{\prime}}=\left\{\begin{array}{l}
i^{aa^{\prime}}Z^a X^{a^{\prime}},\ \ d=2,
\\
\\
\omega^{-\bar{2}aa^{\prime}}Z^a X^{a^{\prime}},\ \ d>2,
\end{array}\right.
\end{equation}
where $\bar{2}$ is a multiplicative inverse of $2$: $2\times\bar{2}\equiv 1 \text{mod}\ d$. Form the generalized Pauli operators one can defined strings of such operators -- Pauli strings:
\begin{equation}
\label{eq:ps}
T_{\mathbf{a}}=T_{a_1a^{\prime}_1}T_{a_2a^{\prime}_2}\ldots T_{a_na^{\prime}_n}.
\end{equation}

The Clifford group is defined as a set of unitary operators $U$ which transform Pauli string into another Pauli string up to a phase:
\begin{equation}
\label{eq:ClDef}
\mathcal{C}=\left\{U\ \ :\ \ UT_{\vec{a}} U^{\dagger}=e^{i\phi} T_{\vec{b}}\right\}.
\end{equation}

Clifford gates -- elements of the Clifford group -- act in the space of the size $d^n$. There are three Clifford gates which generate the whole group. Two of them act in one $d$-dimensional space, namely phase gate and Hadamard gate:
\begin{equation}
\begin{array}{l}
K=\left(\begin{array}{llllll}
1 \\
& \omega \\
&& \omega^2 \\
&&\ldots \\
&&& \omega^{\frac{(d-1)(d-3)}{2}} \\
&&&& \omega^{\frac{(d-1)(d-2)}{2}}
\end{array}\right),
\\ \\
H=\cfrac{1}{\sqrt{d}}\left(\begin{array}{llll}
1 & 1 & \ldots & 1 \\
1 & \omega & \ldots & \omega^{d-1} \\
\ldots & \ldots & \ldots & \ldots \\
1 & \omega^{d-1} & \ldots & \omega^{(d-1)^2}
\end{array}\right).
\end{array}
\end{equation}

Also there is one operator which acts on a pair of $d$-dimensional spaces:
\begin{equation}
\label{eq:Sop}
S=\sum\limits_{ij}|i;i\oplus j><i; j|.
\end{equation}
In the $SU(2)$ case this operator is the CNOT-gate.

Interesting fact is that there is a finite (although quite large) number of elements in the Clifford group.  Due to this fact anything constructed from the Clifford gates can be effectively simulated on a classical computer in polynomial time. Interesting fact is that Clifford group includes the entangling S-gate (CNOT). Thus non-classicality (magic) is in fact a different parameter from the entanglement.

\section{Mana}

To measure the degree of magic of a certain state a quantity called mana was introduced. It is defined through a phase space point operator $A_{\vec{a}}$:
\begin{equation}
\label{eq:PhP}
A_{\vec{a}}=d^{-n}T_{\vec{a}}\sum\limits_{\vec{b}} T_{\vec{b}}T^{\dagger}_{\vec{a}},
\end{equation}
where $T_{\vec{a}}$ are Pauli strings (\ref{eq:ps}). using these operators one can define discrete Wigner function $W_{\rho}(\vec{a})$ of a state with density matrix $\rho$:
\begin{equation}
W_{\rho}(\vec{a})=\frac{1}{d^n} \text{Tr}\ \rho A_{\vec{a}}.
\end{equation}
Phase space point operators form a complete orthonormal basis in the space of $d^n\times d^n$ matrices, thus density matrix can be constructed from the Wigner functions:
\begin{equation}
\rho=\sum\limits_{\vec{a}}W_{\rho}(\vec{a})A_{\vec{a}}.
\end{equation}
This means that for any physical state when $\text{Tr}\ \rho=1$, set of Wigner functions satisfy $\sum W_{\rho}(\vec{a}) = 1$. Mana is defined as a logarithm of negativity of the set of Wigner functions:
\begin{equation}
M(\rho)=\log \sum\limits_{\vec{a}} |W_{\rho}(\vec{a})|.
\end{equation}

Mana possess several interesting properties which describes why this definition was chosen. First, it can be shown that mana is equal to zero if and only if the density matrix $\rho$ is made from Clifford gates. Second, mana is additive:
\begin{equation}
M(\rho_a\otimes\rho_b)=M(\rho_a)+M(\rho_b).
\end{equation}
Third, mana is related to the second Renyi entropy $S_2$:
\begin{equation}
M(\rho)\leq\frac{1}{2}(L \log d -S_2).
\end{equation}

\section{Ambiguity}

The main property of the Clifford group, related to the ease of classical computation is its finiteness. However Clifford group is not the only finite subgroup of the $U(d)^{\otimes n}$ group. From the Clifford group one can easily construct other finite subgroups using just the rotation matrices.

Let us take as an example the case of one unitary group. Then instead of $X$ and $Z$ operators one can take
\begin{equation}
\tilde{X}=WXW^{\dagger},\ \ \  \tilde{Z}=WZW^{\dagger},
\end{equation}
where $W$ is a matrix from the $SU(d)$ group.
These provide another set of matrices instead of the generalized Pauli matrices:
\begin{equation}
\tilde{T}=W T W^{\dagger}.
\end{equation}
Instead of Pauli strings one will also get generalized strings which can be rotated independently for each of the generalized Pauli matrices in the string. These generalized Pauli strings will also be related by some finite group:
\begin{equation}
\tilde{\mathcal{C}}=\left\{\tilde{U}\ \ :\ \ \tilde{U}\tilde{T}_{\vec{a}} \tilde{U}^{\dagger}=e^{i\phi} \tilde{T}_{\vec{b}}\right\}.
\end{equation}
These matrices are related to the Clifford group matrices when there is only one $U(d)$ group by rotation with $W$. Thus instead of Hadamard and Phase gate this group will be generated by
\begin{equation}
\tilde{K}=WKW^{\dagger},\ \ \ \tilde{H}=W H W^{\dagger}.
\end{equation}
When tensor product of several $U(d)$ groups is considered, each group can be rotated independently of others. The $S$ operator from (\ref{eq:Sop}) should be modified accordingly:
\begin{equation}
\tilde{S}=\mathcal{W}_2S\mathcal{W}_2^{\dagger},
\end{equation}
where $\mathcal{W}_n=W_1\otimes W_2\otimes\ldots\otimes W_n$, with matrices $W_1$ and $W_2$ acting on the corresponding pair of two dimensional spaces.

By modifying the Clifford group, the definition of Mana also should be modified. Namely phase space point operator instead of (\ref{eq:PhP}) will be equal to
\begin{equation}
\tilde{A}_{\vec{a}}=\mathcal{W}A_{\vec{a}}\mathcal{W}^{\dagger},
\end{equation}
which also modifies the definition of Wigner function:
\begin{equation}
\begin{array}{r}
\tilde{W}_{\rho}(\vec{a})=\frac{1}{d^n}\text{Tr}\ \rho \tilde{A}_{\vec{a}}=\frac{1}{d^n}\text{Tr}\ \rho \mathcal{W}A_{\vec{a}}\mathcal{W}^{\dagger}=
\\ \\
=\frac{1}{d^n}\text{Tr}\ \mathcal{W}^{\dagger}\rho\mathcal{W} A_{\vec{a}}=W_{\mathcal{W}^{\dagger}\rho\mathcal{W}}(\vec{a}).
\end{array}
\end{equation}

\section{Knot states}

Let us discuss applications of the generalized mana definition to the knot theory. To define mana for knots, one should first define some quantum states, related to knots. This can be done as follows, based on the papers by Atiyah \cite{Ati} and Witten \cite{Witt}. From the physics perspective knot theory is a three-dimensional Chern-Simons theory with action
\begin{equation}
S=\frac{k}{4\pi}\int d^3 x \left(\mathcal{A}\wedge d\mathcal{A}+\frac{2}{3}\mathcal{A}\wedge\mathcal{A}\wedge\mathcal{A}\right),
\end{equation}
where $k$ is an integer, called level of the Chern-Simons theory. Knot polynomials are equal to the Wilson-loop averages of the Chern-Simons theory:
\begin{equation}
J^{\mathcal{K}}_r(q)=\left<\text{Tr}_r\oint\limits_{\mathcal{K}} d\vec{x}\vec{\mathcal{A}}\right>.
\end{equation}
Polynomial depends on the gauge group of the Chern-Simons theory, its level $k$, representation $r$ of the gauge group and contour of integration $\mathcal{K}$. In what follows for the sake of simplicity we will speak only of the $SU(2)$ group, when knot polynomials are Jones polynomials. Contour of integration is a closed curve in three-dimensional space or a knot. Finally the answer for the Wilson-loop average happens to be a Laurent polynomial in variable $q$, related to the level of the Chern-Simons theory:
\begin{equation}
q=e^{\cfrac{\pi i}{k+2}}.
\end{equation}
Jones polynomials are in fact related to the representations of the quantum group $U_q(sl(2))$, rather than $SU(2)$ group. When $k$ is an integer then $q$ is a root of integer and then there is a finite number of highest-weight representations of the $U_q(sl(2))$ quantum group. Namely for even $k$ there are only representations, corresponding to the spin $m/2$ for $m\leq k$ \cite{QG}. This allows us to define a basis of representations for the knot states. Namely, state function can be defined as follows:
\begin{equation}
|\mathcal{K}>=\sum\limits_{j=0}^{k-1} J_j |j>.
\end{equation}
Using such a wave function one can also define a density matrix:
\begin{equation}
\rho_{\mathcal{K}}=|\mathcal{K}><\mathcal{K}|,
\label{Pure}
\end{equation}
and study mana for such states. In the present paper we will discuss only two-strand knots, which polynomials can be calculated using the following formula \cite{RT6}:
\begin{equation}
J^{T[2,n]}_j=\sum\limits_{i=0}^j \frac{q^{2j-2i+1}-q^{-2j-2i-1}}{q^{j+1}-q^{-j-1}} \left((-1)^i q^{i^2-2ji-i+j}\right)^{n},
\label{Jones}
\end{equation}
where $n$ is a number of twists in a two-strand braid. For odd $n$ the closure of braid produces a knot, while for even $n$ the result is a link.

\section{Mana for the Knot states}

Mana for two strand knots is periodical in parameter $n$ with period being equal to $k+2$, but besides that mana for many knots coincide. On the Fig.\ref{f:2Mana} mana for $k=2$ is displayed. Knots appear only if $n$ are odd integers in (\ref{Jones}). However, using formula (\ref{Jones}), one can extend this definition to the other values of $n$, which gives graph on the Fig. \ref{f:2nMana}. This is purely an analytical continuation of the results for knots to an arbitrary $n$ by using (\ref{Jones}).

\begin{figure}[h!]
\centering
\begin{center}
\includegraphics[width=0.45\textwidth]{k=2p.eps}
\end{center}
\caption{1. Mana for two-strand knots for $k=2$ in Clifford group basis \label{f:2Mana}}
\end{figure}

\begin{figure}[h!]
\centering
\begin{center}
\includegraphics[width=0.45\textwidth]{k=2.eps}
\end{center}
\caption{2. Mana for two-strand knots for $k=2$ with continuation to arbitrary real $n$ in Clifford group basis.\label{f:2nMana}}
\end{figure}

Similarly mana for $k=3$ (see Fig.\ref{f:3Mana}) and $k=4$ (see Fig.\ref{f:4Mana}) can be produced.

\begin{figure}[h!]
\centering
\begin{center}
\includegraphics[width=0.45\textwidth]{k=3.eps}
\end{center}
\caption{3. Mana for two-strand knots for $k=3$  in Clifford group basis. \label{f:3Mana}}
\end{figure}

\begin{figure}[h!]
\centering
\begin{center}
\includegraphics[width=0.45\textwidth]{k=4.eps}
\end{center}
\caption{4. Mana for two-strand knots for $k=4$ in Clifford group basis.\label{f:4Mana}}
\end{figure}

For $k=2$, as can be seen from Fig.\ref{f:2Mana}, there are two distinct values of mana (and in fact only two different density matrices). By changing the basis we calculate mana in, one can make for the other half of 2-strand knots to be equal to zero, while making the mana for the other set to be nonzero. This in fact can be done not only for $k=2$, but also for other cases. This means that by changing the basis one can make different sets of knots states to be effectively calculable on a classical computer.

Specific property of knot states is that they are pure, i.e. $\rho$ is a tensor product of ket and bra vectors (\ref{Pure}). The state for unknot ($n=1$ in (\ref{Jones})) is defined by a vector $v=\frac{1}{\sqrt{d}}[1,1,..,1]$ and $\rho=v^{\dagger}\times v$ for any $k$ and always has mana equals to zero. By transitivity of unitary group it is always possible to rotate any other knot state to the same vector $v=U v_1$. This immediately gives rotated density matrix or, in other words, new basis for Clifford group, $\rho=U v_1^{\dagger} \times v_1U^\dagger=v^{\dagger} \times v$ with zero mana.

In $k=2$ case the only distinguished state besides unknot is trefoil knot ($n=3$ in (\ref{Jones})) with state vector $v_1=\frac{1}{\sqrt{3}}[1,-1,1]$. The easies way to rotate it to the unknot state vector $v$ is by using unitary matrix
\begin{equation}
\setlength{\arraycolsep}{0pt}
S=\left(\begin{array}{ccc}
1\phantom{0} & 0\phantom{0}\phantom{0} & 0
\\
0\phantom{0} & 0\phantom{0}\phantom{0} & 1
\\
0\phantom{0} & -1\phantom{0} & 0
\end{array}\right)
\setlength{\arraycolsep}{5pt}
\label{S2}
\end{equation}.
Of course this matrix is not unique due to stability subgroup SU(2) and U(1) unambiguity in definition of $v$ or $v_1$. But this freedom is not enough to rotate both $v$ and $v_1$ to zero-mana states.

\begin{figure}[h!]
\centering
\begin{center}
\includegraphics[width=0.45\textwidth]{k=2shift.eps}

\end{center}
\caption{5. Mana for two-strand knots for $k=2$ in the basis, rotated by matrix $S$ from (\ref{S2}).\label{f:2ManaR}}
\end{figure}

\begin{figure}[h!]
\centering
\begin{center}
\includegraphics[width=0.45\textwidth]{k=3shift.eps}
\end{center}
\caption{6. Mana for two-strand knots for $k=3$ in the basis, rotated by matrix $S$ from (\ref{S3}).\label{f:2ManaR}}
\end{figure}

The $k=3$ case can be analytically solved in a similar way, the rotation matrix between $v$ and $v_1$ can be found for instance in a following way (rounded answer)
\begin{equation}
\setlength{\arraycolsep}{0pt}
S=\left(\begin{array}{cccc}
0.07+0.5 i\phantom{0} & -0.7-0.5 i & 0 & 0
\\
0.7-0.5 i\phantom{0} & 0.07-0.5 i\phantom{0} & 0 & 0
\\
0 & 0 & 0.07-0.5 i\phantom{0} & 0.7-0.5 i
\\
0 & 0 & -0.7-0.5 i\phantom{0} & 0.07+0.5 i
\end{array}\right)
\label{S3}
\setlength{\arraycolsep}{5pt}
\end{equation}

Thus we illustrated for k=2 and k=3 that in different basis different knot states can have zero mana but only separately. There is no basis in which two or more different knot states have zero mana simultaneously. That in fact shows that even if the absolute value of mana for a state is meaningless (it can always be set to zero), the set of values of mana for different states has an invariant property, they can not be set to zero together, thus are not calculable on a classical computer.

\section{Conclusion}

In this paper we described how the definition of magic and mana can be changed, while keeping its main property - finiteness of classical calculations. This can be done by rotating either the basis of Clifford group, or rotating the density matrix, which is in principle the same.

We applied these  principles to the studies of mana of knot states. Knots have deep connections with quantum calculations through Reshetikhin-Turaev formalism, which relies on unitary matrices and thus becomes a natural task for quantum computer. Another important connection is topological quantum computer which in turn relies on knots as its basic algorithms.

We studied the mana of 2-strand knots and we also managed to find rotation matrices which allows to change which knots have zero mana and thus are classically calculable.

The goal of this paper was to make it clear that, while mana is an interesting property of the quantum state, which should measure its ``classicality'', its definition is in fact very subjective. It relies on the specific finite subgroup of the $SU(N)$ group. However one can define a different finite subgroup which also changes the mana definition. This new mana also in fact measures classicality, but in a different basis.

Clifford group is only distinguished by its relation to the $X$ and $Z$ operators (\ref{eq:ZX}), which has historical significance, but is not always connected to the real quantum systems, used to build quantum computers. In fact for many quantum computers there is a certain freedom to choose different basic operations, quantum gates. For example for topological quantum computer natural choice of operations are $\mathcal{R}$-matrices, which has no relation whatsoever to the $X$ and $Z$ operators. Thus it is important to understand which states are ``classical'' in relation to the exact choice of basic operations.

We showed using the example of knot states that depending on the basis different knot states can have zero mana in different bases. This in fact means that mana is a relative rather than absolute property. Calculations on a quantum computer (or on a classical reversible computer, which is its classical counterpart) are made by using unitary matrices (or permutation matrices in the classical case), which transform one state of the quantum system into another. If these two states can have zero mana simultaneously in some basis, then the corresponding calculations can be effectively made on a classical computer. Thus mana definition should be chosen specifically for the problem we are trying to study.

While in this paper we used knot states purely as an example to show that mana can be defined differently, it is very interesting to study the dependence of mana in different bases on the knot, as well as mana for the states, related to other knot polynomials. This however still remains to be studied.

\section*{Acknowledgements}

We are grateful for very useful discussions to N. Kolganov.

This work was supported by Russian Science Foundation grant No 18-71-10073.