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>

22 files changed, 687 insertions, 121 deletions

diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H
index a1ce21a91..e519298ca 100644
--- a/Source/BoundaryConditions/PML.H
+++ b/Source/BoundaryConditions/PML.H
@@ -133,7 +133,7 @@ public:
          int ncell, int delta, amrex::IntVect ref_ratio,
          amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
          int do_moving_window, int pml_has_particles, int do_pml_in_domain,
-         const int J_in_time, const int rho_in_time,
+         const int psatd_solution_type, const int J_in_time, const int rho_in_time,
          const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning,
          const amrex::IntVect& fill_guards_fields,
          const amrex::IntVect& fill_guards_current,
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index c7b773123..5e03822f1 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -548,7 +548,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
           int ncell, int delta, amrex::IntVect ref_ratio,
           Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
           int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain,
-          const int J_in_time, const int rho_in_time,
+          const int psatd_solution_type, const int J_in_time, const int rho_in_time,
           const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning,
           const amrex::IntVect& fill_guards_fields,
           const amrex::IntVect& fill_guards_current,
@@ -724,7 +724,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
 
     if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
 #ifndef WARPX_USE_PSATD
-        amrex::ignore_unused(lev, dt, J_in_time, rho_in_time);
+        amrex::ignore_unused(lev, dt, psatd_solution_type, J_in_time, rho_in_time);
 #   if(AMREX_SPACEDIM!=3)
         amrex::ignore_unused(noy_fft);
 #   endif
@@ -745,7 +745,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
         spectral_solver_fp = std::make_unique<SpectralSolver>(lev, realspace_ba, dm,
             nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero,
             v_comoving_zero, dx, dt, in_pml, periodic_single_box, update_with_rho,
-            fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning);
+            fft_do_time_averaging, psatd_solution_type, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning);
 #endif
     }
 
@@ -852,7 +852,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
             spectral_solver_cp = std::make_unique<SpectralSolver>(lev, realspace_cba, cdm,
                 nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero,
                 v_comoving_zero, cdx, dt, in_pml, periodic_single_box, update_with_rho,
-                fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning);
+                fft_do_time_averaging, psatd_solution_type, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning);
 #endif
         }
     }
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index e8e6a025b..8abfa0e7f 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -549,14 +549,14 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time)
 
     // 3) Deposit rho (in rho_new, since it will be moved during the loop)
     //    (after checking that pointer to rho_fp on MR level 0 is not null)
-    if (rho_fp[0])
+    if (rho_fp[0] && rho_in_time == RhoInTime::Linear)
     {
         // Deposit rho at relative time -dt
         // (dt[0] denotes the time step on mesh refinement level 0)
         mypc->DepositCharge(rho_fp, -dt[0]);
         // Filter, exchange boundary, and interpolate across levels
         SyncRho();
-        // Forward FFT of rho_new
+        // Forward FFT of rho
         PSATDForwardTransformRho(rho_fp, rho_cp, 0, 1);
     }
 
@@ -587,17 +587,14 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time)
     // Loop over multi-J depositions
     for (int i_depose = 0; i_depose < n_loop; i_depose++)
     {
-        // Move J deposited previously, from new to old
-        if (J_in_time == JInTime::Linear)
-        {
-            PSATDMoveJNewToJOld();
-        }
+        // Move J from new to old if J is linear in time
+        if (J_in_time == JInTime::Linear) PSATDMoveJNewToJOld();
 
         const amrex::Real t_depose_current = (J_in_time == JInTime::Linear) ?
             (i_depose-n_depose+1)*sub_dt : (i_depose-n_depose+0.5_rt)*sub_dt;
 
-        // TODO Update this when rho quadratic in time is implemented
-        const amrex::Real t_depose_charge = (i_depose-n_depose+1)*sub_dt;
+        const amrex::Real t_depose_charge = (rho_in_time == RhoInTime::Linear) ?
+            (i_depose-n_depose+1)*sub_dt : (i_depose-n_depose+0.5_rt)*sub_dt;
 
         // Deposit new J at relative time t_depose_current with time step dt
         // (dt[0] denotes the time step on mesh refinement level 0)
@@ -616,14 +613,14 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time)
         // (after checking that pointer to rho_fp on MR level 0 is not null)
         if (rho_fp[0])
         {
-            // Move rho deposited previously, from new to old
-            PSATDMoveRhoNewToRhoOld();
+            // Move rho from new to old if rho is linear in time
+            if (rho_in_time == RhoInTime::Linear) PSATDMoveRhoNewToRhoOld();
 
             // Deposit rho at relative time t_depose_charge
             mypc->DepositCharge(rho_fp, t_depose_charge);
             // Filter, exchange boundary, and interpolate across levels
             SyncRho();
-            // Forward FFT of rho_new
+            // Forward FFT of rho
             PSATDForwardTransformRho(rho_fp, rho_cp, 0, 1);
         }
 
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,
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 3d68e8e52..9df4bb21b 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -280,9 +280,9 @@ void WarpX::PSATDForwardTransformJ (
     {
         Idx = spectral_solver_fp[lev]->m_spectral_index;
 
-        idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx);
-        idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy);
-        idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz);
+        idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid);
+        idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid);
+        idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid);
 
         ForwardTransformVect(lev, *spectral_solver_fp[lev], J_fp[lev], idx_jx, idx_jy, idx_jz);
 
@@ -290,9 +290,9 @@ void WarpX::PSATDForwardTransformJ (
         {
             Idx = spectral_solver_cp[lev]->m_spectral_index;
 
-            idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx);
-            idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy);
-            idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz);
+            idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid);
+            idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid);
+            idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid);
 
             ForwardTransformVect(lev, *spectral_solver_cp[lev], J_cp[lev], idx_jx, idx_jy, idx_jz);
         }
@@ -304,11 +304,23 @@ void WarpX::PSATDForwardTransformJ (
     {
         for (int lev = 0; lev <= finest_level; ++lev)
         {
-            spectral_solver_fp[lev]->ApplyFilter(lev, Idx.Jx, Idx.Jy, Idx.Jz);
+            Idx = spectral_solver_fp[lev]->m_spectral_index;
+
+            idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid);
+            idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid);
+            idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid);
+
+            spectral_solver_fp[lev]->ApplyFilter(lev, idx_jx, idx_jy, idx_jz);
 
             if (spectral_solver_cp[lev])
             {
-                spectral_solver_cp[lev]->ApplyFilter(lev, Idx.Jx, Idx.Jy, Idx.Jz);
+                Idx = spectral_solver_cp[lev]->m_spectral_index;
+
+                idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid);
+                idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid);
+                idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid);
+
+                spectral_solver_cp[lev]->ApplyFilter(lev, idx_jx, idx_jy, idx_jz);
             }
         }
     }
@@ -328,9 +340,11 @@ void WarpX::PSATDBackwardTransformJ (
     {
         Idx = spectral_solver_fp[lev]->m_spectral_index;
 
-        idx_jx = static_cast<int>(Idx.Jx);
-        idx_jy = static_cast<int>(Idx.Jy);
-        idx_jz = static_cast<int>(Idx.Jz);
+        // Note that these backward FFTs are currently called only
+        // with algorithms that do not support J linear in time
+        idx_jx = static_cast<int>(Idx.Jx_mid);
+        idx_jy = static_cast<int>(Idx.Jy_mid);
+        idx_jz = static_cast<int>(Idx.Jz_mid);
 
         BackwardTransformVect(lev, *spectral_solver_fp[lev], J_fp[lev],
                               idx_jx, idx_jy, idx_jz, m_fill_guards_current);
@@ -339,9 +353,11 @@ void WarpX::PSATDBackwardTransformJ (
         {
             Idx = spectral_solver_cp[lev]->m_spectral_index;
 
-            idx_jx = static_cast<int>(Idx.Jx);
-            idx_jy = static_cast<int>(Idx.Jy);
-            idx_jz = static_cast<int>(Idx.Jz);
+            // Note that these backward FFTs are currently called only
+            // with algorithms that do not support J linear in time
+            idx_jx = static_cast<int>(Idx.Jx_mid);
+            idx_jy = static_cast<int>(Idx.Jy_mid);
+            idx_jz = static_cast<int>(Idx.Jz_mid);
 
             BackwardTransformVect(lev, *spectral_solver_cp[lev], J_cp[lev],
                                   idx_jx, idx_jy, idx_jz, m_fill_guards_current);
@@ -359,7 +375,15 @@ void WarpX::PSATDForwardTransformRho (
     const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index;
 
     // Select index in k space
-    const int dst_comp = (dcomp == 0) ? Idx.rho_old : Idx.rho_new;
+    int dst_comp;
+    if (rho_in_time == RhoInTime::Constant)
+    {
+        dst_comp = Idx.rho_mid;
+    }
+    else // rho_in_time == RhoInTime::Linear
+    {
+        dst_comp = (dcomp == 0) ? Idx.rho_old : Idx.rho_new;
+    }
 
     for (int lev = 0; lev <= finest_level; ++lev)
     {
@@ -568,15 +592,15 @@ WarpX::PSATDMoveJNewToJOld ()
 
     for (int lev = 0; lev <= finest_level; ++lev)
     {
-        spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx);
-        spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy);
-        spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz);
+        spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx_old);
+        spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy_old);
+        spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz_old);
 
         if (spectral_solver_cp[lev])
         {
-            spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx);
-            spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy);
-            spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz);
+            spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx_old);
+            spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy_old);
+            spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz_old);
         }
     }
 }
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 0d3e919a0..32cd11dfd 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -480,7 +480,7 @@ WarpX::InitPML ()
                              pml_ncell, pml_delta, amrex::IntVect::TheZeroVector(),
                              dt[0], nox_fft, noy_fft, noz_fft, do_nodal,
                              do_moving_window, pml_has_particles, do_pml_in_domain,
-                             J_in_time, rho_in_time,
+                             psatd_solution_type, J_in_time, rho_in_time,
                              do_pml_dive_cleaning, do_pml_divb_cleaning,
                              amrex::IntVect(0), amrex::IntVect(0),
                              guard_cells.ng_FieldSolver.max(),
@@ -517,7 +517,7 @@ WarpX::InitPML ()
                                    pml_ncell, pml_delta, refRatio(lev-1),
                                    dt[lev], nox_fft, noy_fft, noz_fft, do_nodal,
                                    do_moving_window, pml_has_particles, do_pml_in_domain,
-                                   J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning,
+                                   psatd_solution_type, J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning,
                                    amrex::IntVect(0), amrex::IntVect(0),
                                    guard_cells.ng_FieldSolver.max(),
                                    v_particle_pml,
diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H
index c3160cfad..936fd3d2b 100644
--- a/Source/Utils/WarpXAlgorithmSelection.H
+++ b/Source/Utils/WarpXAlgorithmSelection.H
@@ -84,6 +84,13 @@ struct GatheringAlgo {
     };
 };
 
+struct PSATDSolutionType {
+    enum {
+        FirstOrder = 0,
+        SecondOrder = 1
+    };
+};
+
 struct JInTime {
     enum {
         Constant = 0,
@@ -93,8 +100,8 @@ struct JInTime {
 
 struct RhoInTime {
     enum {
-        Linear = 1,
-        Quadratic = 2
+        Constant = 0,
+        Linear = 1
     };
 };
 
diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp
index b99459b46..c5ef16743 100644
--- a/Source/Utils/WarpXAlgorithmSelection.cpp
+++ b/Source/Utils/WarpXAlgorithmSelection.cpp
@@ -63,6 +63,12 @@ const std::map<std::string, int> gathering_algo_to_int = {
     {"default",             GatheringAlgo::EnergyConserving }
 };
 
+const std::map<std::string, int> psatd_solution_type_to_int = {
+    {"first-order", PSATDSolutionType::FirstOrder},
+    {"second-order", PSATDSolutionType::SecondOrder},
+    {"default", PSATDSolutionType::SecondOrder}
+};
+
 const std::map<std::string, int> J_in_time_to_int = {
     {"constant", JInTime::Constant},
     {"linear", JInTime::Linear},
@@ -70,8 +76,8 @@ const std::map<std::string, int> J_in_time_to_int = {
 };
 
 const std::map<std::string, int> rho_in_time_to_int = {
+    {"constant", RhoInTime::Constant},
     {"linear", RhoInTime::Linear},
-    {"quadratic", RhoInTime::Quadratic},
     {"default", RhoInTime::Linear}
 };
 
@@ -145,6 +151,8 @@ GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){
         algo_to_int = charge_deposition_algo_to_int;
     } else if (0 == std::strcmp(pp_search_key, "field_gathering")) {
         algo_to_int = gathering_algo_to_int;
+    } else if (0 == std::strcmp(pp_search_key, "solution_type")) {
+        algo_to_int = psatd_solution_type_to_int;
     } else if (0 == std::strcmp(pp_search_key, "J_in_time")) {
         algo_to_int = J_in_time_to_int;
     } else if (0 == std::strcmp(pp_search_key, "rho_in_time")) {
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 20d4e4210..2d1124c06 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -187,6 +187,11 @@ public:
      */
     static amrex::Vector<ParticleBoundaryType> particle_boundary_hi;
 
+    //! Integer that corresponds to the order of the PSATD solution
+    //! (whether the PSATD equations are derived from first-order or
+    //! second-order solution)
+    static short psatd_solution_type;
+
     //! Integers that correspond to the time dependency of J (constant, linear)
     //! and rho (linear, quadratic) for the PSATD algorithm
     static short J_in_time;
@@ -1639,7 +1644,7 @@ private:
         const int icomp, const int dcomp, const bool apply_kspace_filter=true);
 
     /**
-     * \brief Copy rho_new to rho_old in spectral space
+     * \brief Copy rho_new to rho_old in spectral space (when rho is linear in time)
      */
     void PSATDMoveRhoNewToRhoOld ();
 
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 579387433..e51197089 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -120,6 +120,7 @@ short WarpX::charge_deposition_algo;
 short WarpX::field_gathering_algo;
 short WarpX::particle_pusher_algo;
 short WarpX::electromagnetic_solver_id;
+short WarpX::psatd_solution_type;
 short WarpX::J_in_time;
 short WarpX::rho_in_time;
 short WarpX::load_balance_costs_update_algo;
@@ -1148,6 +1149,11 @@ WarpX::ReadParameters ()
             WARPX_ALWAYS_ASSERT_WITH_MESSAGE(noz_fft > 0, "PSATD order must be finite unless psatd.periodic_single_box_fft is used");
         }
 
+        // Integer that corresponds to the order of the PSATD solution
+        // (whether the PSATD equations are derived from first-order or
+        // second-order solution)
+        psatd_solution_type = GetAlgorithmInteger(pp_psatd, "solution_type");
+
         // Integers that correspond to the time dependency of J (constant, linear)
         // and rho (linear, quadratic) for the PSATD algorithm
         J_in_time = GetAlgorithmInteger(pp_psatd, "J_in_time");
@@ -1311,13 +1317,6 @@ WarpX::ReadParameters ()
             );
         }
 
-        if (J_in_time == JInTime::Constant)
-        {
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                rho_in_time == RhoInTime::Linear,
-                "psatd.J_in_time=constant supports only psatd.rho_in_time=linear");
-        }
-
         if (J_in_time == JInTime::Linear)
         {
             WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
@@ -1331,6 +1330,14 @@ WarpX::ReadParameters ()
             WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
                 v_comoving_is_zero,
                 "psatd.J_in_time=linear not implemented with comoving PSATD");
+
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                !current_correction,
+                "psatd.current_correction=1 not implemented with psatd.J_in_time=linear");
+
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                current_deposition_algo != CurrentDepositionAlgo::Vay,
+                "algo.current_deposition=vay not implemented with psatd.J_in_time=linear");
         }
 
         for (int dir = 0; dir < AMREX_SPACEDIM; dir++)
@@ -2334,6 +2341,7 @@ void WarpX::AllocLevelSpectralSolver (amrex::Vector<std::unique_ptr<SpectralSolv
                                                 fft_periodic_single_box,
                                                 update_with_rho,
                                                 fft_do_time_averaging,
+                                                psatd_solution_type,
                                                 J_in_time,
                                                 rho_in_time,
                                                 do_dive_cleaning,