diff --git a/adept/_tf1d/modules.py b/adept/_tf1d/modules.py
index bdd53b4e..9816785d 100644
--- a/adept/_tf1d/modules.py
+++ b/adept/_tf1d/modules.py
@@ -59,7 +59,7 @@ def write_units(self):
         wp0 = np.sqrt(n0 * self.ureg.e**2.0 / (self.ureg.m_e * self.ureg.epsilon_0)).to("rad/s")
         tp0 = (1 / wp0).to("fs")
 
-        v0 = np.sqrt(2.0 * T0 / self.ureg.m_e).to("m/s")
+        v0 = np.sqrt(T0 / self.ureg.m_e).to("m/s")
         x0 = (v0 / wp0).to("nm")
         c_light = _Q(1.0 * self.ureg.c).to("m/s") / v0
         beta = (v0 / self.ureg.c).to("dimensionless")
@@ -72,7 +72,7 @@ def write_units(self):
         sim_duration = (self.cfg["grid"]["tmax"] * tp0).to("ps")
 
         # collisions
-        logLambda_ee = 23.5 - np.log(n0.magnitude**0.5 / T0.magnitude**-1.25)
+        logLambda_ee = 23.5 - np.log(n0.magnitude**0.5 * T0.magnitude**-1.25)
         logLambda_ee -= (1e-5 + (np.log(T0.magnitude) - 2) ** 2.0 / 16) ** 0.5
         nuee = _Q(2.91e-6 * n0.magnitude * logLambda_ee / T0.magnitude**1.5, "Hz")
         nuee_norm = nuee / wp0
diff --git a/adept/_vlasov1d/datamodel.py b/adept/_vlasov1d/datamodel.py
index 95f28d1b..13010629 100644
--- a/adept/_vlasov1d/datamodel.py
+++ b/adept/_vlasov1d/datamodel.py
@@ -71,7 +71,6 @@ class GridConfig(BaseModel):
     vmax: float | None = None
     vmin: float | None = None
     parallel: tuple[str, ...] | bool = False
-    c_light: float | None = None
 
     @field_validator("parallel", mode="before")
     @classmethod
diff --git a/adept/_vlasov1d/grid.py b/adept/_vlasov1d/grid.py
index 91b6d09b..b16db67e 100644
--- a/adept/_vlasov1d/grid.py
+++ b/adept/_vlasov1d/grid.py
@@ -52,7 +52,7 @@ def __init__(
         self.tmin = tmin
 
         # Compute dx
-        self.dx = xmax / nx
+        self.dx = (xmax - xmin) / nx
 
         # Override dt for EM wave stability if needed
         if should_override_dt_for_em_waves:
diff --git a/adept/_vlasov1d/helpers.py b/adept/_vlasov1d/helpers.py
index d34e774f..534047ac 100644
--- a/adept/_vlasov1d/helpers.py
+++ b/adept/_vlasov1d/helpers.py
@@ -72,7 +72,10 @@ def _initialize_supergaussian_distribution_(
     # Thermal velocity: v_t = sqrt(T/m)
     v_thermal = np.sqrt(T0 / mass)
 
-    # Alpha factor for supergaussian normalization
+    # Alpha factor for supergaussian normalization. This fixes the moment ratio
+    # <v^4>/<v^2> = 3*T0/mass for every order m; the realized *variance* equals
+    # T0/mass only at m=2 (Maxwellian). For m>2 flat-tops the second-moment
+    # temperature exceeds T0 (x1.24 at m=3, x1.37 at m=4). See docs config.md.
     alpha = np.sqrt(3.0 * gamma_3_over_m(supergaussian_order) / gamma_5_over_m(supergaussian_order))
 
     single_dist = -(np.power(np.abs((vax[None, :] - v0) / (alpha * v_thermal)), supergaussian_order))
diff --git a/adept/_vlasov1d/modules.py b/adept/_vlasov1d/modules.py
index 484c7ecd..c4840e77 100644
--- a/adept/_vlasov1d/modules.py
+++ b/adept/_vlasov1d/modules.py
@@ -118,6 +118,9 @@ def write_units(self) -> dict:
         sim_duration = (grid.tmax * norm.tau).to("ps")
 
         nu_ee = norm.approximate_ee_collision_frequency()
+        # e-e collision rate in code units (1/wp0). This is what the FP/Krook
+        # `baseline` rates should be compared against.
+        nu_ee_norm = (nu_ee * norm.tau).to("").magnitude
 
         beta = 1.0 / norm.speed_of_light_norm()
 
@@ -131,6 +134,7 @@ def write_units(self) -> dict:
             "beta": beta,
             "x0": norm.L0.to("nm"),
             "nuee": nu_ee.to("Hz"),
+            "nuee_norm": nu_ee_norm,
             "logLambda_ee": norm.logLambda_ee(),
             "box_length": box_length,
             "box_width": box_width,
@@ -238,10 +242,12 @@ def get_solver_quantities(self) -> dict:
             cfg_grid["species_grids"][species_name]["one_over_kvr"] = jnp.array(one_over_kvr)
 
             # Build species parameters (charge, mass, charge-to-mass ratio)
+            # T0 is the bulk (first-listed component) temperature, used e.g. by the Krook target
             cfg_grid["species_params"][species_name] = {
                 "charge": species_cfg.charge,
                 "mass": species_cfg.mass,
                 "charge_to_mass": species_cfg.charge / species_cfg.mass,
+                "T0": self.simulation.species_distributions[species_name][0].T0,
             }
 
         cfg_grid["n_prof_total"] = n_prof_total
diff --git a/adept/_vlasov1d/solvers/pushers/field.py b/adept/_vlasov1d/solvers/pushers/field.py
index 76a999af..ff73262c 100644
--- a/adept/_vlasov1d/solvers/pushers/field.py
+++ b/adept/_vlasov1d/solvers/pushers/field.py
@@ -227,7 +227,7 @@ def __call__(self, f_dict: dict, prev_ex: jnp.ndarray, dt: jnp.float64):
 class AmpereSolver:
     """Ampere solver using current density to evolve electric field.
 
-    Solves: a single Euler step of ∂E/∂t = -j, where j = Σ_s (q_s/m_s) ∫v f_s dv
+    Solves: a single Euler step of ∂E/∂t = -j, where j = Σ_s q_s ∫v f_s dv
     is the total current density from all species.
     """
 
diff --git a/adept/_vlasov1d/solvers/pushers/fokker_planck.py b/adept/_vlasov1d/solvers/pushers/fokker_planck.py
index 23166bfd..dfb2a5b1 100644
--- a/adept/_vlasov1d/solvers/pushers/fokker_planck.py
+++ b/adept/_vlasov1d/solvers/pushers/fokker_planck.py
@@ -76,9 +76,9 @@ def compute_vbar(self, f: Array) -> Array:
             f: Distribution function, shape (..., nv)
 
         Returns:
-            vbar: Mean velocity <v>, shape (...)
+            vbar: Mean velocity <v> = ∫v f dv / ∫f dv, shape (...)
         """
-        return jnp.sum(f * self.v, axis=-1) * self.dv
+        return jnp.sum(f * self.v, axis=-1) / jnp.sum(f, axis=-1)
 
     def compute_D(self, f: Array, beta: Array) -> Array:
         """
@@ -242,7 +242,11 @@ def __init__(self, cfg: Mapping[str, Any]):
         self.cfg = cfg
         v = cfg["grid"]["species_grids"]["electron"]["v"]
         dv = cfg["grid"]["species_grids"]["electron"]["dv"]
-        f_mx = np.exp(-(v[None, :] ** 2.0) / 2.0)
+        params = cfg["grid"].get("species_params", {}).get("electron", {})
+        T0 = params.get("T0", 1.0)
+        mass = params.get("mass", 1.0)
+        # Maxwellian with variance T0/m (the species' bulk thermal width)
+        f_mx = np.exp(-(v[None, :] ** 2.0) / (2.0 * T0 / mass))
         self.f_mx = f_mx / np.sum(f_mx, axis=1)[:, None] / dv
         self.dv = dv
 
diff --git a/adept/_vlasov1d/solvers/pushers/vlasov.py b/adept/_vlasov1d/solvers/pushers/vlasov.py
index 384d6f43..828dbf7a 100644
--- a/adept/_vlasov1d/solvers/pushers/vlasov.py
+++ b/adept/_vlasov1d/solvers/pushers/vlasov.py
@@ -28,7 +28,7 @@ def __init__(self, cfg, interp_e):
         """Create interpolation helpers for x-advection and v-advection steps."""
         self.x = cfg["grid"]["x"]
         self.v = cfg["grid"]["v"]
-        self.f_interp = partial(interp2d, x=self.x, y=self.v, period=cfg["grid"]["xmax"])
+        self.f_interp = partial(interp2d, x=self.x, y=self.v, period=cfg["grid"]["xmax"] - cfg["grid"]["xmin"])
         self.dt = cfg["grid"]["dt"]
         self.dummy_x = jnp.ones_like(self.x)
         self.dummy_v = jnp.ones_like(self.v)[None, :]
diff --git a/adept/_vlasov1d/storage.py b/adept/_vlasov1d/storage.py
index 1075ed48..6899e049 100644
--- a/adept/_vlasov1d/storage.py
+++ b/adept/_vlasov1d/storage.py
@@ -137,11 +137,12 @@ def _calc_moment_(inp, _dv=dv):
                 f = y[species_name]
                 species_moments = {}
                 species_moments["n"] = _calc_moment_(f)
-                species_moments["v"] = _calc_moment_(f * v[None, :])
+                species_moments["j"] = _calc_moment_(f * v[None, :])
+                species_moments["v"] = species_moments["j"] / species_moments["n"]
                 v_m_vbar = v[None, :] - species_moments["v"][:, None]
                 species_moments["p"] = _calc_moment_(f * v_m_vbar**2.0)
                 species_moments["q"] = _calc_moment_(f * v_m_vbar**3.0)
-                species_moments["-flogf"] = _calc_moment_(f * jnp.log(jnp.abs(f)))
+                species_moments["-flogf"] = _calc_moment_(-jnp.abs(f) * jnp.log(jnp.abs(f)))
                 species_moments["f^2"] = _calc_moment_(f * f)
 
                 result[species_name] = species_moments
@@ -292,9 +293,11 @@ def save(t, y, args):
         scalars = {}
 
         # Compute scalars for each species
+        mean_kinetic_energy = 0.0
         for species_name in species_names:
             v = species_grids[species_name]["v"][None, :]
             dv = species_grids[species_name]["dv"]
+            mass = cfg["grid"]["species_params"][species_name]["mass"]
 
             def _calc_mean_moment_(inp, _dv=dv):
                 return jnp.mean(jnp.sum(inp, axis=1) * _dv)
@@ -306,12 +309,19 @@ def _calc_mean_moment_(inp, _dv=dv):
             scalars[f"mean_q_{species_name}"] = _calc_mean_moment_(f * v**3.0)
             scalars[f"mean_-flogf_{species_name}"] = _calc_mean_moment_(-jnp.log(jnp.abs(f)) * jnp.abs(f))
             scalars[f"mean_f2_{species_name}"] = _calc_mean_moment_(f * f)
+            mean_kinetic_energy += 0.5 * mass * scalars[f"mean_P_{species_name}"]
 
         # Shared field scalars (not species-specific)
         scalars["mean_de2"] = jnp.mean(y["de"] ** 2.0)
         scalars["mean_e2"] = jnp.mean(y["e"] ** 2.0)
         scalars["mean_pond"] = jnp.mean(-0.5 * jnp.gradient(y["a"] ** 2.0, cfg["grid"]["dx"])[1:-1])
 
+        # Energy conservation monitor (electrostatic): x-averaged kinetic + E-field
+        # energy density. Transverse (EM) field energy is not included.
+        scalars["mean_kinetic_energy"] = mean_kinetic_energy
+        scalars["mean_field_energy"] = 0.5 * scalars["mean_e2"]
+        scalars["mean_total_energy"] = mean_kinetic_energy + 0.5 * scalars["mean_e2"]
+
         return scalars
 
     return save
diff --git a/adept/normalization.py b/adept/normalization.py
index aa9f7714..ae23bbc2 100644
--- a/adept/normalization.py
+++ b/adept/normalization.py
@@ -34,7 +34,7 @@ class PlasmaNormalization:
     def logLambda_ee(self) -> float:
         n0_cc = self.n0.to("1/cc").magnitude
         T0_eV = self.T0.to("eV").magnitude
-        logLambda_ee = 23.5 - jnp.log(n0_cc**0.5 / T0_eV**-1.25)
+        logLambda_ee = 23.5 - jnp.log(n0_cc**0.5 * T0_eV**-1.25)
         logLambda_ee -= (1e-5 + (jnp.log(T0_eV) - 2) ** 2.0 / 16) ** 0.5
         return logLambda_ee
 
@@ -78,9 +78,13 @@ def electron_debye_normalization(n0_str, T0_str):
     """
     Returns the electron thermal normalization for the given density and temperature.
     Unit quantities are:
-        - Debye length
-        - Electron thermal velocity
-        - Langmuir oscillation frequency
+        - L0 = Debye length, sqrt(eps0 T0 / (n0 e^2))
+        - v0 = electron thermal velocity sqrt(T0/m_e) (RMS / standard-deviation
+          convention: a T=1 Maxwellian is exp(-v^2/2) with unit variance)
+        - tau = 1/wp0 (inverse Langmuir oscillation frequency)
+
+    Under this convention a code wavenumber is k*lambda_De and the Bohm-Gross
+    dispersion reads w^2 = 1 + 3 k^2.
     """
     n0 = UREG.Quantity(n0_str)
     T0 = UREG.Quantity(T0_str)
@@ -88,7 +92,7 @@ def electron_debye_normalization(n0_str, T0_str):
     wp0 = ((n0 * UREG.e**2.0 / (UREG.m_e * UREG.epsilon_0)) ** 0.5).to("rad/s")
     tau = 1 / wp0
 
-    v0 = ((2.0 * T0 / UREG.m_e) ** 0.5).to("m/s")
+    v0 = ((T0 / UREG.m_e) ** 0.5).to("m/s")
     x0 = (v0 / wp0).to("nm")
 
     return PlasmaNormalization(m0=UREG.m_e, q0=UREG.e, n0=n0, T0=T0, L0=x0, v0=v0, tau=tau)
diff --git a/configs/vlasov-1d/srs-debug-small.yaml b/configs/vlasov-1d/srs-debug-small.yaml
index c0bfa5ec..5ac73e31 100644
--- a/configs/vlasov-1d/srs-debug-small.yaml
+++ b/configs/vlasov-1d/srs-debug-small.yaml
@@ -48,13 +48,18 @@ mlflow:
   experiment: vlasov1d-srs
   run: debug-small
 
+# The ey driver pair is a backscatter-SRS matched triad in code units
+# (v0 = sqrt(T0/m_e), k in 1/lambda_De, w in wp0, c_hat = c/v0 = 15.984 at 2000 eV):
+#   pump:  w0 = 2.79,    k0 = +sqrt(w0^2 - 1)/c_hat = +0.162949 (rightgoing)
+#   seed:  w0 = 1.700942, k0 = -0.086080 (leftgoing)
+#   EPW:   k lambda_De = 0.249029, w = 1.089058 = w_pump - w_seed
 drivers:
   ex: {}
   ey:
     '0':
       params:
         a0: 1.0
-        k0: 1.0
+        k0: 0.162949
         w0: 2.79
         dw0: 0.
       envelope:
@@ -69,8 +74,8 @@ drivers:
     '1':
       params:
         a0: 1.e-6
-        k0: 1.0
-        w0: 1.63
+        k0: -0.086080
+        w0: 1.700942
         dw0: 0.
       envelope:
         time:
diff --git a/configs/vlasov-1d/wavepacket.yaml b/configs/vlasov-1d/wavepacket.yaml
index 3770512f..5ceba970 100644
--- a/configs/vlasov-1d/wavepacket.yaml
+++ b/configs/vlasov-1d/wavepacket.yaml
@@ -35,7 +35,6 @@ drivers:
           width: 900.0
   ey: {}
 grid:
-  c_light: 10
   dt: 0.1
   nv: 6144
   nx: 4096
diff --git a/docs/source/solvers/vlasov1d/config.md b/docs/source/solvers/vlasov1d/config.md
index 6263565d..a7e38f63 100644
--- a/docs/source/solvers/vlasov1d/config.md
+++ b/docs/source/solvers/vlasov1d/config.md
@@ -48,6 +48,28 @@ units:
   normalizing_density: 1.5e21/cc
 ```
 
+### Normalization convention
+
+All Vlasov-1D quantities are normalized using the following unit system, built from
+`normalizing_density` ($n_0$) and `normalizing_temperature` ($T_0$):
+
+| Unit | Definition | Meaning |
+|------|-----------|---------|
+| time | $\tau = 1/\omega_{p0}$, $\omega_{p0} = \sqrt{n_0 e^2/(\epsilon_0 m_e)}$ | inverse plasma frequency |
+| velocity | $v_0 = \sqrt{T_0/m_e}$ | electron thermal speed (RMS / standard-deviation convention) |
+| length | $L_0 = v_0/\omega_{p0} = \lambda_{De}$ | electron Debye length |
+| wavenumber | $1/L_0$ | a code wavenumber is $k \lambda_{De}$ |
+
+Consequences of the $v_0 = \sqrt{T_0/m_e}$ (σ) convention:
+
+- A Maxwellian at code temperature $T = 1$ is $f \propto e^{-v^2/2}$ (unit variance).
+- The Bohm–Gross dispersion in code units is $\omega^2 = 1 + 3 k^2$.
+- The normalized speed of light is $\hat c = c/\sqrt{T_0/m_e}$ (e.g. $\hat c = 15.98$ at 2000 eV).
+
+Dimensional inputs (strings with units, e.g. `xmax: 100um`) are converted with these
+units; plain numeric inputs are taken to already be in code units and pass through
+unchanged.
+
 ## density
 
 Species and density configuration. You can define multiple "species" by using keys prefixed with `species-`.
@@ -65,9 +87,9 @@ Each species is defined with a key starting with `species-` (e.g., `species-back
 | `noise_seed` | int | Random seed for noise initialization |
 | `noise_type` | string | `"gaussian"` or `"uniform"` |
 | `noise_val` | float | Amplitude of noise |
-| `v0` | float | Drift velocity (normalized) |
-| `T0` | float | Temperature (normalized to `normalizing_temperature`) |
-| `m` | float | Exponent for super-Gaussian distribution. `2.0` is Maxwellian |
+| `v0` | float | Drift velocity in code units of $\sqrt{T_0/m_e}$ (thermal-σ units). Numeric only — dimensional strings are not supported here |
+| `T0` | float | Temperature in units of `normalizing_temperature`. The initialized distribution has velocity variance `T0/mass`. Numeric only |
+| `m` | float | Exponent for super-Gaussian distribution $f \propto \exp[-\|v/(\alpha v_{th})\|^m]$. `2.0` is Maxwellian. Note: $\alpha$ is chosen to fix the moment ratio $\langle v^4\rangle/\langle v^2\rangle = 3\,T_0/\mathrm{mass}$ for all $m$; the *variance* equals `T0/mass` only at `m: 2`. For flat-top distributions (`m > 2`) the second-moment temperature diagnostic will read higher than `T0` (e.g. ×1.24 at `m: 3`, ×1.37 at `m: 4`) |
 | `basis` | string | Spatial profile type (see below) |
 
 #### Basis Types
@@ -535,7 +557,7 @@ Both `fokker_planck` and `krook` use profile objects to define spatiotemporal va
 
 | Field | Type | Description |
 |-------|------|-------------|
-| `baseline` | float | Base collision frequency |
+| `baseline` | float | Base collision frequency in code units of $\omega_{p0}$ (i.e. $\hat\nu = \nu[\mathrm{s}^{-1}]/\omega_{p0}$, a true rate — no $2\pi$). For reference, the NRL e–e rate in these units is logged as `nuee_norm` in `units.yaml`. Note the Dougherty/Lenard–Bernstein $\nu$ agrees with the NRL $\nu_{ee}$ only up to an O(1) factor |
 | `bump_or_trough` | string | `"bump"` or `"trough"` |
 | `center` | float | Profile center |
 | `rise` | float | Transition steepness |
diff --git a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml
index 4134829c..9539860e 100644
--- a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml
+++ b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_array_config.yml
@@ -38,7 +38,7 @@ drivers:
         w0: 1.0
   ey: {}
 grid:
-  beta: 0.088474881879774
+  beta: 0.062561188941867
   dt: 0.1
   dx: 0.785
   ion_charge:
@@ -16507,6 +16507,7 @@ grid:
       vmax: 6.4
   species_params:
     electron:
+      T0: 1.0
       charge: -1.0
       charge_to_mass: -1.0
       mass: 1.0
@@ -17723,17 +17724,18 @@ terms:
 units:
   derived:
     T0: 2000 electron_volt
-    beta: 0.088474881879774
-    box_length: 0.07623654008768617 micron
+    beta: 0.062561188941867
+    box_length: 0.05390737447020297 micron
     box_width: inf
-    c_light: 11.302642950785
-    logLambda_ee: -11.781233290653
+    c_light: 15.984350951661
+    logLambda_ee: 7.221022858202
     n0: 1.5e+21 / cubic_centimeter
-    nuee: -574949910190.1326 hertz
+    nuee: 352401683370.27344 hertz
+    nuee_norm: 0.00016128753860756
     sim_duration: 0.004622577634581562 picosecond
     tp0: 0.45768095391896646 femtosecond
-    v0: 26524102.30999721 meter / second
+    v0: 18755372.608284798 meter / second
     wp0: 2184928150138955.5 radian / second
-    x0: 12.139576447083783 nanometer
+    x0: 8.58397682646544 nanometer
   normalizing_density: 1.5e21/cc
   normalizing_temperature: 2000eV
diff --git a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml
index ebca6209..f83b3fe6 100644
--- a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml
+++ b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_derived_config.yml
@@ -38,7 +38,7 @@ drivers:
         w0: 1.0
   ey: {}
 grid:
-  beta: 0.088474881879774
+  beta: 0.062561188941867
   dt: 0.1
   dx: 0.785
   max_steps: 105
@@ -111,17 +111,18 @@ terms:
 units:
   derived:
     T0: 2000 electron_volt
-    beta: 0.088474881879774
-    box_length: 0.07623654008768617 micron
+    beta: 0.062561188941867
+    box_length: 0.05390737447020297 micron
     box_width: inf
-    c_light: 11.302642950785
-    logLambda_ee: -11.781233290653
+    c_light: 15.984350951661
+    logLambda_ee: 7.221022858202
     n0: 1.5e+21 / cubic_centimeter
-    nuee: -574949910190.1326 hertz
+    nuee: 352401683370.27344 hertz
+    nuee_norm: 0.00016128753860756
     sim_duration: 0.004622577634581562 picosecond
     tp0: 0.45768095391896646 femtosecond
-    v0: 26524102.30999721 meter / second
+    v0: 18755372.608284798 meter / second
     wp0: 2184928150138955.5 radian / second
-    x0: 12.139576447083783 nanometer
+    x0: 8.58397682646544 nanometer
   normalizing_density: 1.5e21/cc
   normalizing_temperature: 2000eV
diff --git a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml
index dca5edde..f1e7ae4b 100644
--- a/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml
+++ b/tests/test_vlasov1d/test_config_regression/fokker_planck_conservation_units.yml
@@ -1,13 +1,14 @@
 T0: 2000 electron_volt
-beta: 0.088474881879774
-box_length: 0.07623654008768617 micron
+beta: 0.062561188941867
+box_length: 0.05390737447020297 micron
 box_width: inf
-c_light: 11.302642950785
-logLambda_ee: -11.781233290653
+c_light: 15.984350951661
+logLambda_ee: 7.221022858202
 n0: 1.5e+21 / cubic_centimeter
-nuee: -574949910190.1326 hertz
+nuee: 352401683370.27344 hertz
+nuee_norm: 0.00016128753860756
 sim_duration: 0.004622577634581562 picosecond
 tp0: 0.45768095391896646 femtosecond
-v0: 26524102.30999721 meter / second
+v0: 18755372.608284798 meter / second
 wp0: 2184928150138955.5 radian / second
-x0: 12.139576447083783 nanometer
+x0: 8.58397682646544 nanometer
diff --git a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml
index 927d9760..e899d7fa 100644
--- a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml
+++ b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_array_config.yml
@@ -29,7 +29,7 @@ drivers:
   ex: {}
   ey: {}
 grid:
-  beta: 0.088474881879774
+  beta: 0.062561188941867
   dt: 0.5
   dx: 9.81746875
   ion_charge:
@@ -28711,10 +28711,12 @@ grid:
       vmax: 0.005
   species_params:
     electron:
+      T0: 1.0
       charge: -1.0
       charge_to_mass: -1.0
       mass: 1.0
     ion:
+      T0: 0.01
       charge: 10.0
       charge_to_mass: 0.00054466230936819
       mass: 18360.0
@@ -38873,17 +38875,18 @@ terms:
 units:
   derived:
     T0: 2000 electron_volt
-    beta: 0.088474881879774
-    box_length: 3.8137571970393944 micron
+    beta: 0.062561188941867
+    box_length: 2.696733575825556 micron
     box_width: inf
-    c_light: 11.302642950785
-    logLambda_ee: -11.781233290653
+    c_light: 15.984350951661
+    logLambda_ee: 7.221022858202
     n0: 1.5e+21 / cubic_centimeter
-    nuee: -574949910190.1326 hertz
+    nuee: 352401683370.27344 hertz
+    nuee_norm: 0.00016128753860756
     sim_duration: 2.288633610071792 picosecond
     tp0: 0.45768095391896646 femtosecond
-    v0: 26524102.30999721 meter / second
+    v0: 18755372.608284798 meter / second
     wp0: 2184928150138955.5 radian / second
-    x0: 12.139576447083783 nanometer
+    x0: 8.58397682646544 nanometer
   normalizing_density: 1.5e21/cc
   normalizing_temperature: 2000eV
diff --git a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml
index e7c61532..9fe52c26 100644
--- a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml
+++ b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_derived_config.yml
@@ -29,7 +29,7 @@ drivers:
   ex: {}
   ey: {}
 grid:
-  beta: 0.088474881879774
+  beta: 0.062561188941867
   dt: 0.5
   dx: 9.81746875
   max_steps: 10005
@@ -120,17 +120,18 @@ terms:
 units:
   derived:
     T0: 2000 electron_volt
-    beta: 0.088474881879774
-    box_length: 3.8137571970393944 micron
+    beta: 0.062561188941867
+    box_length: 2.696733575825556 micron
     box_width: inf
-    c_light: 11.302642950785
-    logLambda_ee: -11.781233290653
+    c_light: 15.984350951661
+    logLambda_ee: 7.221022858202
     n0: 1.5e+21 / cubic_centimeter
-    nuee: -574949910190.1326 hertz
+    nuee: 352401683370.27344 hertz
+    nuee_norm: 0.00016128753860756
     sim_duration: 2.288633610071792 picosecond
     tp0: 0.45768095391896646 femtosecond
-    v0: 26524102.30999721 meter / second
+    v0: 18755372.608284798 meter / second
     wp0: 2184928150138955.5 radian / second
-    x0: 12.139576447083783 nanometer
+    x0: 8.58397682646544 nanometer
   normalizing_density: 1.5e21/cc
   normalizing_temperature: 2000eV
diff --git a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml
index b05d7d96..e018602a 100644
--- a/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml
+++ b/tests/test_vlasov1d/test_config_regression/multispecies_ion_acoustic_units.yml
@@ -1,13 +1,14 @@
 T0: 2000 electron_volt
-beta: 0.088474881879774
-box_length: 3.8137571970393944 micron
+beta: 0.062561188941867
+box_length: 2.696733575825556 micron
 box_width: inf
-c_light: 11.302642950785
-logLambda_ee: -11.781233290653
+c_light: 15.984350951661
+logLambda_ee: 7.221022858202
 n0: 1.5e+21 / cubic_centimeter
-nuee: -574949910190.1326 hertz
+nuee: 352401683370.27344 hertz
+nuee_norm: 0.00016128753860756
 sim_duration: 2.288633610071792 picosecond
 tp0: 0.45768095391896646 femtosecond
-v0: 26524102.30999721 meter / second
+v0: 18755372.608284798 meter / second
 wp0: 2184928150138955.5 radian / second
-x0: 12.139576447083783 nanometer
+x0: 8.58397682646544 nanometer
diff --git a/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml b/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml
index 8912206a..0924eec1 100644
--- a/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml
+++ b/tests/test_vlasov1d/test_config_regression/resonance_array_config.yml
@@ -38,7 +38,7 @@ drivers:
         w0: 1.1598
   ey: {}
 grid:
-  beta: 0.088474881879774
+  beta: 0.062561188941867
   dt: 0.25
   dx: 0.654375
   ion_charge:
@@ -20747,6 +20747,7 @@ grid:
       vmax: 6.4
   species_params:
     electron:
+      T0: 1.0
       charge: -1.0
       charge_to_mass: -1.0
       mass: 1.0
@@ -23327,17 +23328,18 @@ terms:
 units:
   derived:
     T0: 2000 electron_volt
-    beta: 0.088474881879774
-    box_length: 0.25420273080193445 micron
+    beta: 0.062561188941867
+    box_length: 0.17974847474618633 micron
     box_width: inf
-    c_light: 11.302642950785
-    logLambda_ee: -11.781233290653
+    c_light: 15.984350951661
+    logLambda_ee: 7.221022858202
     n0: 1.5e+21 / cubic_centimeter
-    nuee: -574949910190.1326 hertz
+    nuee: 352401683370.27344 hertz
+    nuee_norm: 0.00016128753860756
     sim_duration: 0.21980127811958367 picosecond
     tp0: 0.45768095391896646 femtosecond
-    v0: 26524102.30999721 meter / second
+    v0: 18755372.608284798 meter / second
     wp0: 2184928150138955.5 radian / second
-    x0: 12.139576447083783 nanometer
+    x0: 8.58397682646544 nanometer
   normalizing_density: 1.5e21/cc
   normalizing_temperature: 2000eV
diff --git a/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml b/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml
index 66794372..6e5e2484 100644
--- a/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml
+++ b/tests/test_vlasov1d/test_config_regression/resonance_derived_config.yml
@@ -38,7 +38,7 @@ drivers:
         w0: 1.1598
   ey: {}
 grid:
-  beta: 0.088474881879774
+  beta: 0.062561188941867
   dt: 0.25
   dx: 0.654375
   max_steps: 1925
@@ -119,17 +119,18 @@ terms:
 units:
   derived:
     T0: 2000 electron_volt
-    beta: 0.088474881879774
-    box_length: 0.25420273080193445 micron
+    beta: 0.062561188941867
+    box_length: 0.17974847474618633 micron
     box_width: inf
-    c_light: 11.302642950785
-    logLambda_ee: -11.781233290653
+    c_light: 15.984350951661
+    logLambda_ee: 7.221022858202
     n0: 1.5e+21 / cubic_centimeter
-    nuee: -574949910190.1326 hertz
+    nuee: 352401683370.27344 hertz
+    nuee_norm: 0.00016128753860756
     sim_duration: 0.21980127811958367 picosecond
     tp0: 0.45768095391896646 femtosecond
-    v0: 26524102.30999721 meter / second
+    v0: 18755372.608284798 meter / second
     wp0: 2184928150138955.5 radian / second
-    x0: 12.139576447083783 nanometer
+    x0: 8.58397682646544 nanometer
   normalizing_density: 1.5e21/cc
   normalizing_temperature: 2000eV
diff --git a/tests/test_vlasov1d/test_config_regression/resonance_units.yml b/tests/test_vlasov1d/test_config_regression/resonance_units.yml
index cc18e176..22e9776c 100644
--- a/tests/test_vlasov1d/test_config_regression/resonance_units.yml
+++ b/tests/test_vlasov1d/test_config_regression/resonance_units.yml
@@ -1,13 +1,14 @@
 T0: 2000 electron_volt
-beta: 0.088474881879774
-box_length: 0.25420273080193445 micron
+beta: 0.062561188941867
+box_length: 0.17974847474618633 micron
 box_width: inf
-c_light: 11.302642950785
-logLambda_ee: -11.781233290653
+c_light: 15.984350951661
+logLambda_ee: 7.221022858202
 n0: 1.5e+21 / cubic_centimeter
-nuee: -574949910190.1326 hertz
+nuee: 352401683370.27344 hertz
+nuee_norm: 0.00016128753860756
 sim_duration: 0.21980127811958367 picosecond
 tp0: 0.45768095391896646 femtosecond
-v0: 26524102.30999721 meter / second
+v0: 18755372.608284798 meter / second
 wp0: 2184928150138955.5 radian / second
-x0: 12.139576447083783 nanometer
+x0: 8.58397682646544 nanometer
diff --git a/tests/test_vlasov1d/test_fp_momentum_conservation.py b/tests/test_vlasov1d/test_fp_momentum_conservation.py
new file mode 100644
index 00000000..a50d9f66
--- /dev/null
+++ b/tests/test_vlasov1d/test_fp_momentum_conservation.py
@@ -0,0 +1,83 @@
+#  Copyright (c) Ergodic LLC 2026
+#  research@ergodic.io
+"""
+Momentum conservation of the Dougherty operator where n(x) != 1.
+
+The Dougherty drag must center on the *mean velocity* u(x) = ∫v f dv / n(x).
+A drag centered on the raw first moment ∫v f dv = n*u (the historical bug)
+breaks momentum conservation exactly where the density deviates from 1, e.g.
+for bump-on-tail or large-amplitude density perturbations with collisions on.
+
+This test applies the collision operator in isolation (no advection, no
+fields) to a drifting Maxwellian with a strong density modulation and checks
+per-cell conservation of density, momentum, and energy.
+"""
+
+import jax
+
+jax.config.update("jax_enable_x64", True)
+
+import numpy as np
+import pytest
+from jax import numpy as jnp
+
+from adept._vlasov1d.solvers.pushers.fokker_planck import Collisions
+
+
+def _make_cfg(nx, nv, vmax, fp_type):
+    dv = 2.0 * vmax / nv
+    v = jnp.linspace(-vmax + dv / 2.0, vmax - dv / 2.0, nv)
+    return {
+        "grid": {
+            "species_grids": {"electron": {"v": v, "dv": dv, "nv": nv, "vmax": vmax}},
+            "species_params": {"electron": {"charge": -1.0, "mass": 1.0, "charge_to_mass": -1.0, "T0": 1.0}},
+        },
+        "terms": {
+            "fokker_planck": {"is_on": True, "type": fp_type},
+            "krook": {"is_on": False},
+        },
+    }
+
+
+# Energy tolerance is per-scheme: Chang-Cooper preserves the discrete equilibrium
+# to ~1e-6, plain central differencing only to ~1e-3 at this resolution.
+@pytest.mark.parametrize("fp_type,energy_rtol", [("dougherty", 5e-3), ("chang_cooper_dougherty", 1e-5)])
+def test_dougherty_momentum_conservation_with_density_perturbation(fp_type, energy_rtol):
+    nx, nv, vmax = 16, 512, 6.4
+    cfg = _make_cfg(nx, nv, vmax, fp_type)
+    v = np.array(cfg["grid"]["species_grids"]["electron"]["v"])
+    dv = cfg["grid"]["species_grids"]["electron"]["dv"]
+
+    # Strong density modulation and a finite drift: exactly the regime where
+    # centering the drag on n*u instead of u destroys momentum conservation.
+    x = np.linspace(0, 2 * np.pi, nx, endpoint=False)
+    n_prof = 1.0 + 0.5 * np.sin(x)
+    u0 = 0.5
+    f = n_prof[:, None] * np.exp(-((v[None, :] - u0) ** 2) / 2.0)
+    f = f / (np.sum(np.exp(-((v - u0) ** 2) / 2.0)) * dv)
+
+    f = jnp.array(f)
+    collisions = Collisions(cfg)
+
+    nu = jnp.ones(nx)
+    dt = 0.1
+    f_out = f
+    for _ in range(50):
+        f_out = collisions(nu, None, f_out, dt)
+
+    n0_ = np.sum(np.array(f), axis=1) * dv
+    n1 = np.sum(np.array(f_out), axis=1) * dv
+    p0 = np.sum(np.array(f) * v[None, :], axis=1) * dv
+    p1 = np.sum(np.array(f_out) * v[None, :], axis=1) * dv
+    e0 = np.sum(np.array(f) * v[None, :] ** 2, axis=1) * dv
+    e1 = np.sum(np.array(f_out) * v[None, :] ** 2, axis=1) * dv
+
+    # The historical n*u-centered drag gives O(50%) momentum errors here;
+    # the corrected operator conserves to discretization level (~1e-7).
+    np.testing.assert_allclose(n1, n0_, rtol=1e-10, err_msg="density not conserved")
+    np.testing.assert_allclose(p1, p0, rtol=1e-6, err_msg="momentum not conserved where n(x) != 1")
+    np.testing.assert_allclose(e1, e0, rtol=energy_rtol, err_msg="energy not conserved")
+
+    # And the distribution must still be centered on u0, not n(x)*u0
+    u_final = p1 / n1
+    np.testing.assert_allclose(u_final, u0, atol=1e-6)
diff --git a/tests/test_vlasov1d/test_units_boundary.py b/tests/test_vlasov1d/test_units_boundary.py
new file mode 100644
index 00000000..e7af1c38
--- /dev/null
+++ b/tests/test_vlasov1d/test_units_boundary.py
@@ -0,0 +1,74 @@
+#  Copyright (c) Ergodic LLC 2026
+#  research@ergodic.io
+"""
+Units-boundary regression tests for the Vlasov-1D normalization layer.
+
+Every physics test in this suite works in dimensionless code units, which
+validates the engine but not the dimensional dictionary: numeric inputs bypass
+``normalize()`` entirely. These tests cross the units boundary — dimensional
+input in, dimensional quantity out — and compare against *independent* physical
+references (CODATA constants via scipy, and the NRL formulary), so a convention
+error in ``normalization.py`` cannot hide.
+
+The engine convention under test: v0 = sqrt(T0/m_e) (RMS / standard-deviation
+thermal speed), L0 = v0/wp0 = lambda_De, Maxwellian exp(-v^2/2) at T=1.
+"""
+
+import numpy as np
+from scipy import constants as csts
+
+from adept.normalization import electron_debye_normalization, normalize
+
+N0_STR, T0_STR = "1.5e21/cc", "2000eV"
+N0_SI = 1.5e21 * 1e6  # 1/m^3
+T0_J = 2000.0 * csts.e
+
+
+def test_debye_normalization_against_codata():
+    """v0, L0, tau, and c_hat must match hand-computed CODATA values."""
+    norm = electron_debye_normalization(N0_STR, T0_STR)
+
+    v_th = np.sqrt(T0_J / csts.m_e)  # sigma convention: sqrt(T/m)
+    wp0 = np.sqrt(N0_SI * csts.e**2 / (csts.epsilon_0 * csts.m_e))
+    lambda_de = v_th / wp0
+
+    np.testing.assert_allclose(norm.v0.to("m/s").magnitude, v_th, rtol=1e-6)
+    np.testing.assert_allclose(norm.L0.to("m").magnitude, lambda_de, rtol=1e-6)
+    np.testing.assert_allclose(norm.tau.to("s").magnitude, 1.0 / wp0, rtol=1e-6)
+    np.testing.assert_allclose(norm.speed_of_light_norm(), csts.c / v_th, rtol=1e-6)
+
+
+def test_debye_length_against_nrl_formulary():
+    """lambda_De = 7.43e2 sqrt(T_eV/n_cc) cm (NRL formulary) — a truly external reference."""
+    norm = electron_debye_normalization(N0_STR, T0_STR)
+    lambda_de_nrl_m = 7.43e2 * np.sqrt(2000.0 / 1.5e21) * 1e-2
+    np.testing.assert_allclose(norm.L0.to("m").magnitude, lambda_de_nrl_m, rtol=1e-3)
+
+
+def test_dimensional_string_round_trip():
+    """String inputs must convert with L0 = lambda_De (not sqrt(2) lambda_De)."""
+    norm = electron_debye_normalization(N0_STR, T0_STR)
+
+    v_th = np.sqrt(T0_J / csts.m_e)
+    wp0 = np.sqrt(N0_SI * csts.e**2 / (csts.epsilon_0 * csts.m_e))
+    lambda_de = v_th / wp0
+
+    # x: 100 um -> 100e-6 / lambda_De code units
+    np.testing.assert_allclose(normalize("100um", norm, dim="x"), 100e-6 / lambda_de, rtol=1e-6)
+    # k: 1/um -> lambda_De / 1e-6 code units (k lambda_De)
+    np.testing.assert_allclose(normalize("1/um", norm, dim="k"), lambda_de / 1e-6, rtol=1e-6)
+    # t: 1 ps -> wp0 * 1e-12
+    np.testing.assert_allclose(normalize("1ps", norm, dim="t"), wp0 * 1e-12, rtol=1e-6)
+    # numeric inputs pass through untouched
+    assert normalize(3.25, norm, dim="x") == 3.25
+
+
+def test_logged_collision_frequency_is_physical():
+    """logLambda and nuee must be positive and match the NRL expression."""
+    norm = electron_debye_normalization(N0_STR, T0_STR)
+    log_lambda = float(norm.logLambda_ee())
+    # NRL: 23.5 - ln(n^1/2 T^-5/4) - [1e-5 + (ln T - 2)^2/16]^1/2
+    expected = 23.5 - np.log(np.sqrt(1.5e21) * 2000.0**-1.25) - np.sqrt(1e-5 + (np.log(2000.0) - 2.0) ** 2 / 16.0)
+    np.testing.assert_allclose(log_lambda, expected, rtol=1e-10)
+    assert log_lambda > 0
+    assert norm.approximate_ee_collision_frequency().to("Hz").magnitude > 0