mixing in nonlinear wave physics ( Robinson et al., 2001a; Whitham, 1974 ), and is seen in SSVEP studies using two stimulus frequencies (Regan, 1989 ). To brie ﬂy illustrate how mixing occurs in a general setting, consider a signal y(t)=sin( ω1t)+sin( ω2t) comprising two sinusoids with frequencies ω1andω2. If this signal is passed through a quadratic nonlinearity (i. e., is squared), the resulting signal (after using trigonometric identities) is y2(t)=1-cos(2 ω1t)/2-cos(2 ω2t)/ 2+cos[( ω1-ω2)t]-cos[(ω1+ω2)t], and thus contains second har-monics 2 ω1and 2 ω2, and sum and difference frequencies ω1+ω2and ω1-ω2, respectively. To have high power at frequencies f±requires a sharp alpha peak, so they are not observed in Fig. 1 (a) [except possibly near (0,10)Hz], where the alpha peak is broader and intersubject vari-ability tends to obscure this feature in the group average. Bicoherence (Kim and Powers, 1979 )i ss h o w ni n Fig. 3 (c) for f D=30Hz. Clear peaks at ( f1,f2)=( f D,fα) signify high phase coherence between f D,fα, and f D+fα. Peaks around ( f D,f D), (f D,0), and (0, f D) show phase coher-ence between f Dandf D±f D(2:1 harmonic generation; phase coherence with f1,2=0 represents the limiting case as f1,20). Enhancement is also seen on diagonals where f1+f2=f D,2f D,3f D, showing coupling be-tween frequencies that sum to f D,2f D,a n d3 f D. Entrainment involves more than just suppression of background peaks: enhancements A, B, and C in Fig. 3 (b) are precursors to phase-locked regions entered by increasing the amplitude Φn, as is explored further in Fig. 4. Enhancement A is an incipient 1:4 subharmonic; a small further increase in Φnis suf ﬁcient to enter this regime, yielding high power at f=f D/4, as seen by comparing Figs. 4 (c) and (d). Both B and C have two branches (those at B roughly have a forked shape on the end of the 1:2 line, but are blurred together here into a broad peak), and at each f Dthe sum of their (mean) frequencies equals that of the sharp peak immediately above in f(i. e., the frequencies sum to f Dat B and to f D/2 at C). Under changes in Φnthe branches at B always meet at the left endpoint of the subharmonic regime, shifting the back-ground activity away from f D. This is seen more clearly by comparing Figs. 4 (a)-(c), where increasing Φnleads to competition between the linear purely noise-driven fαresonance and a new resonance near onset of 1:2 entrainment. Such competition can result in incomplete en-trainment, where background activity is shifted in frequency to yield an enhancement at a new frequency different from f D; periodic pulling (Klinger et al., 1995 ) is a similar phenomenon seen in plasma physics. The noise-driven spectrum of the forced system results from interactionbetween background resonances (set by system parameters) and reso-nances due to the drive (set by both system and drive parameters). The system undergoes period doubling bifurcations as Φ nis in-creased for ﬁxed f Din intervals around f≈fα,2fα,....Fig. 3 (d) shows turning points (after transients) of sinusoidally driven time series at f D=20Hz vs. Φnin the absence of noise (i. e., σn=0). We stress that this analysis is of the interaction between the system and the periodic drive, not of the various frequency interactions possible with noise input that are better probed by spectra. For small Φnthe peak-to-peak amplitude (the maximum minus minimum shown) depends linearly onΦn; at larger Φnthe solution jumps to the period doubled branch(subharmonic entrainment). Successive period doubling with increas-ingΦnleads to aperiodic dynamics; further increases yield periodic so-lutions again. The LLE is positive for drive amplitudes in the range Φn≈5-9s-1(apart from narrow periodic windows), as shown in Fig. 3 (f), indicating that the dynamics are chaotic. Since the LLE can be estimated from time series without needing access to the full set of un-derlying dynamical variables, this model prediction should be testableusing experimental EEG time series. Another key feature is the existence of bistability between 1:1 and 1:2 entrainment, shown by coexistence of green and red points around Φ n≈1s-1. Thus both 1:1 and 1:2 entrainments (the latter always ex-hibits 1:1 activity concurrently) exist at the same point in parameter space. The implication is that experimental SSVEPs at given f DandΦn are not necessarily unique, yielding a new source of intertrial variability that warrants study in its own right. We also predict hysteresis in stud-ies where f Dis ramped slowly up and down through a zone where bis-table entrained states exist, which would give different spectra on the up and down paths. This should be relatively straightforward to test; one possible place to look is where the 1:1 and 1:2 entrainment zones meet in Fig. 1. Hysteresis effects would likely be obscured in studies where stimuli at the f Dof interest are presented in random order, as is commonly used to avoid habituation effects ( Regan, 1989 ). Spontane-ous transitions between bistable states during a long recording session might also contribute to intertrial variability. Extending the analysis of period doubling bifurcations in Fig. 3 (d) across a range of drive frequences reveals a rich pattern of entrainment zones. Fig. 3 (e) shows phase-locking regions vs. f DandΦnfor dynamics initially at a stable ﬁxed point. For comparison, Fig. 3 (d) lies at f D=20 Hz, and Fig. 3 (a) lies at Φn=2. 8s-1showing that 1:1 entrain-ment of fαoccurs near the 1:2 regime boundary. The 1:1 region is near-linear for small Φnbut increasingly nonlinear at larger Φn,f o re x-ample yielding regions of N:1 phase-locking at f D≈fα/N. The 1:1 regime bifurcates to 1: Nregimes over ranges of f Daround (and for minimal Φn at)f D≈Nfα, consistent with observed 1:3 and 1:8 entrainments (Herrmann, 2001 ) (in Figs. 5 and 4 therein, respectively). The 1: Nre-gions are always nonlinear, showing nested and overlapping period doublings toward dark blue regions where dynamics are variously of high period ( f D>80 f), quasiperiodic, or chaotic. Gaps in the 1:3 family of period doublings around f D=30Hz and Φn=14s-1reveal the locking regions underneath, hinting at a complicated structure reminis-cent of overlapping Arnold tongues ( Glass, 2001 ). The layout of entrain-ment zones in parameter space is also a new source of parameter constraints, because a given parameter set can be ruled out if it does not exhibit the correct entrainment properties. This is complementary to previous constraints in the linear regime, because satisfying these does not guarantee the same for the nonlinear regime. For example, the requirement that the system lie near a period doubling in the nonlinear dynamics is complementary to constraints previously eluci-dated for linear background EEG spectra ( Robinson et al., 2004 ). The structure in Fig. 3 (d) is for a driven system initially at a stable ﬁxed point, not a pre-existing limit cycle, and the entrainment in f D (Hz)f (Hz)a) 01 0 2 0 3 005101520 f D (Hz)f (Hz)b) 0 1 02 03 005101520 f D (Hz)f (Hz)c) 01 0 2 0 3 005101520 f D (Hz)f (Hz)d) 01 0 2 0 3 005101520 Fig. 4. Spectra vs. f Dfor increasing drive amplitude Φnwith other parameters as in Fig. 3 (a). (a) Φn=1s-1. (b)Φn=2s-1. (c)Φn=3s-1. (d)Φn=4s-1. Dark shading corre-sponds to high values. 1952 J. A. Roberts, P. A. Robinson / Neuro Image 62 (2012) 1947-1955