2.1 Julia implementation

using FFTW, Plots, LinearAlgebra
using BenchmarkTools, Statistics

• FFTW: Fast Fourier Transform library for efficient frequency-domain operations
• Plots: Visualization of results
• LinearAlgebra: Linear algebra utilities
• BenchmarkTools: Performance benchmarking
• Statistics: Statistical functions

2.1.1 Mesh Data Structure

"""
    UniformMesh(start, stop, length)

    1D uniform mesh data structure for storing discretization information.
"""
struct UniformMesh

   start    :: Float64      # Starting point of the domain
   stop     :: Float64      # Ending point of the domain
   length   :: Int64        # Number of grid points
   step     :: Float64      # Grid spacing (Δx or Δv)
   points   :: Vector{Float64}  # Array of grid point coordinates

   function UniformMesh(start, stop, length::Int; endpoint=true)

       # Create uniform grid from start to stop with 'length' points
       # endpoint=false means we exclude the right boundary (periodic domain)
       points = LinRange(start, stop, length+1)[1:end-1]
       
       # Uniform spacing
       step = (stop - start) / length

       new(start, stop, length, step, points)

   end

end

Purpose: This data structure encapsulates all the information needed to describe a 1D uniform grid, making the code cleaner and more organized.

2.1.2 Function to Compute Charge Density

"""

    compute_rho(meshv::UniformMesh, f)

    Compute the charge density from the distribution function:
    ρ(x,t) = ∫ f(x,v,t) dv

    We subtract the mean to remove the constant background charge,
    keeping only the oscillating part that drives the dynamics.

"""
function compute_rho(meshv::UniformMesh, f)

   dv  = meshv.step                      # Velocity grid spacing
   ρ = dv .* vec(sum(real(f), dims=2))   # Integrate f over all velocities
   ρ .- mean(ρ)                          # Subtract the mean
   
end

• We numerically integrate the distribution function over velocity space using the trapezoidal rule (represented by multiplying by \(dv\))
• Subtracting the mean gives us the deviation from the background charge, which is what creates the net electric field
• dims=2 specifies we sum over the velocity dimension (columns), leaving the spatial dimension (rows)

2.1.3 Function to Compute Electric Field

"""
    compute_e(mesh::UniformMesh, ρ)

    Compute the electric field from the charge density using Poisson's equation.
    
    In Fourier space: ∂²E/∂x² = -ρ/ε₀ becomes (-k²) Ê = -ρ̂/ε₀
    
    Solution: Ê_k = ρ̂_k / (i k ε₀)
    
"""
function compute_e(mesh::UniformMesh, ρ)

   n = mesh.length
   k =  2π / (mesh.stop - mesh.start)    # Fundamental wavenumber
   
   # Wavenumber array for FFT output ordering
   modes  = zeros(Float64, n)
   modes .= k * vcat(0:n÷2-1, -n÷2:-1)  # [0, 1, 2, ..., -n/2, ..., -1] × k
   
   modes[1] = 1.0  # Avoid division by zero for the zero mode
   
   ρ̂ = fft(ρ)./modes                    # Fourier transform and divide by k
   
   vec(real(ifft(-1im .* ρ̂)))            # Inverse transform with factor of -i
   
end

1. Fourier transform the charge density: \(\hat{\rho}_k\)
2. For each mode \(k\), solve Poisson equation in Fourier space: \(\hat{E}_k = -i\hat{\rho}_k / k\)
3. The factor \(-i\) comes from integrating \(E = -\partial_x \phi\) in Fourier space
4. Inverse Fourier transform to get \(E(x)\) in physical space

2.1.4 Structure for Advection Operations

"""

    advection! = AmpereAdvection(mesh) 

    This structure encapsulates the advection operators for the Vlasov–Ampere system:
    
    - Advection in velocity: ∂f/∂t + E(x) ∂f/∂v = 0
    - Advection in space:    ∂f/∂t + v ∂f/∂x = 0
    - Field update:          ∂E/∂t = -J = -∫ fv dv

    It precomputes and stores the wavenumber array for efficient Fourier operations.

"""
struct AmpereAdvection 
    
    mesh :: UniformMesh
    kx   :: Vector{Float64}   # Wavenumber array

    function AmpereAdvection(mesh)
        
        nx  = mesh.length
        dx  = mesh.step
        Lx  = mesh.stop - mesh.start
        kx  = zeros(Float64, nx)
        
        # Standard FFT wavenumber ordering: [0, 1, 2, ..., -n/2, ..., -2, -1]
        kx .= 2π/Lx .* fftfreq(nx, nx)
        
        new(mesh, kx)
      
    end

end

Why Store Wavenumbers? Pre-computing the wavenumber array avoids redundant calculations in the inner loop of the simulation.

2.1.5 Advection Function in Velocity Space

"""
    Advection function along velocity (v-direction)

    Solves: ∂f/∂t = E ∂f/∂v
    
    In Fourier space in x, this becomes:
    ∂f̂_k/∂t = E_k ∂f̂_k/∂v
    
    Semi-implicit solution:
    f̂_k^(n+1)(v) = exp(-i dt E_k / k_x) × f̂_k^n(v)
    
    The exponential factor represents how the electric field
    shifts the phase of each Fourier component as a function of velocity.

"""
function (adv::AmpereAdvection)(fᵗ::Array{ComplexF64,2}, 
                                e::Vector{ComplexF64}, 
                                dt::Float64)
    
    # Transform to Fourier space in x (rows are spatial modes, columns are velocity)
    fft!(fᵗ, 1)
    
    # Apply the advection operator: multiply by phase factor
    # -1im * dt * adv.kx gives the exponent for each wavenumber
    # transpose(e) broadcasts the electric field (one value per mode) to all velocities
    fᵗ .= fᵗ .* exp.(-1im * dt * adv.kx * transpose(e))
    
    # Transform back to physical space
    ifft!(fᵗ, 1)

end

• In Fourier space, the operator \(\frac{\partial}{\partial x}\) becomes multiplication by \(ik\)
• The advection equation \(\frac{\partial f}{\partial t} = -v\frac{\partial f}{\partial x}\) in Fourier space becomes \(\frac{\partial \hat{f}_k}{\partial t} = -v \cdot ik \hat{f}_k\)
• The exact solution is \(\hat{f}_k(t) = e^{-ik v t} \hat{f}_k(0)\)

2.1.6 Advection Function in Space

"""
    Advection function along space (x-direction) and field computation

    Solves the coupled system:
    ∂f/∂t = v ∂f/∂x
    ∂E/∂t = -J = -∫ fv dv
    
    Steps:
    1. Transform f to Fourier space in x
    2. Apply spatial advection: f̂_k → exp(-i k v dt) f̂_k
    3. Compute charge density from the advected distribution
    4. Compute electric field from charge density using Poisson equation
    5. Transform everything back to physical space

"""
function (adv::AmpereAdvection)(f::Array{ComplexF64,2}, 
                                e::Vector{ComplexF64}, 
                                v::Vector{Float64}, 
                                dt::Float64)
    
    # Prepare the phase shift factor for each (wavenumber, velocity) pair
    ev = exp.(-1im * dt * adv.kx * transpose(v))    
    
    # Transform f to Fourier space in x
    fft!(f, 1)
    
    # Apply spatial advection: multiply by the phase shift
    f .= f .* ev
    
    # Compute charge density by integrating over velocity
    dv = v[2] - v[1]                    # Velocity spacing
    ρ = dv * vec(sum(f, dims=2))        # ρ_k = Δv ∑_j f̂_k(v_j)
    
    # Solve Poisson equation for each Fourier mode (except k=0)
    for i in 2:length(e)
        e[i] = -1im * ρ[i] ./ adv.kx[i]
    end
    
    # Zero mode remains unchanged
    e[1] = 0.0
    
    # Transform the distribution function back to physical space
    ifft!(f, 1)
    
    # Transform the electric field back to physical space
    ifft!(e)
    
    # The electric field should be real; take the real part
    # (imaginary parts are numerical errors)
    e .= real(e)

end

• The phase shift \(e^{-ikt v}\) encodes how each Fourier mode of the distribution shifts due to spatial advection
• After advection, we immediately compute the charge density and field self-consistently
• This ensures the field responds to the updated distribution function

2.1.7 Initialize the Distribution Function: Landau Damping

"""
    landau(ϵ, kx, x, v)

    Landau damping initialization function
    
    f(x,v) = (1/√(2π)) × (1 + ϵ cos(kₓ x)) × exp(-v²/2)
    
    This represents a perturbed Maxwell-Boltzmann distribution:
    - Base: Gaussian distribution in velocity (thermal equilibrium)
    - Perturbation: Sinusoidal density modulation in space
    
    Parameters:
    - ϵ: Perturbation amplitude (typically small, ~0.001)
    - kx: Wavenumber of the perturbation
    - x: Spatial grid points
    - v: Velocity grid points
    
    Physical Phenomenon:
    In collisionless plasma, this small density perturbation excites
    electrostatic waves. Due to wave-particle interactions (Landau damping),
    the wave amplitude decays exponentially even without collisions,
    with energy transferred to particle heating.
    
    See: https://en.wikipedia.org/wiki/Landau_damping

"""
function landau(ϵ, kx, x, v)
    
    # Create a 2D distribution: rows are spatial points, columns are velocity points
    # cos(kx*x) is a column vector, exp(-v²/2) is a row vector
    # Multiplying them together broadcasts to a 2D array
    (1.0 .+ ϵ * cos.(kx * x)) / sqrt(2π) .* transpose(exp.(-0.5 * v .* v))
    
end

• The Gaussian in velocity represents a thermalized plasma at rest
• The cosine modulation in space creates a density wave
• This wave couples to the plasma and creates oscillations
• Landau damping causes these oscillations to decay logarithmically: \(E(t) \propto e^{-\gamma t}\) with \(\gamma > 0\)