diff --git a/Docs/source/usage/PWFA.png b/Docs/source/usage/PWFA.png
deleted file mode 100644
index f2299ec1ff9..00000000000
Binary files a/Docs/source/usage/PWFA.png and /dev/null differ
diff --git a/Docs/source/usage/examples.rst b/Docs/source/usage/examples.rst
index 0492372b4e6..dec8e9bcfee 100644
--- a/Docs/source/usage/examples.rst
+++ b/Docs/source/usage/examples.rst
@@ -23,7 +23,13 @@ Plasma-Based Acceleration
 
    examples/lwfa/README.rst
    examples/pwfa/README.rst
-   pwfa.rst
+
+Coming soon:
+
+* LWFA: External injection in the boosted frame
+* LWFA: Ionization injection in the lab frame using a LASY data file
+* PWFA: Self-injection in the lab frame
+* MR case?
 
 
 Laser-Plasma Interaction
diff --git a/Docs/source/usage/faq.rst b/Docs/source/usage/faq.rst
index 67cea8d6621..3805b935e2b 100644
--- a/Docs/source/usage/faq.rst
+++ b/Docs/source/usage/faq.rst
@@ -63,6 +63,7 @@ An in-depth discussion of the boosted frame is provided in the :ref:`moving wind
 
 What about Back-transformed diagnostics (BTD)?
 ----------------------------------------------
+.. _faq-btd:
 
 .. figure:: https://user-images.githubusercontent.com/10621396/198702232-9dd595ad-479e-4170-bd25-51e2b72cd50a.png
    :alt: [fig:BTD_features] Minkowski diagram indicating several features of the back-transformed diagnostic (BTD). The diagram explains why the first BTD begins to fill at boosted time :math:`t'=0` but this doesn't necessarily correspond to lab time :math:`t=0`, how the BTD grid-spacing is determined by the boosted time step :math:`\Delta t'`, hence why the snapshot length don't correspond to the grid spacing and length in the input script, and how the BTD snapshots complete when the effective snapshot length is covered in the boosted frame.
diff --git a/Docs/source/usage/python.rst b/Docs/source/usage/python.rst
index 4d53b311764..6f81d6a0504 100644
--- a/Docs/source/usage/python.rst
+++ b/Docs/source/usage/python.rst
@@ -95,6 +95,7 @@ Diagnostics
 .. autoclass:: pywarpx.picmi.Checkpoint
 
 .. autoclass:: pywarpx.picmi.ReducedDiagnostic
+.. _python_lab_frame_diag:
 
 Lab-frame diagnostics diagnostics are used when running boosted-frame simulations.
 
diff --git a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py
deleted file mode 100755
index 296aea48b35..00000000000
--- a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py
+++ /dev/null
@@ -1,80 +0,0 @@
-#!/usr/bin/env python3
-
-from pywarpx import picmi
-
-#from warp import picmi
-
-constants = picmi.constants
-
-nx = 64
-ny = 64
-nz = 64
-
-xmin = -200.e-6
-xmax = +200.e-6
-ymin = -200.e-6
-ymax = +200.e-6
-zmin = -200.e-6
-zmax = +200.e-6
-
-moving_window_velocity = [0., 0., constants.c]
-
-number_per_cell_each_dim = [2, 2, 1]
-
-max_steps = 10
-
-grid = picmi.Cartesian3DGrid(number_of_cells = [nx, ny, nz],
-                             lower_bound = [xmin, ymin, zmin],
-                             upper_bound = [xmax, ymax, zmax],
-                             lower_boundary_conditions = ['periodic', 'periodic', 'open'],
-                             upper_boundary_conditions = ['periodic', 'periodic', 'open'],
-                             lower_boundary_conditions_particles = ['periodic', 'periodic', 'absorbing'],
-                             upper_boundary_conditions_particles = ['periodic', 'periodic', 'absorbing'],
-                             moving_window_velocity = moving_window_velocity,
-                             warpx_max_grid_size=32)
-
-solver = picmi.ElectromagneticSolver(grid=grid, cfl=1)
-
-beam_distribution = picmi.UniformDistribution(density = 1.e23,
-                                              lower_bound = [-20.e-6, -20.e-6, -150.e-6],
-                                              upper_bound = [+20.e-6, +20.e-6, -100.e-6],
-                                              directed_velocity = [0., 0., 1.e9])
-
-plasma_distribution = picmi.UniformDistribution(density = 1.e22,
-                                                lower_bound = [-200.e-6, -200.e-6, 0.],
-                                                upper_bound = [+200.e-6, +200.e-6, None],
-                                                fill_in = True)
-
-beam = picmi.Species(particle_type='electron', name='beam', initial_distribution=beam_distribution)
-plasma = picmi.Species(particle_type='electron', name='plasma', initial_distribution=plasma_distribution)
-
-sim = picmi.Simulation(solver = solver,
-                       max_steps = max_steps,
-                       verbose = 1,
-                       warpx_current_deposition_algo = 'esirkepov',
-                       warpx_use_filter = 0)
-
-sim.add_species(beam, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=number_per_cell_each_dim))
-sim.add_species(plasma, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=number_per_cell_each_dim))
-
-field_diag = picmi.FieldDiagnostic(name = 'diag1',
-                                   grid = grid,
-                                   period = max_steps,
-                                   data_list = ['Ex', 'Ey', 'Ez', 'Jx', 'Jy', 'Jz', 'part_per_cell'],
-                                   write_dir = '.',
-                                   warpx_file_prefix = 'Python_PlasmaAcceleration_plt')
-
-part_diag = picmi.ParticleDiagnostic(name = 'diag1',
-                                     period = max_steps,
-                                     species = [beam, plasma],
-                                     data_list = ['ux', 'uy', 'uz', 'weighting'])
-
-sim.add_diagnostic(field_diag)
-sim.add_diagnostic(part_diag)
-
-# write_inputs will create an inputs file that can be used to run
-# with the compiled version.
-#sim.write_input_file(file_name = 'inputs_from_PICMI')
-
-# Alternatively, sim.step will run WarpX, controlling it from Python
-sim.step()
diff --git a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_3d.py b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_3d.py
new file mode 100755
index 00000000000..12cb6da6bf9
--- /dev/null
+++ b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_3d.py
@@ -0,0 +1,172 @@
+#!/usr/bin/env python3
+
+from pywarpx import picmi
+
+c = picmi.constants.c
+q_e = picmi.constants.q_e
+m_e = picmi.constants.m_e
+m_p = picmi.constants.m_p
+ep0 = picmi.constants.ep0
+
+nx = 64
+nz = 64
+
+xmin = -100.e-6
+xmax = +100.e-6
+zmin = -149.e-6
+zmax = +1.e-6
+
+moving_window_velocity = [0., 0., c]
+gamma_boost = 10.0
+
+npc_x = 1
+npc_z = 1
+
+max_steps = 10
+
+field_bc_z='open'
+field_bc_x='open'
+grid = picmi.Cartesian3DGrid(
+    number_of_cells = [nx, nx, nz],
+    lower_bound = [xmin, xmin, zmin],
+    upper_bound = [xmax, xmax, zmax],
+    lower_boundary_conditions = [field_bc_x, field_bc_x, field_bc_z],
+    upper_boundary_conditions = [field_bc_x, field_bc_x, field_bc_z],
+    lower_boundary_conditions_particles = ['absorbing']*3,
+    upper_boundary_conditions_particles = ['absorbing']*3,
+    moving_window_velocity = moving_window_velocity,
+    warpx_blocking_factor = 32,
+    warpx_max_grid_size=128
+)
+
+smoother = picmi.BinomialSmoother(n_pass=[1]*3)
+solver = picmi.ElectromagneticSolver(
+    method='CKC',
+    grid=grid,
+    cfl=0.99,
+    source_smoother=smoother,
+    warpx_pml_ncell=10)
+
+drive_charge = -10e-11
+driver_gaussian_distribution = picmi.GaussianBunchDistribution(
+    n_physical_particles= abs(drive_charge) / q_e,
+    rms_bunch_size=[2e-6, 2e-6, 4.e-6],
+    rms_velocity=[2.*c, 2.*c, 2.e4*c],
+    centroid_position=[0., 0., -20.e-6],
+    centroid_velocity=[0., 0., 2.e5*c],
+)
+beam_charge = -3.e-13
+beam_z0 = -90e-6
+beam_gaussian_distribution = picmi.GaussianBunchDistribution(
+    n_physical_particles= abs(beam_charge) / q_e,
+    rms_bunch_size=[0.5e-6, 0.5e-6, 1.e-6],
+    rms_velocity=[2.*c, 2.*c, 200.*c],
+    centroid_position=[0., 0., beam_z0],
+    centroid_velocity=[0., 0., 4000.*c],
+)
+
+dens = 1.e23
+lramp = 8.e-3
+plasma_distribution = picmi.AnalyticDistribution(
+    density_expression=
+        f'(z<lramp)*0.5*(1-cos(pi*z/lramp))*dens+(z>lramp)*dens',
+    pi=3.141592653589793,
+    dens=dens,
+    lramp = lramp,
+    lower_bound=[-70.e-6, -70.e-6, 0.],
+    upper_bound=[70.e-6, 70.e-6, 0.2],
+    fill_in=True)
+
+beam = picmi.Species(
+    particle_type='electron',
+    name='beam',
+    initial_distribution=beam_gaussian_distribution)
+driver = picmi.Species(
+    particle_type='electron',
+    name='driver',
+    initial_distribution=driver_gaussian_distribution)
+plasma_e = picmi.Species(
+    particle_type='electron',
+    name='plasma_e',
+    initial_distribution=plasma_distribution)
+plasma_i = picmi.Species(
+    particle_type='proton',
+    name='plasma_p',
+    initial_distribution=plasma_distribution)
+
+sim = picmi.Simulation(
+    solver = solver,
+    max_steps = max_steps,
+    verbose = 1,
+    warpx_use_filter = 0,
+    gamma_boost = gamma_boost,
+    warpx_use_fdtd_nci_corr = True,
+    particle_shape='cubic'
+)
+
+sim.add_species_through_plane(
+    beam,
+    layout=picmi.PseudoRandomLayout(
+        grid=grid,
+        n_macroparticles=1000
+    ),
+    injection_plane_position=0.,
+    injection_plane_normal_vector=[0.,0.,1.]
+)
+sim.add_species_through_plane(
+    driver,
+    layout=picmi.PseudoRandomLayout(
+        grid=grid,
+        n_macroparticles=1000
+    ),
+    injection_plane_position=0.,
+    injection_plane_normal_vector=[0.,0.,1.]
+)
+sim.add_species(
+    plasma_e,
+    layout=picmi.GriddedLayout(
+        grid=grid,
+        n_macroparticle_per_cell=[npc_x, npc_x, npc_z]
+    )
+)
+sim.add_species(
+    plasma_i,
+    layout=picmi.GriddedLayout(
+        grid=grid,
+        n_macroparticle_per_cell=[npc_x, npc_x, npc_z]
+    )
+)
+
+field_diag = picmi.FieldDiagnostic(
+    name = 'Python_PlasmaAcceleration_plt',
+    grid = grid,
+    period = max_steps,
+    data_list = ['Ex', 'Ey', 'Ez', 'Jx', 'Jy', 'Jz', 'part_per_cell'],
+    # write_dir='Python_PlasmaAcceleration_plt',
+    warpx_format='openpmd',
+)
+
+part_diag = picmi.ParticleDiagnostic(
+    name = 'Python_PlasmaAcceleration_plt',
+    period = max_steps,
+    species = [driver, beam, plasma_e],
+    # write_dir='Python_PlasmaAcceleration_plt',
+    warpx_format='openpmd',
+)
+
+beamrel_red_diag = picmi.ReducedDiagnostic(
+    diag_type='BeamRelevant',
+    name='beamrel',
+    species=beam,
+    period=10)
+
+sim.add_diagnostic(field_diag)
+sim.add_diagnostic(part_diag)
+sim.add_diagnostic(beamrel_red_diag)
+
+# write_inputs will create an inputs file that can be used to run
+# with the compiled version.
+sim.write_input_file(file_name = 'inputs_from_PICMI')
+
+# Alternatively, sim.step will run WarpX, controlling it from Python
+sim.step()
diff --git a/Examples/Physics_applications/plasma_acceleration/README.rst b/Examples/Physics_applications/plasma_acceleration/README.rst
index d5775e93aa8..0077a2c4252 100644
--- a/Examples/Physics_applications/plasma_acceleration/README.rst
+++ b/Examples/Physics_applications/plasma_acceleration/README.rst
@@ -11,12 +11,32 @@ On the other hand, for physical situation with strong asymmetries (e.g., non-rou
 
 Additionally, to speed up computation, this example uses the :ref:`boosted frame method <theory-boostedframe>` to effectively model long acceleration lengths.
 
-Alternatively, an other common approximation for PWFAs is quasi-static modeling, e.g., if effects such as self-injection can be ignored.
+Alternatively, another common approximation for PWFAs is quasi-static modeling, useful if effects such as self-injection can be ignored.
 In the Beam, Plasma & Accelerator Simulation Toolkit (BLAST), `HiPACE++ <https://hipace.readthedocs.io>`__ provides such methods.
 
 .. note::
 
-   TODO: The Python (PICMI) input file should use the boosted frame method, like the ``inputs_3d_boost`` file.
+    There are two FDTD Maxwell field solvers that compute the field push implemented in WarpX: the Yee and Cole-Karkkainen solver with Cowan coefficients (CKC) solvers.
+    The later includes a modification that allows the numerical dispersion of light in vacuum to be exact, and that is why we choose CKC for the example.
+
+Lorentz boosted frame
+---------------------
+
+    WarpX simulations can be done in the laboratory or a :ref:`Lorentz-boosted <theory-boostedframe>` frame.
+    In the laboratory frame, there is typically no need to model the plasma ions species, since they are mainly stationary during the short time scales associated with the motion of plasma electrons.
+    In the boosted frame, that argument is no longer valid, as ions have relativistic velocities.
+
+    To avoid spurious effects in the boosted frame, we want the longitudinal cell size to be larger than the transverse one.
+    Translating this condition to the cell transverse (:math:`d_{x}`) and longitudinal dimensions (:math:`d_{z}`) in the laboratory frame leads to: :math:`d_{x} > (d_{z} (1+\beta_{b}) \gamma_{b})`,
+    where :math:`\beta_{b}` is the boosted frame velocity in units of :math:`c`.
+
+    To determine the total number of time steps of the simulation, we could either set the ``<zmax_plasma_to_compute_max_step>`` parameter to the end of the plasma (:math:`z_{\textrm{end}}`), or compute it using:
+
+    * boosted frame edge of the simulation box, :math:`\textrm{corner} = l_{e}/ ((1-\beta_{b}) \gamma_{b})`
+    * time of interaction in the boosted frame, :math:`T = \frac{z_{\textrm{end}}/\gamma_{b}-\textrm{corner}}{c (1+\beta_{b})}`
+    * total number of iterations, :math:`i_{\textrm{max}} = T/dt`
+
+    where :math:`l_{e}` is the position of the left edge of the simulation box (in respect to propagation direction).
 
 
 Run
@@ -24,7 +44,7 @@ Run
 
 This example can be run **either** as:
 
-* **Python** script: ``python3 PICMI_inputs_plasma_acceleration.py`` or
+* **Python** script: ``python3 PICMI_inputs_plasma_acceleration_3d.py`` or
 * WarpX **executable** using an input file: ``warpx.3d inputs_3d_boost``
 
 For `MPI-parallel <https://www.mpi-forum.org>`__ runs, prefix these lines with ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system.
@@ -33,11 +53,7 @@ For `MPI-parallel <https://www.mpi-forum.org>`__ runs, prefix these lines with `
 
    .. tab-item:: Python: Script
 
-      .. note::
-
-         TODO: This input file should use the boosted frame method, like the ``inputs_3d_boost`` file.
-
-      .. literalinclude:: PICMI_inputs_plasma_acceleration.py
+      .. literalinclude:: PICMI_inputs_plasma_acceleration_3d.py
          :language: python3
          :caption: You can copy this file from ``Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py``.
 
@@ -47,17 +63,40 @@ For `MPI-parallel <https://www.mpi-forum.org>`__ runs, prefix these lines with `
          :language: ini
          :caption: You can copy this file from ``Examples/Physics_applications/plasma_acceleration/inputs_3d_boost``.
 
-Analyze
--------
+Visualize
+---------
+
+The plasma has a density profile that consists of a long up-ramp followed by a uniform plateau.
 
-.. note::
+With this configuration the driver expels the plasma electrons, creating a bubble depleted of plasma electrons.
+This bubble is the characteristic shape of the nonlinear 3D plasma wake.
+This creates an accelerating field :math:`E_z`, as can be seen in :numref:`fig_pwfa_wake`, and transverse focusing forces (not shown).
 
-   This section is TODO.
+.. _fig_pwfa_wake:
 
+.. figure:: https://user-images.githubusercontent.com/10621396/292573839-8d6d7296-48f8-4942-94bd-548683fd35fe.png
+   :alt: Electron beam, in black, accelerated by wakefield, in red-blue colormap, driven by drive electron beam, in purple. Plasma particles are plotted in green. Set max_steps to 800 in the input to generate data for this figure.
+   :width: 100%
 
-Visualize
----------
+   Electron beam, in black, accelerated by wakefield, in red-blue colormap, driven by drive electron beam, in purple. Plasma particles are plotted in green. Set ``max_steps`` to 800 in the input to generate data for this figure. Note that the :math:`E_z` field is plotted with an interpolation option and appears smoother than the numerical resolution would suggest.
 
-.. note::
+The beam increases in energy as it propagates in the wake, as shown in :numref:`fig_pwfa_energy`.
+There is a slight inflection point in the evolution of the beam energy, showing how the beam accelerates as the wake forms around the electron beam and then gains energy linearly in time once the wake stabilizes.
+Note that this is the mean energy in the boosted frame.
+:ref:`Lab Frame <python_lab_frame_diag>`  or :ref:`back-transformed <running-cpp-parameters-diagnostics-btd>` diagnostics can be used to analyze beam energy in the lab frame. See also the :ref:`FAQ <faq-btd>`.
+
+.. _fig_pwfa_energy:
+
+.. figure:: https://user-images.githubusercontent.com/10621396/290962801-97d994f9-d48d-4f76-a37e-d14f6781d680.png
+   :alt: Energy gain of electron beam in wake through this simulation (this is the mean beam energy in the boosted frame, NOT the lab frame). Set ``max_steps`` to 800 in the input to generate data for this figure.
+   :width: 80%
+
+   Energy gain of electron beam in wake through this simulation (this is the mean beam energy in the boosted frame, NOT the lab frame). Set ``max_steps`` to 800 in the input to generate data for this figure.
+
+The plots can be generated with the following script:
+
+.. dropdown:: Script ``plot_3d.py``
 
-   This section is TODO.
+   .. literalinclude:: plot_3d.py
+      :language: python3
+      :caption: You can copy this file from ``Examples/Physics_applications/plasma_acceleration/plot_3d.py``.
diff --git a/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost b/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost
index 76dcd3ee286..05b0fcaddaa 100644
--- a/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost
+++ b/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost
@@ -112,5 +112,6 @@ beam.rigid_advance = true
 
 # Diagnostics
 diagnostics.diags_names = diag1
+diag1.format = openpmd
 diag1.intervals = 500
 diag1.diag_type = Full
diff --git a/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost b/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost
index 2264872ec43..5e5f019b986 100644
--- a/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost
+++ b/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost
@@ -1,34 +1,37 @@
 #################################
 ####### GENERAL PARAMETERS ######
 #################################
-stop_time = 3.93151387287e-11
-amr.n_cell = 32 32 256
-amr.max_grid_size = 64
+max_step = 10
+amr.n_cell = 64 64 64
 amr.blocking_factor = 32
+amr.max_grid_size = 128
 amr.max_level = 0
 geometry.dims = 3
-geometry.prob_lo = -0.00015 -0.00015 -0.00012
-geometry.prob_hi = 0.00015 0.00015 1.e-06
+geometry.prob_lo = -100.e-6 -100.e-6 -149.e-6
+geometry.prob_hi =  100.e-6  100.e-6  1e-06
 
 #################################
 ####### Boundary condition ######
 #################################
-boundary.field_lo = periodic periodic pml
-boundary.field_hi = periodic periodic pml
+boundary.field_hi = pml pml pml
+boundary.field_lo = pml pml pml
+boundary.particle_hi = absorbing absorbing absorbing
+boundary.particle_lo = absorbing absorbing absorbing
 
 #################################
 ############ NUMERICS ###########
 #################################
-algo.maxwell_solver = ckc
+algo.maxwell_solver = CKC
 warpx.verbose = 1
-warpx.do_dive_cleaning = 0
 warpx.use_filter = 1
-warpx.cfl = .99
+warpx.pml_ncell = 10
+warpx.cfl = 0.99
 warpx.do_moving_window = 1
-warpx.moving_window_dir = z
-warpx.moving_window_v = 1. # in units of the speed of light
-my_constants.lramp = 8.e-3
-my_constants.dens  = 1e+23
+warpx.boost_direction = z
+warpx.moving_window_v = 1.0
+
+my_constants.dens = 1e+23
+my_constants.lramp = 0.008
 
 # Order of particle shape factors
 algo.particle_shape = 3
@@ -36,113 +39,113 @@ algo.particle_shape = 3
 #################################
 ######### BOOSTED FRAME #########
 #################################
-warpx.gamma_boost = 10.0
-warpx.boost_direction = z
+warpx.gamma_boost =       10.0
+warpx.moving_window_dir = z
 
 #################################
 ############ PLASMA #############
 #################################
-particles.species_names = driver plasma_e plasma_p beam driverback
+particles.species_names = beam driver plasma_e plasma_p
 particles.use_fdtd_nci_corr = 1
-particles.rigid_injected_species = driver beam
+particles.rigid_injected_species = beam driver
+
+beam.charge = -q_e
+beam.initialize_self_fields = 0
+beam.injection_style = gaussian_beam
+beam.mass = m_e
+beam.momentum_distribution_type = gaussian
+beam.npart = 1000
+beam.q_tot = -3e-13
+beam.rigid_advance = 1
+beam.ux_m = 0.0
+beam.ux_th = 2.0
+beam.uy_m = 0.0
+beam.uy_th = 2.0
+beam.uz_m = 4000.0
+beam.uz_th = 200.0
+beam.x_m = 0.0
+beam.x_rms = 5e-07
+beam.y_m = 0.0
+beam.y_rms = 5e-07
+beam.z_m = -9e-05
+beam.z_rms = 1e-06
+beam.zinject_plane = 0.0
 
 driver.charge = -q_e
-driver.mass = 1.e10
-driver.injection_style = "gaussian_beam"
-driver.x_rms = 2.e-6
-driver.y_rms = 2.e-6
-driver.z_rms = 4.e-6
-driver.x_m = 0.
-driver.y_m = 0.
-driver.z_m = -20.e-6
+driver.initialize_self_fields = 0
+driver.injection_style = gaussian_beam
+driver.mass = m_e
+driver.momentum_distribution_type = gaussian
 driver.npart = 1000
-driver.q_tot = -1.e-9
-driver.momentum_distribution_type = "gaussian"
+driver.q_tot = -1e-10
+driver.rigid_advance = 1
 driver.ux_m = 0.0
+driver.ux_th = 2.0
 driver.uy_m = 0.0
-driver.uz_m = 200000.
-driver.ux_th = 2.
-driver.uy_th = 2.
-driver.uz_th = 20000.
-driver.zinject_plane = 0.
-driver.rigid_advance = true
-
-driverback.charge = q_e
-driverback.mass = 1.e10
-driverback.injection_style = "gaussian_beam"
-driverback.x_rms = 2.e-6
-driverback.y_rms = 2.e-6
-driverback.z_rms = 4.e-6
-driverback.x_m = 0.
-driverback.y_m = 0.
-driverback.z_m = -20.e-6
-driverback.npart = 1000
-driverback.q_tot = 1.e-9
-driverback.momentum_distribution_type = "gaussian"
-driverback.ux_m = 0.0
-driverback.uy_m = 0.0
-driverback.uz_m = 200000.
-driverback.ux_th = 2.
-driverback.uy_th = 2.
-driverback.uz_th = 20000.
-driverback.do_backward_propagation = true
+driver.uy_th = 2.0
+driver.uz_m = 200000.0
+driver.uz_th = 20000.0
+driver.x_m = 0.0
+driver.x_rms = 2e-06
+driver.y_m = 0.0
+driver.y_rms = 2e-06
+driver.z_m = -2e-05
+driver.z_rms = 4e-06
+driver.zinject_plane = 0.0
 
 plasma_e.charge = -q_e
+plasma_e.density_function(x,y,z) = (z<lramp)*0.5*(1-cos(pi*z/lramp))*dens+(z>lramp)*dens
+plasma_e.do_continuous_injection = 1
+plasma_e.initialize_self_fields = 0
+plasma_e.injection_style = nuniformpercell
 plasma_e.mass = m_e
-plasma_e.injection_style = "NUniformPerCell"
-plasma_e.zmin = -100.e-6 # 0.e-6
-plasma_e.zmax = 0.2
-plasma_e.xmin = -70.e-6
-plasma_e.xmax =  70.e-6
-plasma_e.ymin = -70.e-6
-plasma_e.ymax =  70.e-6
-# plasma_e.profile = constant
-# plasma_e.density = 1.e23
-plasma_e.profile = parse_density_function
-plasma_e.density_function(x,y,z) = "(z<lramp)*0.5*(1-cos(pi*z/lramp))*dens+(z>lramp)*dens"
+plasma_e.momentum_distribution_type = constant
 plasma_e.num_particles_per_cell_each_dim = 1 1 1
-plasma_e.momentum_distribution_type = "at_rest"
-plasma_e.do_continuous_injection = 1
+plasma_e.profile = parse_density_function
+plasma_e.ux = 0.0
+plasma_e.uy = 0.0
+plasma_e.uz = 0.0
+plasma_e.xmax = 7e-05
+plasma_e.xmin = -7e-05
+plasma_e.ymax = 7e-05
+plasma_e.ymin = -7e-05
+plasma_e.zmax = 0.2
+plasma_e.zmin = 0.0
 
 plasma_p.charge = q_e
+plasma_p.density_function(x,y,z) = (z<lramp)*0.5*(1-cos(pi*z/lramp))*dens+(z>lramp)*dens
+plasma_p.do_continuous_injection = 1
+plasma_p.initialize_self_fields = 0
+plasma_p.injection_style = nuniformpercell
 plasma_p.mass = m_p
-plasma_p.injection_style = "NUniformPerCell"
-plasma_p.zmin = -100.e-6 # 0.e-6
-plasma_p.zmax = 0.2
-# plasma_p.profile = "constant"
-# plasma_p.density = 1.e23
-plasma_p.profile = parse_density_function
-plasma_p.density_function(x,y,z) = "(z<lramp)*0.5*(1-cos(pi*z/lramp))*dens+(z>lramp)*dens"
-plasma_p.xmin = -70.e-6
-plasma_p.xmax =  70.e-6
-plasma_p.ymin = -70.e-6
-plasma_p.ymax =  70.e-6
+plasma_p.momentum_distribution_type = constant
 plasma_p.num_particles_per_cell_each_dim = 1 1 1
-plasma_p.momentum_distribution_type = "at_rest"
-plasma_p.do_continuous_injection = 1
+plasma_p.profile = parse_density_function
+plasma_p.ux = 0.0
+plasma_p.uy = 0.0
+plasma_p.uz = 0.0
+plasma_p.xmax = 7e-05
+plasma_p.xmin = -7e-05
+plasma_p.ymax = 7e-05
+plasma_p.ymin = -7e-05
+plasma_p.zmax = 0.2
+plasma_p.zmin = 0.0
 
-beam.charge = -q_e
-beam.mass = m_e
-beam.injection_style = "gaussian_beam"
-beam.x_rms = .5e-6
-beam.y_rms = .5e-6
-beam.z_rms = 1.e-6
-beam.x_m = 0.
-beam.y_m = 0.
-beam.z_m = -100.e-6
-beam.npart = 1000
-beam.q_tot = -5.e-10
-beam.momentum_distribution_type = "gaussian"
-beam.ux_m = 0.0
-beam.uy_m = 0.0
-beam.uz_m = 2000.
-beam.ux_th = 2.
-beam.uy_th = 2.
-beam.uz_th = 200.
-beam.zinject_plane = .8e-3
-beam.rigid_advance = true
+#################################
+########## DIAGNOSTICS ##########
+#################################
 
-# Diagnostics
 diagnostics.diags_names = diag1
-diag1.intervals = 10000
+
 diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell
+diag1.file_prefix = diags/diag1
+diag1.format = openpmd
+diag1.intervals = 10
+diag1.species = driver beam plasma_e
+diag1.write_species = 1
+
+warpx.reduced_diags_names = beamrel
+beamrel.intervals = 1
+beamrel.species = beam
+beamrel.type = BeamRelevant
diff --git a/Examples/Physics_applications/plasma_acceleration/plot_3d.py b/Examples/Physics_applications/plasma_acceleration/plot_3d.py
new file mode 100644
index 00000000000..7da5efdfb34
--- /dev/null
+++ b/Examples/Physics_applications/plasma_acceleration/plot_3d.py
@@ -0,0 +1,78 @@
+#!/usr/bin/env python3
+
+# Copyright 2023 The WarpX Community
+#
+# This file is part of WarpX.
+#
+# Authors: Ryan Sandberg
+# License: BSD-3-Clause-LBNL
+#
+# This is a script plots the wakefield and energy gain of a PWFA simulation.
+
+
+from matplotlib import pyplot as plt
+
+plt.rcParams.update({'font.size':16})
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+
+ts = OpenPMDTimeSeries('diags/diag1')
+
+iteration=ts.iterations[-1]
+select = {'y':[-3e-6,3e-6]}
+
+pez, pex = ts.get_particle(species='plasma_e',
+                           iteration=iteration,
+                           var_list=['z','x'],
+                           select=select)
+pdz, pdx = ts.get_particle(species='driver',
+                           iteration=iteration,
+                           var_list=['z','x'])
+bz, bx = ts.get_particle(species='beam',
+                         iteration=iteration,
+                         var_list=['z','x'])
+
+fez, info = ts.get_field('E',
+                         coord='z',
+                         iteration=iteration,
+                         slice_across='y'
+                        )
+
+psize = 0.2
+xscale = 1e6
+zscale = 1e2
+ext = info.imshow_extent
+ext = [ext[2]*zscale,ext[3]*zscale,ext[0]*xscale,ext[1]*xscale]
+plt.figure()
+ezplt = plt.imshow(fez.T/1e9,
+                   extent=ext,
+                   origin='lower',
+                   aspect='auto',
+                   cmap='RdBu_r',
+                   vmin=-20,
+                   vmax=20,
+                   interpolation='sinc',
+                  )
+cb = plt.colorbar(ezplt,ax=plt.gca())
+cb.set_label(r'$E_z$ (GV/m)')
+plt.scatter(pez*zscale,pex*xscale,s=psize,c='g')
+plt.scatter(pdz*zscale,pdx*xscale,s=psize,c='tab:purple')
+plt.scatter(bz*zscale,bx*xscale,s=psize,c='k')
+plt.ylim(-60,60)
+
+plt.xlabel(r'Propagation distance z (cm)')
+plt.ylabel(r'Transverse direction x ($\mu$m)')
+plt.tight_layout()
+plt.savefig('pwfa_wakefield.png')
+plt.close()
+
+
+data = np.genfromtxt("diags/reducedfiles/beamrel.txt")
+plt.figure()
+c = 2.998e8
+plt.plot(c*data[:,1]*1e3,data[:,8]*0.511)
+plt.ylabel('Mean energy, boosted frame (MeV)')
+plt.xlabel(r'Propagation distance z (mm)')
+plt.tight_layout()
+plt.savefig('pwfa_beam_energy.png')
+plt.close()
diff --git a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json
index 2860e26ac8f..a9a6e8d8a30 100644
--- a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json
+++ b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json
@@ -1,45 +1,49 @@
 {
+  "lev=0": {
+    "Bx": 0.3939499531818974,
+    "By": 9563.615022656977,
+    "Bz": 0.3591030695783777,
+    "Ex": 1057402535431.6962,
+    "Ey": 144244335.46926767,
+    "Ez": 3010491709181.346,
+    "jx": 3227914372.5244455,
+    "jy": 5068791609.191265,
+    "jz": 102420312561226.16
+  },
   "beam": {
+    "particle_position_x": 0.0003765650634942415,
+    "particle_position_y": 0.0,
+    "particle_position_z": 2.017924300869033,
     "particle_momentum_x": 4.3319528267672493e-19,
     "particle_momentum_y": 4.501000169972693e-19,
     "particle_momentum_z": 2.730833780425186e-17,
-    "particle_position_x": 0.0003765650634942415,
-    "particle_position_y": 2.017924300869033,
     "particle_weight": 3744905444676.4585
   },
   "driver": {
-    "particle_momentum_x": 4.404628898485127e-19,
-    "particle_momentum_y": 4.375137239499453e-19,
-    "particle_momentum_z": 2.7424245305788556e-15,
-    "particle_position_x": 0.0016253419519751734,
-    "particle_position_y": 0.31953805221039305,
+    "particle_position_x": 0.0016253419519739634,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.31953805221039305,
+    "particle_momentum_x": 4.404628902111724e-19,
+    "particle_momentum_y": 4.375137239494432e-19,
+    "particle_momentum_z": 2.74242453056392e-15,
     "particle_weight": 93622636116911.45
   },
-  "lev=0": {
-    "Bx": 0.3939499531816364,
-    "By": 9563.615022656964,
-    "Bz": 0.3591030695780125,
-    "Ex": 1057402535431.6746,
-    "Ey": 144244335.46922195,
-    "Ez": 3010491709181.3467,
-    "jx": 3227914369.7171793,
-    "jy": 5068791609.12788,
-    "jz": 102420312561226.17
-  },
   "plasma_e": {
-    "particle_momentum_x": 1.2615128163589931e-23,
-    "particle_momentum_y": 1.5392234051130659e-27,
-    "particle_momentum_z": 1.3694904132251956e-18,
     "particle_position_x": 0.01771877489382093,
-    "particle_position_y": 0.020620065534187347,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.020620065534187347,
+    "particle_momentum_x": 1.2615128163630441e-23,
+    "particle_momentum_y": 1.5392234055850962e-27,
+    "particle_momentum_z": 1.3694904132251956e-18,
     "particle_weight": 727797380314273.1
   },
   "plasma_p": {
-    "particle_momentum_x": 1.2615099272124483e-23,
-    "particle_momentum_y": 1.53924926768808e-27,
-    "particle_momentum_z": 2.51631809113383e-15,
     "particle_position_x": 0.017718749986451513,
-    "particle_position_y": 0.020620065463180168,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.020620065463180168,
+    "particle_momentum_x": 1.2615099272164996e-23,
+    "particle_momentum_y": 1.5392492681601607e-27,
+    "particle_momentum_z": 2.51631809113383e-15,
     "particle_weight": 727797380314273.1
   }
 }
\ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json
index f26534c6cec..b4aaf46b17b 100644
--- a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json
+++ b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json
@@ -1,58 +1,38 @@
 {
   "beam": {
-    "particle_momentum_x": 4.1784825059093756e-19,
-    "particle_momentum_y": 4.56492260137707e-19,
-    "particle_momentum_z": 2.733972888170628e-17,
-    "particle_position_x": 0.0003995213395426269,
-    "particle_position_y": 0.0004148795632360405,
-    "particle_position_z": 1.971520197241896,
-    "particle_weight": 3120754537.230381
+    "particle_momentum_x": 4.278894587984017e-19,
+    "particle_momentum_y": 4.0908169483354626e-19,
+    "particle_momentum_z": 5.4603362113804313e-17,
+    "particle_position_x": 0.0004073990679742117,
+    "particle_position_y": 0.00040229749956423387,
+    "particle_position_z": 1.7797898827251948,
+    "particle_weight": 1872452.7223382287
   },
   "driver": {
-    "particle_momentum_x": 4.7004364050785617e+21,
-    "particle_momentum_y": 4.6785862113093076e+21,
-    "particle_momentum_z": 2.999596093018443e+25,
-    "particle_position_x": 0.0015811715645607214,
-    "particle_position_y": 0.0016212416120016558,
-    "particle_position_z": 0.3737871668388371,
-    "particle_weight": 6241509074.460762
-  },
-  "driverback": {
-    "particle_momentum_x": 4.813131349021333e+21,
-    "particle_momentum_y": 5.165480740901231e+21,
-    "particle_momentum_z": 3.005830430844926e+25,
-    "particle_position_x": 0.0016489932623583448,
-    "particle_position_y": 0.0016186267675597135,
-    "particle_position_z": 0.42037151381228377,
-    "particle_weight": 6241509074.460762
+    "particle_momentum_x": 4.2537746521796484e-19,
+    "particle_momentum_y": 4.265395425729062e-19,
+    "particle_momentum_z": 2.7386179521483695e-15,
+    "particle_position_x": 0.0015730820359599112,
+    "particle_position_y": 0.001567442518035302,
+    "particle_position_z": 0.3821097659069855,
+    "particle_weight": 624150907.4460764
   },
   "lev=0": {
-    "Bx": 44900.36241701639,
-    "By": 44925.14975609915,
-    "Bz": 6.522149557732535,
-    "Ex": 1740395711918.398,
-    "Ey": 1739383965652.072,
-    "Ez": 15904996412605.766,
-    "jx": 181761783700.78683,
-    "jy": 211922069601.8713,
-    "jz": 1808517835456823.8
+    "Ex": 40124341691.55219,
+    "Ey": 40231025043.63612,
+    "Ez": 767573724673.0515,
+    "jx": 1964457768.6389933,
+    "jy": 2135134278.0101168,
+    "jz": 131704307418024.33,
+    "part_per_cell": 5872.0
   },
   "plasma_e": {
-    "particle_momentum_x": 1.2522307177726677e-25,
-    "particle_momentum_y": 1.255266202984791e-25,
-    "particle_momentum_z": 9.78204821543727e-18,
-    "particle_position_x": 0.12656250028461832,
-    "particle_position_y": 0.12656250028543425,
-    "particle_position_z": 0.06789191622060947,
-    "particle_weight": 4834401345.035692
-  },
-  "plasma_p": {
-    "particle_momentum_x": 1.2522307202410732e-25,
-    "particle_momentum_y": 1.2552662045089724e-25,
-    "particle_momentum_z": 1.7961333879129066e-14,
-    "particle_position_x": 0.12656249999984498,
-    "particle_position_y": 0.12656249999984454,
-    "particle_position_z": 0.06789191621903365,
-    "particle_weight": 4834401345.035692
+    "particle_momentum_x": 7.643918230575959e-27,
+    "particle_momentum_y": 7.663034837733629e-27,
+    "particle_momentum_z": 5.260568129295649e-18,
+    "particle_position_x": 0.06655000001193831,
+    "particle_position_y": 0.06655000001284267,
+    "particle_position_z": 0.013805068230195833,
+    "particle_weight": 116268127.33666785
   }
 }
\ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json b/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json
index 5cad819a7a8..19b04454bc5 100644
--- a/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json
+++ b/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json
@@ -1,29 +1,38 @@
 {
+  "lev=0": {
+    "Ex": 85212422458.16203,
+    "Ey": 84961665803.53134,
+    "Ez": 769698287688.4983,
+    "jx": 677671298.7936147,
+    "jy": 698210615.5401165,
+    "jz": 65852153140593.5,
+    "part_per_cell": 2000.0
+  },
   "beam": {
-    "particle_momentum_x": 5.708397145757236e-20,
-    "particle_momentum_y": 5.708397145757229e-20,
-    "particle_momentum_z": 7.574215683313995e-19,
-    "particle_position_x": 0.009793902444696388,
-    "particle_position_y": 0.00979390244469639,
-    "particle_position_z": 0.1047802881481609,
-    "particle_weight": 7031249999.999999
+    "particle_position_x": 0.00039529289114018036,
+    "particle_position_y": 0.00040531040300041394,
+    "particle_position_z": 1.7639542317504893,
+    "particle_momentum_x": 4.281807877835991e-19,
+    "particle_momentum_y": 4.2619036979363644e-19,
+    "particle_momentum_z": 5.469934524328222e-17,
+    "particle_weight": 1872452.7223382287
   },
-  "lev=0": {
-    "Ex": 7074528555636.087,
-    "Ey": 7074528078305.145,
-    "Ez": 19091561484766.926,
-    "jx": 92775315721769.97,
-    "jy": 92775315721769.94,
-    "jz": 1267816446833707.5,
-    "part_per_cell": 590976.0
+  "driver": {
+    "particle_position_x": 0.0015980853581705075,
+    "particle_position_y": 0.001659518252944162,
+    "particle_position_z": 0.3687834893095841,
+    "particle_momentum_x": 4.1784825059093756e-19,
+    "particle_momentum_y": 4.564922601377071e-19,
+    "particle_momentum_z": 2.734589138454591e-15,
+    "particle_weight": 624150907.4460764
   },
-  "plasma": {
+  "plasma_e": {
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.0,
     "particle_momentum_x": 0.0,
     "particle_momentum_y": 0.0,
     "particle_momentum_z": 0.0,
-    "particle_position_x": 58.98240000000005,
-    "particle_position_y": 58.98239999999999,
-    "particle_position_z": 66.35520000000001,
-    "particle_weight": 360000000000.0
+    "particle_weight": 0.0
   }
 }
\ No newline at end of file
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index f00fb083cbc..117a0f28254 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -2770,7 +2770,8 @@ useOMP = 1
 numthreads = 1
 compileTest = 0
 doVis = 0
-analysisRoutine = Examples/analysis_default_regression.py
+outputFile = PlasmaAccelerationBoost2d_plt
+analysisRoutine = Examples/analysis_default_openpmd_regression.py
 
 [PlasmaAccelerationBoost3d]
 buildDir = .
@@ -2786,7 +2787,8 @@ useOMP = 1
 numthreads = 1
 compileTest = 0
 doVis = 0
-analysisRoutine = Examples/analysis_default_regression.py
+outputFile = PlasmaAccelerationBoost3d_plt
+analysisRoutine = Examples/analysis_default_openpmd_regression.py
 
 [PlasmaAccelerationBoost3d_hybrid]
 buildDir = .
@@ -3702,9 +3704,9 @@ particleTypes = electrons
 
 [Python_PlasmaAcceleration]
 buildDir = .
-inputFile = Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py
+inputFile = Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_3d.py
 runtime_params =
-customRunCmd = python3 PICMI_inputs_plasma_acceleration.py
+customRunCmd = python3 PICMI_inputs_plasma_acceleration_3d.py
 dim = 3
 addToCompileString = USE_PYTHON_MAIN=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_APP=OFF -DWarpX_PYTHON=ON
@@ -3718,7 +3720,8 @@ compileTest = 0
 doVis = 0
 compareParticles = 1
 particleTypes = beam
-analysisRoutine = Examples/analysis_default_regression.py
+outputFile = diags/Python_PlasmaAcceleration_plt
+analysisRoutine = Examples/analysis_default_openpmd_regression.py
 
 [Python_PlasmaAcceleration1d]
 buildDir = .