1. Introduction
The dynamics of small-amplitude electron plasma waves (EPWs) is described by the complex root,
$\omega _L=\omega _{\text{EPW}} + i \nu _{\text{L}}$
, of the linearised dispersion function
$\varepsilon (\omega , k_1; f(\omega )) = 1 - \omega _p^2\int f / (\omega -k_1 v)^2 \, {\rm d}v$
. Assuming a Maxwellian distribution function, EPWs oscillate at a frequency,
$\omega _{\text{EPW}}$
, and damp at the Landau damping rate,
$\nu _{\text{L}}$
(Landau Reference Landau1946; Canosa Reference Canosa1973). In the case of large-amplitude EPWs, Bernstein, Greene & Kruskal (Reference Bernstein, Greene and Kruskal1957) derived a stationary distribution function where the EPW is sustained due to particle trapping in the potential well of the wave. Such a state, often termed a BGK mode, is described by a nearly Maxwell–Boltzmann distribution function and a nonlinear correction due to quasi-linear diffusion (Davidson Reference Davidson2012) at the phase velocity of the EPW. Manheimer & Flynn (Reference Manheimer and Flynn1971), Dewar (Reference Dewar1972) Morales & O’Neil (Reference Morales and O’Neil1972) showed that a BGK mode results in a small negative shift in the wave frequency for a self-consistent solution. In addition to the negative shift in the real part,
$\Delta \omega _{NL}$
, O’Neil (Reference O’Neil1965) shows that the imaginary part, i.e. the damping rate, becomes negligible and the dispersion relation is given by
$\omega _{NL} = \omega _{\text{EPW}} + \Delta \omega _{\text{NL}} + 0i$
. These results rely on constant-amplitude, sudden-turn-on or asymptotic approximations. A more general, non-local, neo-adiabatic framework has since been developed by Bénisti et al. (Reference Bénisti, Strozzi, Gremillet and Morice2009), Bénisti (Reference Bénisti2017a
,
Reference Bénistib
): the damping rate (Bénisti et al. Reference Bénisti, Strozzi, Gremillet and Morice2009), the action distribution function (Bénisti Reference Bénisti2017a
) and the nonlinear frequency shift (Bénisti Reference Bénisti2017b
) of a slowly varying, externally driven EPW are all shown to be history-dependent functionals of the wave amplitude rather than local functions of its instantaneous value. In particular, Bénisti (Reference Bénisti2017b
) shows that this non-locality produces a hysteresis in
$\delta \omega$
over a grow-then-decay amplitude loop, even in a collisionless plasma. A parallel, collisionless adiabatic framework was developed by Dodin & Fisch (Reference Dodin and Fisch2012a
,
Reference Dodin and Fischb
,
Reference Dodin and Fischc
), Dodin (Reference Dodin2014), in which a Whitham-type averaged Lagrangian in oscillation-centre variables yields the nonlinear dispersion relation and frequency shift
$\omega _{\text{NL}}(a)$
directly from the trapped-particle distribution. For smooth trapped distributions, this recovers
$\omega _{\text{NL}} \propto \sqrt {a}$
, but beam-like distributions produce different power laws or a logarithmic nonlinearity and the trapped-particle autoresonance makes
$\omega _{\text{NL}}$
generically non-local. Both the Bénisti and Dodin–Fisch theories are built on the collisionless Vlasov–Poisson system. They describe the Phase I trapping and plateau formation considered later, and, because the slowly-decaying Phase II amplitude keeps the dynamics in their neo-adiabatic regime (the separatrix contraction rate
$\dot v_{\text{tr}}/v_{\text{tr}} \propto \nu _{\text{eff}}$
is small compared with
$\omega _B$
), they continue to contribute a collisionless, adiabatic-recession component to the Phase II trapped-particle response. As shown in the remainder of this paper, however, Phase II also requires a Fokker–Planck contribution to the susceptibility: the net nonlinear frequency shift emerges from a self-consistent balance between the collisionless (Vlasov) and collisional (Fokker–Planck) contributions, and cannot be captured by the collisionless frameworks alone.
Previous works assume idealised, collisionless plasmas. Plasmas accessible in the laboratory, e.g. those in fusion and gas discharge experiments, and in Nature, e.g. in the ionosphere and the interstellar medium, are often weakly collisional such that the ratio of
$\nu _{\textit{ee}} / \omega _{p} \sim 10^{-4} {-} 10^{-10}$
, where
$\nu _{\textit{ee}}$
is the electron–electron (e–e) collision frequency (Thorne & Blandford Reference Thorne and Blandford2017). While it is implicitly known that in collisional systems, BGK modes cannot persist forever and that over time, the wave must exhibit some damping as the particles thermalise, the dynamics of this transition from nonlinear to the linear state is not well understood.
Using asymptotic analysis, Zakharov & Karpman (Reference Zakharov and Karpman1963) (ZK) determined that during the detrapping phase, the effective damping rate in weakly collisional, nonlinear EPW can be given by
$ \nu _{\text{ZK}} \propto ({\nu _{\textit{ee}} \omega _{\text{EPW}}^2}/{\omega _B^3}) \nu _{\text{L}}$
, where
$\nu _{\textit{ee}}$
is the e–e collision frequency,
$\omega _{\text{EPW}}$
is the linear resonance,
$\omega _B$
is the bounce frequency of the trapped particles,
$\nu _L$
is the linear Landau damping rate. While the calculation of this relation contains approximations e.g. fluid Landau damping rate, it correctly suggests that waves can have effective damping rates larger than the base e–e collision rate because the rate is inversely proportional to
$\omega _B/\omega _p$
. Numerical verification of this has been limited to a few simulations using a Particle-In-Cell code where the simulations found qualitative agreement (Valentini Reference Valentini2008).
This paper provides an end-to-end description of the wave dynamics in weakly collisional, large-amplitude EPWs which exhibits three distinct phases. Phase I comprises well-understood collisionless mechanisms of wave excitation, particle trapping and the establishment of a quasi-linear plateau. The physics in this phase can be described by a second-order perturbative expansion of the Vlasov–Poisson system with the electric field,
$E$
, as the small parameter. Phase II is the longest-lived phase during which a quasi-steady state is reached between weak electron–electron collisions and strong wave–electron interactions. Describing Phase II requires a first-order expansion in
$\nu _{\textit{ee}}$
as another small parameter in addition to the second-order expansion in
$E$
. The wave evolution concludes with Phase III when the distribution function has returned to a near-Maxwellian and the wave energy is rapidly damped via Landau damping. Only the first-order expansion in
$E$
is required to describe this phase accurately. While the behaviour in Phases I and III has previously been elucidated, this work describes collisional effects on the wave dispersion relation during Phase II and characterises the lifetime of Phase II. With these quantities, weakly collisional, large-amplitude EPW can be modelled as a piecewise ordinary differential equation.
The remainder of this paper is organised as follows. Section 2 describes the simulation configuration. Section 3 provides an overview of the overall wave dynamics, and introduces the three distinct phases of evolution and highlights the enhancement of the nonlinear frequency shift during the detrapping phase. Section 4 describes the physical mechanism for this enhancement of the nonlinear frequency shift. Sections 5, 6 and 7 describe each phase in detail along with empirical fits for the quantities described in table 1. Finally, § 8 summarises the findings and discusses future work.
Three distinct phases of nonlinear electron plasma wave evolution are framed in terms of perturbative expansions necessary to describe each phase as well as the resulting dispersion relation of each phase. The second and third columns give respectively the lowest order
$m$
in wave amplitude
$E$
and order
$n$
in collision frequency
$\nu _{\textit{ee}}$
at which the dispersion-relation contributions in that phase appear.

Table 1. Long description
The table presents three distinct phases of nonlinear electron plasma wave evolution, detailing the order of wave amplitude and collision frequency for each phase. It has three rows and five columns. The columns are labeled Phase, m in O(E^m), n in O(vee^ee), Re(omega), and Im(omega). Row 1: Phase I, m in O(E^m) is 2, n in O(vee^ee) is 0, Re(omega) is omega subscript EPW plus delta omega subscript NL, Im(omega) is 0. Row 2: Phase II, m in O(E^m) is 2, n in O(vee^ee) is 1, Re(omega) is omega subscript EPW plus delta omega subscript NL times exp(gamma delta omega t), Im(omega) is f subscript NL times v subscript L. Row 3: Phase III, m in O(E^m) is 1, n in O(vee^ee) is 0, Re(omega) is omega subscript EPW, Im(omega) is f subscript L times v subscript L.
2. Simulation configuration
Weakly collisional, large amplitude electron plasma waves are studied using a 1-D1V Eulerian Vlasov–Poisson–Fokker–Planck solver (Joglekar & Levy Reference Joglekar and Levy2020) with periodic boundary conditions. The equations that are solved are given by
While Particle-In-Cell methods are a common tool for studying plasma kinetic behaviour, the slow convergence with particle number is a problem for resolving subtle, kinetic phenomena. In these simulations,
$N_x = 256, N_v = 2048,$
$L_x = 2\pi / k\lambda _D, \Delta t = 0.5, \omega _{pe} t_{\text{max}} = 4000$
. Exponential integrators with spectral derivatives in position and velocity space for the different components of the Vlasov equation are used (Thomas Reference Thomas2016). The simulations step through the phase-space integrators with a 6th-order integrator in time. Poisson’s equation is solved using a spectral solver. The advection–diffusion Fokker–Planck equation (Dougherty Reference Dougherty1964),
$\mathcal{C}(f) = ({\partial }/{\partial v}) (\bar {v} f + \langle v - \bar {v} \rangle ^2 ( {\partial f}/{\partial v}))$
, is solved using a first-order Euler implicit solver with centre differencing of the velocity derivatives. The distribution function is initialised uniformly in space, that is, the density and temperature profiles are uniform. The distribution function is initially Maxwellian.
EPWs are driven using a ponderomotive driver as shown by Joglekar & Thomas (Reference Joglekar and Thomas2023). Because a Hilbert transform of the electric field is necessary for determining the time-dependent oscillation frequency, the driver is started
$450 \omega _p^{-1}$
into the simulation. This prevents the noise that the Hilbert transform introduces into the beginning and end of the transformed signal from obfuscating the changing frequency during the drive phase and Phase I. The parameters for the time envelope are
$t_c = 500 \omega _p^{-1}, t_w = 100 \omega _p^{-1}, t_R = 25 \omega _p^{-1}, \omega =\omega _{\text{EPW}}(k\lambda _D)$
. This results in a driver profile given in figure 1. The amplitude of this driver profile is multiplied by
$\sin (k x - \omega _{\rm EPW} t)$
, where
$\omega _{\rm EPW}$
is given by the solution to the linearised dispersion relation for small amplitude EPWs.
Envelope amplitude of the external forcing term over time. This value is then multiplied by
$\sin (k x - \omega _{\rm EPW} t)$
to give the external forcing term on the space–time grid.

Here, 1600 simulations are performed by scanning over
in 10, 10 and 16 steps, respectively. The parameter ranges were chosen to ensure all three phases are observable within computationally tractable simulation times. For
$k \lambda _D \lt 0.3$
, weak Landau damping requires either very long simulations or smaller amplitudes to observe Phase III. For
$k \lambda _D \gt 0.35$
, strong damping necessitates larger amplitudes, which would push portions of the parameter space into either purely linear or purely nonlinear regimes.
The following section discusses power-law fits to data collected from the 1600 simulations and attempts to explain the dynamics that surface in those data.
3. Overall wave dynamics
Figure 2 shows the three distinct phases of weakly collisional, large-amplitude EPW evolution. Figure 2(a) shows the evolution of the wave amplitude. Phase I is very short-lived, as it encompasses the relatively fast dynamics of wave excitation and particle trapping. The dynamics is primarily collisionless, where wave damping is diminished due to quasi-linear diffusion (Bernstein et al. Reference Bernstein, Greene and Kruskal1957; O’Neil Reference O’Neil1965; Davidson Reference Davidson2012) and there is a corresponding shift in the frequency of oscillation of the wave (Manheimer & Flynn Reference Manheimer and Flynn1971; Dewar Reference Dewar1972; Morales & O’Neil Reference Morales and O’Neil1972; Berger et al. Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013). Figure 2(b) shows the spatially averaged distribution function
$|\hat {f}^{0}|(v)$
in the trapping region. During Phase I, it exhibits flattening at the phase velocity of the wave, as expected via quasi-linear diffusion. It also exhibits heating of the particles in a region
${\approx} \pm\!1.5 v_{\text{tr}}$
from the phase velocity, where the trapping velocity
$v_{\text{tr}} = \sqrt {2E/k}$
and particularly where
$v \gt v_{\text{ph}} + v_{\text{tr}}$
. The two velocity scales play distinct roles:
$v_{\text{ph}} = \omega /k$
is the wave-frame location of the resonant region and is amplitude-independent, while
$v_{\text{tr}}$
is the amplitude-dependent half-width of the separatrix about
$v_{\text{ph}}$
. This modification of the distribution function in the trapping region results in the well-established nonlinear frequency shift where the EPW frequency downshifts by a few percent as shown in figure 2(c).
(a) Wave amplitude, (b) spatially averaged distribution function and (c) wave oscillation frequency are plotted during three distinct phases of wave evolution. The vertical dotted line in panel (b) indicates the EPW phase velocity
$v_{\text{ph}}$
.

During Phase II, figure 2(a) shows that the wave amplitude is weakly damped. Data from simulations indicate that the rate is lower than the Landau damping rate, but higher than the collision frequency. In figure 2(b), it is seen that the distribution function has developed a plateau at the phase velocity. The plateau spans the trapping region, with width
${\approx} 2 v_{\text{tr}}$
, which slowly shrinks during Phase II as the wave damps and
$v_{\text{tr}} = \sqrt {2E/k}$
decreases. The slope
$\partial f/\partial v$
at
$v_{\text{ph}}$
remains small relative to the Maxwellian value but is finite, consistent with the weak, nearly constant damping rate observed in figure 2(a). Compared with the distribution function in Phase I, there are more particles in the trapping region and fewer particles in
$v \gt v_{\text{ph}} + v_{\text{tr}}$
. During this time, the wave frequency clearly decreases, exhibiting a further shift away from the value dictated by linear behaviour,
$\omega _{\rm EPW} = 1.16$
.
After
${\gt}\;\mathcal{O}(100)$
bounce periods, figure 2(a), the damping rate increases rapidly to the order of the Landau damping rate. This indicates the beginning of Phase III. Figure 2(b) shows that the trapped particles have thermalised and the distribution function has become nearly Maxwellian. Consequently, figure 2(c) shows that the EPW frequency reverts to that of a small-amplitude, slightly warmer and strongly damped EPW.
The effects on the dispersion relation during the three phases are summarised in table 1, where
$\omega _{\text{EPW}}$
and
$\nu _L$
are the real and imaginary part of the solution to the linearised dispersion relation, respectively. It is found that for
$0.3 \leqslant k \lambda _D \leqslant 0.35, 0.1 \lt \omega _B \lt 0.25$
and
$10^{-5} \lt \nu _{\textit{ee}} \lt 3 \times 10^{-4}$
, the dynamics can be described by using the set of quantities described in the equations in the table.
Table 2 provides empirical fits for the quantities described in table 1 as power-laws of the form
$Q = C (k\lambda _D)^\alpha (\omega _B)^\beta (\nu _{\textit{ee}})^\gamma$
, where
$C, \alpha , \beta , \gamma$
are fit parameters. Each of the fits has an
$R^2 \gt 0.95$
in real and log space indicating that the empirical power laws work well for the parameter space studied here.
Best fits for the quantities described in table 1 as power-laws of the form
$Q = C (k\lambda _D)^\alpha (\omega _B)^\beta (\nu _{\textit{ee}})^\gamma$
, where
$C, \alpha , \beta , \gamma$
are fit parameters. Uncertainties
${\lt}0.005$
are reported as 0.00.

Table 2. Long description
A table with six rows and five columns comparing empirical fits for various quantities as power-laws. The columns are labeled Quantity, Coefficient, alpha(kλD), beta(ωB), gamma(vece), and R-squared. Row 1: Quantity, ΔωN L, Coefficient, 93.00 ± 1.75, alpha(kλD), 5.92 ± 0.01, beta(ωB), 0.81 ± 0.00, gamma(vece), 0.06 ± 0.00, R-squared, 0.99. Row 2: Quantity, γΔω, Coefficient, 0.03 ± 0.00, alpha(kλD), -0.93 ± 0.06, beta(ωB), -0.88 ± 0.01, gamma(vece), 0.80 ± 0.00, R-squared, 0.97. Row 3: Quantity, fN L, Coefficient, 9.51 ± 0.87, alpha(kλD), -1.43 ± 0.07, beta(ωB), -1.48 ± 0.02, gamma(vece), 1.07 ± 0.00, R-squared, 0.99. Row 4: Quantity, τN L, Coefficient, 0.13 ± 0.01, alpha(kλD), -3.82 ± 0.03, beta(ωB), 1.23 ± 0.01, gamma(vece), -0.86 ± 0.00, R-squared, 0.99. Row 5: Quantity, τf L, Coefficient, 1.73 ± 0.18, alpha(kλD), -2.25 ± 0.06, beta(ωB), 0.78 ± 0.01, gamma(vece), 0.15 ± 0.01, R-squared, 0.92.
The next section focuses on the novel physics in these simulations, the increase in the nonlinear frequency shift due to collisions during Phase II. The subsequent sections describe the dynamics in each phase in detail and provide empirical fits for the quantities described in table 1.
4. Physical mechanism for frequency shift enhancement
To determine the reason for the counter-intuitive increase in the frequency shift, the contribution to the dielectric function from the evolution operators in the Vlasov–Poisson–Fokker–Planck (VPFP) system of equations can be calculated. The evolution equation of quasi-linear theory (Davidson Reference Davidson2012) for the spatially averaged distribution for a one-component collisional, unmagnetised plasma can be written as
\begin{equation} \frac {\partial\!\hat {f}^0}{\partial t} + \frac {q}{m}\left (\hat {E}^{1\star } \frac {\partial\!\hat {f}^{1}}{\partial v} + \hat {E}^{1} \frac {\partial\!\hat {f}^{1\ast }}{\partial v}\right )= \nu _{\textit{ee}} \,\mathcal{C}(\hat {f}^0). \end{equation}
Performing a perturbative expansion assuming
$\nu _{\textit{ee}} \ll \omega _B$
, that is, the collision frequency is much smaller than the bounce frequency, and using
$\hat {f}^0 = \hat {f}^0_0 + \nu _{\textit{ee}} \hat {f}^0_1$
where the subscript indicates the order in
$\nu _{\textit{ee}}$
, (4.1) becomes
\begin{align} \frac {\partial\!\hat {f}^0_0}{\partial t} &+ \frac {q}{m}\left (\hat {E}^{1\star }_0 \frac {\partial\!\hat {f}^{1}_0}{\partial v} + \hat {E}^{1}_0 \frac {\partial\!\hat {f}^{1\star }_0}{\partial v}\right ) = 0 ,\\[-14pt]\nonumber \end{align}
\begin{align} \frac {\partial\!\hat {f}^0_1}{\partial t} &+ \frac {q}{m}\left (\hat {E}^{1\star }_0 \frac {\partial\!\hat {f}^{1}_1}{\partial v} + \hat {E}^{1}_1 \frac {\partial\!\hat {f}^{1\star }_0}{\partial v}\right ) = \hat {C} (\hat {f}^0_0) . \end{align}
Equation (4.2) can be written as the quasilinear diffusion equation that leads to the well-understood collisionless dynamics of the flattening of the distribution function and the nonlinear frequency shift of Phase I. Equation (4.3) is the weakly collisional perturbative correction to the collisionless spatially averaged distribution function
$\hat {f}^0_0$
. It is
$\hat {f}^0_1$
that is responsible for the increase in the frequency shift. That is, the collisional perturbation to the distribution function leads to an enhancement of the nonlinear frequency shift.
This can be formalised by using the dielectric function given by
$ \varepsilon (\omega ;k_1) = 1 - \omega _p^2 \int _{-\infty }^{\infty } {\hat {\tilde {f}}^0}/{(\omega - k_1 v)^2}\,{\rm d}v$
, where
$\hat {\tilde {f}}^0 = \hat {f}^0/n_0$
and
$n_0$
is the plasma density, with
$\omega _p^2 = q^2n_0/m\varepsilon _0$
. If
$\hat {\tilde {f}}^0 = \tilde {f}_M$
, a Maxwellian distribution, the linear wave dispersion relation has the real frequency
$\omega _1$
, which is the solution to
$\varepsilon (\omega _1;k_1) = 0$
. If the dispersion function is expanded about the linear frequency for
$\omega = \omega _1 + \delta \omega$
, to lowest order,
$ \delta \omega = \omega _p^2[\left .({\partial \varepsilon }/{\partial \omega })\right |_{\omega =\omega _1}]^{-1}\int _{-\infty }^{\infty } ({\delta \hat {f}^0})/{(\omega _1 - k_1v)^2}\,{\rm d}v,$
i.e. the nonlinear change in
$\hat {f}^0$
from the coupling of
$f_1$
to the field leads to a well-known frequency shift (Manheimer & Flynn Reference Manheimer and Flynn1971; Dewar Reference Dewar1972; Morales & O’Neil Reference Morales and O’Neil1972).
The time derivative of the frequency shift is
\begin{equation} \frac {\partial }{\partial t}\delta \omega = \omega _p^2\bigg[\left .\frac {\partial \varepsilon }{\partial \omega }\right |_{\omega =\omega _1}\bigg]^{-1}\int _{-\infty }^{\infty } \frac {1}{(\omega _1 - k_1v)^2}\frac {\partial \hat {\tilde {f}}_0}{\partial t}\,{\rm d}v, \end{equation}
with
${\partial \hat {\tilde {f}}_0}/{\partial t}$
given by (4.1) which includes the fast time scale rate of evolution from the normal quasilinear theory (Davidson Reference Davidson2012) and the slow collisional relaxation of the distribution function. Expanding (4.4) in
$\nu _{\textit{ee}}$
gives
\begin{equation} \frac {\partial (\delta \omega _0 + \nu \delta \omega _1)}{\partial t} = \bigg [\!\left .\frac {\partial \varepsilon }{\partial \omega }\right |_{\omega =\omega _1}\!\bigg ]^{-1} \big[\chi \big(\partial _t \hat {f}^0_{0}\big) + \nu \chi \big(\partial _t \hat {f}^0_{1}\big)\big], \end{equation}
where
$\chi (f) = \omega _p^2\int f / (\omega _1 - k_1 v)^2 \,{\rm d}v$
is the susceptibility function.
The first term on each side represents the frequency shift described by Manheimer & Flynn (Reference Manheimer and Flynn1971), Dewar (Reference Dewar1972), Morales & O’Neil (Reference Morales and O’Neil1972). The second term is the collisional contribution and is given by substituting (4.3) for
$\partial _t \hat {f}^0_1$
, which gives
where
$\mathcal{V}(f_{0}^1, f_{1}^1) = E_{0}^{1*} ({\partial f_{1}^1}/{\partial v}) + E_{1}^{1*} ({\partial f_{1}^{1}}/{\partial v}) + \text{c.c.}$
is the linearised contribution from the Vlasov terms in (4.1) and
$E^1 = E_{0}^1 + \nu E_{1}^1$
. By using this relation, the effect on the dielectric function of the Vlasov terms and the Fokker–Planck terms can be isolated and compared. This comparison is performed numerically in fully nonlinear VPFP simulations with a Dougherty collision operator (Dougherty Reference Dougherty1964), given by
$\mathcal{C}(f) = \partial _v (\bar {v} f + \bar {v}_{\text{th}}^2 \partial _v f)$
, where
$\bar {v} = \int {\rm d}v\, v \,f$
and
$\bar {v}_{th}^2=\int {\rm d}v \, (v-\bar {v})^2 \, f$
. This collision operator conserves density, momentum and energy, and performs a linearised approximation of the phase-space drag and diffusion from small-angle collisions.
Evolution of the contributions to the susceptibility function.

Figure 3 compares the change in the dielectric function,
$\Delta \epsilon$
, from the Vlasov terms to
$\Delta \epsilon$
from the collision operator. It shows that during Phase I, the Vlasov contribution dominates as expected from the theory from Manheimer & Flynn (Reference Manheimer and Flynn1971), Dewar (Reference Dewar1972), Morales & O’Neil (Reference Morales and O’Neil1972). However, during Phase II, the quasi-steady equilibrium between the Vlasov and Fokker–Planck term is responsible for a net increase in the frequency shift. The Fokker–Planck term always has a positive contribution while the Vlasov contribution counteracts this behaviour, but only partially, resulting in a net increase in the frequency shift of the EPW. The subsequent increase can be roughly equal in magnitude to the initial collisionless nonlinear frequency shift from Phase I. This behaviour is persistent over all simulations.
Physically, the collision operator continuously scatters particles at the edges of the trapping region. If the particle trapping is not saturated, there will remain a negative slope at the phase velocity of the wave and wave–particle interactions continue to occur. This quasi-steady competition results in a net redistribution of phase space density that contributes positively to the susceptibility function. The precise mechanism by which this balance produces an enhanced frequency shift, rather than simply restoring linear behaviour, remains an open question for future theoretical work.
In the sections that follow, each phase is described in detail along with empirical fits for the quantities described in table 1.
5. Phase I: trapping and plateau formation
Phase I is a short-lived phase that encompasses wave excitation and particle trapping. The modifications to the dispersion relation include the loss of wave damping and a nonlinear frequency shift. In the following sections that describe each phase in detail, four simulations, each with
$k\lambda _D = 0.3, a_0 = 0.1$
and varying collision frequencies are used as examples. These are given in table 3.
Simulations used in the remainder of the manuscript have the same
$k\lambda _D = 0.3$
and drive strength
$a_0$
, and varying collision frequencies
$\nu _{\textit{ee}}$
.

Table 3. Long description
A table comparing simulations with varying collision frequencies. The table has four rows and two columns. The columns are labeled Simulation and v_ec/omega_pe. The row labels are A, B, C, and D. The values in the v_ec/omega_pe column are 10^-5, 10^-4, 1.2 x 10^-4, and 3.6 x 10^-4 respectively.
5.1. Phase space evolution
During Phase I, particles near the phase velocity of the wave are trapped in the potential well. This creates a phase space vortex near the phase velocity of the wave.
Phase space snapshots at
$\omega _{pe} t = 650$
are shown in figure 4(a) for Simulations A and B. Both snapshots show a similarly sized vortex. This is because the effect of weak collisions and Landau damping has not yet surfaced and the wave amplitudes are nearly the same. However, the snapshot from Simulation B exhibits more smoothing throughout and the features are not as sharp, especially near the edges of the trapping region.
Figure 4 also shows that both waves are roughly in phase with one another, as indicated by the fact that the vortex is centred around 10
$\lambda _D$
in both cases.
Phase space portrait at
$t=650\omega _{pe}^{-1}$
, shortly after the drive is turned off, for Simulations A and B.

Evolution of the (a) wave amplitude (b) and instantaneous frequency during Phase I for Simulations A and B. The dashed line in panel (b) indicates the linear frequency.

5.2. Wave evolution and nonlinear frequency shift
The wave evolution early in time in Simulations A and B is plotted in figure 5(a). It shows that the magnitudes of the driven waves in both simulations are within a few percent of one another and only weakly affected by collisions. However, there are two differences that should be highlighted.
The first is that the effect of the weak damping is already evident in the first
$200 \omega _{pe}t$
of the simulation after the drive is turned off around
$600 \omega _{pe} t$
. The difference in the damping rates is clear in figure 5(a), where the wave amplitude in Simulation B decays faster than that in Simulation A. The instantaneous frequency of the wave is plotted in figure 5(b) and shows similar behaviour in that the initial frequency shifts are within a few percent of each other, but the decrease in the frequency as the wave amplitude is damped in Simulation B is evident.
Second, the oscillations in the wave amplitude corresponding to particle trapping are suppressed faster in Simulation B. This is likely due to the additional smoothing of the vortex that is evident in figure 4. Similar to the wave amplitude behaviour, the bounce oscillations in the wave frequency are suppressed more quickly in Simulation B.
The nonlinear frequency shift during Phase I is measured by performing a Hilbert transform of the electric field to extract the instantaneous frequency of oscillation. The frequency shift,
$\Delta \omega = \omega - \omega _{\text{EPW}}$
, is measured at the end of Phase I, which is defined as the time when the wave amplitude reaches its first maximum after the drive has been turned off. Fitting the data from all 1600 simulations results in the empirical relation given by
where
$\omega _B = \sqrt {k E}$
is the bounce frequency of the trapped particles. This fit has an R
$^2$
value of 0.996 indicating excellent agreement with the data.
The nonlinear frequency shift given by (5.1) has a strong dependence on
$k\lambda _D$
because the frequency shift is a direct result of the rearrangement of the phase space density. The phase space density is a strong function of the phase velocity of the wave which is a strong function of
$k\lambda _D$
.
The fact that Phase I comprises collisionless dynamics is also reflected in the effectively negligible dependence on the collision rate in (5.1).
Nonlinear frequency shift versus normaliesd wave amplitude. Data from adept and Berger et al. (Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013) along with the theoretical curve (line) from Dewar (Reference Dewar1972) are compared.

Finally, (5.1) has a weaker dependence on the wave amplitude than that reported by Manheimer & Flynn (Reference Manheimer and Flynn1971), Dewar (Reference Dewar1972), Morales & O’Neil (Reference Morales and O’Neil1972), Berger et al. (Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013), where
$\Delta \omega \propto \omega _B$
. While an analytical derivation of this scaling is outside the scope of this work, comparisons to prior simulations and theory are performed. Figure 6 shows a direct comparison to the plot of Berger et al. (Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013). The data from adept agree well with the data from Berger et al. (Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013) and both datasets are within a few percent of the theoretical curve of Dewar (Reference Dewar1972). This figure suggests that adept is performing similarly to the Vlasov code used by Berger et al. (Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013) and that the discrepancy between the measured scaling and the theoretical scaling is not due to numerical artefacts, but rather due to the approximations used in the derivation of the theory. Further, Berger et al. (Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013) noted that the measured nonlinear frequency shift is sensitive to the specifics of the drive duration.
It is also important to consider that the linear scaling was derived using asymptotic approximations by Manheimer & Flynn (Reference Manheimer and Flynn1971), Dewar (Reference Dewar1972), Morales & O’Neil (Reference Morales and O’Neil1972). Specifically, Morales & O’Neil (Reference Morales and O’Neil1972) stipulate that
$v_{tr} v_{ph} \ll v_{th}^2$
. In these simulations,
$v_{ph} \approx 3.3 v_{th}$
and
$v_tr \gtrapprox 0.5$
. They also stipulate that their theory requires
$v_ph \geqslant 4 v_{th}$
, which leads to
$k \lambda _D \lt 0.3$
. In these simulations,
$0.3 \leqslant k \lambda _D \leqslant 0.35$
. These violations of the assumptions used in the derivation of the theory could explain the discrepancy in the scaling with wave amplitude, while the close agreement with previous simulations as well as the theoretical curve in figure 6 suggest that the simulations are reasonably accurately capturing the nonlinear frequency shift physics.
5.3. Physical mechanism for the nonlinear frequency shift
Figure 7(a) distinguishes the change in the spatially averaged distribution function over the first 600
$\omega _{pe}t$
from the Vlasov terms and the Fokker–Planck terms. It shows that the Vlasov terms cause a substantial increase in density in the trapping region, while the contribution from the Fokker–Planck terms is negligible.
(a) Contribution to
$\Delta |\hat {f}^0(v)|$
integrated over the first
$\omega _{pe} t = 600$
. (b) Corresponding contribution to
$\Delta \chi (v)$
, which is directly proportional to
$\Delta |\hat {f}^0|$
.

Figure 7. Long description
Panel A: A line graph shows the contribution to the integrated value over the first period. The x-axis is labeled v/v_th and ranges from 3 to 5. The y-axis is labeled Δ|f^0| and ranges from -1 to 2 times 10^-6. Two lines are plotted: a blue line representing the Vlasov model and a red line representing the Fokker-Planck model. The blue line shows more significant fluctuations and peaks compared to the red line. Panel B: A line graph shows the corresponding contribution to Δχ_v, which is directly proportional to Δ|f^0|. The x-axis is labeled v/v_th and ranges from 3 to 5. The y-axis is labeled Δχ_v and ranges from 0 to 0.00010. The blue line representing the Vlasov model shows notable peaks and fluctuations, while the red line representing the Fokker-Planck model remains relatively stable with minor fluctuations.
Figure 7(b) shows the corresponding contribution to the susceptibility function,
$\Delta \chi (v) = \omega _{pe}^2 \int _0^v \partial _v' \Delta |\hat {f}^0(v')| / (\omega - k v') \,{\rm d}v'$
, which is directly proportional to the frequency shift. The different contributions are calculated by summing the change in the distribution function from each of the terms in the simulations. It shows that the Vlasov terms are responsible for the nonlinear frequency shift, while the Fokker–Planck terms have a negligible effect. This confirms that the nonlinear frequency shift during Phase I is a collisionless effect arising from wave–particle interactions, as expected from Manheimer & Flynn (Reference Manheimer and Flynn1971), Dewar (Reference Dewar1972), Morales & O’Neil (Reference Morales and O’Neil1972), Berger et al. (Reference Berger, Brunner, Chapman, Divol, Still and Valeo2013).
6. Phase II: quasi-steady collisional phase
Phase II is the longest-lived phase during which a quasi-steady state between weak electron–electron collisions and strong wave–electron interactions governs the dynamics. Figures 2(a) and 2(c) suggest that the wave evolution during this phase is characterised by three quantities; the effective damping rate of the wave, the growth rate of the nonlinear frequency shift and the lifetime of Phase II. Each of these quantities is discussed in the following subsections. First, the phase space evolution during this phase is described.
6.1. Phase space evolution
Figures 8(a) and 8(b) show phase space snapshots during Phase II for Simulations A and Simulation B, respectively, at
$\omega _{pe}t = 1400$
.
Phase space portrait at
$t=1400\omega _{pe}^{-1}$
, long after the drive is turned off, for two different collision frequencies.

The snapshots show that similarly sized vortices persist in both simulations despite the wave amplitude in Simulation B being roughly half that of Simulation A.
Simulation B exhibits substantial phase space diffusion to the extent that the trapped particle orbits are not discernible as they are in Simulation A. This is because the higher collision frequency in Simulation B leads to more rapid phase space diffusion in the deeply trapped region where the gradients in velocity space are large.
It is also evident from figure 8 that there is a phase shift between the two simulations. This is because the wave frequency has experienced a greater downshift in Simulation B due to the effect highlighted in the previous section.
6.2. Wave amplitude and frequency evolution
Figure 9(a) shows the evolution of the electric field amplitude over a longer time scale for all four simulations. The monotonic increase of the effective damping rate with
$\nu _{\textit{ee}}$
has two complementary physical origins. At the linear level, the Lenard–Bernstein treatment of the Vlasov–Fokker–Planck dispersion relation (Lenard & Bernstein Reference Lenard and Bernstein1958) predicts that the bare collision frequency simply augments the Landau rate,
$\nu _{\text{eff}} \approx \nu _L + \mathcal{O}(\nu _{\textit{ee}})$
, so more collisions mean more damping. In the large-amplitude regime considered here, however, trapping would otherwise largely suppress the Landau contribution (O’Neil Reference O’Neil1965); the Phase II damping rate instead reflects the efficiency with which collisions scatter particles across the separatrix and reinstate a Landau-like kinetic response. The near-linear-in-
$\nu _{\textit{ee}}$
scaling of the effective damping rate in the trapped regime was first derived by Zakharov & Karpman (Reference Zakharov and Karpman1963) via an asymptotic boundary-layer analysis of the separatrix (introduced in § 1; the quantitative comparison with our fits is given later in this section). More recently, Berger et al. (Reference Berger, Arrighi, Chapman, Dimits, Banks and Brunner2023) have recast the same trapping–detrapping competition as a binary amplitude-threshold criterion,
$\omega _B \tau _c \sim 2\pi$
: above threshold, bounce motion wins and
$\nu _{\text{eff}} \ll \nu _L$
, while increasing
$\nu _{\textit{ee}}$
at fixed
$\omega _B$
drives the wave towards the threshold and
$\nu _{\text{eff}}$
accordingly grows towards
$\nu _L$
. The Phase II rates in figure 9(a) sit between
$\nu _{\textit{ee}}$
and
$\nu _L$
for all four simulations, consistent with this picture – neither the bare collision rate nor the linear Landau rate describes the dynamics in isolation. Bounce oscillations have been suppressed over this time scale. Of the four simulations, only Simulation D reaches Phase III within the displayed window: its amplitude in figure 9(a) begins to decrease rapidly after
$\omega _{pe}t \approx 2300$
, marking the onset of the collisional return to Landau damping analysed in § 7. This is consistent with (6.3): Simulation D has the largest
$\nu _{\textit{ee}}$
in the scan and therefore the shortest Phase II lifetime, so it is the first to exit Phase II.
Evolution of the (a) wave amplitude and (b) wave frequency for four different collision frequencies. Simulation D, which has the largest
$\nu _{\textit{ee}}$
and therefore the shortest Phase II lifetime (6.3), enters Phase III near
$\omega _{pe}t \approx 2300$
within the displayed window, visible as the rapid drop in amplitude in panel (a) and the rapid rise in frequency back toward
$\omega _{\text{EPW}}$
in panel (b); see § 7.

Figure 9. Long description
Two line graphs depict the evolution of wave amplitude and wave frequency for different collision frequencies. Panel A: The line graph shows the evolution of wave amplitude over time. The horizontal axis is labeled with omega sub pe t, and the vertical axis is labeled with e times absolute value of E sub 1 divided by m sub e times v sub th times omega sub pe. Four different colored lines represent different collision frequencies. Panel B: The line graph shows the evolution of wave frequency over time. The horizontal axis is labeled with omega sub pe t, and the vertical axis is labeled with omega sub EPW divided by omega sub pe. Four different colored lines represent different collision frequencies. The legend indicates the labels A, B, C, and D for the different lines. The graphs show how the wave amplitude and frequency change over time for each collision frequency.
Figure 9(b) shows the corresponding evolution of the instantaneous frequency of the wave. It shows that the frequency can downshift substantially during Phase II and that the rate of downshift increases with increasing collision frequency. Mirroring the amplitude behaviour in panel (a), the frequency in Simulation D begins to rapidly increase after
$\omega _{pe}t \approx 2300$
and returns towards the linear value
$\omega _{\text{EPW}}$
; this is the Phase III signature of the re-Maxwellianisation of the distribution function, and is analysed in detail in § 7. This is the same Simulation D used throughout the manuscript as the reference three-phase evolution shown in figure 2.
6.3. Physical mechanism for frequency shift enhancement
Figure 10(a) distinguishes the change in the spatially averaged distribution function over the first 600
$\omega _{pe}t$
from the Vlasov terms and the Fokker–Planck terms. It shows that the Vlasov terms redistribute the phase space density from slower velocities to faster velocities while the Fokker–Planck terms do the opposite. The net result, not shown here due to space constraints, is a slight increase in phase space density in the trapping region above the phase velocity.
(a) Change in phase space density
$\Delta |\hat {f}^0(v)|$
over Phase II from the different terms. (b) Corresponding contributions to the susceptibility function
$\Delta \chi (v)$
and the sum. Solid, Sim A; dashed, Sim B.

Figure 10(b) shows the corresponding contribution to the susceptibility function. It shows that the Vlasov terms, in blue, are responsible for a net decrease in the frequency shift, while the Fokker–Planck terms, in red, are responsible for a net increase in the frequency shift. The net result, in green, is an increase in the frequency shift during Phase II.
The dashed lines in figure 10 are the results from Simulation B. It shows similar behaviour to Simulation A, but the magnitudes of the contributions are larger due to the higher collision frequency. A larger increase in phase space density in the trapping region results in a larger increase in the frequency shift.
6.4. Summary of Phase II dynamics
In Phase II, the growth rate of the nonlinear frequency shift, the effective damping rate of the wave and the lifetime of Phase II are given by
respectively. Equations (6.1) and (6.2) vanish as
$\nu _{\textit{ee}} \rightarrow 0$
and the collisionless dynamics of Phase I are recovered. Collisionless behaviour prescribes a persistent wave, with a frequency shift
$\Delta \omega _{\text{NL}}$
due to a trapped particle distribution.
The damping rate in (6.2) has a weaker scaling than that predicted by Zakharov & Karpman (Reference Zakharov and Karpman1963). Their derivation results in a
$\omega _B^{-3}$
scaling, while the best fit finds a weaker scaling of
$\omega _B^{-1.48}$
. Due to this, the ZK damping rate significantly underestimates the observed damping rate for large
$\omega _B$
. This could be because ZK uses a perturbative approach that assumes strong trapping such that
$\omega _B \gg \nu _{\textit{ee}} v_{th}^2/v_{tr}^2$
, which may not hold for weaker wave amplitudes. Similar to the nonlinear frequency shift, it might be that
$\omega _B^{-3}$
scaling would be more accurate for much larger wave amplitudes and at smaller wavenumbers.
Simulations show that the EPW persists in Phase II for
${\gt} \mathcal{O}(10)$
trapping oscillations. The lifetime of the quasi-steady state during Phase II can also be characterised and an empirical fit is provided in (6.3).
7. Phase III: return to Landau damping
Phase III begins when the distribution function has returned to a nearly Maxwellian state and the wave rapidly damps via Landau damping. As shown in figure 2(a), the damping rate increases rapidly to approximately the Landau damping rate (or slightly higher) and figure 2(c) shows that the oscillation frequency rapidly returns to a value near the initial
$\omega _{\text{EPW}}$
.
(a)
$\Delta |\hat {f}^0(v)|$
accumulated through Phase I, Phase II and Phase III. (b) Contribution to
$\Delta \chi (v)$
in phase III from each term.

7.1. Distribution function thermalisation
Figure 11(a) shows the total change in the spatially averaged distribution function over all three phases for Simulation C. During Phases I and II, there is an increase in density near the trapping velocity. During Phase III, this density increase has been erased and the additional phase space density has been redistributed to lower velocities as the distribution function returns to a Maxwellian.
Figure 11(b) shows the contribution to the susceptibility function from each term during Phase III. It shows that the Fokker–Planck term is responsible for returning the distribution function to a Maxwellian, while the effect of the Vlasov term is negligible.
7.2. Wave properties and enhanced damping
Finally, the Phase III dynamics are given by
where
$f_L$
gives the damping rate as a multiple of the Landau damping rate in Phase III. It is seen that with more energy absorbed by the particles during Phase II (i.e. larger
$\omega _B$
), the damping rate during Phase III is larger because the final Maxwellian is slightly warmer.
During Phase III,
$\Delta \epsilon$
relaxes to a value smaller than the value at
$t=0$
, suggesting a small frequency upshift. This upshift is likely due to the heating of the distribution function from absorbing the EPW energy. The return to a Maxwellian distribution is gradual during Phase III as the Fokker–Planck contribution asymptotically approaches the necessary value for the slightly warmer EPW.
8. Conclusion
This work presents a comprehensive kinetic study of large-amplitude electron plasma wave evolution in weakly collisional plasmas using 1600 Vlasov–Poisson–Fokker–Planck simulations spanning a broad parameter space in wavenumber, wave amplitude and collision frequency. The simulations reveal that the wave dynamics naturally partition into three distinct phases, each governed by different physical mechanisms and amenable to different perturbative descriptions.
Phase I encompasses the well-understood collisionless dynamics of wave excitation, particle trapping and quasi-linear plateau formation. During this short-lived phase, the distribution function flattens at the phase velocity, wave damping is suppressed and a nonlinear frequency downshift emerges consistent with previous theoretical predictions, though with a weaker amplitude scaling (
$\Delta \omega \propto \omega _B^{0.82}$
) than the linear dependence predicted by asymptotic theory.
The central finding of this work concerns Phase II, during which a quasi-steady equilibrium develops between weak electron–electron collisions and strong wave–electron interactions. Counter-intuitively, collisions are found to enhance the nonlinear frequency shift rather than simply restore linear behaviour. This effect arises because the Fokker–Planck collision operator contributes positively to the susceptibility function, while the Vlasov terms only partially counteract this contribution. The frequency shift can approximately double during Phase II relative to its Phase I value. Empirical scaling relations for the initial frequency shift (5.1), frequency shift growth rate (6.1), effective damping rate (6.2) and phase lifetime (6.3) are provided, enabling the construction of reduced models for weakly collisional EPW dynamics.
A comparison of the damping rate in Phase II with Zakharov–Karpman theory reveals that while the linear scaling with collision frequency is confirmed, the theory underpredicts the damping rate by approximately an order of magnitude and predicts a stronger inverse dependence on bounce frequency (
$\omega _B^{-3}$
) than is observed (
$\omega _B^{-1.48}$
).
Phase III marks the return to linear behaviour as collisions restore a near-Maxwellian distribution and the wave undergoes rapid Landau damping. The transition occurs after a time,
$\tau _{II}$
given by (6.3), and the damping rate during this phase is characterised by (7.1)
Because this study is one-dimensional, it cannot capture transverse escape of trapped electrons from finite-width laser hot-spots (‘sideloss’), often the dominant loss channel in inertial confinement fusion (ICF) geometries. Sideloss is commonly modelled as a Krook relaxation at rate
$\nu _{\text{sl}} \approx v_{\text{th}} / L_\perp$
(Strozzi Reference Strozzi2005). For an NIF-like single speckle (
$L_\perp \approx 2.8\,\mu$
m,
$f/8$
,
$3\omega$
) in plasma with
$T_e = 3$
keV and
$n_e/n_c = 0.1$
,
$\nu _{\text{sl}}/\omega _{pe} \approx 5 \times 10^{-3}$
, roughly an order of magnitude above the Phase II rate
$\nu _{\text{eff}}/\omega _{pe} \sim 6 \times 10^{-4}$
from (6.2) at representative parameters (
$k\lambda _D = 0.3$
,
$\omega _B = 0.15$
,
$\nu _{\textit{ee}} = 10^{-4}$
), and also above the threshold
${\sim}2 \times 10^{-3}\,\omega _{pe}$
at which Strozzi (Reference Strozzi2005) observed sharp disruption of trapping. At the single-speckle scale, sideloss is therefore expected to dominate, and the 1-D Phase II time scales reported here would overestimate the persistence of the nonlinear state. The two rates become comparable for
$L_\perp \gtrsim 20\,\unicode{x03BC}$
m (a few speckle widths), and other geometric or kinetic loss channels may dominate in space, ionospheric and discharge plasmas. The two mechanisms act through different phase-space channels. Krook sideloss damps deviations from a Maxwellian indiscriminately, whereas Fokker–Planck collisions preferentially diffuse the sharp separatrix gradients. A quantitative assessment of their competition is left to future work.
The collisional modification of the EPW dispersion relation reported here parallels recent findings for ion-acoustic waves (IAWs). For linear IAWs in laser-heated, non-Maxwellian plasmas, Capdessus et al. (Reference Capdessus, Ruyer, Debayle, Loiseau and Masson-Laborde2025) showed that binary Coulomb collisions modify both the real frequency and the damping rate. The Langdon-effect contribution and the collision-induced anisotropy contribution act in opposing directions, yielding damping rates either above or below the collisionless Landau rate. The Phase II results here add a nonlinear counterpart: the Fokker–Planck contribution to the susceptibility enhances the nonlinear frequency shift and produces an effective damping rate well above
$\nu _{\textit{ee}}$
. Together, these results suggest that dispersion in weakly collisional plasmas cannot be captured by the collisionless dispersion relation alone, nor by treating collisions as a simple relaxation towards a Maxwellian.
Future work should extend this analysis to explore the role of ion dynamics and investigate the interaction of multiple EPW modes in the weakly collisional regime. Additionally, developing analytic theory that captures the enhanced frequency shift mechanism identified here would provide deeper physical insight into the quasi-steady collisional dynamics of Phase II.
Acknowledgements
A.J. thanks useful discussions with B. Afeyan to motivate this study.
Editor Antoine C. Bret thanks the referees for their advice in evaluating this article.
Funding
This material is based upon work supported by IFE COLoR under U.S. Department of Energy Grant No. DE-SC0024863, US DOE National Nuclear Security Administration (NNSA) Center of Excellence under Cooperative Agreement No. DE-NA0003869, and by the Department of Energy National Nuclear Security Administration under Award Number DE-NA0004144, the University of Rochester, and the New York State Energy Research and Development Authority. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award FES-ERCAP0026741.
Declaration of interests
The authors report no conflict of interest.



m
E
n
νee
sin(kx−ωEPWt)
vph
Q=C(kλD)α(ωB)β(νee)γ
C,α,β,γ
<0.005

kλD=0.3
a0
νee
t=650ωpe−1

Δ|f^0(v)|
ωpet=600
Δχ(v)
Δ|f^0|
t=1400ωpe−1
νee
ωpet≈2300
ωEPW
Δ|f^0(v)|
Δχ(v)
Δ|f^0(v)|
Δχ(v)