1 Introduction
Plasmas can sustain laser intensities many orders of magnitude larger than typical solid state optical components. This makes plasmas attractive for the manipulation and control of intense laser beams, but requires the formation of large and stable electron density structures in the plasma to impact the propagation of the laser light significantly. A typical approach to this problem is via ponderomotive forces using additional driving laser beams. Several applications of plasma photonics based on light scattering off electron density structures have been proposed, including short pulse amplification via the resonant excitation of plasma waves (Malkin, Shvets & Fisch Reference Malkin, Shvets and Fisch1999), transient plasma gratings (Lehmann & Spatschek Reference Lehmann and Spatschek2016), crossed-beam energy transfer for symmetry control in inertial confinement fusion (Michel et al. Reference Michel, Divol, Williams, Weber, Thomas, Callahan, Haan, Salmonson, Dixit and Hinkel2009) and, more recently, plasma-based polarization control (Michel et al. Reference Michel, Divol, Turnbull and Moody2014; Turnbull et al. Reference Turnbull, Michel, Chapman, Tubman, Pollock, Chen, Goyon, Ross, Divol and Woolsey2016, Reference Turnbull, Goyon, Kemp, Pollock, Mariscal, Divol, Ross, Patankar, Moody and Michel2017; Lehmann & Spatschek Reference Lehmann and Spatschek2018). All of these examples involve the formation of travelling plasma structures. Furthermore, efficiency was an important factor in these applications, as the goal was to create the largest amplitude plasma density perturbation using the smallest possible driver intensities. Recently, it was shown that very large amplitude travelling or standing ion acoustic waves (SIAW) can be excited via the autoresonance (AR) approach (Friedland & Shagalov Reference Friedland and Shagalov2014, Reference Friedland and Shagalov2017; Friedland et al. Reference Friedland, Marcus, Wurtele and Michel2019). The latter exploits the salient property of nonlinear waves to stay in resonance with driving perturbations if some parameter in the system (e.g. the driving frequency) varies in time (for a review of several AR applications see Fajans & Friedland (Reference Fajans and Friedland2001) and Friedland (Reference Friedland2009)). Such systems preserve the resonance with the drive via self-adjustment of the driven wave amplitude leading to a large amplitude amplification even when the driver is relatively weak. In this paper, we will discuss the AR approach to excitation of very large amplitude standing plasma waves (SPW) using both the kinetic model and a simplified water bag model. The nonlinear SIAWs (Friedland et al. Reference Friedland, Marcus, Wurtele and Michel2019) and SPWs (see a related analysis of nonlinear Trivelpiece–Gould waves in Dubin & Ashourvan (Reference Dubin and Ashourvan2015)) are much more complex than their travelling wave counterparts. For example, the autoresonant SIAW driven to large amplitudes by a weak, chirped frequency standing ponderomotive wave is a two-phase solution, each phase locked to one of the travelling waves comprising the drive (Friedland et al. Reference Friedland, Marcus, Wurtele and Michel2019). The peak electron density in these waves can reach several times the initial plasma density. We will show that similar giant plasma structures can be excited in the form of autoresonant SPWs.
The paper is organized as follows. In § 2, we will illustrate autoresonant SPWs in Vlasov–Poisson simulations and compare the results to those from the simplified water bag model presented in § 3. We will discuss weakly nonlinear SPWs in § 4 and use their form as an ansatz in developing a weakly nonlinear theory of driven–chirped SPWs based on the Whitham averaged variational principle (Whitham Reference Whitham1974) in § 5. In § 6, we reduce a system of slow coupled amplitude–phase mismatch equations yielding the autoresonance threshold on the driving amplitude for excitation of autoresonant SPWs. In the same section, we will also discuss the control of large amplitude autoresonant SPWs by tailoring the time dependence of the driving frequency. Finally, § 7 will present our conclusions.
2 Standing autoresonant plasma waves in Vlasov–Poison simulations
Our study of autoresonant SPWs is motivated by numerical simulations of the following one-dimensional Vlasov–Poisson (VP) system describing an externally driven plasma wave
$$\begin{eqnarray}f_{t}+uf_{x}+(\unicode[STIX]{x1D711}+\unicode[STIX]{x1D711}_{d})_{x}\,f_{u}=0,\quad \unicode[STIX]{x1D711}_{xx}=\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D711}+\int f\,\text{d}u-1.\end{eqnarray}$$ Here, we assume constant ion density, 
$f(u,x,t)$ and 
$\unicode[STIX]{x1D711}(x,t)$ are the electron velocity distribution and the electric potential, and 
$\unicode[STIX]{x1D711}_{d}=2\unicode[STIX]{x1D700}\cos (kx)\cos \unicode[STIX]{x1D703}_{d}$, where 
$\unicode[STIX]{x1D703}_{d}=kx-\int \unicode[STIX]{x1D714}_{d}\,\text{d}t$ is a small amplitude standing wave-like ponderomotive potential having a slowly varying frequency 
$\unicode[STIX]{x1D714}_{d}(t)$. All dependent and independent variables in (2.1) are dimensionless, such that the position 
$x$, time 
$t$ and velocity 
$u$ are rescaled with respect to the inverse wave vector 
$1/k$, the inverse plasma frequency 
$\unicode[STIX]{x1D714}_{p}^{-1}=\sqrt{m/(4\unicode[STIX]{x03C0}\text{e}^{2}n_{0})}$ (
$n_{0},m,e$ being the ion density, the electron mass and charge) and 
$\unicode[STIX]{x1D714}_{p}/k$. Then, in the dimensionless form of the driving phase, 
$k=1$ and 
$\unicode[STIX]{x1D714}_{d}$ is rescaled by 
$\unicode[STIX]{x1D714}_{p}$. The distribution function and the potentials in (2.1) are rescaled with respect to 
$kn_{0}/\unicode[STIX]{x1D714}_{p}$, and 
$m\unicode[STIX]{x1D714}_{p}^{2}/(ek^{2})$, respectively. We have also added an effective screening term 
$\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D711}$ in the Poisson equation modelling the typical finite radial extent of the driving potential. The term describes for example the radial extent of the driving wave formed by beating two laser beams, or the radial boundary conditions in applications involving driven Trivelpiece–Gould modes (Dubin & Ashourvan Reference Dubin and Ashourvan2015). The screening affects the dispersion of the driven waves, but most importantly, reduces the threshold driving amplitude for the autoresonant excitation (see § 6). We assume 
$2\unicode[STIX]{x03C0}$ periodicity in 
$x$ imposed by the drive and solve the time evolution problem, subject to simple initial equilibrium 
$\unicode[STIX]{x1D711}(x,0)=0$ and 
$f(u,x,0)=(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}^{2})^{-1/2}\exp (-u^{2}/2\unicode[STIX]{x1D70E}^{2})$, where 
$\unicode[STIX]{x1D70E}=k\unicode[STIX]{x1D706}_{D}$ (
$\unicode[STIX]{x1D706}_{D}=u_{e}/\unicode[STIX]{x1D714}_{p}$ being the Debye length and 
$u_{e}$ the initial electron thermal velocity). Note that 
$\unicode[STIX]{x1D70E}$, 
$\unicode[STIX]{x1D705}$ and the driving parameters fully define our rescaled, dimensionless problem. Finally, we are interested in a driving frequency of the order of the plasma frequency and, consequently, assume 
$\unicode[STIX]{x1D70E}\ll 1$ to avoid the kinetic Landau resonance (
$\unicode[STIX]{x1D714}_{d}/k\approx u$) initially.
Figure 1. The colour map of the waveforms of the electron density 
$n$ (a) and fluid velocity 
$v$ (b) in the final time interval 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}=0.125$ of the simulation. The slopes of the white lines are the phase velocities of the two travelling waves comprising the drive.
We have applied our VP code (Friedland, Khain & Shagalov Reference Friedland, Khain and Shagalov2006) to solve this problem numerically and show some results of the simulations in figures 1–4 for
$\unicode[STIX]{x1D70E}=0.1$, 
$\unicode[STIX]{x1D705}=0$ (the effect of the non-zero 
$\unicode[STIX]{x1D705}$ will be discussed later), and driving frequency (dimensionless) 
$\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{0}+\unicode[STIX]{x1D6FC}t$, where 
$\unicode[STIX]{x1D714}_{0}=\sqrt{(1+\unicode[STIX]{x1D705}^{2})^{-1}+3\unicode[STIX]{x1D70E}^{2}}$ is the normalized plasma wave frequency in the linearized warm fluid model (see (4.4)), while 
$\unicode[STIX]{x1D6FC}=0.00005$ and 
$\unicode[STIX]{x1D700}=0.0021$. We present the results of the simulations versus slow time 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6FC}^{1/2}t$ (a convenient representation in AR problems), as the system evolves between initial 
$\unicode[STIX]{x1D70F}_{0}=-5$ and final 
$\unicode[STIX]{x1D70F}_{1}=2.7$ times (the total dimensionless time in this simulation is 
$\unicode[STIX]{x0394}t=7.7/\unicode[STIX]{x1D6FC}^{1/2}=1089$, i.e. 
$173$ periods of plasma oscillations). Figure 1 shows the waveform of the electron density 
$n=\int f\,\text{d}u$ and fluid velocity 
$v=\int uf\,\text{d}u$ in a small time window 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}=0.125\,$ of the evolution just before reaching 
$\unicode[STIX]{x1D70F}_{1}$. One can see a very large amplitude SPW with its peak electron density reaching 
$2.5$ times the initial plasma density (unity in our dimensionless problem). Furthermore, as expected, for a given time, the solutions are 
$2\unicode[STIX]{x03C0}$ periodic in 
$x$. But there also exist two directions shown by white lines in the figure with slopes 
$\text{d}x/\text{d}\unicode[STIX]{x1D70F}=\pm \unicode[STIX]{x1D714}_{d}/\unicode[STIX]{x1D6FC}^{1/2}$ (these slopes are the phase velocities of the two traveling waves comprising the drive) along which the solutions are 
$\unicode[STIX]{x03C0}$ periodic. This suggests that 
$u,n$ are periodic functions of three arguments 
$\unicode[STIX]{x1D6E9}_{1},\unicode[STIX]{x1D6E9}_{2},x$, where 
$\unicode[STIX]{x1D6E9}_{1,2}=x\pm \int \unicode[STIX]{x1D714}_{d}t$ and the solutions are 
$\unicode[STIX]{x03C0}$ periodic in 
$\unicode[STIX]{x1D6E9}_{1,2}$ and 
$2\unicode[STIX]{x03C0}$ periodic in 
$x.$ Similar two-phase solutions, with each phase 
$\unicode[STIX]{x1D6E9}_{1,2}$ locked to one of the phases of the waves comprising the drive, were also observed in autoresonant SIAWs (Friedland et al. Reference Friedland, Marcus, Wurtele and Michel2019). The occurrence of multi-phase solutions is known in the theory of some partial differential equations (Novikov et al. Reference Novikov, Manakov, Pitaevskii and Zakharov1984).
Figure 2. The snapshots of the electron distribution function at two times 
$\unicode[STIX]{x1D70F}=2.69$ (a) and 
$2.7$ (b) during the last period of the wave oscillation in figure 1. The fluid velocity vanishes in (a), while the electron density reaches its maximal value of 
$2.5$. In contrast, in (b), the fluid velocity 
$|v|$ is at its maximum at two locations in 
$x$. The straight white lines are at the locations of the phase velocities 
$\pm \unicode[STIX]{x1D714}_{d}/k$.
Additional results of the simulations are presented in figure 2, showing the snapshots of the electron distribution function at two different times inside 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70F}$. Figure 2(a) corresponds to the time when the fluid velocity vanishes, while the electron density reaches its maximal value of 2.5, while at the time in figure 2(b) the fluid velocity 
$|v|$ is at its maximum at two locations in 
$x$. The white horizontal lines in figure 2 show the locations of the phase velocities 
$\pm \unicode[STIX]{x1D714}_{d}/k$ associated with the travelling waves comprising the drive. One can see in figure 2(b) that a small fraction of the electrons in the tail of the distribution reach the Landau resonance, indicating proximity to the kinetic wave breaking. However, we did not observe any significant effect of this resonance on the excited wave at this stage of the evolution.
Figure 3. The excitation of the autoresonant SPW for driving amplitude above the AR threshold. (a) The envelope of the maxima of the electron fluid velocity 
$v$ versus slow time 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6FC}^{1/2}t$ in VP simulations (blue solid line) are compared to water bag model simulations (dotted red line). The dashed line represents the phase velocity 
$\unicode[STIX]{x1D714}_{d}/k$ of the driving wave. (b) Shows the frequencies of the driving (straight line) and driven (dotted line) waves. The continuous AR frequency locking is seen starting 
$\unicode[STIX]{x1D70F}\approx -3$.
Figure 4. The numerical simulations as in figure 3, but with the driving amplitude below the AR threshold. (a) Illustrates saturation of the excited wave amplitude, while (b) shows that the frequency locking discontinues, leading to saturation shortly after passage through the linear resonance.
The evolution of the system leading to the final stage shown in figures 1 and 2 is illustrated in figure 3, where panel (a) presents the envelope of the maxima of fluid velocity 
$v$ during the evolution versus slow time. In the same figure, we also show a similar envelope (dashed line) obtained by using a simplified water bag model described in the next section. We observe an excellent agreement between the two models even for very large excitations. This agreement is important since the water bag model will be used below in calculating the threshold for AR excitation of SPWs. Also shown in figure 3(a) by the dashed line is the driving wave phase velocity 
$u_{d}=\unicode[STIX]{x1D714}_{d}/k=\unicode[STIX]{x1D714}_{d}$. One can see that at the final times of the excitation the fluid velocity approaches the phase velocity, indicating again the proximity to the Landau resonance. Finally, figure 3(b) shows the driving and driven wave frequencies (the latter is calculated from the time differences between successive peaks of the fluid velocity) versus 
$\unicode[STIX]{x1D70F}$. It illustrates the characteristic signature of autoresonant waves, i.e. the frequency (phase) locking between the driven and driving waves, which starts even prior to passage through the linear resonance 
$\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{0}$ (at 
$\unicode[STIX]{x1D70F}=0$) and continues to the fully nonlinear stage. We complete this section by figure 4, showing the results of the simulations similar to those in figure 3 (the same parameters and initial conditions), but for a smaller driving amplitude 
$\unicode[STIX]{x1D700}=0.0011$. One can see in the figure that the wave excitation saturates in this case (see panel a), as the phase locking between the driven and driving wave discontinues (see panel b) shortly after passage through the linear resonance. We find that the peak electron density in this case reaches 1.7 of the initial density, as compared to 2.5 in the AR case in figure 3. These results illustrate the existence of the characteristic AR threshold 
$\unicode[STIX]{x1D700}_{th}$ on the driving amplitude (
$\unicode[STIX]{x1D700}_{th}=0.0015$ for the parameters of the above simulations, see (6.8) below). The autoresonance threshold is a weakly nonlinear phenomenon studied in many applications (Friedland Reference Friedland2009). We have seen in figures 3 and 4 that a simplified water bag model can be used in describing the weakly nonlinear evolution of the driven SPWs instead of using the full VP system. The next four sections will present the theory of weakly nonlinear autoresonant SPWs based on this model.
Figure 5. The water bag model. The electron distribution is confined between two limiting trajectories 
$u_{1,2}$.
3 The water bag model of driven–chirped plasma waves
The water bag model (DePackh Reference DePackh1962; Feix, Hohl & Staton Reference Feix, Hohl and Staton1969; Hohl Reference Hohl1969; Berk, Nielsen & Roberts Reference Berk, Nielsen and Roberts1970) assumes that the electron distribution function remains constant, 
$f(u,x,t)=1/(2\unicode[STIX]{x1D6E5})$, between two limiting trajectories 
$u_{1,2}(x,t)$ in phase space and vanishes outside these trajectories (see figure 5). Then the electron density is 
$n(x,t)=\int f\,\text{d}u=(u_{1}-u_{2})/(2\unicode[STIX]{x1D6E5})$ and our kinetic SPW problem is governed by the following set of momentum and Poisson equations
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u_{1t}+u_{1}u_{1x}=(\unicode[STIX]{x1D711}+\unicode[STIX]{x1D711}_{d})_{x},\\ u_{2t}+u_{2}u_{2x}=(\unicode[STIX]{x1D711}+\unicode[STIX]{x1D711}_{d})_{x},\\ \displaystyle \unicode[STIX]{x1D711}_{xx}=\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D711}+(u_{1}-u_{2})/(2\unicode[STIX]{x1D6E5})-1,\end{array}\right\}\end{eqnarray}$$ where we use the same normalization for 
$x,t,f,u_{1,2},\unicode[STIX]{x1D711},$ and 
$\unicode[STIX]{x1D711}_{d}$ as in the previous section. We will be solving this system subject to 
$2\unicode[STIX]{x03C0}$ periodicity in 
$x$ for a trivial initial equilibrium 
$\unicode[STIX]{x1D711}=0$, 
$u_{1,2}=\pm \unicode[STIX]{x1D6E5}$. Note that if one also defines 
$v(x,t)=(u_{1}+u_{2})/2$, equations (3.1) yield
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}n_{t}+(vn)_{x}=0,\\ v_{t}+vv_{x}=(\unicode[STIX]{x1D711}+\unicode[STIX]{x1D711}_{d})_{x}-\unicode[STIX]{x1D6E5}^{2}nn_{x},\\ \unicode[STIX]{x1D711}_{xx}=\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D711}+n-1.\end{array}\right\}\end{eqnarray}$$ Therefore, the water bag model is isomorphic to the usual warm fluid limit of the driven plasma waves with the adiabatic electron pressure scaling 
$p\sim n^{3}$ and 
$\unicode[STIX]{x1D6E5}^{2}=3\unicode[STIX]{x1D70E}^{2}$. The results shown by the dotted red lines in figures 3 and 4 were obtained by solving (3.1) numerically and show excellent agreement between the VP and water bag simulations until approach to the Landau resonance. These water bag simulations used a code similar to that for SIAWs (Friedland & Shagalov Reference Friedland and Shagalov2017) and based on a standard spectral method (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988).
In analysing the AR threshold in the problem, we use a Lagrangian approach. We introduce new potentials 
$\unicode[STIX]{x1D713}_{1,2}$ via 
$u_{1,2}=\pm \unicode[STIX]{x1D6E5}+\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713}_{1,2}$ and rewrite (3.1) as
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}_{1tx}+\unicode[STIX]{x0394}\unicode[STIX]{x1D713}_{1xx}+\unicode[STIX]{x1D713}_{1x}\unicode[STIX]{x1D713}_{1xx}=(\unicode[STIX]{x1D711}+\unicode[STIX]{x1D711}_{d})_{x},\\ \unicode[STIX]{x1D713}_{2tx}-\unicode[STIX]{x0394}\unicode[STIX]{x1D713}_{2xx}+\unicode[STIX]{x1D713}_{2x}\unicode[STIX]{x1D713}_{2xx}=(\unicode[STIX]{x1D711}+\unicode[STIX]{x1D711}_{d})_{x},\\ \unicode[STIX]{x1D711}_{xx}=\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D711}+(\unicode[STIX]{x1D713}_{1x}-\unicode[STIX]{x1D713}_{2x})/2\unicode[STIX]{x1D6E5}.\end{array}\right\}\end{eqnarray}$$ This system can be derived from the variation principle 
$\unicode[STIX]{x1D6FF}(\int L\,\text{d}x\,\text{d}t)=0$, with the three-field Lagrangian density
$$\begin{eqnarray}\displaystyle L & = & \displaystyle {\displaystyle \frac{1}{2}}\unicode[STIX]{x1D711}_{x}^{2}+{\displaystyle \frac{\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D711}^{2}}{2}}-{\displaystyle \frac{\unicode[STIX]{x1D713}_{1x}^{2}+\unicode[STIX]{x1D713}_{2x}^{2}}{4}}-{\displaystyle \frac{\unicode[STIX]{x1D713}_{1x}\unicode[STIX]{x1D713}_{1t}-\unicode[STIX]{x1D713}_{2x}\unicode[STIX]{x1D713}_{2t}}{4\unicode[STIX]{x1D6E5}}}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{\unicode[STIX]{x1D713}_{1x}^{3}-\unicode[STIX]{x1D713}_{2x}^{3}}{12\unicode[STIX]{x1D6E5}}}+{\displaystyle \frac{(\unicode[STIX]{x1D713}_{1x}-\unicode[STIX]{x1D713}_{2x})(\unicode[STIX]{x1D711}+\unicode[STIX]{x1D711}_{d})}{2\unicode[STIX]{x1D6E5}}}.\end{eqnarray}$$Our next goal is to apply the Whitham averaged variational principle (Whitham Reference Whitham1974) and derive the weakly nonlinear slow evolution system describing driven–chirped SPWs governed by this Lagrangian. We proceed by discussing weakly nonlinear driven and phase locked, but not chirped SPWs.
4 Weakly nonlinear SPWs
If one starts in the trivial equilibrium 
$\unicode[STIX]{x1D711}=0$, 
$n=1$, 
$v=0$ (
$u_{1,2}=\pm \unicode[STIX]{x1D6E5}$), equations (3.2) after averaging over one spatial period, yield constant-in-time averaged density and fluid velocity 
$\langle n\rangle =1$, 
$\langle v\rangle =0$ (
$\langle u_{1,2}\rangle =\pm \unicode[STIX]{x1D6E5}$). We consider the linear stage of driven SPWs first, i.e. write 
$n=1+\unicode[STIX]{x1D6FF}n$ and linearize (3.2) to get
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}(\unicode[STIX]{x1D6FF}n)_{t}+v_{x}=0,\\ v_{t}=\unicode[STIX]{x1D711}_{x}-\unicode[STIX]{x1D6E5}^{2}\unicode[STIX]{x1D6FF}n_{x}-2\unicode[STIX]{x1D700}\cos \unicode[STIX]{x1D703}_{d}\sin x,\\ \unicode[STIX]{x1D711}_{xx}=\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D711}+\unicode[STIX]{x1D6FF}n.\end{array}\right\}\end{eqnarray}$$ In the case of a constant driving frequency, this set yields phase-locked standing wave solutions of frequency 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{d}$ for all dependent variables
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D6FF}n=a\cos \unicode[STIX]{x1D703}\cos x,\\ v=b\sin \unicode[STIX]{x1D703}\sin x,\\ \unicode[STIX]{x1D711}=c\cos \unicode[STIX]{x1D703}\cos x,\end{array}\right\}\end{eqnarray}$$ where 
$\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{d}=\unicode[STIX]{x1D714}t$,
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}a=-{\displaystyle \frac{2\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D714}^{2}-\unicode[STIX]{x1D714}_{0}^{2}}},\\ b=\unicode[STIX]{x1D714}a,\\ c=-{\displaystyle \frac{a}{1+\unicode[STIX]{x1D705}^{2}}},\end{array}\right\}\end{eqnarray}$$and
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}^{2}={\displaystyle \frac{1}{(1+\unicode[STIX]{x1D705}^{2})}}+\unicode[STIX]{x1D6E5}^{2}.\end{eqnarray}$$ Note that 
$\unicode[STIX]{x1D714}_{0}$ is the natural frequency of a linear SPW in the problem and the case 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{0}$ corresponds to the linear resonance. As usual, this singularity is removed if one adds nonlinearity or dissipation. Equations (4.2) yield the linear solutions
$$\begin{eqnarray}u_{1,2}=v\pm \unicode[STIX]{x1D6E5}(1+\unicode[STIX]{x1D6FF}n)=\pm \unicode[STIX]{x1D6E5}+b\sin \unicode[STIX]{x1D703}\sin x\pm \unicode[STIX]{x0394}a\cos \unicode[STIX]{x1D703}\cos x\end{eqnarray}$$and, thus,
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{1,2}=-b\sin \unicode[STIX]{x1D703}\cos x\pm \unicode[STIX]{x0394}a\cos \unicode[STIX]{x1D703}\sin x.\end{eqnarray}$$ Interestingly, in contrast to 
$n,u$ and 
$\unicode[STIX]{x1D711}$, the linear solutions for 
$u_{1,2}$ and 
$\unicode[STIX]{x1D713}_{1,2}$ are not standing waves.
Our next goal is to include a weak nonlinearity in the problem, but still for a constant driving frequency case. To this end, we use spatial periodicity and write truncated Fourier expansions
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{1,2}=b_{1,2}\cos x+a_{1,2}\sin x+d_{1,2}\cos (2x)+e_{1,2}\sin (2x), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}=c_{1}\cos x+c_{2}\cos (2x), & \displaystyle\end{eqnarray}$$ where 
$a_{1,2}$, 
$b_{1,2}$ and 
$c_{1}$ are time-dependent amplitudes viewed as small first-order objects, while 
$d_{1,2}$, 
$e_{1,2}$ and 
$c_{2}$ are due to the nonlinearity and, thus, are of second order. The time dependence of the first-order amplitudes is assumed to be that of the linear solutions (4.6) and (4.2), i.e. 
$a_{1,2}=A_{1,2}\cos \unicode[STIX]{x1D703}$, 
$b_{1,2}=B_{1,2}\sin \unicode[STIX]{x1D703}$ and 
$c_{1}=C_{1}\cos \unicode[STIX]{x1D703}$. But what is the time dependence of the second-order amplitudes? Instead of working with the original system (3.1) for answering this question, we use a simpler Lagrangian approach in appendix A. The final result is (see (A 3))
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}d_{1,2}=D_{1,2}\sin (2\unicode[STIX]{x1D703}),\\ e_{1,2}=F_{1,2}+E_{1,2}\cos (2\unicode[STIX]{x1D703}),\\ c_{2}=B+C\cos (2\unicode[STIX]{x1D703}),\end{array}\right\}\end{eqnarray}$$ where 
$D_{1,2}$, 
$F_{1,2}$, 
$E_{1,2}$, 
$B$ and 
$C$ are constants. Thus, the potentials characterizing a weakly nonlinear SPW phase locked to the drive are
$$\begin{eqnarray}\displaystyle \hspace{-24.0pt} & \displaystyle \unicode[STIX]{x1D713}_{i}=B_{i}\sin \unicode[STIX]{x1D703}\cos x+A_{i}\cos \unicode[STIX]{x1D703}\sin x+D_{i}\sin (2\unicode[STIX]{x1D703})\cos (2x)+[F_{i}+E_{i}\cos (2\unicode[STIX]{x1D703})]\sin (2x), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \hspace{-24.0pt} & \displaystyle \unicode[STIX]{x1D711}=C_{1}\cos \unicode[STIX]{x1D703}\cos x+[B+C\cos (2\unicode[STIX]{x1D703})]\cos (2x). & \displaystyle\end{eqnarray}$$ Here 
$i=1,2$, the solutions preserve the space–time symmetry of the linear limit, and all the coefficients (amplitudes) are constant and could be related to the driving amplitude 
$\unicode[STIX]{x1D700}$ by using the approach of appendix A. Nevertheless, we will not follow this route because, knowing the form of (4.10) and (4.11) is sufficient for addressing our original driven–chirped problem via the Whitham averaged variational principle (Whitham Reference Whitham1974), as described next.
5 Whitham’s averaged variational principle for driven SPWs
As in the case of the standing IAWs (Friedland et al. Reference Friedland, Marcus, Wurtele and Michel2019), the Whitham approach uses the ansatz of form (4.10) and (4.11) for the solution of our driven–chirped problem, but now all the amplitudes are viewed as slow functions of time and 
$\unicode[STIX]{x1D703}$ is the fast wave phase having a slowly varying frequency 
$\unicode[STIX]{x1D714}(t)=\text{d}\unicode[STIX]{x1D703}/\text{d}t$ generally different from the frequency 
$\unicode[STIX]{x1D714}_{d}(t)$ of the chirped driving wave. This ansatz is then substituted into the Lagrangian density (3.4), where the driver is written as 
$\unicode[STIX]{x1D711}_{d}=2\unicode[STIX]{x1D700}\cos (\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6F7})\sin x$ and the phase mismatch 
$\unicode[STIX]{x1D6F7}(t)=\unicode[STIX]{x1D703}_{d}-\unicode[STIX]{x1D703}$ is viewed as a slow function of time. Then, all the slow variables are frozen in time and the Lagrangian density is averaged over 
$2\unicode[STIX]{x03C0}$ in both 
$x$ and 
$\unicode[STIX]{x1D703}$ (fast scales). This results in a new averaged Lagrangian density 
$\unicode[STIX]{x1D6EC}$, which depends on slow objects only (i.e. 
$A_{1,2}$, 
$B_{1,2}$, 
$C_{1}$ (first-order amplitudes), 
$D_{1,2}$, 
$E_{1,2}$, 
$F_{1,2}$, 
$B$, 
$C$ (second-order amplitudes), the wave frequency 
$\unicode[STIX]{x1D714}$ and the phase mismatch 
$\unicode[STIX]{x1D6F7}$). To fourth order in amplitudes we obtain (via Mathematica (Wolfram Research Inc. 2017)) 
$\unicode[STIX]{x1D6EC}=\unicode[STIX]{x1D6EC}_{2}+\unicode[STIX]{x1D6EC}_{4}+\unicode[STIX]{x1D6EC}_{d}$, where
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}_{2} & = & \displaystyle 2C_{1}^{2}(1+\unicode[STIX]{x1D705}^{2})-A_{1}^{2}-A_{2}^{2}-B_{1}^{2}-B_{2}^{2}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{2}{\unicode[STIX]{x1D6E5}}}[\unicode[STIX]{x1D714}(A_{1}B_{1}-A_{2}B_{2})-C_{1}(A_{1}-A_{2})],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}_{4} & = & \displaystyle -4(2F_{1}^{2}+2F_{2}^{2}+D_{1}^{2}+D_{2}^{2}+E_{1}^{2}+E_{2}^{2})\nonumber\\ \displaystyle & & \displaystyle +\,2(4+\unicode[STIX]{x1D705}^{2})(2B^{2}+C^{2})-{\displaystyle \frac{4}{\unicode[STIX]{x1D6E5}}} [2\unicode[STIX]{x1D714}(D_{1}E_{1}-D_{2}E_{2})\nonumber\\ \displaystyle & & \displaystyle +\,2B(F_{1}-F_{2})+C(E_{1}-E_{2})]-{\displaystyle \frac{1}{\unicode[STIX]{x1D6E5}}}[A_{1}B_{1}D_{1}-A_{2}B_{2}D_{2}]\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{1}{2\unicode[STIX]{x1D6E5}}}[B_{1}^{2}(2F_{1}-E_{1})-A_{1}^{2}(2F_{1}+E_{1})-B_{2}^{2}(2F_{2}-E_{2})+A_{2}^{2}(2F_{2}+E_{2})],\quad\end{eqnarray}$$
$$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{d}={\displaystyle \frac{4\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D6E5}}}(A_{1}-A_{2})\cos \unicode[STIX]{x1D6F7}.\quad\end{eqnarray}$$ The truncation of the averaged Lagrangian at fourth order in amplitudes is necessary for getting the usual lowest-order significant set of wave amplitude equations of weakly nonlinear resonance and second-order nonlinear shift of the wave frequency (see § 6). Note that 
$\unicode[STIX]{x1D6EC}$ is algebraic in terms of all amplitudes with the wave phase 
$\unicode[STIX]{x1D703}$ entering via 
$\unicode[STIX]{x1D6F7}=$
$\unicode[STIX]{x1D703}_{d}-\unicode[STIX]{x1D703}$ and 
$\unicode[STIX]{x1D714}(t)=\text{d}\unicode[STIX]{x1D703}/\text{d}t$. Therefore, Lagrange equations for the amplitudes comprise a set of eight algebraic equations
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6EC}}{\unicode[STIX]{x2202}Z_{i}}}=0,i=1,\ldots ,8,\end{eqnarray}$$ where the 
$Z_{i}$ represent each of the 8 slow amplitudes in the problem. In addition, the Lagrange equation for 
$\unicode[STIX]{x1D703}$ yields an ordinary differential equation (ODE)
$$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6EC}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}\right)=-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6EC}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}}}.\end{eqnarray}$$To lowest order, the last equation becomes
$$\begin{eqnarray}{\displaystyle \frac{\text{d}}{\text{d}t}}(A_{1}B_{1}-A_{2}B_{2})=-2\unicode[STIX]{x1D700}(A_{1}-A_{2})\sin \unicode[STIX]{x1D6F7},\end{eqnarray}$$ and, as expected, describes the slow evolution, provided 
$\unicode[STIX]{x1D700}$ is sufficiently small.
The plan for analysing our slow driven system is as follows. First, we will use the algebraic equations (5.4) for expressing seven of the eight amplitudes and 
$\unicode[STIX]{x1D714}$ in terms of 
$\unicode[STIX]{x1D6F7}$ and the eighth amplitude (chosen to be 
$C_{1}$), i.e. we obtain relations 
$Z_{i}=Z_{i}(C_{1},\unicode[STIX]{x1D6F7})$ and 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}(C_{1},\unicode[STIX]{x1D6F7})$. Then
$$\begin{eqnarray}{\displaystyle \frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}t}}=\unicode[STIX]{x1D714}_{d}(t)-\unicode[STIX]{x1D714}(C_{1},\unicode[STIX]{x1D6F7}),\end{eqnarray}$$which in combination with (5.6) comprise a closed set of two ODEs describing the slow evolution in our driven system. This plan is algebraically complex, but can be performed using Mathematica (Wolfram Research Inc. 2017). We describe the intermediate steps of this calculation in appendix B and here present the final results. To lowest significant order, the amplitudes in (5.6) are as in the linear undriven problem (see (B 2), (B 3))
$$\begin{eqnarray}\displaystyle & \displaystyle A_{1}=-A_{2}\approx -{\displaystyle \frac{C_{1}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2}}}=-C_{1}\unicode[STIX]{x1D6E5}(1+\unicode[STIX]{x1D705}^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle B_{1}=B_{2}\approx {\displaystyle \frac{C_{1}\unicode[STIX]{x1D714}_{0}}{\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2}}}=C_{1}\unicode[STIX]{x1D714}_{0}(1+\unicode[STIX]{x1D705}^{2}), & \displaystyle\end{eqnarray}$$while (see (B 9) and (B 10) in appendix B)
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}(C_{1},\unicode[STIX]{x1D6F7})\approx \unicode[STIX]{x1D714}_{0}^{2}+NC_{1}^{2}+{\displaystyle \frac{2\unicode[STIX]{x1D700}\cos \unicode[STIX]{x1D6F7}}{(1+\unicode[STIX]{x1D705}^{2})C_{1}}},\end{eqnarray}$$where
$$\begin{eqnarray}N={\displaystyle \frac{9\unicode[STIX]{x1D705}^{2}(3+\unicode[STIX]{x1D705}^{2})+R}{96[1+(4+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}^{2}]}},\end{eqnarray}$$ and 
$R=(1+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}^{2}[3(60+56\unicode[STIX]{x1D705}^{2}+11\unicode[STIX]{x1D705}^{4})+20(1+\unicode[STIX]{x1D705}^{2})(4+\unicode[STIX]{x1D705}^{2})(5+2\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}^{2}+16(1+\unicode[STIX]{x1D705}^{2})^{2}(4+\unicode[STIX]{x1D705}^{2})^{2}\unicode[STIX]{x1D6E5}^{4}]$.
6 Dynamics and control of chirped–driven SPWs
At this stage, we discuss the slow evolution in our driven–chirped SPWs system. Upon substitution of (5.8) and (5.9) into (5.6), we obtain
$$\begin{eqnarray}{\displaystyle \frac{\text{d}C_{1}}{\text{d}t}}=-{\displaystyle \frac{\unicode[STIX]{x1D700}}{\unicode[STIX]{x1D714}_{0}(1+\unicode[STIX]{x1D705}^{2})}}\sin \unicode[STIX]{x1D6F7}.\end{eqnarray}$$ Next, for a small frequency deviation in the vicinity of the linear resonance, 
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{0}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}$, equation (5.10) yields
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}\approx {\displaystyle \frac{N}{2\unicode[STIX]{x1D714}_{0}}}C_{1}^{2}+{\displaystyle \frac{\unicode[STIX]{x1D700}\cos \unicode[STIX]{x1D6F7}}{\unicode[STIX]{x1D714}_{0}(1+\unicode[STIX]{x1D705}^{2})C_{1}}},\end{eqnarray}$$ where the two terms on the right represent the wave frequency shifts due to the nonlinearity and interaction, respectively. Then, from (5.7) and using the driving frequency 
$\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{0}+\unicode[STIX]{x1D6FC}t$, we obtain
$$\begin{eqnarray}{\displaystyle \frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}t}}=\unicode[STIX]{x1D6FC}t-{\displaystyle \frac{N}{2\unicode[STIX]{x1D714}_{0}}}C_{1}^{2}-{\displaystyle \frac{\unicode[STIX]{x1D700}\cos \unicode[STIX]{x1D6F7}}{\unicode[STIX]{x1D714}_{0}(1+\unicode[STIX]{x1D705}^{2})C_{1}}}.\end{eqnarray}$$
Figure 6. The rescaled threshold driving amplitude 
$\unicode[STIX]{x1D700}_{th}/\unicode[STIX]{x1D6FC}^{3/4}$ versus 
$\unicode[STIX]{x1D70E}=k\unicode[STIX]{x1D706}_{D}$ for different values of the screening parameter 
$\unicode[STIX]{x1D705}$.
Equations (6.1) and (6.3) comprise a complete set of amplitude–phase mismatch ODEs describing the passage through the linear resonance in our problem. These equations involve several parameters, but the number of parameters can be reduced to just one by rescaling the problem. Indeed, if we use the slow time 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6FC}^{1/2}t$ and define new amplitude 
$Q$ via 
$Q^{2}={\displaystyle \frac{N}{2\unicode[STIX]{x1D714}_{0}\unicode[STIX]{x1D6FC}^{1/2}}}C_{1}^{2}$, our system reduces to
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\text{d}Q}{\text{d}\unicode[STIX]{x1D70F}}}=-\unicode[STIX]{x1D707}\sin \unicode[STIX]{x1D6F7}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\text{d}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D70F}}}=\unicode[STIX]{x1D70F}-Q^{2}-{\displaystyle \frac{\unicode[STIX]{x1D707}\cos \unicode[STIX]{x1D6F7}}{Q}}, & \displaystyle\end{eqnarray}$$where
$$\begin{eqnarray}\unicode[STIX]{x1D707}={\displaystyle \frac{\unicode[STIX]{x1D700}N^{1.2}}{\sqrt{2}\unicode[STIX]{x1D6FC}^{3/4}\unicode[STIX]{x1D714}_{0}^{3/2}(1+\unicode[STIX]{x1D705}^{2})}}.\end{eqnarray}$$ Note that, if one defines a complex variable 
$\unicode[STIX]{x1D6F9}=Q\text{e}^{\text{i}\unicode[STIX]{x1D6F7}}$, our system is further reduced to a single complex ODE
$$\begin{eqnarray}\text{i}\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D70F}}+(\unicode[STIX]{x1D70F}-|\unicode[STIX]{x1D6F9}|^{2})\unicode[STIX]{x1D6F9}=\unicode[STIX]{x1D707}\end{eqnarray}$$ characteristic of AR problems in many different physical systems and studied in numerous applications (Friedland Reference Friedland2009). For example, if initially 
$\unicode[STIX]{x1D6F9}=0$ and one starts at sufficiently large negative 
$\unicode[STIX]{x1D70F}$ (i.e. far from the linear resonance), this equation predicts transition to AR at large positive 
$\unicode[STIX]{x1D70F}$ if 
$\unicode[STIX]{x1D707}$ is above the threshold 
$\unicode[STIX]{x1D707}_{th}\approx 0.41$, or returning to our original parameters,
$$\begin{eqnarray}\unicode[STIX]{x1D700}>\unicode[STIX]{x1D700}_{th}=0.58{\displaystyle \frac{\unicode[STIX]{x1D6FC}^{3/4}\unicode[STIX]{x1D714}_{0}^{3/2}(1+\unicode[STIX]{x1D705}^{2})}{N^{1.2}}}.\end{eqnarray}$$ The threshold (6.8) assumes its simplest form when either 
$\unicode[STIX]{x1D6E5}$ or 
$\unicode[STIX]{x1D705}$ vanish. Indeed, in the cold plasma case, 
$\unicode[STIX]{x1D6E5}=0$ and
$$\begin{eqnarray}N={\displaystyle \frac{3\unicode[STIX]{x1D705}^{2}(3+\unicode[STIX]{x1D705}^{2})}{32}}.\end{eqnarray}$$ Therefore, for small 
$\unicode[STIX]{x1D705}$, 
$\unicode[STIX]{x1D700}_{th}$ scales as 
$\unicode[STIX]{x1D700}_{th}\sim \unicode[STIX]{x1D6FC}^{3/4}/\unicode[STIX]{x1D705}$ and one needs a non-vanishing screening factor 
$\unicode[STIX]{x1D705}$ to get a sufficiently small 
$\unicode[STIX]{x1D700}_{th}$ for the AR excitation. In the case of 
$\unicode[STIX]{x1D705}=0$ and 
$\unicode[STIX]{x1D6E5}\ll 1$,
$$\begin{eqnarray}N\approx {\displaystyle \frac{15\unicode[STIX]{x1D6E5}^{2}}{8}},\end{eqnarray}$$ so 
$\unicode[STIX]{x1D700}_{th}$ scales as 
$\unicode[STIX]{x1D700}_{th}\sim \unicode[STIX]{x1D6FC}^{3/4}/\unicode[STIX]{x1D6E5}$ (recall that 
$\unicode[STIX]{x1D6E5}=3^{1/2}k\unicode[STIX]{x1D706}_{D}$ and measures the electron thermal spread in the problem). We illustrate these results in figure 6 showing 
$\unicode[STIX]{x1D700}_{th}/\unicode[STIX]{x1D6FC}^{3/4}$ versus 
$\unicode[STIX]{x1D70E}$ (this is 
$k\unicode[STIX]{x1D706}_{D}$ in dimensional notations) for several values of 
$\unicode[STIX]{x1D705}$.
Figure 7. The control of the autoresonant SPW by tapering the driving frequency. (a) The envelope of the maxima of the electron fluid velocity 
$v$ versus slow time 
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D6FC}^{1/2}t$. The dashed line represents the phase velocity 
$\unicode[STIX]{x1D714}_{d}/k$ of the driving wave and the dotted red line represents the envelope of the maxima of the upper limiting velocity of the water bag model. Panel (b) shows the frequencies of the driving (red line) and driven (blue line) waves. The waves are frequency locked starting 
$\unicode[STIX]{x1D70F}\approx -3$.
Finally, we address the possibility of autoresonant control of large amplitude SPWs. This control uses a one-to-one correspondence between the amplitude and frequency of nonlinear waves. Since in autoresonance the driven wave frequency is locked to the driving frequency, the wave amplitude can be controlled by simply varying the driving frequency. For example, we have seen previously (see figure 3), that a continuous linear time variation of the driving frequency leads to the excited wave approaching the wave breaking limit. But what if one wants to avoid this limit and excite a wave having some given target amplitude? This goal can be achieved by tailoring the time variation of the driving frequency appropriately. We illustrate such control in figure 7 showing the results of the simulations using the water bag model for the case of 
$\unicode[STIX]{x1D70E}=0.05$, 
$\unicode[STIX]{x1D705}=0.75$, 
$\unicode[STIX]{x1D700}=0.0011$ (
$\unicode[STIX]{x1D700}_{th}=0.00084$ in this case) and the driving frequency 
$\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{0}+\unicode[STIX]{x1D6FC}t$ for 
$t<0$, but 
$\unicode[STIX]{x1D714}_{d}=\unicode[STIX]{x1D714}_{0}+(2\text{d}/\unicode[STIX]{x03C0})\text{arctg}(t/T)$ for 
$t>0$, where 
$T=2d/\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}$. This driving frequency approaches a fixed value of 
$\unicode[STIX]{x1D714}_{0}+d$ at large times (we used 
$d=0.033$ in the simulations). Figure 7(a) shows the envelopes of the maxima of the electron fluid velocity 
$v$ (solid line) and of the upper limiting velocity 
$u_{1}$ of the water bag model (dotted line). In the same figure, we also show the driving wave phase velocity (dashed line). Separate simulations show that if the driving frequency in this example is chirped linearly in time, the upper limiting velocity of the water bag approaches the Landau resonance (
$u_{1}\approx \unicode[STIX]{x1D714}_{d}/k$) at 
$\unicode[STIX]{x1D70F}=6$, where wave breaking is expected. In contrast, for the saturated driving frequency case shown in figure 7, the excited fluid velocity also saturates at larger times and performs slow oscillations around some fixed value. The frequency of these oscillations scales as 
$\unicode[STIX]{x1D700}^{1/2}$ (Friedland Reference Friedland2009) and their presence illustrates stability of the AR evolution. The velocity 
$u_{1}$ also saturates and, thus, the kinetic wave breaking is avoided. This saturation is caused by the autoresonant frequency locking in the system as illustrated in figure 7(b), showing the evolution of the frequencies of the driven and driving waves.
7 Conclusions
We have studied excitation and control of large amplitude SPWs by a chirped frequency driving wave. The process involved passage through the linear resonance in the problem and transition to autoresonant stage of excitation, where the driven SPW self-adjusts its amplitude to stay in a continuous resonance with the drive. The method allowed reaching of extreme regimes, where the electron density developed a sharply peaked spatial profile with the maximum electron density exceeding the initial plasma density significantly (see figure 1). These results were illustrated in both Vlasov–Poisson and water bag simulations. The simpler water bag model (3.1) was used to develop the adiabatic theory using the Whitham averaged variational principle for studying the weakly nonlinear stage of slowly chirped, driven autoresonant SPWs. In this regime, the problem was reduced to the standard set of coupled amplitude–phase mismatch equations (6.4) and (6.5) characteristic of many other autoresonantly driven problems. The reduction also allowed finding of the threshold driving amplitude (see (6.8)) for the transition to autoresonance. By slowly decreasing the chirp rate of the driving frequency and reaching some fixed frequency level, one could arrive at a given target amplitude of the autoresonant SPW and avoid the kinetic wave breaking (see figure 7).
The self-consistent inclusion of the variation of the driving amplitude and of three-dimensional effects in the process of autoresonant excitation of SPWs seem to be important goals for future research. The form of the autoresonant SPW suggests that it comprises a two-phase solution, with each phase locked to one of the travelling waves comprising the drive. A better understanding of such waveforms, analysing other autoresonant multi-phase plasma waves and studying details of the kinetic wave breaking process in application to autoresonant SPWs also comprise interesting goals for the future.
Acknowledgements
This work was supported by the US-Israel Binational Science Foundation grant no. 6079 and the Russian state program AAAA-A18- 118020190095-4. The authors are also grateful to J. S. Wurtele, P. Michel and G. Marcus for helpful comments and suggestions.
Appendix A. Second-order amplitudes
To calculate the second-order amplitudes, we substitute equations (4.7) and (4.8) into the Lagrangian (3.4), write the result to second nonlinear order in amplitudes and average over one spatial period. This spatially averaged Lagrangian governs the time dependence of all the amplitudes. Its variations with respect to the second-order amplitudes yield a system of ODEs for these amplitudes versus time. Algebraically, this reduction is a tedious process done here using Mathematica (Wolfram Research Inc. 2017). The resulting set of equations is
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}(e_{1,2})_{t}={\textstyle \frac{1}{2}}(a_{1,2}b_{1,2}\pm 4\unicode[STIX]{x0394}d_{1,2}),\\[5.0pt] (d_{1,2})_{t}={\textstyle \frac{1}{4}}(4c_{2}-a_{1,2}^{2}+b_{1,2}^{2}\mp 8\unicode[STIX]{x0394}e_{1,2}),\\[5.0pt] c_{2}={\displaystyle \frac{e_{2}-e_{1}}{\unicode[STIX]{x1D6E5}(4+\unicode[STIX]{x1D705}^{2})}},\end{array}\right\}\end{eqnarray}$$ or after substituting the first-order dependencies 
$a_{1,2}\sim \cos \unicode[STIX]{x1D703}$, 
$b_{1,2}\sim \sin \unicode[STIX]{x1D703}$, and 
$c_{1}\sim \cos \unicode[STIX]{x1D703}$
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}(e_{1,2})_{t}=p_{1,2}\sin (2\unicode[STIX]{x1D703})\pm 2\unicode[STIX]{x0394}d_{1,2},\\[5.0pt] (d_{1,2})_{t}=q_{1,2}+r_{1,2}\cos (2\unicode[STIX]{x1D703})+c_{2}\mp 2\unicode[STIX]{x0394}e_{1,2},\\[5.0pt] c_{2}={\displaystyle \frac{e_{2}-e_{1}}{\unicode[STIX]{x1D6E5}(4+\unicode[STIX]{x1D705}^{2})}},\end{array}\right\}\end{eqnarray}$$ where 
$p_{1,2}$, 
$q_{1,2}$, 
$r_{1,2}$ are constants. These equations have the following time periodic solutions
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}d_{1,2}=D_{1,2}\sin (2\unicode[STIX]{x1D703}),\\[5.0pt] e_{1,2}=F_{1,2}+E_{1,2}\cos (2\unicode[STIX]{x1D703}),\\[5.0pt] c_{2}=B+C\cos (2\unicode[STIX]{x1D703}),\end{array}\right\}\end{eqnarray}$$ where 
$D_{1,2}$, 
$F_{1,2}$, 
$E_{1,2}$, 
$B$ and 
$C$ are constants.
Appendix B. The reduction of the variational system
All algebraic manipulations in this appendix use Mathematica (Wolfram Research Inc. 2017). We employ the averaged Lagrangian 
$\unicode[STIX]{x1D6EC}$ (see (5.1)–(5.3)) and proceed from the variations with respect to the first-order amplitudes 
$A_{1,2}$, 
$B_{1,2}$, 
$C_{1}$
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}-2C_{1}+2A_{1}F_{1}+B_{1}D_{1}+A_{1}E_{1}+2B_{1}\unicode[STIX]{x1D714}+2A_{1}\unicode[STIX]{x1D6E5}=4\unicode[STIX]{x1D700}\cos \unicode[STIX]{x1D6F7},\\[5.0pt] -2C_{1}+2A_{2}F_{2}+B_{2}D_{2}+A_{2}E_{2}+2B_{2}\unicode[STIX]{x1D714}-2A_{2}\unicode[STIX]{x1D6E5}=4\unicode[STIX]{x1D700}\cos \unicode[STIX]{x1D6F7},\\[5.0pt] 2F_{1}B_{1}-A_{1}(D_{1}+2\unicode[STIX]{x1D714})-B_{1}(E_{1}+2\unicode[STIX]{x1D6E5})=0,\\[5.0pt] -2F_{2}B_{2}+A_{2}(D_{2}+2\unicode[STIX]{x1D714})+B_{2}(E_{2}-2\unicode[STIX]{x1D6E5})=0,\\[5.0pt] A_{1}-A_{2}+2C_{1}(1+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}=0.\end{array}\right\}\end{eqnarray}$$ The first four equations in this set yield the linear approximation in terms of 
$C_{1}$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{1}^{\prime }=-A_{2}^{\prime }=-{\displaystyle \frac{C_{1}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x1D714}^{2}-\unicode[STIX]{x1D6E5}^{2}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle B_{1}^{\prime }=B_{2}^{\prime }={\displaystyle \frac{C_{1}\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D714}^{2}-\unicode[STIX]{x1D6E5}^{2}}}. & \displaystyle\end{eqnarray}$$Then the fifth equation in (B 1) gives the linear dispersion relation
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}^{2}={\displaystyle \frac{1}{1+\unicode[STIX]{x1D705}^{2}}}+\unicode[STIX]{x1D6E5}^{2}.\end{eqnarray}$$ Next, we take variations with respect to the second-order amplitudes 
$D_{1,2}$, 
$E_{1,2}$, 
$F_{1,2}$, 
$B$, 
$C$ (to lowest significant order using the linear result for the first-order amplitudes and 
$\unicode[STIX]{x1D714}$) to get
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}A_{1}^{\prime }B_{1}^{\prime }+8E_{1}\unicode[STIX]{x1D714}_{0}+8D_{1}\unicode[STIX]{x1D6E5}=0,\\[5.0pt] A_{2}^{\prime }B_{2}^{\prime }+8E_{2}\unicode[STIX]{x1D714}_{0}-8D_{2}\unicode[STIX]{x1D6E5}=0,\\[5.0pt] -8C0+A_{1}^{\prime 2}+B_{1}^{\prime 2}+16D_{1}\unicode[STIX]{x1D714}_{0}+16E_{1}\unicode[STIX]{x1D6E5}=0,\\[5.0pt] -8C0+A_{2}^{\prime 2}+B_{2}^{\prime 2}+16D_{2}\unicode[STIX]{x1D714}_{0}-16E_{2}\unicode[STIX]{x1D6E5}=0,\\[5.0pt] 8B-A_{1}^{\prime 2}+B_{1}^{\prime }-16F_{1}\unicode[STIX]{x1D6E5}=0,\\[5.0pt] 8B-A_{2}^{\prime 2}+B_{2}^{\prime 2}+16F_{2}\unicode[STIX]{x1D6E5}=0,\\[5.0pt] F_{1}-F_{2}+B(4+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}=0,\\[5.0pt] E_{1}-E_{2}+C(4+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}=0.\end{array}\right\}\end{eqnarray}$$ This system is now solved for the second-order amplitudes via 
$C_{1}^{2}$
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}D_{1}={\displaystyle \frac{C_{1}^{2}\unicode[STIX]{x1D714}_{0}[2+\unicode[STIX]{x1D714}_{0}^{2}(4+\unicode[STIX]{x1D705}^{2})+3\unicode[STIX]{x1D6E5}^{2}(4+\unicode[STIX]{x1D705}^{2})]}{16(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})^{2}[1-\unicode[STIX]{x1D714}_{0}^{2}(4+\unicode[STIX]{x1D705}^{2})+\unicode[STIX]{x1D6E5}^{2}(4+\unicode[STIX]{x1D705}^{2})]}},\\[10.0pt] E_{1}=-{\displaystyle \frac{C_{1}^{2}\unicode[STIX]{x1D6E5}[3\unicode[STIX]{x1D714}_{0}^{2}(4+\unicode[STIX]{x1D705}^{2})+\unicode[STIX]{x1D6E5}^{2}(4+\unicode[STIX]{x1D705}^{2})]}{16(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})^{2}[1-\unicode[STIX]{x1D714}_{0}^{2}(4+\unicode[STIX]{x1D705}^{2})+\unicode[STIX]{x1D6E5}^{2}(4+\unicode[STIX]{x1D705}^{2})]}},\\[10.0pt] F_{1}={\displaystyle \frac{C_{1}^{2}(4+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}}{16(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})[1+\unicode[STIX]{x1D6E5}^{2}(4+\unicode[STIX]{x1D705}^{2})]}},\\[10.0pt] D_{2}=D_{1},E_{2}=-E_{1},F_{2}=-F_{1},\\[5.0pt] B=-{\displaystyle \frac{C_{1}^{2}}{8(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})[1+\unicode[STIX]{x1D6E5}^{2}(4+\unicode[STIX]{x1D705}^{2})]}},\\[10.0pt] C={\displaystyle \frac{C_{1}^{2}(3\unicode[STIX]{x1D714}_{0}^{2}+\unicode[STIX]{x1D6E5}^{2})}{8(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})^{2}[1-\unicode[STIX]{x1D714}_{0}^{2}(4+\unicode[STIX]{x1D705}^{2})+-\unicode[STIX]{x1D6E5}^{2}(4+\unicode[STIX]{x1D705}^{2})]}}.\end{array}\right\}\end{eqnarray}$$ Finally, we return to the first two equations in (B 1) and solve these equations for 
$A_{1,2}$ to higher (third) order in 
$C_{1}$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{1}^{\prime \prime }=-{\displaystyle \frac{C_{1}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x1D714}^{2}-\unicode[STIX]{x1D6E5}^{2}}}+{\displaystyle \frac{C_{1}G_{1}-4\unicode[STIX]{x1D700}\unicode[STIX]{x1D6E5}(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})\cos \unicode[STIX]{x1D6F7}}{2(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})^{2}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{2}^{\prime \prime }={\displaystyle \frac{C_{1}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x1D714}^{2}-\unicode[STIX]{x1D6E5}^{2}}}+{\displaystyle \frac{C_{1}G_{2}+4\unicode[STIX]{x1D700}\unicode[STIX]{x1D6E5}(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})\cos \unicode[STIX]{x1D6F7}}{2(\unicode[STIX]{x1D714}_{0}^{2}-\unicode[STIX]{x1D6E5}^{2})^{2}}}, & \displaystyle\end{eqnarray}$$ where 
$G_{1}=\unicode[STIX]{x1D714}_{0}^{2}(2F_{1}-E_{1})+2D_{1}\unicode[STIX]{x1D714}_{0}\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6E5}^{2}(2F_{1}+E_{1})$ and 
$G_{2}=\unicode[STIX]{x1D714}_{0}^{2}(2F_{2}-E_{2})-2D_{2}\unicode[STIX]{x1D714}_{0}\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6E5}^{2}(2F_{2}+E_{2})$. Finally, we use (B 6) in (B 7) and (B 8) and substitute the resulting 
$A_{1,2}^{\prime \prime }$ into the last equation in (B 1). This yields a higher-order approximation for the frequency 
$\unicode[STIX]{x1D714}$ of the wave
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=\unicode[STIX]{x1D714}_{0}^{2}+NC_{1}^{2}+{\displaystyle \frac{2\unicode[STIX]{x1D700}\cos \unicode[STIX]{x1D6F7}}{(1+\unicode[STIX]{x1D705}^{2})C_{1}}},\end{eqnarray}$$where
$$\begin{eqnarray}N={\displaystyle \frac{9\unicode[STIX]{x1D705}^{2}(3+\unicode[STIX]{x1D705}^{2})+R}{96[1+(4+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}^{2}]}},\end{eqnarray}$$ and 
$R=(1+\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}^{2}[3(60+56\unicode[STIX]{x1D705}^{2}+11\unicode[STIX]{x1D705}^{4})+20(1+\unicode[STIX]{x1D705}^{2})(4+\unicode[STIX]{x1D705}^{2})(5+2\unicode[STIX]{x1D705}^{2})\unicode[STIX]{x1D6E5}^{2}+16(1+\unicode[STIX]{x1D705}^{2})^{2}(4+\unicode[STIX]{x1D705}^{2})^{2}\unicode[STIX]{x1D6E5}^{4}]$.
