Why does this Non-Standard FDTD implementation lead to infinite increase in the magnitude of an EM pulse?

I have been working on a Particle-In-Cell Framework in Python and have noticed an issue where the magnitude of a EM pulse increasing infinitely as the simulation updates.

Currently, I am using the Non-Standard Finite Difference algorithm outlined in this paper: J.L.Vay 2008.

Here is the algorithm used to advance both the electric and magnetic field in vacuum:

@n.njit(parallel=True)
def E_Update(E, B, J, ds, dt, R):
C = 2.99792458e8
EPSILON = 8.8541878128e-12
MU = 1.25663706212e-6
for i in n.prange(1,E.shape[0]-1):
for j in n.prange(1,E.shape[1]-1):
for k in n.prange(1,E.shape[2]-1):
# x-component update
# -> Solve DyBz
Bz_y0 = B[i,j-1,k,2]
Bz_y1 = B[i,j,k,2]
DyBz = (Bz_y1-Bz_y0)/ds
# -> Solve DzBy
By_z0 = B[i,j,k-1,1]
By_z1 = B[i,j,k,1]
DzBy = (By_z1-By_z0)/ds
# -> Solve Ex
R[i,j,k,0] = (1/EPSILON)*((DyBz-DzBy) - J[i,j,k,0])*dt + E[i,j,k,0]

# y-component update
# -> Solve DxBz
Bz_x0 = B[i-1,j,k,2]
Bz_x1 = B[i,j,k,2]
DxBz = (Bz_x1-Bz_x0)/ds
# -> Solve DzBx
Bx_z0 = B[i,j,k-1,0]
Bx_z1 = B[i,j,k,0]
DzBx = (Bx_z1-Bx_z0)/ds
# -> Solve Ey
R[i,j,k,1] = (1/EPSILON)*((-1)*(DyBz-DzBy) - J[i,j,k,1])*dt + E[i,j,k,1]

# z-component update
# -> Solve DxBy
By_x0 = B[i-1,j,k,1]
By_x1 = B[i,j,k,1]
DxBy = (By_x1-By_x0)/ds
# -> Solve DyBx
Bx_y0 = B[i,j-1,k,0]
Bx_y1 = B[i,j,k,0]
DyBx = (Bx_y1-Bx_y0)/ds
# -> Solve Ez
R[i,j,k,2] = (1/EPSILON)*((DxBy-DyBx) - J[i,j,k,2])*dt + E[i,j,k,2]

@n.njit(parallel=True)
def B_Update(E, B, ds, dt, a, b, g, R):
C = 2.99792458e8
EPSILON = 8.8541878128e-12
MU = 1.25663706212e-6
C1 = ((1-(dt/(2*MU)))/(1+(dt/(2*MU))))
C2 = (1/(1+(dt/(2*MU))))
for i in n.prange(1,E.shape[0]-1):
for j in n.prange(1,E.shape[1]-1):
for k in n.prange(1,E.shape[2]-1):

# x-component update
# -> Solve DzEy
Ey_z0 = E[i,j,k,1]
Ey_z1 = E[i,j,k+1,1]
DzEy = (Ey_z1-Ey_z0)/ds
# -> Solve S1zEy
S1zEyQ1 = E[i+1,j,k,1] # Ey(i+1,j,k)
S1zEyQ2 = E[i-1,j,k,1] # Ey(i-1,j,k)
S1zEyQ3 = E[i,j+1,k,1] # Ey(i,j+1,k)
S1zEyQ4 = E[i,j-1,k,1] # Ey(i,j-1,k)
S1zEy = S1zEyQ1 + S1zEyQ2 + S1zEyQ3 + S1zEyQ4
# -> Solve S2zEy
S2zEyQ1 = E[i+1,j+1,k,1] # Ey[i+1,j+1,k]
S2zEyQ2 = E[i+1,j-1,k,1] # Ey[i+1,j-1,k]
S2zEyQ3 = E[i-1,j+1,k,1] # Ey[i-1,j+1,k]
S2zEyQ4 = E[i-1,j-1,k,1] # Ey[i-1,j-1,k]
S2zEy = S2zEyQ1 + S2zEyQ2 + S2zEyQ3 + S2zEyQ4
# -> Solve DSzEy
DSzEy = (a+b*S1zEy+g*S2zEy)*(DzEy)
# -> Solve DyEz
Ez_y0 = E[i,j,k,2]
Ez_y1 = E[i,j+1,k,2]
DyEz = (Ez_y1-Ez_y0)/ds
# -> Solve S1yEz
S1yEzQ1 = E[i+1,j,k,2] # Ez(i+1,j,k)
S1yEzQ2 = E[i-1,j,k,2] # Ez(i-1,j,k)
S1yEzQ3 = E[i,j,k+1,2] # Ez(i,j,k+1)
S1yEzQ4 = E[i,j,k-1,2] # Ez(i,j,k-1)
S1yEz = S1yEzQ1 + S1yEzQ2 + S1yEzQ3 + S1yEzQ4
# -> Solve S2yEz
S2yEzQ1 = E[i+1,j,k+1,2] # Ez(i+1,j,k+1)
S2yEzQ2 = E[i+1,j,k-1,2] # Ez(i+1,j,k-1)
S2yEzQ3 = E[i-1,j,k+1,2] # Ez(i-1,j,k+1)
S2yEzQ4 = E[i-1,j,k-1,2] # Ez(i-1,j,k-1)
S2yEz = S2yEzQ1 + S2yEzQ2 + S2yEzQ3 + S2yEzQ4
# -> Solve DSyEz
DSyEz = (a+b*S1yEz+g*S2yEz)*(DyEz)
# -> Solve DtBx
R[i,j,k,0] = (1/MU)*((DSzEy-DSyEz))*dt + B[i,j,k,0]

# y-component update
# -> Solve DzEx
Ex_z0 = E[i,j,k,0]
Ex_z1 = E[i,j,k+1,0]
DzEx = (Ex_z1-Ex_z0)/ds
# -> Solve S1zEx
S1zExQ1 = E[i+1,j,k,0] # Ex(i+1,j,k)
S1zExQ2 = E[i-1,j,k,0] # Ex(i-1,j,k)
S1zExQ3 = E[i,j+1,k,0] # Ex(i,j+1,k)
S1zExQ4 = E[i,j-1,k,0] # Ex(i,j-1,k)
S1zEx = S1zExQ1 + S1zExQ2 + S1zExQ3 + S1zExQ4
# -> Solve S2zEx
S2zExQ1 = E[i+1,j+1,k,0] # Ex(i+1,j+1,k)
S2zExQ2 = E[i+1,j-1,k,0] # Ex(i+1,j-1,k)
S2zExQ3 = E[i-1,j+1,k,0] # Ex(i-1,j+1,k)
S2zExQ4 = E[i-1,j-1,k,0] # Ex(i-1,j-1,k)
S2zEx = S2zExQ1 + S2zExQ2 + S2zExQ3 + S2zExQ4
# -> Solve DSzEx
DSzEx = (a+b*S1zEx+g*S2zEx)*(DzEx)
# -> Solve DxEz
Ez_x0 = E[i,j,k,2]
Ez_x1 = E[i+1,j,k,2]
DxEz = (Ez_x1-Ez_x0)/ds
# -> Solve S1xEz
S1xEzQ1 = E[i,j+1,k,2] # Ez(i,j+1,k)
S1xEzQ2 = E[i,j-1,k,2] # Ez(i,j-1,k)
S1xEzQ3 = E[i,j,k+1,2] # Ez(i,j,k+1)
S1xEzQ4 = E[i,j,k-1,2] # Ez(i,j,k-1)
S1xEz = S1xEzQ1 + S1xEzQ2 + S1xEzQ3 + S1xEzQ4
# -> Solve S2xEz
S2xEzQ1 = E[i,j+1,k+1,2] # Ez(i,j+1,k+1)
S2xEzQ2 = E[i,j+1,k-1,2] # Ez(i,j+1,k-1)
S2xEzQ3 = E[i,j-1,k+1,2] # Ez(i,j-1,k+1)
S2xEzQ4 = E[i,j-1,k-1,2] # Ez(i,j-1,k-1)
S2xEz = S2xEzQ1+S2xEzQ2+S2xEzQ3+S2xEzQ4
# -> Solve DSxEz
DSxEz = (a+b*S1xEz+g*S2xEz)*(DxEz)
# -> Solve DtBy
R[i,j,k,1] = (1/MU)*((DSxEz-DSzEx))*dt + B[i,j,k,1]

# z-component update
# -> Solve DyEx
Ex_y0 = E[i,j,k,0]
Ex_y1 = E[i,j+1,k,0]
DyEx = (Ex_y1-Ex_y0)/ds
# -> Solve S1yEx
S1yExQ1 = E[i+1,j,k,0] # Ex(i+1,j,k)
S1yExQ2 = E[i-1,j,k,0] # Ex(i-1,j,k)
S1yExQ3 = E[i,j,k+1,0] # Ex(i,j,k+1)
S1yExQ4 = E[i,j,k-1,0] # Ex(i,j,k-1)
S1yEx = S1yExQ1 + S1yExQ2 + S1yExQ3 + S1yExQ4
# -> Solve S2yEx
S2yExQ1 = E[i+1,j,k+1,0] # Ex(i+1,j,k+1)
S2yExQ2 = E[i+1,j,k-1,0] # Ex(i+1,j,k-1)
S2yExQ3 = E[i-1,j,k+1,0] # Ex(i-1,j,k+1)
S2yExQ4 = E[i-1,j,k-1,0] # Ex(i-1,j,k-1)
S2yEx = S2yExQ1 + S2yExQ2 + S2yExQ3 + S2yExQ4
# -> Solve DSyEx
DSyEx = (a+b*S1yEx+g*S2yEx)*(DyEx)
# -> Solve DxEy
Ey_x0 = E[i,j,k,1]
Ey_x1 = E[i+1,j,k,1]
DxEy = (Ey_x1-Ey_x0)/ds
# -> Solve S1xEy
S1xEyQ1 = E[i,j+1,k,1] # Ey(i,j+1,k)
S1xEyQ2 = E[i,j-1,k,1] # Ey(i,j-1,k)
S1xEyQ3 = E[i,j,k+1,1] # Ey(i,j,k+1)
S1xEyQ4 = E[i,j,k-1,1] # Ey(i,j,k-1)
S1xEy = S1xEyQ1 + S1xEyQ2 + S1xEyQ3 + S1xEyQ4
# -> Solve S2xEy
S2xEyQ1 = E[i,j+1,k+1,1] # Ey(i,j+1,k+1)
S2xEyQ2 = E[i,j+1,k-1,1] # Ey(i,j+1,k-1)
S2xEyQ3 = E[i,j-1,k+1,1] # Ey(i,j-1,k+1)
S2xEyQ4 = E[i,j-1,k-1,1] # Ey(i,j-1,k-1)
S2xEy = S2xEyQ1 + S2xEyQ2 + S2xEyQ3 + S2xEyQ4
# -> Solve DSxEy
DSxEy = (a+b*S1xEy+g*S2xEy)*(DxEy)
# -> Solve DtBz
R[i,j,k,2] = (1/MU)*((DSyEx-DSxEy))*dt + B[i,j,k,2]


Each update step involves advancing the B-Field by $$\Delta t / 2$$, then the E-Field by $$\Delta t$$, and finally the B-Field once more by $$\Delta t / 2$$.

So far, I've tried to resolve the issue by:

1. Reverting to standard FDTD (i.e. setting a,b, and g in the B_Update method to 1,0, and 0 respectively)
2. Ensuring the CFL Condition is met by keeping the time-step several orders of magnitude lower than traversal time of light across the length of a Yee-Cell.
3. Using natural units for the simulation ($$c = \mu = \epsilon = 1)$$

The only resource I've found that touches on the issue of FDTD simulations "exploding" is a YouTube lecture about the basics of the algorithm. The stability of the FDTD method depends on the existence of certain values used in the update equations at the proper time and spatial coordinates. I could be wrong, but I believe that the update functions I've written meet this requirement for stability.

Any suggested changes or resource recommendations would be greatly appreciated.