diff --git a/analysis/check_distribution_moments.py b/analysis/check_distribution_moments.py
new file mode 100644
index 000000000..123e36389
--- /dev/null
+++ b/analysis/check_distribution_moments.py
@@ -0,0 +1,165 @@
+import os
+import numpy as np
+
+BASE = (
+    "thesis_postprocess_validation/post_processing/"
+    "kinetic_data/hot_elec/distribution_function/v1_v3_density"
+)
+
+TIME_PATH = (
+    "thesis_postprocess_validation/post_processing/t_grid.npy"
+)
+
+t = np.load(TIME_PATH)
+
+v1 = np.load(os.path.join(BASE, "grid_v1.npy"))
+v3 = np.load(os.path.join(BASE, "grid_v3.npy"))
+F = np.load(os.path.join(BASE, "f_binned.npy"))
+
+dv1 = np.mean(np.diff(v1))
+dv3 = np.mean(np.diff(v3))
+
+
+def moments(Fi):
+
+    norm = np.sum(Fi) * dv1 * dv3
+
+    mean1 = (
+        np.sum(v1[:, None] * Fi)
+        * dv1 * dv3 / norm
+    )
+
+    mean3 = (
+        np.sum(v3[None, :] * Fi)
+        * dv1 * dv3 / norm
+    )
+
+    var1 = (
+        np.sum(
+            (v1[:, None] - mean1) ** 2 * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    var3 = (
+        np.sum(
+            (v3[None, :] - mean3) ** 2 * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    covariance = (
+        np.sum(
+            (v1[:, None] - mean1)
+            * (v3[None, :] - mean3)
+            * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    return (
+        norm,
+        mean1,
+        mean3,
+        np.sqrt(var1),
+        np.sqrt(var3),
+        covariance,
+    )
+
+
+results = np.array([moments(Fi) for Fi in F])
+
+norm = results[:, 0]
+mean1 = results[:, 1]
+mean3 = results[:, 2]
+sigma1 = results[:, 3]
+sigma3 = results[:, 4]
+covariance = results[:, 5]
+
+anisotropy = sigma1**2 / sigma3**2
+
+relative_change = np.array(
+    [
+        np.linalg.norm(Fi - F[0])
+        / np.linalg.norm(F[0])
+        for Fi in F
+    ]
+)
+
+marginal_v1_initial = np.sum(F[0], axis=1) * dv3
+marginal_v1_final = np.sum(F[-1], axis=1) * dv3
+
+marginal_v3_initial = np.sum(F[0], axis=0) * dv1
+marginal_v3_final = np.sum(F[-1], axis=0) * dv1
+
+v1_change = (
+    np.linalg.norm(
+        marginal_v1_final - marginal_v1_initial
+    )
+    / np.linalg.norm(marginal_v1_initial)
+)
+
+v3_change = (
+    np.linalg.norm(
+        marginal_v3_final - marginal_v3_initial
+    )
+    / np.linalg.norm(marginal_v3_initial)
+)
+
+print("============================================================")
+print("DISTRIBUTION MOMENT ANALYSIS")
+print("============================================================")
+
+print("\nInitial values:")
+print("integral =", norm[0])
+print("sigma_v1 =", sigma1[0])
+print("sigma_v3 =", sigma3[0])
+print("anisotropy =", anisotropy[0])
+print("covariance =", covariance[0])
+
+print("\nFinal values:")
+print("integral =", norm[-1])
+print("sigma_v1 =", sigma1[-1])
+print("sigma_v3 =", sigma3[-1])
+print("anisotropy =", anisotropy[-1])
+print("covariance =", covariance[-1])
+
+print("\nMarginal changes:")
+print("v1 relative L2 change =", v1_change)
+print("v3 relative L2 change =", v3_change)
+
+print("\nFull 2D distribution change:")
+print("initial -> final =", relative_change[-1])
+print("maximum =", np.max(relative_change))
+print(
+    "time of maximum change =",
+    t[np.argmax(relative_change)],
+)
+
+print("\nRanges over complete run:")
+
+print(
+    "sigma_v1:",
+    np.min(sigma1),
+    np.max(sigma1),
+)
+
+print(
+    "sigma_v3:",
+    np.min(sigma3),
+    np.max(sigma3),
+)
+
+print(
+    "anisotropy:",
+    np.min(anisotropy),
+    np.max(anisotropy),
+)
+
+print(
+    "covariance:",
+    np.min(covariance),
+    np.max(covariance),
+)
+
+print("\nAnalysis complete.")
diff --git a/analysis/plot_e3_density.py b/analysis/plot_e3_density.py
new file mode 100644
index 000000000..baf5b5d60
--- /dev/null
+++ b/analysis/plot_e3_density.py
@@ -0,0 +1,81 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+BASE = (
+    "thesis_postprocess_validation/post_processing/"
+    "kinetic_data/hot_elec/distribution_function/e3_density"
+)
+
+TIME_PATH = (
+    "thesis_postprocess_validation/post_processing/t_grid.npy"
+)
+
+OUT = "validation_plots"
+os.makedirs(OUT, exist_ok=True)
+
+t = np.load(TIME_PATH)
+
+e3 = np.load(os.path.join(BASE, "grid_e3.npy"))
+F = np.load(os.path.join(BASE, "f_binned.npy"))
+DF = np.load(os.path.join(BASE, "delta_f_binned.npy"))
+
+de3 = np.mean(np.diff(e3))
+
+integrals = np.sum(F, axis=1) * de3
+
+print("============================================================")
+print("E3 SPATIAL DENSITY VALIDATION")
+print("============================================================")
+
+print("grid shape =", e3.shape)
+print("distribution shape =", F.shape)
+
+print("\ne3 range:")
+print(np.min(e3), np.max(e3))
+
+print("\nIntegral:")
+print("initial =", integrals[0])
+print("final =", integrals[-1])
+print("minimum =", np.min(integrals))
+print("maximum =", np.max(integrals))
+
+relative_variation = np.max(
+    np.abs(integrals - integrals[0])
+) / abs(integrals[0])
+
+print("maximum relative variation =", relative_variation)
+
+print(
+    "\nDelta-f identically zero =",
+    np.all(DF == 0),
+)
+
+indices = [0, len(t) // 2, len(t) - 1]
+
+plt.figure(figsize=(10, 6))
+
+for i in indices:
+    plt.plot(
+        e3,
+        F[i],
+        label=f"t = {t[i]:.3f}",
+    )
+
+plt.xlabel(r"$e_3$")
+plt.ylabel("f")
+plt.title("Hot-Electron Spatial Density Distribution")
+plt.legend()
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "e3_density_selected_times.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("\nSaved:", path)
diff --git a/analysis/plot_energy_diagnostics.py b/analysis/plot_energy_diagnostics.py
new file mode 100644
index 000000000..45c9bf463
--- /dev/null
+++ b/analysis/plot_energy_diagnostics.py
@@ -0,0 +1,105 @@
+import os
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+
+DATA = (
+    "thesis_postprocess_validation/data/data_proc0.hdf5"
+)
+
+OUT = "validation_plots"
+os.makedirs(OUT, exist_ok=True)
+
+with h5py.File(DATA, "r") as f:
+
+    t = f["time/value"][:]
+
+    en_B = f["scalar/en_B"][:]
+    en_E = f["scalar/en_E"][:]
+    en_J = f["scalar/en_J"][:]
+    en_f = f["scalar/en_f"][:]
+    en_tot = f["scalar/en_tot"][:]
+
+relative_energy_error = (
+    (en_tot - en_tot[0])
+    / en_tot[0]
+)
+
+print("============================================================")
+print("ENERGY DIAGNOSTICS")
+print("============================================================")
+
+print("number of states =", len(t))
+print("time range =", t[0], t[-1])
+
+print("\nInitial energies:")
+print("E_B =", en_B[0])
+print("E_E =", en_E[0])
+print("E_J =", en_J[0])
+print("E_f =", en_f[0])
+print("E_total =", en_tot[0])
+
+print("\nFinal energies:")
+print("E_B =", en_B[-1])
+print("E_E =", en_E[-1])
+print("E_J =", en_J[-1])
+print("E_f =", en_f[-1])
+print("E_total =", en_tot[-1])
+
+print("\nEnergy conservation:")
+print(
+    "final relative change =",
+    relative_energy_error[-1],
+)
+
+print(
+    "maximum absolute relative change =",
+    np.max(np.abs(relative_energy_error)),
+)
+
+plt.figure(figsize=(10, 6))
+
+plt.plot(t, en_B, label=r"$E_B$")
+plt.plot(t, en_E, label=r"$E_E$")
+plt.plot(t, en_J, label=r"$E_J$")
+plt.plot(t, en_f, label=r"$E_f$")
+
+plt.xlabel("Time")
+plt.ylabel("Energy")
+plt.title("Energy Evolution — Validation Run")
+plt.legend()
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "energy_components.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("\nSaved:", path)
+
+
+plt.figure(figsize=(10, 6))
+
+plt.plot(t, relative_energy_error)
+
+plt.xlabel("Time")
+plt.ylabel("Relative total-energy change")
+plt.title("Total Energy Conservation")
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "total_energy_error.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("Saved:", path)
diff --git a/analysis/plot_v1_v3_anisotropy.py b/analysis/plot_v1_v3_anisotropy.py
new file mode 100644
index 000000000..a669886ea
--- /dev/null
+++ b/analysis/plot_v1_v3_anisotropy.py
@@ -0,0 +1,255 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+BASE = (
+    "thesis_postprocess_validation/post_processing/"
+    "kinetic_data/hot_elec/distribution_function/v1_v3_density"
+)
+
+TIME_PATH = (
+    "thesis_postprocess_validation/post_processing/t_grid.npy"
+)
+
+OUT = "validation_plots"
+os.makedirs(OUT, exist_ok=True)
+
+t = np.load(TIME_PATH)
+
+v1 = np.load(os.path.join(BASE, "grid_v1.npy"))
+v3 = np.load(os.path.join(BASE, "grid_v3.npy"))
+F = np.load(os.path.join(BASE, "f_binned.npy"))
+
+dv1 = np.mean(np.diff(v1))
+dv3 = np.mean(np.diff(v3))
+
+
+def moments(Fi):
+
+    norm = np.sum(Fi) * dv1 * dv3
+
+    mean1 = (
+        np.sum(v1[:, None] * Fi)
+        * dv1 * dv3 / norm
+    )
+
+    mean3 = (
+        np.sum(v3[None, :] * Fi)
+        * dv1 * dv3 / norm
+    )
+
+    var1 = (
+        np.sum(
+            (v1[:, None] - mean1) ** 2 * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    var3 = (
+        np.sum(
+            (v3[None, :] - mean3) ** 2 * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    covariance = (
+        np.sum(
+            (v1[:, None] - mean1)
+            * (v3[None, :] - mean3)
+            * Fi
+        )
+        * dv1 * dv3 / norm
+    )
+
+    return (
+        norm,
+        mean1,
+        mean3,
+        np.sqrt(var1),
+        np.sqrt(var3),
+        covariance,
+    )
+
+
+print("============================================================")
+print("V1-V3 ANISOTROPY VALIDATION")
+print("============================================================")
+
+for index in [0, len(t) // 2, len(t) - 1]:
+
+    (
+        norm,
+        mean1,
+        mean3,
+        sigma1,
+        sigma3,
+        covariance,
+    ) = moments(F[index])
+
+    print(f"\nt = {t[index]:.6f}")
+
+    print("integral =", norm)
+    print("mean_v1 =", mean1)
+    print("mean_v3 =", mean3)
+    print("sigma_v1 =", sigma1)
+    print("sigma_v3 =", sigma3)
+    print("variance ratio =", sigma1**2 / sigma3**2)
+    print("covariance =", covariance)
+
+
+def save_distribution(index, filename, title):
+
+    plt.figure(figsize=(9, 7))
+
+    plt.pcolormesh(
+        v3,
+        v1,
+        F[index],
+        shading="auto",
+    )
+
+    plt.xlabel(r"$v_3$")
+    plt.ylabel(r"$v_1$")
+    plt.title(title)
+
+    plt.colorbar(label="Distribution density")
+
+    plt.tight_layout()
+
+    path = os.path.join(OUT, filename)
+
+    plt.savefig(path, dpi=200)
+    plt.close()
+
+    print("Saved:", path)
+
+
+save_distribution(
+    0,
+    "v1_v3_distribution_t0.png",
+    "Hot-Electron v1-v3 Distribution at t = 0",
+)
+
+save_distribution(
+    -1,
+    "v1_v3_distribution_t1.png",
+    "Hot-Electron v1-v3 Distribution at Final Time",
+)
+
+# ------------------------------------------------------------
+# COMPUTE MARGINAL DISTRIBUTIONS
+# ------------------------------------------------------------
+
+# Initial distribution
+F0 = F[0]
+
+# Final distribution
+F1 = F[-1]
+
+# Integrate over the opposite coordinate
+#
+# f(v1) = ∫ f(v1,v3) dv3
+# f(v3) = ∫ f(v1,v3) dv1
+
+f0_v1 = np.sum(F0, axis=1) * dv3
+f0_v3 = np.sum(F0, axis=0) * dv1
+
+f1_v1 = np.sum(F1, axis=1) * dv3
+f1_v3 = np.sum(F1, axis=0) * dv1
+
+# Normalize only for visualization
+f0_v1 /= np.max(f0_v1)
+f0_v3 /= np.max(f0_v3)
+
+f1_v1 /= np.max(f1_v1)
+f1_v3 /= np.max(f1_v3)
+
+# ------------------------------------------------------------
+# MOMENTS
+# ------------------------------------------------------------
+
+sigma_v1_0 = np.sqrt(
+    np.sum((v1**2) * f0_v1) * dv1 /
+    np.sum(f0_v1 * dv1)
+)
+
+sigma_v3_0 = np.sqrt(
+    np.sum((v3**2) * f0_v3) * dv3 /
+    np.sum(f0_v3 * dv3)
+)
+
+sigma_v1_1 = np.sqrt(
+    np.sum((v1**2) * f1_v1) * dv1 /
+    np.sum(f1_v1 * dv1)
+)
+
+sigma_v3_1 = np.sqrt(
+    np.sum((v3**2) * f1_v3) * dv3 /
+    np.sum(f1_v3 * dv3)
+)
+
+# ------------------------------------------------------------
+# PLOT
+# ------------------------------------------------------------
+
+fig, ax = plt.subplots(
+    1,
+    2,
+    figsize=(14,5),
+    sharey=True
+)
+
+# Initial
+
+ax[0].plot(
+    v1,
+    f0_v1,
+    lw=3,
+    label=rf"$v_1$ ($\sigma$={sigma_v1_0:.3f})"
+)
+
+ax[0].plot(
+    v3,
+    f0_v3,
+    lw=3,
+    label=rf"$v_3$ ($\sigma$={sigma_v3_0:.3f})"
+)
+
+ax[0].set_title("Initial Distribution ($t=0$)")
+ax[0].set_xlabel("Velocity")
+ax[0].set_ylabel("Normalized distribution")
+ax[0].grid(True, alpha=0.3)
+ax[0].legend()
+
+# Final
+
+ax[1].plot(
+    v1,
+    f1_v1,
+    lw=3,
+    label=rf"$v_1$ ($\sigma$={sigma_v1_1:.3f})"
+)
+
+ax[1].plot(
+    v3,
+    f1_v3,
+    lw=3,
+    label=rf"$v_3$ ($\sigma$={sigma_v3_1:.3f})"
+)
+
+ax[1].set_title("Final Distribution ($t=1$)")
+ax[1].set_xlabel("Velocity")
+ax[1].grid(True, alpha=0.3)
+ax[1].legend()
+
+plt.tight_layout()
+
+plt.savefig(
+    "validation_plots/v1_v3_distribution_t0_t1.png",
+    dpi=300,
+    bbox_inches="tight",
+)
+
+print()
+print("Saved:")
+print("validation_plots/v1_v3_distribution_t0_t1.png")
diff --git a/analysis/validate_magnetic_field.py b/analysis/validate_magnetic_field.py
new file mode 100644
index 000000000..9d22751b1
--- /dev/null
+++ b/analysis/validate_magnetic_field.py
@@ -0,0 +1,103 @@
+import os
+import pickle
+import numpy as np
+import matplotlib.pyplot as plt
+
+BASE = "thesis_postprocess_validation/post_processing"
+FIELDS = os.path.join(BASE, "fields_data")
+OUT = "validation_plots"
+
+os.makedirs(OUT, exist_ok=True)
+
+with open(os.path.join(FIELDS, "grids_phy.bin"), "rb") as f:
+    grids = pickle.load(f)
+
+with open(
+    os.path.join(FIELDS, "em_fields/b_field_phy.bin"), "rb"
+) as f:
+    B = pickle.load(f)
+
+z = np.asarray(grids[2])[0, 0, :]
+
+times = sorted(B.keys())
+t0 = times[0]
+
+ix = np.asarray(grids[0]).shape[0] // 2
+iy = np.asarray(grids[1]).shape[1] // 2
+
+Bx = np.asarray(B[t0])[0, ix, iy, :]
+By = np.asarray(B[t0])[1, ix, iy, :]
+Bz = np.asarray(B[t0])[2, ix, iy, :]
+
+target = 1e-4 * np.sin(2 * z)
+
+max_abs_error = np.max(np.abs(Bx - target))
+relative_l2 = np.linalg.norm(Bx - target) / np.linalg.norm(target)
+
+basis = np.sin(2 * z)
+A_fit = np.dot(Bx, basis) / np.dot(basis, basis)
+
+fit = A_fit * basis
+
+relative_residual = (
+    np.linalg.norm(Bx - fit) /
+    np.linalg.norm(Bx)
+)
+
+print("============================================================")
+print("INITIAL MAGNETIC FIELD VALIDATION")
+print("============================================================")
+
+print("time =", t0)
+
+print("\nz range:")
+print(np.min(z), np.max(z))
+
+print("\nBx:")
+print("min =", np.min(Bx))
+print("max =", np.max(Bx))
+print("max abs =", np.max(np.abs(Bx)))
+
+print("\nBy max abs =", np.max(np.abs(By)))
+print("Bz max abs =", np.max(np.abs(Bz)))
+
+print("\nComparison with intended 1e-4 sin(2z):")
+print("max absolute error =", max_abs_error)
+print("relative L2 error =", relative_l2)
+
+print("\nBest fit:")
+print("A_fit =", A_fit)
+print("target amplitude =", 1e-4)
+print(
+    "relative amplitude error =",
+    abs(A_fit - 1e-4) / 1e-4,
+)
+print("relative residual =", relative_residual)
+
+plt.figure(figsize=(10, 6))
+
+plt.plot(z, Bx, label=r"Struphy: $B_x(z)$")
+plt.plot(
+    z,
+    target,
+    "--",
+    label=r"Expected: $10^{-4}\sin(2z)$",
+)
+
+plt.xlabel("z")
+plt.ylabel(r"$B_x$")
+plt.title("Initial Magnetic Perturbation Validation")
+plt.legend()
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "magnetic_perturbation_validation.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("\nSaved:", path)
diff --git a/analysis/validate_v3_maxwellian.py b/analysis/validate_v3_maxwellian.py
new file mode 100644
index 000000000..08c6e38ff
--- /dev/null
+++ b/analysis/validate_v3_maxwellian.py
@@ -0,0 +1,115 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+BASE = (
+    "thesis_postprocess_validation/post_processing/"
+    "kinetic_data/hot_elec/distribution_function/v3_density"
+)
+
+TIME_PATH = (
+    "thesis_postprocess_validation/post_processing/t_grid.npy"
+)
+
+OUT = "validation_plots"
+os.makedirs(OUT, exist_ok=True)
+
+t = np.load(TIME_PATH)
+
+v3 = np.load(os.path.join(BASE, "grid_v3.npy"))
+F = np.load(os.path.join(BASE, "f_binned.npy"))
+DF = np.load(os.path.join(BASE, "delta_f_binned.npy"))
+
+dv = np.mean(np.diff(v3))
+
+density = 0.06
+sigma_target = 0.2
+
+analytic = (
+    density
+    / (np.sqrt(2 * np.pi) * sigma_target)
+    * np.exp(
+        -(v3 ** 2)
+        / (2 * sigma_target ** 2)
+    )
+)
+
+f0 = F[0]
+
+integral = np.sum(f0) * dv
+
+mean = np.sum(v3 * f0) * dv / integral
+
+variance = (
+    np.sum((v3 - mean) ** 2 * f0)
+    * dv
+    / integral
+)
+
+sigma = np.sqrt(variance)
+
+relative_l2 = (
+    np.linalg.norm(f0 - analytic)
+    / np.linalg.norm(analytic)
+)
+
+integrals = np.sum(F, axis=1) * dv
+
+print("============================================================")
+print("V3 MAXWELLIAN VALIDATION")
+print("============================================================")
+
+print("grid shape =", v3.shape)
+print("distribution shape =", F.shape)
+
+print("\nInitial integral =", integral)
+print("mean =", mean)
+print("sigma =", sigma)
+print("target sigma =", sigma_target)
+
+print("\nrelative L2 error =", relative_l2)
+
+print(
+    "maximum relative integral variation =",
+    np.max(np.abs(integrals - integrals[0]))
+    / abs(integrals[0]),
+)
+
+print(
+    "Delta-f identically zero =",
+    np.all(DF == 0),
+)
+
+plt.figure(figsize=(10, 6))
+
+plt.step(
+    v3,
+    f0 / density,
+    where="mid",
+    label=r"Loaded $v_3$ distribution",
+)
+
+plt.plot(
+    v3,
+    analytic / density,
+    "--",
+    label=r"Prescribed Maxwellian, $\sigma=0.2$",
+)
+
+plt.xlabel(r"$v_3$")
+plt.ylabel("Probability density")
+plt.title("Validation of Hot-Electron v3 Distribution")
+plt.legend()
+plt.grid()
+
+plt.tight_layout()
+
+path = os.path.join(
+    OUT,
+    "v3_distribution_validation.png",
+)
+
+plt.savefig(path, dpi=200)
+plt.close()
+
+print("\nSaved:", path)
diff --git a/examples/ColdPlasmaVlasov/anisotropic_maxwellian/params.py b/examples/ColdPlasmaVlasov/anisotropic_maxwellian/params.py
new file mode 100644
index 000000000..e4f51f41d
--- /dev/null
+++ b/examples/ColdPlasmaVlasov/anisotropic_maxwellian/params.py
@@ -0,0 +1,197 @@
+# -----------------------------
+# Description of the simulation
+# -----------------------------
+# Please fill in a verbal description of the simulation. 
+# It will be printed at the beginning of the simulation and can be used to keep track of the different runs.
+
+import logging
+import numpy as np
+
+from struphy import (
+    BaseUnits,
+    DerhamOptions,
+    EnvironmentOptions,
+    FieldsBackground,
+    Simulation,
+    Time,
+    domains,
+    equils,
+    grids,
+    perturbations,
+    LoadingParameters,
+    WeightsParameters,
+    BoundaryParameters,
+    SortingParameters,
+    SavingParameters,
+    BinningPlot,
+    maxwellians,
+    set_logging_level,
+)
+
+from struphy.models import ColdPlasmaVlasov
+
+set_logging_level(logging.INFO)
+
+name = "Thesis Figure 4.9 validation"
+
+description = """
+Short validation run for the geometric ColdPlasmaVlasov
+anisotropy-driven instability benchmark from the thesis.
+No control variate. Lie-Trotter splitting.
+"""
+
+# ----------------
+# Model
+# ----------------
+
+model = ColdPlasmaVlasov(
+    base_units=BaseUnits(),
+    thermal_alpha=-2.0,
+    thermal_epsilon=-1.0,
+    hot_epsilon=-1.0,
+)
+
+model.em_fields.e_field.save_data = True
+model.em_fields.b_field.save_data = True
+model.em_fields.phi.save_data = True
+model.thermal_elec.current.save_data = True
+model.hot_elec.var.save_data = True
+
+# ----------------
+# Simulation
+# ----------------
+
+env = EnvironmentOptions(
+    sim_folder="thesis_fig4_9_validation_short",
+    save_step=1,
+)
+
+time_opts = Time(
+    dt=0.0125,
+    Tend=1.0,
+    split_algo="LieTrotter",
+)
+
+# k = 2, hence Lz = 2*pi/k = pi.
+# The other two directions are inactive/minimal.
+
+domain = domains.Cuboid(
+    l1=0.0,
+    r1=1.0,
+    l2=0.0,
+    r2=1.0,
+    l3=0.0,
+    r3=np.pi,
+)
+
+equil = equils.HomogenSlab()
+
+grid = grids.TensorProductGrid(
+    num_elements=(1, 1, 32),
+)
+
+derham_opts = DerhamOptions(
+    degree=(1, 1, 1),
+    bcs=(None, None, None),
+)
+
+sim = Simulation(
+    model=model,
+    name=name,
+    description=description,
+    params_path=__file__,
+    env=env,
+    time_opts=time_opts,
+    domain=domain,
+    equil=equil,
+    grid=grid,
+    derham_opts=derham_opts,
+)
+
+# ----------------
+# Particles
+# ----------------
+
+loading_params = LoadingParameters(
+    Np=100000,
+    loading="pseudo_random",
+    seed=1234,
+    moments=(0.0, 0.0, 0.0, 0.53, 0.53, 0.20),
+)
+
+weights_params = WeightsParameters(
+    control_variate=False,
+)
+
+boundary_params = BoundaryParameters()
+sorting_params = SortingParameters()
+
+binplot = BinningPlot(
+    slice="e3",
+    n_bins=128,
+    ranges=(0.0, 1.0),
+)
+
+saving_params = SavingParameters(
+    binning_plots=(binplot,),
+)
+
+model.hot_elec.set_markers(
+    loading_params=loading_params,
+    weights_params=weights_params,
+    boundary_params=boundary_params,
+    sorting_params=sorting_params,
+    saving_params=saving_params,
+)
+
+# ----------------
+# Propagators
+# ----------------
+
+model.propagators.maxwell.options = model.propagators.maxwell.Options()
+model.propagators.ohm.options = model.propagators.ohm.Options()
+model.propagators.jxb.options = model.propagators.jxb.Options()
+model.propagators.push_eta.options = model.propagators.push_eta.Options()
+model.propagators.push_vxb.options = model.propagators.push_vxb.Options()
+model.propagators.coupling_va.options = model.propagators.coupling_va.Options()
+model.initial_poisson.options = model.initial_poisson.Options()
+
+# ----------------
+# Initial conditions
+# ----------------
+
+# Cold current initially zero.
+model.thermal_elec.current.add_background(
+    FieldsBackground(values=(0.0, 0.0, 0.0))
+)
+
+# Hot anisotropic Maxwellian:
+# parallel direction is v3 (along B0 / z)
+# perpendicular directions are v1 and v2.
+
+hot_background = maxwellians.Maxwellian3D(
+    n=(0.06, None),
+    vth1=(0.53, None),
+    vth2=(0.53, None),
+    vth3=(0.20, None),
+)
+
+model.hot_elec.var.add_background(hot_background)
+
+model.hot_elec.var.add_initial_condition(hot_background)
+# Magnetic perturbation:
+# B_x = 1e-4 sin(2 z)
+
+magnetic_perturbation = perturbations.ModesSin(
+    ls=(0,),
+    ms=(0,),
+    ns=(1,),
+    amps=(1e-4,),
+    given_in_basis="v",
+    comp=0,
+)
+
+model.em_fields.b_field.add_perturbation(magnetic_perturbation)
+
+if __name__ == "__main__":
+    sim.run()
diff --git a/examples/ColdPlasmaVlasov/anisotropic_maxwellian/pproc.py b/examples/ColdPlasmaVlasov/anisotropic_maxwellian/pproc.py
new file mode 100644
index 000000000..3ccf0e924
--- /dev/null
+++ b/examples/ColdPlasmaVlasov/anisotropic_maxwellian/pproc.py
@@ -0,0 +1,12 @@
+from struphy import PostProcessor
+
+def main():
+
+    pp = PostProcessor(
+        path_out="/home/anishojha/Research/struphy/thesis_postprocess_validation"
+    )
+
+    pp.process()
+
+if __name__ == "__main__":
+    main()
diff --git a/src/struphy/models/cold_plasma_vlasov.py b/src/struphy/models/cold_plasma_vlasov.py
index 1c134fa45..d9482fd35 100644
--- a/src/struphy/models/cold_plasma_vlasov.py
+++ b/src/struphy/models/cold_plasma_vlasov.py
@@ -219,6 +219,9 @@ def allocate_helpers(self):
         Propagator.derham.grad.dot(-phi, out=self.em_fields.e_field.spline.vector)
         logger.info("... Done.")
 
+        # reset particle weights
+        particles.weights = particles.weights_at_t0.copy()
+
     ## default parameters
     def generate_default_parameter_file(self, path=None, prompt=True):
         params_path = super().generate_default_parameter_file(path=path, prompt=prompt)