Add background stopping (#2884)

* Added BackgroundStopping

* Added BackgroundStopping CI test

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

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

* For background stopping, allowed parsed background density and temperature

* Updated background stopping CI test

* For background stopping, temperature is specified in Kelvin

* Added documentation for background stopping

* Added ion stopping plus other fixes

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

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

* Made dt level dependent

* Revert "Made dt level dependent"

This reverts commit 1aed9654d24ff1411a7b5fd8558891e0688f0032.
The collisions should be done using the time step of the lowest level.

* Add a comment about dt

* Add const declaration

* Added comment regarding stopped particle

* Several fixes and clean up of the documentation

* Update CI benchmark

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

12 files changed, 771 insertions, 15 deletions

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 32eeec61a..1e70c27b0 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -1353,14 +1353,7 @@ External fields
 Collision initialization
 ------------------------
 
-WarpX provides a relativistic elastic Monte Carlo binary Coulomb collision model,
-following the algorithm given by `Perez et al. (Phys. Plasmas 19, 083104, 2012) <https://doi.org/10.1063/1.4742167>`_.
-When the RZ mode is used, `warpx.n_rz_azimuthal_modes` must be set to 1 at the moment,
-since the current implementation of the collision module assumes axisymmetry.
-A non-relativistic Monte Carlo treatment for particles colliding
-with a neutral, uniform background gas is also available. The implementation follows the so-called
-null collision strategy discussed for example in `Birdsall (IEEE Transactions on
-Plasma Science, vol. 19, no. 2, pp. 65-85, 1991) <https://ieeexplore.ieee.org/document/106800>`_.
+WarpX provides several particle collions models, using varying degrees of approximation.
 
 * ``collisions.collision_names`` (`strings`, separated by spaces)
     The name of each collision type.
@@ -1368,8 +1361,21 @@ Plasma Science, vol. 19, no. 2, pp. 65-85, 1991) <https://ieeexplore.ieee.org/do
     in this documentation we use ``<collision_name>`` as a placeholder.
 
 * ``<collision_name>.type`` (`string`) optional
-    The type of collsion. The types implemented are ``pairwisecoulomb`` for pairwise Coulomb collisions and
-    ``background_mcc`` for collisions between particles and a neutral background. If not specified, it defaults to ``pairwisecoulomb``.
+    The type of collsion. The types implemented are:
+
+    - ``pairwisecoulomb`` for pairwise Coulomb collisions, the default if unspecified.
+      This provides a pair-wise relativistic elastic Monte Carlo binary Coulomb collision model,
+      following the algorithm given by `Perez et al. (Phys. Plasmas 19, 083104, 2012) <https://doi.org/10.1063/1.4742167>`_.
+      When the RZ mode is used, `warpx.n_rz_azimuthal_modes` must be set to 1 at the moment,
+      since the current implementation of the collision module assumes axisymmetry.
+    - ``background_mcc`` for collisions between particles and a neutral background.
+      This is a non-relativistic Monte Carlo treatment for particles colliding
+      with a neutral, uniform background gas. The implementation follows the so-called
+      null collision strategy discussed for example in `Birdsall (IEEE Transactions on
+      Plasma Science, vol. 19, no. 2, pp. 65-85, 1991) <https://ieeexplore.ieee.org/document/106800>`_.
+    - ``background_stopping`` for slowing of ions due to collisions with electrons or ions.
+      This implements the approximate formulae as derived in Introduction to Plasma Physics,
+      from Goldston and Rutherford, section 14.2.
 
 * ``<collision_name>.species`` (`strings`)
     If using ``pairwisecoulomb`` type this should be the names of two species,
@@ -1395,21 +1401,57 @@ Plasma Science, vol. 19, no. 2, pp. 65-85, 1991) <https://ieeexplore.ieee.org/do
     Execute collision every # time steps. The default value is 1.
 
 * ``<collision_name>.background_density`` (`float`)
-    Only for ``background_mcc``. The density of the neutral background gas in :math:`m^{-3}`.
+    Only for ``background_mcc`` and ``background_stopping``. The density of the background in :math:`m^{-3}`.
     Can also provide ``<collision_name>.background_density(x,y,z,t)`` using the parser
-    initialization style for spatially and temporally varying density. If a function
+    initialization style for spatially and temporally varying density. With ``background_mcc``, if a function
     is used for the background density, the input parameter ``<collision_name>.max_background_density``
     must also be provided to calculate the maximum collision probability.
 
 * ``<collision_name>.background_temperature`` (`float`)
-    Only for ``background_mcc``. The temperature of the neutral background gas in Kelvin.
+    Only for ``background_mcc`` and ``background_stopping``. The temperature of the background in Kelvin.
     Can also provide ``<collision_name>.background_temperature(x,y,z,t)`` using the parser
     initialization style for spatially and temporally varying temperature.
 
 * ``<collision_name>.background_mass`` (`float`) optional
-    Only for ``background_mcc``. The mass of the background gas in kg. If not
-    given the mass of the colliding species will be used unless ionization is
+    Only for ``background_mcc`` and ``background_stopping``. The mass of the background gas in kg.
+    With ``background_mcc``, if not given the mass of the colliding species will be used unless ionization is
     included in which case the mass of the product species will be used.
+    With ``background_stopping``, and ``background_type`` set to ``electrons``, if not given defaults to the electron mass. With
+    ``background_type`` set to ``ions``, the mass must be given.
+
+* ``<collision_name>.background_charge_state`` (`float`)
+    Only for ``background_stopping``, where it is required when ``background_type`` is set to ``ions``.
+    This specifies the charge state of the background ions.
+
+* ``<collision_name>.background_type`` (`string`)
+    Only for ``background_stopping``, where it is required, the type of the background.
+    The possible values are ``electrons`` and ``ions``. When ``electrons``, equation 14.12 from Goldston and Rutherford is used.
+    This formula is based on Coulomb collisions with the approximations that :math:`M_b >> m_e` and :math:`V << v_{thermal\_e}`,
+    and the assumption that the electrons have a Maxwellian distribution with temperature :math:`T_e`.
+
+    .. math::
+        \frac{dV}{dt} = - \frac{2^{1/2}n_eZ_b^2e^4m_e^{1/2}\log\Lambda}{12\pi^{3/2}\epsilon_0M_bT_e^{3/2}}V
+
+    where :math:`V` is each velocity component, :math:`n_e` is the background density, :math:`Z_b` is the ion charge state,
+    :math:`e` is the electron charge, :math:`m_e` is the background mass, :math:`\log\Lambda=\log((12\pi/Z_b)(n_e\lambda_{de}^3))`,
+    :math:`\lambda_{de}` is the DeBye length, and :math:`M_b` is the ion mass.
+    The equation is integrated over a time step, giving :math:`V(t+dt) = V(t)*\exp(-\alpha*{dt})`
+    where :math:`\alpha` is the factor multiplying :math:`V`.
+
+    When ``ions``, equation 14.20 is used.
+    This formula is based on Coulomb collisions with the approximations that :math:`M_b >> M` and :math:`V >> v_{thermal\_i}`.
+    The background ion temperature only appears in the :math:`\log\Lambda` term.
+
+    .. math::
+        \frac{dW_b}{dt} = - \frac{2^{1/2}n_iZ^2Z_b^2e^4M_b^{1/2}\log\Lambda}{8\pi\epsilon_0MW_b^{1/2}}
+
+    where :math:`W_b` is the ion energy, :math:`n_i` is the background density,
+    :math:`Z` is the charge state of the background ions, :math:`Z_b` is the ion charge state,
+    :math:`e` is the electron charge, :math:`M_b` is the ion mass, :math:`\log\Lambda=\log((12\pi/Z_b)(n_i\lambda_{di}^3))`,
+    :math:`\lambda_{di}` is the DeBye length, and :math:`M` is the background ion mass.
+    The equation is integrated over a time step, giving :math:`W_b(t+dt) = ((W_b(t)^{3/2}) - 3/2\beta{dt})^{2/3}`
+    where :math:`\beta` is the term on the r.h.s except :math:`W_b`.
+
 
 * ``<collision_name>.scattering_processes`` (`strings` separated by spaces)
     Only for ``background_mcc``. The MCC scattering processes that should be
diff --git a/Examples/Tests/ion_stopping/analysis_ion_stopping.py b/Examples/Tests/ion_stopping/analysis_ion_stopping.py
new file mode 100755
index 000000000..1ddec4a62
--- /dev/null
+++ b/Examples/Tests/ion_stopping/analysis_ion_stopping.py
@@ -0,0 +1,181 @@
+#!/usr/bin/env python3
+
+# Copyright 2022 David Grote
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# This script tests the ion stopping model.
+# It uses the same stopping power formula that
+# is used in the C++ to check the resulting
+# particle energies.
+
+import os
+import re
+import sys
+
+import numpy as np
+from scipy.constants import c, e, epsilon_0, k, m_e, m_p
+import yt
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# Define constants using the WarpX names for the evals below
+q_e = e
+kb = k
+
+# Tolerance on the error in the final energy (in eV)
+tolerance = 1.e-7
+
+last_filename = sys.argv[1]
+
+# Remove trailing '/' from file name, if necessary
+last_filename.rstrip('/')
+# Find last iteration in file name, such as 'test_name_plt000001' (last_it = '000001')
+last_it = re.search('\d+$', last_filename).group()
+# Find output prefix in file name, such as 'test_name_plt000001' (prefix = 'test_name_plt')
+prefix = last_filename[:-len(last_it)]
+
+def stopping_from_electrons(ne, Te, Zb, ion_mass):
+    """Calculate the coefficient in equation 14.13 from
+    "Introduction to Plasma Physics", Goldston and Rutherford.
+    ne: electron density
+    Te: electron temperature (eV)
+    Zb: ion charge state
+    ion_mass: (kg)
+    """
+    vthe = np.sqrt(3.*Te*e/m_e)
+    wpe = np.sqrt(ne*e**2/(epsilon_0*m_e))
+    lambdadb = vthe/wpe
+    loglambda = np.log((12.*np.pi/Zb)*(ne*lambdadb**3))
+    # Goldston's equation 14.13
+    dEdt = - np.sqrt(2.)*ne*Zb**2*e**4*np.sqrt(m_e)*loglambda/(6.*np.pi**1.5*epsilon_0**2*ion_mass*(Te*e)**1.5)
+    return dEdt
+
+def stopping_from_ions(dt, ni, Ti, mi, Zi, Zb, ion_mass, ion_energy):
+    """
+    ni: background ion density
+    Ti: background ion temperature (eV)
+    mi: background ion mass
+    Zi: background ion charge state
+    Zb: ion charge state
+    ion_mass: (kg)
+    ion_energy: (eV)
+    """
+    vthi = np.sqrt(3.*Ti*e/mi)
+    wpi = np.sqrt(ni*e**2/(epsilon_0*mi))
+    lambdadb = vthi/wpi
+    loglambda = np.log((12.*np.pi/Zb)*(ni*lambdadb**3))
+    alpha = np.sqrt(2.)*ni*Zi**2*Zb**2*e**4*np.sqrt(ion_mass)*loglambda/(8.*np.pi*epsilon_0**2*mi)
+    f1 = np.clip((ion_energy*e)**(3./2.) - 3./2.*alpha*dt, 0., None)
+    ion_energy = (f1)**(2./3.)/e
+    return ion_energy
+
+# Fetch background parameters and inital particle data
+ds0 = yt.load(f'{prefix}{len(last_it)*"0"}')
+ad0 = ds0.all_data()
+
+Zb = 1.  # Ion charge state
+
+ne = float(ds0.parameters['stopping_on_electrons_constant.background_density'])
+Te = eval(ds0.parameters['stopping_on_electrons_constant.background_temperature'])*kb/e
+ion_mass12 = m_p
+
+mi = eval(ds0.parameters['stopping_on_ions_constant.background_mass'])
+Zi = float(ds0.parameters['stopping_on_ions_constant.background_charge_state'])
+ni = float(ds0.parameters['stopping_on_ions_constant.background_density'])
+Ti = eval(ds0.parameters['stopping_on_ions_constant.background_temperature'])*kb/e
+ion_mass34 = 4.*m_p
+
+# For ions1, the background parameters are constants
+vx = ad0[('ions1', 'particle_momentum_x')].to_ndarray()/ion_mass12
+vy = ad0[('ions1', 'particle_momentum_y')].to_ndarray()/ion_mass12
+vz = ad0[('ions1', 'particle_momentum_z')].to_ndarray()/ion_mass12
+EE1 = 0.5*ion_mass12*(vx**2 + vy**2 + vz**2)/e
+
+# For ions2, the background parameters are parsed
+xx = ad0[('ions2', 'particle_position_x')].to_ndarray()/ion_mass12
+yy = ad0[('ions2', 'particle_position_y')].to_ndarray()/ion_mass12
+ne2 = np.where(xx > 0., 1.e20, 1.e21)
+Te2 = np.where(yy > 0., 1., 2.)
+
+vx = ad0[('ions2', 'particle_momentum_x')].to_ndarray()/ion_mass12
+vy = ad0[('ions2', 'particle_momentum_y')].to_ndarray()/ion_mass12
+vz = ad0[('ions2', 'particle_momentum_z')].to_ndarray()/ion_mass12
+EE2 = 0.5*ion_mass12*(vx**2 + vy**2 + vz**2)/e
+
+# For ions3, the background parameters are constants
+vx = ad0[('ions3', 'particle_momentum_x')].to_ndarray()/ion_mass34
+vy = ad0[('ions3', 'particle_momentum_y')].to_ndarray()/ion_mass34
+vz = ad0[('ions3', 'particle_momentum_z')].to_ndarray()/ion_mass34
+EE3 = 0.5*ion_mass34*(vx**2 + vy**2 + vz**2)/e
+
+# For ions4, the background parameters are parsed
+xx = ad0[('ions4', 'particle_position_x')].to_ndarray()/ion_mass34
+yy = ad0[('ions4', 'particle_position_y')].to_ndarray()/ion_mass34
+ni4 = np.where(xx > 0., 1.e20, 1.e21)
+Ti4 = np.where(yy > 0., 0.05, 0.10)
+
+vx = ad0[('ions4', 'particle_momentum_x')].to_ndarray()/ion_mass34
+vy = ad0[('ions4', 'particle_momentum_y')].to_ndarray()/ion_mass34
+vz = ad0[('ions4', 'particle_momentum_z')].to_ndarray()/ion_mass34
+EE4 = 0.5*ion_mass34*(vx**2 + vy**2 + vz**2)/e
+
+
+ds = yt.load(last_filename)
+ad = ds.all_data()
+dt = ds.current_time.to_value()/int(last_it)
+
+# Step through the same number of time steps
+a_EE1 = EE1
+a_EE2 = EE2
+a_EE3 = EE3
+a_EE4 = EE4
+for it in range(int(last_it)):
+    dEdt1 = stopping_from_electrons(ne, Te, Zb, ion_mass12)
+    a_EE1 *= np.exp(dEdt1*dt)
+    dEdt2 = stopping_from_electrons(ne2, Te2, Zb, ion_mass12)
+    a_EE2 *= np.exp(dEdt2*dt)
+    a_EE3 = stopping_from_ions(dt, ni, Ti, mi, Zi, Zb, ion_mass34, a_EE3)
+    a_EE4 = stopping_from_ions(dt, ni4, Ti4, mi, Zi, Zb, ion_mass34, a_EE4)
+
+# Fetch the final particle data
+vx = ad[('ions1', 'particle_momentum_x')].to_ndarray()/ion_mass12
+vy = ad[('ions1', 'particle_momentum_y')].to_ndarray()/ion_mass12
+vz = ad[('ions1', 'particle_momentum_z')].to_ndarray()/ion_mass12
+EE1 = 0.5*ion_mass12*(vx**2 + vy**2 + vz**2)/e
+
+vx = ad[('ions2', 'particle_momentum_x')].to_ndarray()/ion_mass12
+vy = ad[('ions2', 'particle_momentum_y')].to_ndarray()/ion_mass12
+vz = ad[('ions2', 'particle_momentum_z')].to_ndarray()/ion_mass12
+EE2 = 0.5*ion_mass12*(vx**2 + vy**2 + vz**2)/e
+
+vx = ad[('ions3', 'particle_momentum_x')].to_ndarray()/ion_mass34
+vy = ad[('ions3', 'particle_momentum_y')].to_ndarray()/ion_mass34
+vz = ad[('ions3', 'particle_momentum_z')].to_ndarray()/ion_mass34
+EE3 = 0.5*ion_mass34*(vx**2 + vy**2 + vz**2)/e
+
+vx = ad[('ions4', 'particle_momentum_x')].to_ndarray()/ion_mass34
+vy = ad[('ions4', 'particle_momentum_y')].to_ndarray()/ion_mass34
+vz = ad[('ions4', 'particle_momentum_z')].to_ndarray()/ion_mass34
+EE4 = 0.5*ion_mass34*(vx**2 + vy**2 + vz**2)/e
+
+error1 = np.abs(EE1 - a_EE1)
+error2 = np.abs(EE2 - a_EE2)
+error3 = np.abs(EE3 - a_EE3)
+error4 = np.abs(EE4 - a_EE4)
+print('stopping on electrons error with constant = ', error1)
+print('stopping on electrons error with parsed = ', error2)
+print('stopping on ions error with constant = ', error3)
+print('stopping on ions error with parsed = ', error4)
+print('tolerance = ', tolerance)
+
+assert np.all(error1 < tolerance)
+assert np.all(error2 < tolerance)
+assert np.all(error3 < tolerance)
+assert np.all(error4 < tolerance)
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, last_filename)
diff --git a/Examples/Tests/ion_stopping/inputs_3d b/Examples/Tests/ion_stopping/inputs_3d
new file mode 100644
index 000000000..9352946b8
--- /dev/null
+++ b/Examples/Tests/ion_stopping/inputs_3d
@@ -0,0 +1,100 @@
+my_constants.max_ion_energy = 100.  # eV
+my_constants.max_ion_velocity12 = sqrt(2.*max_ion_energy*q_e/m_p)/clight
+my_constants.max_ion_velocity34 = sqrt(2.*max_ion_energy*q_e/(4.*m_p))/clight
+
+max_step = 10
+
+amr.n_cell = 8 8 8
+amr.max_level = 0
+geometry.dims = 3
+geometry.prob_lo = -1. -1. -1.
+geometry.prob_hi = +1. +1. +1.
+boundary.field_lo = periodic periodic periodic
+boundary.field_hi = periodic periodic periodic
+boundary.particle_lo = periodic periodic periodic
+boundary.particle_hi = periodic periodic periodic
+algo.particle_shape = 1
+
+particles.species_names = ions1 ions2 ions3 ions4
+
+ions1.species_type = proton
+ions1.injection_style = "MultipleParticles"
+ions1.multiple_particles_pos_x = 0. 0. 0. 0.
+ions1.multiple_particles_pos_y = 0. 0. 0. 0.
+ions1.multiple_particles_pos_z = 0. 0. 0. 0.
+ions1.multiple_particles_vel_x = max_ion_velocity12 0. 0. max_ion_velocity12/2.
+ions1.multiple_particles_vel_y = 0. max_ion_velocity12 0. max_ion_velocity12/2.
+ions1.multiple_particles_vel_z = 0. 0. -max_ion_velocity12 -max_ion_velocity12/2.
+ions1.multiple_particles_weight = 1. 1. 1. 1.
+ions1.do_not_deposit = 1
+
+ions2.species_type = proton
+ions2.injection_style = "MultipleParticles"
+ions2.multiple_particles_pos_x = -0.1 +0.1 -0.1 +0.1
+ions2.multiple_particles_pos_y = -0.1 -0.1 +0.1 +0.1
+ions2.multiple_particles_pos_z = 0. 0. 0. 0.
+ions2.multiple_particles_vel_x = max_ion_velocity12 max_ion_velocity12 max_ion_velocity12 max_ion_velocity12
+ions2.multiple_particles_vel_y = 0. 0. 0. 0.
+ions2.multiple_particles_vel_z = 0. 0. 0. 0.
+ions2.multiple_particles_weight = 1. 1. 1. 1.
+ions2.do_not_deposit = 1
+
+ions3.charge = q_e
+ions3.mass = 4*m_p
+ions3.injection_style = "MultipleParticles"
+ions3.multiple_particles_pos_x = 0. 0. 0. 0.
+ions3.multiple_particles_pos_y = 0. 0. 0. 0.
+ions3.multiple_particles_pos_z = 0. 0. 0. 0.
+ions3.multiple_particles_vel_x = max_ion_velocity34 0. 0. max_ion_velocity34/2.
+ions3.multiple_particles_vel_y = 0. max_ion_velocity34 0. max_ion_velocity34/2.
+ions3.multiple_particles_vel_z = 0. 0. -max_ion_velocity34 -max_ion_velocity34/2.
+ions3.multiple_particles_weight = 1. 1. 1. 1.
+ions3.do_not_deposit = 1
+
+ions4.charge = q_e
+ions4.mass = 4*m_p
+ions4.injection_style = "MultipleParticles"
+ions4.multiple_particles_pos_x = -0.1 +0.1 -0.1 +0.1
+ions4.multiple_particles_pos_y = -0.1 -0.1 +0.1 +0.1
+ions4.multiple_particles_pos_z = 0. 0. 0. 0.
+ions4.multiple_particles_vel_x = max_ion_velocity34 max_ion_velocity34 max_ion_velocity34 max_ion_velocity34
+ions4.multiple_particles_vel_y = 0. 0. 0. 0.
+ions4.multiple_particles_vel_z = 0. 0. 0. 0.
+ions4.multiple_particles_weight = 1. 1. 1. 1.
+ions4.do_not_deposit = 1
+
+collisions.collision_names = stopping_on_electrons_constant stopping_on_electrons_parsed stopping_on_ions_constant stopping_on_ions_parsed
+
+stopping_on_electrons_constant.type = background_stopping
+stopping_on_electrons_constant.species = ions1
+stopping_on_electrons_constant.background_type = electrons
+stopping_on_electrons_constant.background_mass = m_e
+stopping_on_electrons_constant.background_density = 1.e20
+stopping_on_electrons_constant.background_temperature = 1.*q_e/kb # Kelvin
+
+stopping_on_electrons_parsed.type = background_stopping
+stopping_on_electrons_parsed.species = ions2
+stopping_on_electrons_parsed.background_type = electrons
+stopping_on_electrons_parsed.background_density(x,y,z,t) = if(x>0,1.e20,1.e21)
+stopping_on_electrons_parsed.background_temperature(x,y,z,t) = if(y>0,1.,2.)*q_e/kb
+
+stopping_on_ions_constant.type = background_stopping
+stopping_on_ions_constant.species = ions3
+stopping_on_ions_constant.background_type = ions
+stopping_on_ions_constant.background_mass = m_p
+stopping_on_ions_constant.background_charge_state = 1.
+stopping_on_ions_constant.background_density = 1.e20
+stopping_on_ions_constant.background_temperature = 0.05*q_e/kb # Kelvin
+
+stopping_on_ions_parsed.type = background_stopping
+stopping_on_ions_parsed.species = ions4
+stopping_on_ions_parsed.background_type = ions
+stopping_on_ions_parsed.background_mass = m_p
+stopping_on_ions_parsed.background_charge_state = 1.
+stopping_on_ions_parsed.background_density(x,y,z,t) = if(x>0,1.e20,1.e21)
+stopping_on_ions_parsed.background_temperature(x,y,z,t) = if(y>0,0.05,0.10)*q_e/kb
+
+diagnostics.diags_names = diag1
+diag1.diag_type = Full
+diag1.format = plotfile
+diag1.intervals = 1
diff --git a/Regression/Checksum/benchmarks_json/ion_stopping.json b/Regression/Checksum/benchmarks_json/ion_stopping.json
new file mode 100644
index 000000000..feb378914
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/ion_stopping.json
@@ -0,0 +1,57 @@
+{
+  "ions1": {
+    "particle_cpu": 0.0,
+    "particle_id": 10.0,
+    "particle_momentum_x": 3.4560250353184555e-22,
+    "particle_momentum_y": 3.4560250353184555e-22,
+    "particle_momentum_z": 3.4560250353184555e-22,
+    "particle_position_x": 0.0006978687369019061,
+    "particle_position_y": 0.0006978687369019061,
+    "particle_position_z": 0.0006978687369019061,
+    "particle_weight": 4.0
+  },
+  "ions2": {
+    "particle_cpu": 0.0,
+    "particle_id": 26.0,
+    "particle_momentum_x": 9.111625188610108e-22,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 0.4000131356596105,
+    "particle_position_y": 0.4,
+    "particle_position_z": 0.0,
+    "particle_weight": 4.0
+  },
+  "ions3": {
+    "particle_cpu": 0.0,
+    "particle_id": 42.0,
+    "particle_momentum_x": 6.943047764036734e-22,
+    "particle_momentum_y": 6.943047764036734e-22,
+    "particle_momentum_z": 6.943047764036734e-22,
+    "particle_position_x": 0.000349793860769889,
+    "particle_position_y": 0.000349793860769889,
+    "particle_position_z": 0.000349793860769889,
+    "particle_weight": 4.0
+  },
+  "ions4": {
+    "particle_cpu": 0.0,
+    "particle_id": 58.0,
+    "particle_momentum_x": 1.8496124473858214e-21,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 0.4000005258500433,
+    "particle_position_y": 0.4,
+    "particle_position_z": 0.0,
+    "particle_weight": 4.0
+  },
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0,
+    "Bz": 0.0,
+    "Ex": 0.0,
+    "Ey": 0.0,
+    "Ez": 0.0,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0
+  }
+}
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 61d038b8f..22452b898 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -3266,3 +3266,20 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/Tests/scraping/analysis_rz.py
+
+[ion_stopping]
+buildDir = .
+inputFile = Examples/Tests/ion_stopping/inputs_3d
+runtime_params =
+dim = 3
+addToCompileString =
+cmakeSetupOpts =
+restartTest = 0
+useMPI = 1
+numprocs = 1
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 1
+analysisRoutine = Examples/Tests/ion_stopping/analysis_ion_stopping.py
diff --git a/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.H b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.H
new file mode 100644
index 000000000..698943918
--- /dev/null
+++ b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.H
@@ -0,0 +1,77 @@
+/* Copyright 2022 David Grote
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLES_COLLISION_BACKGROUNSTOPPING_H_
+#define WARPX_PARTICLES_COLLISION_BACKGROUNSTOPPING_H_
+
+#include "Particles/MultiParticleContainer.H"
+#include "Particles/Collision/CollisionBase.H"
+
+#include <AMReX_REAL.H>
+
+#include <string>
+
+enum class BackgroundStoppingType {
+    INVALID,
+    ELECTRONS,
+    IONS,
+};
+
+class BackgroundStopping final
+    : public CollisionBase
+{
+public:
+    BackgroundStopping (std::string collision_name);
+
+    virtual ~BackgroundStopping () = default;
+
+    /** Perform the stopping calculation
+     *
+     * @param cur_time Current time
+     * @param mypc Container of species involved
+     *
+     */
+    void doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc) override;
+
+    /** Perform the stopping calculation within a tile for stopping on electrons
+     *
+     * @param pti particle iterator
+     * @param dt time step size
+     * @param t current simulation time
+     * @param species_mass mass of the active species
+     * @param species_charge charge of the active species
+     *
+     */
+    void doBackgroundStoppingOnElectronsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t,
+                                                    amrex::Real species_mass, amrex::Real species_charge);
+
+    /** Perform the stopping calculation within a tile for stopping on ions
+     *
+     * @param pti particle iterator
+     * @param dt time step size
+     * @param t current simulation time
+     * @param species_mass mass of the active species
+     * @param species_charge charge of the active species
+     *
+     */
+    void doBackgroundStoppingOnIonsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t,
+                                               amrex::Real species_mass, amrex::Real species_charge);
+
+private:
+
+    amrex::Real m_background_mass;
+    amrex::Real m_background_charge_state;
+    BackgroundStoppingType m_background_type;
+
+    amrex::Parser m_background_density_parser;
+    amrex::Parser m_background_temperature_parser;
+
+    amrex::ParserExecutor<4> m_background_density_func;
+    amrex::ParserExecutor<4> m_background_temperature_func;
+
+};
+
+#endif // WARPX_PARTICLES_COLLISION_BACKGROUNSTOPPING_H_
diff --git a/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.cpp b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.cpp
new file mode 100644
index 000000000..dd51a3d0d
--- /dev/null
+++ b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.cpp
@@ -0,0 +1,269 @@
+/* Copyright 2022 David Grote
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "BackgroundStopping.H"
+#include "Utils/ParticleUtils.H"
+#include "Utils/WarpXUtil.H"
+#include "Utils/WarpXProfilerWrapper.H"
+#include "WarpX.H"
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_REAL.H>
+
+#include <string>
+
+BackgroundStopping::BackgroundStopping (std::string const collision_name)
+    : CollisionBase(collision_name)
+{
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_species_names.size() == 1,
+                                     "Background stopping must have exactly one species.");
+
+    amrex::ParmParse pp_collision_name(collision_name);
+
+    std::string background_type_str;
+    pp_collision_name.get("background_type", background_type_str);
+    if (background_type_str == "electrons") {
+        m_background_type = BackgroundStoppingType::ELECTRONS;
+    } else if (background_type_str == "ions") {
+        m_background_type = BackgroundStoppingType::IONS;
+    } else {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "background_type must be either electrons or ions");
+    }
+
+    amrex::Real background_density;
+    std::string background_density_str;
+    if (queryWithParser(pp_collision_name, "background_density", background_density)) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(background_density > 0,
+                 "For background stopping, the background density must be greater than 0");
+        m_background_density_parser = makeParser(std::to_string(background_density), {"x", "y", "z", "t"});
+    } else if (pp_collision_name.query("background_density(x,y,z,t)", background_density_str)) {
+        m_background_density_parser = makeParser(background_density_str, {"x", "y", "z", "t"});
+    } else {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                 "For background stopping, the background density must be specified.");
+    }
+
+    amrex::Real background_temperature;
+    std::string background_temperature_str;
+    if (queryWithParser(pp_collision_name, "background_temperature", background_temperature)) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(background_temperature > 0,
+                 "For background stopping, the background temperature must be greater than 0");
+        m_background_temperature_parser = makeParser(std::to_string(background_temperature), {"x", "y", "z", "t"});
+    } else if (pp_collision_name.query("background_temperature(x,y,z,t)", background_temperature_str)) {
+        m_background_temperature_parser = makeParser(background_temperature_str, {"x", "y", "z", "t"});
+    } else {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                 "For background stopping, the background temperature must be specified.");
+    }
+
+    m_background_density_func = m_background_density_parser.compile<4>();
+    m_background_temperature_func = m_background_temperature_parser.compile<4>();
+
+    if (m_background_type == BackgroundStoppingType::ELECTRONS) {
+        m_background_mass = PhysConst::m_e;
+        queryWithParser(pp_collision_name, "background_mass", m_background_mass);
+    } else if (m_background_type == BackgroundStoppingType::IONS) {
+        getWithParser(pp_collision_name, "background_mass", m_background_mass);
+        getWithParser(pp_collision_name, "background_charge_state", m_background_charge_state);
+    }
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_background_mass > 0,
+             "For background stopping, the background mass must be greater than 0");
+
+}
+
+void
+BackgroundStopping::doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc)
+{
+    WARPX_PROFILE("BackgroundStopping::doCollisions()");
+    using namespace amrex::literals;
+
+    // Note that the lowest level time step is used for all levels since the
+    // collision operation will be done only once per step at the lowest level.
+    const amrex::Real dt = WarpX::GetInstance().getdt(0) * m_ndt;
+    if (int(std::floor(cur_time/dt)) % m_ndt != 0) return;
+
+    auto& species = mypc->GetParticleContainerFromName(m_species_names[0]);
+    amrex::Real species_mass = species.getMass();
+    amrex::Real species_charge = species.getCharge();
+
+    BackgroundStoppingType background_type = m_background_type;
+
+    // Loop over refinement levels
+    auto const flvl = species.finestLevel();
+    for (int lev = 0; lev <= flvl; ++lev) {
+
+        auto cost = WarpX::getCosts(lev);
+
+        // loop over particles box by box
+#ifdef _OPENMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (WarpXParIter pti(species, lev); pti.isValid(); ++pti) {
+            if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+            {
+                amrex::Gpu::synchronize();
+            }
+            amrex::Real wt = amrex::second();
+
+            if (background_type == BackgroundStoppingType::ELECTRONS) {
+                doBackgroundStoppingOnElectronsWithinTile(pti, dt, cur_time, species_mass, species_charge);
+            } else if (background_type == BackgroundStoppingType::IONS) {
+                doBackgroundStoppingOnIonsWithinTile(pti, dt, cur_time, species_mass, species_charge);
+            }
+
+            if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+            {
+                amrex::Gpu::synchronize();
+                wt = amrex::second() - wt;
+                amrex::HostDevice::Atomic::Add(&(*cost)[pti.index()], wt);
+            }
+        }
+
+    }
+}
+
+void BackgroundStopping::doBackgroundStoppingOnElectronsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t,
+                                                                    amrex::Real species_mass, amrex::Real species_charge)
+{
+    using namespace amrex::literals;
+
+    // So that CUDA code gets its intrinsic, not the host-only C++ library version
+    using std::sqrt, std::abs, std::log, std::exp;
+
+    // get particle count
+    const long np = pti.numParticles();
+
+    // get background particle mass
+    amrex::Real mass_e = m_background_mass;
+
+    // setup parsers for the background density and temperature
+    auto n_e_func = m_background_density_func;
+    auto T_e_func = m_background_temperature_func;
+
+    // get Struct-Of-Array particle data, also called attribs
+    auto& attribs = pti.GetAttribs();
+    amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+    amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+    amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+
+    // May be needed to evaluate the density and/or temperature functions
+    auto GetPosition = GetParticlePosition(pti);
+
+    amrex::ParallelFor(np,
+        [=] AMREX_GPU_HOST_DEVICE (long ip)
+        {
+
+            amrex::ParticleReal x, y, z;
+            GetPosition.AsStored(ip, x, y, z);
+            amrex::Real const n_e = n_e_func(x, y, z, t);
+            amrex::Real const T_e = T_e_func(x, y, z, t)*PhysConst::kb;
+
+            // This implements the equation 14.12 from Introduction to Plasma Physics,
+            // Goldston and Rutherford, the slowing down of beam ions due to collisions with electrons.
+            // The equation is written as dV/dt = -alpha*V, and integrated to
+            // give V(t+dt) = V(t)*exp(-alpha*dt)
+
+            amrex::Real constexpr pi = MathConst::pi;
+            amrex::Real constexpr ep0 = PhysConst::ep0;
+            amrex::Real constexpr q_e = PhysConst::q_e;
+            amrex::Real constexpr q_e2 = q_e*q_e;
+            amrex::Real constexpr ep02 = ep0*ep0;
+
+            amrex::Real const Zb = abs(species_charge/q_e);
+
+            amrex::Real const vth = sqrt(3._rt*T_e/mass_e);
+            amrex::Real const wp = sqrt(n_e*q_e2/(ep0*mass_e));
+            amrex::Real const lambdadb = vth/wp;
+            amrex::Real const lambdadb3 = lambdadb*lambdadb*lambdadb;
+            amrex::Real const loglambda = log((12._rt*pi/Zb)*(n_e*lambdadb3));
+
+            amrex::Real const pi32 = pi*sqrt(pi);
+            amrex::Real const q2 = species_charge*species_charge;
+            amrex::Real const T32 = T_e*sqrt(T_e);
+
+            amrex::Real const alpha = sqrt(2._rt)*n_e*q2*q_e2*sqrt(mass_e)*loglambda/(12._rt*pi32*ep02*species_mass*T32);
+
+            ux[ip] *= exp(-alpha*dt);
+            uy[ip] *= exp(-alpha*dt);
+            uz[ip] *= exp(-alpha*dt);
+
+        }
+        );
+}
+
+void BackgroundStopping::doBackgroundStoppingOnIonsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t,
+                                                               amrex::Real species_mass, amrex::Real species_charge)
+{
+    using namespace amrex::literals;
+
+    // So that CUDA code gets its intrinsic, not the host-only C++ library version
+    using std::sqrt, std::abs, std::log, std::exp, std::pow;
+
+    // get particle count
+    const long np = pti.numParticles();
+
+    // get background particle mass
+    amrex::Real mass_i = m_background_mass;
+    amrex::Real charge_state_i = m_background_charge_state;
+
+    // setup parsers for the background density and temperature
+    auto n_i_func = m_background_density_func;
+    auto T_i_func = m_background_temperature_func;
+
+    // get Struct-Of-Array particle data, also called attribs
+    auto& attribs = pti.GetAttribs();
+    amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+    amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+    amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+
+    // May be needed to evaluate the density function
+    auto GetPosition = GetParticlePosition(pti);
+
+    amrex::ParallelFor(np,
+        [=] AMREX_GPU_HOST_DEVICE (long ip)
+        {
+
+            amrex::ParticleReal x, y, z;
+            GetPosition.AsStored(ip, x, y, z);
+            amrex::Real const n_i = n_i_func(x, y, z, t);
+            amrex::Real const T_i = T_i_func(x, y, z, t)*PhysConst::kb;
+
+            // This implements the equation 14.20 from Introduction to Plasma Physics,
+            // Goldston and Rutherford, the slowing down of beam ions due to collisions with electrons.
+            // The equation is written with energy, W, as dW/dt = -alpha/W**0.5, and integrated to
+            // give W(t+dt) = (W(t)**1.5 - 3./2.*alpha*dt)**(2/3)
+
+            amrex::Real constexpr pi = MathConst::pi;
+            amrex::Real constexpr q_e = PhysConst::q_e;
+            amrex::Real constexpr q_e2 = q_e*q_e;
+            amrex::Real constexpr ep0 = PhysConst::ep0;
+            amrex::Real constexpr ep02 = ep0*ep0;
+
+            amrex::Real const qi2 = charge_state_i*charge_state_i*q_e2;
+            amrex::Real const qb2 = species_charge*species_charge;
+            amrex::Real const Zb = abs(species_charge/q_e);
+
+            amrex::Real const vth = sqrt(3._rt*T_i/mass_i);
+            amrex::Real const wp = sqrt(n_i*q_e2/(ep0*mass_i));
+            amrex::Real const lambdadb = vth/wp;
+            amrex::Real const lambdadb3 = lambdadb*lambdadb*lambdadb;
+            amrex::Real const loglambda = log((12._rt*pi/Zb)*(n_i*lambdadb3));
+
+            amrex::Real const alpha = sqrt(2._rt)*n_i*qi2*qb2*sqrt(species_mass)*loglambda/(8._rt*pi*ep02*mass_i);
+
+            amrex::Real const W0 = 0.5_rt*species_mass*(ux[ip]*ux[ip] + uy[ip]*uy[ip] + uz[ip]*uz[ip]);
+            amrex::Real const f1 = pow(W0, 1.5_rt) - 1.5_rt*alpha*dt;
+            // If f1 goes negative, the particle has fully stopped, so set W1 to 0.
+            amrex::Real const W1 = pow((f1 > 0._rt ? f1 : 0._rt), 2._rt/3._rt);
+            amrex::Real const vscale = (W0 > 0._rt ? std::sqrt(W1/W0) : 0._rt);
+
+            ux[ip] *= vscale;
+            uy[ip] *= vscale;
+            uz[ip] *= vscale;
+
+        }
+        );
+}
diff --git a/Source/Particles/Collision/BackgroundStopping/CMakeLists.txt b/Source/Particles/Collision/BackgroundStopping/CMakeLists.txt
new file mode 100644
index 000000000..fa7f65a9c
--- /dev/null
+++ b/Source/Particles/Collision/BackgroundStopping/CMakeLists.txt
@@ -0,0 +1,4 @@
+target_sources(WarpX
+  PRIVATE
+    BackgroundStopping.cpp
+)
diff --git a/Source/Particles/Collision/BackgroundStopping/Make.package b/Source/Particles/Collision/BackgroundStopping/Make.package
new file mode 100644
index 000000000..638e6f4b7
--- /dev/null
+++ b/Source/Particles/Collision/BackgroundStopping/Make.package
@@ -0,0 +1,3 @@
+CEXE_sources += BackgroundStopping.cpp
+
+VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/Collision/BackgroundStopping
diff --git a/Source/Particles/Collision/CMakeLists.txt b/Source/Particles/Collision/CMakeLists.txt
index fe7446ca1..d0348fdb4 100644
--- a/Source/Particles/Collision/CMakeLists.txt
+++ b/Source/Particles/Collision/CMakeLists.txt
@@ -7,3 +7,4 @@ target_sources(WarpX
 )
 
 add_subdirectory(BinaryCollision)
+add_subdirectory(BackgroundStopping)
diff --git a/Source/Particles/Collision/CollisionHandler.cpp b/Source/Particles/Collision/CollisionHandler.cpp
index 45d599a9e..1c1fff447 100644
--- a/Source/Particles/Collision/CollisionHandler.cpp
+++ b/Source/Particles/Collision/CollisionHandler.cpp
@@ -11,6 +11,7 @@
 #include "Particles/Collision/BinaryCollision/BinaryCollision.H"
 #include "Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H"
 #include "Particles/Collision/BinaryCollision/ParticleCreationFunc.H"
+#include "Particles/Collision/BackgroundStopping/BackgroundStopping.H"
 #include "Utils/TextMsg.H"
 
 #include <AMReX_ParmParse.H>
@@ -48,6 +49,9 @@ CollisionHandler::CollisionHandler(MultiParticleContainer const * const mypc)
         else if (type == "background_mcc") {
             allcollisions[i] = std::make_unique<BackgroundMCCCollision>(collision_names[i]);
         }
+        else if (type == "background_stopping") {
+            allcollisions[i] = std::make_unique<BackgroundStopping>(collision_names[i]);
+        }
         else if (type == "nuclearfusion") {
             allcollisions[i] =
                std::make_unique<BinaryCollision<NuclearFusionFunc, ParticleCreationFunc>>(
diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package
index f7e073ecf..330e99e08 100644
--- a/Source/Particles/Collision/Make.package
+++ b/Source/Particles/Collision/Make.package
@@ -4,5 +4,6 @@ CEXE_sources += BackgroundMCCCollision.cpp
 CEXE_sources += MCCProcess.cpp
 
 include $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision/Make.package
+include $(WARPX_HOME)/Source/Particles/Collision/BackgroundStopping/Make.package
 
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/Collision