A New Approach to Compute Linear Landau Damping

Problem Statement
The paper addresses the challenge of obtaining exact time-domain solutions for linear Landau damping problems to serve as benchmarks for numerical simulations in regimes where analytical solutions are intractable. Traditional approaches rely on residue-based analyses, which involve computing the complex zeros (poles) of the dielectric function. While effective for late-time asymptotic behavior governed by dominant poles, these methods struggle with the transient response near the initial time. Resolving this transient phase requires an increasingly large number of poles, and determining them numerically can be difficult or impossible if the dielectric function cannot be analytically continued sufficiently far into the lower complex half-plane. Furthermore, standard Fourier inversion of the response function often suffers from poor convergence due to algebraic tails ($\omega^{-1}$) arising from discontinuities at $t=0$.

Methodology
The authors present a semi-analytical method that circumvents the need for pole summation entirely. The core of the methodology involves:

1. Time-Symmetric Decomposition: Instead of solving for a causal, one-sided response starting at $t=0$, the method constructs the solution as a superposition of a forward-evolving component ($f^+_1$, defined for $t \ge 0$) and a backward-evolving component ($f^-_1$, defined for $t \le 0$). Both are uniquely determined by the same initial condition $f_1(x, v, t=0)$.
2. Symmetrization: By combining these components, the method effectively defines a "prehistory" for the distribution function that evolves from $t=-\infty$ to the prescribed initial state at $t=0$. This symmetrization removes the discontinuity (jump) at $t=0$ inherent in one-sided responses.
3. Spectral Properties: The removal of the $t=0$ jump eliminates the $\omega^{-1}$ algebraic tail in the Fourier spectrum. Consequently, the combined spectrum decays super-algebraically (faster than any inverse power of $\omega$), provided the initial condition admits a smooth solution with compact support. This rapid decay ensures that truncation errors in numerical Fourier inversion become negligible outside a finite integration domain.
4. Handling Singularities: For the distribution function response, the method explicitly handles the resonant singularity on the real axis associated with the continuous spectrum of Case-van Kampen modes. This is achieved by applying the Sokhotski-Plemelj theorem to separate the delta-function contribution (evaluated analytically) from the principal-value integral (evaluated numerically using a regularized integrand).
5. Generalization: The authors provide a formal argument extending this logic to linear Vlasov-Maxwell problems with a uniform background magnetic field and Maxwellian equilibrium distributions, demonstrating that the cancellation of source terms and initial-field terms in the frequency domain remains valid.

Key Results
The method was applied to an electrostatic, unmagnetized plasma with kinetic ions and adiabatic electrons, deriving the response of both ion density and the ion distribution function to an initial density perturbation.

• Density Response: The symmetrized density response allows for accurate numerical Fourier inversion. The authors demonstrate that the error is dominated by the truncation of the integration domain rather than quadrature errors, and they provide explicit bounds for these errors. The solution accurately recovers the initial condition at $t=0$ without the Gibbs oscillations seen in non-symmetrized approaches.
• Distribution Function Response: The method successfully resolves the continuous spectrum. The results show that at fixed time, the perturbation vanishes at large velocities due to Gaussian suppression. At fixed velocity, the solution transitions from a perturbed state with amplitude modulations to an asymptotic free-streaming solution after the Landau-damped terms decay.
• Energy Conservation: The solution is verified to be physically consistent by examining the ion free-energy density. The analysis confirms that for the time-symmetric initial condition, free energy is transferred between ions and the adiabatic electron response, with the total free energy conserved.
• Benchmarking: The semi-analytical solution was compared against a fully kinetic 1D1V simulation (BSL6D code). The pointwise relative error was found to be low, with deviations primarily attributed to interpolation errors in the numerical simulation at the edges of the velocity domain.
• Two-Species Extension: An illustrative example of a self-consistent two-species Vlasov-Poisson plasma (electrons and ions) was presented, capturing the coupling between high-frequency Langmuir waves and low-frequency ion-acoustic oscillations.

Significance and Claims
The paper claims that its primary significance lies in providing a robust, semi-analytical framework for calculating exact time-domain solutions without the computational burden and numerical instability associated with finding complex zeros of dielectric functions. By constructing a time-symmetric solution, the method transforms a problem with slow spectral convergence into one with rapid (super-algebraic) convergence, making high-precision numerical Fourier inversion feasible.

The authors emphasize that this approach is not merely a theoretical exercise but a practical tool for validating numerical codes in regimes where indirect diagnostics or simplified limits are insufficient. They explicitly state that while the explicit algebraic construction for magnetized configurations involves significant complexity and is left for future work, the formal argument confirms that the numerical method is, in principle, extensible to such cases provided the background magnetic field is uniform and the equilibrium distribution is Maxwellian. The work does not propose new experimental setups but rather offers a rigorous benchmarking standard for existing and future plasma simulation codes.