diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c
index 2b93be63c..c892ce5e6 100644
--- a/gyrokinetic/apps/gk_species.c
+++ b/gyrokinetic/apps/gk_species.c
@@ -79,7 +79,8 @@ gk_species_omegaH_dt(gkyl_gyrokinetic_app *app, struct gk_species *gks, const st
 
 static double
 gk_species_rhs_dynamic(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
 {
   // Gyroaverage the potential if needed.
   species->gyroaverage(app, species, app->field->phi_smooth, species->gyro_phi);
@@ -92,7 +93,7 @@ gk_species_rhs_dynamic(gkyl_gyrokinetic_app *app, struct gk_species *species,
 
   // Damping term.
   gk_species_damping_advance(app, species, &species->damping, app->field->phi_smooth, fin,
-    species->lte.f_lte, rhs, species->cflrate);
+    fbar_in, species->lte.f_lte, rhs, species->cflrate);
 
   // LBO Collisions.
   gk_species_lbo_rhs(app, species, &species->lbo, fin, rhs);
@@ -182,7 +183,8 @@ gk_species_rhs_implicit_dynamic(gkyl_gyrokinetic_app *app, struct gk_species *sp
 
 static double
 gk_species_rhs_static(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
 {
   double omega_cfl = 1/DBL_MAX;
   return app->cfl/omega_cfl;
@@ -1611,6 +1613,7 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st
   gk_species_anomalous_diff_init(app, gks, &gks->anom_diff);
 
   // Damping term -nu*f on RHS.
+  gks->damping = (struct gk_damping) { };
   gk_species_damping_init(app, gks, &gks->damping);
 
   // Function multiplying df/dt.
@@ -1850,6 +1853,9 @@ gk_species_apply_ic(gkyl_gyrokinetic_app *app, struct gk_species *gks, double t0
   if (gks->info.init_from_file.type == 0)
     gk_species_projection_calc(app, gks, &gks->proj_init, gks->f, t0);
 
+  // Initialize fbar for low-pass filter damping
+  gk_species_damping_set_fbar_to_f(gks, &gks->damping, gks->f);
+
   // We are pre-computing source for now as it is time-independent.
   gk_species_source_calc(app, gks, &gks->src, gks->lte.f_lte, t0);
 }
@@ -1900,9 +1906,10 @@ gk_species_apply_ic_cross(gkyl_gyrokinetic_app *app, struct gk_species *gks_self
 
 double
 gk_species_rhs(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms)
 {
-  return species->rhs_func(app, species, fin, rhs, bflux_moms);
+  return species->rhs_func(app, species, fin, fbar_in, rhs, bflux_moms);
 }
 
 double
diff --git a/gyrokinetic/apps/gk_species_damping.c b/gyrokinetic/apps/gk_species_damping.c
index 16669672d..753da55ec 100644
--- a/gyrokinetic/apps/gk_species_damping.c
+++ b/gyrokinetic/apps/gk_species_damping.c
@@ -1,18 +1,91 @@
 #include <assert.h>
 #include <gkyl_gyrokinetic_priv.h>
-#include <gkyl_loss_cone_mask_gyrokinetic.h>
 #include <gkyl_alloc.h>
-#include <gkyl_dg_basis_ops.h>
+
+// Damping state synchronization helpers.
+
+static void
+gk_species_damping_set_fbar_to_f_cellwise_const(const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *f)
+{
+  gkyl_array_set_offset(damp->fbar, 1.0, f, 0);
+  gkyl_array_set_offset(damp->fbar1, 1.0, f, 0);
+  gkyl_array_set_offset(damp->fbarnew, 1.0, f, 0);
+}
+
+static void
+gk_species_damping_set_fbar_to_f_same_basis(const struct gk_species *gks, struct gk_damping *damp,
+  const struct gkyl_array *f)
+{
+  gkyl_array_set(damp->fbar, 1.0, f);
+  gkyl_array_set(damp->fbar1, 1.0, f);
+  gkyl_array_set(damp->fbarnew, 1.0, f);
+}
 
 void
-gk_species_damping_write_disabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_set_fbar_to_f_disabled(const struct gk_species *gks, struct gk_damping *damp,
+  const struct gkyl_array *f)
 {
 }
 
 void
-gk_species_damping_write_enabled(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping *damp,
+  const struct gkyl_array *f)
+{
+  damp->set_fbar_to_f_func(gks, damp, f);
+}
+
+static enum gkyl_array_rio_status
+gk_species_damping_read_fbar_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  const char *fname)
+{
+  return GKYL_ARRAY_RIO_SUCCESS;
+}
+
+static enum gkyl_array_rio_status
+gk_species_damping_read_fbar_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  const char *fname)
+{
+  struct gk_damping *damp = &gks->damping;
+  enum gkyl_array_rio_status status =
+    gkyl_comm_array_read(gks->comm, &gks->grid, &gks->local, damp->fbar_host, fname);
+
+  if (status == GKYL_ARRAY_RIO_SUCCESS) {
+    if (app->use_gpu)
+      gkyl_array_copy(damp->fbar, damp->fbar_host);
+    gkyl_array_copy(damp->fbar1, damp->fbar);
+    gkyl_array_copy(damp->fbarnew, damp->fbar);
+  }
+
+  return status;
+}
+
+enum gkyl_array_rio_status
+gk_species_damping_read_fbar(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  const char *fname)
+{
+  return gks->damping.read_fbar_func(app, gks, fname);
+}
+
+// Damping diagnostics write helpers.
+
+void
+gk_species_damping_write_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
+  int frame)
+{
+}
+
+static void
+gk_species_damping_write_rate_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+}
+
+static void
+gk_species_damping_write_rate_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
 {
-  struct timespec wst = gkyl_wall_clock();
   // DG metadata for damping rate.
   struct gkyl_msgpack_map_elem mpe_drate[] = {
     { .key = "poly_order", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = 0 },
@@ -21,16 +94,19 @@ gk_species_damping_write_enabled(gkyl_gyrokinetic_app* app, struct gk_species *g
     { .key = "time", .elem_type = GKYL_MP_DOUBLE, .dval = tm },
     { .key = "frame", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = frame },
   };
-  int mpe_drate_len = sizeof(mpe_drate)/sizeof(mpe_drate[0]);
+  int mpe_drate_len = sizeof(mpe_drate) / sizeof(mpe_drate[0]);
   // Package metadata.
-  int io_meta_len[] = {gks->io_meta_basic_len, mpe_drate_len, app->gk_geom->io_meta_basic_len};
-  const struct gkyl_msgpack_map_elem* io_meta[] = {gks->io_meta_basic, mpe_drate, app->gk_geom->io_meta_basic};
-  struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union(sizeof(io_meta_len)/sizeof(int), io_meta_len, io_meta);
+  int io_meta_len[] = { gks->io_meta_basic_len, mpe_drate_len,
+    app->gk_geom->io_meta_basic_len };
+  const struct gkyl_msgpack_map_elem *io_meta[] = { gks->io_meta_basic, mpe_drate,
+    app->gk_geom->io_meta_basic };
+  struct gkyl_msgpack_data *mt = gkyl_msgpack_create_union(sizeof(io_meta_len) / sizeof(int),
+    io_meta_len, io_meta);
 
   // Write out the damping rate.
   const char *fmt = "%s-%s_damping_rate_%d.gkyl";
   int sz = gkyl_calc_strlen(fmt, app->name, gks->info.name, frame);
-  char fileNm[sz+1]; // ensures no buffer overflow
+  char fileNm[sz + 1]; // ensures no buffer overflow
   snprintf(fileNm, sizeof fileNm, fmt, app->name, gks->info.name, frame);
 
   // Copy data from device to host before writing it out.
@@ -40,17 +116,112 @@ gk_species_damping_write_enabled(gkyl_gyrokinetic_app* app, struct gk_species *g
   gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt, gks->damping.rate_host, fileNm);
   app->stat.n_io += 1;
 
-  gkyl_msgpack_data_release(mt); 
+  gkyl_msgpack_data_release(mt);
+}
+
+static void
+gk_species_damping_write_fbar_disabled(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+}
+
+static void
+gk_species_damping_write_fbar_cellwise_const(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+  // Metadata from species/geometry. fbar is always stored as a p0 field.
+  struct gkyl_msgpack_map_elem mpe_fbar_p0[] = {
+    { .key = "poly_order", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = 0 },
+    { .key = "basis_type", .elem_type = GKYL_MP_STRING, .cval = "serendipity" },
+    { .key = "Description", .elem_type = GKYL_MP_STRING,
+      .cval = "Low-pass-filtered distribution function." },
+    { .key = "time", .elem_type = GKYL_MP_DOUBLE, .dval = tm },
+    { .key = "frame", .elem_type = GKYL_MP_UNSIGNED_INT, .uval = frame },
+  };
+  int io_meta_fbar_len[] = {
+    gks->io_meta_basic_len,
+    (int)(sizeof(mpe_fbar_p0) / sizeof(mpe_fbar_p0[0])),
+    app->gk_geom->io_meta_basic_len
+  };
+  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = {
+    gks->io_meta_basic,
+    mpe_fbar_p0,
+    app->gk_geom->io_meta_basic
+  };
+  struct gkyl_msgpack_data *mt_fbar = gkyl_msgpack_create_union(
+    sizeof(io_meta_fbar_len) / sizeof(int), io_meta_fbar_len, io_meta_fbar);
+
+  const char *fmt_fbar = "%s-%s_damping_fbar_%d.gkyl";
+  int sz_fbar = gkyl_calc_strlen(fmt_fbar, app->name, gks->info.name, frame);
+  char fileNm_fbar[sz_fbar + 1]; // ensures no buffer overflow
+  snprintf(fileNm_fbar, sizeof fileNm_fbar, fmt_fbar, app->name, gks->info.name, frame);
+
+  // Copy fbar from device to host before writing it out.
+  if (app->use_gpu)
+    gkyl_array_copy(gks->damping.fbar_host, gks->damping.fbar);
+
+  gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt_fbar, gks->damping.fbar_host,
+    fileNm_fbar);
+  app->stat.n_io += 1;
+
+  gkyl_msgpack_data_release(mt_fbar);
+}
+
+static void
+gk_species_damping_write_fbar_same_basis(gkyl_gyrokinetic_app *app, struct gk_species *gks,
+  double tm, int frame)
+{
+  // Metadata from species/geometry using the species phase-space basis.
+  gkyl_msgpack_map_elem_set_double(gks->io_meta_phase_len, gks->io_meta_phase, "time", tm);
+  gkyl_msgpack_map_elem_set_uint(gks->io_meta_phase_len, gks->io_meta_phase, "frame", frame);
+
+  int io_meta_fbar_len[] = {
+    gks->io_meta_phase_len,
+    app->gk_geom->io_meta_basic_len
+  };
+  const struct gkyl_msgpack_map_elem *io_meta_fbar[] = {
+    gks->io_meta_phase,
+    app->gk_geom->io_meta_basic
+  };
+  struct gkyl_msgpack_data *mt_fbar = gkyl_msgpack_create_union(
+    sizeof(io_meta_fbar_len) / sizeof(int), io_meta_fbar_len, io_meta_fbar);
+
+  const char *fmt_fbar = "%s-%s_damping_fbar_%d.gkyl";
+  int sz_fbar = gkyl_calc_strlen(fmt_fbar, app->name, gks->info.name, frame);
+  char fileNm_fbar[sz_fbar + 1]; // ensures no buffer overflow
+  snprintf(fileNm_fbar, sizeof fileNm_fbar, fmt_fbar, app->name, gks->info.name, frame);
+
+  // Copy fbar from device to host before writing it out.
+  if (app->use_gpu)
+    gkyl_array_copy(gks->damping.fbar_host, gks->damping.fbar);
+
+  gkyl_comm_array_write(gks->comm, &gks->grid, &gks->local, mt_fbar, gks->damping.fbar_host,
+    fileNm_fbar);
+  app->stat.n_io += 1;
+
+  gkyl_msgpack_data_release(mt_fbar);
+}
+
+void
+gk_species_damping_write_enabled(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
+  int frame)
+{
+  struct timespec wst = gkyl_wall_clock();
+  gks->damping.write_rate_func(app, gks, tm, frame);
+  gks->damping.write_fbar_func(app, gks, tm, frame);
   app->stat.species_diag_io_tm += gkyl_time_diff_now_sec(wst);
 }
 
 void
-gk_species_damping_write_init_only(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_write_init_only(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm,
+  int frame)
 {
   gk_species_damping_write_enabled(app, gks, tm, frame);
   gks->damping.write_func = gk_species_damping_write_disabled;
 }
 
+// Damping rate projection helpers.
+
 static void
 proj_on_basis_c2p_phase_func(const double *xcomp, double *xphys, void *ctx)
 {
@@ -59,221 +230,398 @@ proj_on_basis_c2p_phase_func(const double *xcomp, double *xphys, void *ctx)
   gkyl_velocity_map_eval_c2p(c2p_ctx->vel_map, &xcomp[cdim], &xphys[cdim]);
 }
 
+static void
+gk_species_damping_project_phase_rate(const struct gkyl_gyrokinetic_app *app,
+  const struct gk_species *gks, int num_quad,
+  void (*rate_profile)(double t, const double *xn, double *fout, void *ctx),
+  void *rate_profile_ctx, struct gkyl_array *rate_host, struct gkyl_array *rate)
+{
+  struct gk_proj_on_basis_c2p_func_ctx proj_on_basis_c2p_ctx; // c2p function context.
+  proj_on_basis_c2p_ctx.cdim = app->cdim;
+  proj_on_basis_c2p_ctx.vdim = gks->local_vel.ndim;
+  proj_on_basis_c2p_ctx.vel_map = gks->vel_map;
+
+  struct gkyl_basis basis_mult;
+  gkyl_cart_modal_serendip(&basis_mult, gks->basis.ndim, 0);
+
+  gkyl_proj_on_basis *projup = gkyl_proj_on_basis_inew(&(struct gkyl_proj_on_basis_inp) {
+    .grid = &gks->grid,
+    .basis = &basis_mult,
+    .num_quad = num_quad,
+    .num_ret_vals = 1,
+    .eval = rate_profile,
+    .ctx = rate_profile_ctx,
+    .c2p_func = proj_on_basis_c2p_phase_func,
+    .c2p_func_ctx = &proj_on_basis_c2p_ctx,
+  });
+  gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, rate_host);
+  gkyl_proj_on_basis_release(projup);
+  gkyl_array_copy(rate, rate_host);
+
+  if (num_quad == 1)
+    gkyl_array_scale_range(rate, 1.0 / pow(sqrt(2.0), gks->grid.ndim), &gks->local);
+}
+
+// Damping RHS assembly dispatch helpers.
+
+void
+gk_species_damping_advance_disabled(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+}
+
+void
+gk_species_damping_advance_user_input(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  gkyl_array_set(f_buffer, 1.0, fin);
+  gkyl_array_scale_by_cell(f_buffer, damp->rate);
+
+  // Add damping to f and the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
+
+static void
+gk_species_damping_advance_low_pass_filter_cellwise_const(gkyl_gyrokinetic_app *app,
+  const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  // Compute f - fbar, where fbar is stored as a cellwise-constant p0 field.
+  gkyl_array_set(f_buffer, 1.0, fin); // f_buffer = fin
+  gkyl_array_accumulate_offset(f_buffer, -1.0, fbar_in, 0); // f_buffer = fin - fbar
+  gkyl_array_scale_by_cell(f_buffer, damp->rate); // f_buffer = rate * (fin - fbar)
+
+  // Add damping term to RHS: df/dt -= rate * (f - fbar)
+  // Add to the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
+
+static void
+gk_species_damping_advance_low_pass_filter_same_basis(gkyl_gyrokinetic_app *app,
+  const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  // Compute f - fbar, where fbar uses the same basis as f.
+  gkyl_array_set(f_buffer, 1.0, fin); // f_buffer = fin
+  gkyl_array_accumulate(f_buffer, -1.0, fbar_in); // f_buffer = fin - fbar
+  gkyl_array_scale_by_cell(f_buffer, damp->rate); // f_buffer = rate * (fin - fbar)
+
+  // Add damping term to RHS: df/dt -= rate * (f - fbar)
+  // Add to the CFL frequency.
+  gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  gkyl_array_accumulate(cflrate, 1.0, damp->rate);
+}
+
+void
+gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks,
+  struct gk_damping *damp,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
+  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+{
+  struct timespec wst = gkyl_wall_clock();
+
+  damp->advance_func(app, gks, damp, phi, fin, fbar_in, f_buffer, rhs, cflrate);
+
+  app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
+}
+
+void
+gk_species_damping_write(gkyl_gyrokinetic_app *app, struct gk_species *gks, double tm, int frame)
+{
+  gks->damping.write_func(app, gks, tm, frame);
+}
+
+// Low-pass filter fbar RHS helpers.
+
+void
+gk_species_damping_calc_fbar_rhs_disabled(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+}
+
+static void
+gk_species_damping_calc_fbar_rhs_cellwise_const(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  // rhs_fbar = rate * (f - fbar)
+  gkyl_array_set_offset(rhs_fbar, 1.0, fin, 0);
+  gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
+  gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
+}
+
+static void
+gk_species_damping_calc_fbar_rhs_same_basis(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  // rhs_fbar = rate * (f - fbar)
+  gkyl_array_set(rhs_fbar, 1.0, fin);
+  gkyl_array_accumulate(rhs_fbar, -1.0, fbar_in);
+  gkyl_array_scale_by_cell(rhs_fbar, damp->rate);
+}
+
+void
+gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar)
+{
+  damp->calc_fbar_rhs_func(damp, fin, fbar_in, rhs_fbar);
+}
+
+// Low-pass filter stage update helpers.
+
+static void
+gk_species_damping_forward_euler_disabled(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt)
+{
+}
+
+static void
+gk_species_damping_forward_euler_enabled(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt)
+{
+  gk_species_damping_calc_fbar_rhs(&gks->damping, fin, fbar_in, fbar_out);
+  gk_species_step_f(gks, fbar_out, dt, fbar_in);
+}
+
+void
+gk_species_damping_forward_euler(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt)
+{
+  gks->damping.forward_euler_func(gks, fin, fbar_in, fbar_out, dt);
+}
+
+// Low-pass filter RK stage combine/copy helpers.
+
+static void
+gk_species_damping_combine_disabled(struct gk_species *gks, struct gkyl_array *fout,
+  double c1, const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng)
+{
+}
+
+static void
+gk_species_damping_combine_enabled(struct gk_species *gks, struct gkyl_array *fout,
+  double c1, const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng)
+{
+  gk_species_combine(gks, fout, c1, f1, c2, f2, rng);
+}
+
+void
+gk_species_damping_combine(struct gk_species *gks, struct gkyl_array *fout, double c1,
+  const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng)
+{
+  gks->damping.combine_func(gks, fout, c1, f1, c2, f2, rng);
+}
+
+static void
+gk_species_damping_copy_range_disabled(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range)
+{
+}
+
+static void
+gk_species_damping_copy_range_enabled(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range)
+{
+  gk_species_copy_range(gks, fout, fin, range);
+}
+
+void
+gk_species_damping_copy_range(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range)
+{
+  gks->damping.copy_func(gks, fout, fin, range);
+}
+
+// Damping object lifecycle.
+
 void
 gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks,
   struct gk_damping *damp)
 {
   damp->type = gks->info.damping.type;
   damp->evolve = false; // Whether the rate is time dependent.
+  damp->write_rate = gks->info.damping.write_rate;
+  damp->write_fbar = gks->info.damping.write_fbar;
+  damp->cellwise_const = gks->info.damping.cellwise_const;
+  damp->do_not_reset_fbar = gks->info.damping.do_not_reset_fbar;
+  const double rate_const = gks->info.damping.rate_const;
 
   int num_quad = gks->info.damping.num_quad? gks->info.damping.num_quad : 1; // Default is a p=0 mask.
   assert(num_quad == 1); // MF 2025/06/11: Limited to this for now.
 
   // Default function pointers.
   damp->write_func = gk_species_damping_write_disabled;
+  damp->write_rate_func = gk_species_damping_write_rate_disabled;
+  damp->write_fbar_func = gk_species_damping_write_fbar_disabled;
+  damp->read_fbar_func = gk_species_damping_read_fbar_disabled;
+  damp->advance_func = gk_species_damping_advance_disabled;
+  damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_disabled;
+  damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_disabled;
+  damp->forward_euler_func = gk_species_damping_forward_euler_disabled;
+  damp->combine_func = gk_species_damping_combine_disabled;
+  damp->copy_func = gk_species_damping_copy_range_disabled;
+
+  // Default pointers for fbar arrays
+  damp->fbar = 0;
+  damp->fbar1 = 0;
+  damp->fbarnew = 0;
+  damp->fbar_host = 0;
+  damp->fbar_initialized = false;
+
+  if (!damp->type) {
+    return;
+  }
 
-  if (damp->type) {
-    // Allocate rate array.
-    damp->rate = mkarr(app->use_gpu, num_quad==1? 1 : gks->basis.num_basis, gks->local_ext.volume);
-    damp->rate_host = damp->rate;
-    if (app->use_gpu)
-      damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size); 
-
-    if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
-      struct gk_proj_on_basis_c2p_func_ctx proj_on_basis_c2p_ctx; // c2p function context.
-      proj_on_basis_c2p_ctx.cdim = app->cdim;
-      proj_on_basis_c2p_ctx.vdim = gks->local_vel.ndim;
-      proj_on_basis_c2p_ctx.vel_map = gks->vel_map;
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_inew( &(struct gkyl_proj_on_basis_inp) {
-          .grid = &gks->grid,
-          .basis = &gks->basis,
-          .num_quad = num_quad,
-          .num_ret_vals = 1,
-          .eval = gks->info.damping.rate_profile,
-          .ctx = gks->info.damping.rate_profile_ctx,
-          .c2p_func = proj_on_basis_c2p_phase_func,
-          .c2p_func_ctx = &proj_on_basis_c2p_ctx,
-        }
-      );
-      gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, damp->rate_host);
-      gkyl_proj_on_basis_release(projup);
-      gkyl_array_copy(damp->rate, damp->rate_host);
-
-      if (num_quad == 1)
-        gkyl_array_scale_range(damp->rate, 1.0/pow(sqrt(2.0),gks->grid.ndim), &gks->local);
-    }
-    else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
-      damp->evolve = true; // Since the loss cone boundary is proportional to phi(t).
-
-      // Maximum bmag and its location.
-      // NOTE: if the same max bmag occurs at multiple locations,
-      // bmag_max_coord may have different values on different MPI processes.
-      double bmag_max_coord_ho[GKYL_MAX_CDIM];
-      double bmag_max_ho = gkyl_gk_geometry_reduce_arg_bmag(app->gk_geom, GKYL_MAX, bmag_max_coord_ho);
-      double bmag_max_local = bmag_max_ho;
-      double bmag_max_global;
-      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, &bmag_max_local, &bmag_max_global);
-      double bmag_max_coord_local[app->cdim], bmag_max_coord_global[app->cdim];
-      if (fabs(bmag_max_ho - bmag_max_global) < 1e-16) {
-        for (int d=0; d<app->cdim; d++)
-          bmag_max_coord_local[d] = bmag_max_coord_ho[d];
-      }
-      else {
-        for (int d=0; d<app->cdim; d++)
-          bmag_max_coord_local[d] = -DBL_MAX;
-      }
-      gkyl_comm_allreduce_host(app->comm, GKYL_DOUBLE, GKYL_MAX, app->cdim, bmag_max_coord_local, bmag_max_coord_global);
+  if (damp->write_rate) {
+    damp->write_rate_func = gk_species_damping_write_rate_enabled;
+  }
 
-      if (app->use_gpu) {
-        damp->bmag_max = gkyl_cu_malloc(sizeof(double));
-        damp->bmag_max_coord = gkyl_cu_malloc(app->cdim*sizeof(double));
-	gkyl_cu_memcpy(damp->bmag_max, &bmag_max_global, sizeof(double), GKYL_CU_MEMCPY_H2D);
-	gkyl_cu_memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double), GKYL_CU_MEMCPY_H2D);
+  // Allocate rate array.
+  damp->rate = mkarr(app->use_gpu, 1, gks->local_ext.volume);
+  damp->rate_host = damp->rate;
+  if (app->use_gpu) {
+    damp->rate_host = mkarr(false, damp->rate->ncomp, damp->rate->size);
+  }
+
+  if (gks->info.damping.rate_profile) {
+    gk_species_damping_project_phase_rate(app, gks, num_quad, gks->info.damping.rate_profile,
+      gks->info.damping.rate_profile_ctx, damp->rate_host, damp->rate);
+    gkyl_array_scale(damp->rate, rate_const == 0.0 ? 1.0 : rate_const);
+  }
+  else {
+    gkyl_array_clear(damp->rate, rate_const);
+  }
+
+  if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
+    // USER_INPUT supports either a projected profile, or a uniform constant if no profile is provided.
+    damp->advance_func = gk_species_damping_advance_user_input;
+  }
+  else if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+    // LOW_PASS_FILTER supports either a projected phase-space profile or a uniform constant.
+    damp->evolve = true;
+    
+    // Allocate filtered distribution function array in the requested basis.
+    int fbar_ncomp = damp->cellwise_const ? 1 : gks->basis.num_basis;
+    damp->fbar = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
+    damp->fbar1 = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
+    damp->fbarnew = mkarr(app->use_gpu, fbar_ncomp, gks->local_ext.volume);
+    damp->fbar_host = damp->fbar;
+    damp->fbar_initialized = true;
+    gkyl_array_clear(damp->fbar, 0.0);
+
+    if (damp->write_fbar) {
+      damp->read_fbar_func = gk_species_damping_read_fbar_enabled;
+      if (damp->cellwise_const) {
+        damp->write_fbar_func = gk_species_damping_write_fbar_cellwise_const;
       }
       else {
-        damp->bmag_max = gkyl_malloc(sizeof(double));
-        damp->bmag_max_coord = gkyl_malloc(app->cdim*sizeof(double));
-	memcpy(damp->bmag_max, &bmag_max_global, sizeof(double));
-	memcpy(damp->bmag_max_coord, bmag_max_coord_ho, app->cdim*sizeof(double));
+        damp->write_fbar_func = gk_species_damping_write_fbar_same_basis;
       }
-
-      // Electrostatic potential at bmag_max_coord.
       if (app->use_gpu) {
-        damp->phi_m = gkyl_cu_malloc(sizeof(double));
-        damp->phi_m_global = gkyl_cu_malloc(sizeof(double));
+        damp->fbar_host = mkarr(false, damp->fbar->ncomp, damp->fbar->size);
       }
-      else {
-        damp->phi_m = gkyl_malloc(sizeof(double));
-        damp->phi_m_global = gkyl_malloc(sizeof(double));
-      }
-
-      // Operator that projects the loss cone mask.
-      struct gkyl_loss_cone_mask_gyrokinetic_inp inp_proj = {
-        .phase_grid = &gks->grid,
-        .conf_basis = &app->basis,
-        .phase_basis = &gks->basis,
-        .conf_range =  &app->local,
-        .conf_range_ext = &app->local_ext,
-        .vel_range = &gks->local_vel, 
-        .vel_map = gks->vel_map,
-        .bmag = app->gk_geom->geo_int.bmag,
-        .bmag_max = damp->bmag_max,
-        .bmag_max_loc = damp->bmag_max_coord,
-        .mass = gks->info.mass,
-        .charge = gks->info.charge,
-        .num_quad = num_quad,
-        .use_gpu = app->use_gpu,
-      };
-      damp->lcm_proj_op = gkyl_loss_cone_mask_gyrokinetic_inew( &inp_proj );
-
-      // Project the conf-space rate profile provided.
-      struct gkyl_array *scale_prof_high_order = mkarr(app->use_gpu, gks->basis.num_basis, gks->local_ext.volume);
-      struct gkyl_array *scale_prof_high_order_ho = app->use_gpu? mkarr(false, scale_prof_high_order->ncomp, scale_prof_high_order->size)
-                                                     : gkyl_array_acquire(scale_prof_high_order);
-      
-      gkyl_proj_on_basis *projup = gkyl_proj_on_basis_new(&gks->grid, &gks->basis, num_quad, 1, 
-        gks->info.damping.rate_profile, gks->info.damping.rate_profile_ctx);
-      gkyl_proj_on_basis_advance(projup, 0.0, &gks->local, scale_prof_high_order_ho);
-      gkyl_proj_on_basis_release(projup);
-      gkyl_array_copy(scale_prof_high_order, scale_prof_high_order_ho);
-
-      damp->scale_prof = mkarr(app->use_gpu, num_quad == 1? 1 : gks->basis.num_basis, gks->local_ext.volume);
-      gkyl_array_set_offset(damp->scale_prof, pow(sqrt(2.0),gks->grid.ndim), scale_prof_high_order, 0);
-
-      gkyl_array_release(scale_prof_high_order_ho);
-      gkyl_array_release(scale_prof_high_order);
-
-      // Compute the initial damping rate (assuming phi=0 because phi hasn't been computed).
-      // Find the potential at the mirror throat.
-      gkyl_dg_basis_ops_eval_array_at_coord_comp(app->field->phi_smooth, damp->bmag_max_coord,
-        app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
-      gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
-      // Project the loss cone mask.
-      gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
-        app->field->phi_smooth, damp->phi_m_global, damp->rate);
-      // Multiply by the user's scaling profile.
-      gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
     }
 
-    // Set function pointers chosen at runtime.
-    if (damp->evolve) {
-      damp->write_func = gk_species_damping_write_enabled;
+    if (damp->cellwise_const) {
+      damp->advance_func = gk_species_damping_advance_low_pass_filter_cellwise_const;
+      damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_cellwise_const;
+      damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_cellwise_const;
     }
     else {
-      damp->write_func = gk_species_damping_write_init_only;
+      damp->advance_func = gk_species_damping_advance_low_pass_filter_same_basis;
+      damp->set_fbar_to_f_func = gk_species_damping_set_fbar_to_f_same_basis;
+      damp->calc_fbar_rhs_func = gk_species_damping_calc_fbar_rhs_same_basis;
     }
+    damp->forward_euler_func = gk_species_damping_forward_euler_enabled;
+    damp->combine_func = gk_species_damping_combine_enabled;
+    damp->copy_func = gk_species_damping_copy_range_enabled;
+  }
+
+  if (damp->evolve) {
+    damp->write_func = gk_species_damping_write_enabled;
+  }
+  else {
+    damp->write_func = gk_species_damping_write_init_only;
   }
 }
 
 void
-gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
-  struct gkyl_array *rhs, struct gkyl_array *cflrate)
+gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp)
 {
-  if (damp->type) {
-    struct timespec wst = gkyl_wall_clock();
-    if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
-      gkyl_array_set(f_buffer, 1.0, fin);
-      gkyl_array_scale_by_cell(f_buffer, damp->rate);
-      gkyl_array_accumulate(rhs, -1.0, f_buffer);
-    }
-    else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
-      // Find the potential at the mirror throat.
-      gkyl_dg_basis_ops_eval_array_at_coord_comp(phi, damp->bmag_max_coord,
-        app->basis_on_dev, &app->grid, &app->local, damp->phi_m);
-      gkyl_comm_allreduce(app->comm, GKYL_DOUBLE, GKYL_MAX, 1, damp->phi_m, damp->phi_m_global);
-
-      // Project the loss cone mask.
-      gkyl_loss_cone_mask_gyrokinetic_advance(damp->lcm_proj_op, &gks->local, &app->local,
-        phi, damp->phi_m_global, damp->rate);
-
-      // Assemble the damping term -scale_prof * mask * f.
-      gkyl_array_set(f_buffer, 1.0, fin);
-      gkyl_array_scale_by_cell(damp->rate, damp->scale_prof);
-      gkyl_array_scale_by_cell(f_buffer, damp->rate);
-      gkyl_array_accumulate(rhs, -1.0, f_buffer);
+  // Release remembered fbar state even when damping is currently disabled.
+  if (damp->fbar_initialized) {
+    gkyl_array_release(damp->fbar);
+  }
 
+  if (!damp->type) {
+    return;
+  }
+  gkyl_array_release(damp->rate);
+  if (app->use_gpu) {
+    gkyl_array_release(damp->rate_host);
+  }
+  if (damp->type == GKYL_GK_DAMPING_LOW_PASS_FILTER) {
+    gkyl_array_release(damp->fbar1);
+    gkyl_array_release(damp->fbarnew);
+    if (app->use_gpu && damp->write_fbar) {
+      gkyl_array_release(damp->fbar_host);
     }
-
-    // Add the frequency to the CFL frequency.
-    gkyl_array_accumulate(cflrate, 1.0, damp->rate);
-
-    app->stat.species_damp_tm += gkyl_time_diff_now_sec(wst);
   }
 }
 
 void
-gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame)
+gk_species_damping_reset(gkyl_gyrokinetic_app *app, double tm, struct gk_species *gks,
+  struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp)
 {
-  gks->damping.write_func(app, gks, tm, frame);
-}
+  bool has_fbar_backup = false;
+  struct gkyl_array *fbar_backup = 0;
+
+  // Create a reset-local backup if requested and there is restorable state.
+  if (damp_inp.do_not_reset_fbar) {
+    if (damp->fbar_initialized) {
+      has_fbar_backup = true;
+      fbar_backup = mkarr(app->use_gpu, damp->fbar->ncomp, damp->fbar->size);
+      gkyl_array_copy(fbar_backup, damp->fbar);
+    }
+  }
 
-void
-gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp)
-{
-  if (damp->type) {
-    gkyl_array_release(damp->rate);
-    if (app->use_gpu)
-      gkyl_array_release(damp->rate_host);
+  gk_species_damping_release(app, damp);
+
+  gks->info.damping = damp_inp;
+  gk_species_damping_init(app, gks, damp);
 
-    if (damp->type == GKYL_GK_DAMPING_USER_INPUT) {
-      // Nothing to release.
+  gk_species_damping_set_fbar_to_f(gks, damp, gks->f);
+
+  if (has_fbar_backup) {
+    if (damp->fbar_initialized) {
+      gkyl_array_copy(damp->fbar, fbar_backup);
+      gkyl_array_copy(damp->fbar1, fbar_backup);
+      gkyl_array_copy(damp->fbarnew, fbar_backup);
     }
-    else if (damp->type == GKYL_GK_DAMPING_LOSS_CONE) {
-      if (app->use_gpu) {
-        gkyl_cu_free(damp->bmag_max);
-        gkyl_cu_free(damp->bmag_max_coord);
-        gkyl_cu_free(damp->phi_m);
-        gkyl_cu_free(damp->phi_m_global);
-      }
-      else {
-        gkyl_free(damp->bmag_max);
-        gkyl_free(damp->bmag_max_coord);
-        gkyl_free(damp->phi_m);
-        gkyl_free(damp->phi_m_global);
-      }
-      gkyl_loss_cone_mask_gyrokinetic_release(damp->lcm_proj_op);
-      gkyl_array_release(damp->scale_prof);
+    else {
+      damp->fbar = gkyl_array_acquire(fbar_backup);
+      damp->fbar_initialized = true;
     }
+    gkyl_array_release(fbar_backup);
   }
 }
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index e04e1510b..1a59989f8 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -360,20 +360,25 @@ struct gkyl_gyrokinetic_positivity {
 enum gkyl_gyrokinetic_damping_type {
   GKYL_GK_DAMPING_NONE = 0,
   GKYL_GK_DAMPING_USER_INPUT,
-  GKYL_GK_DAMPING_LOSS_CONE,
+  GKYL_GK_DAMPING_LOW_PASS_FILTER,
 };
 
 struct gkyl_gyrokinetic_damping {
   // Add a damping term to the RHS of the gyrokinetic equation
   //   df/dt = - rate(z) * f
   // with the function rate(z) being:
-  //   - a function given by the user (type = GKYL_PROPORTIONAL_TERM_USER_INPUT).
-  //   - I_loss(z) * scale_factor * scale_profile(z), where I_loss(z) is =1 in the loss
-  //     cone and 0 in the confined region (type = GKYL_PROPORTIONAL_TERM_LOSS_CONE).
+  //   - a function given by the user (type = GKYL_GK_DAMPING_USER_INPUT).
+  //   Alternativly, a low-pass filter can be applied (type = GKYL_GK_DAMPING_LOW_PASS_FILTER), which projects f to a lower-dimensional function fbar(z) and applies the damping as df/dt = - rate(z) * (f - fbar).
   enum gkyl_gyrokinetic_damping_type type;
+
+  double rate_const;
   void (*rate_profile)(double t, const double *xn, double *fout, void *ctx);
   void *rate_profile_ctx; // Context for rate_profile function.
   int num_quad; // Number of quadrature points in each direction to use in projecting the rate.
+  bool write_rate; // Whether to write damping-rate diagnostics.
+  bool write_fbar; // For low-pass filter damping, write fbar diagnostics each output frame.
+  bool cellwise_const; // For low-pass filter damping, whether fbar uses a single value per cell.
+  bool do_not_reset_fbar; // On runtime reset, preserve the existing low-pass-filter fbar instead of reinitializing from f.
 };
 
 // Types of df/dt multipliers: M(x,v,t) modifies df/dt -> M * df/dt.
@@ -1489,6 +1494,18 @@ void gkyl_gyrokinetic_app_reset_species_collisionless(gkyl_gyrokinetic_app* app,
 void gkyl_gyrokinetic_app_reset_species_positivity(gkyl_gyrokinetic_app* app, double tm,
   const char *species_name, struct gkyl_gyrokinetic_positivity pos_inp);
 
+/**
+ * Reset the damping for a given species.
+ *
+ * @param app App object.
+ * @param tm Time-stamp.
+ * @param species_name Name of the species to reset.
+ * @param damping_inp Input parameters for damping terms.
+ */
+void
+gkyl_gyrokinetic_app_reset_species_damping(gkyl_gyrokinetic_app* app, double tm,
+  const char *species_name, struct gkyl_gyrokinetic_damping damping_inp);
+
 /**
  * Reset the field solver.
  *
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 561516097..5dafb412f 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -216,7 +216,7 @@ struct gk_rad_drag {
   struct gk_neut_species *collide_with_neut[GKYL_MAX_SPECIES]; // Neutral cross-species we collide with.
   int collide_with_idx[2*GKYL_MAX_SPECIES]; // Index of species we collide with.
   bool is_neut_species[2*GKYL_MAX_SPECIES]; // Flag of whether neutral or Gk species.
-  
+
   // Drag coefficients in vparallel and mu for each species being collided with.
   struct gkyl_gk_rad_drag *vnu_surf;
   struct gkyl_gk_rad_drag *vnu;
@@ -841,15 +841,38 @@ struct gk_source {
 struct gk_damping {
   enum gkyl_gyrokinetic_damping_type type; // Type of damping term.
   bool evolve; // Whether the source is time dependent.
+  bool write_rate; // Whether to write damping-rate diagnostics.
+  bool write_fbar; // Whether to write low-pass-filter fbar diagnostics.
+  bool cellwise_const; // Whether low-pass-filter fbar is stored as one value per cell.
+  bool do_not_reset_fbar; // Whether runtime reset should preserve existing fbar.
   struct gkyl_array *rate; // Damping rate.
   struct gkyl_array *rate_host; // Host copy for use in IO and projecting.
-  struct gkyl_loss_cone_mask_gyrokinetic *lcm_proj_op; // Operator that projects the loss cone mask.
-  double *bmag_max; // Maximum magnetic field amplitude.
-  double *bmag_max_coord; // Location of bmag_max.
-  double *phi_m, *phi_m_global; // Electrostatic potential at bmag_max.
-  struct gkyl_array *scale_prof; // Conf-space scaling factor profile.
+  // Low-pass filter variables
+  bool fbar_initialized; // Whether fbar is allocated/owned and should be released.
+  struct gkyl_array *fbar; // Filtered/averaged distribution function.
+  struct gkyl_array *fbar1, *fbarnew; // SSPRK3 stage arrays for filtered distribution.
+  struct gkyl_array *fbar_host; // Host copy of fbar for use in IO.
   // Functions chosen at runtime.
   void (*write_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
+  void (*write_rate_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
+  void (*write_fbar_func)(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
+  enum gkyl_array_rio_status (*read_fbar_func)(gkyl_gyrokinetic_app* app,
+    struct gk_species *gks, const char *fname);
+  void (*advance_func)(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp, 
+    const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+    struct gkyl_array *f_buffer,
+    struct gkyl_array *rhs, struct gkyl_array *cflrate);
+  void (*set_fbar_to_f_func)(const struct gk_species *gks, struct gk_damping *damp,
+    const struct gkyl_array *f);
+  void (*calc_fbar_rhs_func)(const struct gk_damping *damp,
+    const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
+  void (*forward_euler_func)(struct gk_species *gks, const struct gkyl_array *fin,
+    const struct gkyl_array *fbar_in, struct gkyl_array *fbar_out, double dt);
+  void (*combine_func)(struct gk_species *gks, struct gkyl_array *fout, double c1,
+    const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+    const struct gkyl_range *rng);
+  void (*copy_func)(struct gk_species *gks, struct gkyl_array *fout,
+    const struct gkyl_array *fin, const struct gkyl_range *range);
 };
 
 struct gk_fdot_multiplier_comp {
@@ -1123,7 +1146,8 @@ struct gk_species {
 
   // Pointer to various functions selected at runtime.
   double (*rhs_func)(gkyl_gyrokinetic_app *app, struct gk_species *species,
-    const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
+    const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+    struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
   double (*rhs_implicit_func)(gkyl_gyrokinetic_app *app, struct gk_species *species,
     const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms, double dt);
   void (*bc_func)(gkyl_gyrokinetic_app *app, const struct gk_species *species,
@@ -1514,7 +1538,7 @@ struct gkyl_gyrokinetic_app {
 
   gkyl_dynvec dts; // Record time step over time.
   bool is_first_dt_write_call; // Flag for integrated moments dynvec written first time.
-  
+
   bool is_multib; // Is this a block in a multiblock sim?
 };
 
@@ -2911,6 +2935,18 @@ void gk_species_source_release(const struct gkyl_gyrokinetic_app *app, const str
  */
 void gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species *gks, struct gk_damping *damp);
 
+/**
+ * Initialize filtered distribution function from initial distribution.
+ * Should be called after initial conditions are applied to f.
+ *
+ * @param gks Pointer to species.
+ * @param damp Species damping object.
+ * @param f Initial distribution function.
+ */
+void gk_species_damping_set_fbar_to_f(const struct gk_species *gks, struct gk_damping *damp, 
+  const struct gkyl_array *f);
+
+
 /**
  * Compute species applied source term.
  *
@@ -2919,12 +2955,14 @@ void gk_species_damping_init(struct gkyl_gyrokinetic_app *app, struct gk_species
  * @param damp Species damping object.
  * @param phi Current electrostatic potential.
  * @param fin Current distribution function.
+ * @param fbar_in Filtered distribution corresponding to fin RK stage.
  * @param f_buffer Phase-space buffer.
  * @param rhs df/dt damping term gets added to.
  * @param cflrate CFL frequency in phase space.
  */
 void gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_species *gks, struct gk_damping *damp,
-  const struct gkyl_array *phi, const struct gkyl_array *fin, struct gkyl_array *f_buffer,
+  const struct gkyl_array *phi, const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *f_buffer,
   struct gkyl_array *rhs, struct gkyl_array *cflrate);
 
 /**
@@ -2937,6 +2975,68 @@ void gk_species_damping_advance(gkyl_gyrokinetic_app *app, const struct gk_speci
  */
 void gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks, double tm, int frame);
 
+/**
+ * Read the filtered distribution function used by low-pass-filter damping.
+ *
+ * This is a runtime-dispatched no-op unless fbar output is enabled. On a
+ * successful read, all filtered-distribution RK stage arrays are synchronized
+ * to the restarted state.
+ *
+ * @param app Gyrokinetic app object.
+ * @param gks Species object.
+ * @param fname Name of the fbar restart file.
+ * @return Status of the file read.
+ */
+enum gkyl_array_rio_status gk_species_damping_read_fbar(gkyl_gyrokinetic_app *app,
+  struct gk_species *gks, const char *fname);
+
+/**
+ * Compute filtered distribution RHS for low-pass filter damping:
+ *   d(fbar)/dt = rate * (f - fbar).
+ *
+ * @param damp Species damping object.
+ * @param fin Current distribution function f.
+ * @param fbar_in Current filtered distribution function fbar.
+ * @param rhs_fbar RHS buffer for fbar.
+ */
+void gk_species_damping_calc_fbar_rhs(const struct gk_damping *damp,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in, struct gkyl_array *rhs_fbar);
+
+/**
+ * Take a forward-Euler substep for the filtered distribution.
+ *
+ * This is a runtime-dispatched no-op unless the damping type enables
+ * filtered-state evolution (currently low-pass filter).
+ *
+ * @param gks Species object.
+ * @param fin Current distribution function f.
+ * @param fbar_in Current filtered distribution state.
+ * @param fbar_out Output filtered distribution state after FE substep.
+ * @param dt Time-step for substep.
+ */
+void gk_species_damping_forward_euler(struct gk_species *gks,
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *fbar_out, double dt);
+
+/**
+ * Combine filtered distributions like gk_species_combine.
+ *
+ * This is runtime-dispatched and a no-op unless the damping type enables
+ * filtered-state evolution (currently low-pass filter).
+ */
+void gk_species_damping_combine(struct gk_species *gks, struct gkyl_array *fout, double c1,
+  const struct gkyl_array *f1, double c2, const struct gkyl_array *f2,
+  const struct gkyl_range *rng);
+
+/**
+ * Copy filtered distribution ranges like gk_species_copy_range.
+ *
+ * This is runtime-dispatched and a no-op unless the damping type enables
+ * filtered-state evolution (currently low-pass filter).
+ */
+void gk_species_damping_copy_range(struct gk_species *gks, struct gkyl_array *fout,
+  const struct gkyl_array *fin, const struct gkyl_range *range);
+
 /**
  * Release species damping object.
  *
@@ -2945,6 +3045,18 @@ void gk_species_damping_write(gkyl_gyrokinetic_app* app, struct gk_species *gks,
  */
 void gk_species_damping_release(const struct gkyl_gyrokinetic_app *app, const struct gk_damping *damp);
 
+/**
+ * Reset species damping object.
+ *
+ * @param app gyrokinetic app object.
+ * @param tm Time.
+ * @param gks Species object.
+ * @param damp Species damping object.
+ * @param damp_inp New damping input.
+ */
+void gk_species_damping_reset(gkyl_gyrokinetic_app* app, double tm, struct gk_species *gks,
+  struct gk_damping *damp, struct gkyl_gyrokinetic_damping damp_inp);
+
 /** gk_species_fdot_multiplier API */
 
 /**
@@ -3173,12 +3285,14 @@ void gk_species_apply_ic_cross(gkyl_gyrokinetic_app *app, struct gk_species *spe
  * @param app gyrokinetic app object.
  * @param species Pointer to species.
  * @param fin Input distribution function.
+ * @param fbar_in Input filtered distribution for the current RK stage.
  * @param rhs On output, the RHS from the species object.
  * @param bflux_moms Output boundary flux moments.
  * @return Maximum stable time-step.
  */
 double gk_species_rhs(gkyl_gyrokinetic_app *app, struct gk_species *species,
-  const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
+  const struct gkyl_array *fin, const struct gkyl_array *fbar_in,
+  struct gkyl_array *rhs, struct gkyl_array **bflux_moms);
 
 /**
  * Compute the *implicit* RHS from species distribution function
@@ -4291,6 +4405,7 @@ void gyrokinetic_calc_field_and_apply_bc(gkyl_gyrokinetic_app* app, double tcurr
  * @param tcurr Current simulation time.
  * @param dt Suggested time step.
  * @param fin Input array of charged-species distribution functions.
+ * @param fbar_in Input array of charged-species filtered distributions.
  * @param fout Output array of charged-species change in distribution functions.
  * @param bflux_out Output array of charged-species boundary fluxes.
  * @param fin_neut Input array of neutral-species distribution functions.
@@ -4299,7 +4414,8 @@ void gyrokinetic_calc_field_and_apply_bc(gkyl_gyrokinetic_app* app, double tcurr
  * @param st Time stepping status object.
  */
 void gyrokinetic_rhs(gkyl_gyrokinetic_app* app, double tcurr, double dt,
-  const struct gkyl_array *fin[], struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
+  const struct gkyl_array *fin[], const struct gkyl_array *fbar_in[],
+  struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], struct gkyl_array **bflux_out_neut[],
   struct gkyl_update_status *st); 
 
diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c
index 526aeebc6..c00491364 100644
--- a/gyrokinetic/apps/gyrokinetic.c
+++ b/gyrokinetic/apps/gyrokinetic.c
@@ -2314,7 +2314,8 @@ gkyl_gyrokinetic_app_write(gkyl_gyrokinetic_app* app, double tm, int frame)
 
 void
 gyrokinetic_rhs(gkyl_gyrokinetic_app* app, double tcurr, double dt,
-  const struct gkyl_array *fin[], struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
+  const struct gkyl_array *fin[], const struct gkyl_array *fbar_in[],
+  struct gkyl_array *fout[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], struct gkyl_array **bflux_out_neut[], 
   struct gkyl_update_status *st)
 {
@@ -2362,7 +2363,7 @@ gyrokinetic_rhs(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   // Compute df/dt (not including sources).
   for (int i=0; i<app->num_species; ++i) {
     struct gk_species *gk_s = &app->species[i];
-    double dt1 = gk_species_rhs(app, gk_s, fin[i], fout[i], bflux_out[i]);
+    double dt1 = gk_species_rhs(app, gk_s, fin[i], fbar_in[i], fout[i], bflux_out[i]);
     dtmin = fmin(dtmin, dt1);
   }
   for (int i=0; i<app->num_neut_species; ++i) {
@@ -3331,6 +3332,8 @@ gkyl_gyrokinetic_app_from_file_species(gkyl_gyrokinetic_app *app, int sidx,
       gkyl_array_copy(gk_s->f, gk_s->f_host);
 
     if (rstat.io_status == GKYL_ARRAY_RIO_SUCCESS) {
+      // Use f as a valid fallback filtered state when no fbar file is configured.
+      gk_species_damping_set_fbar_to_f(gk_s, &gk_s->damping, gk_s->f);
       gk_species_source_calc(app, gk_s, &gk_s->src, gk_s->lte.f_lte, 0.0);
       // Read volume and time integrated boundary flux diagnostics.
       gk_species_bflux_read_voltime_integrated_mom(app, gk_s, &gk_s->bflux);
@@ -3384,6 +3387,13 @@ gkyl_gyrokinetic_app_from_frame_species(gkyl_gyrokinetic_app *app, int sidx, int
   struct gkyl_app_restart_status rstat = gkyl_gyrokinetic_app_from_file_species(app, sidx, fileNm.str);
   cstr_drop(&fileNm);
 
+  if (rstat.io_status == GKYL_ARRAY_RIO_SUCCESS) {
+    cstr fbarFileNm = cstr_from_fmt("%s-%s_damping_fbar_%d.gkyl", app->name,
+      gk_s->info.name, frame);
+    rstat.io_status = gk_species_damping_read_fbar(app, gk_s, fbarFileNm.str);
+    cstr_drop(&fbarFileNm);
+  }
+
   // Append to existing integrated diagnostics.
   app->is_first_dt_write_call = false;
   gk_s->is_first_integ_write_call = false;
@@ -3431,13 +3441,19 @@ gkyl_gyrokinetic_app_from_frame_neut_species(gkyl_gyrokinetic_app *app, int sidx
 struct gkyl_app_restart_status
 gkyl_gyrokinetic_app_read_from_frame(gkyl_gyrokinetic_app *app, int frame)
 {
-  struct gkyl_app_restart_status rstat;
+  struct gkyl_app_restart_status rstat = {
+    .io_status = GKYL_ARRAY_RIO_SUCCESS,
+    .frame = frame,
+    .stime = 0.0,
+  };
   for (int i=0; i<app->num_neut_species; i++) {
     if (app->neut_species[i].info.is_static) {
       gk_neut_species_apply_ic(app, &app->neut_species[i], 0.0);
     }
     else {
       rstat = gkyl_gyrokinetic_app_from_frame_neut_species(app, i, frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
   }
   for (int i=0; i<app->num_species; i++) {
@@ -3446,6 +3462,8 @@ gkyl_gyrokinetic_app_read_from_frame(gkyl_gyrokinetic_app *app, int frame)
     }
     else {
       rstat = gkyl_gyrokinetic_app_from_frame_species(app, i, frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
   }
 
@@ -3638,6 +3656,14 @@ gkyl_gyrokinetic_app_reset_species_positivity(gkyl_gyrokinetic_app* app, double
   gk_species_positivity_reset(app, tm, gks, &gks->positivity, pos_inp);
 }
 
+void
+gkyl_gyrokinetic_app_reset_species_damping(gkyl_gyrokinetic_app* app, double tm,
+  const char *species_name, struct gkyl_gyrokinetic_damping damping_inp)
+{
+  struct gk_species *gks = gk_find_species(app, species_name);
+  gk_species_damping_reset(app, tm, gks, &gks->damping, damping_inp);
+}
+
 void
 gkyl_gyrokinetic_app_reset_field(gkyl_gyrokinetic_app* app, double tm,
   struct gkyl_gyrokinetic_field field_inp)
diff --git a/gyrokinetic/apps/gyrokinetic_multib.c b/gyrokinetic/apps/gyrokinetic_multib.c
index de0352dee..99a1558bd 100644
--- a/gyrokinetic/apps/gyrokinetic_multib.c
+++ b/gyrokinetic/apps/gyrokinetic_multib.c
@@ -1067,7 +1067,11 @@ gkyl_gyrokinetic_multib_app_from_file_neut_species(gkyl_gyrokinetic_multib_app *
 struct gkyl_app_restart_status
 gkyl_gyrokinetic_multib_app_read_from_frame(gkyl_gyrokinetic_multib_app *app, int frame)
 {
-  struct gkyl_app_restart_status rstat;
+  struct gkyl_app_restart_status rstat = {
+    .io_status = GKYL_ARRAY_RIO_SUCCESS,
+    .frame = frame,
+    .stime = 0.0,
+  };
   for (int b=0; b<app->num_local_blocks; ++b) {
     struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
     for (int i=0; i<app->num_neut_species; i++) {
@@ -1076,9 +1080,13 @@ gkyl_gyrokinetic_multib_app_read_from_frame(gkyl_gyrokinetic_multib_app *app, in
         neut_frame = 0;
       }
       rstat = gkyl_gyrokinetic_app_from_frame_neut_species(sbapp, i, neut_frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
     for (int i=0; i<app->num_species; i++) {
       rstat = gkyl_gyrokinetic_app_from_frame_species(sbapp, i, frame);
+      if (rstat.io_status != GKYL_ARRAY_RIO_SUCCESS)
+        return rstat;
     }
   }
   
diff --git a/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c b/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c
index 8f079ac88..fa2e1a883 100644
--- a/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c
+++ b/gyrokinetic/apps/gyrokinetic_multib_update_ssp_rk3.c
@@ -2,7 +2,8 @@
 
 static void
 gyrokinetic_multib_forward_euler(struct gkyl_gyrokinetic_multib_app* app, double tcurr, double dt,
-  const struct gkyl_array *fin[], struct gkyl_array *fout[], 
+  const struct gkyl_array *fin[], struct gkyl_array *fout[],
+  const struct gkyl_array *fbar_in[], struct gkyl_array *fbar_out[],
   struct gkyl_array **bflux_in[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], 
   struct gkyl_array **bflux_in_neut[], struct gkyl_array **bflux_out_neut[], 
@@ -22,7 +23,7 @@ gyrokinetic_multib_forward_euler(struct gkyl_gyrokinetic_multib_app* app, double
   for (int b=0; b<app->num_local_blocks; ++b) {
     int li_charged = b * app->num_species;
     int li_neut = b * app->num_neut_species;
-    gyrokinetic_rhs(app->singleb_apps[b], tcurr, dt, &fin[li_charged], &fout[li_charged],
+    gyrokinetic_rhs(app->singleb_apps[b], tcurr, dt, &fin[li_charged], &fbar_in[li_charged], &fout[li_charged],
       &bflux_out[li_charged], &fin_neut[li_neut], &fout_neut[li_neut], &bflux_out_neut[li_neut], st);
     dtmin = fmin(dtmin, st->dt_actual);
   }
@@ -44,6 +45,7 @@ gyrokinetic_multib_forward_euler(struct gkyl_gyrokinetic_multib_app* app, double
     for (int i=0; i<app->num_species; ++i) {
       struct gk_species *gks = &sbapp->species[i];
       gk_species_step_f(gks, fout[li_charged+i], dta, fin[li_charged+i]);
+      gk_species_damping_forward_euler(gks, fin[li_charged+i], fbar_in[li_charged+i], fbar_out[li_charged+i], dta);
       gk_species_bflux_step_f(sbapp, &gks->bflux, bflux_out[li_charged+i], 1.0, bflux_in[li_charged+i]);
     }
     for (int i=0; i<app->num_neut_species; ++i) {
@@ -69,6 +71,8 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
 
   const struct gkyl_array *fin[ns_charged * nblocks_local];
   struct gkyl_array *fout[ns_charged * nblocks_local];
+  const struct gkyl_array *fbar_in[ns_charged * nblocks_local];
+  struct gkyl_array *fbar_out[ns_charged * nblocks_local];
   struct gkyl_array **bflux_in[ns_charged * nblocks_local];
   struct gkyl_array **bflux_out[ns_charged * nblocks_local];
 
@@ -94,6 +98,8 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
             struct gk_species *gks = &sbapp->species[i];
             fin[li_charged+i] = gks->f;
             fout[li_charged+i] = gks->f1;
+            fbar_in[li_charged+i] = gks->damping.fbar;
+            fbar_out[li_charged+i] = gks->damping.fbar1;
             // Boundary fluxes.
             bflux_in[li_charged+i] = gks->bflux.f;
             bflux_out[li_charged+i] = gks->bflux.f1;
@@ -101,14 +107,14 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           for (int i=0; i<ns_neut; ++i) {
             struct gk_neut_species *gkns = &sbapp->neut_species[i];
             fin_neut[li_neut+i] = gkns->f;
-	    fout_neut[li_neut+i] = gkns->f1;
+            fout_neut[li_neut+i] = gkns->f1;
             // Boundary fluxes.
             bflux_in_neut[li_neut+i] = gkns->bflux.f;
             bflux_out_neut[li_neut+i] = gkns->bflux.f1;
           }
         }
 
-        gyrokinetic_multib_forward_euler(app, tcurr, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_multib_forward_euler(app, tcurr, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
         dt = st.dt_actual;
 
@@ -121,11 +127,11 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           int li_charged = b * ns_charged;
           int li_neut = b * ns_neut;
           for (int i=0; i<ns_charged; ++i) {
-	    struct gk_species *gks = &sbapp->species[i]; 
+            struct gk_species *gks = &sbapp->species[i]; 
             gk_species_bflux_accumulate(sbapp, &gks->bflux, bflux_out[li_charged+i], -1.0, bflux_in[li_charged+i]);
-	  }
+          }
           for (int i=0; i<ns_neut; ++i) {
-	    struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
+            struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
             gk_neut_species_bflux_accumulate(sbapp, &gkns->bflux, bflux_out_neut[li_neut+i], -1.0, bflux_in_neut[li_neut+i]);
           }
         }
@@ -159,21 +165,23 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
             struct gk_species *gks = &sbapp->species[i];
             fin[li_charged+i] = gks->f1;
             fout[li_charged+i] = gks->fnew;
+            fbar_in[li_charged+i] = gks->damping.fbar1;
+            fbar_out[li_charged+i] = gks->damping.fbarnew;
             // Boundary fluxes.
             bflux_in[li_charged+i] = gks->bflux.f1;
             bflux_out[li_charged+i] = gks->bflux.fnew;
           }
           for (int i=0; i<ns_neut; ++i) {
             struct gk_neut_species *gkns = &sbapp->neut_species[i];
-	    fin_neut[li_neut+i] = gkns->f1;
-	    fout_neut[li_neut+i] = gkns->fnew;
+            fin_neut[li_neut+i] = gkns->f1;
+            fout_neut[li_neut+i] = gkns->fnew;
             // Boundary fluxes.
             bflux_in_neut[li_neut+i] = gkns->bflux.f1;
             bflux_out_neut[li_neut+i] = gkns->bflux.fnew;
           }
         }
 
-        gyrokinetic_multib_forward_euler(app, tcurr+dt, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_multib_forward_euler(app, tcurr+dt, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
         if (st.dt_actual < dt) {
@@ -207,13 +215,15 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           for (int b=0; b<nblocks_local; ++b) {
             struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
             for (int i=0; i<ns_charged; ++i) {
-	      struct gk_species *gks = &sbapp->species[i]; 
-	      gk_species_combine(gks, gks->f1, 3.0/4.0, gks->f, 1.0/4.0, gks->fnew, &gks->local_ext);
+              struct gk_species *gks = &sbapp->species[i]; 
+              gk_species_combine(gks, gks->f1, 3.0/4.0, gks->f, 1.0/4.0, gks->fnew, &gks->local_ext);
+              gk_species_damping_combine(gks, gks->damping.fbar1, 3.0/4.0, gks->damping.fbar,
+                1.0/4.0, gks->damping.fbarnew, &gks->local_ext);
               gk_species_bflux_set(sbapp, &gks->bflux, gks->bflux.f1, 1.0/4.0, gks->bflux.fnew);
-	    }
+            }
             for (int i=0; i<ns_neut; ++i) {
-	      struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
-	      gk_neut_species_combine(gkns, gkns->f1, 3.0/4.0, gkns->f, 1.0/4.0, gkns->fnew, &gkns->local_ext);
+              struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
+              gk_neut_species_combine(gkns, gkns->f1, 3.0/4.0, gkns->f, 1.0/4.0, gkns->fnew, &gkns->local_ext);
               gk_neut_species_bflux_set(sbapp, &gkns->bflux, gkns->bflux.f1, 1.0/4.0, gkns->bflux.fnew);
             }
           }
@@ -245,24 +255,26 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           int li_charged = b * ns_charged;
           int li_neut = b * ns_neut;
           for (int i=0; i<ns_charged; ++i) {
-	    struct gk_species *gks = &sbapp->species[i]; 
+            struct gk_species *gks = &sbapp->species[i]; 
             fin[li_charged+i] = gks->f1;
             fout[li_charged+i] = gks->fnew;
+            fbar_in[li_charged+i] = gks->damping.fbar1;
+            fbar_out[li_charged+i] = gks->damping.fbarnew;
             // Boundary fluxes.
             bflux_in[li_charged+i] = gks->bflux.f1;
             bflux_out[li_charged+i] = gks->bflux.fnew;
           }
           for (int i=0; i<ns_neut; ++i) {
             struct gk_neut_species *gkns = &sbapp->neut_species[i];
-	    fin_neut[li_neut+i] = sbapp->neut_species[i].f1;
-	    fout_neut[li_neut+i] = sbapp->neut_species[i].fnew;
+            fin_neut[li_neut+i] = sbapp->neut_species[i].f1;
+            fout_neut[li_neut+i] = sbapp->neut_species[i].fnew;
             // Boundary fluxes.
             bflux_in_neut[li_neut+i] = gkns->bflux.f1;
             bflux_out_neut[li_neut+i] = gkns->bflux.fnew;
           }
         }
 
-        gyrokinetic_multib_forward_euler(app, tcurr+dt/2, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_multib_forward_euler(app, tcurr+dt/2, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
         if (st.dt_actual < dt) {
@@ -296,18 +308,22 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
           for (int b=0; b<nblocks_local; ++b) {
             struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
             for (int i=0; i<ns_charged; ++i) {
-	      struct gk_species *gks = &sbapp->species[i]; 
+              struct gk_species *gks = &sbapp->species[i]; 
               // Step f.
-	      gk_species_combine(gks, gks->f1, 1.0/3.0, gks->f, 2.0/3.0, gks->fnew, &gks->local_ext);
-	      gk_species_copy_range(gks, gks->f, gks->f1, &gks->local_ext);
+              gk_species_combine(gks, gks->f1, 1.0/3.0, gks->f, 2.0/3.0, gks->fnew, &gks->local_ext);
+              gk_species_copy_range(gks, gks->f, gks->f1, &gks->local_ext);
+              gk_species_damping_combine(gks, gks->damping.fbar, 1.0/3.0, gks->damping.fbar1,
+                2.0/3.0, gks->damping.fbarnew, &gks->local_ext);
+              gk_species_damping_copy_range(gks, gks->damping.fbar1, gks->damping.fbar,
+                &gks->local_ext);
               // Step boundary fluxes.
               gk_species_bflux_set(sbapp, &gks->bflux, gks->bflux.f, 2.0/3.0, gks->bflux.fnew);
               gk_species_bflux_calc_voltime_integrated_mom(sbapp, gks, &gks->bflux, tcurr);
-	    }
+            }
             for (int i=0; i<ns_neut; ++i) {
-	      struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
-	      gk_neut_species_combine(gkns, gkns->f1, 1.0/3.0, gkns->f, 2.0/3.0, gkns->fnew, &gkns->local_ext);
-	      gk_neut_species_copy_range(gkns, gkns->f, gkns->f1, &gkns->local_ext);
+              struct gk_neut_species *gkns = &sbapp->neut_species[i]; 
+              gk_neut_species_combine(gkns, gkns->f1, 1.0/3.0, gkns->f, 2.0/3.0, gkns->fnew, &gkns->local_ext);
+              gk_neut_species_copy_range(gkns, gkns->f, gkns->f1, &gkns->local_ext);
               // Step boundary fluxes.
               gk_neut_species_bflux_set(sbapp, &gkns->bflux, gkns->bflux.f, 2.0/3.0, gkns->bflux.fnew);
               gk_neut_species_bflux_calc_voltime_integrated_mom(sbapp, gkns, &gkns->bflux, tcurr);
@@ -350,7 +366,7 @@ gyrokinetic_multib_update_ssp_rk3(struct gkyl_gyrokinetic_multib_app* app, doubl
 
           for (int b=0; b<nblocks_local; ++b) {
             struct gkyl_gyrokinetic_app *sbapp = app->singleb_apps[b];
-	    for (int i=0; i<ns_charged; ++i) {
+            for (int i=0; i<ns_charged; ++i) {
               struct gk_species *gks = &sbapp->species[i];
               // Compute moment of f_new to compute moment of df/dt.
               // Need to do it after the fields are updated.
diff --git a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
index 18a1a7e7a..d7f2fce37 100644
--- a/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
+++ b/gyrokinetic/apps/gyrokinetic_update_ssp_rk3.c
@@ -3,6 +3,7 @@
 static void
 gyrokinetic_forward_euler(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   const struct gkyl_array *fin[], struct gkyl_array *fout[], 
+  const struct gkyl_array *fbar_in[], struct gkyl_array *fbar_out[],
   struct gkyl_array **bflux_in[], struct gkyl_array **bflux_out[], 
   const struct gkyl_array *fin_neut[], struct gkyl_array *fout_neut[], 
   struct gkyl_array **bflux_in_neut[], struct gkyl_array **bflux_out_neut[], 
@@ -18,7 +19,8 @@ gyrokinetic_forward_euler(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   app->stat.nfeuler += 1;
 
   // Compute the time rate of change of the distributions, df/dt.
-  gyrokinetic_rhs(app, tcurr, dt, fin, fout, bflux_out, fin_neut, fout_neut, bflux_out_neut, st);
+  gyrokinetic_rhs(app, tcurr, dt, fin, fbar_in, fout, bflux_out, fin_neut, fout_neut,
+    bflux_out_neut, st);
 
   struct timespec wst = gkyl_wall_clock();
   // Complete update of distribution functions.
@@ -26,6 +28,7 @@ gyrokinetic_forward_euler(gkyl_gyrokinetic_app* app, double tcurr, double dt,
   for (int i=0; i<app->num_species; ++i) {
     struct gk_species *gks = &app->species[i];
     gk_species_step_f(gks, fout[i], dta, fin[i]);
+    gk_species_damping_forward_euler(gks, fin[i], fbar_in[i], fbar_out[i], dta);
     gk_species_bflux_accumulate(app, &gks->bflux, bflux_out[i], 1.0, bflux_in[i]);
   }
   for (int i=0; i<app->num_neut_species; ++i) {
@@ -46,6 +49,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
   // from the actual time-step.
   const struct gkyl_array *fin[app->num_species];
   struct gkyl_array *fout[app->num_species];
+  const struct gkyl_array *fbar_in[app->num_species];
+  struct gkyl_array *fbar_out[app->num_species];
   struct gkyl_array **bflux_in[app->num_species];
   struct gkyl_array **bflux_out[app->num_species];
 
@@ -67,6 +72,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           struct gk_species *gks = &app->species[i];
           fin[i] = gks->f;
           fout[i] = gks->f1;
+          fbar_in[i] = gks->damping.fbar;
+          fbar_out[i] = gks->damping.fbar1;
           // Boundary fluxes.
           bflux_in[i] = gks->bflux.f;
           bflux_out[i] = gks->bflux.f1;
@@ -86,7 +93,7 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           gk_species_source_adapt(app, gks, &gks->src, gks->lte.f_lte, bflux_in, tcurr);
         }
 
-        gyrokinetic_forward_euler(app, tcurr, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_forward_euler(app, tcurr, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
         dt = st.dt_actual;
 
@@ -125,6 +132,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           struct gk_species *gks = &app->species[i];
           fin[i] = gks->f1;
           fout[i] = gks->fnew;
+          fbar_in[i] = gks->damping.fbar1;
+          fbar_out[i] = gks->damping.fbarnew;
           // Boundary fluxes.
           bflux_in[i] = gks->bflux.f1;
           bflux_out[i] = gks->bflux.fnew;
@@ -138,7 +147,7 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           bflux_out_neut[i] = gkns->bflux.fnew;
         }
 
-        gyrokinetic_forward_euler(app, tcurr+dt, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_forward_euler(app, tcurr+dt, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
         if (st.dt_actual < dt) {
@@ -168,6 +177,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           for (int i=0; i<app->num_species; ++i) {
             struct gk_species *gks = &app->species[i];
             gk_species_combine(gks, gks->f1, 3.0/4.0, gks->f, 1.0/4.0, gks->fnew, &gks->local_ext);
+            gk_species_damping_combine(gks, gks->damping.fbar1, 3.0/4.0, gks->damping.fbar,
+              1.0/4.0, gks->damping.fbarnew, &gks->local_ext);
             gk_species_bflux_set(app, &gks->bflux, gks->bflux.f1, 1.0/4.0, gks->bflux.fnew);
           }
           for (int i=0; i<app->num_neut_species; ++i) {
@@ -197,6 +208,8 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           struct gk_species *gks = &app->species[i];
           fin[i] = gks->f1;
           fout[i] = gks->fnew;
+          fbar_in[i] = gks->damping.fbar1;
+          fbar_out[i] = gks->damping.fbarnew;
           // Boundary fluxes.
           bflux_in[i] = gks->bflux.f1;
           bflux_out[i] = gks->bflux.fnew;
@@ -210,7 +223,7 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
           bflux_out_neut[i] = gkns->bflux.fnew;
         }
 
-        gyrokinetic_forward_euler(app, tcurr+dt/2, dt, fin, fout, bflux_in, bflux_out,
+        gyrokinetic_forward_euler(app, tcurr+dt/2, dt, fin, fout, fbar_in, fbar_out, bflux_in, bflux_out,
           fin_neut, fout_neut, bflux_in_neut, bflux_out_neut, &st);
 
         if (st.dt_actual < dt) {
@@ -242,6 +255,10 @@ gyrokinetic_update_ssp_rk3(gkyl_gyrokinetic_app* app, double dt0)
             // Step f.
             gk_species_combine(gks, gks->f1, 1.0/3.0, gks->f, 2.0/3.0, gks->fnew, &gks->local_ext);
             gk_species_copy_range(gks, gks->f, gks->f1, &gks->local_ext);
+            gk_species_damping_combine(gks, gks->damping.fbar, 1.0/3.0, gks->damping.fbar1,
+              2.0/3.0, gks->damping.fbarnew, &gks->local_ext);
+            gk_species_damping_copy_range(gks, gks->damping.fbar1, gks->damping.fbar,
+              &gks->local_ext);
             // Step boundary fluxes.
             gk_species_bflux_set(app, &gks->bflux, gks->bflux.f, 2.0/3.0, gks->bflux.fnew);
             gk_species_bflux_calc_voltime_integrated_mom(app, gks, &gks->bflux, tcurr);
diff --git a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
index d6f94e132..91128d814 100644
--- a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_damped_1x2v_p1.c
@@ -704,12 +704,10 @@ int main(int argc, char **argv)
     },
 
     .damping = {
-//      .type = GKYL_GK_DAMPING_USER_INPUT,
-//      .rate_profile = loss_cone_damping_rate_profile,
-//      .rate_profile_ctx = &ctx,
-      .type = GKYL_GK_DAMPING_LOSS_CONE,
+      .type = GKYL_GK_DAMPING_USER_INPUT,
       .rate_profile = loss_cone_damping_rate_scaling,
       .rate_profile_ctx = &ctx,
+      .write_rate = true,
     },
 
     .bcs = {
diff --git a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
index 73f3e90e5..06f025ccb 100644
--- a/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_mirror_boltz_elc_poa_1x2v_p1.c
@@ -27,6 +27,7 @@ struct gk_poa_phase_params {
   bool is_static_field; // Whether to evolve the field.
   bool is_positivity_enabled; // Whether positivity is enabled.
   enum gkyl_gyrokinetic_fdot_multiplier_type fdot_mult_type; // Type of df/dt multipler.
+  enum gkyl_gyrokinetic_damping_type damping_type; // Type of damping.
 };
 
 // Define the context of the simulation. This is basically all the globals
@@ -576,6 +577,7 @@ create_ctx(void)
     poa_phases[2*i].is_static_field = is_static_field_oap;
     poa_phases[2*i].fdot_mult_type = fdot_mult_type_oap;
     poa_phases[2*i].is_positivity_enabled = is_positivity_enabled_oap;
+    poa_phases[2*i].damping_type = GKYL_GK_DAMPING_NONE;
 
     // FDPs.
     poa_phases[2*i+1].phase = GK_POA_FDP;
@@ -585,6 +587,7 @@ create_ctx(void)
     poa_phases[2*i+1].is_static_field = is_static_field_fdp;
     poa_phases[2*i+1].fdot_mult_type = fdot_mult_type_fdp;
     poa_phases[2*i+1].is_positivity_enabled = is_positivity_enabled_fdp;
+    poa_phases[2*i+1].damping_type = GKYL_GK_DAMPING_LOW_PASS_FILTER;
   }
   // Add an extra, longer FDP.
   poa_phases[num_phases-1].phase = GK_POA_FDP;
@@ -594,6 +597,7 @@ create_ctx(void)
   poa_phases[num_phases-1].is_static_field = is_static_field_fdp;
   poa_phases[num_phases-1].fdot_mult_type = fdot_mult_type_fdp;
   poa_phases[num_phases-1].is_positivity_enabled = is_positivity_enabled_fdp;
+  poa_phases[num_phases-1].damping_type = GKYL_GK_DAMPING_LOW_PASS_FILTER;
 
   double write_phase_freq = 0.5; // Frequency of writing phase-space diagnostics (as a fraction of num_frames).
   double int_diag_calc_freq = 5; // Frequency of calculating integrated diagnostics (as a factor of num_frames).
@@ -780,10 +784,19 @@ void run_phase(gkyl_gyrokinetic_app* app, struct gk_mirror_ctx *ctx, double num_
     .type = pparams->is_positivity_enabled? GKYL_GK_POSITIVITY_SHIFT : GKYL_GK_POSITIVITY_NONE,
     .write_diagnostics = pparams->is_positivity_enabled,
   };
+  struct gkyl_gyrokinetic_damping damping_inp = {
+    .type = pparams->damping_type,
+    .cellwise_const = false,
+    .rate_const = 1/5e-6,
+    .write_fbar = true,
+    .write_rate = true,
+    .do_not_reset_fbar = true,
+  };
 
   gkyl_gyrokinetic_app_reset_species_collisionless(app, t_curr, "ion", collisionless_inp);
   gkyl_gyrokinetic_app_reset_species_fdot_multiplier(app, t_curr, "ion", fdot_mult_inp);
   gkyl_gyrokinetic_app_reset_species_positivity(app, t_curr, "ion", positivity_inp);
+  gkyl_gyrokinetic_app_reset_species_damping(app, t_curr, "ion", damping_inp);
   gkyl_gyrokinetic_app_reset_field(app, t_curr, field_inp);
 
   // Compute initial guess of maximum stable time-step.
@@ -1044,8 +1057,8 @@ int main(int argc, char **argv)
     write_data(&trig_write_conf, &trig_write_phase, app, tfs.t_curr, true);
   }
 
-  if (app_args.num_steps != INT_MAX)
-    phase_idx_end = 1;
+  // if (app_args.num_steps != INT_MAX)
+  //   phase_idx_end = 1;
 
   // Loop over number of number of phases;
   for (int pit=phase_idx_init; pit<phase_idx_end; pit++) {
diff --git a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
index ac4da0731..d9910f78e 100644
--- a/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
+++ b/gyrokinetic/creg/rt_gk_wham_1x2v_p1.c
@@ -618,6 +618,14 @@ int main(int argc, char **argv)
       .write_diagnostics = true,
     },
 
+    .damping = {
+      .type = GKYL_GK_DAMPING_LOW_PASS_FILTER,
+      .rate_const = 1e-6,
+      .write_rate = true,
+      .write_fbar = true,
+      .cellwise_const = true,
+    },
+    
     .source = {
       .source_id = GKYL_PROJ_SOURCE,
       .num_sources = 1,