Stability of force-driven shear flows in nonequilibrium molecular simulations with periodic boundaries

Periodic boundary conditions (PBCs) are mainstays of molecular simulations,^1,2 where they are used to construct physically relevant models using material volumes that are often computationally restricted to be many orders of magnitude smaller than those in experiments. In standard PBCs, particles that exit the simulation volume reenter on the opposite face and interact with the nearest periodic images of the other particles. PBCs mitigate surface effects from the small volume, and simulated properties can agree well with experiments while modeling remarkably few particles.^3,4 In addition, PBCs play an especially important role in nonequilibrium molecular simulations of shear rheology.⁵ Although wall-driven shear flows mimicking experimental rheometry can be modeled, surface effects unduly influence these simulations because the gap height between the surfaces that can be modeled is restricted, leading to unrealistic confinement that overemphasizes surface effects and slip compared to experiments. To overcome this limitation, simple or oscillatory shear flow of an unbounded fluid can be mimicked using Lees–Edwards PBCs,⁶ where the images of the fully periodic simulation cell are effectively put in relative motion at a given strain rate. However, determining properties like the shear viscosity by this method remains challenging, particularly at small strain rates, because the fluctuating stress must be carefully sampled.

Motivated by these complications, alternative simulation methods that fix the stress rather than the strain rate have been developed. Here, we focus on one class of these approaches that includes Müller-Plathe’s reverse nonequilibrium simulation (RNES) method⁷ and the periodic Poiseuille flow method.⁸ These methods do not generate flow with forcing at the boundaries but rather impose a spatially varying perturbation (i.e., a distributed body force) on the “unbounded” fluid that is compatible with the PBCs. With thoughtful choice of the implementation and form of the perturbation, rheological properties like the shear viscosity can be extracted from the measured flow field, which usually has rapidly converging statistics.⁷ These methods are convenient to implement and have been widely used to simulate rheological properties of fluids, such as simple liquids,^9,10 water,¹¹ ionic liquids,^12,13 polymer solutions^14,15 and melts,^16–18 and colloidal dispersions.^19–23 However, we were recently surprised to find that the RNES method can fail to produce the desired flow field in certain simulation geometries.²⁴ RNES should generate a periodic parallel shear flow with two opposing Couette-like regions.⁷ We obtained this flow profile [Fig. 1(a)] in a cubic box, but undesired steady vortices [Fig. 1(b)] developed when the box was elongated in the flow direction. We found that the vortex formation depended on the presence of PBCs in the shear-gradient direction, the aspect ratio of the domain, and the shear rate. Moreover, the presence of vortices was not specific to the fluid model or using the RNES algorithm to create the flow. We speculated that this might indicate a general hydrodynamic instability underlying these nonequilibrium simulation techniques, but we were not able to establish a relationship for the conditions under which the instability occurred within the statistical accuracy of the simulations.

Here, we analyze the hydrodynamic stability of parallel viscous shear flows that are generated by nonequilibrium (force-driven) simulation methods within three-dimensional PBCs. Hydrodynamic stability is a long-studied topic,^25,26 but shear flows within a fully periodic domain^27–33 have received less attention than wall-bounded flows. This may be partially due to the fact that experimentally observed flows rarely, if ever, have PBCs; however, PBCs are the norm rather than the exception in molecular simulations. We first consider the general case of a parallel shear flow in PBCs that is driven by a distributed piecewise continuous body force. We derive an expression for the critical Reynolds number for the linear instability as a function of the aspect ratio of the fully periodic domain. We then analyze the stability of three specific flows useful for simulating shear rheology: sinusoidal (Kolmogorov) flow, periodic Poiseuille flow, and a flow like the one generated by RNES that we call periodic Couette flow. We perform complementary particle-based simulations that confirm the analysis and shed valuable insight on the emergence of the instability. For a fixed Reynolds number, the results indicate that many parallel shear flows in PBCs become linearly unstable due to a geometric effect of the periodic domain, which admits unstable disturbances once it becomes sufficiently long in the flow direction. This work provides a simple framework for selecting geometries with PBCs that stabilize a flow, requiring only the desired flow’s Fourier series expansion. It also underscores the need for caution when simulating dynamic processes with PBCs, which can differ in unexpected ways from domains bounded by surfaces because the PBCs do not impose the same constraints as the surfaces.

The rest of this article is organized as follows. We perform the stability analysis in Sec. II, and then, we describe technical details of the particle-based simulations that we used to validate it in Sec. III. We compare the analysis to the simulations in Sec. IV, showing excellent agreement between the two, before concluding in Sec. V.

In the simulation methods of interest, the fluid is modeled as discrete interacting particles within a simulation cell having PBCs in three dimensions. Shear flow is generated by imposing an external force that varies in one dimension on the individual particles. For example, the periodic Poiseuille method applies a constant-magnitude body force to all fluid particles but with a direction that depends on whether the particles are in the upper or lower half of the simulation cell,⁸ while RNES swaps momentum between particles in two slab regions to impose an effective average local force.⁷ In order to analyze the hydrodynamic stability of these particle-based methods, we formulate a continuum description of the force-driven flow that develops.

We consider the incompressible flow u of a Newtonian fluid governed by the standard continuity and Navier–Stokes equations,³⁴ 

with density ρ, pressure p, viscosity μ, and body force per volume f. The velocity field and stresses must be continuous through the periodic boundaries at x = ±L and y = ±H. The third dimension z is also periodic, but its details will not be required for our analysis. In keeping with common practice in force-driven simulations,^7,8 we also require that there is no net acceleration (translation) of the entire domain due to f, i.e., its average value is zero.

Following the standard linear stability analysis procedure,²⁵ we assume that the velocity u = U + δu can be separated into two parts: the steady (time-independent) base flow U and a perturbation δu. (The pressure field, p = P + δp, is expressed similarly.) Both U and u must satisfy the PBCs of the domain. In particular, we are interested in the stability of parallel shear flows U = (U_x(y), 0, 0) generated by steady body forces f = (f_x(y), 0, 0), so it suffices to consider the two-dimensional disturbance,³⁵ δu = e^ωt+ikxv(y) with v = (⁠$vx$⁠, $vy$⁠, 0) and δp = e^ωt+ikxq(y). The specific form of U_x and f_x will depend on the simulation method. Linearization of Eq. (2) ( Appendix A) yields

which has been written in dimensionless form by defining $ỹ$ = y/H and $Ũ$_x = U_x/U, where U is the maximum of U_x, which also naturally gives $ṽ$_y = $vy$/U, $k̃=Hk$⁠, and $ω̃=ωH/U$⁠. The Reynolds number is Re = UH/ν with ν = μ/ρ being the kinematic viscosity. The base flow is unstable when the real part of $ω̃$ is positive for a solution $ṽ$_y that satisfies the PBCs in $ỹ$⁠.

In all but a handful of cases, Eq. (3) must be solved numerically, e.g., by discretization or by expanding $ṽ$_y in a set of orthogonal basis functions.²⁶ Orszag solved Eq. (3) for plane Poiseuille flow using Chebyshev polynomials,³⁶ but other basis functions can be used. Due to the PBCs, it is advantageous to expand the flow fields in complex finite Fourier series that are inherently periodic, $ṽ$_y(⁠$ỹ$⁠) = ∑_n$ṽ$_ne^iπn$ỹ$ and $Ũ$_x(⁠$ỹ$⁠) = ∑_p$Ũ$_pe^iπp$ỹ$. These series can be directly substituted into Eq. (3) using termwise differentiation,

We then use the orthogonality of the basis functions to obtain

Note that the sum in Eq. (5), which resulted from products between $ṽ$_y, $Ũ$_x, and their derivatives in Eq. (3), is a convolution in Fourier space, as expected.

The coefficients $Ũ$_p are determined by the Fourier representation of a given base flow, so Eq. (5) can be written as an eigenvalue problem for $ṽ$_n by truncating the series to −N ≤ n, p ≤ N for sufficiently large N. This problem can be solved for different values of $k̃$ and Re to find modes having unstable $ω̃$⁠. Because of the PBC in x, $k̃$ is restricted in the values it can take by k = πm/L with m being an integer. In dimensionless variables, this gives $k̃=πm/α$ with α = L/H being the aspect ratio of the domain. This constraint introduces the aspect ratio into Eq. (5),

This approach can be considered a generalization of prior analysis for a sinusoidal flow²⁷ to an arbitrary (force-driven) parallel shear flow U_x(y) in a fully periodic domain where the PBC aspect ratio α constrains k. To find the limit of stability for a given α, we numerically determined the eigenvalues of Eq. (6)^36,37 as a function of Re and solved for the critical Reynolds number Re_c that gave the least stable $ω̃$ having real part $R(ω̃)=0$ when m = 1. (Disturbances with larger m are unstable at larger α because they appear as a ratio.) In doing so, we neglected the possibility of instabilities triggered by nonorthogonal eigenmodes.³⁸ 

In order to test the stability analysis [Eq. (6)], we simulated force-driven flows using the computationally efficient multiparticle collision dynamics method.^39–41 We used the stochastic rotation dynamics collision scheme³⁹ with cubic cells of edge length a, fixed 130° rotation angle,⁴² and random grid shifting⁴³; a Maxwell–Boltzmann rescaling thermostat to maintain constant temperature T;⁴⁴ and time 0.1 τ between collisions, where $τ=am/kBT$⁠, m is the particle mass, and k_B is Boltzmann’s constant. The particle density was ρ = 5 m/a³, giving a liquid-like Newtonian fluid with kinematic viscosity ν = 0.79 a²/τ.^45–47 All simulations were performed with hoomd-blue (version 2.6.0)^48–50 using a three-dimensional periodic simulation box with a square cross section (H = 50 a in y and z) and L varied in x to span 0.6 ≤ α ≤ 3. (The total number of particles ranged from 3 × 10⁶ to 15 × 10⁶.) The maximum velocity of the base shear flow was restricted to U ≲ 0.32 a/τ to avoid artificial effects at large Mach numbers,⁵¹ giving Re ≤ 20.

We considered three force-driven shear flows: a sinusoidal (Kolmogorov) flow,²⁷ a periodic Poiseuille flow,⁸ and a periodic Couette-like flow.^7,24 The first is a historically well-studied problem^27–31 that serves as a useful test of the simulations, while flows like the latter two are commonly used to simulate shear rheology. In our simulations, we applied a force per mass F_x(y) to each particle, giving f_x = ρF_x. (All particles had unit mass m, so F_x is also the force per particle.) Figure 2 shows F_x and U_x for the three flows when Re = 20, and details of their functional forms are given in  Appendix B. The periodic Poiseuille flow is the periodic extension of the classic Poiseuille flow, and it has two opposing parabolic regions in the periodic domain. The periodic Couette-like flow, which mimics the flow created by the RNES scheme, differs slightly from the true periodic extension of the classic Couette flow because it has two small parabolic regions that keep the shear stress continuous and F_x bounded, whereas the true periodic Couette flow would have a step change in the shear stress. For simplicity, we will refer to the simulated periodic Couette-like flow as periodic Couette flow.

The simulations were initialized by superimposing the base flow on the thermalized fluid and were run for 10⁵ τ to reach a steady state. The velocity field u(x, y) was then measured by averaging the particle velocities in square bins of size a² every 10 τ for 0.5 × 10⁵ τ. We used the measured u to assess whether the targeted base flow U_x was stable in the simulations. We adopted three empirical criteria to systematically classify the base shear flow as stable if:

1. The deviation of the average total kinetic energy in y from equipartition [⟨E_y⟩ = N_pk_BT/2 for N_p = (ρ/m)LH² particles] was less than 3 standard errors of the mean as computed from its own fluctuations,

2. the absolute deviation of the flow velocity from the base flow, |u_x(y) − U_x(y)|, was larger than 10⁻³ a/τ in fewer than 5% of the bins, and

3. the absolute value of the gradient velocity |u_y(x)| was larger than 10⁻³ a/τ in fewer than 5% of the bins.

If all three criteria were not met, the base flow was classified as unstable; all other cases were unclassified. Criterion (1) is a coarse check for the development of flows directed along y, which are not present in the base parallel shear flow, while criteria (2) and (3) are more detailed tests of the measured flow against the expected flow. The thresholds for criteria (2) and (3) were chosen to be tolerant of the numerical and statistical fluctuations in equilibrium simulations (no flow) and so base flows with small α that were visibly stable were properly classified. Criterion (3) tended to be very stringent and often failed even when criteria (1) and (2) held.

We first considered the stability of the sinusoidal shear flow, $Ũ$_x(⁠$ỹ$⁠) = sin(π$ỹ$) [Fig. 3(a)]. Its Fourier coefficients are $Ũ$_±1 = ∓i/2 and zero otherwise. Secondary flows like those shown in Fig. 1 readily formed in the simulations at sufficiently large Re and α. During the accessible simulation time, the flows remained stationary, although subsequent transitions to chaotic motion have been reported.^28,29 The bottom panel of Fig. 3(a) compares the predictions of the linear stability analysis with the simulation observations, which are in excellent agreement. Both indicate that the base flow remains stable for all α ≤ 1 over this range of Re, as proved theoretically by Meshalkin and Sinai.²⁷ However, Fig. 3(a) also shows that some flows in domains with α > 1 are also stable for sufficiently small Re, and over the range of α tested, all flows having Re ≲ 5 were stable. We confirmed that α and Re are the appropriate dimensionless groups to use for comparison between the simulations and the stability analysis by performing selected simulations with smaller boxes and/or higher fluid viscosity. The flows were completely consistent with the boundaries in Fig. 3(a) when interpreted using α and Re.

We subsequently tested two flows commonly used in nonequilibrium molecular simulations: the periodic extension of plane Poiseuille flow [Fig. 3(b)]⁸ and an approximate periodic extension of plane Couette flow [Fig. 3(c)] that mimics the RNES flow.⁷ Both were generated using constant forces in blocks of dimensionless half-width $d̃=d/H$⁠, where $d̃=1/2$ for the periodic Poiseuille flow and we chose $d̃=1/25$ to approximate the periodic Couette flow [see Fig. 2(a) and  Appendix B]. The corresponding Fourier coefficients are

Because the block forces are only piecewise continuous, $Ũ$_p ∼ 1/p³ and the Fourier series for d²$Ũ$_x/d $ỹ$² converges pointwise rather than uniformly. We accordingly took N = 100 and confirmed that using N = 200 did not significantly affect the computed critical Reynolds numbers.

Similar to the sinusoidal flow, the simulated periodic Poiseuille and Couette flows became unstable at sufficiently large Re and α, in excellent agreement with the stability analysis [Figs. 3(b) and 3(c)]. In our previous study using RNES⁷ to create a periodic Couette flow, we found stable shear flows for α ≲ 1.25;²⁴ this is in complete agreement with our new analysis for the studied parameters. We are unaware of prior reports of an instability in the periodic Poiseuille flow, but its stability was highly similar to the sinusoidal flow. This may not be surprising given that the two flows closely resemble each other [Figs. 3(a) and 3(b)]. The periodic Couette flow [Fig. 3(c)] was stable over a larger parameter space than the other two flows (see below).

Given the similar α dependence of Re_c for the three flows, we posited that the stability might be controlled by only the first (long wavelength) terms in the series expansions and applied a three-mode approximation (N = 1).^28,31 This drastic simplification allows the eigenvalues of Eq. (6) to be computed analytically,

$ω̃0$ and $ω̃+$ will always be negative because Re > 0, but $ω̃−$ can be zero or positive. We note that $Ũ$₁ and $Ũ$₋₁ are complex conjugates, |$Ũ$₁|² = $Ũ$₁$Ũ$₋₁, because $Ũ$_x is real valued. Solving for $ω̃−=0$ gives a relationship between the critical Reynolds number Re_c, α, and |$Ũ$₁|,

There is a minimum critical Reynolds number as α → ∞; this result is well-known for the sinusoidal flow,^27,28,31 where all flows having $Re<π2≈4.4$ are stable. However, Re_c increases as α decreases toward one, expanding the range of Re for which flows are stable. Both are in good agreement with the simulations, and they clarify a point of uncertainty in Ref. 24 about whether there is a flow rate (Reynolds number) below which vortices should not be observed in a simulation box of fixed size.

Equation (10) approximates the solution of Eq. (6) well for all three flows when α ≳ 1.5 (Fig. 3). This explains the similar stability of the periodic Poiseuille (|$Ũ$₁| = 16/π³ ≈ 0.52) and sinusoidal (|$Ũ$₁| = 0.5) flows, as their first Fourier modes are nearly identical, and the stability of the periodic Couette flow to larger Re (|$Ũ$₁| → 4/π² ≈ 0.41 as d → 0). However, there are some discrepancies between Eqs. (6) and (10) at smaller α. In fact, Eq. (6) predicts that the periodic Poiseuille flow becomes unstable for Re ≳ 50 when α = 1, but Eq. (10) diverges as α → 1. Certain shear flows may still be unstable in cubic domains when PBCs are used. Nevertheless, Eq. (10) provides a useful estimate for either (1) choosing a Reynolds number for a fixed aspect ratio or (2) choosing an aspect ratio for a fixed Reynolds number that keeps the base flow stable. It also highlights the generality of the instability, as Eq. (10) indicates similarly shaped flows (to lowest order) become unstable under similar conditions.

Most of the secondary flows that developed in our simulations had two stationary vortices [Fig. 1(b)], although at the largest aspect ratio and Reynolds numbers simulated—α = 3 and Re = 17.5 and 20—we found four stationary vortices for the sinusoidal and periodic Poiseuille flows. We hypothesized that the structure of some of the secondary flows might be connected to the most unstable eigenmode. This eigenmode, w = e^ikxv(y), can be reconstructed using the Fourier coefficients for $vy$ tht comprise the eigenvector having the largest $R(ω)$ and computing $vx$ = (i/k)d$vy$/dy based on the incompressibility of u ( Appendix A). Figure 4(a) shows w for the periodic Couette flow at Re = 15 and α = 1.4, which develops into the flow shown in Fig. 1(b). It has two pairs of counter-rotating vortices with streamlines primarily directed along the shear gradient y.

The final steady flow need not possess the same structure as either the base flow or the most unstable eigenmode. However, we considered as an ansatz that in some cases the simulated flow field might be well-approximated by a linear combination of the two, u ≈ c₁U + c₂w, subject to a shift of the coordinates with respect to the periodic boundaries. As an example, we performed a least-squares regression to the simulated u for the conditions in Fig. 4, determining optimal coefficients c₁ = 0.765 and c₂ = 0.373. The fitted flow field [Fig. 4(b)] bears striking similarity to the simulated flow [Fig. 1(b)], having a dimensionless root-mean-squared error of 0.02 per velocity component.

We also noted that the waiting time τ_w for the secondary flows to emerge in the simulations was shorter for points farther from the curve Re_c(α). This is qualitatively expected because these conditions are less stable based on their eigenvalues and should require smaller fluctuations (and less time) to depart from the base flow. To quantify this, we computed τ_w for the unstable flows by empirically fitting E_y to a hyperbolic tangent during the first 10⁵ τ simulated; we defined τ_w as the first time that E_y increased by 10% of the difference between its initial and final values ( Appendix C). These times can be compared to $1/R(ω)$⁠, which is the typical timescale associated with growth of the most unstable eigenmode. We found that τ_w had a strong linear correlation with $1/R(ω)$ for all three flows (Fig. 5). Moreover, we noted that the instability could take a surprisingly long time to emerge, up to nearly 1.5 × 10⁴ τ in dimensional units. Having this estimate of the timescale for instability is an added benefit of our analysis.

The stability of force-driven shear flows simulated in fully periodic domains seems to stand in stark contrast to those in bounded ones. The plane Couette flow is linearly stable for all infinitesimal disturbances,⁵² while the plane Poiseuille flow is linearly stable up to Re = 5772³⁶; in practice, both become unstable in the range of Re ≈ 10² to Re ≈ 10³.^53–55 However, their periodic extensions [Figs. 3(b) and 3(c)], driven by distributed body forces, are unstable at Re up to two orders of magnitude smaller. The change of stability in PBCs is not a simple consequence of the base flow; we previously showed that an unstable periodic Couette flow could be made stable in particle-based simulations by introducing a no-penetration boundary condition at either a half period or a full period of the flow.²⁴ Instead, the PBCs impose fundamentally different constraints: (1) the velocities and stresses are not obligated to take a certain value at a surface, expanding the types of flows that can be realized, and (2) the PBCs restrict the wavelengths of disturbances to those commensurate with the domain, introducing a strong geometric dependence once the domain admits unstable modes.

We have given a simple recipe for determining the stability of parallel shear flows in spatially periodic domains that requires only the Fourier series expansion of the flows, and we have applied it to understand the flows we observed in nonequilibrium molecular simulations. This recipe can be used to design well-behaved models and simulation methods and to choose appropriate simulation parameters; in the case studied, we encourage performing simulations in domains having α < 1. The presence of a hydrodynamic instability in nonequilibrium simulation methods such as RNES, which we and others had not previously appreciated, also highlights the need for caution when using PBCs to simulate certain dynamic processes, as the PBCs may fundamentally alter the underlying physics in unexpected ways.

The data that support the findings of this study are openly available in Princeton University’s DataSpace at http://arks.princeton.edu/ark:/88435/dsp01zw12z8154.⁵⁶ 

We thank Arash Nikoubashman and Zachary Sherman for helpful comments on this manuscript. M.P.H. and T.M.T. acknowledge support from the Welch Foundation (Grant No. F-1696). A.S. was supported by the Princeton Center for Complex Materials (PCCM), a U.S. National Science Foundation Materials Research Science and Engineering Center (Grant No. DMR-1420541). The simulations were performed using computational resources supported by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University.

Although the derivation of Eq. (3) is widely documented,^25,26 we include details here for completeness. Substituting u and p into Eqs. (1) and (2), using that U must also be a solution of the same, and neglecting terms of $O(δu2)$ yields the linearized equations,

Equation (A1) gives a simple relation between $vx$ and $vy$⁠,

where the prime denotes the ordinary derivative, $vy$′ = d$vy$/dy, while Eq. (A2) gives

Equation (A3) can be inserted in Eq. (A4) and solved for q,

Substituting Eq. (A6) into Eq. (A5) results in Eq. (3), which we note is the well-known Orr–Sommerfeld equation²⁶ after replacing ω by iω.

In this appendix, we give explicit functional forms for the studied periodic flow fields and the body forces that we applied to generate them. All quantities have units consistent with the multiparticle collision dynamics simulations (see Sec. III).

The sinusoidal flow,

was generated by a sinusoidal force applied to each particle based on its position,

We chose U to obtain a certain Re, setting the force amplitude F = (πν)²Re/H³.

Periodic extensions of the plane Poiseuille flow and the plane Couette flow were generated similarly using piecewise constant forces. This functional form ensured that velocities and stresses in the fluid were continuous. We defined two blocks centered at y = ±H/2 of half-width d ≤ H/2 each. Then, a positive, constant force F was applied to all particles in the upper block and a negative, constant force −F was applied to all particles in the lower block,

Linear momentum was conserved on average, but there were instantaneous fluctuations due to the distribution of particles. By solving Eq. (2) piecewise using the symmetries of the problem, the flow field can be deduced,

with

where g₁ is the usual linear (Couette-like) flow, and g₂ is the quadratic (Poiseuille-like) flow. The maximum velocity is

For the periodic plane Poiseuille flow,⁸ the width of the blocks was chosen to cover the full simulation box, d = H/2. Only g₂ contributes to this flow, establishing two opposing parabolic regions, and the force magnitude required for a given Reynolds number is F = 8ν²Re/H³.

By making d much smaller than H/2, the flow is dominated by g₁ and a periodic plane Couette-like flow can be generated having two linear regimes. In practice, d must be sufficiently large that there are enough particles in each block to reliably apply the force and achieve the targeted maximum velocity. The force magnitude required for a given Reynolds number is $F=8ν2Re/H3(1−(1−2d/H)2)$⁠. In the limit d → 0, the shear stress has a step change at y = ±H/2, and the applied forces are necessarily delta functions. For finite d, however, F is bounded, and the flow has two small quadratic regions that keep the shear stress continuous. We chose d = 2 a for our boxes having H = 50 a, which gave predominantly Couette-like flow fields but also ensured the blocks contained sufficient numbers of particles. The resulting flow is highly similar to that generated by the RNES method⁷ even though the origin of the perturbation is different.

In order to estimate the timescale for occurrence of the instability, we empirically fit the kinetic energy E_y to a hyperbolic tangent,

We fixed E₀ = 1/2 because the system was initially thermalized, and so its kinetic energy must obey equipartition. We then fit the final value E₁, the midpoint τ₀, and the width $w$⁠, resulting in fits like those shown in Fig. 6. The waiting time was defined as τ_w = τ₀ − 1.09861$w$ using the “10–90 thickness.” We found it challenging to reliably fit E_y for a handful of conditions close to the stability curves of Fig. 2; we accordingly neglected all points having $1/R(ω̃)>10$ in creating and analyzing Fig. 4.