PSATD: Implement First-Order Equations (#3466)

* Implement First-Order PSATD Equations

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Fix Unused Parameter Warning

* Fix RZ Build

* Fix Normalization of G to Match PML

* Add CI Test: 3D Uniform Plasma

* Cleaning

* Update 2D CI Checksums

* Update 3D CI Checksums

* Add F,G to CI Checksums of `uniform_plasma_multiJ`

* Allow User to Choose First-Order v. Second-Order

* Add WARPX_ALWAYS_ASSERT_WITH_MESSAGE

* Rename New Class `PsatdAlgorithmFirstOrder`

* Remove Inline Comment

* Update RZ CI Checksums

* Fix inline comment

* Use auxiliary variables to avoid divisions

* Use auxiliary variables to avoid divisions

* Make `nci_psatd_stability` dir and merge inputs

* Move all Galilean tests to `nci_psatd_stability`

* Fix CI

* Fix index for backward FFT of J

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

13 files changed, 587 insertions, 70 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
index a370d4b2d..912ed47c4 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
@@ -1,5 +1,6 @@
 target_sources(WarpX
   PRIVATE
+    PsatdAlgorithmFirstOrder.cpp
     PsatdAlgorithmJConstantInTime.cpp
     PsatdAlgorithmJLinearInTime.cpp
     PsatdAlgorithmPml.cpp
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index c798ffb01..40f9d0e9a 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -1,4 +1,5 @@
 CEXE_sources += SpectralBaseAlgorithm.cpp
+CEXE_sources += PsatdAlgorithmFirstOrder.cpp
 CEXE_sources += PsatdAlgorithmJConstantInTime.cpp
 CEXE_sources += PsatdAlgorithmJLinearInTime.cpp
 CEXE_sources += PsatdAlgorithmPml.cpp
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp
index 1d6248f8d..dfd7ffe92 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp
@@ -103,9 +103,9 @@ PsatdAlgorithmComoving::pushSpectralFields (SpectralFieldData& f) const
             const Complex Bz_old = fields(i,j,k,Idx.Bz);
 
             // Shortcuts for the values of J and rho
-            const Complex Jx = fields(i,j,k,Idx.Jx);
-            const Complex Jy = fields(i,j,k,Idx.Jy);
-            const Complex Jz = fields(i,j,k,Idx.Jz);
+            const Complex Jx = fields(i,j,k,Idx.Jx_mid);
+            const Complex Jy = fields(i,j,k,Idx.Jy_mid);
+            const Complex Jz = fields(i,j,k,Idx.Jz_mid);
             const Complex rho_old = fields(i,j,k,Idx.rho_old);
             const Complex rho_new = fields(i,j,k,Idx.rho_new);
 
@@ -447,9 +447,9 @@ void PsatdAlgorithmComoving::CurrentCorrection (SpectralFieldData& field_data)
         amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             // Shortcuts for the values of J and rho
-            const Complex Jx = fields(i,j,k,Idx.Jx);
-            const Complex Jy = fields(i,j,k,Idx.Jy);
-            const Complex Jz = fields(i,j,k,Idx.Jz);
+            const Complex Jx = fields(i,j,k,Idx.Jx_mid);
+            const Complex Jy = fields(i,j,k,Idx.Jy_mid);
+            const Complex Jz = fields(i,j,k,Idx.Jz_mid);
             const Complex rho_old = fields(i,j,k,Idx.rho_old);
             const Complex rho_new = fields(i,j,k,Idx.rho_new);
 
@@ -482,15 +482,15 @@ void PsatdAlgorithmComoving::CurrentCorrection (SpectralFieldData& field_data)
                     const Complex theta = amrex::exp(- I * k_dot_v * dt * 0.5_rt);
                     const Complex den = 1._rt - theta * theta;
 
-                    fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod);
-                    fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod);
-                    fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx.Jx_mid) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx.Jy_mid) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx.Jz_mid) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod);
 
                 } else {
 
-                    fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod);
-                    fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod);
-                    fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx.Jx_mid) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx.Jy_mid) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod);
+                    fields(i,j,k,Idx.Jz_mid) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod);
                 }
             }
         });
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H
new file mode 100644
index 000000000..c90b19e1a
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H
@@ -0,0 +1,100 @@
+/* Copyright 2019
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PSATD_ALGORITHM_FIRST_ORDER_H_
+#define WARPX_PSATD_ALGORITHM_FIRST_ORDER_H_
+
+#include "FieldSolver/SpectralSolver/SpectralFieldData.H"
+#include "FieldSolver/SpectralSolver/SpectralKSpace.H"
+#include "SpectralBaseAlgorithm.H"
+
+#include <AMReX_Array.H>
+#include <AMReX_Config.H>
+#include <AMReX_REAL.H>
+
+#include <AMReX_BaseFwd.H>
+
+#include <array>
+#include <memory>
+
+#if WARPX_USE_PSATD
+/*
+ * \brief Class that updates the fields in spectral space according to the first-order PSATD equations.
+ */
+class PsatdAlgorithmFirstOrder : public SpectralBaseAlgorithm
+{
+    public:
+
+        /**
+         * \brief Constructor of the class PsatdAlgorithmFirstOrder
+         *
+         * \param[in] spectral_kspace spectral space
+         * \param[in] dm distribution mapping
+         * \param[in] spectral_index object containing indices to access data in spectral space
+         * \param[in] norder_x order of the spectral solver along x
+         * \param[in] norder_y order of the spectral solver along y
+         * \param[in] norder_z order of the spectral solver along z
+         * \param[in] nodal whether the E and B fields are defined on a fully nodal grid or a Yee grid
+         * \param[in] dt time step of the simulation
+         * \param[in] div_cleaning whether to use divergence correction for both E and B (thus, F and G)
+         * \param[in] J_in_time time dependency of J (currently supported: constant, linear)
+         * \param[in] rho_in_time time dependency of rho (currently supported: constant, linear)
+         */
+        PsatdAlgorithmFirstOrder (
+            const SpectralKSpace& spectral_kspace,
+            const amrex::DistributionMapping& dm,
+            const SpectralFieldIndex& spectral_index,
+            const int norder_x,
+            const int norder_y,
+            const int norder_z,
+            const bool nodal,
+            const amrex::Real dt,
+            const bool div_cleaning,
+            const int J_in_time,
+            const int rho_in_time);
+
+        /**
+         * \brief Updates E, B, F, and G fields in spectral space,
+         *        according to the first-order PSATD equations
+         *
+         * \param[in,out] f all the fields in spectral space
+         */
+        virtual void pushSpectralFields (SpectralFieldData& f) const override final;
+
+        /**
+         * \brief Virtual function for current correction in Fourier space
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
+         * This function overrides the virtual function \c CurrentCorrection in the
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         */
+        virtual void CurrentCorrection (SpectralFieldData& field_data) override final;
+
+        /**
+         * \brief Virtual function for Vay current deposition in Fourier space
+         * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>).
+         * This function overrides the virtual function \c VayDeposition in the
+         * base class \c SpectralBaseAlgorithm and cannot be overridden by further
+         * derived classes.
+         *
+         * \param[in,out] field_data All fields in Fourier space
+         */
+        virtual void VayDeposition (SpectralFieldData& field_data) override final;
+
+    private:
+
+        SpectralFieldIndex m_spectral_index;
+
+        // Other member variables
+        amrex::Real m_dt;
+        bool m_div_cleaning;
+        int m_J_in_time;
+        int m_rho_in_time;
+};
+#endif // WARPX_USE_PSATD
+#endif // WARPX_PSATD_ALGORITHM_FIRST_ORDER_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp
new file mode 100644
index 000000000..d32604760
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp
@@ -0,0 +1,375 @@
+/* Copyright 2019
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "PsatdAlgorithmFirstOrder.H"
+
+#include "Utils/TextMsg.H"
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "Utils/WarpXConst.H"
+#include "Utils/WarpX_Complex.H"
+
+#include <AMReX_Array4.H>
+#include <AMReX_BLProfiler.H>
+#include <AMReX_BaseFab.H>
+#include <AMReX_BoxArray.H>
+#include <AMReX_GpuComplex.H>
+#include <AMReX_GpuLaunch.H>
+#include <AMReX_GpuQualifiers.H>
+#include <AMReX_IntVect.H>
+#include <AMReX_MFIter.H>
+#include <AMReX_PODVector.H>
+
+#include <cmath>
+
+#if WARPX_USE_PSATD
+
+using namespace amrex::literals;
+
+PsatdAlgorithmFirstOrder::PsatdAlgorithmFirstOrder(
+    const SpectralKSpace& spectral_kspace,
+    const amrex::DistributionMapping& dm,
+    const SpectralFieldIndex& spectral_index,
+    const int norder_x,
+    const int norder_y,
+    const int norder_z,
+    const bool nodal,
+    const amrex::Real dt,
+    const bool div_cleaning,
+    const int J_in_time,
+    const int rho_in_time)
+    // Initializer list
+    : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal),
+    m_spectral_index(spectral_index),
+    m_dt(dt),
+    m_div_cleaning(div_cleaning),
+    m_J_in_time(J_in_time),
+    m_rho_in_time(rho_in_time)
+{}
+
+void
+PsatdAlgorithmFirstOrder::pushSpectralFields (SpectralFieldData& f) const
+{
+    const bool div_cleaning = m_div_cleaning;
+
+    const bool J_constant = (m_J_in_time == JInTime::Constant) ? true : false;
+    const bool J_linear   = (m_J_in_time == JInTime::Linear  ) ? true : false;
+    const bool rho_constant = (m_rho_in_time == RhoInTime::Constant) ? true : false;
+    const bool rho_linear   = (m_rho_in_time == RhoInTime::Linear  ) ? true : false;
+
+    const amrex::Real dt = m_dt;
+    const amrex::Real dt2 = dt*dt;
+
+    const SpectralFieldIndex& Idx = m_spectral_index;
+
+    // Loop over boxes
+    for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi)
+    {
+        const amrex::Box& bx = f.fields[mfi].box();
+
+        // Extract arrays for the fields to be updated
+        amrex::Array4<Complex> fields = f.fields[mfi].array();
+
+        // Extract pointers for the k vectors
+        const amrex::Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+#if defined(WARPX_DIM_3D)
+        const amrex::Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+#endif
+        const amrex::Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr();
+
+        // Loop over indices within one box
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            // Record old values of the fields to be updated
+            const Complex Ex_old = fields(i,j,k,Idx.Ex);
+            const Complex Ey_old = fields(i,j,k,Idx.Ey);
+            const Complex Ez_old = fields(i,j,k,Idx.Ez);
+            const Complex Bx_old = fields(i,j,k,Idx.Bx);
+            const Complex By_old = fields(i,j,k,Idx.By);
+            const Complex Bz_old = fields(i,j,k,Idx.Bz);
+
+            // Shortcuts for the values of J and rho
+            const Complex Jx_mid = (J_constant) ? fields(i,j,k,Idx.Jx_mid) : 0._rt;
+            const Complex Jy_mid = (J_constant) ? fields(i,j,k,Idx.Jy_mid) : 0._rt;
+            const Complex Jz_mid = (J_constant) ? fields(i,j,k,Idx.Jz_mid) : 0._rt;
+            const Complex Jx_old = (J_linear  ) ? fields(i,j,k,Idx.Jx_old) : 0._rt;
+            const Complex Jy_old = (J_linear  ) ? fields(i,j,k,Idx.Jy_old) : 0._rt;
+            const Complex Jz_old = (J_linear  ) ? fields(i,j,k,Idx.Jz_old) : 0._rt;
+            const Complex Jx_new = (J_linear  ) ? fields(i,j,k,Idx.Jx_new) : 0._rt;
+            const Complex Jy_new = (J_linear  ) ? fields(i,j,k,Idx.Jy_new) : 0._rt;
+            const Complex Jz_new = (J_linear  ) ? fields(i,j,k,Idx.Jz_new) : 0._rt;
+
+            const Complex Jx_c0 = (J_constant) ? Jx_mid : Jx_old;
+            const Complex Jy_c0 = (J_constant) ? Jy_mid : Jy_old;
+            const Complex Jz_c0 = (J_constant) ? Jz_mid : Jz_old;
+            const Complex Jx_c1 = (J_linear  ) ? (Jx_new-Jx_old)/dt : 0._rt;
+            const Complex Jy_c1 = (J_linear  ) ? (Jy_new-Jy_old)/dt : 0._rt;
+            const Complex Jz_c1 = (J_linear  ) ? (Jz_new-Jz_old)/dt : 0._rt;
+
+            Complex rho_mid, rho_old, rho_new, F_old, G_old;
+            Complex rho_c0, rho_c1;
+            if (div_cleaning)
+            {
+                rho_mid = (rho_constant) ? fields(i,j,k,Idx.rho_mid) : 0._rt;
+                rho_old = (rho_linear  ) ? fields(i,j,k,Idx.rho_old) : 0._rt;
+                rho_new = (rho_linear  ) ? fields(i,j,k,Idx.rho_new) : 0._rt;
+
+                F_old = fields(i,j,k,Idx.F);
+                G_old = fields(i,j,k,Idx.G);
+
+                rho_c0 = (rho_constant) ? rho_mid : rho_old;
+                rho_c1 = (rho_linear  ) ? (rho_new-rho_old)/dt : 0._rt;
+            }
+
+            // k vector values
+            const amrex::Real kx = modified_kx_arr[i];
+#if defined(WARPX_DIM_3D)
+            const amrex::Real ky = modified_ky_arr[j];
+            const amrex::Real kz = modified_kz_arr[k];
+#else
+            constexpr amrex::Real ky = 0._rt;
+            const     amrex::Real kz = modified_kz_arr[j];
+#endif
+            // Physical constants and imaginary unit
+            constexpr amrex::Real c = PhysConst::c;
+            constexpr amrex::Real c2 = c*c;
+            constexpr amrex::Real inv_c = 1._rt/c;
+            constexpr amrex::Real mu0 = PhysConst::mu0;
+            constexpr Complex I = Complex{0._rt, 1._rt};
+
+            const amrex::Real kx2 = kx*kx;
+            const amrex::Real ky2 = ky*ky;
+            const amrex::Real kz2 = kz*kz;
+
+            const amrex::Real knorm = std::sqrt(kx2 + ky2 + kz2);
+            const amrex::Real knorm2 = knorm*knorm;
+            const amrex::Real knorm4 = knorm2*knorm2;
+
+            // Auxiliary variables
+            const amrex::Real inv_knorm  = 1._rt/knorm;
+            const amrex::Real inv_knorm2 = 1._rt/knorm2;
+            const amrex::Real inv_knorm4 = 1._rt/knorm4;
+
+            const amrex::Real C = std::cos(c*knorm*dt);
+            const amrex::Real S = std::sin(c*knorm*dt);
+
+            // Update equations
+
+            if (knorm == 0._rt)
+            {
+                fields(i,j,k,Idx.Ex) = Ex_old - mu0*c2*dt*Jx_c0 - 0.5_rt*mu0*c2*dt2*Jx_c1;
+                fields(i,j,k,Idx.Ey) = Ey_old - mu0*c2*dt*Jy_c0 - 0.5_rt*mu0*c2*dt2*Jy_c1;
+                fields(i,j,k,Idx.Ez) = Ez_old - mu0*c2*dt*Jz_c0 - 0.5_rt*mu0*c2*dt2*Jz_c1;
+
+                if (div_cleaning)
+                {
+                    fields(i,j,k,Idx.F) = F_old - mu0*c2*dt*rho_c0 - 0.5_rt*mu0*c2*dt2*rho_c1;
+                }
+            }
+            else // knorm != 0
+            {
+                Complex C01, C02, C03, C04, C05, C06, C07, C08,
+                        C09, C10, C11, C12, C13, C14, C15, C16;
+
+                // Ex
+                C01 = (div_cleaning) ? C : (kx2+ky2*C+kz2*C)*inv_knorm2;
+                C02 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2;
+                C03 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2;
+                C04 = 0._rt;
+                C05 = -I*c*kz*S*inv_knorm;
+                C06 =  I*c*ky*S*inv_knorm;
+                C07 = (div_cleaning) ? I*c*kx*S*inv_knorm : 0._rt;
+                C09 = (div_cleaning) ? -mu0*c*S*inv_knorm : -mu0*c*(dt*c*kx2*knorm2+ky2*knorm*S+kz2*knorm*S)*inv_knorm4;
+                C10 = (div_cleaning) ? 0._rt : mu0*c*kx*ky*(knorm*S-dt*c*knorm2)*inv_knorm4;
+                C11 = (div_cleaning) ? 0._rt : mu0*c*kx*kz*(knorm*S-dt*c*knorm2)*inv_knorm4;
+                C12 = 0._rt; // This is not redundant, do not remove this
+                if (J_linear) C12 = (div_cleaning) ? mu0*(C-1._rt)*inv_knorm2 : mu0*(2._rt*ky2*(C-1._rt)+2._rt*kz2*(C-1._rt)-dt2*c2*kx2*knorm2)*inv_knorm4*0.5_rt;
+                C13 = (J_linear && !div_cleaning) ? mu0*kx*ky*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt;
+                C14 = (J_linear && !div_cleaning) ? mu0*kx*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt;
+                C15 = (div_cleaning) ? I*mu0*c2*kx*(C-1._rt)*inv_knorm2 : 0._rt;
+                C16 = (div_cleaning && rho_linear) ? I*mu0*c*kx*(knorm*S-dt*c*knorm2)*inv_knorm4 : 0._rt;
+
+                fields(i,j,k,Idx.Ex) = C01*Ex_old + C02*Ey_old + C03*Ez_old
+                                     + C04*Bx_old + C05*By_old + C06*Bz_old
+                                     + C07*F_old // only with div cleaning
+                                     + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0
+                                     + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time
+                                     + C15*rho_c0  // only with div cleaning
+                                     + C16*rho_c1; // only with div cleaning and rho linear in time
+
+                // Ey
+                C01 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2;
+                C02 = (div_cleaning) ? C : (kx2*C+ky2+kz2*C)*inv_knorm2;
+                C03 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2;
+                C04 = I*c*kz*S*inv_knorm;
+                C05 = 0._rt;
+                C06 = -I*c*kx*S*inv_knorm;
+                C07 = (div_cleaning) ? I*c*ky*S*inv_knorm : 0._rt;
+                C09 = (div_cleaning) ? 0._rt : mu0*c*kx*ky*(knorm*S-dt*c*knorm2)*inv_knorm4;
+                C10 = (div_cleaning) ? -mu0*c*S*inv_knorm : -mu0*c*(dt*c*ky2*knorm2+kx2*knorm*S+kz2*knorm*S)*inv_knorm4;
+                C11 = (div_cleaning) ? 0._rt : mu0*c*ky*kz*(knorm*S-dt*c*knorm2)*inv_knorm4;
+                C12 = (J_linear && !div_cleaning) ? mu0*kx*ky*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt;
+                C13 = 0._rt; // This is not redundant, do not remove this
+                if (J_linear) C13 = (div_cleaning) ? mu0*(C-1._rt)*inv_knorm2 : mu0*(2._rt*kx2*(C-1._rt)+2._rt*kz2*(C-1._rt)-dt2*c2*ky2*knorm2)*inv_knorm4*0.5_rt;
+                C14 = (J_linear && !div_cleaning) ? mu0*ky*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt;
+                C15 = (div_cleaning) ? I*mu0*c2*ky*(C-1._rt)*inv_knorm2 : 0._rt;
+                C16 = (div_cleaning && rho_linear) ? I*mu0*c*ky*(knorm*S-dt*c*knorm2)*inv_knorm4 : 0._rt;
+
+                fields(i,j,k,Idx.Ey) = C01*Ex_old + C02*Ey_old + C03*Ez_old
+                                     + C04*Bx_old + C05*By_old + C06*Bz_old
+                                     + C07*F_old // only with div cleaning
+                                     + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0
+                                     + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time
+                                     + C15*rho_c0  // only with div cleaning
+                                     + C16*rho_c1; // only with div cleaning and rho linear in time
+
+                // Ez
+                C01 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2;
+                C02 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2;
+                C03 = (div_cleaning) ? C : (kx2*C+ky2*C+kz2)*inv_knorm2;
+                C04 = -I*c*ky*S*inv_knorm;
+                C05 = I*c*kx*S*inv_knorm;
+                C06 = 0._rt;
+                C07 = (div_cleaning) ? I*c*kz*S*inv_knorm : 0._rt;
+                C09 = (div_cleaning) ? 0._rt : mu0*c*kx*kz*(knorm*S-dt*c*knorm2)*inv_knorm4;
+                C10 = (div_cleaning) ? 0._rt : mu0*c*ky*kz*(knorm*S-dt*c*knorm2)*inv_knorm4;
+                C11 = (div_cleaning) ? -mu0*c*S*inv_knorm : -mu0*c*(dt*c*kz2*knorm2+kx2*knorm*S+ky2*knorm*S)*inv_knorm4;
+                C12 = (J_linear && !div_cleaning) ? mu0*kx*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt;
+                C13 = (J_linear && !div_cleaning) ? mu0*ky*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt;
+                C14 = 0._rt; // This is not redundant, do not remove this
+                if (J_linear) C14 = (div_cleaning) ? mu0*(C-1._rt)*inv_knorm2 : mu0*(2._rt*kx2*(C-1._rt)+2._rt*ky2*(C-1._rt)-dt2*c2*kz2*knorm2)*inv_knorm4*0.5_rt;
+                C15 = (div_cleaning) ? I*mu0*c2*kz*(C-1._rt)*inv_knorm2 : 0._rt;
+                C16 = (div_cleaning && rho_linear) ? I*mu0*c*kz*(knorm*S-dt*c*knorm2)*inv_knorm4 : 0._rt;
+
+                fields(i,j,k,Idx.Ez) = C01*Ex_old + C02*Ey_old + C03*Ez_old
+                                     + C04*Bx_old + C05*By_old + C06*Bz_old
+                                     + C07*F_old // only with div cleaning
+                                     + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0
+                                     + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time
+                                     + C15*rho_c0  // only with div cleaning
+                                     + C16*rho_c1; // only with div cleaning and rho linear in time
+
+                // Bx
+                C01 = 0._rt;
+                C02 = I*kz*S*inv_knorm*inv_c;
+                C03 = -I*ky*S*inv_knorm*inv_c;
+                C04 = (div_cleaning) ? C : (kx2+ky2*C+kz2*C)*inv_knorm2;
+                C05 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2;
+                C06 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2;
+                C08 = (div_cleaning) ? I*kx*S*inv_knorm*inv_c : 0._rt;
+                C09 = 0._rt;
+                C10 = I*mu0*kz*(C-1._rt)*inv_knorm2;
+                C11 = -I*mu0*ky*(C-1._rt)*inv_knorm2;
+                C12 = 0._rt;
+                C13 = (J_linear) ? I*mu0*kz*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+                C14 = (J_linear) ? -I*mu0*ky*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+
+                fields(i,j,k,Idx.Bx) = C01*Ex_old + C02*Ey_old + C03*Ez_old
+                                     + C04*Bx_old + C05*By_old + C06*Bz_old
+                                     + C08*G_old // only with div cleaning
+                                     + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0
+                                     + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1; // only with J linear in time
+
+                // By
+                C01 = -I*kz*S*inv_knorm*inv_c;
+                C02 = 0._rt;
+                C03 = I*kx*S*inv_knorm*inv_c;
+                C04 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2;
+                C05 = (div_cleaning) ? C : (kx2*C+ky2+kz2*C)*inv_knorm2;
+                C06 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2;
+                C08 = (div_cleaning) ? I*ky*S*inv_knorm*inv_c : 0._rt;
+                C09 = -I*mu0*kz*(C-1._rt)*inv_knorm2;
+                C10 = 0._rt;
+                C11 = I*mu0*kx*(C-1._rt)*inv_knorm2;
+                C12 = (J_linear) ? -I*mu0*kz*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+                C13 = 0._rt;
+                C14 = (J_linear) ? I*mu0*kx*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+
+                fields(i,j,k,Idx.By) = C01*Ex_old + C02*Ey_old + C03*Ez_old
+                                     + C04*Bx_old + C05*By_old + C06*Bz_old
+                                     + C08*G_old // only with div cleaning
+                                     + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0
+                                     + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1; // only with J linear in time
+
+                // Bz
+                C01 = I*ky*S*inv_knorm*inv_c;
+                C02 = -I*kx*S*inv_knorm*inv_c;
+                C03 = 0._rt;
+                C04 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2;
+                C05 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2;
+                C06 = (div_cleaning) ? C : (kx2*C+ky2*C+kz2)*inv_knorm2;
+                C08 = (div_cleaning) ? I*kz*S*inv_knorm*inv_c : 0._rt;
+                C09 = I*mu0*ky*(C-1._rt)*inv_knorm2;
+                C10 = -I*mu0*kx*(C-1._rt)*inv_knorm2;
+                C11 = 0._rt;
+                C12 = (J_linear) ? I*mu0*ky*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+                C13 = (J_linear) ? -I*mu0*kx*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+                C14 = 0._rt;
+
+                fields(i,j,k,Idx.Bz) = C01*Ex_old + C02*Ey_old + C03*Ez_old
+                                     + C04*Bx_old + C05*By_old + C06*Bz_old
+                                     + C08*G_old // only with div cleaning
+                                     + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0
+                                     + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1; // only with J linear in time
+
+                if (div_cleaning)
+                {
+                    // F
+                    C01 = I*kx*S*inv_knorm*inv_c;
+                    C02 = I*ky*S*inv_knorm*inv_c;
+                    C03 = I*kz*S*inv_knorm*inv_c;
+                    C07 = C;
+                    C09 = I*mu0*kx*(C-1._rt)*inv_knorm2;
+                    C10 = I*mu0*ky*(C-1._rt)*inv_knorm2;
+                    C11 = I*mu0*kz*(C-1._rt)*inv_knorm2;
+                    C12 = (J_linear) ? I*mu0*kx*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+                    C13 = (J_linear) ? I*mu0*ky*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+                    C14 = (J_linear) ? I*mu0*kz*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt;
+                    C15 = -mu0*c*S*inv_knorm;
+                    C16 = (rho_linear) ? mu0*(C-1._rt)*inv_knorm2 : 0._rt;
+
+                    fields(i,j,k,Idx.F) = C01*Ex_old + C02*Ey_old + C03*Ez_old
+                                        + C07*F_old
+                                        + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0
+                                        + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time
+                                        + C15*rho_c0
+                                        + C16*rho_c1; // only with rho linear in time
+
+                    // G
+                    C04 = I*c*kx*S*inv_knorm;
+                    C05 = I*c*ky*S*inv_knorm;
+                    C06 = I*c*kz*S*inv_knorm;
+                    C08 = C;
+
+                    fields(i,j,k,Idx.G) = C04*Bx_old + C05*By_old + C06*Bz_old
+                                        + C08*G_old;
+                }
+            }
+        });
+    }
+}
+
+void PsatdAlgorithmFirstOrder::CurrentCorrection (SpectralFieldData& field_data)
+{
+    // Profiling
+    BL_PROFILE("PsatdAlgorithmFirstOrder::CurrentCorrection");
+
+    amrex::ignore_unused(field_data);
+    amrex::Abort(Utils::TextMsg::Err(
+        "Current correction not implemented for first-order PSATD equations"));
+}
+
+void
+PsatdAlgorithmFirstOrder::VayDeposition (SpectralFieldData& field_data)
+{
+    // Profiling
+    BL_PROFILE("PsatdAlgorithmFirstOrder::VayDeposition()");
+
+    amrex::ignore_unused(field_data);
+    amrex::Abort(Utils::TextMsg::Err(
+        "Vay deposition not implemented for first-order PSATD equations"));
+}
+
+#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp
index be850e252..af3e2e336 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp
@@ -103,9 +103,9 @@ PsatdAlgorithmGalileanRZ::pushSpectralFields (SpectralFieldDataRZ & f)
             auto const Bp_m = Idx.Bx + Idx.n_fields*mode;
             auto const Bm_m = Idx.By + Idx.n_fields*mode;
             auto const Bz_m = Idx.Bz + Idx.n_fields*mode;
-            auto const Jp_m = Idx.Jx + Idx.n_fields*mode;
-            auto const Jm_m = Idx.Jy + Idx.n_fields*mode;
-            auto const Jz_m = Idx.Jz + Idx.n_fields*mode;
+            auto const Jp_m = Idx.Jx_mid + Idx.n_fields*mode;
+            auto const Jm_m = Idx.Jy_mid + Idx.n_fields*mode;
+            auto const Jz_m = Idx.Jz_mid + Idx.n_fields*mode;
             auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode;
             auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode;
 
@@ -323,9 +323,9 @@ PsatdAlgorithmGalileanRZ::CurrentCorrection (SpectralFieldDataRZ& field_data)
         [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept
         {
             // All of the fields of each mode are grouped together
-            auto const Jp_m = Idx.Jx + Idx.n_fields*mode;
-            auto const Jm_m = Idx.Jy + Idx.n_fields*mode;
-            auto const Jz_m = Idx.Jz + Idx.n_fields*mode;
+            auto const Jp_m = Idx.Jx_mid + Idx.n_fields*mode;
+            auto const Jm_m = Idx.Jy_mid + Idx.n_fields*mode;
+            auto const Jz_m = Idx.Jz_mid + Idx.n_fields*mode;
             auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode;
             auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp
index 9d1fbf3e1..c52437bf8 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp
@@ -190,9 +190,9 @@ PsatdAlgorithmJConstantInTime::pushSpectralFields (SpectralFieldData& f) const
             const Complex Bz_old = fields(i,j,k,Idx.Bz);
 
             // Shortcuts for the values of J
-            const Complex Jx = fields(i,j,k,Idx.Jx);
-            const Complex Jy = fields(i,j,k,Idx.Jy);
-            const Complex Jz = fields(i,j,k,Idx.Jz);
+            const Complex Jx = fields(i,j,k,Idx.Jx_mid);
+            const Complex Jy = fields(i,j,k,Idx.Jy_mid);
+            const Complex Jz = fields(i,j,k,Idx.Jz_mid);
 
             Complex F_old;
             if (dive_cleaning)
@@ -751,9 +751,9 @@ void PsatdAlgorithmJConstantInTime::CurrentCorrection (SpectralFieldData& field_
         ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             // Shortcuts for the values of J and rho
-            const Complex Jx = fields(i,j,k,Idx.Jx);
-            const Complex Jy = fields(i,j,k,Idx.Jy);
-            const Complex Jz = fields(i,j,k,Idx.Jz);
+            const Complex Jx = fields(i,j,k,Idx.Jx_mid);
+            const Complex Jy = fields(i,j,k,Idx.Jy_mid);
+            const Complex Jz = fields(i,j,k,Idx.Jz_mid);
             const Complex rho_old = fields(i,j,k,Idx.rho_old);
             const Complex rho_new = fields(i,j,k,Idx.rho_new);
 
@@ -787,25 +787,25 @@ void PsatdAlgorithmJConstantInTime::CurrentCorrection (SpectralFieldData& field_
                     const Complex rho_old_mod = rho_old * amrex::exp(I * k_dot_vg * dt);
                     const Complex den = 1._rt - amrex::exp(I * k_dot_vg * dt);
 
-                    fields(i,j,k,Idx.Jx) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
+                    fields(i,j,k,Idx.Jx_mid) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
                         * kx / (k_norm * k_norm);
 
-                    fields(i,j,k,Idx.Jy) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
+                    fields(i,j,k,Idx.Jy_mid) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
                         * ky / (k_norm * k_norm);
 
-                    fields(i,j,k,Idx.Jz) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
+                    fields(i,j,k,Idx.Jz_mid) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den)
                         * kz / (k_norm * k_norm);
                 }
 
                 else
                 {
-                    fields(i,j,k,Idx.Jx) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt)
+                    fields(i,j,k,Idx.Jx_mid) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt)
                         * kx / (k_norm * k_norm);
 
-                    fields(i,j,k,Idx.Jy) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt)
+                    fields(i,j,k,Idx.Jy_mid) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt)
                         * ky / (k_norm * k_norm);
 
-                    fields(i,j,k,Idx.Jz) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt)
+                    fields(i,j,k,Idx.Jz_mid) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt)
                         * kz / (k_norm * k_norm);
                 }
             }
@@ -840,11 +840,11 @@ PsatdAlgorithmJConstantInTime::VayDeposition (SpectralFieldData& field_data)
         ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             // Shortcuts for the values of D
-            const Complex Dx = fields(i,j,k,Idx.Jx);
+            const Complex Dx = fields(i,j,k,Idx.Jx_mid);
 #if defined(WARPX_DIM_3D)
-            const Complex Dy = fields(i,j,k,Idx.Jy);
+            const Complex Dy = fields(i,j,k,Idx.Jy_mid);
 #endif
-            const Complex Dz = fields(i,j,k,Idx.Jz);
+            const Complex Dz = fields(i,j,k,Idx.Jz_mid);
 
             // Imaginary unit
             constexpr Complex I = Complex{0._rt, 1._rt};
@@ -859,18 +859,18 @@ PsatdAlgorithmJConstantInTime::VayDeposition (SpectralFieldData& field_data)
 #endif
 
             // Compute Jx
-            if (kx_mod != 0._rt) fields(i,j,k,Idx.Jx) = I * Dx / kx_mod;
-            else                 fields(i,j,k,Idx.Jx) = 0._rt;
+            if (kx_mod != 0._rt) fields(i,j,k,Idx.Jx_mid) = I * Dx / kx_mod;
+            else                 fields(i,j,k,Idx.Jx_mid) = 0._rt;
 
 #if defined(WARPX_DIM_3D)
             // Compute Jy
-            if (ky_mod != 0._rt) fields(i,j,k,Idx.Jy) = I * Dy / ky_mod;
-            else                 fields(i,j,k,Idx.Jy) = 0._rt;
+            if (ky_mod != 0._rt) fields(i,j,k,Idx.Jy_mid) = I * Dy / ky_mod;
+            else                 fields(i,j,k,Idx.Jy_mid) = 0._rt;
 #endif
 
             // Compute Jz
-            if (kz_mod != 0._rt) fields(i,j,k,Idx.Jz) = I * Dz / kz_mod;
-            else                 fields(i,j,k,Idx.Jz) = 0._rt;
+            if (kz_mod != 0._rt) fields(i,j,k,Idx.Jz_mid) = I * Dz / kz_mod;
+            else                 fields(i,j,k,Idx.Jz_mid) = 0._rt;
         });
     }
 }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp
index bd9df977e..861001c78 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp
@@ -120,9 +120,9 @@ PsatdAlgorithmJLinearInTime::pushSpectralFields (SpectralFieldData& f) const
             const Complex Bz_old = fields(i,j,k,Idx.Bz);
 
             // Shortcuts for the values of J and rho
-            const Complex Jx_old = fields(i,j,k,Idx.Jx);
-            const Complex Jy_old = fields(i,j,k,Idx.Jy);
-            const Complex Jz_old = fields(i,j,k,Idx.Jz);
+            const Complex Jx_old = fields(i,j,k,Idx.Jx_old);
+            const Complex Jy_old = fields(i,j,k,Idx.Jy_old);
+            const Complex Jz_old = fields(i,j,k,Idx.Jz_old);
             const Complex Jx_new = fields(i,j,k,Idx.Jx_new);
             const Complex Jy_new = fields(i,j,k,Idx.Jy_new);
             const Complex Jz_new = fields(i,j,k,Idx.Jz_new);
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
index 55b58821c..effb1cc2b 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
@@ -83,7 +83,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 
     const bool update_with_rho = m_update_with_rho;
     const bool time_averaging = m_time_averaging;
-    const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false;
+    const bool J_linear = (m_J_in_time == JInTime::Linear) ? true : false;
     const bool dive_cleaning = m_dive_cleaning;
     const bool divb_cleaning = m_divb_cleaning;
 
@@ -112,7 +112,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 
         amrex::Array4<const amrex::Real> X5_arr;
         amrex::Array4<const amrex::Real> X6_arr;
-        if (time_averaging && J_in_time_linear)
+        if (time_averaging && J_linear)
         {
             X5_arr = X5_coef[mfi].array();
             X6_arr = X6_coef[mfi].array();
@@ -131,6 +131,9 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
         amrex::ParallelFor(bx, modes,
         [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept
         {
+            int idx_jx = (J_linear) ? static_cast<int>(Idx.Jx_old) : static_cast<int>(Idx.Jx_mid);
+            int idx_jy = (J_linear) ? static_cast<int>(Idx.Jy_old) : static_cast<int>(Idx.Jy_mid);
+            int idx_jz = (J_linear) ? static_cast<int>(Idx.Jz_old) : static_cast<int>(Idx.Jz_mid);
 
             // All of the fields of each mode are grouped together
             int const Ep_m = Idx.Ex + Idx.n_fields*mode;
@@ -139,9 +142,9 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
             int const Bp_m = Idx.Bx + Idx.n_fields*mode;
             int const Bm_m = Idx.By + Idx.n_fields*mode;
             int const Bz_m = Idx.Bz + Idx.n_fields*mode;
-            int const Jp_m = Idx.Jx + Idx.n_fields*mode;
-            int const Jm_m = Idx.Jy + Idx.n_fields*mode;
-            int const Jz_m = Idx.Jz + Idx.n_fields*mode;
+            int const Jp_m = idx_jx + Idx.n_fields*mode;
+            int const Jm_m = idx_jy + Idx.n_fields*mode;
+            int const Jz_m = idx_jz + Idx.n_fields*mode;
             int const rho_old_m = Idx.rho_old + Idx.n_fields*mode;
             int const rho_new_m = Idx.rho_new + Idx.n_fields*mode;
 
@@ -238,7 +241,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
                 G_old = fields(i,j,k,G_m);
             }
 
-            if (J_in_time_linear)
+            if (J_linear)
             {
                 const int Jp_m_new = Idx.Jx_new + Idx.n_fields*mode;
                 const int Jm_m_new = Idx.Jy_new + Idx.n_fields*mode;
@@ -335,7 +338,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f)
 {
     const bool time_averaging = m_time_averaging;
-    const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false;
+    const bool J_linear = (m_J_in_time == JInTime::Linear) ? true : false;
 
     // Fill them with the right values:
     // Loop over boxes and allocate the corresponding coefficients
@@ -356,7 +359,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
 
         amrex::Array4<amrex::Real> X5;
         amrex::Array4<amrex::Real> X6;
-        if (time_averaging && J_in_time_linear)
+        if (time_averaging && J_linear)
         {
             X5 = X5_coef[mfi].array();
             X6 = X6_coef[mfi].array();
@@ -395,7 +398,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
                 X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0);
             }
 
-            if (time_averaging && J_in_time_linear)
+            if (time_averaging && J_linear)
             {
                 constexpr amrex::Real c2 = PhysConst::c;
                 const amrex::Real dt3 = dt * dt * dt;
@@ -450,9 +453,9 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data)
         [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept
         {
             // All of the fields of each mode are grouped together
-            auto const Jp_m = Idx.Jx + Idx.n_fields*mode;
-            auto const Jm_m = Idx.Jy + Idx.n_fields*mode;
-            auto const Jz_m = Idx.Jz + Idx.n_fields*mode;
+            auto const Jp_m = Idx.Jx_mid + Idx.n_fields*mode;
+            auto const Jm_m = Idx.Jy_mid + Idx.n_fields*mode;
+            auto const Jz_m = Idx.Jz_mid + Idx.n_fields*mode;
             auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode;
             auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
index 4ab88f1a3..c7848d731 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -83,18 +83,19 @@ class SpectralFieldIndex
     // Always
     int Ex = -1, Ey = -1, Ez = -1;
     int Bx = -1, By = -1, Bz = -1;
-    int Jx = -1, Jy = -1, Jz = -1;
-    int rho_old = -1, rho_new = -1, divE = -1;
+    int divE = -1;
 
     // Time averaging
     int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1;
     int Bx_avg = -1, By_avg = -1, Bz_avg = -1;
 
-    // J linear in time
+    // J
+    int Jx_old = -1, Jy_old = -1, Jz_old = -1;
+    int Jx_mid = -1, Jy_mid = -1, Jz_mid = -1;
     int Jx_new = -1, Jy_new = -1, Jz_new = -1;
 
-    // rho quadratic in time
-    int rho_mid = -1;
+    // rho
+    int rho_old = -1, rho_mid = -1, rho_new = -1;
 
     // div(E) and div(B) cleaning
     int F = -1, G = -1;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
index 56189e0ff..e1db4e5ab 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -51,10 +51,8 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho,
     {
         Ex = c++; Ey = c++; Ez = c++;
         Bx = c++; By = c++; Bz = c++;
-        Jx = c++; Jy = c++; Jz = c++;
 
         // TODO Allocate rho_old and rho_new only when needed
-        rho_old = c++; rho_new = c++;
 
         // Reuse data corresponding to index Bx = 3 to avoid storing extra memory
         divE = 3;
@@ -69,17 +67,25 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho,
 
         if (divb_cleaning) G = c++;
 
-        if (J_in_time == JInTime::Linear)
+        if (J_in_time == JInTime::Constant)
         {
-            Jx_new = c++;
-            Jy_new = c++;
-            Jz_new = c++;
+            Jx_mid = c++; Jy_mid = c++; Jz_mid = c++;
+        }
+        else if (J_in_time == JInTime::Linear)
+        {
+            Jx_old = c++; Jy_old = c++; Jz_old = c++;
+            Jx_new = c++; Jy_new = c++; Jz_new = c++;
         }
 
-        if (rho_in_time == RhoInTime::Quadratic)
+        if (rho_in_time == RhoInTime::Constant)
         {
             rho_mid = c++;
         }
+        else if (rho_in_time == RhoInTime::Linear)
+        {
+            rho_old = c++;
+            rho_new = c++;
+        }
 
         if (pml_rz)
         {
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 38b242010..da4b9687b 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -58,6 +58,8 @@ class SpectralSolver
          *                                (no domain decomposition)
          * \param[in] update_with_rho whether rho is used in the field update equations
          * \param[in] fft_do_time_averaging whether the time averaging algorithm is used
+         * \param[in] psatd_solution_type whether the PSATD equations are derived
+         *                                from a first-order or second-order model
          * \param[in] J_in_time integer that corresponds to the time dependency of J
          *                      (constant, linear) for the PSATD algorithm
          * \param[in] rho_in_time integer that corresponds to the time dependency of rho
@@ -80,6 +82,7 @@ class SpectralSolver
                         const bool periodic_single_box,
                         const bool update_with_rho,
                         const bool fft_do_time_averaging,
+                        const int psatd_solution_type,
                         const int J_in_time,
                         const int rho_in_time,
                         const bool dive_cleaning,
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index c0d7b8941..f29a0bd03 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -8,10 +8,12 @@
 #include "FieldSolver/SpectralSolver/SpectralFieldData.H"
 #include "SpectralAlgorithms/PsatdAlgorithmComoving.H"
 #include "SpectralAlgorithms/PsatdAlgorithmPml.H"
+#include "SpectralAlgorithms/PsatdAlgorithmFirstOrder.H"
 #include "SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H"
 #include "SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H"
 #include "SpectralKSpace.H"
 #include "SpectralSolver.H"
+#include "Utils/TextMsg.H"
 #include "Utils/WarpXProfilerWrapper.H"
 
 #include <memory>
@@ -30,6 +32,7 @@ SpectralSolver::SpectralSolver(
                 const bool pml, const bool periodic_single_box,
                 const bool update_with_rho,
                 const bool fft_do_time_averaging,
+                const int psatd_solution_type,
                 const int J_in_time,
                 const int rho_in_time,
                 const bool dive_cleaning,
@@ -49,13 +52,13 @@ SpectralSolver::SpectralSolver(
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
 
-    if (pml) // PSATD equations in the PML grids
+    if (pml) // PSATD equations in the PML region
     {
         algorithm = std::make_unique<PsatdAlgorithmPml>(
             k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
             dt, dive_cleaning, divb_cleaning);
     }
-    else // PSATD equations in the regulard grids
+    else // PSATD equations in the regular domain
     {
         // Comoving PSATD algorithm
         if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.)
@@ -64,7 +67,31 @@ SpectralSolver::SpectralSolver(
                 k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
                 v_comoving, dt, update_with_rho);
         }
-        else // PSATD algorithms: standard, Galilean, averaged Galilean, multi-J
+        // Galilean PSATD algorithm (only J constant in time)
+        else if (v_galilean[0] != 0. || v_galilean[1] != 0. || v_galilean[2] != 0.)
+        {
+            algorithm = std::make_unique<PsatdAlgorithmJConstantInTime>(
+                k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
+                v_galilean, dt, update_with_rho, fft_do_time_averaging,
+                dive_cleaning, divb_cleaning);
+        }
+        else if (psatd_solution_type == PSATDSolutionType::FirstOrder)
+        {
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                !fft_do_time_averaging,
+                "psatd.do_time_averaging=1 not supported when psatd.solution_type=first-order");
+
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                (!dive_cleaning && !divb_cleaning) || (dive_cleaning && divb_cleaning),
+                "warpx.do_dive_cleaning and warpx.do_divb_cleaning must be equal when psatd.solution_type=first-order");
+
+            const bool div_cleaning = (dive_cleaning && divb_cleaning);
+
+            algorithm = std::make_unique<PsatdAlgorithmFirstOrder>(
+                k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
+                dt, div_cleaning, J_in_time, rho_in_time);
+        }
+        else if (psatd_solution_type == PSATDSolutionType::SecondOrder)
         {
             if (J_in_time == JInTime::Constant)
             {
@@ -73,7 +100,7 @@ SpectralSolver::SpectralSolver(
                     v_galilean, dt, update_with_rho, fft_do_time_averaging,
                     dive_cleaning, divb_cleaning);
             }
-            else // J linear in time
+            else if (J_in_time == JInTime::Linear)
             {
                 algorithm = std::make_unique<PsatdAlgorithmJLinearInTime>(
                     k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,