SELF-TRAPPING OF DISKOSEISMIC CORRUGATION MODES IN NEUTRON STAR SPACETIMES

The stationary axisymmetric Kerr-spacetime can be characterized by only the central object mass, M, and the dimensionless spin, $j\equiv J/{M}^{2}$, where J is the angular momentum. The Kerr orbital and epicyclic frequencies for free particles near the equatorial plane are given by (Okazaki et al. 1987; Kato 1990),

$$\begin{eqnarray}&&{\rm{\Omega }}=\displaystyle \frac{1}{M({r}_{{\rm{BL}}}^{3/2}+j)},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\perp }={\rm{\Omega }}{\left(1-\displaystyle \frac{4j}{{r}_{{\rm{BL}}}^{3/2}}+\displaystyle \frac{3{j}^{2}}{{r}_{{\rm{BL}}}^{2}}\right)}^{1/2},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\kappa ={\rm{\Omega }}{\left(1-\displaystyle \frac{6}{{r}_{{\rm{BL}}}}+\displaystyle \frac{8j}{{r}_{{\rm{BL}}}^{3/2}}-\displaystyle \frac{3{j}^{2}}{{r}_{{\rm{BL}}}^{2}}\right)}^{1/2},\end{eqnarray} \tag{ 3 }$$

where Ω is the orbital frequency, ${{\rm{\Omega }}}_{\perp }$ and κ are the vertical and radial epicyclic frequencies and ${r}_{{\rm{BL}}}$ is the Boyer–Lindquist radial coordinate divided by the mass.

The inequality between the orbital and vertical-epicyclic frequencies results in Lense–Thirring nodal-precession (Bardeen & Petterson 1975), and allows low-frequency corrugation waves to propagate within thin accretion disks, in the region interior to the inner-vertical resonance (VR), where the wave frequency equals the Lense–Thirring frequency, ${{\rm{\Omega }}}_{{\rm{LT}}}\equiv {\rm{\Omega }}-{{\rm{\Omega }}}_{\perp }$ (see, e.g., Kato 1990; Silbergleit et al. 2001; Tsang & Lai 2009a).

Trapped corrugation waves (c-modes) in thin relativistic disks have frequencies lower than associated inertial-acoustic (p-modes) and inertial modes (g-modes). They have been proposed as an explanation for low-frequency variability for accreting black hole (BH) systems. However, in order for c-modes to exist, a reflecting inner boundary must be present due to, e.g., a strong pressure maximum (Kato 1990), or a truncation of the disk outside the innermost stable circular orbit (ISCO), which will quickly damp modes unless the density transition toward the disk sonic point is extremely sharp (Lai & Tsang 2009; Miranda et al. 2015).

Due to the theoretical uncertainty associated with such reflecting inner-disk boundaries, much of the work in relativistic diskoseismology has focussed on g-modes, which for (non-magnetized, see Fu & Lai 2009) thin BH-disks can become self-trapping due to the general-relativistic turnover of the radial epicyclic frequency, independent of the inner disk boundary condition (see Wagoner 1999; Kato 2001; Lai et al. 2013, for reviews).

In this letter we examine the trapping of disk c-modes in rotating neutron star (NS) spacetimes, that can significantly differ from Kerr in their higher-order multipole moments. We will show that for a certain range of NS parameters, the modification of the orbital and epicyclic frequencies due to the effect of a large spacetime quadrupole (Pappas 2012; Gondek-Rosińska et al. 2014; Pappas 2015; Wiśniewicz et al. 2015) can result in the self-trapping of c-modes due to general-relativistic effects in the disk exterior to the ISCO and the stellar surface. The effects of the quadrupole on geodesics have been of interest to several authors (see e.g., Morsink & Stella 1999; Bauböck et al. 2012; Pachón et al. 2012; Bauböck et al. 2015; Boshkayev et al. 2016). The observational identification of such a self-trapping c-mode in a spun-up accreting neutron star system through (e.g.,) phase-resolved iron-line spectroscopy (Tsang & Butsky 2013) could allow the NS quadrupole-moment to be constrained.

Near NSs additional parameters beyond the mass and spin are needed to characterize the spacetime metric (see e.g., Pappas 2009; Pappas & Apostolatos 2013). Here we expand the appropriate free-particle orbital/epicyclic frequencies to include up to the mass-quadrupole, M₂, and spin-octupole, S₃. We restrict our interest to the characteristic frequencies for particles on the equatorial plane, as the high-spin NSs almost certainly have angular momentum dominated by accretion-induced spin-up.

In a semi-Keplerian approach, these frequencies can be approximated as,

$$\begin{eqnarray}{\rm{\Omega }} & = & \displaystyle \frac{1}{M({r}_{{\rm{circ}}}^{3/2}+j)}\\ & & \times {\left(1-\displaystyle \frac{3\alpha {j}^{2}}{4{r}_{{\rm{circ}}}^{2}}+\displaystyle \frac{(1-4\alpha ){j}^{2}}{2{r}_{{\rm{circ}}}^{3}}-\displaystyle \frac{3(\alpha -2\beta ){j}^{3}}{4{r}_{{\rm{circ}}}^{7/2}}\right)}^{-1},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{{\rm{\Omega }}}_{\perp } & = & {\rm{\Omega }}\left(1-\displaystyle \frac{4j}{{r}_{{\rm{circ}}}^{3/2}}+\displaystyle \frac{3\alpha {j}^{2}}{{r}_{{\rm{circ}}}^{2}}-\displaystyle \frac{6(1-\alpha ){j}^{2}}{{r}_{{\rm{circ}}}^{3}}\right.\\ & & -{\left.\displaystyle \frac{3(4\beta -3\alpha ){j}^{3}}{{r}_{{\rm{circ}}}^{7/2}}\right)}^{1/2},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\kappa & = & {\rm{\Omega }}\left(1-\displaystyle \frac{6}{{r}_{{\rm{circ}}}}+\displaystyle \frac{8j}{{r}_{{\rm{circ}}}^{3/2}}-\displaystyle \frac{3\alpha {j}^{2}}{{r}_{{\rm{circ}}}^{2}}+\displaystyle \frac{3(4-5\alpha ){j}^{2}}{{r}_{{\rm{circ}}}^{3}}\right.\\ & & +{\left.\displaystyle \frac{6(3\beta -2\alpha ){j}^{3}}{{r}_{{\rm{circ}}}^{7/2}}\right)}^{1/2},\end{eqnarray} \tag{ 6 }$$

where ${r}_{{\rm{circ}}}\equiv \sqrt{{g}_{\phi \phi }}/M$ is the circumferential radius over the mass, $\alpha \equiv -{M}_{2}/({j}^{2}\;{M}^{3})$ is the quadrupole rotational-deformability (α-parameter), and $\beta \equiv -{S}_{3}/({j}^{3}\;{M}^{4})$ is the spin-octupole rotational-deformability (β-parameter). M₂ and S₃ are the relativistic multipole-moments (RMMs) (see Quevedo 1990, for review) and are calculated for NSs (Pappas & Apostolatos 2012, see). These frequencies are constructed for a stationary axisymmetric spacetime produced by a series expansion of an Ernst potential in terms of RMMs (similar to Ryan 1995, but resummed and expressed in terms of ${r}_{{\rm{circ}}}$).

These expressions have an exact expansion including up to the leading S₃ contribution, of ${ \mathcal O }({r}_{{\rm{circ}}}^{-5})$ for Ω, and ${ \mathcal O }({r}_{{\rm{circ}}}^{-13/2})$ for ${{\rm{\Omega }}}_{\perp }$ and κ. Expressions up to ${ \mathcal O }({M}_{2})$ are found by setting $\beta =0$. While ${ \mathcal O }({M}_{2})$ would be accurate enough to describe the orbital and the radial epicyclic frequency, ${{\rm{\Omega }}}_{{\rm{LT}}}$ is more sensitive because it is the small difference of two larger quantities. For NSs and Kerr-BHs ${ \mathcal O }({M}_{2})$ expressions underestimate the Lense–Thirring frequency, while ${ \mathcal O }({S}_{3})$ expressions overestimate it. Both capture the qualitative properties for $j\lesssim 0.8$, but a more accurate set of NS-frequencies could be constructed and will be the subject of a future work.

For $\alpha =\beta =1$, the expressions (4)–(6) yield the Kerr frequencies³ for a BH with spin parameter j. Alternatively, we can describe NSs by assigning "NS values" to the parameters α and β.

It was recently shown by Pappas & Apostolatos (2014), Stein et al. (2014), Yagi et al. (2014a, 2014b) that the RMM of NSs that obey realistic equations of state (EOSs) satisfy universal 3-hair relations between higher-order moments and the mass, the spin, and the quadrupole (the three hairs). Thus, the higher-order RMMs of NSs can be expressed in terms of the first three using those relations. For this study, the only higher-order moment of interest is S₃ (β-parameter). Using the relation between S₃ and M₂, β can be expressed in terms of α as (Pappas & Apostolatos 2014),

$$\begin{eqnarray}&&\beta ={(-0.36+1.48{(\sqrt{\alpha })}^{0.65})}^{3},\end{eqnarray} \tag{ 7 }$$

thus resulting in expressions for the NS frequencies that depend only on the parameters j, α, and NS mass M. For NSs, the spin parameter is constrained to be $j\lesssim 0.7$ (the Keplerian breakup limit, which is well within the range of validity of Equations (4)–(6)), while α ranges from ∼1.2 up to ∼10. In Section 5 we will discuss possible ranges for these parameters along with constraints based on NS spin and mass for realistic EOSs.

Following Okazaki et al. (1987), Tsang & Lai (2009a) we consider a thin isothermal accretion disk around an NS, utilizing isothermal Newtonian (flat space) expressions for the disk perturbations, while adopting the relativistic free-particle characteristic frequencies for the disk elements. While this mixed Newtonian/relativistic approach is formally inconsistent, it suffices to capture the qualitative effects of the spacetime on diskoseismic modes. A fully relativistic approach (see, e.g., Ipser & Lindblom 1992; Perez et al. 1997; Silbergleit et al. 2001) can be used in future work for more detailed quantitive study.

For a thin isothermal disk with unperturbed velocity ${{\boldsymbol{u}}}_{o}=(0,r{\rm{\Omega }},0)$ in cylindrical coordinates $(r,\phi ,z)$, the unperturbed vertical density profile can be given as (for small $z\ll r$)

$$\begin{eqnarray}&&{\rho }_{o}(r,z)=\displaystyle \frac{{\rm{\Sigma }}(r)}{\sqrt{2\pi }H}\mathrm{exp}(-{z}^{2}/2\;{H}^{2}).\end{eqnarray} \tag{ 8 }$$

Here ${\rm{\Sigma }}(r)$ is the (vertically integrated) surface density and $H={c}_{s}/{{\rm{\Omega }}}_{\perp }$ is the vertical scale height, where ${c}_{s}\ll r{\rm{\Omega }}$ is the isothermal sound speed.

Perturbing the mass and momentum conservation equations gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\delta \rho +{\rm{\nabla }}\cdot ({\rho }_{o}\delta {\boldsymbol{u}}+{{\boldsymbol{u}}}_{o}\delta \rho )=0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\delta {\boldsymbol{u}}+({{\boldsymbol{u}}}_{o}\cdot {\rm{\nabla }})\delta {\boldsymbol{u}}+(\delta {\boldsymbol{u}}\cdot {\rm{\nabla }}){{\boldsymbol{u}}}_{o}=-{\rm{\nabla }}\delta h,\end{eqnarray} \tag{ 10 }$$

with the enthalpy perturbation $\delta h\equiv \delta P/\rho ={c}_{s}^{2}\delta \rho /\rho$, where we assume that the perturbations are also isothermal. Following Tsang & Lai (2009a) (see also Zhang & Lai 2006) we assume perturbations of the form $\delta h,\delta P,\delta {\boldsymbol{u}},\delta \rho \propto \mathrm{exp}({im}\phi -i\omega t){H}_{n}(z/H)$, where ${H}_{n}(z/H)$ is the Hermite-polynomial of order n, m is the azimuthal wavenumber, and ω is the angular perturbation-frequency. After some manipulation while neglecting terms proportional to ${dH}/{dr}\sim O(1/r)$, the perturbation equations can be combined to give the wave equation for disk enthalpy perturbations,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}}{{{dr}}^{2}}\delta h-\left(\displaystyle \frac{d}{{dr}}\mathrm{ln}\displaystyle \frac{D}{r{\rm{\Sigma }}}\right)\displaystyle \frac{d}{{dr}}\delta h\\ &&\qquad +\left[\displaystyle \frac{2\;m{\rm{\Omega }}}{r\tilde{\omega }}\displaystyle \frac{d}{{dr}}\mathrm{ln}\displaystyle \frac{D}{{\rm{\Omega }}{\rm{\Sigma }}}-\displaystyle \frac{{m}^{2}}{{r}^{2}}-\displaystyle \frac{D({\tilde{\omega }}^{2}-n{{\rm{\Omega }}}_{\perp }^{2})}{{c}_{s}^{2}{\tilde{\omega }}^{2}}\right]\delta h=0,\end{eqnarray} \tag{ 11 }$$

where $D\equiv {\kappa }^{2}-{\tilde{\omega }}^{2}$, and $\tilde{\omega }\equiv \omega -m{\rm{\Omega }}$. Here, to capture the general-relativistic effects, we allow the characteristic frequencies to be given by Equations (4)–(6), with coordinate mapping ${r}_{{\rm{circ}}}=r/M$. This is more self-consistent but quantitatively different from the mapping adopted by, e.g., Kato (1990) or Tsang & Lai (2009a), which take ${r}_{{\rm{BL}}}=r/M$, though the results are qualitatively similar for the Kerr limit.

Far from the Lindblad resonances (where D = 0) and the corotation resonance (where $\tilde{\omega }=0$), we can adopt a WKB approach, assuming radial dependence $\delta h\propto \mathrm{exp}({ikr})$. This gives the WKB dispersion relation (Okazaki et al. 1987)

$$\begin{eqnarray}&&{c}_{s}^{2}{k}^{2}=\displaystyle \frac{({\kappa }^{2}-{\tilde{\omega }}^{2})(n{{\rm{\Omega }}}_{\perp }^{2}-{\tilde{\omega }}^{2})}{{\tilde{\omega }}^{2}},\end{eqnarray} \tag{ 12 }$$

where waves can propagate in regions where ${k}^{2}\gt 0$. Modes with no vertical structure, such that n = 0, are known as inertial-acoustic modes (p-modes) and propagate in the region where $\tilde{\omega }\gt \kappa$ (see e.g., Lai & Tsang 2009; Tsang & Lai 2009b). Waves with $n\gt 0$ can propagate where ${\tilde{\omega }}^{2}\lt {\kappa }^{2}\lt n{{\rm{\Omega }}}_{\perp }^{2}$, the inertial-mode (g-mode) propagation region (see, e.g., Wagoner 2012), or where ${\tilde{\omega }}^{2}\gt n{{\rm{\Omega }}}_{\perp }^{2}\gt {\kappa }^{2}$, the c-mode propagation region. In this work we primarily consider the fundamental c-mode (Kato 1990), with m = 1, n = 1, where the boundary of the propagation region is determined by the Lense–Thirring precession frequency $\omega ={\rm{\Omega }}-{{\rm{\Omega }}}_{\perp }$, though similar results will hold for any $m=\sqrt{n}$.

In order for c-modes in the Kerr metric to be trapped a reflecting inner boundary is required such that a mode cavity can form (Figure 1).

Zoom In Zoom Out Reset image size

For spacetimes with different RMM, the Lense–Thirring frequency can be significantly modified close to the central object. For some NS parameters this can result in a turnover of ${{\rm{\Omega }}}_{{\rm{LT}}}$, allowing c-modes to become self-trapping, reflected on both sides by the inner VR, where $\omega ={{\rm{\Omega }}}_{{\rm{LT}}}$, as shown in the Figure 2. For ${{\rm{\Omega }}}_{{\rm{LT}}}$ we have used Equations (4), (5) and not their $\beta =0$ counterparts, as this gives more conservative estimates for the NS parameter-space, since the peak appears for larger values of j and α.

Zoom In Zoom Out Reset image size

In Figure 3, we show the peak of the Lense–Thirring frequency, ${f}_{{\rm{peak}}}={{\rm{\Omega }}}_{{\rm{LT}}}({r}_{{\rm{peak}}})/2\pi$ (in units of ${\rm{Hz}}\times \;{M}_{\odot }/M$), for the spacetime parameters j and α. We also indicate the j and α where the radius at which ${{\rm{\Omega }}}_{{\rm{LT}}}$ peaks, ${r}_{{\rm{peak}}}$, is equal to the expected equatorial radius of the NS, ${r}_{{\rm{NS}}}$. For NSs with realistic EOSs, the equatorial radius can be approximated by a single function of j and α using the fitting formula⁴ from Appendix A of Pappas (2015).

Zoom In Zoom Out Reset image size

NSs with j and α values above the line where ${r}_{{\rm{peak}}}={r}_{{\rm{NS}}}$ have ${r}_{{\rm{peak}}}$ located outside the star, and thus can induce self-trapping of c-modes. The discrete c-mode spectra and detailed radial structure depend on the disk sound-speed, such that a half-integer number of wavelengths fit in the cavity at the mode frequency. However, all self-trapped c-modes should have a radial trapping region with inner edge distinct from the disk edge, and frequency below $m\times {f}_{{\rm{peak}}}$ (see Figure 2), set by $m\times {{\rm{\Omega }}}_{{\rm{LT}}}$ at both the inner and outer edges of the self-trapping region.

NSs with j and α values at or below this line (or NSs that have disks magnetically truncated outside ${r}_{{\rm{peak}}}$), have maximum value of ${{\rm{\Omega }}}_{{\rm{LT}}}$ occuring at ${r}_{{\rm{NS}}}$ (or the disk inner edge). In these systems c-modes must be trapped by reflection at the disk edge, with frequency set by $m\times {{\rm{\Omega }}}_{{\rm{LT}}}({r}_{{\rm{IVR}}})$ at the outer edge of the trapping region, where the spacetime quadrupole effects are more subtle.

In the previous section we determined the range of the j and α parameters (Figure 3) such that c-modes can become self-trapped. However, in order to convert these j and α ranges into ranges of the more familiar NS properties of spin frequency and mass, an EOS must be selected.

Here we examine NSs with four different standard EOSs: AU, SLy4, APR, and L, in order of increasing stiffness. We have chosen SLy4 and APR to represent mid-range realistic NS EOSs, while AU and L we use as proxies for the soft and stiff extremes for realistic EOSs. In Figure 4 we show NS mass and spin frequency contours as a function of j and α for each of these EOSs. The mass contours are determined using interpolation of numerically generated RNS models (Stergioulas & Friedman 1995). The models are the same used by Pappas & Apostolatos (2014) (see Supplemental material of that paper for details) and cover the full range of NS parameters for each EOS. Pappas (2015) showed in Appendix B that the quantity $M\times {f}_{{\rm{spin}}}$ can be given by a single fitting formula, for the range of realistic EOSs, as a function of j and α that is accurate to within 2% across this parameter range. The frequency contours of Figure 4 are determined using this fitting formula, and the interpolated masses.

Zoom In Zoom Out Reset image size

From these contours we see that to have self-trapping c-modes for mid-range EOSs (SLy4, and APR) we require fast spins (${f}_{{\rm{spin}}}\gtrsim 0.6$ kHz), and low masses $M\lesssim 1.2\;{M}_{\odot }$. Comparing AU and L we see that stiffer EOSs have higher α for a given mass, and higher j for a given spin frequency, such that extremely stiff EOSs, like L, allow for self-trapping to occur for a much wider range of masses and spin frequencies. Extremely soft EOSs, like AU, require unrealistically low-mass and high-spin in order for an NS to have self-trapping c-modes. In general, the possible range of NS parameters for which c-modes can be observed will be different for different EOSs, providing the possibility to distinguish between them.

We have examined the diskoseismic c-modes for thin accretion disks using the relativistic free-particle characteristic frequencies, Ω, ${{\rm{\Omega }}}_{\perp }$, and κ, in spacetimes with more prominent higher-order multipoles. Such spacetimes can arise due to spin-induced deformation of an NS contributing to the mass-quadrupole, spin-octupole, and higher-order RMMs. Expanding up to the lowest-order spin-octupole contribution we have determined the characteristic frequencies for these spacetimes in terms of the mass M, the spin parameter j, the quadrupole α-parameter, and spin-octupole β-parameter. Using these frequencies, and the disk perturbation equations, we have shown that c-modes can become self-trapped in NS spacetimes (Figure 2) due to a turnover of the Lense–Thirring frequency ${{\rm{\Omega }}}_{{\rm{LT}}}={\rm{\Omega }}-{{\rm{\Omega }}}_{\perp }$, for sufficiently high values of the parameters j and α (Figure 3).

The masses and spin-frequencies corresponding to these values of j and α depend on the stiffness of the EOS. In Figure 4 we showed the various mass and spin-frequency ranges above which c-modes can become self-trapping. Mid-range realistic EOSs, like APR or SLy4, require a combination of very high spin frequencies ($\gtrsim 0.6$ kHz), and low masses ($\lesssim 1.2\;{M}_{\odot }$) in order to reach this threshold. The fastest known millisecond pulsar, PSR J1748–2446ad, has spin frequency 716 Hz (Hessels et al. 2006), while the lowest well-constrained NS masses currently measured are $\sim 1.0\;{M}_{\odot }$–$1.1\;{M}_{\odot }$ (Lattimer 2012; Martinez et al. 2015). Thus NS spin and mass combinations sufficient to induce c-mode self-trapping for mid-range EOSs are plausible but likely very rare. This, of course, ignores evolutionary considerations, e.g., that recycled millisecond pulsars require significant accretion in order to be spun up, which may raise the lowest plausible mass by up to $\sim 0.1\;{M}_{\odot }$, depending on magnetic accretion torque history (Alpar et al. 1982; Tauris et al. 2012).

For stiffer EOSs, these mass and spin requirements are loosened. For the extremely stiff L EOS a $1.4\;{M}_{\odot }$ NS with spin period $\gtrsim 400$ Hz can have c-mode self-trapping. However, extremely soft EOSs, represented here by AU, require unrealistically high-spins and low-masses to induce c-mode self-trapping.

If high-spin, low-mass accreting NS candidates can be identified, it may be possible to detect self-trapping c-modes, e.g., through the use of sufficiently high-resolution phase-resolved spectroscopy of the Fe-Kα fluorescence line (Miller & Homan 2005; Tsang & Butsky 2013; Stevens et al. 2014; Ingram & van der Klis 2015). The broadened Fe-Kα line probes the structure of the inner part of the disk, such that the radial propagation region for the c-modes could be identified through the redshift range of the iron line variability at the c-mode frequency (Tsang & Butsky 2013). Such observations may be possible with new or upcoming X-ray missions such as LOFT (Feroci et al. 2012), Astro-H (Takahashi et al. 2008), or Astro-SAT (Singh et al. 2014).

The self-trapping of c-modes in NS accretion disk represents a qualitative difference in the behavior of the system due to the influence of the spacetime quadrupole. Identifying these self-trapped modes can then help to constrain the NS quadrupole, as well as the NS EOS through its rotational deformability.

Finally, we point out that similar effects could also be observed for more "exotic" compact objects. The turnover of ${{\rm{\Omega }}}_{{\rm{LT}}}$ would also be present for quark stars (Gondek-Rosińska et al. 2014), therefore self-trapping c-modes should be expected from these objects as well. Moreover, since quark stars can have smaller masses and radii than NSs as well as higher spin parameters (Gourgoulhon et al. 1999), the range of j and α where such quark stars could develop self-trapping c-modes should be significantly larger.

Compact objects in alternative theories of gravity could also be considered, as NSs and BHs in some of these theories, due to the presence of additional degrees of freedom, can develop additional "hair" that modify the spectrum of multipole moments and the geodesic frequencies (Berti et al. 2015; Pappas & Sotiriou 2015). Self-trapping c-modes might be present in these theories for a quite different range of parameters, distinguishing them from the GR predictions. For example, a rotating BH of an alternative theory that deviates from Kerr by having a larger quadrupole would have self-trapping c-modes in contrast to the usual Kerr BH.

D.T. was supported by the UMCP-Astronomy Center for Theory and Computation Prize Fellowship program. G.P. was supported by NSF CAREER Grant No. PHY-1055103. We thank L.C. Stein, E. Berti, C. Mingarelli, and M.C. Miller for useful discussion.