diff --git a/Examples/Tests/particle_thermal_boundary/analysis_2d.py b/Examples/Tests/particle_thermal_boundary/analysis_2d.py
new file mode 100755
index 00000000000..f5edc775cba
--- /dev/null
+++ b/Examples/Tests/particle_thermal_boundary/analysis_2d.py
@@ -0,0 +1,40 @@
+#!/use/bin/env python3
+#
+# Copyright 2023-2024 Revathi Jambunathan
+#
+# This file is part of WarpX
+#
+# License: BSD-3-Clause-LBNL
+
+"""
+The script checks to see the growth in total field energy and total particle energy
+
+The input file in 2D initializes uniform plasma (electrons,ions) with
+thermal boundary condition. We do not expect the particle energy to increase
+beyond 2% in the time that it takes all particles to cross the domain boundary
+"""
+
+
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(1,'../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+filename = sys.argv[1]
+
+FE_rdiag = '/diags/reducedfiles/EF.txt'
+init_Fenergy = np.loadtxt(FE_rdiag)[1,2]
+final_Fenergy = np.loadtxt(FE_rdiag)[-1,2]
+assert(final_Fenergy/init_Fenergy < 40)
+assert(final_Fenergy < 5.e-5)
+
+PE_rdiag = '/diags/reducedfiles/EN.txt'
+init_Penergy = np.loadtxt(PE_rdiag)[0,2]
+final_Penergy = np.loadtxt(PE_rdiag)[-1,2]
+assert( abs(final_Penergy - init_Penergy)/init_Penergy < 0.02)
+
+test_name = os.path.splot(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename)
diff --git a/Examples/Tests/particle_thermal_boundary/inputs_2d b/Examples/Tests/particle_thermal_boundary/inputs_2d
new file mode 100644
index 00000000000..9c3e18cd35d
--- /dev/null
+++ b/Examples/Tests/particle_thermal_boundary/inputs_2d
@@ -0,0 +1,77 @@
+max_step = 2000
+
+# number of grid points
+amr.n_cell =   16 16
+
+# Maximum allowable size of each subdomain in the problem domain;
+#    this is used to decompose the domain for parallel calculations.
+amr.max_grid_size = 128
+
+# Maximum level in hierarchy (for now must be 0, i.e., one level in total)
+amr.max_level = 0
+
+# Geometry
+geometry.dims = 2
+geometry.prob_lo     = 0.e-6   0.e-6    # physical domain
+geometry.prob_hi     = 2.5e-7   2.5e-7
+
+# Boundary condition
+boundary.field_lo = pml pml
+boundary.field_hi = pml pml
+boundary.particle_lo = thermal thermal
+boundary.particle_hi = thermal thermal
+boundary.electrons.u_th = uth_e
+boundary.C.u_th = uth_C
+warpx.do_dive_cleaning = 0
+
+# Verbosity
+warpx.verbose = 1
+
+# Order of particle shape factors
+algo.particle_shape = 2
+
+# CFL
+warpx.cfl = 0.98
+
+# Density
+my_constants.n0 = 4e26
+my_constants.uth_e = 0.06256112470898544
+my_constants.uth_C = 0.00042148059678527106
+
+# Particles
+particles.species_names = electrons C
+
+electrons.charge = -q_e
+electrons.mass = m_e
+electrons.injection_style = "NUniformPerCell"
+electrons.num_particles_per_cell_each_dim = 8 8
+electrons.profile = constant
+electrons.density = n0
+electrons.momentum_distribution_type = gaussian
+electrons.ux_th = uth_e
+electrons.uy_th = uth_e
+electrons.uz_th = uth_e
+
+C.charge = 6*q_e
+C.mass = 12*m_p
+C.injection_style = "NUniformPerCell"
+C.num_particles_per_cell_each_dim = 8 8
+C.profile = constant
+C.density = n0/6
+C.momentum_distribution_type = gaussian
+C.ux_th = uth_C
+C.uy_th = uth_C
+C.uz_th = uth_C
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 3000
+diag1.write_species = 0
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz rho divE
+diag1.diag_type = Full
+
+warpx.reduced_diags_names = EN EF
+EN.intervals = 10
+EN.type = ParticleEnergy
+EF.intervals = 10
+EF.type = FieldEnergy
diff --git a/Regression/Checksum/benchmarks_json/particle_thermal_boundary.json b/Regression/Checksum/benchmarks_json/particle_thermal_boundary.json
new file mode 100644
index 00000000000..5678a294eff
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/particle_thermal_boundary.json
@@ -0,0 +1,12 @@
+{
+  "lev=0": {
+    "Bx": 229.28366410112736,
+    "By": 242.40221902487008,
+    "Bz": 164.70543954019666,
+    "Ex": 789234700854.9321,
+    "Ey": 15892271613.51144,
+    "Ez": 693631132513.1528,
+    "divE": 5.441439446237374e+19,
+    "rho": 331367728.98455656
+  }
+}
\ No newline at end of file
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 56a7fbc8f10..d9160bebaa7 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -4720,3 +4720,20 @@ compareParticles = 1
 particleTypes = electrons
 outputFile = particle_boundary_interaction_plt
 analysisRoutine = Examples/Tests/particle_boundary_interaction/analysis.py
+
+[particle_thermal_boundary]
+buildDir = .
+inputFile = Examples/Tests/particle_thermal_boundary/inputs_2d
+runtime_params =
+dim = 2
+addToCompileString =
+cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_OPENPMD=ON
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/particle_thermal_boundary/analysis.py