\section{Introduction}

Intermittent locomotion where organisms alternate between active propulsion and a passive phase is frequently encountered in the natural world~\citep{interlocokramer,Gleiss2011}. It is observed in terrestrial, aquatic and aerial modes of locomotion in a diversity of species ranging from unicellular organisms such as ciliates to various reptiles, birds and mammals. Terrestrial organisms that exhibit intermittent locomotion typically come to a complete stop during the passive phase, but organisms in air or water may continue to glide forward~\citep{interlocokramer}. Intermittent locomotion has also been identified in artificial systems. For example, in artificial microswimmers that are powered by chemical activity, the self-interaction of the swimmer with its long-lived chemical wake results in speed bursts~\citep{stopgoswim}. In this paper, we explore the intermittent locomotion that emerges in a system of a self-propelled droplets on a vibrated liquid bath.

Vertically vibrating a bath of silicone oil with frequency $f$ can give rise to steadily walking droplets on the free surface of the liquid~\citep{Couder2005,Couder2005WalkingDroplets}. The walking droplet, also known as a walker, emerges just below the Faraday instability threshold above which the fluid-air interface becomes unstable to standing subharmonic Faraday waves of frequency $f/2$~\citep{Faraday1831a,Molacek2013DropsTheory}. Each bounce of a walker generates a localized, slowly decaying standing wave. The droplet interacts with these self-generated waves on subsequent bounces, resulting in a self-propelled droplet-wave entity on the liquid surface. Such walkers have been shown to mimic several peculiar behaviors that were previously though to be exclusive to the quantum realm~\citep{Bush2015,Bush2020review,Fort17515,harris_bush_2014,Oza2014,Perrard2014b,Perrard2014a,labousse2016,Zeeman,spinstates2018,PhysRevE.88.011001,Giletconfined2016,Saenz2017,Cristea,durey_milewski_wang_2020,Friedal,Eddi2009,tunnelingnachbin,tunneling2020,ValaniHOM,correlationnachbin,Durey2020hqft}. Recently, a new class of walking droplets, coined superwalkers, has been shown to emerge when the bath is vibrated simultaneously with two frequencies, $f$ and $f/2$, along with a relative phase difference $\Delta\phi_0$ between them~\citep{superwalker}. Superwalkers are typically bigger and faster than single-frequency driven walkers and their inter-droplet interactions give rise to novel multi-droplet behaviors.

A given sized walker at a fixed driving amplitude and frequency typically has a fixed walking speed. By contrast, a given sized superwalker at fixed driving amplitudes and frequencies can have a range of walking speeds, dictated by the phase difference $\Delta\phi_0$. Depending on the value of $\Delta\phi_0$, small- to moderate-sized superwalkers ($0.4-0.7$\,mm in radius) can either be in the superwalking or pure bouncing regime while larger superwalkers (bigger than $0.7$\,mm in radius) can also coalesce with the bath for a narrow range of $\Delta\phi_0$ values. \citet{superwalker} observed that a slight detuning of the two driving frequencies to $f$ and $f/2+\epsilon$, results in a \emph{stop-and-go motion} (SGM) for superwalking droplets. In this locomotion, a droplet was observed to periodically switch between active self-propulsion and passive bouncing, see Fig.~\ref{fig: schematic}(a). In this paper, we adopt the theoretical model for superwalkers developed by \citet{superwalkernumerical} and explore the SGM of superwalking droplets through numerical simulations. We start by presenting the theoretical model in Sec.~\ref{theory} followed by a description of the emergence of SGM in Sec.~\ref{emergence}. In Sec.~\ref{PS}, we provide a description of various kinds of SGM observed in the parameter space spanned by the droplet size and the detuning parameter $\epsilon$ and explore them in detail in Sections~\ref{BF}, \ref{FF} and \ref{small detun}. We provide concluding remarks in Sec.~\ref{DC}.

\section{Theoretical formulation}\label{theory}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig_1.pdf}
\caption{Schematic of the droplet system. (a) Snapshots of a superwalker undergoing stop-and-go motion (SGM) showing the passive bouncing (stop) phase and the active self-propulsion (go) phase as observed in experiments. (b) The theoretical setup has a droplet located at a vertical position $z_d$ with an underlying wave field $h(\mathbf{x},t)$ and walking horizontally with velocity $\dot{x}_d$ in the frame of reference of the driven bath. The notation is defined in the main text. The vertical amplitude of the wave field is exaggerated.}
\label{fig: schematic}
\end{figure}

Consider a liquid droplet of mass $m$ and radius $R$ walking on a bath of the same liquid of density $\rho$, viscosity $\nu$ and surface tension $\sigma$, as shown schematically in Fig.~\ref{fig: schematic}(b). The bath is subjected to vertical time varying acceleration $\gamma(t)$. The system geometry is described in the co-moving frame of the bath by horizontal coordinates $\mathbf{x} = (x,y)$ and the vertical coordinate $z$, with the origin chosen to be the undeformed surface of the bath (dashed horizontal line). In the bath's frame of reference, the center of mass of the droplet is located at a horizontal position $\mathbf{x}_d$ and the south pole of the droplet at a vertical position $z_d$ such that $z_d=0$ represents initiation of the droplet's impact with the undeformed surface of the bath. The free surface elevation of the liquid filling the bath is at $z=h(\mathbf{x},t)$.

\citet{superwalkernumerical} developed a theoretical model for superwalkers by considering two-frequency driving of the form $\gamma(t)=\gamma_f\sin(2\pi f t)+\gamma_{f/2}\sin(\pi f t+\Delta\phi_0)$, where $\gamma_f$ is the amplitude of the primary driving frequency, $\gamma_{f/2}$ is the amplitude of the subharmonic frequency and $\Delta\phi_0$ is the constant relative phase difference between the two. In the experiments of \citet{superwalker}, SGM was observed when the bath was driven at frequencies $f$ and $f/2+\epsilon$ giving rise to time varying acceleration: $\gamma_f\sin(2\pi f t)+\gamma_{f/2+\epsilon}\sin(\pi f t+\Delta\phi_0+2\pi\epsilon t)$.
This can be rewritten as
\begin{equation}\label{drivingform}
\gamma(t)=\gamma_f\sin(2\pi f t)+\gamma_{f/2}\sin(\pi f t+\Delta\phi(t)),
\end{equation}
where $\Delta\phi(t)=\Delta\phi_0+2\pi\epsilon t$, and $\gamma_{f/2} \approx \gamma_{f/2+\epsilon}$ for $\epsilon\ll f/2$. Hence, two frequency driving at $f$ and $f/2+\epsilon$ can be interpreted as driving at $f$ and $f/2$ with a continuously changing phase difference $\Delta\phi(t)$ that varies linearly with time. Assuming that $\epsilon\ll f/2$ so that $\Delta\phi(t)$ evolves slowly, we can use a `quasi-static' approximation and use the model developed by \citet{superwalkernumerical} for a constant phase difference $\Delta\phi_0$ and apply it to the case of time varying phase difference $\Delta\phi(t)$. Hence, we will use the driving acceleration prescribed by Eq.~\eqref{drivingform} in the model of \citet{superwalkernumerical} to study the SGM. We proceed by reviewing the key equations governing the vertical dynamics, horizontal dynamics and wave form in the superwalker model of \citet{superwalkernumerical}.

\subsection{Vertical dynamics}
The vertical motion of the droplet is governed by
\begin{equation}\label{eq: vertical}
    m\ddot{z}_d=-m[\text{g}+\gamma(t)]+{F_N}(t).
\end{equation}
In this equation, the first term on the right hand side is the effective gravitational force acting on the droplet in the oscillating frame of the bath, with g the constant acceleration due to gravity. The second term on the right hand side is the normal force imparted to the droplet during contact with the liquid surface. This contact force is calculated by modeling the bath as a spring and a damper~\cite{molacek_bush_2013},
\begin{equation}
{F_N}(t)=H(-\bar{z}_d) \,\max\left(-k\bar{z}_d-b\dot{\bar{z}}_d,0\right).
\label{springdamp}
\end{equation}
In Eq.~(\ref{springdamp}), $H(\cdot)$ stands for the Heaviside step function and $\bar{z}_d=z_d-h(\mathbf{x}_d,t)$ is the vertical position of the droplet above the free surface of the bath. The constants $k$ and $b$ are the spring constant and damping force coefficient, respectively.

\subsection{Wave field}
The free surface elevation $z=h(\mathbf{x},t)$ is calculated by adding individual waves generated by the droplet on each bounce:
\begin{equation}\label{eq: wave}
\textstyle h(\mathbf{x},t)=\sum_{n} h_n(\mathbf{x},\mathbf{x}_n,t,t_n)\, ,
\end{equation}
where $h_n(\mathbf{x},\mathbf{x}_n,t,t_n)$ is the wave field generated by bounce $n$ at location $\mathbf{x}_n$ and time $t_n$.
The individual waves generated by the droplet on each bounce are localized decaying Faraday waves, which for the case of the bath being driven at frequencies $f$ and $f/2$ take the form
\begin{widetext}
\begin{align}\label{sw_wave}
h_n(\mathbf{x},t)= A_{f/2}\frac{\cos(\pi f t +\theta_{F/2}^{+})}{\sqrt{t-t_n}} \,\text{J}_0(k_{F/2}|\mathbf{x}-{\mathbf{x}}_n|)\,\exp{\left[ -\frac{t-{t}_n}{T_F \text{Me}_{f/2}} - \frac{T_F |\mathbf{x}-{\mathbf{x}}_n|^2}{8 \pi D_{f/2} (t-{t}_n)} \right]}& \nonumber\\
+ A_{f/4}\frac{\cos(\pi f t/2 + \theta_{F/4}^{+})}{\sqrt{t-t_n}} \,\text{J}_0(k_{F/4}|\mathbf{x}-{\mathbf{x}}_n|)\,\exp{\left[ -\frac{t-{t}_n}{T_F\text{Me}_{f/4}} - \frac{T_F |\mathbf{x}-{\mathbf{x}}_n|^2}{8 \pi D_{f/4} (t-t_n)} \right]}&,
\end{align}
\end{widetext}
where the impact location $\mathbf{x}_n$ and the time of impact $t_n$ are given by
\begin{equation*}
    \mathbf{x}_n=\int_{t_n^i}^{t_n^c}\mathbf{x}_d(t')F_N(t')\,\text{d}t'\Big{/}\int_{t_n^i}^{t_n^c}F_N(t')\,\text{d}t',
\end{equation*}
and
\begin{equation*}
    t_n=\int_{t_n^i}^{t_n^c}t'F_N(t')\,\text{d}t'\Big{/}\int_{t_n^i}^{t_n^c}F_N(t')\,\text{d}t',
\end{equation*}
where $t_n^i$ and $t_n^c$ are the time of initiation and completion of the $n$th impact. The interpretation of Eq.~(\ref{sw_wave}) is that a droplet bouncing under the prescribed two-frequency driving excites two dominant subharmonic standing waves with amplitude coefficients $A_{f/2}$ and $A_{f/4}$ and wavenumbers $k_{F/2}$ and $k_{F/4}$, corresponding to  frequencies of $f/2$ and $f/4$. These waves decay in time with a time scale $T_F \text{Me}_{f/2}$ and $T_F \text{Me}_{f/4}$, where $T_F=2/f$ is the period of $f/2$ Faraday waves. The spatial structure of each component of the wave takes the form of a Bessel function of the first kind and zeroth order, $\text{J}_0(\cdot)$.
The waves also spread diffusively with diffusion coefficients $D_{f/2}$ and $D_{f/4}$, and have phase shifts $\theta_{F/2}^{+}$ and $\theta_{F/4}^{+}$ with respect to the driving signal. Due to the evolving phase difference $\Delta\phi(t)$, the phase shift $\theta_{F/4}^{+}(t)$ will vary with time. We refer the reader to \citet{superwalkernumerical} for explicit equations for these parameters. We note that although we have considered both $f/2$ and $f/4$ waves in the wave field of the superwalker, the dominant contribution to the wave field comes from the $f/2$ waves in the parameter regime where superwalking is realized~\citep{superwalkernumerical}.

\subsection{Horizontal dynamics}

The horizontal dynamics of the droplet is modeled using the following equation~\citep{Molacek2013DropsTheory}:
\begin{equation}\label{eq: hor}
    m\ddot{\mathbf{x}}_d+D_{tot}(t)\dot{\mathbf{x}}_d=-{F_N}(t) \boldsymbol{\nabla} h(\mathbf{x}_d,t),
\end{equation}
where $D_{tot}(t)=D_{mom}(t)+D_{air}$. Here, $D_{mom}(t)=C\sqrt{\frac{\rho R}{\sigma}}{F_N}(t)$ is the drag comprising of momentum loss during contact with the bath and $D_{air}=6\pi R \mu_a$ is the air drag. Here $C$ is the contact drag coefficient and $\mu_a$ is the dynamic viscosity of air. The force on the right hand side is the horizontal component of the contact force arising from the small slope $|\nabla h(\mathbf{x}_d,t)|\ll 1$ of the underlying wave field.

\subsection{Numerical implementation}

We restrict the simulations to one horizontal dimension since even if two-dimensional planar dynamics were implemented, the SGM still proceeds along a line.
Eqs.~\eqref{eq: vertical} and \eqref{eq: hor}  corresponding to vertical and horizontal equations of motions respectively are solved using the leap-frog method~\citep{Sprott}, a modified version of the Euler method where the new horizontal and vertical positions of the droplet are calculated using the old velocities and then the new velocities are calculated using the new positions. Converting the second order differential equation for the vertical dynamics in Eq.~(\ref{eq: vertical}) into a system of two first order ordinary differential equations and discretizing using the leap-frog method yields,
\begin{equation*}\label{eq: Vertical disc 1}
    z_d(t_{i+1})=z_d(t_{i})+\Delta t\,v_d(t_{i}),
\end{equation*}
and
\begin{equation*}\label{eq: Vertical disc 2}
    v_d(t_{i+1})=v_d(t_{i})+\frac{\Delta t}{m} \left[-m(\text{g}+\gamma(t_{i+1}))+F_N(t_{i+1})\right],
\end{equation*}
where $v_d(t)=\dot{z}_d(t)$ and
\begin{equation*}\label{eq: Vertical disc 3}
{F_N}(t_{i+1})=H(-\bar{z}_d(t_{i+1})) \,\max\left(-k\bar{z}_d(t_{i+1})-b\bar{v}_d(t_{i}),0\right).
\end{equation*}
Here $\bar{z}_d(t_{i+1})=z_d(t_{i+1})-h(x_d(t_{i+1}),t_{i+1})$ and $\bar{v}_d(t_{i})=v_d(t_{i})-\frac{\partial h}{\partial t}(x_d(t_{i+1}),t_{i+1})$. The total wave height beneath the droplet $h(x_d(t_{i+1}),t_{i+1})$ is calculated using Eq.~\eqref{sw_wave} by keeping the waves from the last $N$ impacts of the droplet. For the integral required to calculate the location of impact $\mathbf{x}_n$, the time of impact $t_n$ and the amplitudes $A_{f/2}$ and $A_{f/4}$ we used the MATLAB inbuilt trapezoid function.

Similarly, Eq.~(\ref{eq: hor}) governing the horizontal dynamics takes the following form,
\begin{equation*}\label{eq: Horizontal disc 1}
    x_d(t_{i+1})=x_d(t_{i})+\Delta t\,u_d(t_{i}),
\end{equation*}
and
\begin{align*}\label{eq: Horizontal disc 2}
    u_d(t_{i+1})=u_d(t_{i})+\frac{\Delta t}{m}\Big[-D_{tot}(t_{i+1})u_d(t_{i})&\nonumber\\
    -{F_N}(t_{i+1}) \frac{\partial h}{\partial x}(x_d(t_{i+1}),t_{i+1})\Big]&,
\end{align*}
where $u_d(t)=\dot{x}_d(t)$.

The physical parameters were fixed to match the experiments of \citet{superwalker}: $\rho=950$\,$\text{kg}/\text{m}^3$, $\nu=20$\,cSt and $\sigma=20.6$\,$\text{mN}/\text{m}$. We also
fix the driving frequencies to $f=80$\,Hz and $f/2=40$\,Hz, and the acceleration amplitudes to $\gamma_{f}=3.8$\,g and $\gamma_{f/2}=0.6$\,g with an initial phase difference $\Delta\phi_0=0^{\circ}$. There are three adjustable parameters in the model: the spring constant of the bath $k$, the damping coefficient of the bath $b$ and the dimensionless contact drag coefficient $C$. The corresponding dimensionless parameters are given by $K=k/m\omega_d^2$ and $B=b/m\omega_d$, where $\omega_d=\sqrt{\sigma/\rho R^3}$ is the droplet's characteristic oscillation frequency~\cite{molacek_bush_2013}. We choose $C=0.17$, $K=0.70$ and $B=0.60$ as these values give a good fit to the experimental data for small- to moderate-sized superwalkers~\citep{superwalkernumerical}.

To increase computational speed, we only stored the waves generated by the $N=100$ most recent bounces of the droplet and discarded the earlier ones, which have typically decayed to below $10^{-5}$ of their initial amplitude for the chosen parameters. A time step of $\Delta t =T_F/100$ was used. The simulations were initialized with $x_d=0\,\text{mm}$, $u_d=1\,\text{mm/s}$, $v_d=0$\,mm/s and three different vertical positions $z_d= (0,5,10)R$. Multiple initial conditions were used so that different modes of walking and bouncing existing for the same parameter values are likely to be captured.

\section{Emergence of Stop-and-go motion (SGM)}\label{emergence}

\begin{figure*}
\centering
\includegraphics[width=2\columnwidth]{Fig_2.pdf}
\caption{Droplet dynamics as a function of the relative phase difference. (a) Steady walking speed $u_d$ of a superwalker as a function of a fixed phase difference $\Delta\phi_0$ for a droplet of radius $R=0.60$\,mm from experiments of \citet{superwalker} (markers) and simulations of \citet{superwalkernumerical} (multi-colored curve). The experimental data to the right of the vertical dashed line at $\Delta\phi_0=180^\circ$ is repeated. The style of marker indicates the bouncing modes observed in experiments: $(1,2,2)$ are black circles $\bullet$, $(1,2,1)$ are black triangles {\protect \scalebox{0.8}{$\blacktriangle$}}, transition between a $(1,2,1)$ and a $(1,2,2)$ mode are gray squares {\protect \scalebox{0.5}{$\blacksquare$}}. The colors on the solid curve indicate different vertical bouncing mode observed in simulations: $(1,2,1)^{\text{L}}${} in navy blue, $(1,2,2)$ in red, and chaotic modes that arise at the transition between $(1,2,1)^{\text{L}}${} and $(1,2,2)$ are shown in gray. The light purple region indicates the bouncing regime where a zero superwalking speed, in other words pure bouncing, is observed. The dotted multi-colored horizontal line in the superwalking regime indicates pure bouncing states that are unstable to horizontal perturbations. In this panel, a time-varying phase difference $\Delta\phi(t)=2\pi \epsilon t$ would correspond to traversing the $\Delta\phi_0$ axis from either left to right $(\epsilon>0)$ or right to left $(\epsilon<0)$. Panel (b) shows the vertical periodic motion of the bath, $\mathcal{B}(t)=-(\gamma_f/(2\pi f)^2)\sin(2\pi f t)-(\gamma_{f/2}/(\pi f)^2)\sin(\pi f t+\Delta\phi_0)$, in the lab frame at different constant phase differences $\Delta\phi_0$. Panel (c) shows vertical bouncing modes obtained for different values of $\Delta\phi_0$ that are marked by a black line segment at the top of panel (a). In this panel, the solid black curves indicate the bath motion, $\mathcal{B}(t)$, the multi-colored curves represent the vertical motion of the south pole of the droplet, $z_d(t)+\mathcal{B}(t)$, and the filled blue regions illustrate the vertical motion of the liquid surface directly beneath the droplet, $h(x_d,t)+\mathcal{B}(t)$, all in the lab frame. The gray regions indicate times at which the droplet is in contact with the bath.}
\label{Fig: stopgo1}
\end{figure*}

A given sized superwalker at fixed driving amplitudes and frequencies can have a range of walking speeds $u_d$ depending on the phase difference $\Delta\phi_0$. This walking speed dependence on the phase difference of superwalkers from experiments (markers) and numerical simulations (multi-colored curve) is shown for a typical droplet of radius $R=0.60$\,mm in Fig.~\ref{Fig: stopgo1}(a). We find two different regimes for such a droplet: a bouncing regime where the droplet bounces in place with no horizontal motion, and a superwalking regime where the droplet also self-propels horizontally while bouncing. However, we note that in the superwalking regime, there exist purely bouncing states (dotted multi-colored horizontal line in the superwalking regime) that are unstable to horizontal perturbations. The two regimes realized for $R=0.60$\,mm droplet is typical for small- to moderate-sized superwalkers~\citep{superwalker}. Note that due to the form of driving considered here at frequencies $f$ and $f/2$, the driving signal, see Fig.~\ref{Fig: stopgo1}(b), and hence the walking speed dependence of the phase difference in the region $180^{\circ}\leq \Delta\phi_0 \leq 360^{\circ}$ is a repeat from the region $0^{\circ}\leq \Delta\phi_0 \leq 180^{\circ}$.

To describe the vertical dynamics of the droplet, we follow \citet{superwalker} and use the notation $(l,m,n)$ to indicate that the droplet impacts the surface $n$ times during $m$ oscillation periods of the bath at frequency $f$, which equals $l$ oscillation periods of the bath at frequency $f/2$. In the superwalking regime, we find that the droplet is typically bouncing in a $(1,2,1)^\text{L}$ mode, where the superscript `L' denotes a low-bouncing, long-contact mode compared to a high-bouncing, short-contact mode that is also observed for smaller superwalkers~\citep{superwalkernumerical}. The $(1,2,1)$ bouncing mode is crucial for walking as the droplet in this mode is bouncing at the same frequency as the frequency of the subharmonic Faraday waves that emerge beyond the Faraday instability threshold. Thus, the droplet's bouncing is in resonance with the damped Faraday waves it generates and with which it interacts. In the bouncing regime and at the start of the superwalking regime, we also find a $(1,2,2)$ bouncing mode where the droplet contacts the bath twice, typically a high bounce and a low bounce, every two up-and-down cycles of the bath. The bouncing modes from simulations match well with experiments, however, for some values of $\Delta\phi_0$, we find $(1,2,1)^{\text{L}}${} mode where a $(1,2,2)$ is observed in experiments. We note that it is difficult to distinguish between a $(1,2,1)^{\text{L}}${} and a $(1,2,2)$ mode in experiments and hence it is not clear whether all the $(1,2,2)$ modes observed in experiments of \citet{superwalker} are truly $(1,2,2)$ or some in fact may even be $(1,2,1)^{\text{L}}${}. The predicted vertical dynamics of the droplet for a selection of $\Delta\phi_0$ values is shown in Fig.~\ref{Fig: stopgo1}(c).

By allowing the phase difference $\Delta\phi_0$ to vary slowly according to $\Delta\phi(t)=2\pi\epsilon t$, we can continuously traverse the $\Delta\phi_0$ axis in Fig.~\ref{Fig: stopgo1}(a) either from left to right ($\epsilon>0$) or right to left ($\epsilon<0$). This gives rise to periodic traversals of the bouncing regime and the superwalking regime resulting in the emergence of SGM. The rate of traversal is controlled by $\epsilon$ with the period of one cycle of traversal, in other words one SGM cycle, given by $1/|\epsilon|$ corresponding to phases sweeping from $0^{\circ}$ to $360^{\circ}$.

Once the phase difference is allowed to evolve with time, the system is subjected to three characteristic time scales set by the input parameters: (i) the time scale, $T_F$, of the bath driving or the oscillation of the underlying dominant $f/2$ standing waves generated by the droplet on each bounce, (ii) the memory time scale, $T_F \text{Me}_{f/2}$, of temporal decay of the dominant $f/2$ waves generated by the droplet and (iii) the SGM time scale, $1/|\epsilon|$, introduced by the evolving phase difference $\Delta\phi(t)$. The interactions of these time scales along with the time scale for the inertial response of the droplet, give rise to the emergent droplet dynamics where the vertical bouncing dynamics of the droplet occurs on the time scale (i), and the evolution of the horizontal walking dynamics of the droplet occurs on the time scale (iii).

 \begin{figure*}
\centering
\includegraphics[width=2\columnwidth]{Fig_3.pdf}
\caption{Various types of SGM observed from simulations in the $(\epsilon,R)$ parameter space. Four distinct types of dynamics are observed: (i) No walking ($\times$, beige region) where the droplet bounces vertically but does not propel horizontally, (ii) irregular SGM ($\circ$, purple region) where the droplet erratically switches its walking direction and/or traverses unequal distances after every half SGM cycle, (iii) back-and-forth SGM ($\Delta$, blue region) where the droplet switches direction after every half SGM cycle and (iv) forth-and-forth SGM (+, maroon region) where the droplet moves in the same direction after every half SGM cycle. The markers indicate the parameter values where simulations were performed while the background colors represent the interpolated region.}
\label{Fig: stopgo_PS}
\end{figure*}

\section{Parameter-space description}\label{PS}

\begin{figure*}
\centering
\includegraphics[width=2\columnwidth]{Fig_4.pdf}
\caption{Back-and-forth SGM: (a) Horizontal walking speed $u_d$ (blue curve) as a function of the evolving phase difference $\Delta\phi(t)=2\pi\epsilon t$ for one cycle of the SGM. The walking speed from constant phase difference superwalker simulations from Fig.~\ref{Fig: stopgo1}(a) is shown as a multi-colored curve. The light purple region indicates the pure bouncing regime. Panel (b) shows the vertical velocity $v_d$ of the droplet at impact (lower frame) and lift-off (upper frame) from constant phase superwalker simulations (background blue and red multi-colored markers with $\bullet$ indicating the superwalking solution branch and $\times$ indicating the non-walking, pure bouncing solution branch in the superwalking regime) and the back-and-forth SGM (blue $\bullet$ indicating impact and blue {\protect \scalebox{0.8}{$\blacktriangle$}} indicating lift-off). Panel (c) shows the vertical motion of the bath (black curve), wave height beneath the droplet (blue shaded region) and the vertical location of the droplet's `south pole' (blue curve) during back-and-forth SGM. The gray shaded region represents contact between the droplet and the underlying wave. Panel (d) shows the evolution of the horizontal wave force $F_{W}=-{F_N}(t)\partial h/\partial x$ (blue), horizontal drag force $F_{D}=-D_{tot}(t)\dot{x}_d$ (red) and the total horizontal force $F_{H}=F_{W}+F_{D}$ (black) during the back-and-forth SGM. Panel (e) shows the evolution of the gradient of the wave beneath the droplet $\partial h/\partial x$ (black solid curve) with the dashed line indicating zero gradient. The droplet size is fixed to $R=0.60$\,mm and the detuning is fixed to $\epsilon=-0.5$\,Hz.}
\label{Fig: backforth1}
\end{figure*}

We have simulated droplets in the parameter space formed by the detuning parameter $\epsilon$ that prescribed the rate at which the different phase differences are traversed and the droplet radius $R$. For the chosen parameters, the time scale for the bouncing dynamics or the oscillations of the waves generated by the droplet is $T_F=2/f=0.025$\,s, while the memory time scale is $T_F \text{Me}_{f/2} \approx 0.22$\,s. We explore the parameter space in the range $-1\leq \epsilon \leq 1$\,Hz and $0.30 \leq R \leq 0.65$\,mm. We have chosen $|\epsilon|\leq1$\,Hz to ensure $\epsilon\ll f/2$ so that our quasi-static approximation remains valid. The minimum droplet size is set by the start of the walking regime of droplets as droplets smaller than this are typically unable to walk at these parameters~\citep{superwalker,superwalkernumerical}. We restrict the maximum droplet size to $R=0.65$\,mm as the theoretical model used here for superwalkers fails to capture experimentally observed characteristics for larger droplets~\citep{superwalkernumerical}. The different dynamical behaviors observed in the parameter-space plot are shown in Fig.~\ref{Fig: stopgo_PS}. We have observed four qualitatively different types of dynamics in the $(\epsilon,R)$ parameter space: (i) no walking, (ii) irregular SGM, (iii) back-and-forth SGM and (iv) forth-and-forth SGM. In the irregular SGM, we find that the droplet erratically reverses its walking direction and/or travels a variable distance after every half cycle of the SGM. In the back-and-forth SGM, intriguingly, the droplet reverses its walking direction after every half cycle of the SGM resulting in the droplet oscillating back-and-forth about a fixed mean position or the droplet performs back-and-forth oscillations with a drift in the mean position. Conversely, in the forth-and-forth SGM, the droplet maintains the same walking direction after every half cycle of the SGM.

The parameter-space diagram in Fig.~\ref{Fig: stopgo_PS} shows that relatively small droplets ($0.30\lesssim R \lesssim 0.50$\,mm) are either unable to walk under this prescribed driving or undergo irregular SGM for small negative detunings. Relatively larger droplets ($0.50\lesssim R \lesssim 0.65$\,mm) show a wide range of stop-and-go behavior including back-and-forth and forth-and-forth motion. For these we typically observe only back-and-forth SGM for relatively large magnitude negative detunings ($-1.0 \lesssim \epsilon \lesssim -0.2\,$Hz) and typically forth-and-forth SGM for relatively large magnitude positive detunings ($0.2 \lesssim \epsilon \lesssim 0.7\,$Hz). For both positive and negative detunings that are small in magnitude, we see a mix of back-and-forth, forth-and-forth and irregular SGM. Moreover, for very small magnitudes of detunings $|\epsilon|\leq0.04$\,Hz, we find that the SGM ceases. To understand these different behaviors in more detail, we study a fixed droplet size of $R=0.60$\,mm, which shows a diversity of behaviors. We start by exploring the back-and-forth SGM in Sec.~\ref{BF}, followed by forth-and-forth SGM in Sec.~\ref{FF} and the small detuning regime in Sec.~\ref{small detun}.

\begin{figure*}
\centering
\includegraphics[width=2\columnwidth]{Fig_5.pdf}
\caption{Forth-and-forth SGM: (a) Horizontal walking speed $u_d$ (maroon curve) as a function of the evolving phase difference $\Delta\phi(t)=2\pi\epsilon t$ for one cycle of the SGM. The walking speed from constant phase difference superwalker simulations from Fig.~\ref{Fig: stopgo1}(a) is shown as a multi-colored curve. The light purple region indicates the bouncing regime. Panel (b) shows the vertical velocity $v_d$ of the droplet at impact (lower frame) and lift-off (upper frame) from constant phase superwalker simulations (background blue and red multi-colored markers with $\bullet$ indicating the superwalking solution branch and $\times$ indicating the non-walking, pure bouncing solution branch in the superwalking regime) and the forth-and-forth SGM (maroon $\bullet$ indicating impact and maroon {\protect \scalebox{0.8}{$\blacktriangle$}} indicating lift-off). Panel (c) shows the vertical motion of the bath (black curve), wave height beneath the droplet (blue shaded region) and the vertical location of the droplet's `south pole' (maroon curve) during forth-and-forth SGM. The gray shaded region represents contact between the droplet and the underlying wave. Panel (d) shows the evolution of the horizontal wave force $F_{W}=-{F_N}(t)\partial h/\partial x$ (blue), horizontal drag force $F_{D}=-D_{tot}(t)\dot{x}_d$ (red) and the total horizontal force $F_{H}=F_{W}+F_{D}$ (black) during the forth-and-forth SGM. Panel (e) shows the evolution of the gradient of the wave beneath the droplet $\partial h/\partial x$ (black solid curve) with the dashed line indicating zero gradient. The droplet size is fixed to $R=0.60$\,mm and the detuning is fixed to $\epsilon=0.5$\,Hz.}
\label{Fig: forthforth1}
\end{figure*}

\section{Back-and-forth SGM}\label{BF}

For relatively large droplets and typically negative detunings in Fig.~\ref{Fig: stopgo_PS}, we find that in the SGM, the droplet reverses its walking direction after every half SGM cycle resulting in a back-and-forth motion. One cycle of a typical back-and-forth SGM for $R=0.60$\,mm and $\epsilon=-0.5$\,Hz is shown in Fig.~\ref{Fig: backforth1} (see also Supplementary Material~\citep{supplement} for a video). During the acceleration period of this SGM, the droplet bounces in a $(1,2,1)^{\text{L}}$ mode where the droplet is efficiently able to generate a slowly decaying localized Faraday waves on each bounce and thus build up the wave field that propels it, see Fig.~\ref{Fig: backforth1}(c) and (e). In this period, during contact the droplet receives two horizontal impulses in the direction of its motion and one horizontal impulse in the opposite direction such that the overall horizontal force, sum of the force from the underlying wave field and the net drag force, is in the direction of the droplet's motion, see Fig.~\ref{Fig: backforth1}(d).
Since the overall horizontal force is imparted in the direction of motion of the droplet, the droplet accelerates. Note from Fig.~\ref{Fig: backforth1}(a) that during this acceleration stage, the droplet's instantaneous velocity $u_d$ (blue curve) is lower than the corresponding steady superwalking velocity for a corresponding constant phase difference driving (multi-colored curve).

Once the instantaneous speed of the droplet exceeds the corresponding steady superwalking speed for a constant phase difference driving, further evolution in this bouncing mode results in an overall horizontal force in the direction opposite to $u_d$ causing the droplet to decelerate.
This is soon followed by the droplet transitioning to a $(1,2,2)$ mode where the droplet impacts the bath twice in two up-and-down motions of the bath. In this mode, during one period of the bath, the droplet generates two waves of opposite phase that interfere destructively and thus their contribution to building up the overall wave field is small. This results in the decay of the previously built-up wave field. At this time in the deceleration stage, the droplet is in the bouncing regime and two competing effects are taking place simultaneously: (i) the horizontal force is acting opposite to $u_d$ and will slow down the droplet and eventually reverse its walking direction, and (ii) the magnitude of the horizontal force is decreasing due to the decay of the wave field. In the back-and-forth SGM, the effect of the horizontal force acting opposite to $u_d$ dominates the decay of the wave field and the droplet reverses its walking direction before the wave field has decayed sufficiently. Thus, the droplet is now walking in the opposite direction at a low speed. Eventually the wave field has decayed significantly and the small horizontal force does not have much further impact on the motion of the droplet.

As $\Delta \phi(t)$ continues to evolve, the droplet transitions back to a $(1,2,1)^{\text{L}}$ mode again and it enters the superwalking regime. The droplet begins to build up the wave field and accelerate again. However, since the droplet was already traveling at a low speed in the opposite direction before the acceleration re-started, the droplet continues to accelerate in the opposite direction. This cycle repeats periodically and results in the back-and-forth SGM.

We note that after reversing the walking direction, the droplet does not immediately accelerate even though it has entered the stage where the constant $\Delta\phi_0$ equilibrium solution has a large superwalking speed. By inspecting the plot of vertical velocity of the droplet $v_d$ at impact and lift-off as a function of the evolving phase difference $\Delta\phi(t)$ in Fig.~\ref{Fig: backforth1}(b), we see that during this time, the droplet is following the equilibrium solution branch corresponding to non-walking solution in the superwalking regime and the droplet begins to accelerate only once it has left this solution branch. This may be due to the droplet losing its wave field during the deceleration stage and hence it takes time for the droplet to build up its wave field again and become unstable to the pure bouncing state.

For a typical droplet size of $R=0.60$\,mm, such a back-and-forth SGM is realized for a range of negative detunings $-1\leq\epsilon\lesssim -0.2\,$Hz. The horizontal walking velocity $u_d$ and the vertical velocity $v_d$ at impact and lift-off for different detunings are shown in Figs.~\ref{Fig: sGM diff detun}(a) and (c) respectively. We find that as the magnitude of detuning gets smaller in the range $-1\leq\epsilon\lesssim -0.2$\,Hz, the droplet follows the constant phase difference superwalking branch more closely during the deceleration stage.

For small negative detunings $(-0.2\lesssim\epsilon< 0\,\text{Hz})$, the SGM is qualitatively different and we see a mixture of back-and-forth, forth-and-forth and irregular SGM. These are explored in detail in Sec.~\ref{small detun}. Here we note that in the back-and-forth SGM in the small negative detuning regime, see for example the $\epsilon=-0.07$\,Hz curves in Figs.~\ref{Fig: sGM diff detun}(a) and (c), we find that in addition to the droplet following the constant phase difference superwalking branch more closely during deceleration, the droplet only starts accelerating near a phase difference of $\Delta\phi=135^{\circ}$. This also corresponds to a sharp transition from the non-walking solution branch to the superwalking solution branch in Fig.~\ref{Fig: sGM diff detun}(c).

\section{Forth-and-forth SGM}\label{FF}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig_6.pdf}
\caption{Back-and-forth and forth-and-forth SGM for different detunings. Horizontal walking speed $u_d$ as a function of the evolving phase difference $\Delta\phi(t)$ for one cycle of the SGM showing (a) back-and-forth SGM at different negative detunings and (b) forth-and-forth SGM at different positive detunings. Panels (c) and (d) show the vertical velocity of the droplet $v_d$ at impact (lower frame) and lift-off (upper frame) from the back-and-forth and the forth-and-forth SGM corresponding to panels (a) and (b) respectively. The constant phase difference superwalker simulation results are shown as multi-colored curves in the background in these panels.}
\label{Fig: sGM diff detun}
\end{figure}

\begin{figure*}
\centering
\includegraphics[width=2\columnwidth]{Fig_7.pdf}
\caption{Back-and-forth and forth-and-forth SGM for small detunings for a droplet of radius $R=0.60$\,mm. (a) Walking velocity $u_d$ as a function of time for $\epsilon=-0.10$\,Hz showing forth-and-forth SGM without noise (black) and SGM with erratic direction reversals due to added noise (gray). Panel (b) shows a zoomed in section of panel (a) with the vertical axis stretched and the green background highlighting that the droplet continues to walk in the same direction as the previous cycle when no noise is added. (c) Walking velocity $u_d$ as a function of time for $\epsilon=0.08$\,Hz showing back-and-forth SGM without noise (black) and SGM with erratic direction reversals due to added noise (gray). Panel (d) shows a zoomed in section of panel (c) with the vertical axis stretched and the pink and green regions highlighting that the droplet reverses its walking direction compared to the previous cycle when no noise is added. The multi-colored curves represent the corresponding results for constant phase difference simulations.}
\label{Fig: nowalking1}
\end{figure*}

For relatively large droplets and typically positive detunings in Fig.~\ref{Fig: stopgo_PS}, we find that in the SGM, the droplet maintains its walking direction after every half SGM cycle resulting in a forth-and-forth SGM. One cycle of a typical forth-and-forth SGM for $R=0.60$\,mm and $\epsilon=0.5$\,Hz is shown in Fig.~\ref{Fig: forthforth1} (see also Supplementary Material~\citep{supplement} for a video). During the acceleration stage of this SGM, we also find that the droplet bounces in a $(1,2,1)^{\text{L}}$ mode. However, this $(1,2,1)^{\text{L}}$ is qualitatively different from the one observed for back-and-forth SGM, see Fig.~\ref{Fig: backforth1}(c). Here the droplet's vertical motion has three turning points when the droplet is in contact with the bath as compared to a single turning point for the $(1,2,1)^{\text{L}}$ mode in back-and-forth SGM. This difference can be attributed to the fact that in this SGM, as shown in Fig.~\ref{Fig: forthforth1}(b), the vertical velocity of the droplet is unable to completely leave the non-walking branch and transition to the superwalking branch where a $(1,2,1)^{\text{L}}$ mode with only a single turning point during contact would be observed. Instead, the droplet's vertical velocity during the acceleration stage oscillates about the non-walking branch where a qualitatively different $(1,2,1)^{\text{L}}$ mode is realized. However, similar to the back-and-forth SGM, the droplet receives two impulses in the direction of motion and one impulse in the direction opposite to its motion such that the overall horizontal force acting on the droplet is in the direction of its motion allowing it to accelerate. Moreover, during acceleration, the instantaneous velocity of the droplet is smaller than the steady superwalking velocity at constant phase differences.

Once the droplet's instantaneous speed exceeds the corresponding equilibrium speed for constant phase differences, the droplet starts to decelerate. In contrast to the back-and-forth SGM, here we find that the droplet's wave field decays significantly during the initial stages of the deceleration period in the $(1,2,1)^{\text{L}}$ mode. Hence, the horizontal component of the force arising from the gradient of the underlying wave does not a have a significant impact on the droplet. At this time the droplet also transitions to a $(1,2,2)$ mode in which the droplet is unable to build its wave field. Since the droplet is already near the peak velocity at the start of the deceleration stage, it experiences large drag forces that eventually slow it down. But unlike the back-and-forth SGM, the droplet's walking direction does not get reversed as the drag force can only slow the droplet down and not reverse its walking direction. Hence, the droplet continues to walk with a decaying speed in the same direction as the $\Delta\phi(t)$ passes through the bouncing regime. Eventually, the droplet enters the phase differences for the superwalking regime where the droplet's steady superwalking velocity exceeds its instantaneous velocity and the droplet also transitions to a $(1,2,1)^{\text{L}}$ mode, where it can build up its wave field and start to accelerate again. Since the droplet already had an initial perturbation in the direction of motion before it entered the acceleration stage, the motion continues in the same direction resulting in forth-and-forth SGM.

We find that after the velocity of the droplet has decayed significantly and is very small, the droplet does not immediately accelerate even though it has reached the phase differences where the equilibrium solution has a large superwalking speed. Here we find that, in comparison to the back-and-forth SGM, the droplet takes a longer time to build up the wave field and accelerate and thus the peak speed it can reach before it starts to decelerate is smaller compared to back-and-forth SGM. This again can be attributed to the droplet never completely transiting to the superwalking branch in Fig.~\ref{Fig: backforth1}(b), and thus being unable to accelerate efficiently.

For a typical droplet size of $R=0.60$\,mm, such a forth-and-forth SGM is realized for a range of positive detunings $0.2\lesssim\epsilon\leq0.7$\,Hz. The horizontal walking velocity $u_d$ and the vertical velocity $v_d$ at impact and lift-off for different positive detunings is shown in Figs.~\ref{Fig: sGM diff detun}(b) and (d) respectively. We find that as the magnitude of detuning increases from $0.2$\,Hz to $0.7$\,Hz, the droplet is able to achieve slightly higher peak speed during the `go' phase. This is also correlated with the increase in the magnitude of oscillations around the non-walking branch in the superwalking regime as the detuning increases.

For small positive detunings $(0<\epsilon\lesssim 0.2\,\text{Hz})$, the SGM is again qualitatively different and we see a mixture of back-and-forth, forth-and-forth and irregular SGM. It is explored in detail in Sec.~\ref{small detun}. Here we note that in the forth-and-forth motion that is realized for small positive detunings (see $\epsilon=0.05$\,Hz curves in Figs.~\ref{Fig: sGM diff detun}(b) and (d)), we find that the droplet is able to completely leave the non-walking branch and transition to the superwalking branch and thus achieve a significantly higher superwalking speed in the forth-and-forth SGM compared to the ones described above.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{Fig_8.pdf}
\caption{No walking for very small detunings $(|\epsilon|<0.04)$ and the effect of added noise for a droplet of radius $R=0.60$\,mm. Panels (a) and (c) show the walking velocity $u_d$ as a function of time without noise (black) and with noise (gray) for $\epsilon=-0.03$\,Hz and $\epsilon=0.03$\,Hz respectively. Panels (b) and (d) shows the vertical velocity of the droplet $v_d$ at impact and lift-off without noise (black) and with noise (gray) corresponding to panels (a) and (c) respectively. The constant phase superwalker simulations results are shown as multi-colored curves in these panels.}
\label{Fig: nowalking2}
\end{figure}

\section{SGM for small detuning}\label{small detun}

For both positive and negative detunings that are relatively small in magnitude ($-0.2\lesssim\epsilon\lesssim-0.05$\,Hz and $0.05\lesssim \epsilon \lesssim 0.2$\,Hz), we observe a mixture of back-and-forth, forth-and-forth and irregular SGM for a typical droplet size of $R=0.60$\,mm. However, we find that these SGMs are qualitatively different from the ones that emerge at larger magnitudes of positive and negative detunings.

In the irregular SGM that emerges in this regime for small detunings, the droplet travels the same distance every half SGM cycle while the reversals of direction are erratic. In contrast, the irregular SGM that arises at large detunings in the parameter space for relatively small droplets, the droplet traverses unequal path lengths as well as performing erratic direction reversals every half cycle of the SGM.

The back-and-forth SGM discussed in Sec.~\ref{BF} reverses its walking direction only once every half cycle and then continues moving in the opposite direction while the forth-and-forth SGM discussed in Sec.~\ref{FF} never reverses its walking direction. As shown in Figs.~\ref{Fig: nowalking1}(a-d), for both back-and-forth and forth-and-forth SGM at small detunings, we find that the droplet reverses its direction multiple times when traversing the bouncing regime. For this regime of SGM, when the droplet first reaches the phase differences corresponding to the bouncing regime, the droplet decelerates. Since the detuning is small, the bouncing regime is traversed slowly and the droplet first overshoots the zero velocity while decelerating and reverses it walking direction. But the droplet cannot accelerate yet in the opposite direction since it is still in the bouncing regime, and therefore the droplet approaches the zero velocity solution again and overshoots, resulting in a second reversal of the walking direction. These reversals of direction continue until the droplet reaches the start of the superwalking regime where it can begin to build up the wave field again and then accelerate. Whichever direction the droplet happens to move when it encounters the start of the superwalking regime, it continues to walk and accelerate in that direction. If the number of zero velocity crossings are even, then a forth-and-forth SGM emerges and if the number of crossings are odd, then back-and-forth SGM is realized.

By performing simulations with added noise in the form of a horizontal random force sampled from a uniform distribution in the range $F_{noise}=[-10^{-7},10^{-7}]$\,N, much smaller in magnitude than the typical forces experienced by the droplet from the underlying wave, we find that we are able to destroy these correlated back-and-forth and forth-and-forth motions at small detunings (see gray curves in Fig.~\ref{Fig: nowalking1}). In these simulations, an SGM with irregular direction reversal emerges due to the random external noise added to the system. In typical experiments where small noise sources are to be expected, back-and-forth and forth-and-forth SGM may not be realized at small detunings. We also note that for these simulations with added noise, the droplet's vertical velocity is able to leave the non-walking branch more easily and following the superwalking branch closely.

For very small detunings $|\epsilon|\leq 0.04$\,Hz, we find that the SGM ceases to exist and the droplet remains stationary. Fig.~\ref{Fig: nowalking2}(a-d) shows the walking velocity $u_d$ as a function of time and the vertical velocity of the droplet $v_d$ at impact and lift-off as a function of the evolving phase difference. Here we see that the droplet always follows the non-walking branch for all phase differences, which means that the droplet is not able to transition to the superwalking branch and build up the wave field it needs to propel itself. However, we find that adding noise in the horizontal dynamics, similar to Fig.~\ref{Fig: nowalking1}, destroys this fragile non-walking state and we again observe SGM, with random switches in the walking direction resulting from the noise.

\section{Conclusions}\label{DC}

We have numerically investigated a new type of intermittent locomotion, the stop-and-go motion (SGM) of superwalking droplets, that emerges when a liquid bath is driven simultaneously at frequencies $f$ and $f/2$ with an evolving phase difference $\Delta\phi(t)=2\pi \epsilon t$. The SGM is a complex nonlinear phenomenon with multiple timescales coming into play such as the bouncing time scale of the droplet, the memory time scale associated with decay of Faraday waves, the even longer time scale introduced by the detuning and the time scale of the inertial response of the droplet. We have observed three qualitatively different kinds of SGM in simulations: back-and-forth, forth-and-forth and irregular. We note that motion similar to irregular SGM reported here has also been observed in experiments with two frequency driving with frequencies $80$\,Hz and $64$\,Hz~\citep{PhysRevE.94.053112}. In preliminary experiments with the setup described in \citet{superwalker}, we have frequently observed back-and-forth SGM. We have also observed irregular SGM often, and less frequently what appears to be forth-and-forth SGM. However, it is difficult to differentiate between a forth-and-forth and irregular SGM in such an experimental setup as the horizontal extent of the domain is limited and a droplet undergoing a forth-and-forth SGM would encounter the boundary after only a few cycles of the SGM making it difficult to conclude whether it is truly forth-and-forth motion. A detailed experimental study of SGM is needed to reveal the different types of SGM that may be realized in laboratory.  An annular domain might be appropriate to study forth-and-forth SGM experimentally.

Investigating the back-and-forth SGM for negative detunings that are relatively large in magnitude reveals the mechanism for the switch in the walking direction that occurs after every half SGM cycle. In this SGM, the droplet receives persistent impulses in the direction opposite to its motion during its deceleration stage. These persistent impulses eventually slow down the droplet and reverse its walking direction before the underlying wave field decays significantly. Hence, the droplet already has an initial perturbation in the direction opposite to the direction of motion when it encounters the superwalking regime again and it is able to accelerate in the opposite direction for the next cycle of SGM. We note that reversals of motion have also been demonstrated in the system of single-frequency driven walkers where a pulse in the driving signal was engineered to abruptly introduce a $\pi$ shift in the bouncing phase of the droplet, which results in the droplet retracing its previous trajectory~\citep{perrard2016}. Here we observe this behavior arising as an emergent phenomenon.

Investigating the forth-and-forth SGM for large positive detunings reveals the mechanism for the droplet to maintain its walking direction after every half SGM cycle. In this SGM, the droplet's wave field decays significantly during the initial stages of the deceleration period in the bouncing regime and hence for the remainder of the deceleration stage, only the drag force is acting on the droplet. The drag slows down the droplet but cannot reverse the walking direction so when the droplet encounters the superwalking regime again, it accelerates in the same direction. It is worth noting that in both the back-and-forth and forth-and-forth SGM, the periodic reversal or maintenance of the walking direction is not directly due to the memory of the waves that are generated by the droplet since the wave field decays significantly during each cycle for both kinds of SGM.

For small positive and negative values of detunings, we found a mixture of back-and-forth, forth-and-forth and irregular SGM. The back-and-forth and the forth-and-forth here are qualitatively different from the ones at large magnitude detunings. Here we find that the droplet switches its direction multiple times while traversing across the bouncing regime. Moreover, these SGM also tend to be less robust and adding a small noise force destroy the back-and-forth and forth-and-forth SGM in this regime.

In this work, we considered SGM that is inspired by the slight detuning of the two driving frequencies in superwalker experiments that results in the phase difference varying linearly in time. Using the back-and-forth and forth-and-forth SGM as building blocks, one might be able to control and program the droplet's motion by engineering any slowly varying time-dependent functions $\Delta\phi(t)$. Moreover, walkers and superwalkers are known to exhibit rich multi-droplet behaviors. It would be interesting to investigate such multi-droplet interactions in the SGM regime where the continuously changing phase difference will result in periodically varying inter-droplet interactions that may give rise to novel collective behaviors.

\acknowledgements
We acknowledge financial support from an Australian Government Research
Training Program (RTP) Scholarship (R.V.) and the Australian Research Council via the Future Fellowship Project No.\ FT180100020 (T.S.).


\section{The Functional $\mathcal{F}_N$}

In this section we define a functional $\mathcal{F}_N : \mathbb{C}[z,\bar{z}]\rightarrow \mathbb{C}$ and use it to represent matrix products.
The goal is to transfer the matrix positivity conditions into block Toeplitz positivity conditions.

\subsection{Definition and Basic Properties}
\begin{defn}
 For $f(z,\overline{z}) \in \mathbb{C}[z,\overline{z}]$, define:
 \begin{equation}
   \label{eq:52}
 \mathcal{F}_N(f) :=
 \frac{1}{2\pi} \int_{-\pi}^\pi [ \frac{1}{N^2} \sum_{j,k=0}^{N-1} f( \omega^j e^{i\theta} , \overline{\omega^k e^{i\theta}} ) ] d\theta
 \end{equation}

\end{defn}

Observe that the integrand is the average of the entries of $\textrm{Gram}(f)[e^{i\theta},\dots,\omega^{N-1}e^{i\theta}]$.

We require the following computation for the next proposition.

\begin{lem}\label{cha:main-results:circulant-average-lemma}
  \begin{equation}
    \label{eq:51}
  \frac{1}{N^2} \sum_{j,k=0}^{N-1} \omega^{\ell (j-k)}
  =
  \begin{cases}
   1 & \text{if } \ell \mid N \\
   0 & \text{if } \ell \nmid N
  \end{cases}
  \end{equation}

\end{lem}
\proof

Define the symbol:

\begin{equation}
  \label{eq:6f}
    \mu(\ell) :=
  \begin{cases}
    N & \textrm{ if } N \mid \ell \\
    0 & \textrm{ if } N \nmid \ell \\
  \end{cases}
\end{equation}
Then
\begin{align}
  \sum_{j=0}^{N-1} \omega^{\ell j} &= 1 + \omega^\ell + \omega^{2\ell} + \dots + \omega^{\ell(N-1)} \\
  &= \mu(\ell)
\end{align}
which gives
\begin{align}
  \sum_{j,k=0}^{N-1} \omega^{\ell(j-k)} &= \sum_{j=0}^{N-1} \sum_{k=0}^{N-1} \omega^{\ell j}\omega^{-\ell k} \\
                                        &= \sum_{j=0}^{N-1} \omega^{\ell j} ( \sum_{k=0}^{N-1} \omega^{-\ell k} ) \\
                                        &= \sum_{j=0}^{N-1} \omega^{\ell j} \mu(\ell) \\
  &= \mu(\ell)^2
\end{align}
\endproof

The key idea is that $\mathcal{F}_N(f)$ computes the sum of the diagonal entries $a_{\ell \ell}$ of the coefficient matrix of $f$ such that $\ell$ is a multiple of $N$.

\begin{prop}\label{sec:defin-basic-prop-1:computation_of_F}
 Let $f(z,\bar{z}) = \sum a_{jk} \bar{z}^j z^k \in \mathbb{C}_[z,\bar{z}]$.
 Then:
 \begin{equation}
   \label{eq:53}
  \mathcal{F}_N( f ) = a_{0,0} + a_{N,N} + a_{2N,2N} + \cdots =  \sum_{\substack{\ell = 0 \\ N \mid \ell}}^m a_{\ell \ell}
 \end{equation}

\end{prop}
\begin{proof}
\begin{align}
\mathcal{F}_N(f)
   &= \frac{1}{2\pi} \int_{-\pi}^\pi [ \frac{1}{N^2} \sum_{j,k=0}^{N-1} f(\omega^j e^{i\theta},\overline{\omega^k e^{i\theta}} ) ] d\theta \\
   &= \frac{1}{N^2} \sum_{j,k=0} [ \frac{1}{2\pi} \int_{-\pi}^\pi f(\omega^j e^{i\theta}, \overline{\omega^k d^{i\theta}}) d\theta ] \\
   &= \frac{1}{N^2} \sum_{j,k=0}^{N-1} [ \sum_{\ell = 0}^m a_{\ell \ell} \omega^{\ell(j-k)} ] \\
   &= \sum_{\ell=0}^m a_{\ell \ell} [ \frac{1}{N^2} \sum_{j,k=0}^{N-1} \omega^{\ell(j-k)} ] \\
   &= \sum_{\substack{\ell = 0 \\ N \mid \ell}}^m a_{\ell \ell}
\end{align}
where we use Lemma \ref{cha:main-results:diagonals-lemma} and Lemma \ref{cha:main-results:circulant-average-lemma}.
\end{proof}

\begin{cor} $\mathcal{F}_N$ satisfies the following properties:
  \begin{enumerate}
  \item $\mathcal{F}_N( z^N \bar{z}^N f(z,\bar{z})) = \mathcal{F}_N(f(z,\bar{z}))$ for all $f \in \mathbb{C}_h[z,\bar{z}]$
  \item $\mathcal{F}_N$ is a $\mathbb{C}$-linear map $\mathbb{C}[z,\bar{z}] \rightarrow \mathbb{C}$
  \item $\mathcal{F}_N$ is an $\mathbb{R}$-linear map $\mathbb{C}_h[z,\bar{z}] \rightarrow \mathbb{R}$
  \end{enumerate}
\end{cor}
\proof
Follows from Proposition \ref{sec:defin-basic-prop-1:computation_of_F}.
\endproof

\begin{prop}\label{sec:defin-basic-prop:functional_positivity}
 Let $f \in \mathbb{C}_h[z,\bar{z}]$.
 Suppose $\text{Gram}(f)[e^{i\theta},\omega e^{i\theta},\dots,\omega^{N-1}e^{i\theta}] \geq 0$ for all $\theta \in [0,2\pi]$.
 Then:

 \begin{equation}
   \label{eq:7f}
  \mathcal{F}_N(f(z,\bar{z})) \geq 0
 \end{equation}

 Furthermore:
 \begin{equation}
   \label{eq:8f}
  \mathcal{F}_N( |h(z)|^2 f(z,\bar{z}) ) \geq 0
 \end{equation}

 for all holomorphic polynomials $h(z) \in \mathbb{C}[z]$.
\end{prop}
\begin{proof}

The average of the entries of a positive matrix is positive.  Thus:
\begin{equation}
  \label{eq:9f}
  \frac{1}{N^2} \sum_{j,k=0}^{N-1} f( \omega^j e^{i\theta} , \overline{\omega^k e^{i\theta}} )  \geq 0 \qquad \text{ for all } \theta \in [0,2\pi]
\end{equation}

The average of positive numbers is positive.  Thus:
\begin{equation}
  \label{eq:10f}
  \mathcal{F}_N(f) :=
  \int_{-\pi}^\pi [ \frac{1}{N^2} \sum_{j,k=0}^{N-1} f( \omega^j e^{i\theta} , \overline{\omega^k e^{i\theta}} ) ] \frac{d\theta}{2\pi} \geq 0
\end{equation}

For a Hermitian polynomial $g \in \mathbb{C}_h[z,\bar{z}]$, for ease of notation, let $\text{Gram}(g)$ denote the matrix $\text{Gram}(g)[e^{i\theta},\omega e^{i\theta},\dots,\omega^{N-1}e^{i\theta}] $.
We know $\text{Gram}(|h(z)|^2) \geq 0$ for all holomorphic polynomials $h(z) \in \mathbb{C}[z]$.
Then
\begin{equation}
  \label{eq:11f}
 \text{Gram}(|h(z)|^2 f(z,\bar{z}) ) = \text{Gram}(|h(z)|^2) \circ \text{Gram}(f(z,\bar{z})) \geq 0
\end{equation}
by the Schur product theorem (where $\circ$ denotes the Hadamard product).

Then $\mathcal{F}_N(|h(z)|^2f(z,\bar{z})) \geq 0$ follows from the first part of the proof.
\end{proof}

\subsection{Representing Matrix Products with $\mathcal{F}_N$}

Given $v,w \in \mathbb{C}^N$ and $A \in M_N$, our goal is to construct a polynomial $f$ such that $\mathcal{F}_N(f) = v^\ast A w$.
Let us expand the expression $v^\ast A w$ so that we may recognize its appearance later.
 Write
 \[ v = \begin{bmatrix} v_0 & \cdots & v_{N-1} \end{bmatrix}^T \]
 \[ w = \begin{bmatrix} w_0 & \cdots & w_{N-1} \end{bmatrix}^T \]
 \[ A = [ a_{jk} ]_{j,k=0}^{N-1} \]
 Then
 \begin{equation}
   \label{eq:54}
  v^\ast A w = \sum_{j,k=0}^{N-1} a_{jk} \bar{v}_j w_k
 \end{equation}

The following example demonstrates the construction for the case $N=2$.

\begin{ex}
 Let $v,w \in \mathbb{C}^2$ and $A \in \mathbb{C}^{2 \times 2}$.  Write:
 \[ v = \begin{bmatrix} v_0 \\ v_1 \end{bmatrix} \qquad
    w = \begin{bmatrix} w_0 \\ w_1 \end{bmatrix} \qquad
    A = \begin{bmatrix} a_{00} & a_{01} \\ a_{10} & a_{11} \end{bmatrix}
 \]
 Define
 \[  f(z,\bar{z}) = \sum_{j,k=0}^1 a_{jk} \bar{z}^j z^k \]
 \[  \tilde{v}(z) = v_0 z^2 + v_1 z \]
 \[  \tilde{w}(z) = w_0 z^2 + w_1 z \]
 Then
 \[
   \overline{ \tilde{v}(z)} \tilde{w}(z) f(z,\bar{z}) =
   a_{00} \bar{v}_0 w_0 z^2 \bar{z}^2 + za_{01}\bar{v}_0w_1 z \bar{z}^2 + \bar{z} a_{10} \bar{v}_1 w_0 z^2 \bar{z} + z\bar{z} a_{11} \bar{v}_1 w_1 z \bar{z}
 \]
 Then
 \[
  \mathcal{F}_2( \overline{ \tilde{v}(z) } \tilde{w}(z) f(z,\bar{z} )) = a_{00} \bar{v}_0 w_0 + a_{01} \bar{v}_0 w_1 + a_{10} \bar{v}_1 w_0 + a_{11}\bar{w}_1 w_1 = v^\ast A w
 \]
 by Proposition \ref{sec:defin-basic-prop-1:computation_of_F}.

 Written in matrix notation: (consider the coefficient of $z^2 \bar{z}^2$ before applying $\mathcal{F}_2$)
 \begin{equation}
   \label{eq:55}
   \mathcal{F}_2 (
   \begin{bmatrix} z^2 \\ z \end{bmatrix}^\ast
   \begin{bmatrix} \bar{v}_0 \\ \bar{v}_1 \end{bmatrix}
   \begin{bmatrix} \bar{w}_0 \\ \bar{w}_1 \end{bmatrix}^\ast
   \begin{bmatrix} z^2 \\ z \end{bmatrix}
   \begin{bmatrix} 1 \\ z \end{bmatrix}^\ast
   \begin{bmatrix} a_{00} & a_{01}\\
   a_{10} & a_{11} \end{bmatrix}
 \begin{bmatrix} 1 \\ z \end{bmatrix})
 =
 \begin{bmatrix} v_0 \\ v_1 \end{bmatrix}^\ast
    \begin{bmatrix} a_{00} & a_{01}\\
      a_{10} & a_{11} \end{bmatrix}
     \begin{bmatrix} w_0 \\ w_1 \end{bmatrix}
 \end{equation}
\end{ex}

We generalize the previous example to obtain the desired construction:

\begin{prop}
  Let:
  \[ v = \begin{bmatrix} v_0 & \cdots & v_{N-1} \end{bmatrix}^T \in \mathbb{C}^N   \]
  \[ w = \begin{bmatrix} w_0 & \cdots & w_{N-1} \end{bmatrix}^T \in \mathbb{C}^N   \]
  \[ A = [ a_{jk} ]_{j,k=0}^{N-1} \in \mathbb{C}^{N\times N} \]
  Define:
  \begin{equation}
    \label{eq:14f}
  \tilde{v}(z) = \sum_{j=0}^{N-1} v_j z^{N-j} \in \mathbb{C}[z]
  \end{equation}
  \begin{equation}
    \label{eq:21f}
   \tilde{w}(z) = \sum_{j=0}^{N-1} w_j z^{N-j} \in \mathbb{C}[z]
  \end{equation}

  \begin{equation}
    \label{eq:24f}
   f(z,\bar{z}) = \sum_{j,k=0}^{N-1} a_{jk} \bar{z}^j z^k \in \mathbb{C}[z,\bar{z}]
  \end{equation}
  Then:
  \begin{equation}
    \label{eq:12f}
   \mathcal{F}_N( \overline{ \tilde{v}(z) } \tilde{w}(z) f(z,\bar{z} )) = v^\ast A w
  \end{equation}
  Furthermore, for integers $s,t \geq 0$ we have:
 \begin{equation}
   \label{eq:13f}
   \mathcal{F}_N( \bar{z}^{Ns} z^{Nt} \overline{ \tilde{v}(z) } \tilde{w}(z) f(z,\bar{z} )) = \delta_{st} v^\ast A w
 \end{equation}

\end{prop}
\begin{proof}
 We have
 \[
  \overline{\tilde{v}(z)} \tilde{w}(z) = \sum_{j,k=0}^{N-1} \bar{v}_j w_k \bar{z}^{N-j} z^{N-k}
 \]
 and
 \[
  f(z,\bar{z}) = \sum_{c,d=0}^{N-1} a_{cd} \bar{z}^c z^d
 \]
 Thus
 \[
  \overline{\tilde{v}(z)} \tilde{w}(z) f(z,\bar{z}) = \sum_{j,k,c,d=0}^{N-1} a_{cd} \bar{v}_j w_k \bar{z}^{N-j+c} z^{N-k+d}
 \]
Use the $\mathbb{C}$-linearity of $\mathcal{F}_N$:
\[
 \mathcal{F}_N(\overline{\tilde{v}(z)} \tilde{w}(z) f(z,\bar{z})) = \sum_{j,k,c,d=0}^{N-1} a_{cd} \bar{v}_j w_k \mathcal{F}_N(\bar{z}^{N-j+c} z^{N-k+d})
\]
Since $0\leq j,k,c,d \leq N-1$, applying Proposition we get
\[
\mathcal{F}_N(\bar{z}^{N-j+c} z^{N-k+d}) = \delta_{jc}\delta_{kd}
\]
Therefore only terms with $j=c$ and $k=d$ will survive $\mathcal{F}_N$:
\[
 \mathcal{F}_N(\overline{\tilde{v}(z)} \tilde{w}(z) f(z,\bar{z})) = \sum_{j,k=0}^{N-1} a_{jk}\bar{v}_j w_k = v^\ast A w
\]

In order to prove the second part, for a Hermitian polynomial $f$, let $\textrm{Mon}(f) \subset \mathbb{N} \times \mathbb{N}$ denote the monomial support (the set of monomials corresponding to nonzero coefficients of $f$).
Note $\textrm{Mon}(fg) = \textrm{Mon}(f) + \textrm{Mon}(g)$.
Say a monomial in $\mathbb{N} \times \mathbb{N}$ is \emph{$N$-divisible} if both components are divisible by $N$ and say the monomial is \emph{diagonal} if both entries are equal.
We know that $\mathcal{F}_N(f)$ is the sum of the coefficients corresponding to monomials of $f$ which are both diagonal and $N$-divisible.

For integers $a,b \in \mathbb{Z}$ let $[a,b] := \{ j \in \mathbb{Z} | a \leq j \leq b \}$.
Then:
\begin{equation}
  \label{eq:57}
  \textrm{Mon}( \overline{ \tilde{v}(z)} \tilde{w}(z) ) = [1,N] \times [1,N]
\end{equation}
\begin{equation}
  \label{eq:58}
  \textrm{Mon} ( f(z,\bar{z}) = [0,N-1] \times [0,N-1]
\end{equation}
Hence
\begin{equation}
  \label{eq:59}
  \textrm{Mon} (  \overline{ \tilde{v}(z)} \tilde{w}(z) f(z,\bar{z}) ) =[1,2N-1]\times[1,2N-1]
\end{equation}
whose only diagonal $N$-divisible element is $(N,N)$.

We have
\begin{equation}
  \label{eq:61}
  \textrm{Mon} (  \overline{ \tilde{v}(z)} \tilde{w}(z) f(z,\bar{z}) \bar{z}^{Ns} z^{Nt} )
  =
  [1+Ns, 2N-1+Ns] \times [ 1+Nt, 2N-1+Nt]
\end{equation}

For convenience, denote $Q(s,t) = Q(s,t)(z,\bar{z}) =  \overline{ \tilde{v}(z)} \tilde{w}(z) f(z,\bar{z}) \bar{z}^{Ns} z^{Nt}$.

As we vary $(s,t) \in \mathbb{Z}_{\geq 0} \times \mathbb{Z}_{\geq 0}$ we obtain a partition of $\mathbb{Z}_{>0} \times \mathbb{Z}_{>0}$ by squares $\textrm{Mon}(Q(s,t))$.  Only the bottom right monomial of each tile is $N$-divisible.
If $s=t$, then the bottom right element of $\textrm{Mon}(Q(s,t))$ is both $N$-divisible and diagonal.
If $s\neq t$, then no elements of $\textrm{Mon} (Q(s,t))$ are both $N$-divisible and diagonal.
From these observations the result follows.

\end{proof}

\section{Introduction}

In this paper we study the problem of writing a Hermitian polynomial as a Hermitian sum of squares modulo the Hermitian ideal $(z^N \bar{z}^N -1)$, where $N$ is a positive integer.
We investigate an idea of Putinar and Scheiderer from \cite{MR2729629} which constructs a counterexample to a question of D'Angelo.
The idea provides necessary conditions on the Hermitian polynomial, and in our situation we show that they are sufficient.
We begin by stating the problem, discussing the context, describing the idea to obtain the matrix positivity conditions, and formulating the main result.

\subsection{Problem Statement}

Let $\mathbb{C}[z,\bar{z}]$ denote the polynomial algebra in the variables $z=(z_1,\dots,z_n)$ and $\bar{z}=(\bar{z}_1,\dots,\bar{z}_n)$ with coefficients in $\mathbb{C}$.
We can write each element $f \in \mathbb{C}[z,\bar{z}]$ using standard multinomial notation:
\begin{equation}
  \label{eq:1}
  f(z,\bar{z}) = \sum_{\alpha,\beta } a_{\alpha \beta} \bar{z}^\alpha z^\beta
\end{equation}
Or in dimension one ($n=1$) using matrix notation:
\begin{equation}
  \label{eq:3}
  f(z,\bar{z}) =
  \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^d
  \end{bmatrix}^\ast
  \begin{bmatrix}
    a_{00} & a_{01} & \cdots & a_{0d} \\
    a_{10} & a_{11} & \cdots & a_{1d} \\
    \vdots & \vdots & \ddots & \vdots \\
    a_{d0} & a_{d1} & \cdots & a_{dd}
  \end{bmatrix}
  \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^d
  \end{bmatrix}
  = \psi_d(z)^\ast A \psi_d(z)
\end{equation}
where
\begin{equation}
  \label{eq:4}
  \psi_d(z) =
  \begin{bmatrix}
    1 & z & \cdots & z^d
  \end{bmatrix}^T
\end{equation}
denotes the \emph{tautological monomial map} and $A=[a_{jk}]$ is an $(d+1)\times (d+1)$ matrix with complex entries.
We often write $\psi(z)$ for $\psi_d(z)$ when $d$ is understood.

Here we use standard matrix notation and operations.  If $A$ is an $n \times m$ matrix, then $A^T$ denotes the transpose, $\overline{A}$ denotes the conjugate, and $A^*$ denotes the conjugate transpose.  We think of elements $v \in \mathbb{C}^n$ as column vectors and $\langle v , w \rangle = v^\ast w$ denotes the standard inner product on complex Euclidean space, with the convention being conjugate linear in the first component.

A matrix is \emph{Hermitian} if $A=A^\ast$ and positive semidefinite if $v^\ast A v \geq 0$ for all vectors $v$.  We use the term \emph{positive matrix} when $A$ is positive semidefinite and write $A \geq 0$.

We have an involution $f \mapsto f^\ast$ given by conjugation:
\begin{equation}
  \label{eq:2}
  f^\ast(z,\bar{z}) := \overline{f(z,\bar{z})} = \psi(z)^\ast A^\ast \psi(z)
\end{equation}
Say $f$ is \emph{Hermitian} if $f=f^\ast$.
Hermitian polynomials are real-valued on $\mathbb{C}^n$.
Let $\mathbb{C}_h[z,\bar{z}]$ denote the collection of Hermitian polynomials.
An ideal $I$ in $\mathbb{C}_h[z,\bar{z}]$ is called a \emph{Hermitian ideal}.
If $h(z) \in \mathbb{C}[z]$ is a holomorphic polynomial, then
\begin{equation}
  \label{eq:5}
  |h(z)|^2 = h(z)\overline{h(z)} \in \mathbb{C}_h[z,\bar{z}]
\end{equation}
is a \emph{Hermitian square}.
Let $\Sigma^2_h \subset \mathbb{C}_h[z,\bar{z}]$ denote the collection of finite sums of Hermitian squares.
We have the basic test for a Hermitian polynomial $f$ to be a Hermitian sum of squares: $f \in \Sigma^2_h$ if and only if the Hermitian coefficient matrix is positive.

We now state the main problem:
  \begin{prob}
    Suppose $f$ is a Hermitian polynomial and $I$ is a Hermitian ideal.
    Under what conditions on $f$ and $I$ does there exist an identity:
\begin{equation}
      f(z,\bar{z}) = \sum_{j=1}^\ell |h_j(z)|^2 + g(z,\bar{z})
\end{equation}
    where $h_1,\dots,h_\ell$ are holomorphic polynomials and $g \in I$.
  \end{prob}

  We easily obtain the \emph{trivial necessary condition},  point-wise positivity on the zero set of $I$:
\begin{equation}
    f(p,\bar{p}) \geq 0 \qquad \forall p \in \mathcal{Z}(I)
\end{equation}
  where $\mathcal{Z} (I)$ denotes the zero-set of $I$:
\begin{equation}
    \mathcal{Z} (I) = \{ p \in \mathbb{C}^n \ | \ g(p,\bar{p}) = 0 \ \forall g \in \mathcal{Z}( I) \}
\end{equation}

This condition is not sufficient in general.  For example, consider the ideal $I=(0)$ and the polynomial $f(z,\bar{z}) = (z+\bar{z})^2$.
Since $z+\bar{z} = 2\textrm{Re}(z)$, we see $f(p,\bar{p}) \geq 0$ for all $p \in \mathbb{C}$.
Writing $f$ in matrix notation:
\begin{equation}
  \label{eq:7b}
  (z+\bar{z})^2 =     \begin{bmatrix}
      1 \\ z \\ z^2
    \end{bmatrix}^\ast
    \begin{bmatrix}
      0 & 0 & 1 \\
      0 & 2 & 0 \\
      1 & 0 & 0
    \end{bmatrix}
    \begin{bmatrix}
      1 \\ z \\ z^2
    \end{bmatrix}
\end{equation}
The coefficient matrix is not positive, hence $f \notin \Sigma^2_h$.

\subsection{Context}

In dimension one, for the ideal $I=(z\bar{z}-1)$, the trivial necessary condition is sufficient.  The idea goes back to the classical Riesz-Fejer lemma concerning positive trigonometric polynomials on the circle \cite{MR1580923} \cite{MR1580922}.
The first multivariable result appears in \cite{MR233770} showing that the trivial necessary condition is sufficient for the odd-dimensional sphere in $\mathbb{C}^n$.
This theorem was rediscovered in \cite{MR1386836} by D'Angelo and collaborators.
At the 2006 AIM Conference ``CR Complexity Theory'' in Palo Alto, D'Angelo asks the following generalization of Quillen's theorem:

\emph{    If $f$ is a Hermitian polynomial on a pseudoconvex hypersurface, then does $f$ agree with a Hermitian sum of squares along the hypersurface?}

  The answer to this paper lies in the paper \cite{MR2729629}.
  Two main ideas of this paper:
  \begin{enumerate}
  \item An obstruction to the question is introduced and a counterexample is constructed.  The obstruction provides necessary matrix positivity conditions for the Hermitian polynomial $f$ beyond pointwise positivity given by the trivial necessary condition.
  \item A characterization of the Hermitian ideals for which every positive Hermitian polynomial is a Hermitian sum of squares is obtained using the Archimedean Positivstellensatz of real algebra.
  \end{enumerate}

  The first point is explored further in the paper \cite{MR2925474}, and the second point is explored further in the paper \cite{MR3426234}.

  Our goal is to continue the idea of the matrix positivity conditions for Hermitian polynomials, and obtain equivalent characterizations for Hermitian sums of squares modulo Hermitian ideals, thereby going beyond the Archimedean positivstellensatz and pointwise-positivity conditions.

\subsection{Matrix Positivity Conditions for Hermitian Polynomials}
  We describe the idea of matrix positivity conditions used in \cite{MR2729629} to construct the counterexample to D'Angelo's question.  Two remarks are needed before we begin.

  First, it is notationally convenient to collect the terms of a Hermitian sum of squares into the components of a holomorphic polynomial map.
  A \emph{holomorphic polynomial map} is a map $h: \mathbb{C}^n \rightarrow \mathbb{C}^\ell$ whose components are holomorphic polynomials.
  Every Hermitian sum of squares may be written as the squared norm of a holomorphic polynomial map:
  \begin{equation}
    \label{eq:15}
    |h_1(z)|^2 + \dots + |h_\ell(z)|^2 = \langle h(z), h(z) \rangle
  \end{equation}
  where $h(z) = (h_1(z),\dots,h_\ell(z))$ is a holomorphic polynomial map.

  Second, for a polynomial $f(z,\bar{z})\in \mathbb{C}[z,\bar{z}]$, we may polarize and treat $z$,$\bar{z}$ as independent variables.  If $f$ is Hermitian, then the polarization satisfies the Hermitian symmetric condition:
  \begin{equation}
    \label{eq:6b}
    f(z,\bar{w}) = \overline{f(w,\bar{z})}
  \end{equation}

  Now given a Hermitian polynomial $f$ and a Hermitian ideal $I$, we show how to obtain necessary matrix positivity conditions which must be satisfied by $f$.

  Suppose $f$ is a Hermitian sum of squares modulo $I$:
  \begin{equation}
    \label{eq:16}
    f(z,\bar{z}) = \langle h(z), h(z) \rangle + g(z,\bar{z})
  \end{equation}
  where $h$ is a holomorphic polynomial map and $g \in I$.

  Suppose further that $p_1,\dots,p_\ell \in \mathbb{C}^n$ are points such that
  \begin{equation}
    \label{eq:17}
    g(p_j, \bar{p}_k) = 0 \qquad \forall g \in I \qquad \forall j,k=1,\dots,\ell
  \end{equation}
  Then
  \begin{equation}
    \label{eq:18}
    f(p_j,\bar{p}_k) = \langle h(p_k), h(p_j) \rangle \qquad \forall j,k=1,\dots,\ell
  \end{equation}
  The matrix determined by the right side is a matrix of pairwise inner products (a \emph{Gram matrix}), which is always positive, hence
  \begin{equation}
    \label{eq:19}
    [  f(p_j,\bar{p}_k)]_{j,k=1}^\ell \geq 0
  \end{equation}
  Thus we conclude:
  \begin{lem}[Matrix Positivity Conditions]\label{MATRIX-POSITIVITY-CONDITIONS}
    Let $f$ be a Hermitian polynomial and $I$ a Hermitian ideal.
    If $f \in \Sigma^2_h+I$, then $f$ satisfies \eqref{eq:19} for all collections of points $p_1,\dots,p_\ell$ satisfying \eqref{eq:17}
  \end{lem}

  Using this idea we obtain necessary conditions for $f \in \Sigma^2_h + I$.
  Our main question of interest is whether these conditions are sufficient?

  With these ideas in mind, given a Hermitian polynomial $f$ and points $p_1,\dots,p_\ell$, we define the \emph{Gram matrix of $f$ on $p_1,\dots,p_\ell$} as the matrix of pairwise evaluations:
  \begin{equation}
    \label{eq:8b}
    \textrm{Gram}(f)[p_1,\dots,p_\ell] := [ f(p_j,\bar{p_k})]_{j,k=1}^\ell
  \end{equation}

  As an example, consider $f=(z+\bar{z})^2$ and $I=(0)$.
  Choose $p_1 = 0$ and $p_2 = 1$.
  Then consider the Gram matrix of $f$ on $p_1,p_2$:
  \begin{equation}
    \label{eq:20}
    \textrm{Gram}(f)[p_1,p_2] =
    \begin{bmatrix}
      f(p_1,\bar{p}_1) & f(p_1,\bar{p}_2) \\
      f(p_2,\bar{p}_1) & f(p_2,\bar{p}_2)
    \end{bmatrix}
    =
    \begin{bmatrix}
      0 & 1 \\
      1 & 4
    \end{bmatrix}
  \end{equation}
  The matrix has determinant $-1$, hence is not positive.
  Therefore $f \notin \Sigma^2_h$.

  \subsection{Main Results}

  Our main result is to prove sufficiency of the matrix positivity conditions for the ideal $    I = (z^N \bar{z}^N - 1 )$, where $N$ is a positive integer.
  The conditions are obtained as follows:
  Let $P(z,\bar{z}) = z^N \bar{z}^N -1$ and $\omega = e^{2\pi i / N}$.
  Then:
  \begin{equation}
    \label{eq:22}
    P ( \omega^j \xi , \overline{\omega^k \xi} ) = 0 \qquad \forall \xi \in \mathbb{T} \qquad \forall j,k=0,\dots,N-1
  \end{equation}
  Therefore by Lemma \ref{MATRIX-POSITIVITY-CONDITIONS}, if $f \in \Sigma^2_h + (z^N \bar{z}^N - 1)$, then
  \begin{equation}
    \label{eq:23}
    [ f(  \omega^j \xi , \overline{\omega^k \xi} ) ]_{j,k=0}^{N-1} \geq 0 \qquad \forall \xi \in \mathbb{T}
  \end{equation}
  We will show that these conditions are sufficient.
  The case $N=1$ is the classical Riesz-Fejer lemma.

The next example provides evidence to the investigation.

\renewcommand{\labelenumi}{(\roman{enumi})}
\begin{ex} ($N=2$)
We demonstrate an example of a Hermitian polynomial $f$ such that:
  \begin{enumerate}
  \item $f \notin \Sigma^2_h$
  \item $\textrm{Gram}(f)[e^{i\theta},-e^{i\theta}] \geq 0$ for all $\theta$
  \item $f \in \Sigma^2_h + (z^2\bar{z}^2-1)$
  \end{enumerate}

  Consider:
  \begin{align}
    f(z,\bar{z}) &= 10 + 2z + 2\bar{z} + 10z\bar{z} -2z^2 \bar{z} -2z\bar{z}^2 \\
    &=
      \begin{bmatrix}
        1 \\ z \\ z^2
      \end{bmatrix}^\ast
    \begin{bmatrix}
      10 & 2 & 0 \\
      2 & 10 & -2 \\
      0 & -2 & 0
    \end{bmatrix}
       \begin{bmatrix}
        1 \\ z \\ z^2
      \end{bmatrix}
  \end{align}
  \begin{enumerate}
  \item By Sylvester's criterion, the Hermitian coefficient matrix is not positive, hence $f \notin \Sigma^2_h$.
  \item We compute $\textrm{Gram}(f)[e^{i\theta},-e^{i\theta}]$:
    \begin{equation}
      \label{eq:3a}
      \textrm{Gram}(f)[e^{i\theta},-e^{i\theta}] =
      \begin{bmatrix}
        20 & 8i\sin(\theta) \\
        -8i\sin(\theta) & 20
      \end{bmatrix}
    \end{equation}
    We have:
    \begin{equation}
      \label{eq:2a}
      \det \textrm{Gram} (f)[e^{i\theta},-e^{i\theta}] = 400 -64\sin^2(\theta) > 0
    \end{equation}
    Thus by Sylvester's criterion:
    \begin{equation}
      \label{eq:111}
      \textrm{Gram}(f)[e^{i\theta},-e^{i\theta}] \geq 0 \qquad \forall \theta \in [0,2\pi]
    \end{equation}

  \item Now
    \begin{equation}
      \label{eq:4a}
      f(z,\bar{z})+ 5(z^2\bar{z}^2-1) =
       \begin{bmatrix}
        1 \\ z \\ z^2
      \end{bmatrix}^\ast
      \begin{bmatrix}
        5 & 2 & 0 \\
        2 & 10 & -2 \\
        0 & -2 & 5
      \end{bmatrix}
 \begin{bmatrix}
        1 \\ z \\ z^2
      \end{bmatrix}
    \end{equation}
    By Sylvester's criterion, the Hermitian coefficient matrix is positive.
    Thus:
    \begin{equation}
      \label{eq:5a}
      f \in \Sigma^2_h + (z^2 \bar{z}^2-1)
    \end{equation}

  \end{enumerate}
\end{ex}

We now describe the main result.
Consider a Hermitian polynomial $g(z,\bar{z})$ and assume $[g(\omega^j \xi, \overline{\omega^k \xi})]_{j,k=0}^{N-1} \geq 0$ for all $\xi \in \mathbb{T}$.  We may write $g$ in block matrix form:
\begin{equation}
    g(z,\bar{z}) =
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}^\ast
  \begin{bmatrix}
    B_{00} & B_{01} & \cdots & B_{0m} \\
    B_{10} & B_{11} & \cdots & B_{1m} \\
    \vdots & \vdots & \ddots & \vdots \\
    B_{m0} & B_{m1} & \cdots & B_{mm}
  \end{bmatrix}
    \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}
\end{equation}
Reducing coefficients modulo $(z^N\bar{z}^N-1)$, we have that $g$ is equivalent to a unique polynomial of the form
\begin{equation}
    f(z,\bar{z}) =
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}^\ast
  \begin{bmatrix}
   A_0 & A_1 & \cdots & A_m \\
   A_{-1} & 0 &  \cdots & 0 \\
   \vdots & \vdots & \ddots & \vdots \\
   A_{-m} & 0 & \cdots & 0
  \end{bmatrix}
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}
\end{equation}
without changing the positivity condition.
Our plan is to show that the matrix positivity conditions imply positivity of the associated block Toeplitz form:
\begin{equation}
     \textrm{Toep}(A_0,\dots,A_m) :=
  \begin{bmatrix}
    A_0 & A_1 & A_2 & \cdots & A_m \\
    A_{-1} & A_0 & A_1 & \ddots & \vdots \\
    A_{-2} & A_{-1} & A_0 & \ddots & \vdots \\
    \vdots & \ddots & \ddots & \ddots & A_1 \\
    A_{-m} & \cdots & \cdots & A_{-1} & A_0
  \end{bmatrix}
\end{equation}
We then use the block Toeplitz positivity to invoke the operator-valued Riesz-Fejer theorem \cite{MR2743422} to obtain $f \in \Sigma^2_h +(z^N\bar{z}^N-1)$, and hence $g \in \Sigma^2_h + (z^N \bar{z}^N-1)$ as well.  The result is formulated as Theorem \ref{MAIN-THEOREm}.

\subsection{Outline}

The paper is organized as follows:
In section 2 we prove some basic lemmas and computations involving orthogonal polynomials on the circle inspired by the $N=1$ case.  We also introduce the block trace parametrization of trigonometric polynomials and recall the operator-valued RF theorem for our purposes.
In section 3 we define a functional $\mathcal{F}_N$, then discuss basic properties and how to use the functional to represent matrix inner products.
In section 4 we utilize the developed tools to prove the full characterization of Hermitian sums of squares modulo the Hermitian ideal $I=(z^N\bar{z}^N-1)$.

\section{Preliminaries}

In this section we discuss preliminary ideas for the main result: basic orthogonal polynomial computations on the circle, block trace parametrization of trigonometric polynomials, and the operator-valued Riesz-Fejer theorem.

\subsection{Basic Computations on the Circle}

To begin, let $h_j(z) = z^j$ for $j \in \mathbb{Z}$.
We have the classical formulae:
  \begin{equation}
    \label{eq:48}
    \frac{1}{2\pi} \int_{-\pi}^\pi h_j(e^{i\theta}) d\theta = \delta_{j0}
  \end{equation}
and
\begin{equation}
  \label{eq:49}
    \frac{1}{2\pi} \int_{-\pi}^\pi h_j(e^{i\theta}) \overline{h_k(e^{i\theta})} d\theta = \delta_{jk}
\end{equation}

  where $\delta_{jk}$ is the standard Kronecker symbol:
  \[
    \delta_{jk} =
    \begin{cases}
      1 & \textrm{ if } j=k \\
      0 & \textrm{ if } j\neq k \\
    \end{cases}
  \]

  We also note the identities:
  \begin{equation}
    \label{eq:6a}
    \overline{h_j(e^{i\theta})} = h_{-j}(e^{i\theta})
  \end{equation}
  \begin{equation}
    \label{eq:7a}
    h_j(z)h_k(z) = h_{j+k}(z)
  \end{equation}

Given $f \in \mathbb{C}[z,\bar{z}]$, we can integrate around the circle to obtain the trace of the coefficient matrix:

  \begin{lem}\label{cha:main-results:diagonals-lemma}
 Let $f(z,\bar{z}) = \sum_{j,k=0}^m a_{jk} \bar{z}^j z^k \in \mathbb{C}[z,\bar{z}]$.
 Then:
 \begin{equation}
   \label{eq:50}
  \frac{1}{2\pi} \int_{-\pi}^\pi f( e^{i\theta}, \overline{e^{i\theta}} ) d\theta = \sum_{\ell = 0}^m a_{\ell \ell}
 \end{equation}

\end{lem}
\proof
\begin{align}
  \frac{1}{2\pi} \int_{-\pi}^\pi f(e^{i\theta}, \overline{e^{i\theta}}) d \theta &=  \frac{1}{2\pi} \int_{-\pi}^\pi ( \sum_{j,k=0}^m a_{jk} e^{-ij\theta}e^{ik\theta}  ) d\theta \\
                                                                                &= \sum_{j,k=0}^m a_{jk} (  \frac{1}{2\pi} \int_{-\pi}^\pi e^{i\theta(k-j)} d\theta ) \\
                                                                                &= \sum_{j,k=0}^m a_{jk} \delta_{jk} \\
  &= \sum_{\ell=0}^m a_{\ell \ell}
\end{align}
\endproof

Consider a Hermitian polynomial $g(z,\bar{z})$ written in matrix notation:
\begin{equation}
  \label{eq:8a}
  g(z,\bar{z}) =
  \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^m
  \end{bmatrix}^\ast
  \begin{bmatrix}
    b_{00} & b_{01} & \cdots & b_{0m} \\
    b_{10} & b_{11} & \cdots & b_{1m} \\
    \vdots & \vdots & \ddots & \vdots \\
    b_{m0} & b_{m1} & \cdots & b_{mm}
  \end{bmatrix}
   \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^m
  \end{bmatrix}
\end{equation}

By reducing the coefficients modulo $(z\bar{z}-1)$, we see that $g$ is congruent modulo $(z\bar{z}-1)$ to a unique polynomial of the form:
\begin{equation}
  \label{eq:9a}
  f(z,\bar{z}) =
   \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^m
  \end{bmatrix}^\ast
  \begin{bmatrix}
    a_0 & a_1 & \cdots & a_m \\
    a_{-1} & 0 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    a_{-m} & 0 & \cdots & 0
  \end{bmatrix}
   \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^m
  \end{bmatrix}
\end{equation}
where $a_{-j} = \bar{a}_j$.
A Hermitian polynomial of this form is called a \emph{trigonometric polynomial} with data $(a_0,\dots,a_m)$.

Since $g(z,\bar{z}) = f(z,\bar{z}) + q(z,\bar{z})(z\bar{z}-1)$ for some Hermitian polynomial $q(z,\bar{z})$, we have the following observations:
\begin{enumerate}
\item $f(e^{i\theta},\overline{e^{i\theta}}) \geq 0$ for all $\theta \in [0,2\pi]$ if and only if $g(e^{i\theta},\overline{e^{i\theta}}) \geq 0$ for all $\theta \in [0,2\pi]$
\item $f \in \Sigma^2_h + (z\bar{z}-1)$ if and only if $g \in \Sigma^2_h + (z\bar{z}-1)$
\end{enumerate}

It is known \cite[pg. 17]{MR890515} that the condition
\begin{equation}
  \label{eq:10a}
  f(e^{i\theta},\overline{e^{i\theta}}) \geq 0 \qquad \forall \theta \in [0,2\pi]
\end{equation}
is equivalent to positivity of the associated Toeplitz matrix
\begin{equation}
  \label{eq:11a}
  \textrm{Toep}(a_0,\dots,a_m) :=
  \begin{bmatrix}
    a_0 & a_1 & a_2 & \cdots & a_m \\
    a_{-1} & a_0 & a_1 & \ddots & \vdots \\
    a_{-2} & a_{-1} & a_0 & \ddots & \vdots \\
    \vdots & \ddots & \ddots & \ddots & a_1 \\
    a_{-m} & \cdots & \cdots & a_{-1} & a_0
  \end{bmatrix}
\end{equation}
A matrix is \emph{Toeplitz} if the entries are constant along the diagonals.
We collect the information in the following lemma:
\begin{lem}\cite[pg. 17]{MR890515}
  Suppose $f$ is a trigonometric polynomial with data $(a_0,\dots,a_m)$.
  Let $w=
  \begin{bmatrix}
    w_0 & \cdots & w_m
  \end{bmatrix}^T \in \mathbb{C}^{m+1}$.
  Define $w(z) = w_0 + w_1 z + \dots + w_m z^m \in \mathbb{C}[z]$.
  Then:
  \begin{equation}
    \label{eq:12a}
    \frac{1}{2\pi}\int_{-\pi}^\pi |w(e^{i\theta})|^2 f(e^{i\theta},\overline{e^{i\theta}}) d\theta
    =
    w^\ast \textrm{Toep}(a_0,\dots,a_m) w
  \end{equation}
  In particular, if $f(e^{i\theta},\overline{e^{i\theta}}) \geq 0$ for all $\theta \in [0,2\pi]$,
  then $\textrm{Toep}(a_0,\dots,a_m) \geq 0$.
\end{lem}

Now consider a holomorphic polynomial $h(z) = h_0 + h_1z + \dots + h_m z^m \in \mathbb{C}[z]$.
Then
\begin{equation}
  \label{eq:13a}
  |h(z)|^2 =
  \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^m
  \end{bmatrix}^\ast
  \begin{bmatrix}
    h_0 \bar{h}_0 & h_1\bar{h}_0 & \cdots &  h_m \bar{h}_0 \\
    h_0 \bar{h}_1 & h_1 \bar{h}_1 & \cdots & h_m \bar{h}_1 \\
    \vdots & \vdots & \ddots & \vdots \\
    h_0 \bar{h}_m & h_1 \bar{h}_m & \cdots & h_m \bar{h}_m
  \end{bmatrix}
  \begin{bmatrix}
    1 \\ z \\ \vdots \\ z^m
  \end{bmatrix}
\end{equation}
Suppose further that $f$ is a trigonometric polynomial with data $(a_0,\dots,a_m)$.
Then, after reducing coefficients modulo $(z\bar{z}-1)$, we see that the following conditions are equivalent:
\begin{enumerate}
\item $f(z,\bar{z}) \equiv |h(z)|^2 \bmod (z \bar{z}-1)$
\item $a_k = \sum_{j=k}^m h_j \bar{h}_{j-k}$ for $k=0,\dots,m$
\end{enumerate}

From which we can conclude that the following conditions are equivalent:
\begin{enumerate}
\item $f \in \Sigma^2_h + (z\bar{z}-1)$
\item There exist $h_0,\dots,h_m \in \mathbb{C}$ such that $a_k = \sum_{j=k}^m h_j \bar{h}_{j-k}$ for $k=0,\dots,m$
\end{enumerate}

These equations characterize the coefficients of the holomorphic polynomial $h(z)$ (see \cite[pg. 22]{MR890515}), and can be solved by the method of spectral factorization.

Our plan for the case $N>1$ is to develop block analogues of these ideas and invoke the operator-valued Riesz-Fejer lemma.

\subsection{Block Trace Parametrization}\label{sec:block-trace-param}

In this section we discuss a block analog of the trace parametrization technique utilized throughout \cite{MR3618747}.

Now consider a Hermitian polynomial $g(z,\bar{z})$.
We may write $g$ in block matrix form:
\begin{equation}
  \label{eq:14a}
  g(z,\bar{z}) =
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}^\ast
  \begin{bmatrix}
    B_{00} & B_{01} & \cdots & B_{0m} \\
    B_{10} & B_{11} & \cdots & B_{1m} \\
    \vdots & \vdots & \ddots & \vdots \\
    B_{m0} & B_{m1} & \cdots & B_{mm}
  \end{bmatrix}
    \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}
\end{equation}
where $B_{jk} \in \mathbb{C}^{N \times N}$ and $B_{kj} = B_{jk}^\ast$.
By reducing the coefficients modulo $(z^N \bar{z}^N - 1)$ we see that $g$ is congruent modulo $(z^N \bar{z}^N-1)$ to a unique polynomial of the form

\begin{equation}
  \label{eq:7d}
    f(z,\bar{z}) =
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}^\ast
  \begin{bmatrix}
   A_0 & A_1 & \cdots & A_m \\
   A_{-1} & 0 &  \cdots & 0 \\
   \vdots & \vdots & \ddots & \vdots \\
   A_{-m} & 0 & \cdots & 0
  \end{bmatrix}
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}
\end{equation}

where $A_j \in \mathbb{C}^{N \times N}$, $A_{-j} = A_j^\ast$ and $\psi(z) = \begin{bmatrix} 1 & z & \cdots & z^{N-1} \end{bmatrix}^T$.
 If a Hermitian polynomial $f$ has the above form, then we say $f$ is \emph{trigonometric modulo $(z^N\bar{z}^N-1)$ of degree $m$ with data $(A_0,\dots,A_m)$}.

 Since
 \begin{equation}
   \label{eq:6b}
   g(z,\bar{z}) = f(z,\bar{z}) + q(z,\bar{z})(z^N \bar{z}^N -1)
 \end{equation}
for some Hermitian polynomial $q$, we see that
 \begin{equation}
   \label{eq:15a}
   \textrm{Gram}(f)[\xi,\omega \xi, \dots, \omega^{N-1} \xi] = \textrm{Gram}(g)[\xi,\omega \xi, \dots, \omega^{N-1}\xi]
 \end{equation}
 Hence the Gram matrices are simultaneously positive.

 Furthermore, we have that $f \in \Sigma^2_h + (z^N \bar{z}^N -1)$ if and only if $g \in \Sigma^2_h + (z^N\bar{z}^N-1)$.

 If $f$ is trigonometric modulo $(z^N \bar{z}^N -1)$ with data $(A_0,\dots,A_m)$, then our goal is to show that the condition
 \begin{equation}
   \label{eq:16}
   \textrm{Gram}(f)[\xi,\omega \xi,\dots,\omega^{N-1}\xi] \geq 0 \qquad \forall \xi \in \mathbb{T}
 \end{equation}
implies positivity of the associated block Toeplitz matrix
\begin{equation}
  \label{eq:17a}
   \textrm{Toep}(A_0,\dots,A_m) :=
  \begin{bmatrix}
    A_0 & A_1 & A_2 & \cdots & A_m \\
    A_{-1} & A_0 & A_1 & \ddots & \vdots \\
    A_{-2} & A_{-1} & A_0 & \ddots & \vdots \\
    \vdots & \ddots & \ddots & \ddots & A_1 \\
    A_{-m} & \cdots & \cdots & A_{-1} & A_0
  \end{bmatrix}
\end{equation}
which allows us to apply the operator-valued RF theorem to obtain $f \in \Sigma^2_h + (z^N \bar{z}^N -1)$.

The idea of the trace parametrization depends on the following observation:
Let $T_k$ denote the elementary block Toeplitz matrix with $I$ on the $k$-th diagonal and $0$ elsewhere, where the main diagonal is counted as $k=0$ and positive diagonals count to the right.
Let $\textrm{Trace}[Q]$ denote the sum of the diagonal blocks of $Q=[Q_{jk}]_{j,k=0}^m$.
Then
\begin{equation}
  \label{eq:21a}
  \textrm{Trace}[T_{-k}Q] = \sum_{j=k}^m Q_{j-k,j}
\end{equation}
is the sum of the $k$-th diagonal of $Q$.

Suppose $h(z)$ is a holomorphic polynomial of degree $N(m+1)$.
We can write $|h(z)|^2$ in block matrix form:
\begin{equation}
  \label{eq:18a}
  |h(z)|^2 =
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}^\ast
  \begin{bmatrix}
    Q_{00} & Q_{01} & \cdots & Q_{0m} \\
    Q_{10} & Q_{11} & \cdots & Q_{1m} \\
    \vdots & \vdots & \ddots & \vdots \\
    Q_{m0} & Q_{m1} & \cdots & Q_{mm}
  \end{bmatrix}
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}
\end{equation}
where $Q = [Q_{jk}]_{j,k=0}^m$ is a positive Hermitian block matrix with $Q_{jk} \in \mathbb{C}^{N \times N}$ and $Q_{kj} = Q_{jk}^\ast$.

Suppose further that $f$ is trigonometric modulo $(z^N \bar{z}^N -1)$ with data $(A_0,\dots,A_m)$.
Then, after by considering the reduction of coefficients modulo $(z^N\bar{z}^N-1)$, we see that the following conditions are equivalent:
\begin{enumerate}
\item $f(z,\bar{z}) \equiv |h(z)|^2 \bmod (z^N \bar{z}^N -1)$
\item $A_k = \textrm{Trace}[T_{-k} Q]$
\end{enumerate}

From which we conclude that the following conditions are equivalent:
\begin{enumerate}
\item $f \in \Sigma^2_h + (z^N \bar{z}^N -1)$
\item There exists a positive block matrix $Q=[Q_{jk}]_{j,k=0}^m$ such that $A_k = \textrm{Trace}[T_{-k}Q]$, $k=0,\dots,m$.
\end{enumerate}

\subsection{Operator-Valued RF Theorem}

We recall the operator-valued Riesz-Fejer theorem:
\begin{thm}\cite[Theorem 2.1]{MR2743422}
  Let $A(z) = \sum_{k=-m}^m A_k z^k$ be a Laurent polynomial with matrix coefficients $A_k \in \mathbb{C}^{N \times N}$.
  The following conditions are equivalent:
  \begin{enumerate}
  \item $A(\xi)\geq 0$ for all $\xi \in \mathbb{T}$
  \item $\textrm{Toep}(A_0,\dots,A_m) \geq 0$
  \item There exists $P(z) = P_0 + P_1 z + \dots + P_m z^M$ with matrix coefficients $P_k \in \mathbb{C}^{N \times N}$ such that $A(z) = P(z)^\ast P(z)$
  \end{enumerate}
\end{thm}

Let $Q_{jk} = P_j^\ast P_k$ and let $Q= [ Q_{jk} ]_{j,k=0}^m$.
If $A(z) = P(z)^\ast P(z)$, then we can equate coefficients to get $A_k = \textrm{Trace}[T_{-k}Q]$.
Thus we get:

\begin{cor}\label{sec:operator-valued-rf:OPERATOR-RF} Let $A_0,\dots,A_m \in \mathbb{C}^{N \times N}$.
  The following conditions are equivalent:
  \begin{enumerate}
  \item $\textrm{Toep}(A_0,\dots,A_m) \geq 0$
  \item There exists a positive block matrix $Q = [Q_{jk}]_{j,k=0}^m$ such that $A_k = \textrm{Trace}[T_{-k}Q]$.
  \end{enumerate}
\end{cor}

With these considerations in mind, starting with $f$ trigonometric modulo $(z^N \bar{z}^N -1)$ with data $(A_0,\dots,A_m)$, we will assume
\begin{equation}
  \label{eq:19a}
  \textrm{Gram}(f)[\xi, \omega \xi, \dots, \omega^{N-1} \xi] \geq 0 \qquad \forall \xi \in \mathbb{T}
\end{equation}
and show that this implies
\begin{equation}
  \label{eq:20a}
  \textrm{Toep}(A_0,\dots,A_m) \geq 0
\end{equation}
from which we can obtain $f \in \Sigma^2_h + (z^N \bar{z}^N -1)$.

\section{Proof of Main Results}

  We are now in position to utilize all the tools developed so far and show that the matrix positivity conditions imply positivity of the associated block Toeplitz matrix.
  We then prove the theorem characterizing Hermitian sums of squares modulo $(z^N \bar{z}^N-1)$.

\begin{prop}\label{sec:repr-matr-prod:gram-to-toep}
 Let $\omega = e^{2 \pi i / N}$.

 Let $f \in \mathbb{C}_h[z,\bar{z}]$ be trigonometric mod $(z^N \bar{z}^N - 1)$ with data $(A_0,\dots,A_m)$.

 Suppose $\text{Gram}(f)[e^{i\theta},\omega e^{i\theta},\dots,\omega^{N-1}e^{i\theta}] \geq 0$ for all $\theta \in [0,2\pi]$.

 Then $\text{Toep}(A_0,\dots,A_m) \geq 0$.
\end{prop}
\begin{proof}

Let $v^{(0)},\dots,v^{(m)} \in \mathbb{C}^N$ be arbitrary.

Write $v^{(j)} = \begin{bmatrix} v_0^{(j)} & \cdots & v_{N-1}^{(j)} \end{bmatrix}^T$ with $v_k^{(j)} \in \mathbb{C}$.

Let $v = \begin{bmatrix} v^{(0)} & \cdots & v^{(m)} \end{bmatrix}^T \in (\mathbb{C}^N)^{m+1}$.

It suffices to show $v^\ast \text{Toep}(A_0,\dots,A_m) v \geq 0$.

Observe that
\[
 v^\ast \text{Toep}(A_0,\dots,A_m) v = \sum_{j,k=0}^m v^{(j)\ast} A_{k-j} v^{(k)}
\]

For $j=0,\dots,m$ define $\tilde{v}^{(j)}(z) = \sum_{k=0}^{N-1} v_k^{(j)} z^{N-k} \in \mathbb{C}[z]$.

Define $v(z) = \sum_{j=0}^m z^{Nj} \tilde{v}^{(j)}(z) \in \mathbb{C}[z]$.

Let $\psi(z) = \phi_N(z) = \begin{bmatrix} 1 & z & \cdots & z^{N-1} \end{bmatrix}^T$.

Claim: $\mathcal{F}_N( |v(z)|^2 f(z,\bar{z}) ) = v^\ast \text{Toep}(A_0,\dots,A_m) v \geq 0$.

Introduce double-index notation for matrix coefficients of $f$ in order to simplify computations:
 \[
  f(z,\bar{z}) =
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}^\ast
  \begin{bmatrix}
    A_{00} & A_{01} & \cdots & A_{0m} \\
    A_{10} & A_{11} & \cdots & A_{1m} \\
    \vdots & \vdots & \ddots & \vdots \\
    A_{m0} & A_{m1} & \cdots & A_{mm}
  \end{bmatrix}
  \begin{bmatrix} \psi(z) \\ z^N \psi(z) \\ \vdots \\ z^{Nm}\psi(z) \end{bmatrix}
 \]
where $A_{jk} = 0$ if $\min(k,j) \neq 0$ and $A_{jk} = A_{k-j}$ otherwise.

Then:
\[
 f(z,\bar{z}) = \sum_{j,k=0}^m \bar{z}^{jN} z^{kN} \psi(z)^\ast A_{jk} \psi(z)
\]

We have:
\[
 |v(z)|^2 = \sum_{j,k=0}^m \bar{z}^{Nj} z^{Nk} \overline{ \tilde{v}^{(j)}(z) } \tilde{v}^{(k)}(z)
\]
Hence:
\[
 |v(z)|^2 f(z,\bar{z})
 =
 \sum_{j,k,c,d=0}^m \bar{z}^{N(c+j)} z^{N(d+k)} \overline{ \tilde{v}^{(j)}(z) } \tilde{v}^{(k)}(z) \psi(z)^\ast A_{cd} \psi(z)
\]
Then using Proposition \ref{sec:defin-basic-prop-1:computation_of_F} and passing between the double and single index coefficients for $f$:
\begin{align}
\mathcal{F}_N( |v(z)|^2 f(z,\bar{z}) ) &= \sum_{j,k,c,d=0}^m \mathcal{F}_N ( \bar{z}^{N(c+j)} z^{N(d+k)} \overline{ \tilde{v}^{(j)}(z) } \tilde{v}^{(k)}(z) \psi(z)^\ast A_{cd} \psi(z) ) \\
&= \sum_{j,k,c,d=0}^m \delta_{c+j,d+k} v^{(j)\ast} A_{cd} v^{(k)} \\
&= \sum_{\substack{j,k,c,d=0 \\ c+j=d+k}}^m v^{(j)\ast} A_{cd} v^{(k)} \\
&= \sum_{\substack{j,k,c,d=0 \\ j-k=d-c}}^m v^{(j)\ast} A_{d-c} v^{(k)} \\
&= \sum_{j,k=0}^m v^{(j)\ast} A_{j-k} v^{(k)} \\
&= v^\ast \text{Toep}(A_0,\dots,A_m)^T v
\end{align}

By Proposition \ref{sec:defin-basic-prop:functional_positivity}  we have
\[
 v^\ast \text{Toep}(A_0,\dots,A_m)^T v \geq 0 \qquad \forall v \in (\mathbb{C}^N)^{m+1}
\]
Therefore $\text{Toep}(A_0,\dots,A_m)^T \geq 0$, and hence $\text{Toep}(A_0,\dots,A_m) \geq 0$.

\end{proof}

\begin{thm}\label{MAIN-THEOREm}
 Let $\omega = e^{2 \pi i / N}$.

 Let $f \in \mathbb{C}_h[z,\bar{z}]$ be trigonometric mod $(z^N \bar{z}^N - 1)$ with data $(A_0,\dots,A_m)$.

 The following conditions are equivalent:
 \begin{enumerate}
  \item $f \in \Sigma^2_h + (z^N \bar{z}^N - 1)$
  \item $\text{Gram}(f)[e^{i\theta},\omega e^{i\theta},\dots,\omega^{N-1}e^{i\theta}] \geq 0$ for all $\theta \in [0,2\pi]$
  \item $\text{Toep}(A_0,\dots,A_m) \geq 0$
 \end{enumerate}
\end{thm}
\begin{proof}

( (i) $\implies$ (ii) ) By Lemma \ref{MATRIX-POSITIVITY-CONDITIONS}.

( (ii) $\implies$ (iii) ) By Proposition \ref{sec:repr-matr-prod:gram-to-toep}.

( (iii) $\implies$ (i) ) By Corollary \ref{sec:operator-valued-rf:OPERATOR-RF} and the conclusion of Section \ref{sec:block-trace-param}.

\end{proof}