diff --git a/epoch1d/Makefile b/epoch1d/Makefile
index b7941021a..a9077b27c 100644
--- a/epoch1d/Makefile
+++ b/epoch1d/Makefile
@@ -206,6 +206,9 @@ DEFINES := $(DEFINE)
 # Use fifth order particle weighting (default is third order).
 #DEFINES += $(D)PARTICLE_SHAPE_BSPLINE3
 
+# Use WT scheme for field interpolation
+#DEFINES += $(D)WT_INTERPOLATION
+
 # Include a unique global particle ID. The first flag defines the ID using
 # an 8-byte integer, the second uses 4-bytes.
 #DEFINES += $(D)PARTICLE_ID
diff --git a/epoch1d/src/constants.F90 b/epoch1d/src/constants.F90
index e1f5170eb..482f31778 100644
--- a/epoch1d/src/constants.F90
+++ b/epoch1d/src/constants.F90
@@ -151,6 +151,7 @@ MODULE constants
   INTEGER, PARAMETER :: c_maxwell_solver_lehe_z = 4
   INTEGER, PARAMETER :: c_maxwell_solver_cowan = 5
   INTEGER, PARAMETER :: c_maxwell_solver_pukhov = 6
+  INTEGER, PARAMETER :: c_maxwell_solver_m4 = 7
 
   ! Particle distribution type codes
   INTEGER, PARAMETER :: c_ic_df_thermal = 1
@@ -266,6 +267,7 @@ MODULE constants
   INTEGER(i8), PARAMETER :: c_def_use_isatty = 2**24
   INTEGER(i8), PARAMETER :: c_def_use_mpi3 = 2**25
   INTEGER(i8), PARAMETER :: c_def_bremsstrahlung = 2**26
+  INTEGER(i8), PARAMETER :: c_def_wt_interpolation = 2**27
 
   ! Stagger types
   INTEGER, PARAMETER :: c_stagger_ex = c_stagger_face_x
@@ -448,11 +450,12 @@ MODULE constants
   INTEGER, PARAMETER :: c_const_maxwell_solver_lehe_z = 104
   INTEGER, PARAMETER :: c_const_maxwell_solver_cowan = 105
   INTEGER, PARAMETER :: c_const_maxwell_solver_pukhov = 106
-  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_m4 = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 108
 
-  INTEGER, PARAMETER :: c_const_px = 108
-  INTEGER, PARAMETER :: c_const_py = 109
-  INTEGER, PARAMETER :: c_const_pz = 110
+  INTEGER, PARAMETER :: c_const_px = 109
+  INTEGER, PARAMETER :: c_const_py = 110
+  INTEGER, PARAMETER :: c_const_pz = 111
 
   ! Custom constants
   INTEGER, PARAMETER :: c_const_deck_lowbound = 4096
diff --git a/epoch1d/src/deck/deck_control_block.F90 b/epoch1d/src/deck/deck_control_block.F90
index 3f2f9dab1..2f26986c3 100644
--- a/epoch1d/src/deck/deck_control_block.F90
+++ b/epoch1d/src/deck/deck_control_block.F90
@@ -367,6 +367,7 @@ FUNCTION control_block_handle_element(element, value) RESULT(errcode)
           .AND. maxwell_solver /= c_maxwell_solver_lehe_x &
           .AND. maxwell_solver /= c_maxwell_solver_cowan &
           .AND. maxwell_solver /= c_maxwell_solver_pukhov &
+          .AND. maxwell_solver /= c_maxwell_solver_m4 &
           .AND. maxwell_solver /= c_maxwell_solver_custom) THEN
         errcode = c_err_bad_value
       END IF
diff --git a/epoch1d/src/fields.f90 b/epoch1d/src/fields.f90
index f485cabb9..d7647160a 100644
--- a/epoch1d/src/fields.f90
+++ b/epoch1d/src/fields.f90
@@ -59,6 +59,11 @@ SUBROUTINE set_maxwell_solver
       dx_cdt = dx / (c * dt)
       deltax = 0.25_num * (1.0_num - dx_cdt**2 * SIN(0.5_num * pi / dx_cdt)**2)
       alphax = 1.0_num - 3.0_num * deltax
+
+    ELSE IF (maxwell_solver == c_maxwell_solver_m4) THEN
+      ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+      deltax = ((c * dt / dx)**2 - 1.0) / 12.0_num
+      alphax = 1.0_num - 3.0_num * deltax
     END IF
 
     IF (rank == 0 .AND. maxwell_solver /= c_maxwell_solver_yee) THEN
diff --git a/epoch1d/src/housekeeping/setup.F90 b/epoch1d/src/housekeeping/setup.F90
index d7dc92239..f2b65006d 100644
--- a/epoch1d/src/housekeeping/setup.F90
+++ b/epoch1d/src/housekeeping/setup.F90
@@ -582,6 +582,8 @@ SUBROUTINE set_dt        ! sets CFL limited step
       ! R. Lehe, PhD Thesis (2014)
       dt = dx / c
 
+    ELSE IF (maxwell_solver == c_maxwell_solver_m4) THEN
+      dt = dx / c
     END IF
 
     IF (any_open) THEN
@@ -621,6 +623,16 @@ SUBROUTINE set_dt        ! sets CFL limited step
 
     dt = dt_multiplier * dt
 
+#ifdef WT_INTERPOLATION
+    IF (c * dt / dx > 0.5_num) THEN
+      IF (rank == 0) THEN
+        PRINT*, '*** ERROR ***'
+        PRINT*, 'Cannot use WT beacause c*dt/dx>0.5'
+      END IF
+      CALL abort_code(c_err_bad_setup)
+    END IF
+#endif
+
     IF (.NOT. any_average) RETURN
 
     DO io = 1, n_io_blocks
diff --git a/epoch1d/src/housekeeping/welcome.F90 b/epoch1d/src/housekeeping/welcome.F90
index 9e179f2a1..5696c1126 100644
--- a/epoch1d/src/housekeeping/welcome.F90
+++ b/epoch1d/src/housekeeping/welcome.F90
@@ -99,6 +99,9 @@ SUBROUTINE compiler_directives
 #ifdef PARTICLE_SHAPE_TOPHAT
     found = .TRUE.
 #endif
+#ifdef WT_INTERPOLATION
+    found = .TRUE.
+#endif
 #ifdef PER_SPECIES_WEIGHT
     found = .TRUE.
 #endif
@@ -186,6 +189,10 @@ SUBROUTINE compiler_directives
     defines = IOR(defines, c_def_particle_shape_tophat)
     WRITE(*,*) 'Top-hat particle shape -DPARTICLE_SHAPE_TOPHAT'
 #endif
+#ifdef WT_INTERPOLATION
+    defines = IOR(defines, c_def_wt_interpolation)
+    WRITE(*,*) 'WT interpolation scheme -DWT_INTERPOLATION'
+#endif
 #ifdef PER_SPECIES_WEIGHT
     defines = IOR(defines, c_def_per_particle_weight)
     WRITE(*,*) 'Per species weighting -DPER_SPECIES_WEIGHT'
diff --git a/epoch1d/src/include/bspline3/gx_wt.inc b/epoch1d/src/include/bspline3/gx_wt.inc
new file mode 100644
index 000000000..7103acafd
--- /dev/null
+++ b/epoch1d/src/include/bspline3/gx_wt.inc
@@ -0,0 +1,15 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/24
+        wt_var1 = wt_facx * (wt_dtx + cell_frac_x)**4
+        wt_var2 = wt_facx * (wt_dtx - cell_frac_x)**4
+        wt_var3 = SIGN(wt_var1, wt_dtx + cell_frac_x) &
+            + SIGN(wt_var2, wt_dtx - cell_frac_x)
+        wt_var4 = ((3.0_num - cell_frac_x) * cell_frac_x + 3.0_num - wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        wt_var5 = ((3.0_num + cell_frac_x) * cell_frac_x - 3.0_num + wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        gx(-2) = wt_var3 + wt_var1 - wt_var2
+        gx(-1) = 4.0_num * (wt_var4 - wt_var3)
+        gx( 0) = 24.0_num + 6.0_num * wt_var3 - 4.0_num * (wt_var4 + wt_var5)
+        gx( 1) = 4.0_num * (wt_var5 - wt_var3)
+        gx( 2) = wt_var3 + wt_var2 - wt_var1
diff --git a/epoch1d/src/include/tophat/gx_wt.inc b/epoch1d/src/include/tophat/gx_wt.inc
new file mode 100644
index 000000000..f76a1c8fd
--- /dev/null
+++ b/epoch1d/src/include/tophat/gx_wt.inc
@@ -0,0 +1,4 @@
+        wt_var1 = wt_facx * (ABS(wt_dtx + cell_frac_x) &
+            - ABS(wt_dtx - cell_frac_x))
+        gx( 0) = 0.5_num + wt_var1
+        gx( 1) = 0.5_num - wt_var1
diff --git a/epoch1d/src/include/triangle/gx_wt.inc b/epoch1d/src/include/triangle/gx_wt.inc
new file mode 100644
index 000000000..2ed74ffe4
--- /dev/null
+++ b/epoch1d/src/include/triangle/gx_wt.inc
@@ -0,0 +1,8 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/2
+        wt_var1 = wt_facx * ( &
+            ABS(wt_dtx + cell_frac_x) * (wt_dtx + cell_frac_x) &
+            + ABS(wt_dtx - cell_frac_x) * (wt_dtx - cell_frac_x))
+        gx(-1) = wt_var1 + cell_frac_x
+        gx( 0) = 2.0_num - 2.0_num * wt_var1
+        gx( 1) = wt_var1 - cell_frac_x
diff --git a/epoch1d/src/parser/evaluator_blocks.F90 b/epoch1d/src/parser/evaluator_blocks.F90
index a37cb6999..bd4ab2bb9 100644
--- a/epoch1d/src/parser/evaluator_blocks.F90
+++ b/epoch1d/src/parser/evaluator_blocks.F90
@@ -634,6 +634,11 @@ SUBROUTINE do_constant(opcode, simplify, parameters, err)
       RETURN
     END IF
 
+    IF (opcode == c_const_maxwell_solver_m4) THEN
+      CALL push_on_eval(REAL(c_maxwell_solver_m4, num))
+      RETURN
+    END IF
+
     IF (opcode == c_const_maxwell_solver_custom) THEN
       CALL push_on_eval(REAL(c_maxwell_solver_custom, num))
       RETURN
diff --git a/epoch1d/src/parser/tokenizer_blocks.f90 b/epoch1d/src/parser/tokenizer_blocks.f90
index aba1d582a..8f5187c0c 100644
--- a/epoch1d/src/parser/tokenizer_blocks.f90
+++ b/epoch1d/src/parser/tokenizer_blocks.f90
@@ -302,6 +302,9 @@ FUNCTION as_constant(name)
     ELSE IF (str_cmp(name, 'pukhov')) THEN
       as_constant = c_const_maxwell_solver_pukhov
 
+    ELSE IF (str_cmp(name, 'm4')) THEN
+      as_constant = c_const_maxwell_solver_m4
+
     ELSE IF (str_cmp(name, 'lehe_x')) THEN
       as_constant = c_const_maxwell_solver_lehe_x
 
diff --git a/epoch1d/src/particles.F90 b/epoch1d/src/particles.F90
index 52eae7f11..6773fcc17 100644
--- a/epoch1d/src/particles.F90
+++ b/epoch1d/src/particles.F90
@@ -82,6 +82,14 @@ SUBROUTINE push_particles
     ! This is to deal with the grid stagger
     REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx
 
+    ! For WT scheme, the weight factors used to weight particle properties onto
+    ! grid to calculate J is same as bspline3, triangle or tophat.  But the
+    ! particle weight factors to weight fields onto particles are not bspline3,
+    ! triangle or tophat. We use hx0 to store the former weight factors at Xi0*.
+#ifdef WT_INTERPOLATION
+    REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx0
+#endif
+
     ! Fields at particle location
     REAL(num) :: ex_part, ey_part, ez_part, bx_part, by_part, bz_part
 
@@ -125,6 +133,15 @@ SUBROUTINE push_particles
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_facx
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 #ifdef DELTAF_METHOD
     REAL(num) :: weight_back
@@ -141,6 +158,9 @@ SUBROUTINE push_particles
     jz = 0.0_num
 
     gx = 0.0_num
+#ifdef WT_INTERPOLATION
+    hx0 = 0.0_num
+#endif
 
     ! Unvarying multiplication factors
 
@@ -153,6 +173,14 @@ SUBROUTINE push_particles
     idtf = idt * fac
     idxf = idx * fac
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_facx = 0.25_num / wt_dtx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+#endif
+#endif
+
     DO ispecies = 1, n_species
       current => species_list(ispecies)%attached_list%head
       IF (species_list(ispecies)%immobile) CYCLE
@@ -267,6 +295,21 @@ SUBROUTINE push_particles
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+
+        ! For WT scheme, use hx0 to store the weight factors used to weigh
+        ! particle properties onto grid at Xi0*. And calculate gx for WT.
+        ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+        ! NOTE: These weights require an additional multiplication factor!
+#ifdef WT_INTERPOLATION
+        hx0 = gx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -431,6 +474,9 @@ SUBROUTINE push_particles
 
           ! Now change Xi1* to be Xi1*-Xi0*. This makes the representation of
           ! the current update much simpler
+#ifdef WT_INTERPOLATION
+          gx = hx0
+#endif
           hx = hx - gx
 
           ! Remember that due to CFL condition particle can never cross more
diff --git a/epoch1d/src/physics_packages/ionise.F90 b/epoch1d/src/physics_packages/ionise.F90
index 04ca51b07..8841d8c78 100644
--- a/epoch1d/src/physics_packages/ionise.F90
+++ b/epoch1d/src/physics_packages/ionise.F90
@@ -374,12 +374,29 @@ SUBROUTINE multiphoton_tunnelling_bsi
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_facx
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
 
     dfac = fac**2 / dt / dx
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_facx = 0.25_num / wt_dtx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -420,6 +437,7 @@ SUBROUTINE multiphoton_tunnelling_bsi
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
@@ -427,6 +445,16 @@ SUBROUTINE multiphoton_tunnelling_bsi
 #else
 #include "triangle/gx.inc"
 #endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
+#endif
+
 
         ! Now redo shifted by half a cell due to grid stagger.
         ! Use shifted version for ex in X, ey in Y, ez in Z
@@ -635,12 +663,29 @@ SUBROUTINE multiphoton_tunnelling
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_facx
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
 
     dfac = fac**2 / dt / dx
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_facx = 0.25_num / wt_dtx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -681,12 +726,22 @@ SUBROUTINE multiphoton_tunnelling
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -882,12 +937,29 @@ SUBROUTINE tunnelling_bsi
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_facx
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
 
     dfac = fac**2 / dt / dx
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_facx = 0.25_num / wt_dtx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -928,12 +1000,22 @@ SUBROUTINE tunnelling_bsi
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -1119,12 +1201,29 @@ SUBROUTINE tunnelling
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_facx
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
 
     dfac = fac**2 / dt / dx
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_facx = 0.25_num / wt_dtx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -1165,12 +1264,22 @@ SUBROUTINE tunnelling
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
diff --git a/epoch1d/src/physics_packages/photons.F90 b/epoch1d/src/physics_packages/photons.F90
index de0752e37..8fb092342 100644
--- a/epoch1d/src/physics_packages/photons.F90
+++ b/epoch1d/src/physics_packages/photons.F90
@@ -779,6 +779,24 @@ SUBROUTINE field_at_particle(part_x, e_at_part, b_at_part)
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
 #endif
 
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_facx
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
+#endif
+
+    ! Unvarying factors for WT scheme
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_facx = 0.25_num / wt_dtx
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+#endif
+#endif
     ! Grid cell position as a fraction.
 #ifdef PARTICLE_SHAPE_TOPHAT
     cell_x_r = part_x / dx - 0.5_num
@@ -796,12 +814,22 @@ SUBROUTINE field_at_particle(part_x, e_at_part, b_at_part)
     ! Also used to weight particle properties onto grid, used later
     ! to calculate J
     ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
     ! Now redo shifted by half a cell due to grid stagger.
diff --git a/epoch2d/Makefile b/epoch2d/Makefile
index 471f95233..d2dfbb65a 100644
--- a/epoch2d/Makefile
+++ b/epoch2d/Makefile
@@ -206,6 +206,9 @@ DEFINES := $(DEFINE)
 # Use fifth order particle weighting (default is third order).
 #DEFINES += $(D)PARTICLE_SHAPE_BSPLINE3
 
+# Use WT scheme for field interpolation
+#DEFINES += $(D)WT_INTERPOLATION
+
 # Include a unique global particle ID. The first flag defines the ID using
 # an 8-byte integer, the second uses 4-bytes.
 #DEFINES += $(D)PARTICLE_ID
diff --git a/epoch2d/src/constants.F90 b/epoch2d/src/constants.F90
index d43b10d39..d90495c50 100644
--- a/epoch2d/src/constants.F90
+++ b/epoch2d/src/constants.F90
@@ -151,6 +151,7 @@ MODULE constants
   INTEGER, PARAMETER :: c_maxwell_solver_lehe_z = 4
   INTEGER, PARAMETER :: c_maxwell_solver_cowan = 5
   INTEGER, PARAMETER :: c_maxwell_solver_pukhov = 6
+  INTEGER, PARAMETER :: c_maxwell_solver_m4 = 7
 
   ! Particle distribution type codes
   INTEGER, PARAMETER :: c_ic_df_thermal = 1
@@ -266,6 +267,7 @@ MODULE constants
   INTEGER(i8), PARAMETER :: c_def_use_isatty = 2**24
   INTEGER(i8), PARAMETER :: c_def_use_mpi3 = 2**25
   INTEGER(i8), PARAMETER :: c_def_bremsstrahlung = 2**26
+  INTEGER(i8), PARAMETER :: c_def_wt_interpolation = 2**27
 
   ! Stagger types
   INTEGER, PARAMETER :: c_stagger_ex = c_stagger_face_x
@@ -448,11 +450,12 @@ MODULE constants
   INTEGER, PARAMETER :: c_const_maxwell_solver_lehe_z = 104
   INTEGER, PARAMETER :: c_const_maxwell_solver_cowan = 105
   INTEGER, PARAMETER :: c_const_maxwell_solver_pukhov = 106
-  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_m4 = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 108
 
-  INTEGER, PARAMETER :: c_const_px = 108
-  INTEGER, PARAMETER :: c_const_py = 109
-  INTEGER, PARAMETER :: c_const_pz = 110
+  INTEGER, PARAMETER :: c_const_px = 109
+  INTEGER, PARAMETER :: c_const_py = 110
+  INTEGER, PARAMETER :: c_const_pz = 111
 
   ! Custom constants
   INTEGER, PARAMETER :: c_const_deck_lowbound = 4096
diff --git a/epoch2d/src/deck/deck_control_block.F90 b/epoch2d/src/deck/deck_control_block.F90
index 3822733e3..f2fa976df 100644
--- a/epoch2d/src/deck/deck_control_block.F90
+++ b/epoch2d/src/deck/deck_control_block.F90
@@ -396,6 +396,7 @@ FUNCTION control_block_handle_element(element, value) RESULT(errcode)
           .AND. maxwell_solver /= c_maxwell_solver_lehe_y &
           .AND. maxwell_solver /= c_maxwell_solver_cowan &
           .AND. maxwell_solver /= c_maxwell_solver_pukhov &
+          .AND. maxwell_solver /= c_maxwell_solver_m4 &
           .AND. maxwell_solver /= c_maxwell_solver_custom) THEN
         errcode = c_err_bad_value
       END IF
diff --git a/epoch2d/src/fields.f90 b/epoch2d/src/fields.f90
index efd7a0ee2..57893d904 100644
--- a/epoch2d/src/fields.f90
+++ b/epoch2d/src/fields.f90
@@ -85,6 +85,15 @@ SUBROUTINE set_maxwell_solver
       deltay = 0.0_num
       alphax = 1.0_num - 2.0_num * betaxy
       alphay = 1.0_num - 2.0_num * betayx
+
+    ELSE IF (maxwell_solver == c_maxwell_solver_m4) THEN
+      ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+      betayx = (c * dt / dx)**2 / 12.0_num
+      betaxy = (c * dt / dy)**2 / 12.0_num
+      deltax = ((c * dt / dx)**2 - 1.0) / 12.0_num
+      deltay = ((c * dt / dy)**2 - 1.0) / 12.0_num
+      alphax = 1.0_num - 2.0_num * betaxy - 3.0_num * deltax
+      alphay = 1.0_num - 2.0_num * betayx - 3.0_num * deltay
     END IF
 
     IF (rank == 0 .AND. maxwell_solver /= c_maxwell_solver_yee) THEN
diff --git a/epoch2d/src/housekeeping/setup.F90 b/epoch2d/src/housekeeping/setup.F90
index 613ace2ce..af63b4b3c 100644
--- a/epoch2d/src/housekeeping/setup.F90
+++ b/epoch2d/src/housekeeping/setup.F90
@@ -690,6 +690,16 @@ SUBROUTINE set_dt        ! sets CFL limited step
 
     dt = dt_multiplier * dt
 
+#ifdef WT_INTERPOLATION
+    IF (c * dt / MIN(dx, dy) > 0.5_num) THEN
+      IF (rank == 0) THEN
+        PRINT*, '*** ERROR ***'
+        PRINT*, 'Cannot use WT beacause c*dt/min(dx,dy)>0.5'
+      END IF
+      CALL abort_code(c_err_bad_setup)
+    END IF
+#endif
+
     IF (.NOT. any_average) RETURN
 
     DO io = 1, n_io_blocks
diff --git a/epoch2d/src/housekeeping/welcome.F90 b/epoch2d/src/housekeeping/welcome.F90
index 9e179f2a1..5696c1126 100644
--- a/epoch2d/src/housekeeping/welcome.F90
+++ b/epoch2d/src/housekeeping/welcome.F90
@@ -99,6 +99,9 @@ SUBROUTINE compiler_directives
 #ifdef PARTICLE_SHAPE_TOPHAT
     found = .TRUE.
 #endif
+#ifdef WT_INTERPOLATION
+    found = .TRUE.
+#endif
 #ifdef PER_SPECIES_WEIGHT
     found = .TRUE.
 #endif
@@ -186,6 +189,10 @@ SUBROUTINE compiler_directives
     defines = IOR(defines, c_def_particle_shape_tophat)
     WRITE(*,*) 'Top-hat particle shape -DPARTICLE_SHAPE_TOPHAT'
 #endif
+#ifdef WT_INTERPOLATION
+    defines = IOR(defines, c_def_wt_interpolation)
+    WRITE(*,*) 'WT interpolation scheme -DWT_INTERPOLATION'
+#endif
 #ifdef PER_SPECIES_WEIGHT
     defines = IOR(defines, c_def_per_particle_weight)
     WRITE(*,*) 'Per species weighting -DPER_SPECIES_WEIGHT'
diff --git a/epoch2d/src/include/bspline3/gx_wt.inc b/epoch2d/src/include/bspline3/gx_wt.inc
new file mode 100644
index 000000000..98ba952ea
--- /dev/null
+++ b/epoch2d/src/include/bspline3/gx_wt.inc
@@ -0,0 +1,29 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/24
+        wt_var1 = wt_facx * (wt_dtx + cell_frac_x)**4
+        wt_var2 = wt_facx * (wt_dtx - cell_frac_x)**4
+        wt_var3 = SIGN(wt_var1, wt_dtx + cell_frac_x) &
+            + SIGN(wt_var2, wt_dtx - cell_frac_x)
+        wt_var4 = ((3.0_num - cell_frac_x) * cell_frac_x + 3.0_num - wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        wt_var5 = ((3.0_num + cell_frac_x) * cell_frac_x - 3.0_num + wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        gx(-2) = wt_var3 + wt_var1 - wt_var2
+        gx(-1) = 4.0_num * (wt_var4 - wt_var3)
+        gx( 0) = 24.0_num + 6.0_num * wt_var3 - 4.0_num * (wt_var4 + wt_var5)
+        gx( 1) = 4.0_num * (wt_var5 - wt_var3)
+        gx( 2) = wt_var3 + wt_var2 - wt_var1
+
+        wt_var1 = wt_facy * (wt_dty + cell_frac_y)**4
+        wt_var2 = wt_facy * (wt_dty - cell_frac_y)**4
+        wt_var3 = SIGN(wt_var1, wt_dty + cell_frac_y) &
+            + SIGN(wt_var2, wt_dty - cell_frac_y)
+        wt_var4 = ((3.0_num - cell_frac_y) * cell_frac_y + 3.0_num - wt_dty2) &
+            * cell_frac_y + 1.0_num + wt_dty2
+        wt_var5 = ((3.0_num + cell_frac_y) * cell_frac_y - 3.0_num + wt_dty2) &
+            * cell_frac_y + 1.0_num + wt_dty2
+        gy(-2) = wt_var3 + wt_var1 - wt_var2
+        gy(-1) = 4.0_num * (wt_var4 - wt_var3)
+        gy( 0) = 24.0_num + 6.0_num * wt_var3 - 4.0_num * (wt_var4 + wt_var5)
+        gy( 1) = 4.0_num * (wt_var5 - wt_var3)
+        gy( 2) = wt_var3 + wt_var2 - wt_var1
diff --git a/epoch2d/src/include/tophat/gx_wt.inc b/epoch2d/src/include/tophat/gx_wt.inc
new file mode 100644
index 000000000..5bae9d300
--- /dev/null
+++ b/epoch2d/src/include/tophat/gx_wt.inc
@@ -0,0 +1,9 @@
+        wt_var1 = wt_facx * (ABS(wt_dtx + cell_frac_x) &
+            - ABS(wt_dtx - cell_frac_x))
+        gx( 0) = 0.5_num + wt_var1
+        gx( 1) = 0.5_num - wt_var1
+
+        wt_var1 = wt_facy * (ABS(wt_dty + cell_frac_y) &
+            - ABS(wt_dty - cell_frac_y))
+        gy( 0) = 0.5_num + wt_var1
+        gy( 1) = 0.5_num - wt_var1
diff --git a/epoch2d/src/include/triangle/gx_wt.inc b/epoch2d/src/include/triangle/gx_wt.inc
new file mode 100644
index 000000000..db1649a8e
--- /dev/null
+++ b/epoch2d/src/include/triangle/gx_wt.inc
@@ -0,0 +1,15 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/2
+        wt_var1 = wt_facx * ( &
+            ABS(wt_dtx + cell_frac_x) * (wt_dtx + cell_frac_x) &
+            + ABS(wt_dtx - cell_frac_x) * (wt_dtx - cell_frac_x))
+        gx(-1) = wt_var1 + cell_frac_x
+        gx( 0) = 2.0_num - 2.0_num * wt_var1
+        gx( 1) = wt_var1 - cell_frac_x
+
+        wt_var1 = wt_facy * ( &
+            ABS(wt_dty + cell_frac_y) * (wt_dty + cell_frac_y) &
+            + ABS(wt_dty - cell_frac_y) * (wt_dty - cell_frac_y))
+        gy(-1) = wt_var1 + cell_frac_y
+        gy( 0) = 2.0_num - 2.0_num * wt_var1
+        gy( 1) = wt_var1 - cell_frac_y
diff --git a/epoch2d/src/parser/evaluator_blocks.F90 b/epoch2d/src/parser/evaluator_blocks.F90
index 345926956..6005f7e05 100644
--- a/epoch2d/src/parser/evaluator_blocks.F90
+++ b/epoch2d/src/parser/evaluator_blocks.F90
@@ -654,6 +654,11 @@ SUBROUTINE do_constant(opcode, simplify, parameters, err)
       RETURN
     END IF
 
+    IF (opcode == c_const_maxwell_solver_m4) THEN
+      CALL push_on_eval(REAL(c_maxwell_solver_m4, num))
+      RETURN
+    END IF
+
     IF (opcode == c_const_maxwell_solver_custom) THEN
       CALL push_on_eval(REAL(c_maxwell_solver_custom, num))
       RETURN
diff --git a/epoch2d/src/parser/tokenizer_blocks.f90 b/epoch2d/src/parser/tokenizer_blocks.f90
index 4450d28b9..60df4705a 100644
--- a/epoch2d/src/parser/tokenizer_blocks.f90
+++ b/epoch2d/src/parser/tokenizer_blocks.f90
@@ -339,6 +339,9 @@ FUNCTION as_constant(name)
     ELSE IF (str_cmp(name, 'pukhov')) THEN
       as_constant = c_const_maxwell_solver_pukhov
 
+    ELSE IF (str_cmp(name, 'm4')) THEN
+      as_constant = c_const_maxwell_solver_m4
+
     ELSE IF (str_cmp(name, 'lehe_x')) THEN
       as_constant = c_const_maxwell_solver_lehe_x
 
diff --git a/epoch2d/src/particles.F90 b/epoch2d/src/particles.F90
index de257e55d..6972a730b 100644
--- a/epoch2d/src/particles.F90
+++ b/epoch2d/src/particles.F90
@@ -85,6 +85,14 @@ SUBROUTINE push_particles
     ! This is to deal with the grid stagger
     REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx, hy
 
+    ! For WT scheme, the weight factors used to weight particle properties onto
+    ! grid to calculate J is same as bspline3, triangle or tophat.  But the
+    ! particle weight factors to weight fields onto particles are not bspline3,
+    ! triangle or tophat. We use hx0 to store the former weight factors at Xi0*.
+#ifdef WT_INTERPOLATION
+    REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx0, hy0
+#endif
+
     ! Fields at particle location
     REAL(num) :: ex_part, ey_part, ez_part, bx_part, by_part, bz_part
 
@@ -129,6 +137,15 @@ SUBROUTINE push_particles
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_facx, wt_facy
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 #ifdef DELTAF_METHOD
     REAL(num) :: weight_back
@@ -146,6 +163,10 @@ SUBROUTINE push_particles
 
     gx = 0.0_num
     gy = 0.0_num
+#ifdef WT_INTERPOLATION
+    hx0 = 0.0_num
+    hy0 = 0.0_num
+#endif
 
     ! Unvarying multiplication factors
 
@@ -161,6 +182,17 @@ SUBROUTINE push_particles
     idtx = idt * idx * fac
     idxy = idx * idy * fac
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+#endif
+#endif
+
     DO ispecies = 1, n_species
       current => species_list(ispecies)%attached_list%head
       IF (species_list(ispecies)%immobile) CYCLE
@@ -286,6 +318,22 @@ SUBROUTINE push_particles
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+
+        ! For WT scheme, use hx0 to store the weight factors used to weigh
+        ! particle properties onto grid at Xi0*. And calculate gx for WT.
+        ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+        ! NOTE: These weights require an additional multiplication factor!
+#ifdef WT_INTERPOLATION
+        hx0 = gx
+        hy0 = gy
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -469,6 +517,10 @@ SUBROUTINE push_particles
 
           ! Now change Xi1* to be Xi1*-Xi0*. This makes the representation of
           ! the current update much simpler
+#ifdef WT_INTERPOLATION
+          gx = hx0
+          gy = hy0
+#endif
           hx = hx - gx
           hy = hy - gy
 
diff --git a/epoch2d/src/physics_packages/ionise.F90 b/epoch2d/src/physics_packages/ionise.F90
index 0e3c4f5fa..4c1fb103b 100644
--- a/epoch2d/src/physics_packages/ionise.F90
+++ b/epoch2d/src/physics_packages/ionise.F90
@@ -411,6 +411,15 @@ SUBROUTINE multiphoton_tunnelling_bsi
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_facx, wt_facy
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -418,6 +427,17 @@ SUBROUTINE multiphoton_tunnelling_bsi
 
     dfac = fac**2 / dt / dx / dy
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -465,12 +485,22 @@ SUBROUTINE multiphoton_tunnelling_bsi
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -693,6 +723,15 @@ SUBROUTINE multiphoton_tunnelling
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_facx, wt_facy
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -700,6 +739,17 @@ SUBROUTINE multiphoton_tunnelling
 
     dfac = fac**2 / dt / dx / dy
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -747,12 +797,22 @@ SUBROUTINE multiphoton_tunnelling
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -961,6 +1021,15 @@ SUBROUTINE tunnelling_bsi
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_facx, wt_facy
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -968,6 +1037,17 @@ SUBROUTINE tunnelling_bsi
 
     dfac = fac**2 / dt / dx / dy
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -1037,12 +1117,22 @@ SUBROUTINE tunnelling_bsi
         dcellx = 0
         dcelly = 0
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/hx_dcell.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/hx_dcell.inc"
 #else
 #include "triangle/hx_dcell.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! These are the electric and magnetic fields interpolated to the
@@ -1219,6 +1309,15 @@ SUBROUTINE tunnelling
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_facx, wt_facy
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -1226,6 +1325,17 @@ SUBROUTINE tunnelling
 
     dfac = fac**2 / dt / dx / dy
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -1273,12 +1383,22 @@ SUBROUTINE tunnelling
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
diff --git a/epoch2d/src/physics_packages/photons.F90 b/epoch2d/src/physics_packages/photons.F90
index f43978bd2..2873bb10b 100644
--- a/epoch2d/src/physics_packages/photons.F90
+++ b/epoch2d/src/physics_packages/photons.F90
@@ -779,6 +779,27 @@ SUBROUTINE field_at_particle(part_x, part_y, e_at_part, b_at_part)
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_facx, wt_facy
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
+#endif
+
+    ! Unvarying factors for WT scheme
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+#endif
 #endif
 
     ! Grid cell position as a fraction.
@@ -804,12 +825,22 @@ SUBROUTINE field_at_particle(part_x, part_y, e_at_part, b_at_part)
     ! Also used to weight particle properties onto grid, used later
     ! to calculate J
     ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
     ! Now redo shifted by half a cell due to grid stagger.
diff --git a/epoch3d/Makefile b/epoch3d/Makefile
index afbd8b51c..4e4842eae 100644
--- a/epoch3d/Makefile
+++ b/epoch3d/Makefile
@@ -206,6 +206,9 @@ DEFINES := $(DEFINE)
 # Use fifth order particle weighting (default is third order).
 #DEFINES += $(D)PARTICLE_SHAPE_BSPLINE3
 
+# Use WT scheme for field interpolation
+#DEFINES += $(D)WT_INTERPOLATION
+
 # Include a unique global particle ID. The first flag defines the ID using
 # an 8-byte integer, the second uses 4-bytes.
 #DEFINES += $(D)PARTICLE_ID
diff --git a/epoch3d/src/constants.F90 b/epoch3d/src/constants.F90
index 5ecd7f6a2..655d6756e 100644
--- a/epoch3d/src/constants.F90
+++ b/epoch3d/src/constants.F90
@@ -151,6 +151,7 @@ MODULE constants
   INTEGER, PARAMETER :: c_maxwell_solver_lehe_z = 4
   INTEGER, PARAMETER :: c_maxwell_solver_cowan = 5
   INTEGER, PARAMETER :: c_maxwell_solver_pukhov = 6
+  INTEGER, PARAMETER :: c_maxwell_solver_m4 = 7
 
   ! Particle distribution type codes
   INTEGER, PARAMETER :: c_ic_df_thermal = 1
@@ -266,6 +267,7 @@ MODULE constants
   INTEGER(i8), PARAMETER :: c_def_use_isatty = 2**24
   INTEGER(i8), PARAMETER :: c_def_use_mpi3 = 2**25
   INTEGER(i8), PARAMETER :: c_def_bremsstrahlung = 2**26
+  INTEGER(i8), PARAMETER :: c_def_wt_interpolation = 2**27
 
   ! Stagger types
   INTEGER, PARAMETER :: c_stagger_ex = c_stagger_face_x
@@ -448,11 +450,12 @@ MODULE constants
   INTEGER, PARAMETER :: c_const_maxwell_solver_lehe_z = 104
   INTEGER, PARAMETER :: c_const_maxwell_solver_cowan = 105
   INTEGER, PARAMETER :: c_const_maxwell_solver_pukhov = 106
-  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_m4 = 107
+  INTEGER, PARAMETER :: c_const_maxwell_solver_custom = 108
 
-  INTEGER, PARAMETER :: c_const_px = 108
-  INTEGER, PARAMETER :: c_const_py = 109
-  INTEGER, PARAMETER :: c_const_pz = 110
+  INTEGER, PARAMETER :: c_const_px = 109
+  INTEGER, PARAMETER :: c_const_py = 110
+  INTEGER, PARAMETER :: c_const_pz = 111
 
   ! Custom constants
   INTEGER, PARAMETER :: c_const_deck_lowbound = 4096
diff --git a/epoch3d/src/deck/deck_control_block.F90 b/epoch3d/src/deck/deck_control_block.F90
index 002749cf8..ef66479d4 100644
--- a/epoch3d/src/deck/deck_control_block.F90
+++ b/epoch3d/src/deck/deck_control_block.F90
@@ -415,6 +415,7 @@ FUNCTION control_block_handle_element(element, value) RESULT(errcode)
           .AND. maxwell_solver /= c_maxwell_solver_lehe_z &
           .AND. maxwell_solver /= c_maxwell_solver_cowan &
           .AND. maxwell_solver /= c_maxwell_solver_pukhov &
+          .AND. maxwell_solver /= c_maxwell_solver_m4 &
           .AND. maxwell_solver /= c_maxwell_solver_custom) THEN
         errcode = c_err_bad_value
       END IF
diff --git a/epoch3d/src/fields.f90 b/epoch3d/src/fields.f90
index 6d49781a7..51ca6645a 100644
--- a/epoch3d/src/fields.f90
+++ b/epoch3d/src/fields.f90
@@ -162,6 +162,30 @@ SUBROUTINE set_maxwell_solver
       alphax = 1.0_num - 2.0_num * betaxy - 2.0_num * betaxz
       alphay = 1.0_num - 2.0_num * betayx - 2.0_num * betayz
       alphaz = 1.0_num - 2.0_num * betazx - 2.0_num * betazy
+
+    ELSE IF (maxwell_solver == c_maxwell_solver_m4) THEN
+      ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+      c1 = (c * dt / dx)**2
+      c2 = (c * dt / dy)**2
+      c3 = (c * dt / dz)**2
+      betayx = c1 / 12.0_num
+      betaxy = c2 / 12.0_num
+      betaxz = c3 / 12.0_num
+      betazx = betayx
+      betazy = betaxy
+      betayz = betaxz
+      deltax = (c1 - 1.0) / 12.0_num
+      deltay = (c2 - 1.0) / 12.0_num
+      deltaz = (c3 - 1.0) / 12.0_num
+      gammax = 0.0_num
+      gammay = 0.0_num
+      gammaz = 0.0_num
+      alphax = 1.0_num - 2.0_num * betaxy - 2.0_num * betaxz &
+                       - 4.0_num * gammax - 3.0_num * deltax
+      alphay = 1.0_num - 2.0_num * betayx - 2.0_num * betayz &
+                       - 4.0_num * gammay - 3.0_num * deltay
+      alphaz = 1.0_num - 2.0_num * betazx - 2.0_num * betazy &
+                       - 4.0_num * gammaz - 3.0_num * deltaz
     END IF
 
     IF (rank == 0 .AND. maxwell_solver /= c_maxwell_solver_yee) THEN
diff --git a/epoch3d/src/housekeeping/setup.F90 b/epoch3d/src/housekeeping/setup.F90
index 4b282b591..aacaa9eeb 100644
--- a/epoch3d/src/housekeeping/setup.F90
+++ b/epoch3d/src/housekeeping/setup.F90
@@ -721,6 +721,7 @@ SUBROUTINE set_dt        ! sets CFL limited step
       dt = MIN(dz, dx * dy / SQRT(dx**2 + dy**2)) / c
 
     ELSE IF (maxwell_solver == c_maxwell_solver_cowan &
+        .OR. maxwell_solver == c_maxwell_solver_m4 &
         .OR. maxwell_solver == c_maxwell_solver_pukhov) THEN
       ! Cowan et al., Phys. Rev. ST Accel. Beams 16, 041303 (2013)
       ! A. Pukhov, Journal of Plasma Physics 61, 425-433 (1999)
@@ -764,6 +765,16 @@ SUBROUTINE set_dt        ! sets CFL limited step
 
     dt = dt_multiplier * dt
 
+#ifdef WT_INTERPOLATION
+    IF (c * dt / MIN(dx, dy, dx) > 0.5_num) THEN
+      IF (rank == 0) THEN
+        PRINT*, '*** ERROR ***'
+        PRINT*, 'Cannot use WT beacause c*dt/min(dx,dy,dx)>0.5'
+      END IF
+      CALL abort_code(c_err_bad_setup)
+    END IF
+#endif
+
     IF (.NOT. any_average) RETURN
 
     DO io = 1, n_io_blocks
diff --git a/epoch3d/src/housekeeping/welcome.F90 b/epoch3d/src/housekeeping/welcome.F90
index 9e179f2a1..5696c1126 100644
--- a/epoch3d/src/housekeeping/welcome.F90
+++ b/epoch3d/src/housekeeping/welcome.F90
@@ -99,6 +99,9 @@ SUBROUTINE compiler_directives
 #ifdef PARTICLE_SHAPE_TOPHAT
     found = .TRUE.
 #endif
+#ifdef WT_INTERPOLATION
+    found = .TRUE.
+#endif
 #ifdef PER_SPECIES_WEIGHT
     found = .TRUE.
 #endif
@@ -186,6 +189,10 @@ SUBROUTINE compiler_directives
     defines = IOR(defines, c_def_particle_shape_tophat)
     WRITE(*,*) 'Top-hat particle shape -DPARTICLE_SHAPE_TOPHAT'
 #endif
+#ifdef WT_INTERPOLATION
+    defines = IOR(defines, c_def_wt_interpolation)
+    WRITE(*,*) 'WT interpolation scheme -DWT_INTERPOLATION'
+#endif
 #ifdef PER_SPECIES_WEIGHT
     defines = IOR(defines, c_def_per_particle_weight)
     WRITE(*,*) 'Per species weighting -DPER_SPECIES_WEIGHT'
diff --git a/epoch3d/src/include/bspline3/gx_wt.inc b/epoch3d/src/include/bspline3/gx_wt.inc
new file mode 100644
index 000000000..8901664b3
--- /dev/null
+++ b/epoch3d/src/include/bspline3/gx_wt.inc
@@ -0,0 +1,43 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/24
+        wt_var1 = wt_facx * (wt_dtx + cell_frac_x)**4
+        wt_var2 = wt_facx * (wt_dtx - cell_frac_x)**4
+        wt_var3 = SIGN(wt_var1, wt_dtx + cell_frac_x) &
+            + SIGN(wt_var2, wt_dtx - cell_frac_x)
+        wt_var4 = ((3.0_num - cell_frac_x) * cell_frac_x + 3.0_num - wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        wt_var5 = ((3.0_num + cell_frac_x) * cell_frac_x - 3.0_num + wt_dtx2) &
+            * cell_frac_x + 1.0_num + wt_dtx2
+        gx(-2) = wt_var3 + wt_var1 - wt_var2
+        gx(-1) = 4.0_num * (wt_var4 - wt_var3)
+        gx( 0) = 24.0_num + 6.0_num * wt_var3 - 4.0_num * (wt_var4 + wt_var5)
+        gx( 1) = 4.0_num * (wt_var5 - wt_var3)
+        gx( 2) = wt_var3 + wt_var2 - wt_var1
+
+        wt_var1 = wt_facy * (wt_dty + cell_frac_y)**4
+        wt_var2 = wt_facy * (wt_dty - cell_frac_y)**4
+        wt_var3 = SIGN(wt_var1, wt_dty + cell_frac_y) &
+            + SIGN(wt_var2, wt_dty - cell_frac_y)
+        wt_var4 = ((3.0_num - cell_frac_y) * cell_frac_y + 3.0_num - wt_dty2) &
+            * cell_frac_y + 1.0_num + wt_dty2
+        wt_var5 = ((3.0_num + cell_frac_y) * cell_frac_y - 3.0_num + wt_dty2) &
+            * cell_frac_y + 1.0_num + wt_dty2
+        gy(-2) = wt_var3 + wt_var1 - wt_var2
+        gy(-1) = 4.0_num * (wt_var4 - wt_var3)
+        gy( 0) = 24.0_num + 6.0_num * wt_var3 - 4.0_num * (wt_var4 + wt_var5)
+        gy( 1) = 4.0_num * (wt_var5 - wt_var3)
+        gy( 2) = wt_var3 + wt_var2 - wt_var1
+
+        wt_var1 = wt_facz * (wt_dtz + cell_frac_z)**4
+        wt_var2 = wt_facz * (wt_dtz - cell_frac_z)**4
+        wt_var3 = SIGN(wt_var1, wt_dtz + cell_frac_z) &
+            + SIGN(wt_var2, wt_dtz - cell_frac_z)
+        wt_var4 = ((3.0_num - cell_frac_z) * cell_frac_z + 3.0_num - wt_dtz2) &
+            * cell_frac_z + 1.0_num + wt_dtz2
+        wt_var5 = ((3.0_num + cell_frac_z) * cell_frac_z - 3.0_num + wt_dtz2) &
+            * cell_frac_z + 1.0_num + wt_dtz2
+        gz(-2) = wt_var3 + wt_var1 - wt_var2
+        gz(-1) = 4.0_num * (wt_var4 - wt_var3)
+        gz( 0) = 24.0_num + 6.0_num * wt_var3 - 4.0_num * (wt_var4 + wt_var5)
+        gz( 1) = 4.0_num * (wt_var5 - wt_var3)
+        gz( 2) = wt_var3 + wt_var2 - wt_var1
diff --git a/epoch3d/src/include/tophat/gx_wt.inc b/epoch3d/src/include/tophat/gx_wt.inc
new file mode 100644
index 000000000..97371a922
--- /dev/null
+++ b/epoch3d/src/include/tophat/gx_wt.inc
@@ -0,0 +1,14 @@
+        wt_var1 = wt_facx * (ABS(wt_dtx + cell_frac_x) &
+            - ABS(wt_dtx - cell_frac_x))
+        gx( 0) = 0.5_num + wt_var1
+        gx( 1) = 0.5_num - wt_var1
+
+        wt_var1 = wt_facy * (ABS(wt_dty + cell_frac_y) &
+            - ABS(wt_dty - cell_frac_y))
+        gy( 0) = 0.5_num + wt_var1
+        gy( 1) = 0.5_num - wt_var1
+
+        wt_var1 = wt_facz * (ABS(wt_dtz + cell_frac_z) &
+            - ABS(wt_dtz - cell_frac_z))
+        gz( 0) = 0.5_num + wt_var1
+        gz( 1) = 0.5_num - wt_var1
diff --git a/epoch3d/src/include/triangle/gx_wt.inc b/epoch3d/src/include/triangle/gx_wt.inc
new file mode 100644
index 000000000..19f878b7b
--- /dev/null
+++ b/epoch3d/src/include/triangle/gx_wt.inc
@@ -0,0 +1,22 @@
+! IMPORTANT NOTE
+! These weight need to be multiplied by 1/2
+        wt_var1 = wt_facx * ( &
+            ABS(wt_dtx + cell_frac_x) * (wt_dtx + cell_frac_x) &
+            + ABS(wt_dtx - cell_frac_x) * (wt_dtx - cell_frac_x))
+        gx(-1) = wt_var1 + cell_frac_x
+        gx( 0) = 2.0_num - 2.0_num * wt_var1
+        gx( 1) = wt_var1 - cell_frac_x
+
+        wt_var1 = wt_facy * ( &
+            ABS(wt_dty + cell_frac_y) * (wt_dty + cell_frac_y) &
+            + ABS(wt_dty - cell_frac_y) * (wt_dty - cell_frac_y))
+        gy(-1) = wt_var1 + cell_frac_y
+        gy( 0) = 2.0_num - 2.0_num * wt_var1
+        gy( 1) = wt_var1 - cell_frac_y
+
+        wt_var1 = wt_facz * ( &
+            ABS(wt_dtz + cell_frac_z) * (wt_dtz + cell_frac_z) &
+            + ABS(wt_dtz - cell_frac_z) * (wt_dtz - cell_frac_z))
+        gz(-1) = wt_var1 + cell_frac_z
+        gz( 0) = 2.0_num - 2.0_num * wt_var1
+        gz( 1) = wt_var1 - cell_frac_z
diff --git a/epoch3d/src/parser/evaluator_blocks.F90 b/epoch3d/src/parser/evaluator_blocks.F90
index 9efab12ab..473085f65 100644
--- a/epoch3d/src/parser/evaluator_blocks.F90
+++ b/epoch3d/src/parser/evaluator_blocks.F90
@@ -669,6 +669,11 @@ SUBROUTINE do_constant(opcode, simplify, parameters, err)
       RETURN
     END IF
 
+    IF (opcode == c_const_maxwell_solver_m4) THEN
+      CALL push_on_eval(REAL(c_maxwell_solver_m4, num))
+      RETURN
+    END IF
+
     IF (opcode == c_const_maxwell_solver_custom) THEN
       CALL push_on_eval(REAL(c_maxwell_solver_custom, num))
       RETURN
diff --git a/epoch3d/src/parser/tokenizer_blocks.f90 b/epoch3d/src/parser/tokenizer_blocks.f90
index 29f434e4c..dbb351552 100644
--- a/epoch3d/src/parser/tokenizer_blocks.f90
+++ b/epoch3d/src/parser/tokenizer_blocks.f90
@@ -382,6 +382,9 @@ FUNCTION as_constant(name)
     ELSE IF (str_cmp(name, 'pukhov')) THEN
       as_constant = c_const_maxwell_solver_pukhov
 
+    ELSE IF (str_cmp(name, 'm4')) THEN
+      as_constant = c_const_maxwell_solver_m4
+
     ELSE IF (str_cmp(name, 'lehe_x')) THEN
       as_constant = c_const_maxwell_solver_lehe_x
 
diff --git a/epoch3d/src/particles.F90 b/epoch3d/src/particles.F90
index 97c2f2ced..d6a294f5c 100644
--- a/epoch3d/src/particles.F90
+++ b/epoch3d/src/particles.F90
@@ -88,6 +88,14 @@ SUBROUTINE push_particles
     ! This is to deal with the grid stagger
     REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx, hy, hz
 
+    ! For WT scheme, the weight factors used to weight particle properties onto
+    ! grid to calculate J is same as bspline3, triangle or tophat.  But the
+    ! particle weight factors to weight fields onto particles are not bspline3,
+    ! triangle or tophat. We use hx0 to store the former weight factors at Xi0*.
+#ifdef WT_INTERPOLATION
+    REAL(num), DIMENSION(sf_min-1:sf_max+1) :: hx0, hy0, hz0
+#endif
+
     ! Fields at particle location
     REAL(num) :: ex_part, ey_part, ez_part, bx_part, by_part, bz_part
 
@@ -133,6 +141,15 @@ SUBROUTINE push_particles
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_dtz, wt_facx, wt_facy, wt_facz
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2, wt_dtz2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 #ifdef DELTAF_METHOD
     REAL(num) :: weight_back
@@ -151,6 +168,11 @@ SUBROUTINE push_particles
     gx = 0.0_num
     gy = 0.0_num
     gz = 0.0_num
+#ifdef WT_INTERPOLATION
+    hx0 = 0.0_num
+    hy0 = 0.0_num
+    hz0 = 0.0_num
+#endif
 
     ! Unvarying multiplication factors
 
@@ -167,6 +189,20 @@ SUBROUTINE push_particles
     idtxz = idt * idx * idz * fac
     idtxy = idt * idx * idy * fac
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_dtz = c * dt / dz
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+    wt_facz = 0.25_num / wt_dtz
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+    wt_dtz2 = wt_dtz**2
+#endif
+#endif
+
     DO ispecies = 1, n_species
       current => species_list(ispecies)%attached_list%head
       IF (species_list(ispecies)%immobile) CYCLE
@@ -301,6 +337,23 @@ SUBROUTINE push_particles
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+
+        ! For WT scheme, use hx0 to store the weight factors used to weigh
+        ! particle properties onto grid at Xi0*. And calculate gx for WT.
+        ! Y. Lu et al., J. Comput. Phys 413, 109388 (2020)
+        ! NOTE: These weights require an additional multiplication factor!
+#ifdef WT_INTERPOLATION
+        hx0 = gx
+        hy0 = gy
+        hz0 = gz
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -499,6 +552,11 @@ SUBROUTINE push_particles
 
           ! Now change Xi1* to be Xi1*-Xi0*. This makes the representation of
           ! the current update much simpler
+#ifdef WT_INTERPOLATION
+          gx = hx0
+          gy = hy0
+          gz = hz0
+#endif
           hx = hx - gx
           hy = hy - gy
           hz = hz - gz
diff --git a/epoch3d/src/physics_packages/ionise.F90 b/epoch3d/src/physics_packages/ionise.F90
index e5b076ab6..d0e8bae58 100644
--- a/epoch3d/src/physics_packages/ionise.F90
+++ b/epoch3d/src/physics_packages/ionise.F90
@@ -448,6 +448,15 @@ SUBROUTINE multiphoton_tunnelling_bsi
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_dtz, wt_facx, wt_facy, wt_facz
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2, wt_dtz2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -456,6 +465,20 @@ SUBROUTINE multiphoton_tunnelling_bsi
 
     dfac = fac**2 / dt / dx / dy / dz
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_dtz = c * dt / dz
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+    wt_facz = 0.25_num / wt_dtz
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+    wt_dtz2 = wt_dtz**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -510,12 +533,22 @@ SUBROUTINE multiphoton_tunnelling_bsi
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -751,6 +784,15 @@ SUBROUTINE multiphoton_tunnelling
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_dtz, wt_facx, wt_facy, wt_facz
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2, wt_dtz2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -759,6 +801,20 @@ SUBROUTINE multiphoton_tunnelling
 
     dfac = fac**2 / dt / dx / dy / dz
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_dtz = c * dt / dz
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+    wt_facz = 0.25_num / wt_dtz
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+    wt_dtz2 = wt_dtz**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -813,12 +869,22 @@ SUBROUTINE multiphoton_tunnelling
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -1040,6 +1106,15 @@ SUBROUTINE tunnelling_bsi
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_dtz, wt_facx, wt_facy, wt_facz
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2, wt_dtz2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -1048,6 +1123,20 @@ SUBROUTINE tunnelling_bsi
 
     dfac = fac**2 / dt / dx / dy / dz
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_dtz = c * dt / dz
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+    wt_facz = 0.25_num / wt_dtz
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+    wt_dtz2 = wt_dtz**2
+#endif
+#endif
+ 
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -1102,12 +1191,22 @@ SUBROUTINE tunnelling_bsi
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
@@ -1319,6 +1418,15 @@ SUBROUTINE tunnelling
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_dtz, wt_facx, wt_facy, wt_facz
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2, wt_dtz2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
     idx = 1.0_num / dx
@@ -1327,6 +1435,20 @@ SUBROUTINE tunnelling
 
     dfac = fac**2 / dt / dx / dy / dz
 
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_dtz = c * dt / dz
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+    wt_facz = 0.25_num / wt_dtz
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+    wt_dtz2 = wt_dtz**2
+#endif
+#endif
+
     ! Stores ionised species until close of ionisation run. Main purpose of this
     ! method is to ensure proper statistics (i.e. prevent ionisation rate being
     ! calculated for full dt when particle has already ionised in time step)
@@ -1381,12 +1503,22 @@ SUBROUTINE tunnelling
         ! Also used to weight particle properties onto grid, used later
         ! to calculate J
         ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
         ! Now redo shifted by half a cell due to grid stagger.
diff --git a/epoch3d/src/physics_packages/photons.F90 b/epoch3d/src/physics_packages/photons.F90
index e1b7aaeb8..e317469f2 100644
--- a/epoch3d/src/physics_packages/photons.F90
+++ b/epoch3d/src/physics_packages/photons.F90
@@ -781,8 +781,31 @@ SUBROUTINE field_at_particle(part_x, part_y, part_z, e_at_part, b_at_part)
 #else
     REAL(num) :: cf2
     REAL(num), PARAMETER :: fac = (0.5_num)**c_ndims
+#endif
+    ! Factors for WT scheme
+#ifdef WT_INTERPOLATION
+    REAL(num) :: wt_dtx, wt_dty, wt_dtz, wt_facx, wt_facy, wt_facz
+    REAL(num) :: wt_var1
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    REAL(num) :: wt_dtx2, wt_dty2, wt_dtz2
+    REAL(num) :: wt_var2, wt_var3, wt_var4, wt_var5
+#endif
 #endif
 
+    ! Unvarying factors for WT scheme
+#ifdef WT_INTERPOLATION
+    wt_dtx = c * dt / dx
+    wt_dty = c * dt / dy
+    wt_dtz = c * dt / dz
+    wt_facx = 0.25_num / wt_dtx
+    wt_facy = 0.25_num / wt_dty
+    wt_facz = 0.25_num / wt_dtz
+#ifdef PARTICLE_SHAPE_BSPLINE3
+    wt_dtx2 = wt_dtx**2
+    wt_dty2 = wt_dty**2
+    wt_dtz2 = wt_dtz**2
+#endif
+#endif
     ! Grid cell position as a fraction.
 #ifdef PARTICLE_SHAPE_TOPHAT
     cell_x_r = part_x / dx - 0.5_num
@@ -812,12 +835,22 @@ SUBROUTINE field_at_particle(part_x, part_y, part_z, e_at_part, b_at_part)
     ! Also used to weight particle properties onto grid, used later
     ! to calculate J
     ! NOTE: These weights require an additional multiplication factor!
+#ifndef WT_INTERPOLATION
 #ifdef PARTICLE_SHAPE_BSPLINE3
 #include "bspline3/gx.inc"
 #elif  PARTICLE_SHAPE_TOPHAT
 #include "tophat/gx.inc"
 #else
 #include "triangle/gx.inc"
+#endif
+#else
+#ifdef PARTICLE_SHAPE_BSPLINE3
+#include "bspline3/gx_wt.inc"
+#elif  PARTICLE_SHAPE_TOPHAT
+#include "tophat/gx_wt.inc"
+#else
+#include "triangle/gx_wt.inc"
+#endif
 #endif
 
     ! Now redo shifted by half a cell due to grid stagger.