\section{Introduction}\label{sec:intro}
DAVs exhibit multi-mode oscillations.\footnote{Photometry dominates the discovery and investigation of DAVs.} The longest period modes have $P\sim 2\,\mathrm{min}$ near the blue edge and $P\sim 20\,\mathrm{min}$ near the red edge of the DAV instability strip. Longer-period modes in cooler DAVs exhibit larger photometric variations and greater amplitude and frequency variability than shorter period modes in hotter ones \citep{Clemens-1993, Mukadam}.

Brickhill \citep{Brickhill-1983, Brickhill-1990, Brickhill-1991a, Brickhill-1991b} recognized that the convective flux in DAVs responds on a much shorter timescale than the oscillation periods of their overstable g-modes. He showed that g-modes with $\omega\gtrsim \tau_c^{-1}$ are overstable due to a mechanism he named `convective driving'. Here $\tau_c$ is the effective thermal timescale of the convection zone (CVZ).\footnote{Our paper adopts terminology used by \cite{Goldreich-Wu-I} in their investigation of convective driving.} Brickhill also realized that the amplitude of the flux perturbation that emerges from the photosphere, $\Delta F_{ph}$, is diminished by the factor $(1+(\omega\tau_c)^2)^{-1/2}$ relative to the one that enters the bottom of the CVZ, $\Delta F_b$. His twin insights explain why overstable modes shift to longer periods as a DAV cools and why longer period modes usually show larger photometric variations. These trends are well-established among DAVs above the red edge.

By itself,  the overstability criterion, $\omega\tau_c>1$ cannot account for the red edge of the DAV instability strip.  An essential step toward understanding its origin was taken by \cite{Dolez-Vauclair}  in a paper that predated Brickhill's work on convective driving. Their linear nonadiabatic calculations found $\omega>L_{\ell,b}$ to be a necessary condition for the overstability of g-modes in cool DAVs.\footnote{$L_{\ell,b}$ is the Lamb frequency at the bottom of the surface CVZ for modes of angular degree $\ell$.} However, \cite{Dolez-Vauclair} mistakenly attributed overstability to the $\kappa$-mechanism. Later, \cite{Pesnell} noted that the overstability they found was caused by `convective blocking' due to their incorrect assumption of frozen convection.\footnote{That the assumption of frozen convection could lead to mode excitation was well-established prior to its being named convective blocking \citep{Cox1}.}  Convective blocking is the antithesis of convective driving. It is the latter which is responsible for the excitation of overstable g-modes across the entire DAV instabililty strip. The identification by \cite{Dolez-Vauclair} of $\omega>L_{\ell,b}$ as a necessary condition for overstability remains valid under convective driving as presented in the numerical results of \citep{Wu-Goldreich-II}. Thus together, $\omega>1/\tau_c$ and $\omega>L_{\ell,b}$ define the low frequency limit, $\omega_\mathrm{crit}$, for overstable g-modes at each $T_{\mathrm{eff}}$.

Inspired by \cite{Dolez-Vauclair}, Brickhill appears to have accepted $\omega> L_{\ell,b}$ as a necessary condition for overstability. However, we are unsure why.  We have been unable to locate a statement in any of Brickhill's papers on DAVs in which he clearly attributes this condition to radiative damping exceeding convective excitation.  Instead, his focus seems to have been on turbulent damping in the convective overshoot region below $z_b$ \citep{Brickhill-1991b}. This is a puzzle, since the calculations by \cite{Dolez-Vauclair} only accounted for radiative damping.

In our paper, we explain that modes with $\omega<L_{\ell,b}$ suffer enhanced radiative damping that  overwhelms convective driving. We demonstrate that the quasi-adiabatic approximation suffices to account for the enhancement in radiative damping and thus provides a minimally complicated conceptual understanding of the red edge, although it would be ideal to have fully non-adiabatic calculations including turbulent viscosity.

As a DAV cools, $1/\tau_c$ steeply declines and $L_{\ell,b}$ gradually rises; $\omega_\mathrm{crit}$ reaches a minimum value close to the red edge as defined by photometric surveys of DAVs. Well below the red edge, photospheric flux variations associated with overstable modes are severely diminished; at $\omega_\mathrm{crit}$, $\Delta F_{ph}/\Delta F_b \sim [(1+(L_{\ell,b}\tau_c)^2]^{-1}\ll 1$.  Moreover, almost nothing is known about DAVs below the red edge because photometry is a poor technique for identifying them.  Although more challenging, spectroscopic observations of bright DAs just below the red edge might succeed in detecting doppler shifts due to the largely horizontal velocities associated with overstable modes. These do not suffer significant reductions between $z_b$ and $z_{ph}$ \citep{van-Kerkwijk}.

The discovery of outbursts in photometric data obtained by the {\it Kepler} satellite from cool DAVs is an exciting recent development \citep{Bell-2015, Hermes, Bell-2016a, Bell-2016b, Hermes-2017}. Individual outbursts recur sporadically with separations of days and emit $10^{33}-10^{34}\,\mathrm{erg}$ on timescales of hours.  \cite{Hermes} suggested that outbursts may involve nonlinear mode couplings as described by \cite{Wu-Goldreich-IV}. To lowest order, relevant nonlinearities couple an overstable parent mode to a pair of lower frequency stable daughter modes. Photometric amplitudes observed across the DAV instability strip are generally in line with those predicted by calculations of saturation based on 3-mode couplings \citep{Wu-Goldreich-IV}.

We concur with \cite{Hermes} and associate outbursts with limit cycles exhibited by 3-mode couplings satisfying $\delta\omega<|\gamma_d|$.  Here $\delta\omega\equiv |\omega_p-\omega_{d_1}-\omega_{d_2}|$ is the frequency mismatch between an overstable parent mode (labeled by $p$) and two stable daughter modes (labeled by $d_1$ and $d_2$), and $\gamma_d\equiv (\gamma_{d_1}+\gamma_{d_2})$ is the average damping rate of daughter modes \citep{Wu-Goldreich-IV}. Consistent with observations, this condition is most likely to be satisfied in cool DAVs near the red edge of the instability strip.  The red edge is where the lowest frequency overstable parent modes are found, where the density of these modes per unit frequency is greatest, and where candidate daughter modes have the largest damping rates, $\gamma_d\sim -(n_d\tau_{th,b})^{-1}$.\footnote{Here $n_d$ is the number of radial nodes of the daughter modes and $\tau_{th,b}\sim \tau_c/5$ is the thermal timescale at the bottom of the CVZ, $z_b$.}  Away from the red edge, toward both higher and lower $T_{\mathrm{eff}}$, the lowest frequency for an overstable mode, $\omega_\mathrm{crit}$, increases, and $\delta\omega$ for candidate 3-mode couplings also increases accordingly because density of g-modes per unit frequency grows with higher frequency. Consequently, $\delta\omega<|\gamma_d|$, is most easily achieved near the red edge.

Our paper is arranged as follows. Section~\ref{sec:red-edge} recapitulates the overstability criterion, $\omega>\max(\tau_c^{-1},L_{\ell,b})$, focusing on the physical explanation for the second threshold $\omega>L_{\ell,b}$ which essentially yields the photometric red edge. Section~\ref{sec:burst} argues that parametric instability involving 3-mode couplings that are sufficiently resonant to trigger limit cycles can quantitatively account for the outbursts. We summarize and discuss unresolved problems in Section~\ref{sec:conclusion}. All figures except Figure~\refnew{fig:omegatauc-T} are based on a WD model with $T_{\mathrm{eff}}\approx 10860\,\,\mathrm{K}$ and $g\approx 10^8\,\mathrm{cm\, s^{-2}}$ that is described in Section~\ref{sec:rededge}.

\section{Red edge}\label{sec:red-edge}

\subsection{Order-of-Magnitude Results}

Gravity waves propagate where $\omega$ is smaller than both the Brunt-Vaisala frequency, $N$, and the Lamb frequency, $L_{\ell}$ \citep[e.g.][]{Cox}. We denote by $z_\omega\sim \omega^2 r^2/(g \ell(\ell+1))$ the depth at which $\omega=L_{\ell}$, and by $z_b$ the depth at the bottom of the surface CVZ. The upper lid of a g-mode's propagation cavity is detached from the CVZ if $z_\omega>z_b$, or equivalently, $\omega>L_{\ell,b}\equiv L_{\ell}(z_b)$. Otherwise the propagation cavity terminates at $z_b$.

Consider a g-mode with $z_\omega\gg z_b$.  The mode is evanescent in both the CVZ and the first scale height below $z_b$, and each stores a similar amount of mode energy. Upon compression of the evanescent region, the first scale height, $H_1$, loses heat (entropy) by radiative diffusion, and thus contributes to mode damping, $\gamma_\mathrm{rad}$. Note that $H_1$ dominates $\gamma_\mathrm{rad}$, because the local damping rate is inversely proportional to the local thermal timescale $\tau_\mathrm{th}(z)\sim p H_p/F\appropto z^6$ which grows steeply with depth. Meanwhile entropy lost by the radiative zone gets absorbed into the CVZ and does not escape until after a timescale $\sim \tau_c$. Thus the CVZ contributes to mode excitation.  According to \cite{Goldreich-Wu-I}, the ratio of convective driving to radiative damping is
\begin{equation}
{\gamma_{\mathrm{cvz}}\over\left|\gamma_{\mathrm{rad}}\right|}\sim{2(\omega\tau_c)^2\over 1+(\omega\tau_c)^2}\, ,\,\, (\omega\gg L_{\ell,b})\, .\label{eq:drive-damp-ratio}
\end{equation}

Next we consider a g-mode with $z_\omega\ll z_b$, or equivalently $\omega\ll L_{\ell,b}$. This mode propagates up to $z_b$, and its effective wavelength, $\lambda_1$, becomes smaller than $H_1$. Radiative diffusion rate per unit volume is proportional to the Laplacian of the temperature perturbation associated with the mode. Therefore $\gamma_\mathrm{rad}$ gets enhanced by a factor $(H_1/\lambda_1)^2 \sim (k_z H_1)^2 \sim (k_h N H_1/\omega)^2\sim (k_h c_s/\omega)^2 \sim (L_{\ell,b}/\omega)^2$.\footnote{We substitute Equation (18) into the second line of Equation (50) in \cite{Goldreich-Wu-I} to get this result. For the purpose of order-of-magnitude estimate, we only retain terms on the order of $k_z^2$ in $\gamma_\mathrm{rad}$ and neglect terms on the order of $k_z /H_1$ and $H_1^{-2}$. We keep all the terms in the quasi-adiabatic calculations described in the subsection \ref{sec:quasi-adiabatic}.} Consequently, radiative damping is enhanced by a factor $(L_{\ell,b}/\omega)^2$ relative to convective driving and
\begin{equation}
{\gamma_{\mathrm{cvz}}\over\left|\gamma_{\mathrm{rad}}\right|}\sim  {2(\omega\tau_c)^2\over 1+(\omega\tau_c)^2}\left(\omega\over L_{\ell,b}\right)^2\, ,\,\, (\omega\ll L_{\ell,b})\, ,\label{eq:drive-damp-ratio-smallw}
\end{equation}

Hot DAVs have shallow CVZs such that $L_{\ell,b} \ll 1/\tau_c$ and modes with $\omega\ll L_{\ell,b}$ are damped according to Equation~\refnew{eq:drive-damp-ratio-smallw}. Equation~\refnew{eq:drive-damp-ratio} sets the boundary between overstable and damped modes at
\begin{equation}
\omega_\mathrm{crit} \sim {1\over \tau_c}\, , \,\, (L_{\ell,b}\ll 1/\tau_c)\, .\label{eq:omegacrit-hot}
\end{equation}
Cool DAVs have deep CVZs such that $L_{\ell,b}\gg 1/\tau_c$ and modes with $\omega\gg L_{\ell,b}$ are overstable according to Equation~\refnew{eq:drive-damp-ratio}.  Equation~\refnew{eq:drive-damp-ratio-smallw} sets the boundary between overstable and damped modes at
\begin{equation}
\omega_\mathrm{crit}\sim {L_{\ell,b}\over \sqrt{2}}\, ,\,\, (L_{\ell,b}\gg 1/\tau_c)\, .\label{eq:omegacrit-cool}
\end{equation}

The flux perturbation that emerges from the photosphere is diminished relative to that entering the base of the CVZ following the relation
\begin{equation}\label{eq:dilution}
{\Delta F_{\mathrm{ph}}\over F}\approx {1\over 1-i\omega\tau_c}\left({\Delta F\over F}\right)_{z_b}\, .
\end{equation}
In cool WDs, overstable modes with $\omega\gtrsim \omega_\mathrm{crit}$ all have $\omega\tau_c\gg 1$. Thus they escape detection in photometric surveys. This explains the red edge of the DAV instability strip. It reflects the transition between the two overstability criterions which results from the deepening of the surface CVZ as a WD cools. Existence of overstable modes at $T_{\mathrm{eff}}$ well below the photometric red edge is the key to test our explanation of the red edge. High-quality spectroscopy of a nearby cool WD may reveal velocities associated with an overstable g-mode since velocities do not suffer the reduction effect like the flux perturbation, although such test seems very challenging in practice.\footnote{Based on private communications with Chris Clemens and JJ Hermes.}

\subsection{Comparison with Quasi-Adiabatic Calculations}\label{sec:quasi-adiabatic}
{\it MESA} \citep{Paxton} is employed to produce WD models with total mass, $M\approx 0.60M_\odot$, surface gravity, $g\approx 10^8\,\mathrm{cm\, s^{-2}}$, and hydrogen envelope mass, $M_H\sim 10^{-4}M_\mathrm{WD}$. We adopt the ML2 version of the mixing length prescription, and adjust the mixing length, $\alpha$, such that entropy levels at the bottom of the CVZ agree with those in Figure~(1) of \cite{Tremblay-MLT}.  Our WD models are input to {\it gyre} \citep{Townsend} to calculate adiabatic eigen-frequencies and eigen-functions for g-modes. The quasi-adiabatic technique as described in \cite{Goldreich-Wu-I} is used to calculate convective driving and radiative damping rates. We adopt the full expression for the quasi-adiabatic flux perturbation in Equation (10) of \cite{Goldreich-Wu-I}, and then we apply their equation (50) to calculate the driving and damping rates. Here we focus on comparing the order-of-magnitude results with those obtained from quasi-adiabatic calculations.

Figure~\refnew{fig:omegacrit-Teff} shows that angular frequencies of marginally overstable g-modes with $\ell=1$ as a function of $T_{\mathrm{eff}}$ exhibit a gradual transition between the two limiting cases, Equations~\refnew{eq:omegacrit-hot} and \refnew{eq:omegacrit-cool}.
\begin{figure}
\includegraphics[width=0.99\linewidth]{omegacrit-Teff.pdf}
\caption{\label{fig:omegacrit-Teff}  The black dots show the radian frequency for marginally overstable modes, $\omega_\mathrm{crit}$, versus $T_{\mathrm{eff}}$, based on quasi-adiabatic calculations. The solid curve shows $1/\tau_c$, the overstability criterion suitable for hot DAVs. The dashed one shows $L_{\ell,b}/\sqrt{2}$, the overstability criterion suitable for cool WDs.    }
\end{figure}

Figure~\refnew{fig:omegatauc-T} plots the reduction factor for the photometric variation, $\omega_\mathrm{crit}\tau_c$. It is on the order of a few above the observational red edge, and increases exponentially with lower $T_{\mathrm{eff}}$ below the red edge. The dashed line show our fitting formula,
\begin{equation}\label{eq:omegatauc-Teff-fit}
\omega_\mathrm{crit}\tau_c\sim \exp\left[(11012\,\mathrm{K} - T_{\mathrm{eff}})\over 111\,\mathrm{K}\right]\, ,\, (\ell=1)\, .
\end{equation}
Therefore, overstable modes become invisible below the red edge.
\begin{figure}
\includegraphics[width=0.99\linewidth]{omegatauc-T.pdf}
\caption{\label{fig:omegatauc-T}  The reduction factor for the photometric flux perturbation for marginally overstable modes, $\omega_\mathrm{crit}\tau_c$, versus $T_{\mathrm{eff}}$. The black dots are based on our quasi-adiabatic calculations using WD models made by MESA as detailed in Section~\ref{sec:rededge}. The dashed line plots the fitting formula, Equation~\refnew{eq:omegatauc-Teff-fit}. The observational photometric red edge is $\sim 10800\,\mathrm{K}$ where $\omega_\mathrm{crit}\tau_c\sim 10$. The reduction factor increases exponentially with lower temperature, making overstable modes invisible below the red edge.
\end{figure}

$T_{\mathrm{eff}}$ at the red edge depends sensitively upon $z_b$, the depth the surface CVZ. But $z_b$ has a sensitive dependence on the mixing length parameter, $\alpha$.  This warrents a short detour to explain how $\alpha$ affects $z_b$. Except for the top scale height right below the photosphere, the bulk of the CVZ is nearly isentropic, i.e., $ds/dz\sim 0$, so the specific entropy is well-approximated by $s_b$. Given $T_{\mathrm{eff}}$ and $g$ at the photosphere, a smaller $\alpha$ corresponds to a smaller convection efficiency and thus yields a smaller $ds/dz$ in the top scale height, and therefore a lower $s_b$. At fixed density, smaller entropy corresponds to lower temperature and higher opacity, and thus a deeper surface CVZ. WD models at the same $T_{\mathrm{eff}}$ but with different choices of $\alpha$ could have completely different depths of CVZs and thus would predict very different $\omega_\mathrm{crit}$. The $\alpha$ in our WD models at different $T_{\mathrm{eff}}$ are chosen so that the $s_b$ agrees with those in \cite{Tremblay-MLT} who calibrated $s_b$ by observational spectra of WDs.

\subsection{At the Red Edge}\label{sec:rededge}
Material in this subsection is based on a particular DA WD model at $T_{\mathrm{eff}}\approx 10860\,\mathrm{K}$ and $g\approx 10^8\,\mathrm{cm\, s^{-2}}$ with the longest period overstable dipole modes at $P\approx 20\,\mathrm{min}$. At the base of the CVZ, the specific entropy and pressure are $s_b\sim 2.6\times 10^9\,\mathrm{erg\, g^{-1}\, K^{-1}}$ and $p_b\approx 7.5\times 10^8\,\mathrm{dyn\, cm^{-2}}$, the conventional thermal timescale, $\tau_\mathrm{th,b}\equiv \int^{z_b} dz {c_p \rho T/F}\approx 5.13\,\mathrm{min}$ as compared to the effective thermal timescale, $\tau_c\approx 15.6\,\mathrm{min}$. Radiative damping is concentrated within the first scale height below the CVZ, $z_b<z<z_b+H_1$. At the red edge, the lowest frequency of overstable modes satisfies $\omega_\mathrm{crit}\tau_\mathrm{th}\approx 1.5$ at $z_b$, and $\omega_\mathrm{crit}\tau_\mathrm{th}\sim 7.2$ at $z_b+H_1$. Therefore, the quasi-adiabatic condition, $\omega\tau_\mathrm{th}\gg 1$, is marginal for overstable modes at the red edge. Because $\omega_\mathrm{crit}\tau_c$ grows exponentially with lower $T_{\mathrm{eff}}$ below the red edge and $\tau_c\sim 5\tau_\mathrm{th}(z_b)$, the quasi-adiabatic condition is well-satisfied for overstable modes in DAVs that are cooler than those at the red edge. Ideally, non-adiabatic calculations that include turbulent viscosity should be carried out to determine the theoretical red edge temperature, $T_\mathrm{red}$. However, as described in the previous subsection, the choice of the mixing length $\alpha$ already inserts considerable uncertainty in $T_\mathrm{red}$. Quasi-adiabatic approximations have the virtue of technical simplicity while being adequate to explain the existence of the red edge.

Work integrals for an overstable dipole mode with $z_\omega\ggz_b$ and a damped diple mode with $z_\omega\ll z_b$ are displayed in Figure~\refnew{fig:WI-compare}.
\begin{figure}
\includegraphics[width=0.9\linewidth]{WI-over-stable.eps}
\includegraphics[width=0.9\linewidth]{WI-stable.eps}
\caption{\label{fig:WI-compare} Work integrals for two $\ell=1$ g-modes in our MESA WD model with $T_{\mathrm{eff}}\approx 10860\,\mathrm{K}$. The upper panel and lower panels demonstrate that the mode with $z_\omega\ggz_b$ is overstable, whereas the mode with $z_\omega\ll z_b$ is damped. Note that $\gamma$ at pressure $p$ is the net contribution to damping and driving from the region between the stellar center up to pressure $p$. The total damping or driving rate is given by the left end of the black solid curve in each panel. Note that quasi-adiabatic approximation for the mode shown in the lower panel is marginal.}
\end{figure}
Results from quasi-adiabatic calculations of linear growth and damping rates for g-modes with $\ell=1,\,\, 2,\,\, 3$ are summarized in Figure~\refnew{fig:timescale}. We will take use of these results for calculating the 3-mode couplings in the next part of the paper.
\begin{figure}
\includegraphics[width=0.99\linewidth]{timescale-1.pdf}
\includegraphics[width=0.99\linewidth]{timescale-2.pdf}
\includegraphics[width=0.99\linewidth]{timescale-3.pdf}
\caption{\label{fig:timescale} Linear growth and damping rates for g-modes with $\ell=1,\,\, 2,\,\, 3$. Solid curves correspond to overstable modes with $\gamma>0$, whereas dashed curves apply to damped modes with $\gamma<0$. Dotted lines indicate $\omega/(\pi n)$, the reciprocal of the time a g-wave takes to travel across its cavity.  A mode transitions to a traveling wave if its damping rate is greater than $\omega/(\pi n)$.
In the same WD model at $T_{\mathrm{eff}}\approx 10860\,\mathrm{K}$, the lowest frequency of an overstable mode rises with increasing $\ell$. Corresponding values of $\omega_\mathrm{crit}\tau_c$ are approximately, $5$ for $\ell =1$, $6$ for $\ell =2$, and $7$ for $\ell =3$. Peak driving rates and associated frequencies, $(\gamma_\mathrm{pk},\nu_\mathrm{pk}, \omega_\mathrm{pk}\tau_c)$, are $(1.4\times 10^{-6}\,\,\mathrm{s^{-1}}, 1048\,\,\mu\mathrm{Hz}, 6.2)$ at $\ell=1$, $(5.1\times 10^{-6}\,\,\mathrm{s^{-1}}, 1376\,\,\mu\mathrm{Hz}, 8.1)$ at $\ell=2$, and $(1.1\times 10^{-5}\,\mathrm{s^{-1}}, 1604\,\mu\mathrm{Hz}, 9.4)$ at $\ell=3$. Horizontal grid lines label $1/(\mathrm{day})$ and $1/(\mathrm{hour})$. Dashed vertical lines mark the overstable mode with maximal growth rate, $\gamma_\mathrm{pk}$. }
\end{figure}

\subsection{Comparison with Previous Work}\label{sec:sub-comparison}
As described in Section \ref{sec:intro}, nonadiabatic calculations by \cite{Dolez-Vauclair} first identified $\omega> L_{\ell,b}$ as a necessary condition for the overstability of g-modes in DAVs. Not long after, \cite{Brickhill-1983}, realized that \cite{Dolez-Vauclair} erred when they attributed mode excitation to the $\kappa$-mechanism. Because the response time of convection in DAVs is shorter than the periods of overstable modes, he proposed a new excitation mechanism that he named `convective driving'.  Brickhill seems to have accepted $\omega> L_{\ell,b}$ as a necessary condition for mode overstability although it is not obvious to us why he did so. In subsequent papers, Brickhill investigated damping by turbulent viscosity in the region of convective overshoot below $z_b$ but to our knowledge, he never clearly implicated radiative damping as the culprit.

A quasi-adiabatic stability analysis by \cite{Goldreich-Wu-I} confirmed Brickhill's results regarding convective excitation and radiative damping for modes with $\omega\ggL_{\ell,b}$.  More significantly, non-adiabatic
stability calculations by \cite{Wu-Goldreich-II} that include convective driving obtained a similar behavior for $\omega_\mathrm{crit}$ as a function of $T_{\mathrm{eff}}$ to that observed by \cite{Dolez-Vauclair}. The results of the current investigation resemble those obtained by \cite{Wu-Goldreich-II} in the sense that as a DAV cools $\omega_\mathrm{crit}$ shows a steep decline above the red edge and a gradual rise below it. The principal difference between \cite{Wu-Goldreich-II} and our current ones is that the former find the theoretical red edge to occur at a $T_{\mathrm{eff}}$ approximately $1000\,\mathrm{K}$ greater than we do. It is likely that this discrepancy stems from differences in the parameterization of convection.

\section{Outbursts}\label{sec:burst}
Photometric outbursts from 6 cool DAVs were discovered in data collected by the Kepler satellite \citep{Bell-2015,Hermes, Bell-2016a, Bell-2016b}.\footnote{According to Hermes (private communication), the number of DAVs exhibiting outbursts currently stands at 12.}  Individual outbursts release $\sim 10^{33}-10^{34}\,\mathrm{erg}$ over several hours and recur aperiodically with separations of days. \cite{Hermes} suggested that outbursts may result from parametric instability \citep{Dziembowski,Wu-Goldreich-IV}.

The simplest case of parametric instability involves the nonlinear coupling of an overstable parent mode to two damped daughter modes. It requires the parent mode's amplitude, $A_p$, to exceed the threshold value
\begin{equation}
|A_p|^2=\frac{\gamma_{d_1}\gamma_{d_2}}{18\kappa^2\omega_{d_1}\omega_{d_2}}\left[1+\left(\frac{2\delta\omega}{\omega_{d_1}+\omega_{d_2}}\right)^2\right]\, ,
\label{eq:ampth}
\end{equation}
where $\delta\omega\equiv |\omega_p-\omega_{d_1}-\omega_{d_2}|$ is the frequency mismatch and $\kappa$ is the nonlinear coupling coefficient as defined in Equation~(9) of \cite{Wu-Goldreich-IV}. Here subscripts $p$, $d_1$, $d_2$ denote quantities pertaining to the parent and daughter modes. For $\delta\omega> |\gamma_d|\equiv |(\gamma_{d_1}+\gamma_{d_2})/2|$, parametric instability settles down on an equilibrium state in which the parent and daughter modes maintain constant amplitudes and the star's luminosity is slightly elevated. More stringent resonant couplings with $\delta\omega<|\gamma_d|$ yield limit cycles. We propose that they generate photometric outbursts.

Individual cycles consist of three stages. In the first, the parent mode's energy exponentiates. The second stage begins once the parent mode's amplitude has grown sufficiently such that the rate at which energy is transferred to its daughter modes overcomes their linear damping rates and their amplitudes begin to grow. At this point, the parent mode's energy is
\begin{equation}\label{eq:Eth}
E_{\mathrm{th}}\sim \frac{\gamma_{d_1}\gamma_{d_2}}{18 \kappa^2\omega_{d_1}\omega_{d_2}}\, .
\end{equation}
Subsequently, the parent mode's energy continues growing until eventually it drains into the daughter modes
as fast as convective driving can replenish it.\footnote{At maximum, the parent mode's energy may significantly exceed $E_\mathrm{th}$ \citep{Kumar-Goodman}.} This initiates the third and final stage in the limit cycle during which the parent mode's energy undergoes a precipitous decline. Once it drops below the threshold, the parent mode becomes too weak to feed the daughter modes, and thus the daughter modes'  mechanical energies also decay, since radiative damping converts them into thermal energy. The net result is the conversion of the peak mechanical energy of the parent mode to a pulse of thermal energy.  An outburst arises as the pulse of thermal energy leaks through the photosphere. This completes one cycle.

Growth of the parent mode's energy in the first and second stages occurs on timescale $\sim\gamma_p^{-1}$, and may not be visible; the parent mode's propagation cavity could be deeply buried below the CVZ such that $\omega_p\tau_c\gg 1$ and its photometric flux perturbation is severely reduced  (see equation~\ref{eq:dilution}), or the parent mode might have $\ell>1$ such that the angularly integrated flux perturbation is reduced. These stages corresponds to the quiescent interval between consecutive outbursts. The third stage, conversion of the parent mode's peak energy into heat by the daughter modes lasts for $\sim \gamma_d^{-1}$. The duration of an outburst is set by the maximum of $\gamma_d^{-1}$, $\tau_\mathrm{th}$ at the depth where the thermal pulse is deposited, and $\tau_c$.  The peak energy of the parent mode sets an upper limit to the energy emitted in an outburst.

As shown by Equation~(12) in \cite{Wu-Goldreich-IV},  magnitudes of 3-mode coupling coefficients are bounded from above by
\begin{equation}\label{eq:kappa-max}
|\kappa|_{\mathrm{max}}\sim (n_p^3 L_* \tau_{\omega_p})^{-1/2}\, .
\end{equation}
Here $n_p$ is the number of radial nodes in the parent mode, $L_*$ the stellar luminosity, and $\tau_{\omega_p}$ the thermal timescale at depth $z_{\omega_p}$. Thus the lower bound for the threshold energy is
\begin{equation}\label{eq:Ep-min}
E_\mathrm{th,min}\sim \frac{\gamma_{d_1}\gamma_{d_2}}{18 \omega_{d_1}\omega_{d_2}}\, n_p^3\, L_* \,\tau_{\omega_p}\, .
\end{equation}

Nonlinear effects constrain $\delta p/p<1$ everywhere throughout the star.  Since $\delta p/p$ increases outward as the background density drops, the upper bound on the parent mode's energy corresponds to $(\delta p/p)_{\mathrm{ph}}=1$.\footnote{In reality, it is difficult to determine the energy of the parent mode, because more descendant modes probably also participate in the energy cascade \citep{Kumar-Goodman}. \cite{Arras} propose that energy cascades from parent mode to descendant modes in a way similar to the Kolmogorov turbulent energy cascade. However, \cite{Brink} conclude that a small fraction of 3-mode couplings dominate the energy cascade and that the energy spectrum differs from Kolmogorov.} Figure~\refnew{fig:E-range} bounds the energy of prospective parent modes between an upper limit corresponding to $(\delta p/p)_{\mathrm{ph}}=1$ and the lower limit $E_\mathrm{th,min}$. Note that $(\delta p/p)^2_{\mathrm{ph}}=1$ corresponds to $E_p\sim n_p\tau_{\omega_p}L_*$  \citep{Goldreich-Wu-I} so $E_p$ rises with increasing $\omega_p$ because $\tau_{\omega_p}$ grows much faster than $n_p$ decreases.

\begin{figure}
\includegraphics[width=0.99\linewidth]{Eng-1.eps}
\includegraphics[width=0.99\linewidth]{Eng-2.eps}
\includegraphics[width=0.99\linewidth]{Eng-3.eps}
\caption{\label{fig:E-range} Energy range for prospective parent modes with $\ell=1,\,\, 2,\,\, 3$. The vertical grid lines label the frequencies of the prospective overstable parent modes. These are the modes with maximal growth rates at each $\ell$ (cf. Figure~\refnew{fig:timescale}). Horizontal lines label the typical energy released per outburst.}
\end{figure}

Next, we examine timescales for limit cycles to assess how they compare to recurrence rates and durations of outbursts. Figure~\refnew{fig:timescale} plots driving and damping rates, $\gamma$, for overstable and damped modes with $\ell =1,\, 2,\, 3$.  For each $\ell$, we define $\omega_\mathrm{pk}$ as the frequency at which the peak driving rate, $\gamma_\mathrm{pk}$, occurs.  Note that $\omega_\mathrm{pk}$ is only slightly larger than the minimum frequency for an overstable mode, and also that  $\omega_\mathrm{pk}\tau_c$ increases slowly with $\ell$, running from $6$ at $\ell=1$ through $8$ at $\ell=2$ to $9$ at $\ell=3$. More importantly, values of peak growth rates, $\gamma_\mathrm{pk}$ rise steeply with increasing $\ell$ as do the damping rates of prospective daughter modes.

We select a parent mode with $\omega\approx\omega_\mathrm{pk}$ along with a matched pair of near resonant daughter modes with $\omega\approx \omega_\mathrm{pk}/2$. To quantitatively account for observed outbursts, these must satisfy
\begin{eqnarray}
\gamma_\mathrm{pk} &\sim&(\mathrm{day})^{-1}\, ,  \label{eq:gamma-p-d}\\
\gamma_d&\sim &(\mathrm{hour})^{-1}\, , \nonumber\\
\delta\omega&<&\gamma_d\, .\nonumber
\end{eqnarray}
Horizontal lines in Figure~\refnew{fig:timescale} label levels at $1/\mathrm{day}$ and $1/\mathrm{hour}$.

Modes with $\omega\approx \omega_\mathrm{pk}$ are particularly well-suited for driving limit cycles. Radial orders of high $n$ g-modes roughly obey \footnote{See Figure 4 in \cite{Goldreich-Wu-I}.}
\begin{equation}
n \propto {[\ell(\ell+1)]^{1/2}\over \omega}\, .
\end{equation}
Thus the density of modes per unit frequency, $dn/d\omega\propto [\ell(\ell+1)]^{1/2}/\omega^2$, is larger at higher $\ell$ and lower $\omega$.  This is beneficial for finding combinations of near resonant parent and daughter modes with $\delta\omega<|\gamma_d|$ that are suitable for giving rise to limit cycles and thus producing outbursts.

Choosing potential parent modes likely to be responsible for the outbursts involves an interesting set of tradeoffs. The price for short limit cycles and short pulses is a smaller upper bound on peak parent mode energy. Currently, outbursts are only observed through flux increases and flux increases due to frequent short pulses may rival those from rarer long pulses although the latter have larger fluences.

Integrating Equations (1), (2) and (3) in \cite{Wu-Goldreich-IV},\footnote{We filter out the rapidly oscillating part of each modal amplitude by substituting a modified amplitude $A_j\exp(-i\omega_j t)$ ($j=p\, d_1\, d_2$) in these equations.} we present in Figures~\refnew{fig:outburst-22} and \refnew{fig:outburst-23} two examples of limit cycles involving the same $\ell=2$ parent mode coupled to different pairs of daughter modes.  In the former, the daughters are also $\ell=2$ and the frequency mismatch, $\delta\omega\ll \gamma_d$ is even smaller than $\gamma_p$.  In the latter, the daughters are identical twins with $\ell=3$ and although $\delta\omega<\gamma_d$, it is substantially larger than $\gamma_p$. In each example, recurrence times are in excess of 10 days. The limit cycles with two quadrupole daughter modes appear more irregular, a characteristic we observed in other cases with $|\delta\omega|<\gamma_p$. On the other hand, the limit cycle with two octupole daughter modes has $|\delta\omega|>\gamma_p$, and although not periodic, is more regular.\footnote{\cite{Wersinger} explored analogous behavior in a simple model of 3-mode coupling.}

\begin{figure}
\includegraphics[width=0.99\linewidth]{outburst-22.eps}
\includegraphics[width=0.99\linewidth]{dF-22.eps}
\caption{\label{fig:outburst-22} Limit cycles involving a quadrupole parent mode and two quadrupole daughter modes with oscillation periods $P_p\approx 12.1145\,\mathrm{min}$, $P_{d_1}\approx 23.8249\,\mathrm{min}$, and $P_{d_2}\approx 24.6537\,\mathrm{min}$. Radial orders are $n_p=28$, $n_{d_1}=58$ and $n_{d_2}=59$. Rates of driving and damping are $\gamma_p\approx 5.1017\times 10^{-6}\,\mathrm{s^{-1}}$, $\gamma_{d_1}\approx -9.6923\times 10^{-5}\,\mathrm{s^{-1}}$, and $\gamma_{d_2}\approx - 1.1678\times 10^{-4}\,\mathrm{s^{-1}}$. The frequency mismatch, $\delta\omega\equiv \omega_p-\omega_{d_1}-\omega_{d_2} \approx 1.17736\times 10^{-6} \,\mathrm{rad\, s^{-1}}$. The modal frequencies, radial orders, and driving/damping rates are calculated based on the WD model at $T_{\mathrm{eff}}\approx 10860\,\mathrm{K}$ and $g\approx 10^8\,\mathrm{cm\, s^{-1}}$. The upper panel plots the dimensionless $|\kappa A|$ for the parent mode (solid line) and the first daughter mode (dashed line). The lower panel plots their photometric flux perturbations, assuming the peak $|\kappa A_p|\approx 0.015$ in the upper panel corresponds to $\approx 10^{33}\,\mathrm{erg}$ for the parent mode.  }
\end{figure}

\begin{figure}
\includegraphics[width=0.99\linewidth]{outburst-23.eps}
\includegraphics[width=0.99\linewidth]{dF-23.eps}
\caption{\label{fig:outburst-23} Limit cycles involving a quadrupole parent mode coupled to identical octupole daughter modes with oscillation periods $P_p\approx 12.1145\,\mathrm{min}$, and $P_{d_1}=P_{d_2}\approx 24.2999\,\mathrm{min}$. Radial orders are $n_p=28$ and $n_{d_1}=n_{d_2}=82$. Rates of driving and damping are $\gamma_p\approx 1.1082\times 10^{-5}\,\mathrm{s^{-1}}$ and $\gamma_{d_1}=\gamma_{d_2}\approx -3.3758\times 10^{-4}\,\mathrm{s^{-1}}$. The frequency mismatch, $\delta\omega\equiv \omega_p-\omega_{d_1}-\omega_{d_2} \approx 2.5239\times 10^{-5} \,\mathrm{rad\, s^{-1}}$. The modes are calculated based on the same WD model as indicated in the caption of Figure~\refnew{fig:outburst-22}. The upper panel plots the dimensionless $|\kappa A|$ for the parent mode (solid line) and the first daughter mode (dashed line). The lower panel plots their photometric fluctuations, assuming the peak $|\kappa A_p|\approx 0.04$ corresponds to $10^{33}\,\mathrm{erg}$ for the parent mode.}
\end{figure}

\section{Conclusion and Discussion}\label{sec:conclusion}
Our paper clarifies the physics behind the red edge of the DAV instability strip.  It also provides an
explanation for the photometric outbursts recently discovered from cool DAVs.

Regarding the red edge, we recapitulate the necessary and sufficient condition for overstable modes, $\omega>\max(\tau_c^{-1},\, L_{\ell,b})$. Modes with angular frequencies below $L_{\ell,b}$ suffer enhanced radiative damping that overwhelms convective driving and are damped.
Below the red edge, all overstable modes have $\omega>L_{\ell,b}>1$ which renders them invisible to photometric surveys.  This completes the conceptual understanding of the observational red edge. High-quality spectroscopic observations of bright DAs with $T_{\mathrm{eff}}$ well below the red edge might reveal overstable modes through their surface velocities. This would be a valuable test for our explanation of the red edge.

We demonstrate that limit cycles due to sufficiently resonant 3-mode couplings quantitatively account for the energies, durations and repetition times of the outbursts observed by {\it Kepler}. Outbursting DAVs cluster close to the red edge because the condition for triggering limit cycles is preferentially satisfied there.

\section*{Acknowledgements}
We thank Yanqin Wu for her generous help, especially providing her programs and models for comparison with our work. We thank JJ Hermes, Chris Clemens, and Keaton Bell for comments on interpreting observations. Jing Luan is supported by the Theoretical Astronomy Center and Center for Integrative Planetary Science at University of California at Berkeley.


\section{Introduction}

The Kumbh Mela in India \cite{baranwal2015, barnett2016} and the Hajj in Arabia \cite{hughes2003, Kim2015, curtis2013virtual} are the two biggest periodic human gathering events on earth. Estimates have shown that $\sim 10^6$ to $10^7$ people gather into a confined space during these events. The dynamics of such large crowds and its possible spontaneous transition to a crush have perplexed researchers for over thirty years. Gatherings at a carnival, in a sports stadium or at a train station are no less susceptible to crowd disasters \cite{kok2016}. Such transitions are a matter of grave concern to both law enforcement and public health. We study these phase transitions from an orderly movement to a crush in large confined mobile crowds with the motive of proposing a control strategy to delay the onset of a spontaneous transition to a crush or even reverse this state back to ordered motion.

The social force based model introduced by Helbing and Moln\`ar \cite{helbing1995} and its variants have been shown to be capable of describing the emergent dynamics in human crowds \cite{vicsek2000nature,kwak2013,Helbing2006prl,ma2013,cirillo2013}. Most of these previous studies have relied on a combination of theory and agent-based simulations to study escape dynamics and evacuation efficiency of crowds through narrow openings \cite{helbing1995}. Such approaches are now being complemented by on-site computer vision studies \cite{kok2016,Kim2015}, cognitive science \cite{moussaid2011} and data analytics \cite{helbing2007pre,Silverberg2013}. For example, deployment of authority figures in escaping crowds in a metro station \cite{Song2017}, placement of obstacles near the escaping door \cite{Lin2017,Zuriguelpre2016}, mixing individualistic and herding behavior \cite{vicsek2000nature} are being proposed as means to enhance the evacuation efficiency. Of course, analysis of these specific situations and mechanisms to avoid crowd disasters \cite{kwak2013,vicsek2000nature,ma2013,Chraibi2010} have led us to an understanding of the underlying dynamics. However, it is much more promising if we can generalize the physics and use this knowledge to devise strategies to control crowd behavior.

\begin{figure}
\centering
\includegraphics[height=3.2in]{figs/seperatrix-crop.pdf}\label{fig:seperatrix
\caption{Bifurcation diagram showing the states of dynamical order and disorder in a human crowd. As the coordination coefficient, $\mu$, between agents decreases, a uniformly moving crowd with a finite magnitude of the averaged velocity $\lvert \langle \mathbf{v} \rangle \rvert$ loses order and changes to a state of crush, with $\lvert\langle \mathbf{v} \rangle \rvert \approx 0$. Insets (a) and (b) show two different choices of location of game changers as red points among other people marked as yellow. Inset (c) shows a measure of success (finite $\lvert\langle \mathbf{v} \rangle \rvert$ as $t \rightarrow \infty$) and failure ($\lvert\langle \mathbf{v} \rangle \rvert \rightarrow 0$ as $t \rightarrow \infty$) for a given impulse $\mathcal{I}$ compared to the time-averaged momentum in the ordered state $\mathcal{P}_o$ for the value of $\mu$ marked with $\times$. The separatrices delineating domains of attraction towards the ordered and disordered states are shown by dashed lines.}

\label{fig:figure2}
\end{figure}

Indeed, in this \textit{Letter}, we show that such a generalized approach towards crowd control is possible.
Fig.~\ref{fig:figure2} summarizes the main findings of this work.
We model human crowd in the frame work of active matter \cite{marchetti2013rmp1,vicsek1995prl,vicsek2000nature} and obtain two distinct crowd states, namely an ordered state of collective motion and a disordered state of crush. The ordered state with a finite magnitude of the averaged velocity $\lvert\langle \mathbf{v} \rangle \rvert$ is realized when coordination amongst the agents is high. Here, $\langle \hspace{0.05cm} \cdot \hspace{0.05cm} \rangle$ denotes average over all agents and over time. A transition to a disordered state where $\lvert\langle \mathbf{v} \rangle \rvert \approx 0$ is observed as a result of breakdown in agent coordination. We will show that transition from collective ordered motion to a state of crush happens as a non-equilibrium phase transition. More importantly we find that this transition can be reversed, without requiring a change in the overall level of coordination, by imparting an impulse to a small fraction of randomly chosen agents deemed as \textit{game changers}. As we will show, game changers work best when they are placed in regions of maximum crowd speed (or lowest local panic). For example, the optimal location of game changers is on a ring located at 70\% of the radius in a circular domain. We now describe the model, the implications of these findings for the safety of a human crowd and the general rationale for optimal game changer placement.

\textit{Model description}:
We use a dynamical variant of the well-studied agent based Vicsek model \cite{vicsek2000nature, vicsek1995prl} to simulate crowd dynamics.
This is an off-lattice model, where each agent is modeled as having a mass $m$ and occupying an area of a soft disc of diameter $d$. $N$ such agents are confined to move in a circular boundary. Each agent responds to three isotropic `interaction forces'. Thus the momentum balance for an agent $i$ can be written as
\begin{align}
m \frac{d\mathbf{v}_{i}}{dt} = \mathbf{F}^{pp}_{i} + \mathbf{F}^{sp}_{i} + \mathbf{F}^{D}_{i}
\label{eqn:forcebal}
\end{align}
where $\mathbf{F}^{pp}_{i}$, $\mathbf{F}^{sp}_{i}$ and $\mathbf{F}^{D}_{i}$ represent, respectively, the sum of all repulsive forces on the $i^{th}$ agent due to binary interactions with neighboring agents, a self propelling force generated by the $i^{th}$ agent and an alignment force on the $i^{th}$ agent due to agent-neighbor interactions. Agent inertia, which is not part of the classical Vicsek model \cite{vicsek1995prl} is included in Eq.~\eqref{eqn:forcebal}. The repulsive interaction force between the agent $i$ and a neighboring agent $j$ occurs due to space exclusion and is modeled as a linear soft spring: $\mathbf{F}^{pp}_{i}=\sum_j -k_n \boldsymbol{\delta}_{ij} \mathcal{H}\left( |\mathbf{r}_{i} - \mathbf{r}_{j}|-d \right)$.
Here, $\mathcal{H}(\cdot)$ is Heaviside function. If the position vector $\mathbf{r}_i$ points to the center of the $i^{th}$ agent, the extent of compression, $\boldsymbol{\delta}_{ij}$, is the vectorial distance along the separation vector between two agents $i$ and $j$,
$\boldsymbol{\delta}_{ij} = \left(|\mathbf{r}_{i} - \mathbf{r}_{j}|-d\right) \frac{\mathbf{r}_{i} - \mathbf{r}_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|}$. The coefficient $k_n$ determines the strength of the space exclusion force drawing from the analogous Herztian contact theory. The self propelling force $\mathbf{F}^{sp}_{i}$, produced by an agent has two components: one component aligned with the instantaneous velocity direction $\hat{\mathbf{v}}_{i}$ of strength $\beta$ and a second along a pre-determined motive direction $\hat{\mathbf{v}}^{m}_{i}$ of strength $\gamma$. Mathematically, $\mathbf{F}_{i}^{sp} =  m(\beta\hat{\mathbf{v}}_{i} + \gamma\hat{\mathbf{v}}^m_{i}) -\alpha |\mathbf{v}_{i}|\hat{\mathbf{v}}_{i}$ where $|\hat{\mathbf{v}}_{i}| = |\hat{\mathbf{v}}^m_{i}| = 1$. The parameter $\alpha$ is responsible for limiting the agent  to a terminal speed \cite{hinz2014motility}. The last force in the list is the collective crowd influence force experienced by each agent and is given by $\mathbf{F}^{D}_{i} =- \mu d (\mathbf{v}_{i} - {\mathbf{v}}_\mathbf{c})$, where $\mu$ is a co-ordination coefficient that controls the coupling strength between the $i^{th}$ agent and its neighborhood crowd. We calculate ${\mathbf{v}}_\mathbf{c}$ from a Gaussian weighted average of the velocities of all neighborhood agents in a radius $h$ \cite{pavan2015granular,Mahapatra2017}, chosen appropriately for dense crowds \cite{curtis2013virtual,Bechinger2016}.

We consider $N=6120$ active agents of diameter $d=0.5~m$ confined in a domain of circular shape of radius $R=22.5~m$ and mass $m=60~kg$, which results in a crowd density $\rho \approx 4~$persons$/m^2$ as observed in typical crowd conditions \cite{Parisi2009}. We select $k_n = 3 \times 10^{6}~N/m $ and $\beta = 1~m/s^{2}$. We have set $\alpha=\gamma=0$ without loss of generality of the conclusions. The conclusions also remain unaltered for systems with polydispersed agents \cite{sm1}. Typical agent speeds in the simulations match with those observed in literature \cite{Baglietto2011,helbing2007pre}. The influence radius for the neighborhood is set as $h=5d$ with the walls modeled as fixed agents.

\begin{figure}
\centering
\subfigure[]{\includegraphics[height=0.93in]{figs/half3p6-crop.pdf}\label{fig:half3p6}}
\subfigure[]{\includegraphics[height=0.93in]{figs/panic3p6-crop.pdf}\label{fig:panic3p6}}
\subfigure[]{\includegraphics[height=0.95in]{figs/speed3p6-crop.pdf}\label{fig:speed3p6}}\\
\centering
\subfigure[]{\includegraphics[height=0.93in]{figs/half3p4-crop.pdf}\label{fig:half3p4}}
\subfigure[]{\includegraphics[height=0.93in]{figs/panic3p4-crop.pdf}\label{fig:panic3p4}}
\subfigure[]{\includegraphics[height=0.95in]{figs/speed3p4-crop.pdf}\label{fig:speed3p4}}
\caption{Instantaneous fields of (a) velocity of agents, left half shows direction and right half shows the magnitude (b) panic factor and the (c) radial distribution of azimuthally averaged radial and angular components of velocity and panic factor $P$ for a system exhibiting ordered motion at $\mu=540~(Pa\cdot{s})$.  As $\mu$ decreases to $510~(Pa\cdot{s})$, the system transits to a disordered state. Correspondingly (d) velocity, (e) panic factor $P$ and (f) averaged components of velocity and panic factor. The inset of (f) shows the instantaneous velocity profiles.
}\label{fig:figure3}
\end{figure}

\textit{Dynamic states and phase transition}:
Two dynamically stable states of the system are obtained from our simulations. They are described by Fig~\ref{fig:figure3} (see SM \cite{sm1} for videos).
At high values of the coordination coefficient $\mu$, active agents organize into an ordered velocity field (see Fig.~\ref{fig:half3p6}). The local orientation of this ordered velocity field is determined by the geometry of the confinement. In our simulations it manifests as a rotating velocity field induced by the circular confinement, a preferred state observed in other active systems as well \cite{wioland2013prl,lushi2014, czirok1996formation}. In periodic domains the ordered velocity field corresponds to a flocking state \cite{sm1} consistent with previous research \cite{chate2008collective}. From this state, as $\mu$ is reduced, collective motion disappears, the order in the velocity field is lost and we obtain a disordered state
(see Fig.\ref{fig:half3p4}).

The two distinct states also  show different velocity profiles. In Fig.~\ref{fig:speed3p6}, the velocity field is dominated by its azimuthal component which increases from zero at the center to a maximum and reduces to zero again near the walls. On the other hand in Fig.~\ref{fig:speed3p4}, which corresponds to a disordered state, both the radial and azimuthal components of velocity fluctuate around zero mean and are small when compared to those observed under collective motion. This state where agents' motion is dictated by individual choice is referred to as the state of crush. In this state, the thrust force by each agent is balanced by the drag force, $m\beta \sim \mu d |\mathbf{v}_{i}|$
resulting in individual agent speed $|\mathbf{v}_{i}| \sim \frac{m \beta}{\mu d}$.

A brief discussion of the order-disorder transition is presented next.
Unlike previous variants of the Vicsek model \cite{vicsek1995prl,marchetti2013rmp1}
where disorder is brought in via a random force, physical or virtual collisions (generated by elastic repulsive forces) acting in consonance with agent inertia is the source of disorder in our model. Activity further amplifies this noise in the system. An increase in coordination between agents reduces this deterministic noise \cite{sm1} and beyond a threshold value, we obtain ordered crowd motion. On the other hand, when noise dominates, the system goes into a state of crush

In order to describe such a transition from an individual choice to a crowd dominated choice Helbing et al. \cite{vicsek2000nature} have postulated a \textit{panic factor}, $P$. Casting their \cite{vicsek2000nature} equations in our form, we identify the panic factor as
\begin{align}
P =\frac{m \beta}{m \beta+\mu d \lvert \mathbf{v_c} \rvert}
\label{eqn:panic}
\end{align}
Thus panic factor is a measure of the individualistic behavior ratioed to the crowd behavior. In a panicked state, as observed in a crush, selfish active agents exhibit individualistic behavior and self propelling force will exceed the coordination force since $\mathbf{v_c} \approx 0$. Then $m\beta >> \mu d \lvert \mathbf{v_c} \rvert$ and $P \approx 1$. On the other hand, when the crowd is in an ordered state, the self-propelling force of individual agents would be just sufficient to match the forces exerted by the neighboring crowd: $m\beta \approx \mu d \lvert \mathbf{v_c}\rvert$ and $P \approx \frac{1}{2}$. Since $\mathbf{v_c}$ is a local quantity, Eq.~\eqref{eqn:panic} is  a generalization of the constant panic factor postulated by \cite{vicsek2000nature}. Thus, our model allows the crowd to have spatially varying choices. The panic factor corresponding to the ordered and disordered states are shown in Fig.~\ref{fig:panic3p6} and \ref{fig:panic3p4} respectively. As expected, $P\approx 1$ in the disordered state as opposed to $P \approx \frac{1}{2}$ in spatial regions of high order. The azimuthally averaged panic factor corresponding to the two states are also shown in Fig.~\ref{fig:speed3p6} and \ref{fig:speed3p4}. Notably, $\frac{1}{2} \le P \le 1$ in the Fig.\ref{fig:speed3p6}.

In order to quantitatively describe these two states, we define the magnitude of the averaged velocity of the crowd, $\lvert\langle \mathbf{v} \rangle \rvert =  \frac{1}{T}{}\int_{0}^{T} \frac{1}{N}\sum_{i} \mathbf{v}_{i} dt$ as the order parameter.
Fig.~\ref{fig:figure2} shows the variation of $\lvert\langle \mathbf{v} \rangle \rvert$ as a function of $\mu$ for a fixed active forcing ($\beta =$ constant) and crowd density $\rho$. At large values of $\mu$, where collective motion dominates, the average velocity $\lvert\langle \mathbf{v} \rangle \rvert$, of the crowd is weakly dependent on $\mu$.
The collective motion persists even when $\mu$ is lowered but only upto a critical value, $\mu_1$. At $\mu_1$ there is a sharp decrease in the order parameter indicating a phase transition. The collective motion disappears and the agents exhibit diffusive motion. Further reduction in $\mu$ maintains $\lvert\langle \mathbf{v} \rangle \rvert \approx 0$. Reversing the experiment by increasing the value of $\mu$ gives a different $\mu_2$ for the reverse transition, namely order from a crush.

\textit{Safety - Recovery from a crush via an induced transition}: So far, we analyzed the order - disorder transition of a crowd using $\mu$ as the control parameter. However, as Fig.~\ref{fig:figure2} indicates, the critical values for the transitions are different, $\mu_1 \ne \mu_2$. In other words, the system exhibits hysteresis. An interesting outcome of this observation is that two distinct stable dynamical states exist for $\mu_1 \le \mu \le \mu_2$. This presents an opportunity to investigate recovering order from a crush.

We start by considering a system exhibiting a state of crush for a particular $\mu$ in the hysteresis region, for example, the point marked as `$\times$' in Fig.~\ref{fig:figure2}. As can be seen, there also exists a state of order at this same level of coordination between the agents. We now empower game changers in the system and study whether the system can be driven to an ordered state. Game changers can be people either  employed at these locations a priori or they can be promoted agents. Game changers differ from the rest of the crowd only in their ability to control their velocity for a short duration of time. The distribution of game changers can be selected in different ways; we discuss two of them here: game changers located in a stripe at various radial locations (Fig.~\ref{fig:figure2}(a)) and uniformly dispersed throughout the domain (Fig.~\ref{fig:figure2}(b)).

In short, an impulse $\mathcal{I}=N_g m \gamma \tau$ is imparted to a state of crush through $N_g (=10\%N)$ game changers by assigning them a tangential motive force $m\gamma$ for a certain interval of time $\tau$. For the rest of the crowd, $\gamma=0$. The force is then removed and the system is allowed to evolve. The response of the system to varying impulse $\mathcal{I}$ corresponding to the point marked `$\times$' in Fig.~\ref{fig:figure2} is shown as an inset (c) in the same figure. It may be seen that $\lvert\langle \mathbf{v} \rangle \rvert$ steadily increases for $t > \tau$ when ${\mathcal{I}}/{\mathcal{P}_{o}}\geq 0.67$, indicating spontaneous restoration of order in the system after the forcing has been turned off. Here $\mathcal{P}_o=Nm{\lvert\langle \mathbf{v} \rangle \rvert}_o$ is the total momentum in the ordered state and $\lvert\langle \mathbf{v} \rangle \rvert_o$ is the average speed in the ordered state.
On the other hand, when ${\mathcal{I}}/{\mathcal{P}_{o}}< 0.67$, $\lvert\langle \mathbf{v} \rangle \rvert$ reduces towards zero when the forcing is removed (for $t>\tau$). This implies that the system fails to recover from the crush despite the impulsive forcing.

Of course, successful recovery of the system from a state of crush depends upon (i) the strength of the impulse ($\mathcal{I}$) provided by game changers and (ii) the distribution of game changers - location, geometry and dimensions of the mode (for e.g.\ stripe vs dispersed). For $\mu_1 \le \mu \le \mu_2$, we have systematically identified the minimum value of $\mathcal{I}$ for which transition to order is achieved, for any possibility in $\frac{N_g}{N}$ and distribution of game changers. Thus we obtain the separatrices, in Fig.~\ref{fig:figure2} shown as dotted lines that the game changers have to move the system to, before spontaneous transition to the ordered state sets in.
When the system is impulsively forced to an order parameter value above (below) this curve, game changers succeed (fail) in recovering order. As $\mu$ increases, $\lvert\langle \mathbf{v} \rangle \rvert$ (and thus the corresponding impulse) required for the system to transit to an ordered state decreases. We also note that the separatrix that lies closest to the disordered state is for the case where the game changers are distributed on a ring located at $\approx 70\%$ radius of the domain. Correspondingly,  $\lvert\langle \mathbf{v} \rangle \rvert$ (and thus the corresponding impulse) required is always a minimum for this particular spatial distribution of game changers

In other words, locating game changers at $70\%$ of the domain radius is the optimal way to recover from a crush in a circular domain without changing the level of coordination between the agents. This is explained by studying the radial velocity distribution as well as panic factor in the ordered state. At approximately $0.7R$, the azimuthal velocity ($|V_\theta|$) of the agents\, and therefore $\mathbf{v_c}$, are both at their maxima, as seen in Fig.~\ref{fig:speed3p6}. Hence the power input to the system ($\mathcal{I}|\mathbf{v_c}|$) through the organized motion of the game changers is maximum. Therefore, the impulse required to transform the system is minimum when they are placed at the locations of maximum crowd speed. Maximum $|\mathbf{v_c}|$ at this location also implies that the panic factor (and thus individualism) is at a minimum.

This is an important result with regard to crowd safety. In order to recover from a crush without changing the level of coordination between the agents while imparting the least impulse, game changers would have to be placed at the locations of maximum crowd speed which correspond to locations of lowest panic in the domain. In case of an accident like crush, these game changers can impart power to the crowd in identified directions, thus helping align agent motion throughout the domain.

We have also studied the transition from an ordered state to a disordered state. Here the role of game changers is observed to be similar to the role of a small fraction of rogue elements in flocking models studied earlier \cite{Ariel2015, yllanes2017many}. We find that breaking order via rogue elements is easier than its restoration.

\begin{figure}
\centering
\subfigure[]{{\includegraphics[scale=0.17]{figs/new-crop.pdf}\label{fig:orderp}}}
\subfigure[]{\includegraphics[scale=0.16]{figs/binder-crop.pdf}\label{fig:binder}}
\caption{(a) During the order-disorder transition in a crowd, the order parameter $\lvert\langle \mathbf{v} \rangle \rvert$ shows a discontinuous change at $\mu = \mu_1, \mu_2$. Different sets correspond to different system sizes. (b) Change in the value of Binder cumulant from $1/3$ to $2/3$ through a dip at $\mu = \mu_2$ during the transition.}

\end{figure}

\textit{First order phase transition and hysteresis:} We now discuss the reasons that lead to the basic premise that the recovery from a crush via game changers is possible and we show that this possibility exists due to the underlying nature of the phase transition.
The sharp change in the order parameter in Fig.~\ref{fig:orderp} at $\mu=\mu_1, \mu_2$ indicates that the phase transition from an ordered state to a crush is first order. We confirm this by calculating the Binder cumulant of the order parameter through the transition while moving to a state of order from a state of crush. We define the Binder cumulant $G$ following Chat\'e et al. \cite{chate2008collective} and plot it as a function of the control parameter $\mu$.
Fig.~\ref{fig:binder} shows this plot for various domain sizes, $R$, for a fixed $\rho$. $G$ varies from  $2/3$ to $1/3$ from an ordered state to a disordered state as is known for two dimensional systems.
During the transition, Binder cumulant undergoes a sharp dip at the transition point proving that crowd crush is indeed a first order phase transition from a uniformly moving crowd. Moreover the size and the sharpness of the dip also increases with an increase in the system size, another feature of a first order phase transition \cite{chate2008collective}. Thus, our simulations underscore the discontinuous nature of order-disorder transition seen in Vicsek-like models of active systems \cite{chate2008collective}, additionally under confinement.

Hysteresis is a hallmark of first order phase transitions. Therefore, near the transition point,  both ordered and disordered states are solutions at the same $\mu$. This existence of multiple solutions is the reason that opens a pathway to reverse a crush state through mobilization of game-changers. It appears that the role played by game changers in restoring order is analogous to nucleation in equilibrium phase transitions \cite{binder1987theory,hohenberg1995metastability,csernai1992nucleation}. Our analysis also shows that it would not be possible to recover the state through game changers for $\mu < \mu_1$. On the other hand, the width of the hysteresis loop increases with an increase in domain size (see Fig.~\ref{fig:orderp}), providing better scope for recovery in case of an accident in large gatherings.

\begin{acknowledgments}
The authors would like to gratefully acknowledge discussions with Srikanth Vedantam, Silke Henkes and Hugues Chat\'e.
\end{acknowledgments}


\section{Introduction}
Despite a history stretching over five decades, magnetic colloids continue to find new applications. Paramagnetic beads are highly attractive bases for dynamic materials due to the relative ease and precision with which researchers can control magnetic fields in many media. Magnetic colloids are being used for medical tasks such as drug delivery \cite{drugdelivery}, tissue scaffolding \cite{tissue2}, image contrast \cite{contrastagents}, and tumor reduction \cite{hyperthermia}. More exotic functions include self-assembling swimmers \cite{swimmers1,swimmers3}, walkers \cite{walkingribbon1,walkingribbon2}, grabbers \cite{asters}, and self-healing membranes \cite{magicangle}.
\\
Both medical \cite{bettermed} and non-medical \cite{lumpactuator,wireflow1,wireflow2} applications  benefit from using particles with extended aspect ratios. 1D chains of magnetic colloids, or magnetic filaments, are of particular interest on account of their rich behavior arising from the interplay of their elastic and magnetic properties \cite{wei2016}. Chain synthesis techniques have advanced steadily over almost two decades \cite{firstchains}, and a wide variety of bead-linking methods are now available to scientists \cite{squishchains, chainbundle, onepot, blockingtemp, ferromagnetic, hollowchain, branching, fattyassembly, anisotropicassembly, streptavidinlink}. In recent years a variety of research teams have studied chain behaviors such as actuation via bending \cite{bending, magnetictorque, coreshell}, buckling transitions \cite{overview1}, desynchronization in slowly rotating fields \cite{rotation}, swimming \cite{swimmers2}, and pumping \cite{pump, pumptheory}.
\\
Most work to date has focused on static \cite{Lusebrink2016, Cerda2016} or relatively slow field changes with single filaments \cite{overview2}. In this work, we combine analytics and simulations to study the behavior of isolated and bulk filaments in dynamic magnetic fields. Our goal is not to characterize the equilibrium ensemble of bulk filaments, as thermodynamic results will depend strongly on the microscopic details of the magnetic filaments. Instead, we focus on exploiting dynamic effects common to many types of filaments in order to produce desired behaviors. In the case of single filaments under fast precession, we find that for large precession angles the filaments naturally form helices with harmonic potentials. These helices can be used for well-controlled contraction at the microscale. For multi-chain ensembles, we show that aggregation behavior is a function of two effects. The dominant effect is the effective interaction of beads in the bulk, which is determined by the precession angle. The second effect is a frequency-dependent residual torque on chain ends. Combining these effects with more complex field choices, such as precession with a dynamically changing precession angle, enables a much greater variety of behaviors. We illustrate these possibilities by driving chains to twist together into an effectively cross-linked gel.

\section{Model}
For computational purposes, we treat the chains as collections of hard, spherical magnetic dipoles bound to their neighbors by flexible but inextensible bonds. The flexibility of the chain is governed by the bending rigidity $\epsilon$ through the harmonic angle potential
\begin{equation}
	U_B (\delta\theta) = \frac{\epsilon}{2} (\delta\theta)^2 \,,
	\label{eq:theta}
\end{equation}
where $\delta\theta$ is the angle formed by a bead and its two neighbors. In addition to hard-sphere interactions, beads also interact via the dipole-dipole potential
\begin{equation}
	U_{ij} = \frac{1}{r_{ij}^3} \left( \bm{\mu}_i \cdot \bm{\mu}_j - \frac{3}{r_{ij}^2} (\bm{\mu}_i \cdot \mathbf{r}_{ij}) (\bm{\mu}_j \cdot \mathbf{r}_{ij}) \right) \,,
	\label{eq:dipole_general}
\end{equation}
where $\bm{\mu}_i$ is the dipole moment of a bead and $\mathbf{r}_{ij}$ the displacement between beads. In MD simulations we take the beads to be freely rotating spheres with dipole moments of fixed magnitude, appropriate to systems such as Ref. \citen{hollowchain} in a saturating external field. We refer to this below as the computational model of the filaments. In the computational model, all magnetic interactions between dipoles less than 10$\sigma$ apart are considered.

The basic behavior of the filaments is readily described. For the regimes considered here we find that magnetization remains uniform. The dominant component of Eq. (\ref{eq:dipole_general}) is then
\begin{equation}
	U(\mathbf{r}, \bm{\mu}) = \frac{\mu^2}{r^3} \left( 1 - 3 (\hat{\bm \mu} \cdot \mathbf{\hat{r}})^2 \right) \,.
	\label{eq:dipole}
\end{equation}
For the case of a static field, free filaments will align with the field. Parallel filaments have an $r^{-3}$ repulsive interaction when side by side, and an oscillatory attraction that decays exponentially with a length scale set by bead size. Consequently the filaments will eventually aggregate into ribbons and columns, but on time scales too long to interest us here.
\\
Rapidly changing the field orientation frustrates the filament's alignment with the field. Consider a field precessing around the $z$ axis so that the magnetic moment directions are given by $\hat{\bm \mu} = \cos \omega t \, \sin \beta\, \mathbf{\hat{i}} + \sin \omega t \, \sin \beta\, \mathbf{\hat{j}} + \cos\beta\, \mathbf{\hat{k}}$. If the field precesses with a short period compared to the characteristic time scales for bead translation, we are justified in treating the conformation of the filament as quasistatic with respect to the precession. We may average Eq. (\ref{eq:dipole}) over the precession period to obtain the quasistatic magnetic energy
\begin{equation}
	U(r, \beta, \alpha) = \frac{\mu^2}{r^3} \left (1 - 3 \left( \cos^2{\beta}\cos^2{\alpha} + \frac{1}{2} \sin^2{\beta} \sin^2{\alpha}\right) \right)\,,
	\label{eq:precess}
\end{equation}
where $\beta$ is the angle of the magnetic field with its precession axis $z$, and $\alpha$ is the angle of the displacement vector with $z$. For small values of $\beta$, Eq. (\ref{eq:precess}) drives configurations similar to a static field aligned with $z$. For $\beta$ values near $\pi/2$, dipoles prefer to form sheets in the $x-y$ plane, and repel in the $z$ direction \cite{Martin2000}. At the magic angle $\beta_m = \cos^{-1}\sqrt{\sfrac{1}{3}}$, the $r^{-3}$ potential vanishes and interactions are dominated by dipole correlations that scale as $r^{-6}$. Reference \citen{magicangle} discusses the behavior of free colloids at the magic angle in depth. In the present work, we verified that filaments fail to aggregate or align at the magic angle for the parameters used in our system. Instead, we observed dominant hard-sphere interactions. We conclude that the $r^{-6}$ potential is weak compared to thermal forces and does not appreciably contribute to conformations or aggregate formation in the regimes studied here
\\
We assume that the bending modulus $\epsilon$ is large compared to the thermal energy $kT$ and magnetic energy $\mu^2\sigma^{-3}$, so that the chain follows a smooth curve with small $\delta\theta$ everywhere. With this assumption we can now develop a continuum model of the filaments, which is complementary to the computational model above. Let the unit vector $\mathbf{T}(s)$ be the tangent of the filament curve at bead $s$, or equivalently the displacement vector pointing from $s$ to the next bead $s+ds$. Denoting bead diameter with $\sigma$, we write the interaction between neighboring beads as
\begin{align}
	U_{NN} = \frac{M(\beta)}{2}T_z^2\,, \quad
	M(\beta) = -\frac{9\mu^2}{\sigma^3} \left( \cos^2 \beta -\frac{1}{3} \right) \,,
	\label{eq:NN}
\end{align}
where $T_z = \mathbf{T} \cdot \hat{\mathbf{z}} = \cos \alpha$. The next-nearest neighbor interactions have leading prefactor $\sfrac{1}{8}$. An expansion in $\delta \theta$ shows that the leading term is simply a rescaling of the nearest-neighbor energy by $1+\sfrac{1}{8}$. The next term is proportional to $\sfrac{3}{16} \mu^2 \sigma^{-3} (\delta \theta)^2$. For the regimes considered here, this term is negligible compared to $U_B$. Accordingly we neglect interactions beyond nearest-neighbor for the continuum model.
\\
We may now write the zero-$T$ Hamiltonian for a single static strand of length $L$ forming a continuous curve $\mathbf{x}(s)$:
\begin{equation} \label{eq:continuum}
	H = \int_0^L ds\, \left( \frac{\epsilon}{2} (\mathbf{T}')^2 + \frac{M(\beta)}{2} (T_z)^2 \right) \,,
\end{equation}
where the first term gives the bending energy and the second gives the magnetic energy.In the following, we rescale lengths in terms of the bead diameter $\sigma$, v.g. $s \rightarrow s/\sigma$ and the total length $L$ corresponds to the number of beads in the filament. Equation (\ref{eq:continuum}) is similar to previous continuum studies of magnetic filaments, but generalized to the case of fast precession \cite{overview1}. Despite their  simplifications, continuum models well describe experimental results for magnetic filaments in static fields, including U and S-shaped hairpin turns, \cite{overview2}. In the following sections we will use MD simulations of the computational model with a Langevin thermostat to both complement and verify the continuum model in the case of rapid precession and to bound its applicability.

\section{Results}

\begin{figure}
	\centering
	\includegraphics[width=.95\columnwidth]{fig1.pdf}
	\caption{Helical contraction. (a) Schematic of the boundary conditions on a contracting helix. (b) Frame from MD simulation using LAMMPS with a 40 bead magnetic filament connecting two inert test masses. Under the influence of a precessing field with $M > 0$, the helix attempts to contract to the $x-y$ plane. Here $\epsilon = 40$, $\mu^2\sigma^{-3} = 4 kT$, and $\beta=\pi / 2$. (c) If $M < 0$, the filament attempts to align with the $z$ axis and forms planar solutions. Both bending and magnetic energy costs are localized to hairpin turns, so the state is metastable.}
\end{figure}
The combination of the magnetic and bending energies can be exploited for controllable actuation. Suppose that each end of the filament is anchored to two test loads a distance $d$ apart. In Appendix \ref{sec:Helices} we show that there are two families of solutions for this Hamiltonian with fixed boundaries: helices about the $z$ axis [Figs. 1a and 1b], and planar curves that include the $z$ axis [Fig. 1c]. We find from simulations that the one-coil helices are the stable solutions for positive $M$ ($\beta > \beta_m$). For negative $M$, ($\beta < \beta_m$) the preferred solutions are planar with hairpins that concentrate bending energy, similar to the hairpins that occur for free ends in Ref. \citen{overview1}.
\\
The total energy of the helix is given for arbitrary initial conditions in Appendix \ref{sec:Helices}. Let us take the ideal case, where the start- and end-points of the chain are aligned with the $z$ axis. For a helix, $T_z = d/L$ and $(\mathbf{T}')^2 = 4\pi^2L^{-2}(1-d^2/L^2)$. Integrating Eq. (\ref{eq:continuum}) we get
\begin{equation}\label{eq:spring}
	H = \frac{k}{2} d^2 + \frac{2\pi^2\epsilon}{L} \,, \quad 	k = \frac{M(\beta)}{L} - \frac{4\pi^2\epsilon}{L^3}\,.
\end{equation}
This potential is harmonic in the separation of the end points $d$ for any (physical) value of $d$.  Thus, the helix behaves as a spring with a force constant $k$ that can be adjusted by simply by changing the precession angle $\beta$. If $k$ is positive, the harmonic favors contraction. This offers remarkably precise control over actuation for a given angle.
\\
We verified the formation of harmonic contractile helices with MD simulations in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package. Details regarding simulation methods can be found in Appendix \ref{sec:mdsim}. Movie 1 in Supplemental Material illustrates the results.

\begin{figure*}
	\centering
	\includegraphics[width=.99\textwidth]{fig2.pdf}
	\caption{Effects of precession frequency and bending modulus on aggregation for 80-bead filaments with $\mu^2 \sigma^{-3} = 4kT$ and $\gamma \sigma=10$ (Lennard-Jones units). For sufficiently low frequencies and bending moduli, the torque on filament ends forces filaments to aggregate in metastable spirals instead of branching networks.}
\end{figure*}
\vskip0pc \noindent
We turn to the case that the filament ends are not fixed, so that it is possible for multiple filaments to aggregate. If $\beta < \beta_m$, the filaments align with the precession axis in a manner similar to free dipolar colloids. If $\beta > \beta_m$, the filaments arrange themselves in the plane orthogonal to precession and readily aggregate. In this second case, their conformation and aggregation behavior depend on the speed of precession $\omega$. Here we consider the high-(but finite-) frequency regime and determine the minimum necessary frequency for the validity of the Eq. (\ref{eq:precess}).
\\
For high frequencies the static term given by Eq. (\ref{eq:precess}) remains the dominant energy for individual colloids, but additional force terms arise when the chain is permitted to move in response to finite-frequency precession. For simplicity, consider the case that $\beta=\sfrac{\pi}{2}$, {\it i.e.} that the magnetic field rotates entirely in the plane of the filament. The path of the filament is parametrized by the angle $\psi(s)$ its tangent forms with some axis in the $x-y$ plane. The angle between the tangent and the field is instantaneously $\omega t -\psi$, and we find the nearest-neighbor force $\mathbf{F}_{NN} = -\nabla U(\mathbf{r}, \mathbf{\mu})$ in cylindrical coordinates:
\begin{equation}
\mathbf{F}_{NN} = 3\frac{\mu^2}{\sigma^4} \left[ \sin (2\omega t -2\psi)\mathbf{N} + (1-3 \cos^2 (\omega t - \psi))\mathbf{T} \right] \,,
	\label{eq:FNN}
\end{equation}
where $\mathbf{N}$ is the unit vector in the plane perpendicular to $\mathbf{T}$. For beads in the middle of the filament, the net force is proportional to $\delta \psi$ and amounts to an oscillating adjustment to the bending rigidity $\epsilon$. However, the two end beads have a finite force. Discarding the force in the $\mathbf{T}$ direction, and including the viscosity of the medium, the net force on the ends of a straight filament is given by
\begin{equation}
	\mathbf{F}_{end} = 3\frac{\mu^2}{\sigma^4} \sin (2(\omega - \dot{\psi})t)\mathbf{N} -\gamma \dot{\psi} \mathbf{N} \,,
	\label{eq:Fend}
\end{equation}
where $\gamma$ is the viscous drag coefficient. The magnetic force oscillates back and forth as the field rotates. We are interested in the cycle-averaged motion
\begin{equation}
	\langle \dot{\psi} \rangle = \frac{\omega}{\pi} \int_0^\frac{\pi}{\omega} dt \dot{\psi} \,,
\end{equation}
In a Brownian regime, the instantaneous total force $\mathbf{F_{end}}$ is zero. Further, in the high-frequency regime the motion of the filament in a single precession cycle is small. We write $\psi=\psi_0+\delta \psi$, and set $\psi_0=0$ without loss of generality. Expanding the force to linear order in $\delta \psi$, we find
\begin{equation}
	\dot{\psi} = \frac{1}{\tau} \big ( \sin (2\omega t) - 2\delta \psi \cos (2\omega t) \big)
	\label{eq:psidot}
\end{equation}
where we have introduced the magnetoviscous timescale $\tau = 1/3  \gamma \,\mu^{-2} \sigma^{4}$. On averaging over one period of the motion, the left term vanishes. The right term may be integrated by parts:
\begin{equation}
	\langle \dot{\psi} \rangle = \frac{ 1}{ \pi \tau} \int_0^\frac{\pi}{\omega} dt  \dot{\psi} \sin (2\omega t)
	\label{eq:recursion}
\end{equation}
One may now recursively substitute Eq. (\ref{eq:psidot}) into Eq. (\ref{eq:recursion}), integrate the first term, and integrate the second term by parts. Concretely, the first such iteration generates
\begin{align}
	\langle \dot{\psi} \rangle &= \frac{ 1}{ \pi \tau} \int_0^\frac{\pi}{\omega} dt \frac{1}{\tau} \big ( \sin (2\omega t) - 2\delta \psi \cos (2\omega t) \big) \sin (2\omega t) \nonumber \\
	&= \frac{ 1}{ \pi \tau} \left( \frac{\pi}{2\omega\tau} + \frac{1}{4\omega \tau} \int_0^\frac{\pi}{\omega} dt  \dot{\psi} \sin (4\omega t) \right) \nonumber \\
	&= \frac{1}{\tau} \left( \frac{1}{2\omega \tau} + \mathcal{O} (\omega \tau) ^ {-2} \right)
	\label{eq:first_recurse}
\end{align}
The net motion of the filament ends against the drag force of the solvent may be identified with an average magnetic force $\langle \mathbf{F}_{NN} \rangle = \gamma \dot{\psi} \mathbf{N}$, which vanishes with $(\omega \tau)^{-1}$. This net force generates a magnetic torque that tends to wind the two ends of the filament in curves of opposite chirality. The magnetic torque competes with the bending modulus to determine the conformation of the filament aggregates. For large values of $\omega$, the bending modulus dominates and filaments aggregate in branching networks. For smaller values, the magnetic end force succeeds in forcing the filament ends to curve until they make contact with the body of the filament. Magnetic attractions between beads then cause the filament to roll into a spiral, and separate spirals aggregate in compact clusters.  Figure 2 illustrates the resulting aggregate metastable ``phases'', characterized in terms of the inverse magnetoelastic parameter $\xi = \epsilon / (ML^2)$ and magnetically normalized end force $\chi = M / (\gamma \omega \sigma)$, with $M= 3 \mu^2/\sigma^3$, as generated by the MD model.
\\
For large values of the bending modulus the end force never rolls the filament in a spiral, instead rotating the entire filament coherently (top right). For small bending moduli, the filaments tend to naturally roll and fold even at high frequencies (bottom left). At low frequencies, the quasistatic approximation for the filament bulk fails and the filaments form rapidly rotating compact shapes (bottom right). For intermediate values of $\xi$, there is a continuous transition from branching to spiral conformations in the region $500 \le \chi/\xi \le 800$.  In all cases, filaments were initially thermalized with no magnetic interactions. Note that since the aggregates are metastable, alternative initial conditions will alter the diagram.
\\
It is now possible to bound the validity of Eq. (\ref{eq:precess}) at $\chi \le 1$. For 100 nm-diameter filaments in room temperature water with the parameters quoted in Fig. 2, this requires a precession frequency of 1 kHz. The required frequency falls with increasing diameter as $\sigma^{-3}$.
\begin{figure*}
	\centering
	\includegraphics[width=.99\textwidth]{fig3.pdf}
	\caption{Achieving twisting with a magnetic drive. (a) Path traced by the driving field with $\omega=2.5$, $\omega_2=1.0$, and $\beta_0=0.5$ over $t \in [0, 4\pi]i$ (top) and $[0, \pi]$ (bottom). (b) The twisted cross-links created by this drive with the parameters $\mu=2.0$ and $\epsilon = 80.0$. Further simulation details are given in the SM. (c) Increasing chain length and system density creates a system-spanning network of cross-linked filaments. The succeeding panels are the results of altering the magnetic drive parameters as noted: (d) $\beta_0 = 0.4$; (e) $\beta_0 = 0.6$; (f) $\omega_2 = 2.0$; (g) $\omega_2 = 0.5$; (h) $\omega = 25$, $\omega_2 = 10$; (i) $\omega = 0.25$, $\omega_2 = 0.1$.}
	\label{Fig:3}
\end{figure*}
More complex aggregates are possible with a more general form of driving field than simple precession. A truly arbitrary field opens an enormous design space. Here, we confine ourselves to illustrating how one may design the field to achieve a desired end state. We select as our target material a self-assembling 3D gel in which strands have twisted around one another to effectively cross-link. Such a gel is much stronger than one that relies on magnetic attractions alone, and is of special interest due to its dynamic nature and potential for self-healing.
\\
For an isotropic gel, we require $\langle M \rangle \approx 0$ while retaining a nonzero $r^{-3}$ attraction between filaments. To achieve this we use a precessing field as before, but with a time-dependent precession angle:
\begin{equation}
	\beta(t) = \frac{\pi}{2} - \left( \frac{\pi}{2} - \beta_0 \right) \cos^2 \omega_2 t
	\label{eq:beta}
\end{equation}
and set $\beta_0$ so that the precession angle roughly averages to the magic angle $\beta_m$, or $\beta_0 \approx 0.5$. Obtaining the desired twisting behavior is more challenging. The end points of the filaments should trace a circular pattern around the filament axis, and should do so regardless of how the filament is aligned in the laboratory frame. Therefore, the field must trace open loops on the unit sphere at frequencies low enough for the filament ends to respond. Guided by the preceding frequency analysis of the filament response to finite frequencies, we choose $\xi=0.001$ and $\chi \approx 1$. In order to ensure an isotropic gel, each filament end must respond to the precession regardless of filament orientation. We therefore choose $\omega = 2.5$, $\omega_2 = 1$, so that the field forms a series of open loops with different orientations. Figure 3(a) illustrates the path the field traces on the unit sphere with these parameters.
\\
We test these parameters in MD simulations of 80-bead filaments with $\epsilon=40$, $\mu=2$, and volume fraction is $0.001$. Figure 3(b) illustrates the successful ``spun'' cross-links, but the system density is too low to form a system-spanning network. We achieve such a network in Fig. 3(c) by increasing chain length to 160 and volume fraction to 0.0018.
\\
We also verify the effects of changing the field parameters in the low-density system. Setting $\beta_0$ too high or too low creates a net bias in the chain bulk interactions towards planar or linear conformations, preventing gelation [Figs. 3(d)-3(e)]. If $\omega_2/\omega$ is too large or too small, the chain ends do not exhibit strong circular motion and chains do not spin together into a gel [Fig. 3(f)-3(g)]. Large values for $\omega$ and $\omega_2$ lead to weak torques on chain ends and no chain aggregation [Fig. 3(h)], while small values violate our assumption that we can time-average magnetic interactions in the bulk. In this case, chains form compact dynamic shapes as the entire chain responds to changes in the field [Fig. 3i].

\section{Conclusion}
In this paper, we presented analytic and computational studies of superparamagnetic filaments under the influence of a rapidly precessing field. Single filaments with their ends attached to loads can be used as microscopic actuators, with the strength of contraction determined by the angle of precession. Collections of free chains are capable of complex configurations. We found that these configurations are determined primarily by two effects: the desired orientation of beads in the chain bulk, which is fixed by the precession angle, and the residual torque on chain ends which vanishes as the inverse square of the precession frequency. We can exploit the combination of these effects to generate particular desired materials and behavior. Our work demonstrates the power of magnetic filaments as the basis for dynamic designer materials, in which desired properties can be created and destroyed in real time. We note that in the case of driven systems, metastable structures may depend on the initial configuration, a feature worth of exploring. Extending the present discussion with exact equilibrium configurations for filaments in multiple precession regimes is a subject for future research.

\section*{Acknowledgements}

This work was supported by the Center for Bio-Inspired Energy Science (CBES), which is an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0000989.

\begin{appendix}

\section{Helices in a precessing field} \label{sec:Helices}

Here we show that helices minimize the energy of the curve given by the elastic and magnetic contributions. Let $\mathbf{Y}(s)$ be the position vector of the curve parametrized by arc-length $s$ and of total length $L$, with given start and end points for the chain, $\mathbf{Y}(0) = \mathbf{Y}_0$ and $\mathbf{Y}(L)=\mathbf{Y}_L$, subject to the physical constraint that they are not more than a distance $L$ apart, i.e., $d = |\mathbf{Y}_L - \mathbf{Y}_0| < L$.
The Frenet-Serret (FS) frame adapted to the curve is composed by $\mathbf{T}(s)= \mathbf{Y}'(s)$ the tangent vector, $\mathbf{N}(s)$ the principal normal
and $\mathbf{B}(s) = \mathbf{T}(s) \times \mathbf{N}(s)$ the binormal. The FS formulas
describe how the frame rotates along the curve:
\begin{equation}
	\begin{pmatrix} \mathbf{T} ' \\ \mathbf{N} ' \\ \mathbf{B} ' \end{pmatrix} =
	\begin{pmatrix}
		0 & \kappa & 0 \\
		-\kappa & 0 & \tau \\
		0 & -\tau & 0
	\end{pmatrix}
	\begin{pmatrix} \mathbf{T} \\ \mathbf{N} \\ \mathbf{B} \end{pmatrix}\,,
\end{equation}
where $\kappa = \mathbf{T}' \cdot \mathbf{N}$ is the curvature and $\tau = \mathbf{N}' \cdot \mathbf{B}$ is the torsion. If $\kappa=0$, the curve is a straight line; if $\tau=0$, the curve lies in a plane.
\\
In order to minimize the total energy $H$ defined in Eq. (\ref{eq:continuum}), it is convenient to consider the effective energy, \cite{Auxiliary} :
\begin{equation} \label{eq:efcontinuumham}
	H_E = H + \int_0^L ds\, \left( \mathbf{F} \cdot (\mathbf{T}-\mathbf{\mathbf{Y}}') + \frac{\lambda}{2}(\mathbf{T}^2-1) \right) \,,
\end{equation}
where $\mathbf{F}$ is an auxiliary vector field acting as a Lagrange multiplier implementing the definition of the tangent vector as the arc-length derivative of the position vector \cite{Sphereconf, Envbias}, whereas $\lambda$ is another Lagrange multiplier implementing rigid bonds, thus constraining the distance between beads.
\\
The Euler-Lagrange equation for $\mathbf{Y}$ yields the conservation law $\mathbf{F}'=0$.
Using the FS equations and decomposing $\hat{\mathbf{z}}$ in the FS basis, the Euler-Lagrange equation for $\mathbf{T}$ yields $\mathbf{F}$ spanned in the FS basis
\begin{eqnarray}
	\mathbf{F} &=& -\left(\epsilon \kappa^2 + M T_z^2 +\lambda \right)  \mathbf{T} + \left(\epsilon \kappa' - M T_z N_z \right) \mathbf{N} \nonumber \\
	&&+	\left(\epsilon \kappa \tau - M T_z B_z \right) \mathbf{B} \,.
\end{eqnarray}
Differentiating $\mathbf{F}$, using again the FS formulas and projecting onto the tangent, we get
\begin{equation}
\mathbf{F}' \cdot \mathbf{T} =  -	\left(\epsilon \kappa^2 + M T_z^2 + \lambda \right)' - \kappa (\epsilon \kappa'- M T_zN_z)
	\label{eq:FpT}
\end{equation}
Recognizing that $\kappa N_z = T_z'$, we can write the second term as a total derivative--a consequence of the reparametrization invariance of the energy \cite{Hamforcurves}---allowing us to integrate and solve for $\lambda$:
\begin{equation}
	\lambda = -\frac{1}{2}(3 \epsilon \kappa^2 + M T_z^2) + c \,,
	\label{eq:lambda}
\end{equation}
where $c$ is a constant of integration, which can be regarded as a line tension controlling the total length.
We can now write
\begin{eqnarray}
\mathbf{F} &=& \left(\frac{\epsilon}{2} \kappa^2 + \frac{M}{2} T_z^2 - c \right)  \mathbf{T} + \left(\epsilon \kappa' - M T_z N_z \right) \mathbf{N} \nonumber \\
&&+	\left(\epsilon \kappa \tau - M T_z B_z \right) \mathbf{B}\,.
\end{eqnarray}
$-\mathbf{F}$ is the external force on the curve. The EL equations correspond to the projections onto the two normals
\begin{subequations} \label{eq:FpNB}
\begin{eqnarray}
\mathbf{F}' \cdot \mathbf{N} &=& \epsilon \left( \kappa'' + \kappa \left(\frac{\kappa^2}{2}-\tau^2\right) \right) \nonumber \\
&& + \kappa\left(M \left(\frac{T_z^2}{2}- N_z^2 \right) - c \right) =0 \,, \label{eq:FpN}\\
\mathbf{F}' \cdot \mathbf{B} &=&	\epsilon\left( \kappa \tau' + 2 \tau  \kappa' \right) - M \kappa N_z B_z =0\,. \label{eq:FpB}
\end{eqnarray}
\end{subequations}
Equation (\ref{eq:FpB}) can be satisfied when $\tau$ and $B_z$ vanish, which corresponds to configurations lying on a single plane containing the $z$ axis. A detailed analysis of these planar curves, determined by solving Eq. (\ref{eq:FpN}) under appropriate boundary conditions, will be discussed elsewhere.
\\
Since $\mathbf{F}$ is conserved, its magnitude is constant
\begin{eqnarray}
	F^2 &=& \left( \frac{\epsilon}{2} \kappa^2 + \frac{M}{2} T_z^2 - c\right)^2 + (\epsilon \kappa' - MT_zN_z)^2 \nonumber \\
	&&+(\epsilon \kappa \tau-MT_zB_z)^2\,.
	\label{eq:F2}
\end{eqnarray}
This provides a first integral of Eqs. (\ref{eq:FpNB}). The trivial way to satisfy this equation is to set all the quantities on the right hand side to be constant. Moreover, for constant $\kappa$, $\tau$, $T_z$ and $B_z$, Eq. (\ref{eq:FpB}) is satisfied if $N_z=0$, which corresponds to a helix winding around the $z$ axis. In this case Eq. (\ref{eq:FpN}) determines the constant $c$.
\\
In cylindrical coordinates $(\rho,\phi,z)$, the position vector of a $Z$-aligned helix centered on the origin with radius $\rho$ and pitch $p = 2 \pi \xi$\footnote{$\xi$ can be expressed as $\xi = \rho \tan \gamma$, with $\gamma$ the angle that the tangent makes with the azimuthal direction $\hat{\bm \phi}$, so it measures the slanting of the helix.} can be spanned with respect to the cylindrical basis $\{\hat{\bm \rho}, \hat{\bm \phi}, \hat{\bf z}\}$ as $\mathbf{Y}(s) =  \rho \hat{\bm \rho} + \xi \phi \hat{\bf z}$. The tangent vector is $\mathbf{T} = (\rho \hat{\bm \phi} + \xi \hat{\bf z})/\sqrt{\rho^2 + \xi^2}$. The FS curvature and torsion of the helix are $\kappa = \rho /(\rho^2 + \xi^2)$ and $\tau = \xi /(\rho^2 + \xi^2)$.
\\
To complete the characterization of the helical segment, besides $\rho$ and $\xi$, the total azimuthal angle $\Phi$ between the two ends has to be specified, so $\phi \in (0, \Phi]$. In this manner also the height difference is determined, $h= \xi \Phi$ and $z \in (0, h]$. The total length $L$ of the helix segment and distance $d$ between the two ends are given in terms of these three parameters by
\begin{equation} \label{eq:defLd}
L = \Phi\sqrt{\rho^2 + \xi^2}\,, \quad d = \sqrt{4 \rho^2 \sin \frac{\Phi}{2} + \xi^2 \Phi^2}
\end{equation}
Alternatively, the helical segment can be described in terms of $L$, $d$ and the angle $\eta$ that the displacement vector between the two ends makes with the helical axis. To this end we start by expressing $\rho$ and $\xi$ through the relations
\begin{equation}  \label{eq:xietarel}
 2 \rho \sin \frac{\Phi}{2} = d \sin \eta \,, \quad h = \xi \Phi = d \cos \eta  \,.
\end{equation}
Combining these expression with those for $L$ and $d$ in Eq. (\ref{eq:defLd}), we find that the total angle $\Phi$ is given in terms of $L$, $d$ and $\eta$ by the equation
\begin{equation}
	\mathrm{sinc}{\frac{\Phi}{2}} = \frac{\sin \eta}{\sqrt{\left(\frac{L}{d}\right)^2 - \cos^2{\eta}}}\,
	\label{eq:phiL}
\end{equation}
where $\mathrm{sinc} \, x = \sin x/x$. Since this is a transcendental equation it has to be solved numerically. Once $\Phi$ is known, $\rho$ and $\xi$ can be determined from Eqs. (\ref{eq:xietarel}).
\\
The bending and magnetic energies of the helical segment are given by
\begin{equation} \label{eq:EBEMhlx}
H_B =  \frac{\epsilon}{2} \frac{\Phi^2}{L}\left(1 - \left(\frac{h}{L}\right)^2\right)  \,, \quad
H_M = \frac{M(\beta)}{2}\, \frac{h^2}{L}\,.
\end{equation}
Clearly, the bending energy is minimized by reducing $\Phi$ and increasing $h$, so that the helix with one or fewer windings and higher pitch is preferred. By contrast the magnetic energy is minimized by reducing $h$, or equivalently the pitch of the helix, tending to form approximately circular loops on the plane orthogonal to the direction of precession. We see that the competition of the bending and magnetic energies can be exploited for controllable, repeatable actuation: the first tends to straighten the helix, whereas the latter tends to collapse it.
\\
In the special case that the helix completes one period $\Phi = 2\pi$, ($\eta=0$ and $h=d$) the total energy of the chain reduces to Eq. (\ref{eq:spring}).

\section{MD Simulation} \label{sec:mdsim}

We used the open source code LAMMPS to perform molecular dynamics simulations of beads. System units were chosen such that bead diameters, bead masses, the magnetic constant $\mu_0/(4\pi)$, and thermal energies were unitary. Beads in a field were modeled as point dipoles of fixed magnitude $2.0$, with hard sphere interactions given by truncated Lennard-Jones potentials, linear bonds with force constant $k=500$, and a Langevin thermostat with damping set to 0.1 so that the drag constant on each bead is ten times its mass. This ensures an overdamped behavior. Stiffness was set by a harmonic angular potential as described in the text. The external field was of magnitude $2000$, so its interaction strength with the dipoles was $4000 kT$.
\\
For the helices in Fig. 1(b) and the ESI movie, chains of length 40 beads were anchored to fixed points on the surface of load spheres with diameter 10 and mass 1000. We attempted both fixed and unfixed tangent boundary conditions. The fixed tangent boundary condition is enforced by applying the angular potential to the angle between the chain tangent at the surface of the load sphere and the radial vector from the center of the load sphere to the chain attachment point. Switching between boundary conditions did not observably alter the resulting dynamics. Figure 2(b) and the ESI movie use the fixed tangent condition.
\\
The conformations of Fig. 2 were generated with 100 chains of length 80 beads. The volume density for the system was 0.002 in system units.
\\
The conformations of Fig. 3 use the same system parameters as Fig. 2, except that the two cases the demonstrate twisting (Figs. \ref{Fig:3}(b) and \ref{Fig:3}(c)) were repeated with 400 chains of 120 beads to ensure the twist-linked gel remains intact at larger system sizes. In this case the density was 0.0035.
\\
Various conditions, including the spirals of Fig. 2(a) and the twisted gel of Fig. 3(b), were tested with lattice-Boltzmann hydrodynamics. Lattice spacing was set to bead size, and viscosity and density were set to 1 in system units. We did not observe any significant changes, but cannot guarantee that sufficiently strong hydrodynamic effects will not disrupt the results reported here.

\end{appendix}