diff --git a/Docs/source/latex_theory/AMR/AMR.tex b/Docs/source/latex_theory/AMR/AMR.tex
index 0e9eace8390..3ca94d95436 100644
--- a/Docs/source/latex_theory/AMR/AMR.tex
+++ b/Docs/source/latex_theory/AMR/AMR.tex
@@ -22,7 +22,7 @@ \section{Mesh refinement}
 In addition, for some implementations where the field that is computed at a given level is affected by the solution at finer levels, there are cases where the procedure violates the integral of Gauss' Law around the refined patch, leading to long range errors \cite{Vaylpb2002,Colellajcp2010}. As will be shown below, in the procedure that has been developed in WarpX, the field at a given refinement level is not affected by the solution at finer levels, and is thus not affected by this type of error.
 
 \subsection{Electrostatic}
-A cornerstone of the Particle-In-Cell method is that assuming a particle lying in a hypothetical infinite grid, then if the grid is regular and symmetrical, and if the order of field gathering matches the order of charge (or current) deposition, then there is no self-force of the particle acting on itself: a) anywhere if using the so-called ``momentum conserving'' gathering scheme; b) on average within one cell if using the ``energy conserving'' gathering scheme \cite{Birdsalllangdon}. A breaking of the regularity and/or symmetry in the grid, whether it is from the use of irregular meshes or mesh refinement, and whether one uses finite difference, finite volume or finite elements, results in a net spurious self-force (which does not average to zero over one cell)  for a macroparticle close to the point of irregularity (mesh refinement interface for the current purpose) \cite{Vaylpb2002,Colellajcp2010}.
+A cornerstone of the Particle-In-Cell method is that, given a particle lying in a hypothetical infinite grid, if the grid is regular and symmetrical, and if the order of field gathering matches the order of charge (or current) deposition, then there is no self-force of the particle acting on itself: a) anywhere if using the so-called ``momentum conserving'' gathering scheme; b) on average within one cell if using the ``energy conserving'' gathering scheme \cite{Birdsalllangdon}. A breaking of the regularity and/or symmetry in the grid, whether it is from the use of irregular meshes or mesh refinement, and whether one uses finite difference, finite volume or finite elements, results in a net spurious self-force (which does not average to zero over one cell)  for a macroparticle close to the point of irregularity (mesh refinement interface for the current purpose) \cite{Vaylpb2002,Colellajcp2010}.
 
 A sketch of the implementation of mesh refinement in WarpX is given in Figure~\ref{fig:ESAMR} (left). Given the solution of the electric potential at a refinement level $L_n$, it is interpolated onto the boundaries of the grid patch(es) at the next refined level $L_{n+1}$. The electric potential is then computed at level $L_{n+1}$ by solving the Poisson equation. This procedure necessitates the knowledge of the charge density at every level of refinement. For efficiency, the macroparticle charge is deposited on the highest level patch that contains them, and the charge density of each patch is added recursively to lower levels, down to the lowest.
 
diff --git a/Docs/source/latex_theory/allbibs.bib b/Docs/source/latex_theory/allbibs.bib
index a3c355aae25..b3475c5a81b 100644
--- a/Docs/source/latex_theory/allbibs.bib
+++ b/Docs/source/latex_theory/allbibs.bib
@@ -2181,7 +2181,6 @@ @article{LehePRE2016
 volume = {94},
 year = {2016}
 }
-
 @book{godfrey1985iprop,
 author = {Godfrey, B. B.},
 publisher = {Defense Technical Information Center},
diff --git a/Docs/source/latex_theory/input_output/input_output.tex b/Docs/source/latex_theory/input_output/input_output.tex
index e013e236445..446f3f2b5d6 100644
--- a/Docs/source/latex_theory/input_output/input_output.tex
+++ b/Docs/source/latex_theory/input_output/input_output.tex
@@ -30,7 +30,7 @@ \subsection{Inputs and outputs in a boosted frame simulation}
 
 \subsubsection{Input in a boosted frame simulation}
 \paragraph{Particles - }
-Particles are launched through a plane using a technique that is generic and applies to Lorentz boosted frame simulations in general, including plasma acceleration, and is illustrated using the case of a positively charged particle beam propagating through a background of cold electrons in an assumed continuous transverse focusing system, leading to a well-known growing transverse ``electron cloud'' instability  \cite{Vayprl07}. In the laboratory frame, the electron background is initially at rest and a moving window is used to follow the beam progression. Traditionally, the beam macroparticles are initialized all at once in the window, while background electron macroparticles are created continuously in front of the beam on a plane that is perpendicular to the beam velocity. In a frame moving at some fraction of the beam velocity in the laboratory frame, the beam initial conditions at a given time in the calculation frame are generally unknown and one must initialize the beam differently. However, it can be taken advantage of the fact that the beam initial conditions are often known for a given plane in the laboratory, either directly, or via simple calculation or projection from the conditions at a given time in the labortory frame. Given the position and velocity $\{x,y,z,v_x,v_y,v_z\}$ for each beam macroparticle at time $t=0$ for a beam moving at the average velocity $v_b=\beta_b c$ (where $c$ is the speed of light) in the laboratory, and using the standard synchronization ($z=z'=0$ at $t=t'=0$) between the laboratory and the calculation frames, the procedure for transforming the beam quantities for injection in a boosted frame moving at velocity $\beta c$ in the laboratory is as follows (the superscript $'$ relates to quantities known in the boosted frame while the superscript $^*$ relates to quantities that are know at a given longitudinal position $z^*$ but different times of arrival):
+Particles are launched through a plane using a technique that is generic and applies to Lorentz boosted frame simulations in general, including plasma acceleration, and is illustrated using the case of a positively charged particle beam propagating through a background of cold electrons in an assumed continuous transverse focusing system, leading to a well-known growing transverse ``electron cloud'' instability \cite{Vayprl07}. In the laboratory frame, the electron background is initially at rest and a moving window is used to follow the beam progression. Traditionally, the beam macroparticles are initialized all at once in the window, while background electron macroparticles are created continuously in front of the beam on a plane that is perpendicular to the beam velocity. In a frame moving at some fraction of the beam velocity in the laboratory frame, the beam initial conditions at a given time in the calculation frame are generally unknown and one must initialize the beam differently. However, it can be taken advantage of the fact that the beam initial conditions are often known for a given plane in the laboratory, either directly, or via simple calculation or projection from the conditions at a given time in the labortory frame. Given the position and velocity $\{x,y,z,v_x,v_y,v_z\}$ for each beam macroparticle at time $t=0$ for a beam moving at the average velocity $v_b=\beta_b c$ (where $c$ is the speed of light) in the laboratory, and using the standard synchronization ($z=z'=0$ at $t=t'=0$) between the laboratory and the calculation frames, the procedure for transforming the beam quantities for injection in a boosted frame moving at velocity $\beta c$ in the laboratory is as follows (the superscript $'$ relates to quantities known in the boosted frame while the superscript $^*$ relates to quantities that are know at a given longitudinal position $z^*$ but different times of arrival):
 
 \begin{enumerate}
 \item project positions at $z^*=0$ assuming ballistic propagation
diff --git a/Docs/source/theory/amr.rst b/Docs/source/theory/amr.rst
index 730593beea4..d83b1d7db9d 100644
--- a/Docs/source/theory/amr.rst
+++ b/Docs/source/theory/amr.rst
@@ -3,158 +3,64 @@
 Mesh refinement
 ===============
 
-.. raw:: latex
-
-   \centering
+.. _fig_ESAMR:
 
 .. figure:: ICNSP_2011_Vay_fig1.png
    :alt: Sketches of the implementation of mesh refinement in WarpX with the electrostatic (left) and electromagnetic (right) solvers. In both cases, the charge/current from particles are deposited at the finest levels first, then interpolated recursively to coarser levels. In the electrostatic case, the potential is calculated first at the coarsest level :math:`L_0`, the solution interpolated to the boundaries of the refined patch :math:`r` at the next level :math:`L_{1}` and the potential calculated at :math:`L_1`. The procedure is repeated iteratively up to the highest level. In the electromagnetic case, the fields are computed independently on each grid and patch without interpolation at boundaries. Patches are terminated by absorbing layers (PML) to prevent the reflection of electromagnetic waves. Additional coarse patch :math:`c` and fine grid :math:`a` are needed so that the full solution is obtained by substitution on :math:`a` as :math:`F_{n+1}(a)=F_{n+1}(r)+I[F_n( s )-F_{n+1}( c )]` where :math:`F` is the field, and :math:`I` is a coarse-to-fine interpolation operator. In both cases, the field solution at a given level :math:`L_n` is unaffected by the solution at higher levels :math:`L_{n+1}` and up, allowing for mitigation of some spurious effects (see text) by providing a transition zone via extension of the patches by a few cells beyond the desired refined area (red & orange rectangles) in which the field is interpolated onto particles from the coarser parent level only.
-   :name: fig:ESAMR
-   :width: 15cm
+   :width: 95%
 
    Sketches of the implementation of mesh refinement in WarpX with the electrostatic (left) and electromagnetic (right) solvers. In both cases, the charge/current from particles are deposited at the finest levels first, then interpolated recursively to coarser levels. In the electrostatic case, the potential is calculated first at the coarsest level :math:`L_0`, the solution interpolated to the boundaries of the refined patch :math:`r` at the next level :math:`L_{1}` and the potential calculated at :math:`L_1`. The procedure is repeated iteratively up to the highest level. In the electromagnetic case, the fields are computed independently on each grid and patch without interpolation at boundaries. Patches are terminated by absorbing layers (PML) to prevent the reflection of electromagnetic waves. Additional coarse patch :math:`c` and fine grid :math:`a` are needed so that the full solution is obtained by substitution on :math:`a` as :math:`F_{n+1}(a)=F_{n+1}(r)+I[F_n( s )-F_{n+1}( c )]` where :math:`F` is the field, and :math:`I` is a coarse-to-fine interpolation operator. In both cases, the field solution at a given level :math:`L_n` is unaffected by the solution at higher levels :math:`L_{n+1}` and up, allowing for mitigation of some spurious effects (see text) by providing a transition zone via extension of the patches by a few cells beyond the desired refined area (red & orange rectangles) in which the field is interpolated onto particles from the coarser parent level only.
 
-The mesh refinement methods that have been implemented in WarpX were developed following the following principles: i) avoidance of spurious effects from mesh refinement, or minimization of such effects; ii) user controllability of the spurious effects’ relative magnitude; iii) simplicity of implementation. The two main generic issues that were identified are: a) spurious self-force on macroparticles close to the mesh refinement interface (J. Vay et al. 2002; Colella and Norgaard 2010); b) reflection (and possible amplification) of short wavelength electromagnetic waves at the mesh refinement interface (Vay 2001). The two effects are due to the loss of translation invariance introduced by the asymmetry of the grid on each side of the mesh refinement interface.
+The mesh refinement methods that have been implemented in WarpX were developed following the following principles: i) avoidance of spurious effects from mesh refinement, or minimization of such effects; ii) user controllability of the spurious effects’ relative magnitude; iii) simplicity of implementation. The two main generic issues that were identified are: a) spurious self-force on macroparticles close to the mesh refinement interface :cite:p:`amr-Vaylpb2002,amr-Colellajcp2010`; b) reflection (and possible amplification) of short wavelength electromagnetic waves at the mesh refinement interface :cite:p:`amr-Vayjcp01`. The two effects are due to the loss of translation invariance introduced by the asymmetry of the grid on each side of the mesh refinement interface.
 
-In addition, for some implementations where the field that is computed at a given level is affected by the solution at finer levels, there are cases where the procedure violates the integral of Gauss’ Law around the refined patch, leading to long range errors (J. Vay et al. 2002; Colella and Norgaard 2010). As will be shown below, in the procedure that has been developed in WarpX, the field at a given refinement level is not affected by the solution at finer levels, and is thus not affected by this type of error.
+In addition, for some implementations where the field that is computed at a given level is affected by the solution at finer levels, there are cases where the procedure violates the integral of Gauss’ Law around the refined patch, leading to long range errors :cite:p:`amr-Vaylpb2002,amr-Colellajcp2010`. As will be shown below, in the procedure that has been developed in WarpX, the field at a given refinement level is not affected by the solution at finer levels, and is thus not affected by this type of error.
 
 Electrostatic
 -------------
 
-A cornerstone of the Particle-In-Cell method is that assuming a particle lying in a hypothetical infinite grid, then if the grid is regular and symmetrical, and if the order of field gathering matches the order of charge (or current) deposition, then there is no self-force of the particle acting on itself: a) anywhere if using the so-called “momentum conserving” gathering scheme; b) on average within one cell if using the “energy conserving” gathering scheme (Birdsall and Langdon 1991). A breaking of the regularity and/or symmetry in the grid, whether it is from the use of irregular meshes or mesh refinement, and whether one uses finite difference, finite volume or finite elements, results in a net spurious self-force (which does not average to zero over one cell) for a macroparticle close to the point of irregularity (mesh refinement interface for the current purpose) (J. Vay et al. 2002; Colella and Norgaard 2010).
-
-A sketch of the implementation of mesh refinement in WarpX is given in Figure \ `[fig:ESAMR] <#fig:ESAMR>`__ (left). Given the solution of the electric potential at a refinement level :math:`L_n`, it is interpolated onto the boundaries of the grid patch(es) at the next refined level :math:`L_{n+1}`. The electric potential is then computed at level :math:`L_{n+1}` by solving the Poisson equation. This procedure necessitates the knowledge of the charge density at every level of refinement. For efficiency, the macroparticle charge is deposited on the highest level patch that contains them, and the charge density of each patch is added recursively to lower levels, down to the lowest.
+A cornerstone of the Particle-In-Cell method is that given a particle lying in a hypothetical infinite grid, if the grid is regular and symmetrical, and if the order of field gathering matches the order of charge (or current) deposition, then there is no self-force of the particle acting on itself: a) anywhere if using the so-called “momentum conserving” gathering scheme; b) on average within one cell if using the “energy conserving” gathering scheme :cite:p:`amr-Birdsalllangdon`. A breaking of the regularity and/or symmetry in the grid, whether it is from the use of irregular meshes or mesh refinement, and whether one uses finite difference, finite volume or finite elements, results in a net spurious self-force (which does not average to zero over one cell) for a macroparticle close to the point of irregularity (mesh refinement interface for the current purpose) :cite:p:`amr-Vaylpb2002,amr-Colellajcp2010`.
 
-.. raw:: latex
+A sketch of the implementation of mesh refinement in WarpX is given in :numref:`fig_ESAMR`. Given the solution of the electric potential at a refinement level :math:`L_n`, it is interpolated onto the boundaries of the grid patch(es) at the next refined level :math:`L_{n+1}`. The electric potential is then computed at level :math:`L_{n+1}` by solving the Poisson equation. This procedure necessitates the knowledge of the charge density at every level of refinement. For efficiency, the macroparticle charge is deposited on the highest level patch that contains them, and the charge density of each patch is added recursively to lower levels, down to the lowest.
 
-   \centering
+.. _fig_ESselfforce:
 
 .. figure:: ICNSP_2011_Vay_fig2.png
    :alt: Position history of one charged particle attracted by its image induced by a nearby metallic (dirichlet) boundary. The particle is initialized at rest. Without refinement patch (reference case), the particle is accelerated by its image, is reflected specularly at the wall, then decelerates until it reaches its initial position at rest. If the particle is initialized inside a refinement patch, the particle is initially accelerated toward the wall but is spuriously reflected before it reaches the boundary of the patch whether using the method implemented in WarpX or the MC method. Providing a surrounding transition region 2 or 4 cells wide in which the potential is interpolated from the parent coarse solution reduces significantly the effect of the spurious self-force.
-   :name: fig:ESselfforce
-   :width: 15cm
+   :width: 95%
 
    Position history of one charged particle attracted by its image induced by a nearby metallic (dirichlet) boundary. The particle is initialized at rest. Without refinement patch (reference case), the particle is accelerated by its image, is reflected specularly at the wall, then decelerates until it reaches its initial position at rest. If the particle is initialized inside a refinement patch, the particle is initially accelerated toward the wall but is spuriously reflected before it reaches the boundary of the patch whether using the method implemented in WarpX or the MC method. Providing a surrounding transition region 2 or 4 cells wide in which the potential is interpolated from the parent coarse solution reduces significantly the effect of the spurious self-force.
 
-The presence of the self-force is illustrated on a simple test case that was introduced in (J. Vay et al. 2002) and also used in (Colella and Norgaard 2010): a single macroparticle is initialized at rest within a single refinement patch four cells away from the patch refinement boundary. The patch at level :math:`L_1` has :math:`32\times32` cells and is centered relative to the lowest :math:`64\times64` grid at level :math:`L_0` (“main grid”), while the macroparticle is centered in one direction but not in the other. The boundaries of the main grid are perfectly conducting, so that the macroparticle is attracted to the closest wall by its image. Specular reflection is applied when the particle reaches the boundary so that the motion is cyclic. The test was performed with WarpX using either linear or quadratic interpolation when gathering the main grid solution onto the refined patch boundary. It was also performed using another method from P. McCorquodale et al (labeled “MC” in this paper) based on the algorithm given in (Mccorquodale et al. 2004), which employs a more elaborate procedure involving two-ways interpolations between the main grid and the refined patch. A reference case was also run using a single :math:`128\times128` grid with no refined patch, in which it is observed that the particle propagates toward the closest boundary at an accelerated pace, is reflected specularly at the boundary, then slows down until it reaches its initial position at zero velocity. The particle position histories are shown for the various cases in Fig. `[fig:ESselfforce] <#fig:ESselfforce>`__. In all the cases using the refinement patch, the particle was spuriously reflected near the patch boundary and was effectively trapped in the patch. We notice that linear interpolation performs better than quadratic, and that the simple method implemented in WarpX performs better than the other proposed method for this test (see discussion below).
-
-.. raw:: latex
+The presence of the self-force is illustrated on a simple test case that was introduced in :cite:t:`amr-Vaylpb2002` and also used in :cite:t:`amr-Colellajcp2010`: a single macroparticle is initialized at rest within a single refinement patch four cells away from the patch refinement boundary. The patch at level :math:`L_1` has :math:`32\times32` cells and is centered relative to the lowest :math:`64\times64` grid at level :math:`L_0` ("main grid"), while the macroparticle is centered in one direction but not in the other. The boundaries of the main grid are perfectly conducting, so that the macroparticle is attracted to the closest wall by its image. Specular reflection is applied when the particle reaches the boundary so that the motion is cyclic. The test was performed with WarpX using either linear or quadratic interpolation when gathering the main grid solution onto the refined patch boundary. It was also performed using another method from P. McCorquodale et al (labeled "MC" in this paper) based on the algorithm given in :cite:t:`amr-Mccorquodalejcp2004`, which employs a more elaborate procedure involving two-ways interpolations between the main grid and the refined patch. A reference case was also run using a single :math:`128\times128` grid with no refined patch, in which it is observed that the particle propagates toward the closest boundary at an accelerated pace, is reflected specularly at the boundary, then slows down until it reaches its initial position at zero velocity. The particle position histories are shown for the various cases in :numref:`fig_ESselfforce`. In all the cases using the refinement patch, the particle was spuriously reflected near the patch boundary and was effectively trapped in the patch. We notice that linear interpolation performs better than quadratic, and that the simple method implemented in WarpX performs better than the other proposed method for this test (see discussion below).
 
-   \centering
+.. _fig_ESselfforcemap:
 
 .. figure:: ICNSP_2011_Vay_fig3.png
    :alt: (left) Maps of the magnitude of the spurious self-force :math:`\epsilon` in arbitrary units within one quarter of the refined patch, defined as :math:`\epsilon=\sqrt{(E_x-E_x^{ref})^2+(E_y-E_y^{ref})^2}`, where :math:`E_x` and :math:`E_y` are the electric field components within the patch experienced by one particle at a given location and :math:`E_x^{ref}` and :math:`E_y^{ref}` are the electric field from a reference solution. The map is given for the WarpX and the MC mesh refinement algorithms and for linear and quadratic interpolation at the patch refinement boundary. (right) Lineouts of the maximum (taken over neighboring cells) of the spurious self-force. Close to the interface boundary (x=0), the spurious self-force decreases at a rate close to one order of magnitude per cell (red line), then at about one order of magnitude per six cells (green line).
-   :name: fig:ESselfforcemap
-   :width: 15cm
+   :width: 95%
 
    (left) Maps of the magnitude of the spurious self-force :math:`\epsilon` in arbitrary units within one quarter of the refined patch, defined as :math:`\epsilon=\sqrt{(E_x-E_x^{ref})^2+(E_y-E_y^{ref})^2}`, where :math:`E_x` and :math:`E_y` are the electric field components within the patch experienced by one particle at a given location and :math:`E_x^{ref}` and :math:`E_y^{ref}` are the electric field from a reference solution. The map is given for the WarpX and the MC mesh refinement algorithms and for linear and quadratic interpolation at the patch refinement boundary. (right) Lineouts of the maximum (taken over neighboring cells) of the spurious self-force. Close to the interface boundary (x=0), the spurious self-force decreases at a rate close to one order of magnitude per cell (red line), then at about one order of magnitude per six cells (green line).
 
-The magnitude of the spurious self-force as a function of the macroparticle position was mapped and is shown in Fig. `[fig:ESselfforcemap] <#fig:ESselfforcemap>`__ for the WarpX and MC algorithms using linear or quadratic interpolations between grid levels. It is observed that the magnitude of the spurious self-force decreases rapidly with the distance between the particle and the refined patch boundary, at a rate approaching one order of magnitude per cell for the four cells closest to the boundary and about one order of magnitude per six cells beyond. The method implemented in WarpX offers a weaker spurious force on average and especially at the cells that are the closest to the coarse-fine interface where it is the largest and thus matters most.
+The magnitude of the spurious self-force as a function of the macroparticle position was mapped and is shown in :numref:`fig_ESselfforcemap` for the WarpX and MC algorithms using linear or quadratic interpolations between grid levels. It is observed that the magnitude of the spurious self-force decreases rapidly with the distance between the particle and the refined patch boundary, at a rate approaching one order of magnitude per cell for the four cells closest to the boundary and about one order of magnitude per six cells beyond. The method implemented in WarpX offers a weaker spurious force on average and especially at the cells that are the closest to the coarse-fine interface where it is the largest and thus matters most.
 We notice that the magnitude of the spurious self-force depends strongly on the distance to the edge of the patch and to the nodes of the underlying coarse grid, but weakly on the order of deposition and size of the patch.
 
-A method was devised and implemented in WarpX for reducing the magnitude of spurious self-forces near the coarse-fine boundaries as follows. Noting that the coarse grid solution is unaffected by the presence of the patch and is thus free of self-force, extra “transition” cells are added around the “effective” refined area.
-Within the effective area, the particles gather the potential in the fine grid. In the extra transition cells surrounding the refinement patch, the force is gathered directly from the coarse grid (an option, which has not yet been implemented, would be to interpolate between the coarse and fine grid field solutions within the transition zone so as to provide continuity of the force experienced by the particles at the interface). The number of cells allocated in the transition zones is controllable by the user in WarpX, giving the opportunity to check whether the spurious self-force is affecting the calculation by repeating it using different thicknesses of the transition zones. The control of the spurious force using the transition zone is illustrated in Fig. \ `[fig:ESselfforce] <#fig:ESselfforce>`__, where the calculation with WarpX using linear interpolation at the patch interface was repeated using either two or four cells transition regions (measured in refined patch cell units). Using two extra cells allowed for the particle to be free of spurious trapping within the refined area and follow a trajectory that is close to the reference one, and using four extra cells improved further to the point where the resulting trajectory becomes indistinguishable from the reference one.
-We note that an alternative method was devised for reducing the magnitude of self-force near the coarse-fine boundaries for the MC method, by using a special deposition procedure near the interface (Colella and Norgaard 2010).
+A method was devised and implemented in WarpX for reducing the magnitude of spurious self-forces near the coarse-fine boundaries as follows. Noting that the coarse grid solution is unaffected by the presence of the patch and is thus free of self-force, extra "transition" cells are added around the "effective" refined area.
+Within the effective area, the particles gather the potential in the fine grid. In the extra transition cells surrounding the refinement patch, the force is gathered directly from the coarse grid (an option, which has not yet been implemented, would be to interpolate between the coarse and fine grid field solutions within the transition zone so as to provide continuity of the force experienced by the particles at the interface). The number of cells allocated in the transition zones is controllable by the user in WarpX, giving the opportunity to check whether the spurious self-force is affecting the calculation by repeating it using different thicknesses of the transition zones. The control of the spurious force using the transition zone is illustrated in :numref:`fig_ESselfforce`, where the calculation with WarpX using linear interpolation at the patch interface was repeated using either two or four cells transition regions (measured in refined patch cell units). Using two extra cells allowed for the particle to be free of spurious trapping within the refined area and follow a trajectory that is close to the reference one, and using four extra cells improved further to the point where the resulting trajectory becomes indistinguishable from the reference one.
+We note that an alternative method was devised for reducing the magnitude of self-force near the coarse-fine boundaries for the MC method, by using a special deposition procedure near the interface :cite:p:`amr-Colellajcp2010`.
 
 Electromagnetic
 ---------------
 
-The method that is used for electrostatic mesh refinement is not directly applicable to electromagnetic calculations. As was shown in section 3.4 of (Vay 2001), refinement schemes relying solely on interpolation between coarse and fine patches lead to the reflection with amplification of the short wavelength modes that fall below the cutoff of the Nyquist frequency of the coarse grid. Unless these modes are damped heavily or prevented from occurring at their source, they may affect particle motion and their effect can escalate if trapped within a patch, via multiple successive reflections with amplification.
+The method that is used for electrostatic mesh refinement is not directly applicable to electromagnetic calculations. As was shown in section 3.4 of :cite:t:`amr-Vayjcp01`, refinement schemes relying solely on interpolation between coarse and fine patches lead to the reflection with amplification of the short wavelength modes that fall below the cutoff of the Nyquist frequency of the coarse grid. Unless these modes are damped heavily or prevented from occurring at their source, they may affect particle motion and their effect can escalate if trapped within a patch, via multiple successive reflections with amplification.
 
-To circumvent this issue, an additional coarse patch (with the same resolution as the parent grid) is added, as shown in Fig. \ `[fig:ESAMR] <#fig:ESAMR>`__-right and described in (Vay, Adam, and Heron 2004). Both the fine and the coarse grid patches are terminated by Perfectly Matched Layers, reducing wave reflection by orders of magnitude, controllable by the user (Berenger 1996; J.-L. Vay 2002). The source current resulting from the motion of charged macroparticles within the refined region is accumulated on the fine patch and is then interpolated onto the coarse patch and added onto the parent grid. The process is repeated recursively from the finest level down to the coarsest. The Maxwell equations are then solved for one time interval on the entire set of grids, by default for one time step using the time step of the finest grid. The field on the coarse and fine patches only contain the contributions from the particles that have evolved within the refined area but not from the current sources outside the area. The total contribution of the field from sources within and outside the refined area is obtained by adding the field from the refined grid :math:`F(r)`, and adding an interpolation :math:`I` of the difference between the relevant subset :math:`s` of the field in the parent grid :math:`F(s)` and the field of the coarse grid :math:`F( c )`, on an auxiliary grid :math:`a`, i.e. :math:`F(a)=F(r)+I[F(s)-F( c )]`. The field on the parent grid subset :math:`F(s)` contains contributions from sources from both within and outside of the refined area. Thus, in effect, there is substitution of the coarse field resulting from sources within the patch area by its fine resolution counterpart. The operation is carried out recursively starting at the coarsest level up to the finest.
+To circumvent this issue, an additional coarse patch (with the same resolution as the parent grid) is added, as shown in :numref:`fig_ESAMR` and described in :cite:t:`amr-Vaycpc04`. Both the fine and the coarse grid patches are terminated by Perfectly Matched Layers, reducing wave reflection by orders of magnitude, controllable by the user :cite:p:`amr-Berengerjcp96,amr-Vayjcp02`. The source current resulting from the motion of charged macroparticles within the refined region is accumulated on the fine patch and is then interpolated onto the coarse patch and added onto the parent grid. The process is repeated recursively from the finest level down to the coarsest. The Maxwell equations are then solved for one time interval on the entire set of grids, by default for one time step using the time step of the finest grid. The field on the coarse and fine patches only contain the contributions from the particles that have evolved within the refined area but not from the current sources outside the area. The total contribution of the field from sources within and outside the refined area is obtained by adding the field from the refined grid :math:`F(r)`, and adding an interpolation :math:`I` of the difference between the relevant subset :math:`s` of the field in the parent grid :math:`F(s)` and the field of the coarse grid :math:`F( c )`, on an auxiliary grid :math:`a`, i.e. :math:`F(a)=F(r)+I[F(s)-F( c )]`. The field on the parent grid subset :math:`F(s)` contains contributions from sources from both within and outside of the refined area. Thus, in effect, there is substitution of the coarse field resulting from sources within the patch area by its fine resolution counterpart. The operation is carried out recursively starting at the coarsest level up to the finest.
 An option has been implemented in which various grid levels are pushed with different time steps, given as a fixed fraction of the individual grid Courant conditions (assuming same cell aspect ratio for all grids and refinement by integer factors). In this case, the fields from the coarse levels, which are advanced less often, are interpolated in time.
 
 The substitution method has two potential drawbacks due to the inexact cancellation between the coarse and fine patches of : (i) the remnants of ghost fixed charges created by the particles entering and leaving the patches (this effect is due to the use of the electromagnetic solver and is different from the spurious self-force that was described for the electrostatic case); (ii) if using a Maxwell solver with a low-order stencil, the electromagnetic waves traveling on each patch at slightly different velocity due to numerical dispersion.
 The first issue results in an effective spurious multipole field whose magnitude decreases very rapidly with the distance to the patch boundary, similarly to the spurious self-force in the electrostatic case. Hence, adding a few extra transition cells surrounding the patches mitigates this effect very effectively.
-The tunability of WarpX’s electromagnetic finite-difference and pseudo-spectral solvers provides the means to optimize the numerical dispersion so as to minimize the second effect for a given application, which has been demonstrated on the laser-plasma interaction test case presented in (Vay, Adam, and Heron 2004).
-Both effects and their mitigation are described in more detail in (Vay, Adam, and Heron 2004).
+The tunability of WarpX’s electromagnetic finite-difference and pseudo-spectral solvers provides the means to optimize the numerical dispersion so as to minimize the second effect for a given application, which has been demonstrated on the laser-plasma interaction test case presented in :cite:t:`amr-Vaycpc04`.
+Both effects and their mitigation are described in more detail in :cite:t:`amr-Vaycpc04`.
 
 Caustics are supported anywhere on the grid with an accuracy that is set by the local resolution, and will be adequately resolved if the grid resolution supports the necessary modes from their sources to the points of wavefront crossing. The mesh refinement method that is implemented in WarpX has the potential to provide higher efficiency than the standard use of fixed gridding, by offering a path toward adaptive gridding following wavefronts.
 
-.. raw:: html
-
-   <div id="refs" class="references">
-
-.. raw:: html
-
-   <div id="ref-Berengerjcp96">
-
-Berenger, Jp. 1996. “Three-Dimensional Perfectly Matched Layer for the Absorption of Electromagnetic Waves.” *Journal of Computational Physics* 127 (2): 363–79.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Birdsalllangdon">
-
-Birdsall, C K, and A B Langdon. 1991. *Plasma Physics via Computer Simulation*. Adam-Hilger.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Colellajcp2010">
-
-Colella, Phillip, and Peter C Norgaard. 2010. “Controlling Self-Force Errors at Refinement Boundaries for Amr-Pic.” *Journal of Computational Physics* 229 (4): 947–57. https://doi.org/10.1016/J.Jcp.2009.07.004.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Mccorquodalejcp2004">
-
-Mccorquodale, P, P Colella, Dp Grote, and Jl Vay. 2004. “A Node-Centered Local Refinement Algorithm For Poisson’s Equation In Complex Geometries.” *Journal of Computational Physics* 201 (1): 34–60. https://doi.org/10.1016/J.Jcp.2004.04.022.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Vayjcp01">
-
-Vay, J.-L. 2001. “An Extended Fdtd Scheme for the Wave Equation: Application to Multiscale Electromagnetic Simulation.” *Journal of Computational Physics* 167 (1): 72–98.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Vayjcp02">
-
-———. 2002. “Asymmetric Perfectly Matched Layer for the Absorption of Waves.” *Journal of Computational Physics* 183 (2): 367–99. https://doi.org/10.1006/Jcph.2002.7175.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Vaycpc04">
-
-Vay, J.-L., J.-C. Adam, and A Heron. 2004. “Asymmetric Pml for the Absorption of Waves. Application to Mesh Refinement in Electromagnetic Particle-in-Cell Plasma Simulations.” *Computer Physics Communications* 164 (1-3): 171–77. https://doi.org/10.1016/J.Cpc.2004.06.026.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Vaylpb2002">
-
-Vay, Jl, P Colella, P Mccorquodale, B Van Straalen, A Friedman, and Dp Grote. 2002. “Mesh Refinement for Particle-in-Cell Plasma Simulations: Applications to and Benefits for Heavy Ion Fusion.” *Laser and Particle Beams* 20 (4): 569–75. https://doi.org/10.1017/S0263034602204139.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   </div>
+.. bibliography::
+   :keyprefix: amr-
diff --git a/Docs/source/theory/input_output.rst b/Docs/source/theory/input_output.rst
index 75427ad9c6d..21a5f5c8d2c 100644
--- a/Docs/source/theory/input_output.rst
+++ b/Docs/source/theory/input_output.rst
@@ -3,7 +3,7 @@
 Inputs and Outputs
 ==================
 
-Initialization of the plasma columns and drivers (laser or particle beam) is performed via the specification of multidimensional functions that describe the initial state with, if needed, a time dependence, or from reconstruction of distributions based on experimental data. Care is needed when initializing quantities in parallel to avoid double counting and ensure smoothness of the distributions at the interface of computational domains. When the sum of the initial distributions of charged particles is not charge neutral, initial fields are computed using generally a static approximation with Poisson solves accompanied by proper relativistic scalings (Vay 2008; Cowan et al. 2013).
+Initialization of the plasma columns and drivers (laser or particle beam) is performed via the specification of multidimensional functions that describe the initial state with, if needed, a time dependence, or from reconstruction of distributions based on experimental data. Care is needed when initializing quantities in parallel to avoid double counting and ensure smoothness of the distributions at the interface of computational domains. When the sum of the initial distributions of charged particles is not charge neutral, initial fields are computed using generally a static approximation with Poisson solves accompanied by proper relativistic scalings :cite:p:`io-Vaypop2008, io-CowanPRSTAB13`.
 
 Outputs include dumps of particle and field quantities at regular intervals, histories of particle distributions moments, spectra, etc, and plots of the various quantities. In parallel simulations, the diagnostic subroutines need to handle additional complexity from the domain decomposition, as well as large amount of data that may necessitate data reduction in some form before saving to disk.
 
@@ -12,16 +12,17 @@ Simulations in a Lorentz boosted frame require additional considerations, as des
 Inputs and outputs in a boosted frame simulation
 ------------------------------------------------
 
-.. _Fig_inputoutput:
+.. _fig_inputoutput:
+
 .. figure:: Input_output.png
    :alt: (top) Snapshot of a particle beam showing “frozen" (grey spheres) and “active" (colored spheres) macroparticles traversing the injection plane (red rectangle). (bottom) Snapshot of the beam macroparticles (colored spheres) passing through the background of electrons (dark brown streamlines) and the diagnostic stations (red rectangles). The electrons, the injection plane and the diagnostic stations are fixed in the laboratory plane, and are thus counter-propagating to the beam in a boosted frame.
-   :width: 120mm
+   :width: 100%
 
    (top) Snapshot of a particle beam showing “frozen" (grey spheres) and “active" (colored spheres) macroparticles traversing the injection plane (red rectangle). (bottom) Snapshot of the beam macroparticles (colored spheres) passing through the background of electrons (dark brown streamlines) and the diagnostic stations (red rectangles). The electrons, the injection plane and the diagnostic stations are fixed in the laboratory plane, and are thus counter-propagating to the beam in a boosted frame.
 
 The input and output data are often known from, or compared to, experimental data. Thus, calculating in
 a frame other than the laboratory entails transformations of the data between the calculation frame and the laboratory
-frame. This section describes the procedures that have been implemented in the Particle-In-Cell framework Warp (Grote et al. 2005) to handle the input and output of data between the frame of calculation and the laboratory frame (Vay et al. 2011). Simultaneity of events between two frames is valid only for a plane that is perpendicular to the relative motion of the frame. As a result, the input/output processes involve the input of data (particles or fields) through a plane, as well as output through a series of planes, all of which are perpendicular to the direction of the relative velocity between the frame of calculation and the other frame of choice.
+frame. This section describes the procedures that have been implemented in the Particle-In-Cell framework Warp :cite:p:`io-Warp` to handle the input and output of data between the frame of calculation and the laboratory frame :cite:p:`io-Vaypop2011`. Simultaneity of events between two frames is valid only for a plane that is perpendicular to the relative motion of the frame. As a result, the input/output processes involve the input of data (particles or fields) through a plane, as well as output through a series of planes, all of which are perpendicular to the direction of the relative velocity between the frame of calculation and the other frame of choice.
 
 Input in a boosted frame simulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -29,17 +30,17 @@ Input in a boosted frame simulation
 Particles -
 ^^^^^^^^^^^^
 
-Particles are launched through a plane using a technique that is generic and applies to Lorentz boosted frame simulations in general, including plasma acceleration, and is illustrated using the case of a positively charged particle beam propagating through a background of cold electrons in an assumed continuous transverse focusing system, leading to a well-known growing transverse “electron cloud” instability (Vay 2007). In the laboratory frame, the electron background is initially at rest and a moving window is used to follow the beam progression. Traditionally, the beam macroparticles are initialized all at once in the window, while background electron macroparticles are created continuously in front of the beam on a plane that is perpendicular to the beam velocity. In a frame moving at some fraction of the beam velocity in the laboratory frame, the beam initial conditions at a given time in the calculation frame are generally unknown and one must initialize the beam differently. However, it can be taken advantage of the fact that the beam initial conditions are often known for a given plane in the laboratory, either directly, or via simple calculation or projection from the conditions at a given time in the labortory frame. Given the position and velocity :math:`\{x,y,z,v_x,v_y,v_z\}` for each beam macroparticle at time :math:`t=0` for a beam moving at the average velocity :math:`v_b=\beta_b c` (where :math:`c` is the speed of light) in the laboratory, and using the standard synchronization (:math:`z=z'=0` at :math:`t=t'=0`) between the laboratory and the calculation frames, the procedure for transforming the beam quantities for injection in a boosted frame moving at velocity :math:`\beta c` in the laboratory is as follows (the superscript :math:`'` relates to quantities known in the boosted frame while the superscript :math:`^*` relates to quantities that are know at a given longitudinal position :math:`z^*` but different times of arrival):
+Particles are launched through a plane using a technique that is generic and applies to Lorentz boosted frame simulations in general, including plasma acceleration, and is illustrated using the case of a positively charged particle beam propagating through a background of cold electrons in an assumed continuous transverse focusing system, leading to a well-known growing transverse “electron cloud” instability :cite:p:`io-Vayprl07`. In the laboratory frame, the electron background is initially at rest and a moving window is used to follow the beam progression. Traditionally, the beam macroparticles are initialized all at once in the window, while background electron macroparticles are created continuously in front of the beam on a plane that is perpendicular to the beam velocity. In a frame moving at some fraction of the beam velocity in the laboratory frame, the beam initial conditions at a given time in the calculation frame are generally unknown and one must initialize the beam differently. However, it can be taken advantage of the fact that the beam initial conditions are often known for a given plane in the laboratory, either directly, or via simple calculation or projection from the conditions at a given time in the labortory frame. Given the position and velocity :math:`\{x,y,z,v_x,v_y,v_z\}` for each beam macroparticle at time :math:`t=0` for a beam moving at the average velocity :math:`v_b=\beta_b c` (where :math:`c` is the speed of light) in the laboratory, and using the standard synchronization (:math:`z=z'=0` at :math:`t=t'=0`) between the laboratory and the calculation frames, the procedure for transforming the beam quantities for injection in a boosted frame moving at velocity :math:`\beta c` in the laboratory is as follows (the superscript :math:`'` relates to quantities known in the boosted frame while the superscript :math:`^*` relates to quantities that are know at a given longitudinal position :math:`z^*` but different times of arrival):
 
 #. project positions at :math:`z^*=0` assuming ballistic propagation
 
    .. math::
 
       \begin{aligned}
-          t^* &=& \left(z-\bar{z}\right)/v_z \label{Eq:t*}\\
-          x^* &=& x-v_x t^* \label{Eq:x*}\\
-          y^* &=& y-v_y t^* \label{Eq:y*}\\
-          z^* &=& 0 \label{Eq:z*}\end{aligned}
+          t^* &= \left(z-\bar{z}\right)/v_z \label{Eq:t*}\\
+          x^* &= x-v_x t^* \label{Eq:x*}\\
+          y^* &= y-v_y t^* \label{Eq:y*}\\
+          z^* &= 0 \label{Eq:z*}\end{aligned}
 
    the velocity components being left unchanged,
 
@@ -48,13 +49,13 @@ Particles are launched through a plane using a technique that is generic and app
    .. math::
 
       \begin{aligned}
-          t'^* &=& -\gamma t^* \label{Eq:tp*}\\
-          x'^* &=& x^* \label{Eq:xp*}\\
-          y'^* &=& y^* \label{Eq:yp*}\\
-          z'^* &=& \gamma\beta c t^* \label{Eq:zp*}\\
-          v'^*_x&=&\frac{v_x^*}{\gamma\left(1-\beta \beta_b\right)} \label{Eq:vxp*}\\
-          v'^*_y&=&\frac{v_y^*}{\gamma\left(1-\beta \beta_b\right)} \label{Eq:vyp*}\\
-          v'^*_z&=&\frac{v_z^*-\beta c}{1-\beta \beta_b} \label{Eq:vzp*}\end{aligned}
+          t'^* &= -\gamma t^* \label{Eq:tp*}\\
+          x'^* &= x^* \label{Eq:xp*}\\
+          y'^* &= y^* \label{Eq:yp*}\\
+          z'^* &= \gamma\beta c t^* \label{Eq:zp*}\\
+          v'^*_x&=\frac{v_x^*}{\gamma\left(1-\beta \beta_b\right)} \label{Eq:vxp*}\\
+          v'^*_y&=\frac{v_y^*}{\gamma\left(1-\beta \beta_b\right)} \label{Eq:vyp*}\\
+          v'^*_z&=\frac{v_z^*-\beta c}{1-\beta \beta_b} \label{Eq:vzp*}\end{aligned}
 
    where :math:`\gamma=1/\sqrt{1-\beta^2}`. With the knowledge of the time at which each beam macroparticle crosses the plane into consideration, one can inject each beam macroparticle in the simulation at the appropriate location and time.
 
@@ -63,11 +64,11 @@ Particles are launched through a plane using a technique that is generic and app
    .. math::
 
       \begin{aligned}
-          z' &=& z'^*-\bar{v}'^*_z t'^* \label{Eq:zp}\end{aligned}
+          z' &= z'^*-\bar{v}'^*_z t'^* \label{Eq:zp}\end{aligned}
 
-   This additional step is needed for setting the electrostatic or electromagnetic fields at the plane of injection. In a Particle-In-Cell code, the three-dimensional fields are calculated by solving the Maxwell equations (or static approximation like Poisson, Darwin or other (Vay 2008)) on a grid on which the source term is obtained from the macroparticles distribution. This requires generation of a three-dimensional representation of the beam distribution of macroparticles at a given time before they cross the injection plane at :math:`z'^*`. This is accomplished by expanding the beam distribution longitudinally such that all macroparticles (so far known at different times of arrival at the injection plane) are synchronized to the same time in the boosted frame. To keep the beam shape constant, the particles are “frozen” until they cross that plane: the three velocity components and the two position components perpendicular to the boosted frame velocity are kept constant, while the remaining position component is advanced at the average beam velocity. As particles cross the plane of injection, they become regular “active” particles with full 6-D dynamics.
+   This additional step is needed for setting the electrostatic or electromagnetic fields at the plane of injection. In a Particle-In-Cell code, the three-dimensional fields are calculated by solving the Maxwell equations (or static approximation like Poisson, Darwin or other :cite:p:`io-Vaypop2008`) on a grid on which the source term is obtained from the macroparticles distribution. This requires generation of a three-dimensional representation of the beam distribution of macroparticles at a given time before they cross the injection plane at :math:`z'^*`. This is accomplished by expanding the beam distribution longitudinally such that all macroparticles (so far known at different times of arrival at the injection plane) are synchronized to the same time in the boosted frame. To keep the beam shape constant, the particles are “frozen” until they cross that plane: the three velocity components and the two position components perpendicular to the boosted frame velocity are kept constant, while the remaining position component is advanced at the average beam velocity. As particles cross the plane of injection, they become regular “active” particles with full 6-D dynamics.
 
-Figure :numref:`Fig_inputoutput` (top) shows a snapshot of a beam that has passed partly through the injection plane. As the frozen beam macroparticles pass through the injection plane (which moves opposite to the beam in the boosted frame), they are converted to “active" macroparticles. The charge or current density is accumulated from the active and the frozen particles, thus ensuring that the fields at the plane of injection are consistent.
+A snapshot of a beam that has passed partly through the injection plane in shown in :numref:`fig_inputoutput` (top). As the frozen beam macroparticles pass through the injection plane (which moves opposite to the beam in the boosted frame), they are converted to “active" macroparticles. The charge or current density is accumulated from the active and the frozen particles, thus ensuring that the fields at the plane of injection are consistent.
 
 Laser -
 ^^^^^^^^
@@ -89,93 +90,28 @@ If, for convenience, the injection plane is moving at constant velocity :math:`\
 .. math::
 
    \begin{aligned}
-   E_\perp\left(x,y,t\right)&=&\left(1-\beta_s\right)E_0 f\left(x,y,t\right)\nonumber \\
-   &\times& \sin\left[\left(1-\beta_s\right)\omega t+\phi\left(x,y,\omega\right)\right].\end{aligned}
+   E_\perp\left(x,y,t\right)&=\left(1-\beta_s\right)E_0 f\left(x,y,t\right)\sin\left[\left(1-\beta_s\right)\omega t+\phi\left(x,y,\omega\right)\right]
+   \end{aligned}
 
 The injection of a laser of frequency :math:`\omega` is considered for a simulation using a boosted frame moving at :math:`\beta c` with respect to the laboratory. Assuming that the laser is injected at a plane that is fixed in the laboratory, and thus moving at :math:`\beta_s=-\beta` in the boosted frame, the injection in the boosted frame is given by
 
 .. math::
 
    \begin{aligned}
-   E_\perp\left(x',y',t'\right)&=&\left(1-\beta_s\right)E'_0 f\left(x',y',t'\right)\nonumber \\
-   &\times&\sin\left[\left(1-\beta_s\right)\omega' t'+\phi\left(x',y',\omega'\right)\right]\\
-   &=&\left(E_0/\gamma\right) f\left(x',y',t'\right) \nonumber\\
-   &\times&\sin\left[\omega t'/\gamma+\phi\left(x',y',\omega'\right)\right]\end{aligned}
+   E_\perp\left(x',y',t'\right)&=\left(1-\beta_s\right)E'_0 f\left(x',y',t'\right)\sin\left[\left(1-\beta_s\right)\omega' t'+\phi\left(x',y',\omega'\right)\right]
+   \\
+   &=\left(E_0/\gamma\right) f\left(x',y',t'\right)\sin\left[\omega t'/\gamma+\phi\left(x',y',\omega'\right)\right]
+   \end{aligned}
 
 since :math:`E'_0/E_0=\omega'/\omega=1/\left(1+\beta\right)\gamma`.
 
-The electric field is then converted into currents that get injected via a 2D array of macro-particles, with one positive and one dual negative macro-particle for each array cell in the plane of injection, whose weights and motion are governed by :math:`E_\perp\left(x',y',t'\right)`. Injecting using this dual array of macroparticles offers the advantage of automatically including the longitudinal component that arises from emitting into a boosted frame, and to automatically verify the discrete Gauss’ law thanks to using charge conserving (e.g. Esirkepov) current deposition scheme (Esirkepov 2001).
+The electric field is then converted into currents that get injected via a 2D array of macro-particles, with one positive and one dual negative macro-particle for each array cell in the plane of injection, whose weights and motion are governed by :math:`E_\perp\left(x',y',t'\right)`. Injecting using this dual array of macroparticles offers the advantage of automatically including the longitudinal component that arises from emitting into a boosted frame, and to automatically verify the discrete Gauss’ law thanks to using charge conserving (e.g. Esirkepov) current deposition scheme :cite:p:`io-Esirkepovcpc01`.
 
 Output in a boosted frame simulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Some quantities, e.g. charge or dimensions perpendicular to the boost velocity, are Lorentz invariant.
-Those quantities are thus readily available from standard diagnostics in the boosted frame calculations. Quantities that do not fall in this category are recorded at a number of regularly spaced “stations", immobile in the laboratory frame, at a succession of time intervals to record data history, or averaged over time. A visual example is given on Fig. :numref:`Fig_inputoutput` (bottom). Since the space-time locations of the diagnostic grids in the laboratory frame generally do not coincide with the space-time positions of the macroparticles and grid nodes used for the calculation in a boosted frame, some interpolation is performed at runtime during the data collection process. As a complement or an alternative, selected particle or field quantities can be dumped at regular intervals and quantities are reconstructed in the laboratory frame during a post-processing phase. The choice of the methods depends on the requirements of the diagnostics and particular implementations.
-
-.. raw:: html
-
-   <div id="refs" class="references">
-
-.. raw:: html
-
-   <div id="ref-CowanPRSTAB13">
-
-Cowan, Benjamin M, David L Bruhwiler, John R Cary, Estelle Cormier-Michel, and Cameron G R Geddes. 2013. “Generalized algorithm for control of numerical dispersion in explicit time-domain electromagnetic simulations.” *Physical Review Special Topics-Accelerators and Beams* 16 (4). https://doi.org/10.1103/PhysRevSTAB.16.041303.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Esirkepovcpc01">
-
-Esirkepov, Tz. 2001. “Exact Charge Conservation Scheme for Particle-in-Cell Simulation with an Arbitrary Form-Factor.” *Computer Physics Communications* 135 (2): 144–53.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Warp">
-
-Grote, D P, A Friedman, J.-L. Vay, and I Haber. 2005. “The Warp Code: Modeling High Intensity Ion Beams.” In *Aip Conference Proceedings*, 55–58. 749.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Vaypop2008">
-
-Vay, J L. 2008. “Simulation of Beams or Plasmas Crossing at Relativistic Velocity.” *Physics of Plasmas* 15 (5): 56701. https://doi.org/10.1063/1.2837054.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Vayprl07">
-
-Vay, J.-L. 2007. “Noninvariance of Space- and Time-Scale Ranges Under A Lorentz Transformation and the Implications for the Study of Relativistic Interactions.” *Physical Review Letters* 98 (13): 130405/1–4.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
-
-   <div id="ref-Vaypop2011">
-
-Vay, J -L., C G R Geddes, E Esarey, C B Schroeder, W P Leemans, E Cormier-Michel, and D P Grote. 2011. “Modeling of 10 Gev-1 Tev Laser-Plasma Accelerators Using Lorentz Boosted Simulations.” *Physics of Plasmas* 18 (12). https://doi.org/10.1063/1.3663841.
-
-.. raw:: html
-
-   </div>
-
-.. raw:: html
+Those quantities are thus readily available from standard diagnostics in the boosted frame calculations. Quantities that do not fall in this category are recorded at a number of regularly spaced “stations", immobile in the laboratory frame, at a succession of time intervals to record data history, or averaged over time. A visual example is given on :numref:`fig_inputoutput` (bottom). Since the space-time locations of the diagnostic grids in the laboratory frame generally do not coincide with the space-time positions of the macroparticles and grid nodes used for the calculation in a boosted frame, some interpolation is performed at runtime during the data collection process. As a complement or an alternative, selected particle or field quantities can be dumped at regular intervals and quantities are reconstructed in the laboratory frame during a post-processing phase. The choice of the methods depends on the requirements of the diagnostics and particular implementations.
 
-   </div>
+.. bibliography::
+   :keyprefix: io-
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 61f1d8fc155..ec360135f6f 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -1389,12 +1389,12 @@ Laser initialization
 
         E(\boldsymbol{x},t) \propto Re\left[ \exp\left(  -\frac{(t-t_{peak})^2}{\tau^2 + 2i\phi^{(2)}} + i\omega_0 (t-t_{peak}) + i\phi_0 \right) \right]
 
-    where :math:`\tau` is given by ``<laser_name>.laser_duration`` and represents the
+    where :math:`\tau` is given by ``<laser_name>.profile_duration`` and represents the
     Fourier-limited duration of the laser pulse. Thus, the actual duration of the chirped laser pulse is:
 
     .. math::
 
-        \tau' = \sqrt{ \tau^2 + 4 \phi^{(2)}/\tau^2 }
+        \tau' = \sqrt{ \tau^2 + 4 (\phi^{(2)})^2/\tau^2 }
 
 * ``<laser_name>.do_continuous_injection`` (`0` or `1`) optional (default `0`).
     Whether or not to use continuous injection.
diff --git a/Docs/source/usage/python.rst b/Docs/source/usage/python.rst
index c9d846b41e7..2cf7982c69c 100644
--- a/Docs/source/usage/python.rst
+++ b/Docs/source/usage/python.rst
@@ -13,7 +13,7 @@ In the input file, instances of classes are created defining the various aspects
 A variable of type :py:class:`pywarpx.picmi.Simulation` is the central object to which all other options are passed, defining the simulation time, field solver, registered species, etc.
 
 Once the simulation is fully configured, it can be used in one of two modes.
-**Interactive** use is the most common and can be :ref:`extended with custom runtime functionality <usage-picmi-extend>`:
+**Interactive** use is the most common and can be :ref:`extended with custom runtime functionality <usage-python-extend>`:
 
 .. tab-set::
 
@@ -25,6 +25,9 @@ Once the simulation is fully configured, it can be used in one of two modes.
 
       :py:meth:`~pywarpx.picmi.Simulation.write_input_file`: create an :ref:`inputs file for a WarpX executable <running-cpp-parameters>`
 
+When run directly from Python, one can also extend WarpX with further custom user logic.
+See the :ref:`detailed workflow page <usage-python-extend>` on how to extend WarpX from Python.
+
 
 .. _usage-picmi-parameters:
 
@@ -137,118 +140,3 @@ Laser profiles can be used to initialize laser pulses in the simulation.
 Laser injectors control where to initialize laser pulses on the simulation grid.
 
 .. autoclass:: pywarpx.picmi.LaserAntenna
-
-
-.. _usage-picmi-extend:
-
-Extending a Simulation from Python
-----------------------------------
-
-When running WarpX directly from Python it is possible to interact with the simulation.
-
-For instance, with the :py:meth:`~pywarpx.picmi.Simulation.step` method of the simulation class, one could run ``sim.step(nsteps=1)`` in a loop:
-
-.. code-block:: python3
-
-   # Preparation: set up the simulation
-   #   sim = picmi.Simulation(...)
-   #   ...
-
-   steps = 1000
-   for _ in range(steps):
-       sim.step(nsteps=1)
-
-       # do something custom with the sim object
-
-As a more flexible alternative, one can install `callback functions <https://en.wikipedia.org/wiki/Callback_(computer_programming)>`__, which will execute a given Python function at a
-specific location in the WarpX simulation loop.
-
-.. automodule:: pywarpx.callbacks
-   :members: installcallback, uninstallcallback, isinstalled
-
-Data Access
-^^^^^^^^^^^
-
-While the simulation is running, callbacks can access the WarpX simulation data *in situ*.
-
-An important object for data access is ``Simulation.extension.warpx``, which is available only during the simulation run.
-This object is the Python equivalent to the C++ ``WarpX`` simulation class and provides access to field ``MultiFab`` and ``ParticleContainer`` data.
-
-.. py:function:: pywarpx.picmi.Simulation.extension.warpx.getistep()
-
-.. py:function:: pywarpx.picmi.Simulation.extension.warpx.gett_new()
-
-.. py:function:: pywarpx.picmi.Simulation.extension.warpx.evolve()
-
-.. autofunction:: pywarpx.picmi.Simulation.extension.finalize()
-
-These and other classes are provided through `pyAMReX <https://github.com/AMReX-Codes/pyamrex>`__.
-After the simulation is initialized, pyAMReX can be accessed via
-
-.. code-block:: python
-
-   from pywarpx import picmi, libwarpx
-
-   # ... simulation definition ...
-
-   # equivalent to
-   #   import amrex.space3d as amr
-   # for a 3D simulation
-   amr = libwarpx.amr  # picks the right 1d, 2d or 3d variant
-
-.. py:function:: amr.ParallelDescriptor.NProcs()
-
-.. py:function:: amr.ParallelDescriptor.MyProc()
-
-.. py:function:: amr.ParallelDescriptor.IOProcessor()
-
-.. py:function:: amr.ParallelDescriptor.IOProcessorNumber()
-
-Particles can be added to the simulation at specific positions and with specific attribute values:
-
-.. code-block:: python
-
-   from pywarpx import particle_containers, picmi
-
-   # ...
-
-   electron_wrapper = particle_containers.ParticleContainerWrapper("electrons")
-
-.. autofunction:: pywarpx.particle_containers.ParticleContainerWrapper.add_particles
-
-Properties of the particles already in the simulation can be obtained with various functions.
-
-.. autofunction:: pywarpx.particle_containers.ParticleContainerWrapper.get_particle_count
-
-.. autofunction:: pywarpx.particle_containers.ParticleContainerWrapper.get_particle_structs
-
-.. autofunction:: pywarpx.particle_containers.ParticleContainerWrapper.get_particle_arrays
-
-The ``get_particle_structs()`` and ``get_particle_arrays()`` functions are called
-by several utility functions of the form ``get_particle_{comp_name}`` where
-``comp_name`` is one of ``x``, ``y``, ``z``, ``r``, ``theta``, ``id``, ``cpu``,
-``weight``, ``ux``, ``uy`` or ``uz``.
-
-New components can be added via Python.
-
-.. autofunction:: pywarpx.particle_containers.ParticleContainerWrapper.add_real_comp
-
-Various diagnostics are also accessible from Python.
-This includes getting the deposited or total charge density from a given species as well as accessing the scraped particle buffer.
-See the example in ``Examples/Tests/ParticleBoundaryScrape`` for a reference on how to interact with scraped particle data.
-
-.. autofunction:: pywarpx.particle_containers.ParticleContainerWrapper.get_species_charge_sum
-
-.. autofunction:: pywarpx.particle_containers.ParticleContainerWrapper.deposit_charge_density
-
-.. autofunction:: pywarpx.particle_containers.ParticleBoundaryBufferWrapper.get_particle_boundary_buffer_size
-
-.. autofunction:: pywarpx.particle_containers.ParticleBoundaryBufferWrapper.get_particle_boundary_buffer_structs
-
-.. autofunction:: pywarpx.particle_containers.ParticleBoundaryBufferWrapper.get_particle_boundary_buffer
-
-.. autofunction:: pywarpx.particle_containers.ParticleBoundaryBufferWrapper.clear_buffer
-
-The embedded boundary conditions can be modified when using the electrostatic solver.
-
-.. py:function:: pywarpx.picmi.Simulation.extension.warpx.set_potential_on_eb()
diff --git a/Docs/source/usage/workflows.rst b/Docs/source/usage/workflows.rst
index 7a93d081c6b..e68ec9391be 100644
--- a/Docs/source/usage/workflows.rst
+++ b/Docs/source/usage/workflows.rst
@@ -8,11 +8,12 @@ This section collects typical user workflows and best practices for WarpX.
 .. toctree::
    :maxdepth: 2
 
+   workflows/python_extend
    workflows/domain_decomposition
+   workflows/plot_distribution_mapping
    workflows/debugging
    workflows/libensemble
    workflows/plot_timestep_duration
-   workflows/plot_distribution_mapping
    workflows/psatd_stencil
    workflows/archiving
    workflows/ml_dataset_training
diff --git a/Docs/source/usage/workflows/python_extend.rst b/Docs/source/usage/workflows/python_extend.rst
new file mode 100644
index 00000000000..6c9286c02ce
--- /dev/null
+++ b/Docs/source/usage/workflows/python_extend.rst
@@ -0,0 +1,279 @@
+.. _usage-python-extend:
+
+Extend a Simulation with Python
+===============================
+
+When running WarpX directly :ref:`from Python <usage-picmi>` it is possible to interact with the simulation.
+
+For instance, with the :py:meth:`~pywarpx.picmi.Simulation.step` method of the simulation class, one could run ``sim.step(nsteps=1)`` in a loop:
+
+.. code-block:: python3
+
+   # Preparation: set up the simulation
+   #   sim = picmi.Simulation(...)
+   #   ...
+
+   steps = 1000
+   for _ in range(steps):
+       sim.step(nsteps=1)
+
+       # do something custom with the sim object
+
+As a more flexible alternative, one can install `callback functions <https://en.wikipedia.org/wiki/Callback_(computer_programming)>`__, which will execute a given Python function at a
+specific location in the WarpX simulation loop.
+
+.. automodule:: pywarpx.callbacks
+   :members: installcallback, uninstallcallback, isinstalled
+
+
+pyAMReX
+-------
+
+Many of the following classes are provided through `pyAMReX <https://github.com/AMReX-Codes/pyamrex>`__.
+After the simulation is initialized, the pyAMReX module can be accessed via
+
+.. code-block:: python
+
+   from pywarpx import picmi, libwarpx
+
+   # ... simulation definition ...
+
+   # equivalent to
+   #   import amrex.space3d as amr
+   # for a 3D simulation
+   amr = libwarpx.amr  # picks the right 1d, 2d or 3d variant
+
+
+Full details for pyAMReX APIs are `documented here <https://pyamrex.readthedocs.io/en/latest/usage/api.html>`__.
+Important APIs include:
+
+* `amr.ParallelDescriptor <https://pyamrex.readthedocs.io/en/latest/usage/api.html#amrex.space3d.ParallelDescriptor.IOProcessor>`__: MPI-parallel rank information
+* `amr.MultiFab <https://pyamrex.readthedocs.io/en/latest/usage/api.html#amrex.space3d.MultiFab>`__: MPI-parallel field data
+* `amr.ParticleContainer_* <https://pyamrex.readthedocs.io/en/latest/usage/api.html#amrex.space3d.ParticleContainer_1_1_2_1_default>`__: MPI-parallel particle data for a particle species
+
+
+Data Access
+-----------
+
+While the simulation is running, callbacks can have read and write access the WarpX simulation data *in situ*.
+
+An important object in the ``pywarpx.picmi`` module for data access is ``Simulation.extension.warpx``, which is available only during the simulation run.
+This object is the Python equivalent to the C++ ``WarpX`` simulation class.
+
+.. py:class:: WarpX
+
+   .. py:method:: getistep(lev: int)
+
+      Get the current step on mesh-refinement level ``lev``.
+
+   .. py:method:: gett_new(lev: int)
+
+      Get the current physical time on mesh-refinement level ``lev``.
+
+   .. py:method:: getdt(lev: int)
+
+      Get the current physical time step size on mesh-refinement level ``lev``.
+
+   .. py:method:: multifab(multifab_name: str)
+
+      Return MultiFabs by name, e.g., ``"Efield_aux[x][level=0]"``, ``"Efield_cp[x][level=0]"``, ...
+
+      The physical fields in WarpX have the following naming:
+
+      - ``_fp`` are the "fine" patches, the regular resolution of a current mesh-refinement level
+      - ``_aux`` are temporary (auxiliar) patches at the same resolution as ``_fp``.
+        They usually include contributions from other levels and can be interpolated for gather routines of particles.
+      - ``_cp`` are "coarse" patches, at the same resolution (but not necessary values) as the ``_fp`` of ``level - 1``
+        (only for level 1 and higher).
+
+   .. py:method:: multi_particle_container
+
+   .. py:method:: get_particle_boundary_buffer
+
+   .. py:method:: set_potential_on_eb(potential: str)
+
+      The embedded boundary (EB) conditions can be modified when using the electrostatic solver.
+      This set the EB potential string and updates the function parser.
+
+   .. py:method:: evolve(numsteps=-1)
+
+      Evolve the simulation the specified number of steps.
+
+   .. autofunction:: pywarpx.picmi.Simulation.extension.finalize
+
+.. py::def:: pywarpx.picmi.Simulation.extension.get_instance
+
+   Return a reference to the WarpX object.
+
+
+The :py:class:`WarpX` also provides read and write access to field ``MultiFab`` and ``ParticleContainer`` data, shown in the following examples.
+
+Fields
+^^^^^^
+
+This example accesses the :math:`E_x(x,y,z)` field at level 0 after every time step and calculate the largest value in it.
+
+.. code-block:: python3
+
+   from pywarpx import picmi
+   from pywarpx.callbacks import callfromafterstep
+
+   # Preparation: set up the simulation
+   #   sim = picmi.Simulation(...)
+   #   ...
+
+
+   @callfromafterstep
+   def set_E_x():
+       warpx = sim.extension.warpx
+
+       # data access
+       E_x_mf = warpx.multifab(f"Efield_fp[x][level=0]")
+
+       # compute
+       # iterate over mesh-refinement levels
+       for lev in range(warpx.finest_level + 1):
+           # grow (aka guard/ghost/halo) regions
+           ngv = E_x_mf.n_grow_vect
+
+           # get every local block of the field
+           for mfi in E_x_mf:
+               # global index space box, including guards
+               bx = mfi.tilebox().grow(ngv)
+               print(bx)  # note: global index space of this block
+
+               # numpy representation: non-copying view, including the
+               # guard/ghost region;     .to_cupy() for GPU!
+               E_x_np = E_x_mf.array(mfi).to_numpy()
+
+               # notes on indexing in E_x_np:
+               # - numpy uses locally zero-based indexing
+               # - layout is F_CONTIGUOUS by default, just like AMReX
+
+               # notes:
+               # Only the next lines are the "HOT LOOP" of the computation.
+               # For efficiency, use numpy array operation for speed on CPUs.
+               # For GPUs use .to_cupy() above and compute with cupy or numba.
+               E_x_np[()] = 42.0
+
+
+   sim.step(nsteps=100)
+
+For further details on how to `access GPU data <https://pyamrex.readthedocs.io/en/latest/usage/zerocopy.html>`__ or compute on ``E_x``, please see the `pyAMReX documentation <https://pyamrex.readthedocs.io/en/latest/usage/compute.html#fields>`__.
+
+
+High-Level Field Wrapper
+""""""""""""""""""""""""
+
+.. note::
+
+   TODO
+
+.. note::
+
+   TODO: What are the benefits of using the high-level wrapper?
+   TODO: What are the limitations (e.g., in memory usage or compute scalability) of using the high-level wrapper?
+
+
+Particles
+^^^^^^^^^
+
+.. code-block:: python3
+
+   from pywarpx import picmi
+   from pywarpx.callbacks import callfromafterstep
+
+   # Preparation: set up the simulation
+   #   sim = picmi.Simulation(...)
+   #   ...
+
+   @callfromafterstep
+   def my_after_step_callback():
+       warpx = sim.extension.warpx
+
+       # data access
+       multi_pc = warpx.multi_particle_container()
+       pc = multi_pc.get_particle_container_from_name("electrons")
+
+       # compute
+       # iterate over mesh-refinement levels
+       for lvl in range(pc.finest_level + 1):
+           # get every local chunk of particles
+           for pti in pc.iterator(pc, level=lvl):
+               # default layout: AoS with positions and cpuid
+               aos = pti.aos().to_numpy()
+
+               # additional compile-time and runtime attributes in SoA format
+               soa = pti.soa().to_numpy()
+
+               # notes:
+               # Only the next lines are the "HOT LOOP" of the computation.
+               # For efficiency, use numpy array operation for speed on CPUs.
+               # For GPUs use .to_cupy() above and compute with cupy or numba.
+
+               # write to all particles in the chunk
+               # note: careful, if you change particle positions, you need to
+               #       redistribute particles before continuing the simulation step
+               # aos[()]["x"] = 0.30
+               # aos[()]["y"] = 0.35
+               # aos[()]["z"] = 0.40
+
+               for soa_real in soa.real:
+                   soa_real[()] = 42.0
+
+               for soa_int in soa.int:
+                   soa_int[()] = 12
+
+
+   sim.step(nsteps=100)
+
+For further details on how to `access GPU data <https://pyamrex.readthedocs.io/en/latest/usage/zerocopy.html>`__ or compute on ``electrons``, please see the `pyAMReX documentation <https://pyamrex.readthedocs.io/en/latest/usage/compute.html#particles>`__.
+
+
+High-Level Particle Wrapper
+"""""""""""""""""""""""""""
+
+.. note::
+
+   TODO: What are the benefits of using the high-level wrapper?
+   TODO: What are the limitations (e.g., in memory usage or compute scalability) of using the high-level wrapper?
+
+Particles can be added to the simulation at specific positions and with specific attribute values:
+
+.. code-block:: python
+
+   from pywarpx import particle_containers, picmi
+
+   # ...
+
+   electron_wrapper = particle_containers.ParticleContainerWrapper("electrons")
+
+
+.. autoclass:: pywarpx.particle_containers.ParticleContainerWrapper
+   :members:
+
+The ``get_particle_structs()`` and ``get_particle_arrays()`` functions are called
+by several utility functions of the form ``get_particle_{comp_name}`` where
+``comp_name`` is one of ``x``, ``y``, ``z``, ``r``, ``theta``, ``id``, ``cpu``,
+``weight``, ``ux``, ``uy`` or ``uz``.
+
+
+Diagnostics
+-----------
+
+Various diagnostics are also accessible from Python.
+This includes getting the deposited or total charge density from a given species as well as accessing the scraped particle buffer.
+See the example in ``Examples/Tests/ParticleBoundaryScrape`` for a reference on how to interact with scraped particle data.
+
+
+.. autoclass:: pywarpx.particle_containers.ParticleBoundaryBufferWrapper
+   :members:
+
+
+Modify Solvers
+--------------
+
+From Python, one can also replace numerical solvers in the PIC loop or add new physical processes into the time step loop.
+Examples:
+
+* :ref:`Capacitive Discharge <examples-capacitive-discharge>`: replaces the Poisson solver of an electrostatic simulation (default: MLMG) with a python function that uses `superLU <https://portal.nersc.gov/project/sparse/superlu/>`__ to directly solve the Poisson equation.
diff --git a/Python/pywarpx/callbacks.py b/Python/pywarpx/callbacks.py
index 1c191e607f4..8f3c1dc5e2c 100644
--- a/Python/pywarpx/callbacks.py
+++ b/Python/pywarpx/callbacks.py
@@ -7,7 +7,7 @@
 # License: BSD-3-Clause-LBNL
 
 """Callback Locations
-==================
+------------------
 
 These are the functions which allow installing user created functions so that
 they are called at various places along the time step.
@@ -23,9 +23,6 @@
 instance method as an argument. Note that if an instance method is used, an
 extra reference to the method's object is saved.
 
-The install can be done using a decorator, which has the prefix ``callfrom``. See
-example below.
-
 Functions can be called at the following times:
 
 * ``beforeInitEsolve``: before the initial solve for the E fields (i.e. before the PIC loop starts)
@@ -46,27 +43,37 @@
 * ``particleinjection``: called when particle injection happens, after the position
   advance and before deposition is called, allowing a user
   defined particle distribution to be injected each time step
-* ``appliedfields``: allows directly specifying any fields to be applied to the particles
-  during the advance
 
-To use a decorator, the syntax is as follows. This will install the function
-``myplots`` to be called after each step.
+Example that calls the Python function ``myplots`` after each step:
 
 .. code-block:: python3
 
-   @callfromafterstep
+   from pywarpx.callbacks import installcallback
+
    def myplots():
-       ppzx()
+       # do something here
+
+   installcallback('afterstep', myplots)
 
-This is equivalent to the following:
+   # run simulation
+   sim.step(nsteps=100)
+
+The install can also be done using a `Python decorator <https://docs.python.org/3/glossary.html#term-decorator>`__, which has the prefix ``callfrom``.
+To use a decorator, the syntax is as follows. This will install the function ``myplots`` to be called after each step.
+The above example is quivalent to the following:
 
 .. code-block:: python3
 
+   from pywarpx.callbacks import callfromafterstep
+
+   @callfromafterstep
    def myplots():
-       ppzx()
+       # do something here
 
-   installcallback('afterstep', myplots)
+   # run simulation
+   sim.step(nsteps=100)
 """
+
 import copy
 import sys
 import time
diff --git a/Source/FieldSolver/ElectrostaticSolver.H b/Source/FieldSolver/ElectrostaticSolver.H
index 7bd982bb4e0..7fc6c9bf6da 100644
--- a/Source/FieldSolver/ElectrostaticSolver.H
+++ b/Source/FieldSolver/ElectrostaticSolver.H
@@ -52,7 +52,7 @@ namespace ElectrostaticSolver {
         void buildParsers ();
         void buildParsersEB ();
 
-        /* \brief Sets the EB potential string and updates the parsers
+        /* \brief Sets the EB potential string and updates the function parser
          *
          * \param [in] potential The string value of the potential
          */
diff --git a/Source/Parallelization/CMakeLists.txt b/Source/Parallelization/CMakeLists.txt
index 48f0d3a49bd..d38f8d6bbf8 100644
--- a/Source/Parallelization/CMakeLists.txt
+++ b/Source/Parallelization/CMakeLists.txt
@@ -5,5 +5,6 @@ foreach(D IN LISTS WarpX_DIMS)
         GuardCellManager.cpp
         WarpXComm.cpp
         WarpXRegrid.cpp
+        WarpXSumGuardCells.cpp
     )
 endforeach()
diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package
index 629cfafea62..7930118a1d2 100644
--- a/Source/Parallelization/Make.package
+++ b/Source/Parallelization/Make.package
@@ -1,5 +1,6 @@
 CEXE_sources += WarpXComm.cpp
 CEXE_sources += WarpXRegrid.cpp
 CEXE_sources += GuardCellManager.cpp
+CEXE_sources += WarpXSumGuardCells.cpp
 
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Parallelization
diff --git a/Source/Parallelization/WarpXSumGuardCells.H b/Source/Parallelization/WarpXSumGuardCells.H
index 13d6dce4af2..260ebb3871a 100644
--- a/Source/Parallelization/WarpXSumGuardCells.H
+++ b/Source/Parallelization/WarpXSumGuardCells.H
@@ -8,10 +8,6 @@
 #ifndef WARPX_SUM_GUARD_CELLS_H_
 #define WARPX_SUM_GUARD_CELLS_H_
 
-#include "Utils/WarpXAlgorithmSelection.H"
-
-#include <ablastr/utils/Communication.H>
-
 #include <AMReX_MultiFab.H>
 
 /** \brief Sum the values of `mf`, where the different boxes overlap
@@ -26,21 +22,10 @@
  *    updates both the *valid* cells and *guard* cells. (This is because a
  *    spectral solver requires the value of the sources over a large stencil.)
  */
-inline void
+void
 WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
                    const amrex::IntVect& src_ngrow,
-                   const int icomp=0, const int ncomp=1)
-{
-    amrex::IntVect n_updated_guards;
-
-    // Update both valid cells and guard cells
-    if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
-        n_updated_guards = mf.nGrowVect();
-    } else {  // Update only the valid cells
-        n_updated_guards = amrex::IntVect::TheZeroVector();
-    }
-    ablastr::utils::communication::SumBoundary(mf, icomp, ncomp, src_ngrow, n_updated_guards, WarpX::do_single_precision_comms, period);
-}
+                   int icomp=0, int ncomp=1);
 
 /** \brief Sum the values of `src` where the different boxes overlap
  * (i.e. in the guard cells) and copy them into `dst`
@@ -57,24 +42,10 @@ WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
  * Note: `i_comp` is the component where the results will be stored in `dst`;
  *       The component from which we copy in `src` is always 0.
  */
-inline void
+void
 WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src,
                    const amrex::Periodicity& period,
                    const amrex::IntVect& src_ngrow,
-                   const int icomp=0, const int ncomp=1)
-{
-    amrex::IntVect n_updated_guards;
-
-    // Update both valid cells and guard cells
-    if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
-        n_updated_guards = dst.nGrowVect();
-    } else {  // Update only the valid cells
-        n_updated_guards = amrex::IntVect::TheZeroVector();
-    }
-
-    dst.setVal(0., icomp, ncomp, n_updated_guards);
-//    ablastr::utils::communication::ParallelAdd(dst, src, 0, icomp, ncomp, src_ngrow, n_updated_guards, period);
-    dst.ParallelAdd(src, 0, icomp, ncomp, src_ngrow, n_updated_guards, period);
-}
+                   int icomp=0, int ncomp=1);
 
 #endif // WARPX_SUM_GUARD_CELLS_H_
diff --git a/Source/Parallelization/WarpXSumGuardCells.cpp b/Source/Parallelization/WarpXSumGuardCells.cpp
new file mode 100644
index 00000000000..20ac73cab50
--- /dev/null
+++ b/Source/Parallelization/WarpXSumGuardCells.cpp
@@ -0,0 +1,51 @@
+/* Copyright 2019 Maxence Thevenet, Remi Lehe, Weiqun Zhang
+ *
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "WarpXSumGuardCells.H"
+
+#include "Utils/WarpXAlgorithmSelection.H"
+
+#include "WarpX.H"
+
+#include <ablastr/utils/Communication.H>
+
+void
+WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
+                   const amrex::IntVect& src_ngrow,
+                   const int icomp, const int ncomp)
+{
+    amrex::IntVect n_updated_guards;
+
+    // Update both valid cells and guard cells
+    if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
+        n_updated_guards = mf.nGrowVect();
+    } else {  // Update only the valid cells
+        n_updated_guards = amrex::IntVect::TheZeroVector();
+    }
+    ablastr::utils::communication::SumBoundary(mf, icomp, ncomp, src_ngrow, n_updated_guards, WarpX::do_single_precision_comms, period);
+}
+
+
+void
+WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src,
+                   const amrex::Periodicity& period,
+                   const amrex::IntVect& src_ngrow,
+                   const int icomp, const int ncomp)
+{
+    amrex::IntVect n_updated_guards;
+
+    // Update both valid cells and guard cells
+    if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
+        n_updated_guards = dst.nGrowVect();
+    } else {  // Update only the valid cells
+        n_updated_guards = amrex::IntVect::TheZeroVector();
+    }
+
+    dst.setVal(0., icomp, ncomp, n_updated_guards);
+    dst.ParallelAdd(src, 0, icomp, ncomp, src_ngrow, n_updated_guards, period);
+}
diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp
index 2ceb17be116..e42b8268fe1 100644
--- a/Source/Python/WarpX.cpp
+++ b/Source/Python/WarpX.cpp
@@ -86,7 +86,15 @@ void init_WarpX (py::module& m)
             "Evolve the simulation the specified number of steps"
         )
 
-        // from AmrCore->AmrMesh
+        // from amrex::AmrCore / amrex::AmrMesh
+        .def_property_readonly("max_level",
+            [](WarpX const & wx){ return wx.maxLevel(); },
+            "The maximum mesh-refinement level for the simulation."
+        )
+        .def_property_readonly("finest_level",
+            [](WarpX const & wx){ return wx.finestLevel(); },
+            "The currently finest level of mesh-refinement used. This is always less or equal to max_level."
+        )
         .def("Geom",
             //[](WarpX const & wx, int const lev) { return wx.Geom(lev); },
             py::overload_cast< int >(&WarpX::Geom, py::const_),
@@ -112,7 +120,15 @@ void init_WarpX (py::module& m)
             },
             py::arg("multifab_name"),
             py::return_value_policy::reference_internal,
-            "Return MultiFabs by name, e.g., 'Efield_aux[x][l=0]', 'Efield_cp[x][l=0]', ..."
+            R"doc(Return MultiFabs by name, e.g., ``\"Efield_aux[x][level=0]\"``, ``\"Efield_cp[x][level=0]\"``, ...
+
+The physical fields in WarpX have the following naming:
+
+- ``_fp`` are the "fine" patches, the regular resolution of a current mesh-refinement level
+- ``_aux`` are temporary (auxiliar) patches at the same resolution as ``_fp``.
+  They usually include contributions from other levels and can be interpolated for gather routines of particles.
+- ``_cp`` are "coarse" patches, at the same resolution (but not necessary values) as the ``_fp`` of ``level - 1``
+  (only for level 1 and higher).)doc"
         )
         .def("multi_particle_container",
             [](WarpX& wx){ return &wx.GetPartContainer(); },
@@ -140,22 +156,26 @@ void init_WarpX (py::module& m)
         // Expose functions to get the current simulation step and time
         .def("getistep",
             [](WarpX const & wx, int lev){ return wx.getistep(lev); },
-            py::arg("lev")
+            py::arg("lev"),
+            "Get the current step on mesh-refinement level ``lev``."
         )
         .def("gett_new",
             [](WarpX const & wx, int lev){ return wx.gett_new(lev); },
-            py::arg("lev")
+            py::arg("lev"),
+            "Get the current physical time on mesh-refinement level ``lev``."
         )
         .def("getdt",
             [](WarpX const & wx, int lev){ return wx.getdt(lev); },
-            py::arg("lev")
+            py::arg("lev"),
+            "Get the current physical time step size on mesh-refinement level ``lev``."
         )
 
         .def("set_potential_on_eb",
             [](WarpX& wx, std::string potential) {
                 wx.m_poisson_boundary_handler.setPotentialEB(potential);
             },
-            py::arg("potential")
+            py::arg("potential"),
+            "Sets the EB potential string and updates the function parser."
         )
     ;