diff --git a/lasy/backend.py b/lasy/backend.py
new file mode 100644
index 00000000..f47858fc
--- /dev/null
+++ b/lasy/backend.py
@@ -0,0 +1,10 @@
+try:
+    import cupy as xp
+
+    use_cupy = True
+except ImportError:
+    import numpy as xp
+
+    use_cupy = False
+
+__all__ = ["use_cupy", "xp"]
diff --git a/lasy/laser.py b/lasy/laser.py
index 7f4f2ef8..68b4f14f 100644
--- a/lasy/laser.py
+++ b/lasy/laser.py
@@ -10,6 +10,8 @@
 )
 from lasy.utils.openpmd_output import write_to_openpmd_file
 
+from .backend import use_cupy, xp
+
 
 class Laser:
     """
@@ -87,7 +89,7 @@ class Laser:
     >>>     extent[2:] *= 1e6
     >>>     extent[:2] *= 1e12
     >>>     tmin, tmax, rmin, rmax = extent
-    >>>     vmax = np.abs(E_rt).max()
+    >>>     vmax = xp.abs(E_rt).max()
     >>>     axes[step].imshow(
     ...         E_rt,
     ...         origin="lower",
@@ -113,11 +115,11 @@ def __init__(
         # Get the spectral axis
         dt = self.grid.dx[time_axis_indx]
         Nt = self.grid.shape[time_axis_indx]
-        self.omega_1d = 2 * np.pi * np.fft.fftfreq(Nt, dt) + profile.omega0
+        self.omega_1d = 2 * xp.pi * xp.fft.fftfreq(Nt, dt) + profile.omega0
 
         # Create the grid on which to evaluate the laser, evaluate it
         if self.dim == "xyt":
-            x, y, t = np.meshgrid(*self.grid.axes, indexing="ij")
+            x, y, t = xp.meshgrid(*self.grid.axes, indexing="ij")
             self.grid.set_temporal_field(profile.evaluate(x, y, t))
         elif self.dim == "rt":
             if n_theta_evals is None:
@@ -126,17 +128,17 @@ def __init__(
                 n_theta_evals = 2 * self.grid.n_azimuthal_modes - 1
             # Make sure that there are enough points to resolve the azimuthal modes
             assert n_theta_evals >= 2 * self.grid.n_azimuthal_modes - 1
-            theta1d = 2 * np.pi / n_theta_evals * np.arange(n_theta_evals)
-            theta, r, t = np.meshgrid(theta1d, *self.grid.axes, indexing="ij")
-            x = r * np.cos(theta)
-            y = r * np.sin(theta)
+            theta1d = 2 * xp.pi / n_theta_evals * xp.arange(n_theta_evals)
+            theta, r, t = xp.meshgrid(theta1d, *self.grid.axes, indexing="ij")
+            x = r * xp.cos(theta)
+            y = r * xp.sin(theta)
             # Evaluate the profile on the generated grid
             envelope = profile.evaluate(x, y, t)
             # Perform the azimuthal decomposition
-            azimuthal_modes = np.fft.ifft(envelope, axis=0)
+            azimuthal_modes = xp.fft.ifft(envelope, axis=0)
             field = azimuthal_modes[:n_azimuthal_modes]
             if n_azimuthal_modes > 1:
-                field = np.concatenate(
+                field = xp.concatenate(
                     (field, azimuthal_modes[-n_azimuthal_modes + 1 :])
                 )
             self.grid.set_temporal_field(field)
@@ -179,7 +181,8 @@ def apply_optics(self, optical_element):
         # Apply optical element
         spectral_field = self.grid.get_spectral_field()
         if self.dim == "rt":
-            r, omega = np.meshgrid(self.grid.axes[0], self.omega_1d, indexing="ij")
+            r, omega = xp.meshgrid(self.grid.axes[0], self.omega_1d, indexing="ij")
+
             # The line below assumes that amplitude_multiplier
             # is cylindrically symmetric, hence we pass
             # `r` as `x` and 0 as `y`
@@ -194,7 +197,7 @@ def apply_optics(self, optical_element):
             for i_m in range(self.grid.azimuthal_modes.size):
                 spectral_field[i_m, :, :] *= multiplier
         else:
-            x, y, omega = np.meshgrid(
+            x, y, omega = xp.meshgrid(
                 self.grid.axes[0], self.grid.axes[1], self.omega_1d, indexing="ij"
             )
             spectral_field *= optical_element.amplitude_multiplier(
@@ -202,7 +205,7 @@ def apply_optics(self, optical_element):
             )
         self.grid.set_spectral_field(spectral_field)
 
-    def propagate(self, distance, nr_boundary=None, backend="NP", show_progress=True):
+    def propagate(self, distance, nr_boundary=None, show_progress=True):
         """
         Propagate the laser pulse by the distance specified.
 
@@ -216,16 +219,14 @@ def propagate(self, distance, nr_boundary=None, backend="NP", show_progress=True
             will be attenuated (to assert proper Hankel transform).
             Only used for ``'rt'``.
 
-        backend : string (optional)
-            Backend used by axiprop (see axiprop documentation).
         show_progress : bool (optional)
             Whether to show a progress bar when performing the computation
         """
         # apply boundary "absorption" if required
         if nr_boundary is not None:
             assert type(nr_boundary) is int and nr_boundary > 0
-            absorb_layer_axis = np.linspace(0, np.pi / 2, nr_boundary)
-            absorb_layer_shape = np.cos(absorb_layer_axis) ** 0.5
+            absorb_layer_axis = xp.linspace(0, xp.pi / 2, nr_boundary)
+            absorb_layer_shape = xp.cos(absorb_layer_axis) ** 0.5
             absorb_layer_shape[-1] = 0.0
             field = self.grid.get_temporal_field()
             if self.dim == "rt":
@@ -237,16 +238,27 @@ def propagate(self, distance, nr_boundary=None, backend="NP", show_progress=True
                 field[:, :nr_boundary, :] *= absorb_layer_shape[::-1][None, :, None]
             self.grid.set_temporal_field(field)
 
+        # Select backend
+        if use_cupy:
+            backend = "CU"
+        else:
+            backend = "NP"
+
+        k = self.omega_1d / c
         if self.dim == "rt":
             # Construct the propagator (check if exists)
             if not hasattr(self, "prop"):
                 spatial_axes = (self.grid.axes[0],)
                 self.prop = []
+                if use_cupy:
+                    # Move quantities to CPU to create propagator
+                    k = xp.asnumpy(k)
+                    spatial_axes = (xp.asnumpy(spatial_axes[0]),)
                 for m in self.grid.azimuthal_modes:
                     self.prop.append(
                         PropagatorResampling(
                             *spatial_axes,
-                            self.omega_1d / c,
+                            k,
                             mode=m,
                             backend=backend,
                             verbose=False,
@@ -255,14 +267,14 @@ def propagate(self, distance, nr_boundary=None, backend="NP", show_progress=True
             # Propagate the spectral image
             spectral_field = self.grid.get_spectral_field()
             for i_m in range(self.grid.azimuthal_modes.size):
-                transform_data = np.transpose(spectral_field[i_m]).copy()
+                transform_data = xp.transpose(spectral_field[i_m]).copy()
                 self.prop[i_m].step(
                     transform_data,
                     distance,
                     overwrite=True,
                     show_progress=show_progress,
                 )
-                spectral_field[i_m, :, :] = np.transpose(transform_data).copy()
+                spectral_field[i_m, :, :] = xp.transpose(transform_data).copy()
             self.grid.set_spectral_field(spectral_field)
         else:
             # Construct the propagator (check if exists)
@@ -279,11 +291,11 @@ def propagate(self, distance, nr_boundary=None, backend="NP", show_progress=True
                 )
             # Propagate the spectral image
             spectral_field = self.grid.get_spectral_field()
-            transform_data = np.moveaxis(spectral_field, -1, 0).copy()
+            transform_data = xp.moveaxis(spectral_field, -1, 0).copy()
             self.prop.step(
                 transform_data, distance, overwrite=True, show_progress=show_progress
             )
-            spectral_field = np.moveaxis(transform_data, 0, -1).copy()
+            spectral_field = xp.moveaxis(transform_data, 0, -1).copy()
 
         # Choose the time translation assuming propagation at v=c
         translate_time = distance / c
@@ -293,7 +305,7 @@ def propagate(self, distance, nr_boundary=None, backend="NP", show_progress=True
         # propagators, so it needs to be added by hand.
         # Note: subtracting by omega0 is only a global phase convention,
         # that derives from the definition of the envelope in lasy.
-        spectral_field *= np.exp(
+        spectral_field *= xp.exp(
             -1j * (self.omega_1d[None, None, :] - self.profile.omega0) * translate_time
         )
         self.grid.set_spectral_field(spectral_field)
@@ -349,7 +361,8 @@ def show(self, **kw):
         ----------
         **kw : additional arguments to be passed to matplotlib's imshow command
         """
-        temporal_field = self.grid.get_temporal_field()
+        # Get field on CPU
+        temporal_field = self.grid.get_temporal_field(to_cpu=True)
         if self.dim == "rt":
             # Show field in the plane y=0, above and below axis, with proper sign for each mode
             E = [
diff --git a/lasy/optical_elements/optical_element.py b/lasy/optical_elements/optical_element.py
index cd7e5d50..92f08d58 100644
--- a/lasy/optical_elements/optical_element.py
+++ b/lasy/optical_elements/optical_element.py
@@ -1,6 +1,6 @@
 from abc import ABC, abstractmethod
 
-import numpy as np
+from lasy.backend import xp
 
 
 class OpticalElement(ABC):
@@ -48,4 +48,4 @@ def amplitude_multiplier(self, x, y, omega, omega0):
         """
         # The base class only defines dummy multiplier
         # (This should be replaced by any class that inherits from this one.)
-        return np.ones_like(x, dtype="complex128")
+        return xp.ones_like(x, dtype="complex128")
diff --git a/lasy/optical_elements/parabolic_mirror.py b/lasy/optical_elements/parabolic_mirror.py
index 54ede494..143bf766 100644
--- a/lasy/optical_elements/parabolic_mirror.py
+++ b/lasy/optical_elements/parabolic_mirror.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.constants import c
 
+from lasy.backend import xp
+
 from .optical_element import OpticalElement
 
 
@@ -47,4 +48,4 @@ def amplitude_multiplier(self, x, y, omega, omega0):
             Contains the value of the multiplier at the specified points.
             This array has the same shape as the array omega.
         """
-        return np.exp(-1j * omega * (x**2 + y**2) / (2 * c * self.f))
+        return xp.exp(-1j * omega * (x**2 + y**2) / (2 * c * self.f))
diff --git a/lasy/profiles/from_array_profile.py b/lasy/profiles/from_array_profile.py
index 1586581f..d3cd161d 100644
--- a/lasy/profiles/from_array_profile.py
+++ b/lasy/profiles/from_array_profile.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.interpolate import RegularGridInterpolator
 
+from lasy.backend import xp
+
 from .profile import Profile
 
 
@@ -56,13 +57,13 @@ def __init__(self, wavelength, pol, array, dim, axes, axes_order=["x", "y", "t"]
             assert axes_order in [["x", "y", "t"], ["t", "y", "x"]]
 
             if axes_order == ["t", "y", "x"]:
-                self.array = np.swapaxes(array, 0, 2)
+                self.array = xp.swapaxes(array, 0, 2)
             else:
                 self.array = array
 
             self.combined_field_interp = RegularGridInterpolator(
                 (axes["x"], axes["y"], axes["t"]),
-                np.abs(array) + 1.0j * np.unwrap(np.angle(array), axis=-1),
+                xp.abs(array) + 1.0j * xp.unwrap(xp.angle(array), axis=-1),
                 bounds_error=False,
                 fill_value=0.0,
             )
@@ -70,21 +71,21 @@ def __init__(self, wavelength, pol, array, dim, axes, axes_order=["x", "y", "t"]
             assert axes_order in [["r", "t"], ["t", "r"]]
 
             if axes_order == ["t", "r"]:
-                self.array = np.swapaxes(array, 0, 2)
+                self.array = xp.swapaxes(array, 0, 2)
             else:
                 self.array = array
 
             # If the first point of radial axis is not 0, we "mirror" it,
             # to make correct interpolation within the first cell
             if axes["r"][0] != 0.0:
-                r = np.concatenate(([-axes["r"][0]], axes["r"]))
-                array = np.concatenate(([array[0]], array))
+                r = xp.concatenate(([-axes["r"][0]], axes["r"]))
+                array = xp.concatenate(([array[0]], array))
             else:
                 r = axes["r"]
 
             self.combined_field_interp = RegularGridInterpolator(
                 (r, axes["t"]),
-                np.abs(array) + 1.0j * np.unwrap(np.angle(array), axis=-1),
+                xp.abs(array) + 1.0j * xp.unwrap(xp.angle(array), axis=-1),
                 bounds_error=False,
                 fill_value=0.0,
             )
@@ -94,10 +95,10 @@ def evaluate(self, x, y, t):
         if self.dim == "xyt":
             combined_field = self.combined_field_interp((x, y, t))
         else:
-            combined_field = self.combined_field_interp((np.sqrt(x**2 + y**2), t))
+            combined_field = self.combined_field_interp((xp.sqrt(x**2 + y**2), t))
 
-        envelope = np.abs(np.real(combined_field)) * np.exp(
-            1.0j * np.imag(combined_field)
+        envelope = xp.abs(xp.real(combined_field)) * xp.exp(
+            1.0j * xp.imag(combined_field)
         )
 
         return envelope
diff --git a/lasy/profiles/from_insight_file.py b/lasy/profiles/from_insight_file.py
index bc832335..aa09047f 100644
--- a/lasy/profiles/from_insight_file.py
+++ b/lasy/profiles/from_insight_file.py
@@ -1,7 +1,8 @@
 import h5py
-import numpy as np
 from scipy.constants import c
 
+from lasy.backend import xp
+
 from .from_array_profile import FromArrayProfile
 
 
@@ -30,10 +31,10 @@ class FromInsightFile(FromArrayProfile):
     def __init__(self, file_path, pol, omega0="barycenter"):
         # read the data from H5 filed
         with h5py.File(file_path, "r") as hf:
-            data = np.asanyarray(hf["data/Exyt_0"][()], dtype=np.complex128, order="C")
-            t = np.asanyarray(hf["scales/t"][()], dtype=np.float64, order="C")
-            x = np.asanyarray(hf["scales/x"][()], dtype=np.float64, order="C")
-            y = np.asanyarray(hf["scales/y"][()], dtype=np.float64, order="C")
+            data = xp.asanyarray(hf["data/Exyt_0"][()], dtype=xp.complex128, order="C")
+            t = xp.asanyarray(hf["scales/t"][()], dtype=xp.float64, order="C")
+            x = xp.asanyarray(hf["scales/x"][()], dtype=xp.float64, order="C")
+            y = xp.asanyarray(hf["scales/y"][()], dtype=xp.float64, order="C")
 
         # convert data and axes to SI units
         t *= 1e-15
@@ -42,29 +43,29 @@ def __init__(self, file_path, pol, omega0="barycenter"):
 
         # get the field on axis and local frequencies
         field_onaxis = data[data.shape[0] // 2, data.shape[1] // 2, :]
-        omega_array = -np.gradient(np.unwrap(np.angle(field_onaxis)), t)
+        omega_array = -xp.gradient(xp.unwrap(xp.angle(field_onaxis)), t)
 
         # choose the central frequency
         if omega0 == "peak":
             # using peak field frequency
-            omega0 = omega_array[np.abs(field_onaxis).argmax()]
+            omega0 = omega_array[xp.abs(field_onaxis).argmax()]
         elif omega0 == "barycenter":
             # or "center of mass" frequency
-            omega0 = np.average(omega_array, weights=np.abs(field_onaxis) ** 2)
+            omega0 = xp.average(omega_array, weights=xp.abs(field_onaxis) ** 2)
         else:
             assert type(omega0) == float
 
         # check the complex field convention and correct if needed
         if omega0 < 0:
             omega0 *= -1
-            data = np.conj(data)
+            data = xp.conj(data)
             print("Warning: input field will be conjugated")
 
         # remove the envelope frequency
-        data *= np.exp(1j * omega0 * t[None, None, :])
+        data *= xp.exp(1j * omega0 * t[None, None, :])
 
         # created LASY profile using FromArrayProfile class
-        wavelength = 2 * np.pi * c / omega0
+        wavelength = 2 * xp.pi * c / omega0
         dim = "xyt"
         axes = {"x": x, "y": y, "t": t}
         super().__init__(
diff --git a/lasy/profiles/from_openpmd_profile.py b/lasy/profiles/from_openpmd_profile.py
index a44b39b3..2aaa74e5 100644
--- a/lasy/profiles/from_openpmd_profile.py
+++ b/lasy/profiles/from_openpmd_profile.py
@@ -1,8 +1,8 @@
-import numpy as np
 import openpmd_api as io
 from openpmd_viewer import OpenPMDTimeSeries
 from scipy.constants import c
 
+from lasy.backend import xp
 from lasy.utils.laser_utils import create_grid, field_to_envelope
 from lasy.utils.openpmd_input import reorder_array
 
@@ -85,7 +85,7 @@ def __init__(
 
         # If `is_envelope` is not given, assume that complex arrays are envelopes.
         if is_envelope is None:
-            is_envelope = np.iscomplexobj(F)
+            is_envelope = xp.iscomplexobj(F)
 
         if theta is None:  # Envelope obtained from the full 3D array
             dim = "xyt"
@@ -109,7 +109,7 @@ def __init__(
             omg0 = it.meshes["laserEnvelope"].get_attribute("angularFrequency")
             array = F
 
-        wavelength = 2 * np.pi * c / omg0
+        wavelength = 2 * xp.pi * c / omg0
         if verbose:
             print(
                 "Wavelength used in the definition of the envelope (nm):",
diff --git a/lasy/profiles/gaussian_profile.py b/lasy/profiles/gaussian_profile.py
index 81e56a04..1d0df69e 100644
--- a/lasy/profiles/gaussian_profile.py
+++ b/lasy/profiles/gaussian_profile.py
@@ -89,7 +89,7 @@ class GaussianProfile(CombinedLongitudinalTransverseProfile):
     >>> extent[2:] *= 1e6
     >>> extent[:2] *= 1e15
     >>> tmin, tmax, rmin, rmax = extent
-    >>> vmax = np.abs(E_rt).max()
+    >>> vmax = xp.abs(E_rt).max()
     >>> plt.imshow(
     ...     E_rt,
     ...     origin="lower",
diff --git a/lasy/profiles/longitudinal/cosine_profile.py b/lasy/profiles/longitudinal/cosine_profile.py
index 108ddff7..d2fb195d 100644
--- a/lasy/profiles/longitudinal/cosine_profile.py
+++ b/lasy/profiles/longitudinal/cosine_profile.py
@@ -1,4 +1,4 @@
-import numpy as np
+from lasy.backend import xp
 
 from .longitudinal_profile import LongitudinalProfile
 
@@ -62,10 +62,10 @@ def evaluate(self, t):
         tn = (t - self.t_peak) / self.tau_fwhm
 
         envelope = (
-            np.cos(0.5 * np.pi * tn)
+            xp.cos(0.5 * xp.pi * tn)
             * (tn > -1)
             * (tn < 1)
-            * np.exp(+1.0j * (self.cep_phase + self.omega0 * self.t_peak))
+            * xp.exp(+1.0j * (self.cep_phase + self.omega0 * self.t_peak))
         )
 
         return envelope
diff --git a/lasy/profiles/longitudinal/gaussian_profile.py b/lasy/profiles/longitudinal/gaussian_profile.py
index 3d52131b..a05d4f90 100644
--- a/lasy/profiles/longitudinal/gaussian_profile.py
+++ b/lasy/profiles/longitudinal/gaussian_profile.py
@@ -1,4 +1,4 @@
-import numpy as np
+from lasy.backend import xp
 
 from .longitudinal_profile import LongitudinalProfile
 
@@ -55,7 +55,7 @@ def evaluate(self, t):
             Contains the value of the longitudinal envelope at the
             specified points. This array has the same shape as the array t.
         """
-        envelope = np.exp(
+        envelope = xp.exp(
             -((t - self.t_peak) ** 2) / self.tau**2
             + 1.0j * (self.cep_phase + self.omega0 * self.t_peak)
         )
diff --git a/lasy/profiles/longitudinal/longitudinal_profile.py b/lasy/profiles/longitudinal/longitudinal_profile.py
index 68a49c5f..19e64aaf 100644
--- a/lasy/profiles/longitudinal/longitudinal_profile.py
+++ b/lasy/profiles/longitudinal/longitudinal_profile.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.constants import c, pi
 
+from lasy.backend import xp
+
 
 class LongitudinalProfile(object):
     """
@@ -31,4 +32,4 @@ def evaluate(self, t):
         """
         # The base class only defines dummy fields
         # (This should be replaced by any class that inherits from this one.)
-        return np.zeros(t.shape, dtype="complex128")
+        return xp.zeros(t.shape, dtype="complex128")
diff --git a/lasy/profiles/longitudinal/longitudinal_profile_from_data.py b/lasy/profiles/longitudinal/longitudinal_profile_from_data.py
index 7047f313..3f89a91f 100644
--- a/lasy/profiles/longitudinal/longitudinal_profile_from_data.py
+++ b/lasy/profiles/longitudinal/longitudinal_profile_from_data.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.constants import c
 
+from lasy.backend import xp
+
 from .longitudinal_profile import LongitudinalProfile
 
 
@@ -62,58 +63,58 @@ def __init__(self, data, lo, hi):
         if data["datatype"] == "spectral":
             # First find central frequency
             wavelength = data["axis"]
-            assert np.all(
-                np.diff(wavelength) > 0
+            assert xp.all(
+                xp.diff(wavelength) > 0
             ), 'data["axis"] must be in monotonously increasing order.'
             spectral_intensity = data["intensity"]
             if data.get("phase") is None:
-                spectral_phase = np.zeros_like(wavelength)
+                spectral_phase = xp.zeros_like(wavelength)
             else:
                 spectral_phase = data["phase"]
             dt = data["dt"]
-            cwl = np.sum(spectral_intensity * wavelength) / np.sum(spectral_intensity)
+            cwl = xp.sum(spectral_intensity * wavelength) / xp.sum(spectral_intensity)
             cfreq = c / cwl
             # Determine required sampling frequency for desired dt
             sample_freq = 1 / dt
             # Determine number of points in temporal domain. This is the number of
             # points required to maintain the input spectral resolution while spanning
             # enough spectrum to achieve the desired temporal resolution.
-            indx = np.argmin(np.abs(wavelength - cwl))
-            dfreq = np.abs(c / wavelength[indx] - c / wavelength[indx + 1])
+            indx = xp.argmin(xp.abs(wavelength - cwl))
+            dfreq = xp.abs(c / wavelength[indx] - c / wavelength[indx + 1])
             N = int(sample_freq / dfreq)
-            freq = np.linspace(cfreq - sample_freq / 2, cfreq + sample_freq / 2, N)
+            freq = xp.linspace(cfreq - sample_freq / 2, cfreq + sample_freq / 2, N)
             # interpolate the spectrum onto this new array
-            freq_intensity = np.interp(
+            freq_intensity = xp.interp(
                 freq, c / wavelength[::-1], spectral_intensity[::-1], left=0, right=0
             )
-            freq_phase = np.interp(
+            freq_phase = xp.interp(
                 freq, c / wavelength[::-1], spectral_phase[::-1], left=0, right=0
             )
 
-            freq_amplitude = np.sqrt(freq_intensity)
+            freq_amplitude = xp.sqrt(freq_intensity)
 
             # Inverse Fourier Transform to the time domain
             t_amplitude = (
-                np.fft.fftshift(
-                    np.fft.ifft(
-                        np.fft.ifftshift(freq_amplitude * np.exp(-1j * freq_phase))
+                xp.fft.fftshift(
+                    xp.fft.ifft(
+                        xp.fft.ifftshift(freq_amplitude * xp.exp(-1j * freq_phase))
                     )
                 )
                 / dt
             )
-            time = np.linspace(-dt * N / 2, dt * N / 2 - dt, N)
+            time = xp.linspace(-dt * N / 2, dt * N / 2 - dt, N)
 
             # Extract intensity and phase
-            temporal_intensity = np.abs(t_amplitude) ** 2
-            temporal_intensity /= np.max(temporal_intensity)
-            temporal_phase = np.unwrap(-np.angle(t_amplitude))
-            temporal_phase -= temporal_phase[np.argmin(np.abs(time))]
+            temporal_intensity = xp.abs(t_amplitude) ** 2
+            temporal_intensity /= xp.max(temporal_intensity)
+            temporal_phase = xp.unwrap(-xp.angle(t_amplitude))
+            temporal_phase -= temporal_phase[xp.argmin(xp.abs(time))]
 
         elif data["datatype"] == "temporal":
             time = data["axis"]
             temporal_intensity = data["intensity"]
             if data.get("phase") is None:
-                temporal_phase = np.zeros_like(time)
+                temporal_phase = xp.zeros_like(time)
             else:
                 temporal_phase = data["phase"]
             cwl = data["wavelength"]
@@ -126,8 +127,8 @@ def __init__(self, data, lo, hi):
         # Finally crop the temporal domain to the physical domain
         # of interest
 
-        tIndLo = np.argmin(np.abs(time - lo))
-        tIndHi = np.argmin(np.abs(time - hi))
+        tIndLo = xp.argmin(xp.abs(time - lo))
+        tIndHi = xp.argmin(xp.abs(time - hi))
 
         self.time = time[tIndLo:tIndHi]
         self.temporal_intensity = temporal_intensity[tIndLo:tIndHi]
@@ -148,9 +149,9 @@ def evaluate(self, t):
             Contains the value of the longitudinal envelope at the
             specified points. This array has the same shape as the array t.
         """
-        intensity = np.interp(t, self.time, self.temporal_intensity)
-        phase = np.interp(t, self.time, self.temporal_phase)
+        intensity = xp.interp(t, self.time, self.temporal_intensity)
+        phase = xp.interp(t, self.time, self.temporal_phase)
 
-        envelope = np.sqrt(intensity) * np.exp(-1j * phase)
+        envelope = xp.sqrt(intensity) * xp.exp(-1j * phase)
 
         return envelope
diff --git a/lasy/profiles/longitudinal/super_gaussian_profile.py b/lasy/profiles/longitudinal/super_gaussian_profile.py
index 8668d47d..a9f6026c 100644
--- a/lasy/profiles/longitudinal/super_gaussian_profile.py
+++ b/lasy/profiles/longitudinal/super_gaussian_profile.py
@@ -1,4 +1,4 @@
-import numpy as np
+from lasy.backend import xp
 
 from .longitudinal_profile import LongitudinalProfile
 
@@ -61,8 +61,8 @@ def evaluate(self, t):
             Contains the value of the longitudinal envelope at the
             specified points. This array has the same shape as the array t.
         """
-        envelope = np.exp(
-            -np.power(((t - self.t_peak) ** 2) / self.tau**2, self.n_order / 2)
+        envelope = xp.exp(
+            -xp.power(((t - self.t_peak) ** 2) / self.tau**2, self.n_order / 2)
             + 1.0j * (self.cep_phase + self.omega0 * self.t_peak)
         )
 
diff --git a/lasy/profiles/profile.py b/lasy/profiles/profile.py
index d64e8ddb..e25d6777 100644
--- a/lasy/profiles/profile.py
+++ b/lasy/profiles/profile.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.constants import c
 
+from lasy.backend import xp
+
 
 class Profile(object):
     """
@@ -29,10 +30,10 @@ class Profile(object):
 
     def __init__(self, wavelength, pol):
         assert len(pol) == 2
-        norm_pol = np.sqrt(np.abs(pol[0]) ** 2 + np.abs(pol[1]) ** 2)
-        self.pol = np.array([pol[0] / norm_pol, pol[1] / norm_pol])
+        norm_pol = xp.sqrt(xp.abs(pol[0]) ** 2 + xp.abs(pol[1]) ** 2)
+        self.pol = xp.array([pol[0] / norm_pol, pol[1] / norm_pol])
         self.lambda0 = wavelength
-        self.omega0 = 2 * np.pi * c / self.lambda0
+        self.omega0 = 2 * xp.pi * c / self.lambda0
 
     def evaluate(self, x, y, t):
         """
@@ -52,7 +53,7 @@ def evaluate(self, x, y, t):
         """
         # The base class only defines dummy fields
         # (This should be replaced by any class that inherits from this one.)
-        return np.zeros_like(x, dtype="complex128")
+        return xp.zeros_like(x, dtype="complex128")
 
     def __add__(self, other):
         """Return the sum of two profiles."""
@@ -90,12 +91,12 @@ def __init__(self, *profiles):
         lambda0s = [p.lambda0 for p in self.profiles]
         pols = [p.pol for p in self.profiles]
         # Check that all wavelengths are the same
-        assert np.allclose(
+        assert xp.allclose(
             lambda0s, lambda0s[0]
         ), "Added profiles must have the same wavelength."
         lambda0 = profiles[0].lambda0
         # Check that all polarizations are the same
-        assert np.allclose(
+        assert xp.allclose(
             pols, pols[0]
         ), "Added profiles must have the same polarization."
         pol = profiles[0].pol
diff --git a/lasy/profiles/speckle_profile.py b/lasy/profiles/speckle_profile.py
index 6dddc5eb..4d68f587 100644
--- a/lasy/profiles/speckle_profile.py
+++ b/lasy/profiles/speckle_profile.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.constants import c
 
+from lasy.backend import xp
+
 from .profile import Profile
 
 
@@ -24,18 +25,18 @@ def gen_gaussian_time_series(t_num, dt, fwhm, rms_mean):
         A time series array of complex numbers with shape [t_num]
     """
     if fwhm == 0.0:
-        temporal_amplitude = np.zeros(t_num, dtype=np.complex128)
+        temporal_amplitude = xp.zeros(t_num, dtype=xp.complex128)
     else:
-        omega = np.fft.fftshift(np.fft.fftfreq(t_num, d=dt))
-        psd = np.exp(-np.log(2) * 0.5 * np.square(omega / fwhm * 2 * np.pi))
-        spectral_amplitude = np.array(psd) * (
-            np.random.normal(size=t_num) + 1j * np.random.normal(size=t_num)
+        omega = xp.fft.fftshift(xp.fft.fftfreq(t_num, d=dt))
+        psd = xp.exp(-xp.log(2) * 0.5 * xp.square(omega / fwhm * 2 * xp.pi))
+        spectral_amplitude = xp.array(psd) * (
+            xp.random.normal(size=t_num) + 1j * xp.random.normal(size=t_num)
         )
-        temporal_amplitude = np.fft.ifftshift(
-            np.fft.fft(np.fft.fftshift(spectral_amplitude))
+        temporal_amplitude = xp.fft.ifftshift(
+            xp.fft.fft(xp.fft.fftshift(spectral_amplitude))
         )
-        temporal_amplitude *= rms_mean / np.sqrt(
-            np.mean(np.square(np.abs(temporal_amplitude)))
+        temporal_amplitude *= rms_mean / xp.sqrt(
+            xp.mean(xp.square(xp.abs(temporal_amplitude)))
         )
     return temporal_amplitude
 
@@ -189,8 +190,8 @@ def __init__(
         super().__init__(wavelength, pol)
         self.laser_energy = laser_energy
         self.focal_length = focal_length
-        self.beam_aperture = np.array(beam_aperture, dtype="float")
-        self.n_beamlets = np.array(n_beamlets, dtype="int")
+        self.beam_aperture = xp.array(beam_aperture, dtype="float")
+        self.n_beamlets = n_beamlets
         self.temporal_smoothing_type = temporal_smoothing_type
         self.laser_bandwidth = relative_laser_bandwidth
 
@@ -198,22 +199,22 @@ def __init__(
         self.dt_update = 1 / self.laser_bandwidth / 50
         self.do_include_transverse_envelope = do_include_transverse_envelope
 
-        self.x_lens_list = np.linspace(
+        self.x_lens_list = xp.linspace(
             -0.5 * (self.n_beamlets[0] - 1),
             0.5 * (self.n_beamlets[0] - 1),
             num=self.n_beamlets[0],
         )
-        self.y_lens_list = np.linspace(
+        self.y_lens_list = xp.linspace(
             -0.5 * (self.n_beamlets[1] - 1),
             0.5 * (self.n_beamlets[1] - 1),
             num=self.n_beamlets[1],
         )
-        self.Y_lens_matrix, self.X_lens_matrix = np.meshgrid(
+        self.Y_lens_matrix, self.X_lens_matrix = xp.meshgrid(
             self.y_lens_list, self.x_lens_list
         )
-        self.Y_lens_index_matrix, self.X_lens_index_matrix = np.meshgrid(
-            np.arange(self.n_beamlets[1], dtype=float),
-            np.arange(self.n_beamlets[0], dtype=float),
+        self.Y_lens_index_matrix, self.X_lens_index_matrix = xp.meshgrid(
+            xp.arange(self.n_beamlets[1], dtype=float),
+            xp.arange(self.n_beamlets[0], dtype=float),
         )
 
         if "SSD" in self.temporal_smoothing_type.upper():
@@ -226,7 +227,7 @@ def __init__(
             self.ssd_transverse_bandwidth_distribution = (
                 ssd_transverse_bandwidth_distribution
             )
-            ssd_normalization = np.sqrt(
+            ssd_normalization = xp.sqrt(
                 self.ssd_transverse_bandwidth_distribution[0] ** 2
                 + self.ssd_transverse_bandwidth_distribution[1] ** 2
             )
@@ -260,7 +261,7 @@ def __init__(
             SpeckleProfile.supported_smoothing
         )
         assert relative_laser_bandwidth > 0, "laser_bandwidth must be greater than 0"
-        assert np.size(n_beamlets) == 2, "has to be a size 2 array"
+        assert len(n_beamlets) == 2, "has to be a size 2 array"
         if "SSD" in self.temporal_smoothing_type.upper():
             assert (
                 ssd_number_color_cycles is not None
@@ -276,7 +277,7 @@ def __init__(
                 ssd_transverse_bandwidth_distribution,
                 ssd_phase_modulation_amplitude,
             ):
-                assert np.size(q) == 2, "has to be a size 2 array"
+                assert xp.size(q) == 2, "has to be a size 2 array"
                 assert q[0] > 0 or q[1] > 0, "cannot be all zeros"
 
     def init_gaussian_time_series(
@@ -315,36 +316,36 @@ def init_gaussian_time_series(
         if "SSD" in self.temporal_smoothing_type.upper():
             pm_phase0 = gen_gaussian_time_series(
                 series_time.size
-                + int(np.sum(self.ssd_time_delay) / self.dt_update)
+                + int(xp.sum(self.ssd_time_delay) / self.dt_update)
                 + 2,
                 self.dt_update,
-                2 * np.pi * self.ssd_phase_modulation_frequency[0],
+                2 * xp.pi * self.ssd_phase_modulation_frequency[0],
                 self.ssd_phase_modulation_amplitude[0],
             )
             pm_phase1 = gen_gaussian_time_series(
                 series_time.size
-                + int(np.sum(self.ssd_time_delay) / self.dt_update)
+                + int(xp.sum(self.ssd_time_delay) / self.dt_update)
                 + 2,
                 self.dt_update,
-                2 * np.pi * self.ssd_phase_modulation_frequency[1],
+                2 * xp.pi * self.ssd_phase_modulation_frequency[1],
                 self.ssd_phase_modulation_amplitude[1],
             )
-            time_interp = np.arange(
+            time_interp = xp.arange(
                 start=0,
-                stop=series_time[-1] + np.sum(self.ssd_time_delay) + 3 * self.dt_update,
+                stop=series_time[-1] + xp.sum(self.ssd_time_delay) + 3 * self.dt_update,
                 step=self.dt_update,
             )[: pm_phase0.size]
             return (
                 time_interp,
                 [
-                    (np.real(pm_phase0) + np.imag(pm_phase0)) / np.sqrt(2),
-                    (np.real(pm_phase1) + np.imag(pm_phase1)) / np.sqrt(2),
+                    (xp.real(pm_phase0) + xp.imag(pm_phase0)) / xp.sqrt(2),
+                    (xp.real(pm_phase1) + xp.imag(pm_phase1)) / xp.sqrt(2),
                 ],
             )
         elif "ISI" in self.temporal_smoothing_type.upper():
-            complex_amp = np.stack(
+            complex_amp = xp.stack(
                 [
-                    np.stack(
+                    xp.stack(
                         [
                             gen_gaussian_time_series(
                                 series_time.size,
@@ -393,37 +394,37 @@ def beamlets_complex_amplitude(
         if any(
             rpp_type in temporal_smoothing_type.upper() for rpp_type in ["RPP", "CPP"]
         ):
-            return np.ones_like(self.X_lens_matrix)
+            return xp.ones_like(self.X_lens_matrix)
         if temporal_smoothing_type.upper() == "FM SSD":
-            phase_t = self.ssd_phase_modulation_amplitude[0] * np.sin(
+            phase_t = self.ssd_phase_modulation_amplitude[0] * xp.sin(
                 self.ssd_x_y_dephasing[0]
                 + 2
-                * np.pi
+                * xp.pi
                 * self.ssd_phase_modulation_frequency[0]
                 * (
                     t_now
                     - self.X_lens_matrix * self.ssd_time_delay[0] / self.n_beamlets[0]
                 )
-            ) + self.ssd_phase_modulation_amplitude[1] * np.sin(
+            ) + self.ssd_phase_modulation_amplitude[1] * xp.sin(
                 self.ssd_x_y_dephasing[1]
                 + 2
-                * np.pi
+                * xp.pi
                 * self.ssd_phase_modulation_frequency[1]
                 * (
                     t_now
                     - self.Y_lens_matrix * self.ssd_time_delay[1] / self.n_beamlets[1]
                 )
             )
-            return np.exp(1j * phase_t)
+            return xp.exp(1j * phase_t)
         elif temporal_smoothing_type.upper() == "GP RPM SSD":
-            phase_t = np.interp(
+            phase_t = xp.interp(
                 t_now
                 + self.X_lens_index_matrix
                 * self.ssd_time_delay[0]
                 / self.n_beamlets[0],
                 series_time,
                 time_series[0],
-            ) + np.interp(
+            ) + xp.interp(
                 t_now
                 + self.Y_lens_index_matrix
                 * self.ssd_time_delay[1]
@@ -431,7 +432,7 @@ def beamlets_complex_amplitude(
                 series_time,
                 time_series[1],
             )
-            return np.exp(1j * phase_t)
+            return xp.exp(1j * phase_t)
         elif temporal_smoothing_type.upper() == "GP ISI":
             return time_series[:, :, int(round(t_now / self.dt_update))]
         else:
@@ -475,21 +476,21 @@ def generate_speckle_pattern(
         Y_focus_matrix = y[:, :, 0] / lambda_fnum[1]
         x_focus_list = X_focus_matrix[:, 0]
         y_focus_list = Y_focus_matrix[0, :]
-        x_phase_focus_matrix = np.exp(
+        x_phase_focus_matrix = xp.exp(
             -2
-            * np.pi
+            * xp.pi
             * 1j
             / self.n_beamlets[0]
-            * self.x_lens_list[:, np.newaxis]
-            * x_focus_list[np.newaxis, :]
+            * self.x_lens_list[:, xp.newaxis]
+            * x_focus_list[xp.newaxis, :]
         )
-        y_phase_focus_matrix = np.exp(
+        y_phase_focus_matrix = xp.exp(
             -2
-            * np.pi
+            * xp.pi
             * 1j
             / self.n_beamlets[1]
-            * self.y_lens_list[:, np.newaxis]
-            * y_focus_list[np.newaxis, :]
+            * self.y_lens_list[:, xp.newaxis]
+            * y_focus_list[xp.newaxis, :]
         )
 
         bca = self.beamlets_complex_amplitude(
@@ -498,15 +499,15 @@ def generate_speckle_pattern(
             time_series=time_series,
             temporal_smoothing_type=self.temporal_smoothing_type,
         )
-        speckle_amp = np.einsum(
+        speckle_amp = xp.einsum(
             "jk,jl->kl",
-            np.einsum("ij,ik->jk", bca * exp_phase_plate, x_phase_focus_matrix),
+            xp.einsum("ij,ik->jk", bca * exp_phase_plate, x_phase_focus_matrix),
             y_phase_focus_matrix,
         )
         if self.do_include_transverse_envelope:
             speckle_amp = (
-                np.sinc(X_focus_matrix / self.n_beamlets[0])
-                * np.sinc(Y_focus_matrix / self.n_beamlets[1])
+                xp.sinc(X_focus_matrix / self.n_beamlets[0])
+                * xp.sinc(Y_focus_matrix / self.n_beamlets[1])
                 * speckle_amp
             )
         return speckle_amp
@@ -533,30 +534,30 @@ def evaluate(self, x, y, t):
 
         # Calculate auxiliary parameters
         if "RPP" == self.temporal_smoothing_type.upper():
-            phase_plate = np.random.choice([0, np.pi], self.n_beamlets)
+            phase_plate = xp.random.choice([0, xp.pi], self.n_beamlets)
         elif any(
             cpp_smoothing_type in self.temporal_smoothing_type.upper()
             for cpp_smoothing_type in ["CPP", "SSD"]
         ):
-            phase_plate = np.random.uniform(
-                -np.pi, np.pi, size=self.n_beamlets[0] * self.n_beamlets[1]
+            phase_plate = xp.random.uniform(
+                -xp.pi, xp.pi, size=self.n_beamlets[0] * self.n_beamlets[1]
             ).reshape(self.n_beamlets)
         elif "ISI" in self.temporal_smoothing_type.upper():
-            phase_plate = np.zeros(self.n_beamlets)  # ISI does not require phase plates
+            phase_plate = xp.zeros(self.n_beamlets)  # ISI does not require phase plates
         else:
             raise NotImplementedError
-        exp_phase_plate = np.exp(1j * phase_plate)
+        exp_phase_plate = xp.exp(1j * phase_plate)
         if self.temporal_smoothing_type.upper() == "FM SSD":
-            self.ssd_x_y_dephasing = np.random.standard_normal(2) * np.pi
+            self.ssd_x_y_dephasing = xp.random.standard_normal(2) * xp.pi
 
-        series_time = np.arange(0, t_max + self.dt_update, self.dt_update)
+        series_time = xp.arange(0, t_max + self.dt_update, self.dt_update)
 
         if "GP" in self.temporal_smoothing_type.upper():
             new_series_time, time_series = self.init_gaussian_time_series(series_time)
         else:
             new_series_time, time_series = series_time, None
 
-        envelope = np.zeros(x.shape, dtype=complex)
+        envelope = xp.zeros(x.shape, dtype=complex)
         for i, t_i in enumerate(t_norm):
             envelope[:, :, i] = self.generate_speckle_pattern(
                 t_i,
diff --git a/lasy/profiles/transverse/gaussian_profile.py b/lasy/profiles/transverse/gaussian_profile.py
index 9905f2ca..4b3ddc88 100644
--- a/lasy/profiles/transverse/gaussian_profile.py
+++ b/lasy/profiles/transverse/gaussian_profile.py
@@ -1,4 +1,4 @@
-import numpy as np
+from lasy.backend import xp
 
 from .transverse_profile import TransverseProfile
 
@@ -50,7 +50,7 @@ def __init__(self, w0, wavelength=None, z_foc=0):
             assert (
                 wavelength is not None
             ), "You need to pass the wavelength, when `z_foc` is non-zero."
-            self.z_foc_over_zr = z_foc * wavelength / (np.pi * w0**2)
+            self.z_foc_over_zr = z_foc * wavelength / (xp.pi * w0**2)
 
     def _evaluate(self, x, y):
         """
@@ -73,6 +73,6 @@ def _evaluate(self, x, y):
         # Calculate the argument of the complex exponential
         exp_argument = -(x**2 + y**2) / (self.w0**2 * diffract_factor)
         # Get the transverse profile
-        envelope = np.exp(exp_argument) / diffract_factor
+        envelope = xp.exp(exp_argument) / diffract_factor
 
         return envelope
diff --git a/lasy/profiles/transverse/hermite_gaussian_profile.py b/lasy/profiles/transverse/hermite_gaussian_profile.py
index 8511361c..67ffa504 100644
--- a/lasy/profiles/transverse/hermite_gaussian_profile.py
+++ b/lasy/profiles/transverse/hermite_gaussian_profile.py
@@ -1,8 +1,9 @@
 from math import factorial
 
-import numpy as np
 from scipy.special import hermite
 
+from lasy.backend import xp
+
 from .transverse_profile import TransverseProfile
 
 
@@ -63,7 +64,7 @@ def __init__(self, w0, n_x, n_y, wavelength=None, z_foc=0):
             assert (
                 wavelength is not None
             ), "You need to pass the wavelength, when `z_foc` is non-zero."
-            self.z_foc_over_zr = z_foc * wavelength / (np.pi * w0**2)
+            self.z_foc_over_zr = z_foc * wavelength / (xp.pi * w0**2)
 
     def _evaluate(self, x, y):
         """
@@ -84,15 +85,15 @@ def _evaluate(self, x, y):
         # Term for wavefront curvature, waist and Gouy phase
         diffract_factor = 1.0 - 1j * self.z_foc_over_zr
         w = self.w0 * abs(diffract_factor)
-        psi = np.angle(diffract_factor)
+        psi = xp.angle(diffract_factor)
 
         envelope = (
-            np.sqrt(2 / np.pi)
-            * np.sqrt(1 / (2 ** (self.n_x) * factorial(self.n_x) * self.w0))
-            * np.sqrt(1 / (2 ** (self.n_y) * factorial(self.n_y) * self.w0))
-            * hermite(self.n_x)(np.sqrt(2) * x / w)
-            * hermite(self.n_y)(np.sqrt(2) * y / w)
-            * np.exp(
+            xp.sqrt(2 / xp.pi)
+            * xp.sqrt(1 / (2 ** (self.n_x) * factorial(self.n_x) * self.w0))
+            * xp.sqrt(1 / (2 ** (self.n_y) * factorial(self.n_y) * self.w0))
+            * hermite(self.n_x)(xp.sqrt(2) * x / w)
+            * hermite(self.n_y)(xp.sqrt(2) * y / w)
+            * xp.exp(
                 -(x**2 + y**2) / (self.w0**2 * diffract_factor)
                 - 1.0j * (self.n_x + self.n_y) * psi
             )
diff --git a/lasy/profiles/transverse/jinc_profile.py b/lasy/profiles/transverse/jinc_profile.py
index 3eb92eed..aea79177 100644
--- a/lasy/profiles/transverse/jinc_profile.py
+++ b/lasy/profiles/transverse/jinc_profile.py
@@ -1,6 +1,7 @@
-import numpy as np
 import scipy.special as scispe
 
+from lasy.backend import xp
+
 from .transverse_profile import TransverseProfile
 
 
@@ -42,11 +43,11 @@ def _evaluate(self, x, y):
             Contains the value of the envelope at the specified points
             This array has the same shape as the arrays x, y
         """
-        r_over_w0 = np.sqrt(x**2 + y**2) / self.w0
+        r_over_w0 = xp.sqrt(x**2 + y**2) / self.w0
 
-        envelope = np.ones_like(r_over_w0)
+        envelope = xp.ones_like(r_over_w0)
         # Avoid dividing by zero
-        np.divide(
+        xp.divide(
             2.0 * scispe.jv(1, r_over_w0),
             r_over_w0,
             out=envelope,
diff --git a/lasy/profiles/transverse/laguerre_gaussian_profile.py b/lasy/profiles/transverse/laguerre_gaussian_profile.py
index 9df47223..450b1c55 100644
--- a/lasy/profiles/transverse/laguerre_gaussian_profile.py
+++ b/lasy/profiles/transverse/laguerre_gaussian_profile.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.special import genlaguerre
 
+from lasy.backend import xp
+
 from .transverse_profile import TransverseProfile
 
 
@@ -61,7 +62,7 @@ def __init__(self, w0, p, m, wavelength=None, z_foc=0):
             assert (
                 wavelength is not None
             ), "You need to pass the wavelength, when `z_foc` is non-zero."
-            self.z_foc_over_zr = z_foc * wavelength / (np.pi * w0**2)
+            self.z_foc_over_zr = z_foc * wavelength / (xp.pi * w0**2)
 
     def _evaluate(self, x, y):
         """
@@ -82,7 +83,7 @@ def _evaluate(self, x, y):
         # Term for wavefront curvature, waist and Gouy phase
         diffract_factor = 1.0 - 1j * self.z_foc_over_zr
         w = self.w0 * abs(diffract_factor)
-        psi = np.angle(diffract_factor)
+        psi = xp.angle(diffract_factor)
         # complex_position corresponds to r e^{+/-i\theta}
         if self.m > 0:
             complex_position = x - 1j * y
@@ -92,7 +93,7 @@ def _evaluate(self, x, y):
         envelope = (
             complex_position ** abs(self.m)
             * genlaguerre(self.p, abs(self.m))(2 * radius**2 / w**2)
-            * np.exp(
+            * xp.exp(
                 -(radius**2) / (self.w0**2 * diffract_factor)
                 - 1.0j * (2 * self.p + self.m) * psi
             )  # Additional Gouy phase
diff --git a/lasy/profiles/transverse/super_gaussian_profile.py b/lasy/profiles/transverse/super_gaussian_profile.py
index 395d565c..17e2629a 100644
--- a/lasy/profiles/transverse/super_gaussian_profile.py
+++ b/lasy/profiles/transverse/super_gaussian_profile.py
@@ -1,4 +1,4 @@
-import numpy as np
+from lasy.backend import xp
 
 from .transverse_profile import TransverseProfile
 
@@ -45,6 +45,6 @@ def _evaluate(self, x, y):
             Contains the value of the envelope at the specified points
             This array has the same shape as the arrays x, y
         """
-        envelope = np.exp(-np.power((x**2 + y**2) / self.w0**2, self.n_order / 2))
+        envelope = xp.exp(-xp.power((x**2 + y**2) / self.w0**2, self.n_order / 2))
 
         return envelope
diff --git a/lasy/profiles/transverse/transverse_profile.py b/lasy/profiles/transverse/transverse_profile.py
index 77a30704..a771deb2 100644
--- a/lasy/profiles/transverse/transverse_profile.py
+++ b/lasy/profiles/transverse/transverse_profile.py
@@ -1,4 +1,4 @@
-import numpy as np
+from lasy.backend import xp
 
 
 class TransverseProfile(object):
@@ -32,7 +32,7 @@ def _evaluate(self, x, y):
         """
         # The base class only defines dummy fields
         # (This should be replaced by any class that inherits from this one.)
-        return np.zeros(x.shape, dtype="complex128")
+        return xp.zeros(x.shape, dtype="complex128")
 
     def __add__(self, other):
         """Return the sum of two transverse profiles."""
diff --git a/lasy/profiles/transverse/transverse_profile_from_data.py b/lasy/profiles/transverse/transverse_profile_from_data.py
index 68bbe0b4..ebdca3a6 100644
--- a/lasy/profiles/transverse/transverse_profile_from_data.py
+++ b/lasy/profiles/transverse/transverse_profile_from_data.py
@@ -1,6 +1,6 @@
-import numpy as np
 from scipy.interpolate import RegularGridInterpolator
 
+from lasy.backend import xp
 from lasy.utils.exp_data_utils import find_center_of_mass
 
 from .transverse_profile import TransverseProfile
@@ -48,35 +48,35 @@ def __init__(self, intensity_data, lo, hi, center_data=True):
 
         intensity_data = intensity_data.astype("float64")
 
-        n_y, n_x = np.shape(intensity_data)
+        n_y, n_x = xp.shape(intensity_data)
 
         dx = (hi[0] - lo[0]) / n_x
         dy = (hi[1] - lo[1]) / n_y
 
         if center_data:
-            x_range = np.abs(hi[0] - lo[0])
-            y_range = np.abs(hi[1] - lo[1])
-            x_data = np.linspace(-x_range / 2, x_range / 2, n_x)
-            y_data = np.linspace(-y_range / 2, y_range / 2, n_y)
+            x_range = xp.abs(hi[0] - lo[0])
+            y_range = xp.abs(hi[1] - lo[1])
+            x_data = xp.linspace(-x_range / 2, x_range / 2, n_x)
+            y_data = xp.linspace(-y_range / 2, y_range / 2, n_y)
 
             n_x0, n_y0 = find_center_of_mass(intensity_data)
-            intensity_data = np.roll(
-                np.roll(intensity_data, -int(n_x0 - n_x / 2), axis=1),
+            intensity_data = xp.roll(
+                xp.roll(intensity_data, -int(n_x0 - n_x / 2), axis=1),
                 -int(n_y0 - n_y / 2),
                 axis=0,
             )
 
         else:
-            x_data = np.linspace(lo[0], hi[0], n_x)
-            y_data = np.linspace(lo[1], hi[1], n_y)
+            x_data = xp.linspace(lo[0], hi[0], n_x)
+            y_data = xp.linspace(lo[1], hi[1], n_y)
 
         # Normalise the profile such that its squared integeral == 1
-        intensity_data /= np.sum(intensity_data) * dx * dy
+        intensity_data /= xp.sum(intensity_data) * dx * dy
 
         # Note here we use the square root of intensity to get the 'field'
         self.field_interp = RegularGridInterpolator(
             (y_data, x_data),
-            np.sqrt(intensity_data),
+            xp.sqrt(intensity_data),
             bounds_error=False,
             fill_value=0.0,
         )
diff --git a/lasy/utils/exp_data_utils.py b/lasy/utils/exp_data_utils.py
index ef861c40..6e69365e 100644
--- a/lasy/utils/exp_data_utils.py
+++ b/lasy/utils/exp_data_utils.py
@@ -1,4 +1,4 @@
-import numpy as np
+from lasy.backend import xp
 
 
 def find_center_of_mass(img):
@@ -17,14 +17,14 @@ def find_center_of_mass(img):
         and the vertical. The units are in pixels.
 
     """
-    rows, cols = np.shape(img)
-    x = np.linspace(0, cols - 1, cols)
-    y = np.linspace(0, rows - 1, rows)
+    rows, cols = xp.shape(img)
+    x = xp.linspace(0, cols - 1, cols)
+    y = xp.linspace(0, rows - 1, rows)
 
     # find the beam center using COM
-    img_tot = np.sum(img)
-    x0 = np.sum(np.dot(img, x)) / img_tot
-    y0 = np.sum(np.dot(img.T, y)) / img_tot
+    img_tot = xp.sum(img)
+    x0 = xp.sum(xp.dot(img, x)) / img_tot
+    y0 = xp.sum(xp.dot(img.T, y)) / img_tot
 
     return x0, y0
 
@@ -44,14 +44,14 @@ def find_d4sigma(img):
     D4sigX : The D4sigma along the first (x) axis
     D4sigY : The D4sigma along the second (y) axis
     """
-    rows, cols = np.shape(img)
-    x = np.linspace(0, cols - 1, cols)
-    y = np.linspace(0, rows - 1, rows)
+    rows, cols = xp.shape(img)
+    x = xp.linspace(0, cols - 1, cols)
+    y = xp.linspace(0, rows - 1, rows)
 
     x0, y0 = find_center_of_mass(img)
 
-    img_tot = np.sum(img)
-    D4sigX = 4 * np.sqrt(np.sum(np.dot(img, (x - x0) ** 2)) / img_tot)
-    D4sigY = 4 * np.sqrt(np.sum(np.dot(img.T, (y - y0) ** 2)) / img_tot)
+    img_tot = xp.sum(img)
+    D4sigX = 4 * xp.sqrt(xp.sum(xp.dot(img, (x - x0) ** 2)) / img_tot)
+    D4sigY = 4 * xp.sqrt(xp.sum(xp.dot(img.T, (y - y0) ** 2)) / img_tot)
 
     return D4sigX, D4sigY
diff --git a/lasy/utils/grid.py b/lasy/utils/grid.py
index 54a3453d..9629314f 100644
--- a/lasy/utils/grid.py
+++ b/lasy/utils/grid.py
@@ -1,5 +1,7 @@
 import numpy as np
 
+from lasy.backend import use_cupy, xp
+
 time_axis_indx = -1
 
 
@@ -44,7 +46,7 @@ def __init__(self, dim, lo, hi, npoints, n_azimuthal_modes=None):
         self.axes = []
         self.dx = []
         for i in range(ndims):
-            self.axes.append(np.linspace(lo[i], hi[i], npoints[i]))
+            self.axes.append(xp.linspace(lo[i], hi[i], npoints[i]))
             self.dx.append(self.axes[i][1] - self.axes[i][0])
 
         if dim == "rt":
@@ -61,10 +63,9 @@ def __init__(self, dim, lo, hi, npoints, n_azimuthal_modes=None):
             # 0, 1, 2, ..., n_azimuthal_modes-1, -n_azimuthal_modes+1, ..., -1
             ncomp = 2 * self.n_azimuthal_modes - 1
             self.shape = (ncomp, self.npoints[0], self.npoints[1])
-
-        self.temporal_field = np.zeros(self.shape, dtype="complex128")
+        self.temporal_field = xp.zeros(self.shape, dtype="complex128")
         self.temporal_field_valid = False
-        self.spectral_field = np.zeros(self.shape, dtype="complex128")
+        self.spectral_field = xp.zeros(self.shape, dtype="complex128")
         self.spectral_field_valid = False
 
     def set_temporal_field(self, field):
@@ -78,6 +79,8 @@ def set_temporal_field(self, field):
         """
         assert field.shape == self.temporal_field.shape
         assert field.dtype == "complex128"
+        if use_cupy and type(field) == np.ndarray:
+            field = xp.asarray(field)  # Copy to GPU
         self.temporal_field[:, :, :] = field
         self.temporal_field_valid = True
         self.spectral_field_valid = False  # Invalidates the spectral field
@@ -93,11 +96,13 @@ def set_spectral_field(self, field):
         """
         assert field.shape == self.spectral_field.shape
         assert field.dtype == "complex128"
+        if use_cupy and type(field) == np.ndarray:
+            field = xp.asarray(field)  # Copy to GPU
         self.spectral_field[:, :, :] = field
         self.spectral_field_valid = True
         self.temporal_field_valid = False  # Invalidates the temporal field
 
-    def get_temporal_field(self):
+    def get_temporal_field(self, to_cpu=False):
         """
         Return a copy of the temporal field.
 
@@ -109,37 +114,41 @@ def get_temporal_field(self):
         field : ndarray of complexs
             The temporal field.
         """
-        # We return a copy, so that the user cannot modify
-        # the original field, unless get_temporal_field is called
-        if self.temporal_field_valid:
-            return self.temporal_field.copy()
-        elif self.spectral_field_valid:
+        if not self.temporal_field_valid:
             self.spectral2temporal_fft()
-            return self.temporal_field.copy()
+        # Return a copy of the field, either on CPU or GPU, so that the user
+        # cannot modify the original field, unless set_spectral_field is called
+        if to_cpu and use_cupy:
+            return xp.asnumpy(self.temporal_field)
         else:
-            raise ValueError("Both temporal and spectral fields are invalid")
+            return self.temporal_field.copy()
 
-    def get_spectral_field(self):
+    def get_spectral_field(self, to_cpu=False):
         """
         Return a copy of the spectral field.
 
         (Modifying the returned object will not modify the original field stored
         in the Grid object ; one must use set_spectral_field to do so.)
 
+        Parameters
+        ----------
+        to_cpu : bool
+            If True, the returned field is always returned as a numpy array on CPU
+            (even when the lasy backend is cupy)
+
         Returns
         -------
         field : ndarray of complexs
             The spectral field.
         """
-        # We return a copy, so that the user cannot modify
-        # the original field, unless set_spectral_field is called
-        if self.spectral_field_valid:
-            return self.spectral_field.copy()
-        elif self.temporal_field_valid:
+        if not self.spectral_field_valid:
             self.temporal2spectral_fft()
-            return self.spectral_field.copy()
+        # Return a copy of the field, either on CPU or GPU, so that the user
+        # cannot modify the original field, unless set_spectral_field is called
+        if to_cpu and use_cupy:
+            return xp.asnumpy(self.spectral_field)
         else:
-            raise ValueError("Both temporal and spectral fields are invalid")
+            return self.spectral_field.copy()
 
     def temporal2spectral_fft(self):
         """
@@ -148,7 +157,7 @@ def temporal2spectral_fft(self):
         (Only along the time axis, not along the transverse spatial coordinates.)
         """
         assert self.temporal_field_valid
-        self.spectral_field = np.fft.ifft(
+        self.spectral_field = xp.fft.ifft(
             self.temporal_field, axis=time_axis_indx, norm="backward"
         )
         self.spectral_field_valid = True
@@ -160,7 +169,7 @@ def spectral2temporal_fft(self):
         (Only along the time axis, not along the transverse spatial coordinates.)
         """
         assert self.spectral_field_valid
-        self.temporal_field = np.fft.fft(
+        self.temporal_field = xp.fft.fft(
             self.spectral_field, axis=time_axis_indx, norm="backward"
         )
         self.temporal_field_valid = True
diff --git a/lasy/utils/laser_utils.py b/lasy/utils/laser_utils.py
index 76e8fef3..c48ae04e 100644
--- a/lasy/utils/laser_utils.py
+++ b/lasy/utils/laser_utils.py
@@ -1,10 +1,11 @@
-import numpy as np
 from axiprop.containers import ScalarFieldEnvelope
 from axiprop.lib import PropagatorFFT2, PropagatorResampling
 from scipy.constants import c, e, epsilon_0, m_e
 from scipy.interpolate import interp1d
 from scipy.signal import hilbert
 
+from lasy.backend import use_cupy, xp
+
 from .grid import Grid
 
 
@@ -47,7 +48,7 @@ def compute_laser_energy(dim, grid):
         energy = ((dV * epsilon_0 * 0.5) * abs(envelope) ** 2).sum()
     else:  # dim == "rt":
         energy = (
-            dV[np.newaxis, :, np.newaxis]
+            dV[xp.newaxis, :, xp.newaxis]
             * epsilon_0
             * 0.5
             * abs(envelope[:, :, :]) ** 2
@@ -103,7 +104,7 @@ def normalize_peak_field_amplitude(amplitude, grid):
     """
     if amplitude is not None:
         field = grid.get_temporal_field()
-        field_max = np.abs(field).max()
+        field_max = xp.abs(field).max()
         if field_max == 0.0:
             print("Field is zero everywhere, normalization will be skipped")
         else:
@@ -125,12 +126,12 @@ def normalize_peak_intensity(peak_intensity, grid):
     """
     if peak_intensity is not None:
         field = grid.get_temporal_field()
-        intensity = np.abs(epsilon_0 * field**2 / 2 * c)
+        intensity = xp.abs(epsilon_0 * field**2 / 2 * c)
         input_peak_intensity = intensity.max()
         if input_peak_intensity == 0.0:
             print("Field is zero everywhere, normalization will be skipped")
         else:
-            grid.set_temporal_field(np.sqrt(peak_intensity / input_peak_intensity))
+            grid.set_temporal_field(xp.sqrt(peak_intensity / input_peak_intensity))
 
 
 def get_full_field(laser, theta=0, slice=0, slice_axis="x", Nt=None):
@@ -161,13 +162,13 @@ def get_full_field(laser, theta=0, slice=0, slice_axis="x", Nt=None):
     time_axis = laser.grid.axes[-1]
 
     if laser.dim == "rt":
-        azimuthal_phase = np.exp(-1j * laser.grid.azimuthal_modes * theta)
+        azimuthal_phase = xp.exp(-1j * laser.grid.azimuthal_modes * theta)
         env_upper = env * azimuthal_phase[:, None, None]
         env_upper = env_upper.sum(0)
-        azimuthal_phase = np.exp(1j * laser.grid.azimuthal_modes * theta)
+        azimuthal_phase = xp.exp(1j * laser.grid.azimuthal_modes * theta)
         env_lower = env * azimuthal_phase[:, None, None]
         env_lower = env_lower.sum(0)
-        env = np.vstack((env_lower[::-1][:-1], env_upper))
+        env = xp.vstack((env_lower[::-1][:-1], env_upper))
     elif slice_axis == "x":
         Nx_middle = env.shape[0] // 2 - 1
         Nx_slice = int((1 + slice) * Nx_middle)
@@ -181,24 +182,24 @@ def get_full_field(laser, theta=0, slice=0, slice_axis="x", Nt=None):
 
     if Nt is not None:
         Nr = env.shape[0]
-        time_axis_new = np.linspace(laser.grid.lo[-1], laser.grid.hi[-1], Nt)
-        env_new = np.zeros((Nr, Nt), dtype=env.dtype)
+        time_axis_new = xp.linspace(laser.grid.lo[-1], laser.grid.hi[-1], Nt)
+        env_new = xp.zeros((Nr, Nt), dtype=env.dtype)
 
         for ir in range(Nr):
-            interp_fu_abs = interp1d(time_axis, np.abs(env[ir]))
+            interp_fu_abs = interp1d(time_axis, xp.abs(env[ir]))
             slice_abs = interp_fu_abs(time_axis_new)
-            interp_fu_angl = interp1d(time_axis, np.unwrap(np.angle(env[ir])))
+            interp_fu_angl = interp1d(time_axis, xp.unwrap(xp.angle(env[ir])))
             slice_angl = interp_fu_angl(time_axis_new)
-            env_new[ir] = slice_abs * np.exp(1j * slice_angl)
+            env_new[ir] = slice_abs * xp.exp(1j * slice_angl)
 
         time_axis = time_axis_new
         env = env_new
 
-    env *= np.exp(-1j * omega0 * time_axis[None, :])
-    env = np.real(env)
+    env *= xp.exp(-1j * omega0 * time_axis[None, :])
+    env = xp.real(env)
 
     if laser.dim == "rt":
-        ext = np.array(
+        ext = xp.array(
             [
                 laser.grid.lo[-1],
                 laser.grid.hi[-1],
@@ -207,7 +208,7 @@ def get_full_field(laser, theta=0, slice=0, slice_axis="x", Nt=None):
             ]
         )
     else:
-        ext = np.array(
+        ext = xp.array(
             [
                 laser.grid.lo[-1],
                 laser.grid.hi[-1],
@@ -301,8 +302,8 @@ def get_spectrum(
         Array with the angular frequencies of the spectrum.
     """
     # Get the frequencies of the fft output.
-    freq = np.fft.fftfreq(grid.shape[-1], d=(grid.axes[-1][1] - grid.axes[-1][0]))
-    omega = 2 * np.pi * freq
+    freq = xp.fft.fftfreq(grid.shape[-1], d=(grid.axes[-1][1] - grid.axes[-1][0]))
+    omega = 2 * xp.pi * freq
 
     # Get on axis or full field.
     field = grid.get_temporal_field()
@@ -318,10 +319,10 @@ def get_spectrum(
         # Assume that the FFT of the envelope and the FFT of the complex
         # conjugate of the envelope do not overlap. Then we only need
         # one of them.
-        spectrum = 0.5 * np.fft.fft(field) * grid.dx[-1]
+        spectrum = 0.5 * xp.fft.fft(field) * grid.dx[-1]
         omega = omega0 - omega
         # Sort frequency array (and the spectrum accordingly).
-        i_sort = np.argsort(omega)
+        i_sort = xp.argsort(omega)
         omega = omega[i_sort]
         spectrum = spectrum[..., i_sort]
         # Keep only positive frequencies.
@@ -329,7 +330,7 @@ def get_spectrum(
         omega = omega[i_keep]
         spectrum = spectrum[..., i_keep]
     else:
-        spectrum = np.fft.fft(field) * grid.dx[-1]
+        spectrum = xp.fft.fft(field) * grid.dx[-1]
         # Keep only positive frequencies.
         i_keep = spectrum.shape[-1] // 2
         omega = omega[:i_keep]
@@ -337,21 +338,21 @@ def get_spectrum(
 
     # Convert to spectral energy density (J/(m^2 rad Hz)).
     if method != "raw":
-        spectrum = np.abs(spectrum) ** 2 * epsilon_0 * c / np.pi
+        spectrum = xp.abs(spectrum) ** 2 * epsilon_0 * c / xp.pi
 
     # Integrate transversely.
     if method == "sum":
         dV = get_grid_cell_volume(grid, dim)
         dz = grid.dx[-1] * c
         if dim == "xyt":
-            spectrum = np.sum(spectrum * dV / dz, axis=(0, 1))
+            spectrum = xp.sum(spectrum * dV / dz, axis=(0, 1))
         else:
-            spectrum = np.sum(spectrum[0] * dV[:, np.newaxis] / dz, axis=0)
+            spectrum = xp.sum(spectrum[0] * dV[:, xp.newaxis] / dz, axis=0)
 
     # If the user specified a frequency range, interpolate into it.
     if method in ["sum", "on_axis"] and range is not None:
-        omega_interp = np.linspace(*range, bins)
-        spectrum = np.interp(omega_interp, omega, spectrum)
+        omega_interp = xp.linspace(*range, bins)
+        spectrum = xp.interp(omega_interp, omega, spectrum)
         omega = omega_interp
 
     return spectrum, omega
@@ -431,16 +432,16 @@ def get_frequency(
     # Assumes t is last dimension!
     if is_envelope:
         assert omega0 is not None
-        phase = np.unwrap(np.angle(field))
-        omega = omega0 + np.gradient(-phase, grid.axes[-1], axis=-1, edge_order=2)
-        central_omega = np.average(omega, weights=np.abs(field))
+        phase = xp.unwrap(xp.angle(field))
+        omega = omega0 + xp.gradient(-phase, grid.axes[-1], axis=-1, edge_order=2)
+        central_omega = xp.average(omega, weights=xp.abs(field))
     else:
         assert dim in ["xyt", "rt"]
         if dim == "xyt" and phase_unwrap_nd:
             print("WARNING: using 3D phase unwrapping, this can be expensive")
 
         h = field if is_hilbert else hilbert_transform(grid)
-        h = np.squeeze(field)
+        h = xp.squeeze(field)
         if phase_unwrap_nd:
             try:
                 from skimage.restoration import unwrap_phase
@@ -452,22 +453,22 @@ def get_frequency(
                 "scikit-image must be install for nd phase unwrapping.",
                 "Please install scikit-image or use phase_unwrap_nd=False.",
             )
-            phase = unwrap_phase(np.angle(h))
+            phase = unwrap_phase(xp.angle(h))
         else:
-            phase = np.unwrap(np.angle(h))
-        omega = np.gradient(-phase, grid.axes[-1], axis=-1, edge_order=2)
+            phase = xp.unwrap(xp.angle(h))
+        omega = xp.gradient(-phase, grid.axes[-1], axis=-1, edge_order=2)
 
         if dim == "xyt":
-            weights = np.abs(h)
+            weights = xp.abs(h)
         else:
             r = grid.axes[0].reshape((grid.axes[0].size, 1))
-            weights = r * np.abs(h)
-        central_omega = np.average(omega, weights=weights)
+            weights = r * xp.abs(h)
+        central_omega = xp.average(omega, weights=weights)
 
     # Filter out too small frequencies
-    omega = np.maximum(omega, lower_bound * central_omega)
+    omega = xp.maximum(omega, lower_bound * central_omega)
     # Filter out too large frequencies
-    omega = np.minimum(omega, upper_bound * central_omega)
+    omega = xp.minimum(omega, upper_bound * central_omega)
 
     return omega, central_omega
 
@@ -491,11 +492,11 @@ def get_duration(grid, dim):
     dV = get_grid_cell_volume(grid, dim)
     field = grid.get_temporal_field()
     if dim == "xyt":
-        weights = np.abs(field) ** 2 * dV
+        weights = xp.abs(field) ** 2 * dV
     else:  # dim == "rt":
-        weights = np.abs(field) ** 2 * dV[np.newaxis, :, np.newaxis]
+        weights = xp.abs(field) ** 2 * dV[xp.newaxis, :, xp.newaxis]
     # project weights to longitudinal axes
-    weights = np.sum(weights, axis=(0, 1))
+    weights = xp.sum(weights, axis=(0, 1))
     return weighted_std(grid.axes[-1], weights)
 
 
@@ -549,7 +550,7 @@ def vector_potential_to_field(grid, omega0, direct=True):
     field = grid.get_temporal_field()
     if direct:
         A = (
-            -np.gradient(field, grid.axes[-1], axis=-1, edge_order=2)
+            -xp.gradient(field, grid.axes[-1], axis=-1, edge_order=2)
             + 1j * omega0 * field
         )
         return m_e * c / e * A
@@ -593,7 +594,7 @@ def field_to_envelope(grid, dim, phase_unwrap_nd=False):
         is_hilbert=True,
         phase_unwrap_nd=phase_unwrap_nd,
     )
-    field *= np.exp(1j * omg0_h * grid.axes[-1])
+    field *= xp.exp(1j * omg0_h * grid.axes[-1])
     grid.set_temporal_field(field)
 
     return grid, omg0_h
@@ -637,7 +638,7 @@ def get_grid_cell_volume(grid, dim):
         r = grid.axes[0]
         dr = grid.dx[0]
         # 1D array that computes the volume of radial cells
-        dV = np.pi * ((r + 0.5 * dr) ** 2 - (r - 0.5 * dr) ** 2) * dz
+        dV = xp.pi * ((r + 0.5 * dr) ** 2 - (r - 0.5 * dr) ** 2) * dz
     return dV
 
 
@@ -656,8 +657,8 @@ def weighted_std(values, weights=None):
     -------
     A float with the value of the standard deviation
     """
-    mean_val = np.average(values, weights=weights)
-    std = np.sqrt(np.average((values - mean_val) ** 2, weights=weights))
+    mean_val = xp.average(values, weights=weights)
+    std = xp.sqrt(xp.average((values - mean_val) ** 2, weights=weights))
     return std
 
 
@@ -684,17 +685,17 @@ def create_grid(array, axes, dim):
         hi = (axes["x"][-1], axes["y"][-1], axes["t"][-1])
         npoints = (axes["x"].size, axes["y"].size, axes["t"].size)
         grid = Grid(dim, lo, hi, npoints)
-        assert np.all(grid.axes[0] == axes["x"])
-        assert np.all(grid.axes[1] == axes["y"])
-        assert np.allclose(grid.axes[2], axes["t"], rtol=1.0e-14)
+        assert xp.all(grid.axes[0] == axes["x"])
+        assert xp.all(grid.axes[1] == axes["y"])
+        assert xp.allclose(grid.axes[2], axes["t"], rtol=1.0e-14)
         grid.set_temporal_field(array)
     else:  # dim == "rt":
         lo = (axes["r"][0], axes["t"][0])
         hi = (axes["r"][-1], axes["t"][-1])
         npoints = (axes["r"].size, axes["t"].size)
         grid = Grid(dim, lo, hi, npoints, n_azimuthal_modes=1)
-        assert np.all(grid.axes[0] == axes["r"])
-        assert np.allclose(grid.axes[1], axes["t"], rtol=1.0e-14)
+        assert xp.all(grid.axes[0] == axes["r"])
+        assert xp.allclose(grid.axes[1], axes["t"], rtol=1.0e-14)
         grid.set_temporal_field(array)
     return grid
 
@@ -737,6 +738,9 @@ def export_to_z(dim, grid, omega0, z_axis=None, z0=0.0, t0=0.0, backend="NP"):
     time_axis_indx = -1
 
     t_axis = grid.axes[time_axis_indx]
+    if use_cupy:
+        t_axis = xp.asnumpy(t_axis)
+
     if z_axis is None:
         z_axis = t_axis * c
 
@@ -758,18 +762,18 @@ def export_to_z(dim, grid, omega0, z_axis=None, z0=0.0, t0=0.0, backend="NP"):
                 )
             )
 
-        field_z = np.zeros(
+        field_z = xp.zeros(
             (field.shape[0], field.shape[1], z_axis.size),
             dtype=field.dtype,
         )
 
         # Convert the spectral image to the spatial field representation
         for i_m in range(grid.azimuthal_modes.size):
-            FieldAxprp.import_field(np.transpose(field[i_m]).copy())
+            FieldAxprp.import_field(xp.transpose(field[i_m]).copy())
 
             field_z[i_m] = prop[i_m].t2z(FieldAxprp.Field_ft, z_axis, z0=z0, t0=t0).T
 
-            field_z[i_m] *= np.exp(-1j * (z_axis / c + t0) * omega0)
+            field_z[i_m] *= xp.exp(-1j * (z_axis / c + t0) * omega0)
     else:
         # Construct the propagator
         Nx, Ny, Nt = field.shape
@@ -783,10 +787,10 @@ def export_to_z(dim, grid, omega0, z_axis=None, z0=0.0, t0=0.0, backend="NP"):
             verbose=False,
         )
         # Convert the spectral image to the spatial field representation
-        FieldAxprp.import_field(np.moveaxis(field, -1, 0).copy())
+        FieldAxprp.import_field(xp.moveaxis(field, -1, 0).copy())
         field_z = prop.t2z(FieldAxprp.Field_ft, z_axis, z0=z0, t0=t0)
-        field_z = np.moveaxis(field_z, 0, -1)
-        field_z *= np.exp(-1j * (z_axis / c + t0) * omega0)
+        field_z = xp.moveaxis(field_z, 0, -1)
+        field_z *= xp.exp(-1j * (z_axis / c + t0) * omega0)
 
     return field_z
 
@@ -831,10 +835,10 @@ def import_from_z(dim, grid, omega0, field_z, z_axis, z0=0.0, t0=0.0, backend="N
     Nz = z_axis.size
 
     # Transform the field from spatial to wavenumber domain
-    field_fft = np.fft.fft(field_z, axis=z_axis_indx, norm="forward")
+    field_fft = xp.fft.fft(field_z, axis=z_axis_indx, norm="forward")
 
     # Create the axes for wavenumbers, and for corresponding frequency
-    omega = 2 * np.pi * np.fft.fftfreq(Nz, dz / c) + omega0
+    omega = 2 * xp.pi * xp.fft.fftfreq(Nz, dz / c) + omega0
     k_z = omega / c
 
     if dim == "rt":
@@ -852,12 +856,12 @@ def import_from_z(dim, grid, omega0, field_z, z_axis, z0=0.0, t0=0.0, backend="N
             )
 
         # Convert the spectral image to the spatial field representation
-        field = np.zeros(grid.shape, dtype=np.complex128)
+        field = xp.zeros(grid.shape, dtype=xp.complex128)
         for i_m in range(grid.azimuthal_modes.size):
-            transform_data = np.transpose(field_fft[i_m]).copy()
-            transform_data *= np.exp(-1j * z_axis[0] * (k_z[:, None] - omega0 / c))
+            transform_data = xp.transpose(field_fft[i_m]).copy()
+            transform_data *= xp.exp(-1j * z_axis[0] * (k_z[:, None] - omega0 / c))
             field[i_m] = prop[i_m].z2t(transform_data, t_axis, z0=z0, t0=t0).T
-            field[i_m] *= np.exp(1j * (z0 / c + t_axis) * omega0)
+            field[i_m] *= xp.exp(1j * (z0 / c + t_axis) * omega0)
         grid.set_temporal_field(field)
     else:
         # Construct the propagator
@@ -872,8 +876,8 @@ def import_from_z(dim, grid, omega0, field_z, z_axis, z0=0.0, t0=0.0, backend="N
             verbose=False,
         )
         # Convert the spectral image to the spatial field representation
-        transform_data = np.moveaxis(field_fft, -1, 0).copy()
-        transform_data *= np.exp(-1j * z_axis[0] * (k_z[:, None, None] - omega0 / c))
-        field = np.moveaxis(prop.z2t(transform_data, t_axis, z0=z0, t0=t0), 0, -1)
-        field *= np.exp(1j * (z0 / c + t_axis) * omega0)
+        transform_data = xp.moveaxis(field_fft, -1, 0).copy()
+        transform_data *= xp.exp(-1j * z_axis[0] * (k_z[:, None, None] - omega0 / c))
+        field = xp.moveaxis(prop.z2t(transform_data, t_axis, z0=z0, t0=t0), 0, -1)
+        field *= xp.exp(1j * (z0 / c + t_axis) * omega0)
         grid.set_temporal_field(field)
diff --git a/lasy/utils/mode_decomposition.py b/lasy/utils/mode_decomposition.py
index b88fe7c0..75e01101 100644
--- a/lasy/utils/mode_decomposition.py
+++ b/lasy/utils/mode_decomposition.py
@@ -1,7 +1,6 @@
 import math
 
-import numpy as np
-
+from lasy.backend import xp
 from lasy.profiles.transverse.hermite_gaussian_profile import (
     HermiteGaussianTransverseProfile,
 )
@@ -71,9 +70,9 @@ def hermite_gauss_decomposition(laserProfile, n_x_max=12, n_y_max=12, res=1e-6):
     N_pts_y = int((hi[1] - lo[1]) / res)
 
     # Define spatial arrays
-    x = np.linspace(lo[0], hi[0], N_pts_x)
-    y = np.linspace(lo[1], hi[1], N_pts_y)
-    X, Y = np.meshgrid(x, y)
+    x = xp.linspace(lo[0], hi[0], N_pts_x)
+    y = xp.linspace(lo[1], hi[1], N_pts_y)
+    X, Y = xp.meshgrid(x, y)
     dx = x[2] - x[1]
     dy = y[2] - y[1]
 
@@ -88,8 +87,8 @@ def hermite_gauss_decomposition(laserProfile, n_x_max=12, n_y_max=12, res=1e-6):
     for i in range(n_x_max):
         for j in range(n_y_max):
             HGMode = HermiteGaussianTransverseProfile(w0, i, j)
-            coef = np.real(
-                np.sum(field * HGMode.evaluate(X, Y)) * dx * dy
+            coef = xp.real(
+                xp.sum(field * HGMode.evaluate(X, Y)) * dx * dy
             )  # modalDecomposition
             if math.isnan(coef):
                 coef = 0
@@ -125,23 +124,23 @@ def estimate_best_HG_waist(x, y, field):
     dy = y[1] - y[0]
     assert dx == dy
 
-    X, Y = np.meshgrid(x, y)
+    X, Y = xp.meshgrid(x, y)
 
-    D4SigX, D4SigY = find_d4sigma(np.abs(field) ** 2)
-    w0Est = np.mean((D4SigX, D4SigY)) / 2 * dx  # convert this to a 1/e^2 width
+    D4SigX, D4SigY = find_d4sigma(xp.abs(field) ** 2)
+    w0Est = xp.mean((D4SigX, D4SigY)) / 2 * dx  # convert this to a 1/e^2 width
 
     # Scan around the waist obtained from the D4sigma calculation,
     # and keep the waist for which this HG mode has the highest scalar
     # product with the input profile.
-    waistTest = np.linspace(w0Est / 2, w0Est * 1.5, 30)
-    coeffTest = np.zeros_like(waistTest)
+    waistTest = xp.linspace(w0Est / 2, w0Est * 1.5, 30)
+    coeffTest = xp.zeros_like(waistTest)
 
     for i, wTest in enumerate(waistTest):
         # create a gaussian
         HGMode = HermiteGaussianTransverseProfile(wTest, 0, 0)
         profile = HGMode.evaluate(X, Y)
-        coeffTest[i] = np.real(np.sum(profile * field))
-    w0 = waistTest[np.argmax(coeffTest)]
+        coeffTest[i] = xp.real(xp.sum(profile * field))
+    w0 = waistTest[xp.argmax(coeffTest)]
 
     print("Estimated w0 = %.2f microns" % (w0Est * 1e6))
     return w0
diff --git a/lasy/utils/openpmd_input.py b/lasy/utils/openpmd_input.py
index c683dd3b..a37b8a89 100644
--- a/lasy/utils/openpmd_input.py
+++ b/lasy/utils/openpmd_input.py
@@ -1,6 +1,7 @@
-import numpy as np
 from scipy.constants import c
 
+from lasy.backend import xp
+
 
 def reorder_array(array, md, dim):
     """Reorder an openPMD array to the lasy representation.
@@ -57,7 +58,7 @@ def reorder_array_xyt(array, md):
     if "z" in md.axes.values():
         t = (md.z - md.z[0]) / c
         # Flip to get complex envelope in t assuming z = -c*t
-        array = np.flip(array, axis=-1)
+        array = xp.flip(array, axis=-1)
     else:
         t = md.t
     axes = {"x": md.x, "y": md.y, "t": t}
@@ -94,7 +95,7 @@ def reorder_array_rt(array, md):
     if "z" in md.axes.values():
         t = (md.z - md.z[0]) / c
         # Flip to get complex envelope in t assuming z = -c*t
-        array = np.flip(array, axis=-1)
+        array = xp.flip(array, axis=-1)
     else:
         t = md.t
     r = md.r[md.r.size // 2 :]
@@ -102,6 +103,6 @@ def reorder_array_rt(array, md):
 
     array = 0.5 * (
         array[array.shape[0] // 2 :, :]
-        + np.flip(array[: array.shape[0] // 2, :], axis=0)
+        + xp.flip(array[: array.shape[0] // 2, :], axis=0)
     )
     return array, axes
diff --git a/lasy/utils/openpmd_output.py b/lasy/utils/openpmd_output.py
index b210ed5a..6f5648dc 100644
--- a/lasy/utils/openpmd_output.py
+++ b/lasy/utils/openpmd_output.py
@@ -59,7 +59,8 @@ def write_to_openpmd_file(
         Whether the envelope is converted to normalized vector potential
         before writing to file.
     """
-    array = grid.get_temporal_field()
+    # Get field on CPU
+    array = grid.get_temporal_field(to_cpu=True)
 
     # Create file
     full_filepath = os.path.join(
diff --git a/lasy/utils/phase_retrieval.py b/lasy/utils/phase_retrieval.py
index 4d9f915d..3221bf29 100644
--- a/lasy/utils/phase_retrieval.py
+++ b/lasy/utils/phase_retrieval.py
@@ -1,6 +1,6 @@
 import copy
 
-import numpy as np
+from lasy.backend import xp
 
 
 def gerchberg_saxton_algo(
@@ -51,10 +51,10 @@ def gerchberg_saxton_algo(
     """
     laser1 = copy.deepcopy(laserPos1)
     laser2 = copy.deepcopy(laserPos2)
-    amp1 = np.abs(laser1.grid.get_temporal_field())
-    amp1_summed = np.sum(amp1)
-    amp2 = np.abs(laser2.grid.get_temporal_field())
-    phase1 = np.zeros_like(amp1)
+    amp1 = xp.abs(laser1.grid.get_temporal_field())
+    amp1_summed = xp.sum(amp1)
+    amp2 = xp.abs(laser2.grid.get_temporal_field())
+    phase1 = xp.zeros_like(amp1)
 
     if condition == "max_iterations":
         breakout = lambda i: i < max_iterations
@@ -65,17 +65,17 @@ def gerchberg_saxton_algo(
 
     i = 0
     while breakout(cond):
-        laser1.grid.set_temporal_field(amp1 * np.exp(1j * phase1))
+        laser1.grid.set_temporal_field(amp1 * xp.exp(1j * phase1))
         laser1.propagate(dz, show_progress=False)
 
-        phase2 = np.angle(laser1.grid.get_temporal_field())
-        laser2.grid.set_temporal_field(amp2 * np.exp(1j * phase2))
+        phase2 = xp.angle(laser1.grid.get_temporal_field())
+        laser2.grid.set_temporal_field(amp2 * xp.exp(1j * phase2))
         laser2.propagate(-dz, show_progress=False)
 
         field2 = laser2.grid.get_temporal_field()
-        phase1 = np.angle(field2)
-        amp1_calc = np.abs(field2)
-        amp_error_summed = np.sum(np.abs(amp1_calc) - amp1)
+        phase1 = xp.angle(field2)
+        amp1_calc = xp.abs(field2)
+        amp_error_summed = xp.sum(xp.abs(amp1_calc) - amp1)
         if debug:
             i += 1
             print(
diff --git a/lasy/utils/zernike.py b/lasy/utils/zernike.py
index c6946508..2259ea13 100644
--- a/lasy/utils/zernike.py
+++ b/lasy/utils/zernike.py
@@ -2,6 +2,8 @@
 
 import numpy as np
 
+from lasy.backend import xp
+
 
 def get_zernike_nm(j):
     """
@@ -20,7 +22,7 @@ def get_zernike_nm(j):
     n,m : ints
         The standard Zernike Polynomial Indexes
     """
-    n = int(np.ceil((-3 + np.sqrt(9 + 8 * j)) / 2))
+    n = int(xp.ceil((-3 + xp.sqrt(9 + 8 * j)) / 2))
     m = 2 * j - n * (n + 2)
     return int(m), int(n)
 
@@ -50,8 +52,8 @@ def zernike(x, y, pupilCoords, j):
     """
     # Setup
     (cgx, cgy, r) = pupilCoords
-    rho = np.sqrt((x - cgx) ** 2 + (y - cgy) ** 2) / r
-    theta = np.arctan2(y - cgy, x - cgx)
+    rho = xp.sqrt((x - cgx) ** 2 + (y - cgy) ** 2) / r
+    theta = xp.arctan2(y - cgy, x - cgx)
 
     m, n = get_zernike_nm(j)
 
@@ -60,18 +62,18 @@ def zernike(x, y, pupilCoords, j):
 
     # Now multiply by the azimuthal part
     if m < 0:
-        Z = R * np.sin(-m * theta)
+        Z = R * xp.sin(-m * theta)
     else:
-        Z = R * np.cos(m * theta)
+        Z = R * xp.cos(m * theta)
 
     # Normalization
     if n == 0:
         scaling = 1
     else:
         if m == 0:
-            scaling = np.sqrt((n + 1))
+            scaling = xp.sqrt((n + 1))
         else:
-            scaling = np.sqrt(2 * (n + 1))
+            scaling = xp.sqrt(2 * (n + 1))
     Z = Z * scaling
 
     return Z
@@ -97,23 +99,23 @@ def RmnGenerator(n, m, rho):
     if n == 0:
         try:
             (r,) = rho.shape
-            Rmn = np.ones(
+            Rmn = xp.ones(
                 r,
             )
         except:
             r, c = rho.shape
-            Rmn = np.ones((r, c))
+            Rmn = xp.ones((r, c))
     elif (n - m) % 2 == 0:
         # Even, Rmn is not 0
         k = np.linspace(0, int((n - m) / 2), int((n - m) / 2) + 1).astype(int)
         try:
             (r,) = rho.shape
-            Rmn = np.zeros(
+            Rmn = xp.zeros(
                 r,
             )
         except:
             r, c = rho.shape
-            Rmn = np.zeros((r, c))
+            Rmn = xp.zeros((r, c))
         for i in k:
             Rmn = Rmn + ((-1) ** i * math.factorial(n - i)) / (
                 math.factorial(i)
@@ -124,11 +126,11 @@ def RmnGenerator(n, m, rho):
     else:
         try:
             (r,) = rho.shape
-            Rmn = np.zeros(
+            Rmn = xp.zeros(
                 r,
             )
         except:
             r, c = rho.shape
-            Rmn = np.zeros((r, c))
+            Rmn = xp.zeros((r, c))
 
     return Rmn
diff --git a/requirements.txt b/requirements.txt
index 6325cdaa..633a6ac0 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-axiprop
+axiprop @ git+https://github.com/hightower8083/axiprop@add-step-on-device
 numpy
 openpmd-api
 openpmd-viewer
diff --git a/tests/test_gaussian_propagator.py b/tests/test_gaussian_propagator.py
index e669eff6..cd89269c 100644
--- a/tests/test_gaussian_propagator.py
+++ b/tests/test_gaussian_propagator.py
@@ -3,6 +3,7 @@
 import numpy as np
 import pytest
 
+from lasy.backend import use_cupy, xp
 from lasy.laser import Laser
 from lasy.profiles.gaussian_profile import GaussianProfile
 
@@ -23,14 +24,18 @@ def gaussian():
 
 def get_w0(laser):
     # Calculate the laser waist
-    field = laser.grid.get_temporal_field()
+    field = laser.grid.get_temporal_field(to_cpu=True)
     if laser.dim == "xyt":
         Nx, Ny, Nt = field.shape
         A2 = (np.abs(field[Nx // 2 - 1, :, :]) ** 2).sum(-1)
         ax = laser.grid.axes[1]
+        if use_cupy:
+            ax = xp.asnumpy(ax)
     else:
         A2 = (np.abs(field[0, :, :]) ** 2).sum(-1)
         ax = laser.grid.axes[0]
+        if use_cupy:
+            ax = xp.asnumpy(ax)
         if ax[0] > 0:
             A2 = np.r_[A2[::-1], A2]
             ax = np.r_[-ax[::-1], ax]
diff --git a/tests/test_gsalgo.py b/tests/test_gsalgo.py
index bd913800..26cad070 100644
--- a/tests/test_gsalgo.py
+++ b/tests/test_gsalgo.py
@@ -5,6 +5,7 @@
 import numpy as np
 import pytest
 
+from lasy.backend import xp
 from lasy.laser import Laser
 from lasy.profiles.gaussian_profile import GaussianProfile
 from lasy.utils.phase_retrieval import gerchberg_saxton_algo
@@ -37,20 +38,20 @@ def test_3D_case(gaussian):
 
     # Add a phase aberration
     # CALCULATE THE REQUIRED PHASE ABERRATION
-    x = np.linspace(lo[0], hi[0], npoints[0])
-    y = np.linspace(lo[1], hi[1], npoints[1])
-    X, Y = np.meshgrid(x, y)
+    x = xp.linspace(lo[0], hi[0], npoints[0])
+    y = xp.linspace(lo[1], hi[1], npoints[1])
+    X, Y = xp.meshgrid(x, y)
     pupilRadius = 2 * w0
     phase = -0.2 * zernike(X, Y, (0, 0, pupilRadius), 3)
 
-    R = np.sqrt(X**2 + Y**2)
-    phaseMask = np.ones_like(phase)
+    R = xp.sqrt(X**2 + Y**2)
+    phaseMask = xp.ones_like(phase)
     phaseMask[R > pupilRadius] = 0
 
     # NOW ADD THE PHASE TO EACH SLICE OF THE FOCUS
-    phase3D = np.repeat(phase[:, :, np.newaxis], npoints[2], axis=2)
+    phase3D = xp.repeat(phase[:, :, xp.newaxis], npoints[2], axis=2)
     field = laser.grid.get_temporal_field()
-    laser.grid.set_temporal_field(np.abs(field) * np.exp(1j * phase3D))
+    laser.grid.set_temporal_field(xp.abs(field) * xp.exp(1j * phase3D))
 
     # PROPAGATE THE FIELD FIELD FOWARDS AND BACKWARDS BY 1 MM
     propDist = 2e-3
diff --git a/tests/test_parabolic_mirror.py b/tests/test_parabolic_mirror.py
index 1c60a4b1..45e1baa3 100644
--- a/tests/test_parabolic_mirror.py
+++ b/tests/test_parabolic_mirror.py
@@ -9,6 +9,7 @@
 
 import numpy as np
 
+from lasy.backend import use_cupy, xp
 from lasy.laser import Laser
 from lasy.optical_elements import ParabolicMirror
 from lasy.profiles.gaussian_profile import GaussianProfile
@@ -26,14 +27,18 @@
 
 def get_w0(laser):
     # Calculate the laser waist
-    field = laser.grid.get_temporal_field()
+    field = laser.grid.get_temporal_field(to_cpu=True)
     if laser.dim == "xyt":
-        Nx = field.shape[0]
+        Nx, Ny, Nt = field.shape
         A2 = (np.abs(field[Nx // 2 - 1, :, :]) ** 2).sum(-1)
         ax = laser.grid.axes[1]
+        if use_cupy:
+            ax = xp.asnumpy(ax)
     else:
         A2 = (np.abs(field[0, :, :]) ** 2).sum(-1)
         ax = laser.grid.axes[0]
+        if use_cupy:
+            ax = xp.asnumpy(ax)
         if ax[0] > 0:
             A2 = np.r_[A2[::-1], A2]
             ax = np.r_[-ax[::-1], ax]