diff --git a/docs/sphinx/components.rst b/docs/sphinx/components.rst
index 59efd819..98b839e2 100644
--- a/docs/sphinx/components.rst
+++ b/docs/sphinx/components.rst
@@ -2214,6 +2214,11 @@ electromagnetic
    (``electromagnetic:magnetic_flutter=true``).  They use an
    ``Apar_flutter`` field, not the ``Apar`` field that is calculated
    from the induction terms.
+5. If using ``vorticity:phi_boundary_relax`` to evolve the radial
+   boundary of the electrostatic potential, the timescale
+   ``phi_boundary_timescale`` should be set much longer than the
+   Alfven wave period or unphysical instabilities may grow from the
+   boundaries.
 
 This component modifies the definition of momentum of all species, to
 include the contribution from the electromagnetic potential
@@ -2289,5 +2294,12 @@ parallel derivative terms in the ``evolve_density``, ``evolve_pressure``, ``evol
 ``evolve_momentum`` components. Parallel flow terms are modified, and parallel heat
 conduction.
 
+.. math::
+
+   \begin{aligned}\mathbf{b}\cdot\nabla f =& \mathbf{b}_0\cdot\nabla f + \delta\mathbf{b}\cdot\nabla f \\
+   =& \mathbf{b}_0\cdot\nabla f + \frac{1}{B}\nabla\times\left(\mathbf{b}A_{||}\right)\cdot\nabla f \\
+   \simeq& \mathbf{b}_0\cdot\nabla f + \frac{1}{B_0}\left[A_{||}\nabla\times\mathbf{b} + \left(\nabla A_{||}\right)\times\mathbf{b}_0\right]\cdot\nabla f \\
+   \simeq& \mathbf{b}_0\cdot\nabla f + \frac{1}{B_0}\mathbf{b}_0\times \nabla A_{||} \cdot \nabla f\end{aligned}
+
 .. doxygenstruct:: Electromagnetic
    :members:
diff --git a/src/evolve_momentum.cxx b/src/evolve_momentum.cxx
index 46c768af..07a5a3ad 100644
--- a/src/evolve_momentum.cxx
+++ b/src/evolve_momentum.cxx
@@ -150,6 +150,14 @@ void EvolveMomentum::finally(const Options &state) {
         }
         ddt(NV) += Z * Apar * dndt;
       }
+      if (state["fields"].isSet("Apar_flutter")) {
+        // Magnetic flutter term
+        const Field3D Apar_flutter = get<Field3D>(state["fields"]["Apar_flutter"]);
+
+        // Using the approximation for small delta-B/B
+        // b dot Grad(phi) = Grad_par(phi) + [phi, Apar]
+        ddt(NV) -= Z * N * bracket(phi, Apar_flutter, BRACKET_ARAKAWA);
+      }
     } else {
       ddt(NV) = 0.0;
     }
diff --git a/src/evolve_pressure.cxx b/src/evolve_pressure.cxx
index 3f4d131d..d888827e 100644
--- a/src/evolve_pressure.cxx
+++ b/src/evolve_pressure.cxx
@@ -383,6 +383,8 @@ void EvolvePressure::finally(const Options& state) {
       Field3D db_dot_T = bracket(T, Apar_flutter, BRACKET_ARAKAWA);
       Field3D b0_dot_T = Grad_par(T);
       mesh->communicate(db_dot_T, b0_dot_T);
+      db_dot_T.applyBoundary("neumann");
+      b0_dot_T.applyBoundary("neumann");
       ddt(P) += (2. / 3) * (Div_par(kappa_par * db_dot_T) -
                             Div_n_g_bxGrad_f_B_XZ(kappa_par, db_dot_T + b0_dot_T, Apar_flutter));
     }
diff --git a/src/vorticity.cxx b/src/vorticity.cxx
index 496f7296..4ff3ab6a 100644
--- a/src/vorticity.cxx
+++ b/src/vorticity.cxx
@@ -634,6 +634,7 @@ void Vorticity::transform(Options& state) {
 
 void Vorticity::finally(const Options& state) {
   AUTO_TRACE();
+  auto coord = mesh->getCoordinates();
 
   phi = get<Field3D>(state["fields"]["phi"]);
 
@@ -670,7 +671,7 @@ void Vorticity::finally(const Options& state) {
     }
   }
 
-  if (state.isSection("fields") and state["fields"].isSet("DivJextra")) {
+  if (state["fields"].isSet("DivJextra")) {
     auto DivJextra = get<Field3D>(state["fields"]["DivJextra"]);
 
     // Parallel current is handled here, to allow different 2D or 3D closures
@@ -695,7 +696,18 @@ void Vorticity::finally(const Options& state) {
     const BoutReal A = get<BoutReal>(species["AA"]);
 
     // Note: Using NV rather than N*V so that the cell boundary flux is correct
-    ddt(Vort) += Div_par((Z / A) * NV);
+    const Field3D jpar = (Z / A) * NV;
+    ddt(Vort) += Div_par(jpar);
+
+    if (state["fields"].isSet("Apar_flutter")) {
+      // Magnetic flutter term
+      const Field3D Apar_flutter = get<Field3D>(state["fields"]["Apar_flutter"]);
+
+      // Div_par(jpar) = B * Grad_par(jpar / B)
+      // Using the approximation for small delta-B/B
+      // b dot Grad(jpar) = Grad_par(jpar) + [jpar, Apar]
+      ddt(Vort) += coord->Bxy * bracket(jpar / coord->Bxy, Apar_flutter, BRACKET_ARAKAWA);
+    }
   }
 
   // Viscosity