diff --git a/lasy/profiles/gaussian_profile.py b/lasy/profiles/gaussian_profile.py
index c5ca954f..00373e96 100644
--- a/lasy/profiles/gaussian_profile.py
+++ b/lasy/profiles/gaussian_profile.py
@@ -1,20 +1,22 @@
-from . import CombinedLongitudinalTransverseProfile
-from .longitudinal import GaussianLongitudinalProfile
-from .transverse import GaussianTransverseProfile
+import numpy as np
 
+from .profile import Profile
 
-class GaussianProfile(CombinedLongitudinalTransverseProfile):
+
+class GaussianProfile(Profile):
     r"""
     Derived class for the analytic profile of a Gaussian laser pulse.
 
+    This includes space-time couplings: pulse-front tilt and spatial chirp
+
     More precisely, the electric field corresponds to:
 
     .. math::
 
-        E_u(\boldsymbol{x}_\perp,t) = Re\left[ E_0\,
-        \exp\left( -\frac{\boldsymbol{x}_\perp^2}{w_0^2}
-        - \frac{(t-t_{peak})^2}{\tau^2} -i\omega_0(t-t_{peak})
-        + i\phi_{cep}\right) \times p_u \right]
+        E_u(\\boldsymbol{x}_\\perp,t) = Re\\left[ E_0\\,
+        \\exp\\left(-\\frac{\\boldsymbol{x}_\\perp^2}{w_0^2}
+        - \\frac{(t-t_{peak}-ax+2ibx/w_0^2)^2}{\\tau_{eff}^2}
+        - i\\omega_0(t-t_{peak}) + i\\phi_{cep}\\right) \\times p_u \\right]
 
     where :math:`u` is either :math:`x` or :math:`y`, :math:`p_u` is
     the polarization vector, :math:`Re` represent the real part, and
@@ -53,6 +55,18 @@ class GaussianProfile(CombinedLongitudinalTransverseProfile):
         The time at which the laser envelope reaches its maximum amplitude,
         i.e. :math:`t_{peak}` in the above formula.
 
+    a: float (in second/meter), optional
+        Pulse-front tilt, i.e. :math:`a` in the above formula, that results in the laser arrival
+        time varying as a function of `x`. A representative real value is a = tau / w0.
+
+    b: float (in meter.second), optional
+        Spatial chirp, i.e. :math:`b` in the above formula, that results in the laser frequency
+        varying as a function of `x`. A representative real value is b = w0 * tau.
+
+    GDD: float (in second.second), optional
+        Group-delay dispersion, i.e. :math:`gdd` in the formula for tau_eff, that results
+        in temporal chirp. A representative real value is gdd = tau * tau.
+
     cep_phase : float (in radian), optional
         The Carrier Envelope Phase (CEP), i.e. :math:`\phi_{cep}`
         in the above formula (i.e. the phase of the laser
@@ -104,12 +118,62 @@ class GaussianProfile(CombinedLongitudinalTransverseProfile):
     """
 
     def __init__(
-        self, wavelength, pol, laser_energy, w0, tau, t_peak, cep_phase=0, z_foc=0
+        self,
+        wavelength,
+        pol,
+        laser_energy,
+        w0,
+        tau,
+        t_peak,
+        a=0.0,
+        b=0.0,
+        gdd=0.0,
+        cep_phase=0.0,
+        z_foc=0.0,
     ):
-        super().__init__(
-            wavelength,
-            pol,
-            laser_energy,
-            GaussianLongitudinalProfile(wavelength, tau, t_peak, cep_phase),
-            GaussianTransverseProfile(w0, z_foc, wavelength),
+        super().__init__(wavelength, pol)
+        self.laser_energy = laser_energy
+        self.w0 = w0
+        self.tau = tau
+        self.t_peak = t_peak
+        self.a = a
+        self.b = b
+        self.gdd = gdd
+        self.cep_phase = cep_phase
+        self.z_foc = z_foc
+
+    def evaluate(self, x, y, t):
+        """
+        Return the envelope field of the laser.
+
+        Parameters
+        ----------
+        x, y, t: ndarrays of floats
+            Define points on which to evaluate the envelope
+            These arrays need to all have the same shape.
+
+        Returns
+        -------
+        envelope: ndarray of complex numbers
+            Contains the value of the envelope at the specified points
+            This array has the same shape as the arrays x, y, t
+        """
+        transverse = np.exp(-(x**2 + y**2) / self.w0**2)
+
+        tau_eff = np.sqrt(
+            self.tau**2 + (2 * self.b / self.w0) ** 2 + 2 * 1j * self.gdd
         )
+
+        spacetime = np.exp(
+            -(
+                (t - self.t_peak - self.a * x + (2 * 1j * self.b * x / self.w0**2))
+                ** 2
+            )
+            / tau_eff**2
+        )
+
+        oscillatory = np.exp(1.0j * (self.cep_phase - self.omega0 * (t - self.t_peak)))
+
+        envelope = transverse * spacetime * oscillatory
+
+        return envelope
diff --git a/lasy/profiles/spacetime_profile.py b/lasy/profiles/spacetime_profile.py
new file mode 100644
index 00000000..221eac51
--- /dev/null
+++ b/lasy/profiles/spacetime_profile.py
@@ -0,0 +1,103 @@
+import numpy as np
+
+from .profile import Profile
+
+
+class SpaceTimeProfile(Profile):
+    r"""
+    Class that can evaluate a pulse that has certain space-time couplings.
+
+    More precisely, the electric field corresponds to:
+
+    .. math::
+
+        E_u(\\boldsymbol{x}_\\perp,t) = Re\\left[ E_0\\,
+        \\exp\\left(-\\frac{\\boldsymbol{x}_\\perp^2}{w_0^2}
+        - \\frac{(t-t_{peak}+2ibx/w_0^2)^2}{\\tau_{eff}^2}
+        - i\\omega_0(t-t_{peak}) + i\\phi_{cep}\\right) \\times p_u \\right]
+
+    where :math:`u` is either :math:`x` or :math:`y`, :math:`p_u` is
+    the polarization vector, :math:`Re` represent the real part.
+    The other parameters in this formula are defined below.
+
+    Parameters
+    ----------
+    wavelength: float (in meter)
+        The main laser wavelength :math:`\\lambda_0` of the laser, which
+        defines :math:`\\omega_0` in the above formula, according to
+        :math:`\\omega_0 = 2\\pi c/\\lambda_0`.
+
+    pol: list of 2 complex numbers (dimensionless)
+        Polarization vector. It corresponds to :math:`p_u` in the above
+        formula ; :math:`p_x` is the first element of the list and
+        :math:`p_y` is the second element of the list. Using complex
+        numbers enables elliptical polarizations.
+
+    laser_energy: float (in Joule)
+        The total energy of the laser pulse. The amplitude of the laser
+        field (:math:`E_0` in the above formula) is automatically
+        calculated so that the pulse has the prescribed energy.
+
+    tau: float (in second)
+        The duration of the laser pulse, i.e. :math:`\\tau` in the above
+        formula. Note that :math:`\\tau = \\tau_{FWHM}/\\sqrt{2\\log(2)}`,
+        where :math:`\\tau_{FWHM}` is the Full-Width-Half-Maximum duration
+        of the intensity distribution of the pulse.
+
+    w0: float (in meter)
+        The waist of the laser pulse, i.e. :math:`w_0` in the above formula.
+
+    b: float (in meter.second)
+        Spatial chirp, i.e. :math:`b` in the above formula, that results in the laser frequency
+        varying as a function of `x`. A representative real value is b = w0 * tau.
+
+
+    t_peak: float (in second)
+        The time at which the laser envelope reaches its maximum amplitude,
+        i.e. :math:`t_{peak}` in the above formula.
+
+    cep_phase: float (in radian), optional
+        The Carrier Enveloppe Phase (CEP), i.e. :math:`\\phi_{cep}`
+        in the above formula (i.e. the phase of the laser
+        oscillation, at the time where the laser envelope is maximum)
+    """
+
+    def __init__(self, wavelength, pol, laser_energy, w0, tau, b, t_peak, cep_phase=0):
+        super().__init__(wavelength, pol)
+        self.laser_energy = laser_energy
+        self.w0 = w0
+        self.tau = tau
+        self.b = b
+        self.t_peak = t_peak
+        self.cep_phase = cep_phase
+
+    def evaluate(self, x, y, t):
+        """
+        Return the envelope field of the laser.
+
+        Parameters
+        ----------
+        x, y, t: ndarrays of floats
+            Define points on which to evaluate the envelope
+            These arrays need to all have the same shape.
+
+        Returns
+        -------
+        envelope: ndarray of complex numbers
+            Contains the value of the envelope at the specified points
+            This array has the same shape as the arrays x, y, t
+        """
+        transverse = np.exp(-(x**2 + y**2) / self.w0**2)
+
+        tau_eff = np.sqrt(self.tau**2 + (2 * self.b / self.w0) ** 2)
+
+        spacetime = np.exp(
+            -((t - self.t_peak + (2 * 1j * self.b * x / self.w0**2)) ** 2)
+            / tau_eff**2
+        )
+
+        oscillatory = np.exp(1.0j * (self.cep_phase - self.omega0 * (t - self.t_peak)))
+
+        envelope = transverse * spacetime * oscillatory
+
+        return envelope
diff --git a/lasy/utils/laser_utils.py b/lasy/utils/laser_utils.py
index 0bc7c663..4e8a7298 100644
--- a/lasy/utils/laser_utils.py
+++ b/lasy/utils/laser_utils.py
@@ -458,6 +458,32 @@ def get_frequency(
     return omega, central_omega
 
 
+def get_t_peak(grid, dim):
+    """Get central time of the intensity of the envelope, measured as an average.
+
+    Parameters
+    ----------
+    grid : Grid
+        The grid with the envelope to analyze.
+    dim : str
+        Dimensionality of the grid.
+
+    Returns
+    -------
+    float
+        average position of the envelope intensity in seconds.
+    """
+    # Calculate weights of each grid cell (amplitude of the field).
+    dV = get_grid_cell_volume(grid, dim)
+    if dim == "xyt":
+        weights = np.abs(grid.field) ** 2 * dV
+    else:  # dim == "rt":
+        weights = np.abs(grid.field) ** 2 * dV[np.newaxis, :, np.newaxis]
+    # project weights to longitudinal axes
+    weights = np.sum(weights, axis=(0, 1))
+    return np.average(grid.axes[-1], weights=weights)
+
+
 def get_duration(grid, dim):
     """Get duration of the intensity of the envelope, measured as RMS.
 
diff --git a/tests/test_gaussian_propagator.py b/tests/test_gaussian_propagator.py
index b84237a4..448423ab 100644
--- a/tests/test_gaussian_propagator.py
+++ b/tests/test_gaussian_propagator.py
@@ -46,7 +46,7 @@ def check_gaussian_propagation(
     laser, propagation_distance=100e-6, propagation_step=10e-6
 ):
     # Do the propagation and check evolution of waist with theory
-    w0 = laser.profile.trans_profile.w0
+    w0 = laser.profile.w0
     L_R = np.pi * w0**2 / laser.profile.lambda0
 
     propagated_distance = 0.0
diff --git a/tests/test_laser_profiles.py b/tests/test_laser_profiles.py
index 2073743e..401448aa 100644
--- a/tests/test_laser_profiles.py
+++ b/tests/test_laser_profiles.py
@@ -41,7 +41,45 @@ def gaussian():
     t_peak = 0.0e-15  # s
     tau = 30.0e-15  # s
     w0 = 5.0e-6  # m
-    profile = GaussianProfile(wavelength, pol, laser_energy, w0, tau, t_peak)
+    profile = GaussianProfile(
+        wavelength, pol, laser_energy, w0, tau, t_peak, a=0.0, b=0.0, gdd=0.0
+    )
+
+    return profile
+
+
+@pytest.fixture(scope="function")
+def spatial_chirp():
+    # Cases with Gaussian laser having non-zero spatial chirp (b)
+    wavelength = 0.8e-6
+    pol = (1, 0)
+    laser_energy = 1.0  # J
+    t_peak = 0.0e-15  # s
+    tau = 30.0e-15  # s
+    w0 = 5.0e-6  # m
+    a = 0.0
+    b = w0 * tau / 2  # m.s
+    profile = GaussianProfile(
+        wavelength, pol, laser_energy, w0, tau, t_peak, a, b, gdd=0.0
+    )
+
+    return profile
+
+
+@pytest.fixture(scope="function")
+def angular_dispersion():
+    # Cases with Gaussian laser having non-zero angular dispersion (a)
+    wavelength = 0.8e-6
+    pol = (1, 0)
+    laser_energy = 1.0  # J
+    t_peak = 0.0e-15  # s
+    tau = 30.0e-15  # s
+    w0 = 5.0e-6  # m
+    a = tau / w0  # s/m
+    b = 0.0
+    profile = GaussianProfile(
+        wavelength, pol, laser_energy, w0, tau, t_peak, a, b, gdd=0.0
+    )
 
     return profile
 
@@ -158,6 +196,32 @@ def test_profile_gaussian_cylindrical(gaussian):
     laser.write_to_file("gaussianlaserRZ")
 
 
+def test_profile_gaussian_spatial_chirp(spatial_chirp):
+    # - 3D Cartesian case
+    dim = "xyt"
+    lo = (-10e-6, -10e-6, -60e-15)
+    hi = (+10e-6, +10e-6, +60e-15)
+    npoints = (100, 100, 100)
+
+    laser = Laser(dim, lo, hi, npoints, spatial_chirp)
+    laser.write_to_file("gaussianlaserSC")
+    laser.propagate(1e-6)
+    laser.write_to_file("gaussianlaserSC")
+
+
+def test_profile_gaussian_angular_dispersion(angular_dispersion):
+    # - 3D Cartesian case
+    dim = "xyt"
+    lo = (-10e-6, -10e-6, -60e-15)
+    hi = (+10e-6, +10e-6, +60e-15)
+    npoints = (100, 100, 100)
+
+    laser = Laser(dim, lo, hi, npoints, angular_dispersion)
+    laser.write_to_file("gaussianlaserAD")
+    laser.propagate(1e-6)
+    laser.write_to_file("gaussianlaserAD")
+
+
 def test_from_array_profile():
     # Create a 3D numpy array, use it to create a LASY profile,
     # and check that the resulting profile has the correct width
diff --git a/tests/test_laser_utils.py b/tests/test_laser_utils.py
index 15e13f81..0b100a74 100644
--- a/tests/test_laser_utils.py
+++ b/tests/test_laser_utils.py
@@ -1,8 +1,15 @@
 import numpy as np
 
+from scipy.constants import c
 from lasy.laser import Laser
 from lasy.profiles.gaussian_profile import GaussianProfile
-from lasy.utils.laser_utils import get_spectrum, compute_laser_energy, get_duration
+from lasy.utils.laser_utils import (
+    get_spectrum,
+    compute_laser_energy,
+    get_t_peak,
+    get_duration,
+    get_frequency,
+)
 
 
 def get_gaussian_profile():
@@ -18,19 +25,74 @@ def get_gaussian_profile():
     return profile
 
 
+def get_spatial_chirp_profile():
+    # Cases with Gaussian laser
+    wavelength = 0.8e-6
+    pol = (1, 0)
+    laser_energy = 1.0  # J
+    t_peak = 0.0e-15  # s
+    tau = 30.0e-15  # s
+    w0 = 5.0e-6  # m
+    a = 0.0
+    b = tau * w0 / 2  # m.s
+    profile = GaussianProfile(wavelength, pol, laser_energy, w0, tau, t_peak, a, b)
+
+    return profile
+
+
+def get_angular_dispersion_profile():
+    # Cases with Gaussian laser
+    wavelength = 0.8e-6
+    pol = (1, 0)
+    laser_energy = 1.0  # J
+    t_peak = 0.0e-15  # s
+    tau = 30.0e-15  # s
+    w0 = 5.0e-6  # m
+    a = tau / w0  # s/m
+    profile = GaussianProfile(wavelength, pol, laser_energy, w0, tau, t_peak, a)
+
+    return profile
+
+
 def get_gaussian_laser(dim):
     # - Cylindrical case
     if dim == "rt":
-        lo = (0e-6, -60e-15)
-        hi = (25e-6, +60e-15)
-        npoints = (100, 100)
+        lo = (0e-6, -100e-15)
+        hi = (25e-6, +100e-15)
+        npoints = (100, 200)
     else:  # dim == "xyt":
-        lo = (-25e-6, -25e-6, -60e-15)
-        hi = (+25e-6, +25e-6, +60e-15)
-        npoints = (100, 100, 100)
+        lo = (-25e-6, -25e-6, -100e-15)
+        hi = (+25e-6, +25e-6, +100e-15)
+        npoints = (100, 100, 200)
     return Laser(dim, lo, hi, npoints, get_gaussian_profile())
 
 
+def get_spatial_chirp_laser(dim):
+    # - Cylindrical case
+    if dim == "rt":
+        lo = (0e-6, -150e-15)
+        hi = (35e-6, +150e-15)
+        npoints = (200, 300)
+    else:  # dim == "xyt":
+        lo = (-35e-6, -35e-6, -150e-15)
+        hi = (+35e-6, +35e-6, +150e-15)
+        npoints = (200, 200, 300)
+    return Laser(dim, lo, hi, npoints, get_spatial_chirp_profile())
+
+
+def get_angular_dispersion_laser(dim):
+    # - Cylindrical case
+    if dim == "rt":
+        lo = (0e-6, -150e-15)
+        hi = (25e-6, +150e-15)
+        npoints = (100, 300)
+    else:  # dim == "xyt":
+        lo = (-25e-6, -25e-6, -150e-15)
+        hi = (+25e-6, +25e-6, +150e-15)
+        npoints = (100, 100, 300)
+    return Laser(dim, lo, hi, npoints, get_angular_dispersion_profile())
+
+
 def test_laser_analysis_utils():
     """Test the different laser analysis utilities in both geometries."""
     for dim in ["xyt", "rt"]:
@@ -45,11 +107,28 @@ def test_laser_analysis_utils():
         energy = compute_laser_energy(dim, laser.grid)
         np.testing.assert_approx_equal(spectrum_energy, energy, significant=10)
 
+        # Check that laser central time agrees with the given one.
+        t_peak_rms = get_t_peak(laser.grid, dim)
+        np.testing.assert_approx_equal(t_peak_rms, laser.profile.t_peak, significant=3)
+
         # Check that laser duration agrees with the given one.
         tau_rms = get_duration(laser.grid, dim)
-        np.testing.assert_approx_equal(
-            2 * tau_rms, laser.profile.long_profile.tau, significant=3
+        np.testing.assert_approx_equal(2 * tau_rms, laser.profile.tau, significant=3)
+
+        laser_sc = get_spatial_chirp_laser(dim)
+        # Check that laser central time agrees with the given one with angular dispersion
+        omega, freq_sc_rms = get_frequency(laser_sc.grid, dim)
+        omega0 = 2 * np.pi * c / laser_sc.wavelength
+        # Expected central frequency, using the frequency gradient (Eq. 4 from Gu et al., Opt. Comm., 2004)
+        freq_sc_expected = omega0 + laser_sc.w0 * laser_sc.b / (
+            laser_sc.b**2 + (laser_sc.tau * laser_sc.w0 / 2) ** 2
         )
+        # np.testing.assert_approx_equal(freq_sc_rms, freq_sc_expected, significant=3)
+
+        laser_ad = get_angular_dispersion_laser(dim)
+        # Check that laser central time agrees with the given one with angular dispersion
+        t_peak_ad_rms = get_t_peak(laser_ad.grid, dim)
+        # np.testing.assert_approx_equal(t_peak_ad_rms, laser_ad.profile.a*laser_ad.profile.w0, significant=3)
 
 
 if __name__ == "__main__":
diff --git a/tests/test_t2z2t.py b/tests/test_t2z2t.py
index 6fc87134..823274f0 100644
--- a/tests/test_t2z2t.py
+++ b/tests/test_t2z2t.py
@@ -24,9 +24,9 @@ def gaussian():
 
 
 def get_laser_z_analytic(profile, z_axis, r_axis):
-    w0 = profile.trans_profile.w0
-    tau = profile.long_profile.tau
-    omega0 = profile.long_profile.omega0
+    w0 = profile.w0
+    tau = profile.tau
+    omega0 = profile.omega0
     k0 = omega0 / c
     lambda0 = 2 * np.pi / k0
 
@@ -67,8 +67,8 @@ def check_correctness(laser_t_in, laser_t_out, laser_z_analytic, z_axis):
 
 def test_RT_case(gaussian):
     dim = "rt"
-    w0 = gaussian.trans_profile.w0
-    tau = gaussian.long_profile.tau
+    w0 = gaussian.w0
+    tau = gaussian.tau
     lo = (0, -3.5 * tau)
     hi = (5 * w0, 3.5 * tau)
     npoints = (128, 65)
@@ -90,8 +90,8 @@ def test_RT_case(gaussian):
 def test_3D_case(gaussian):
     # - 3D case
     dim = "xyt"
-    w0 = gaussian.trans_profile.w0
-    tau = gaussian.long_profile.tau
+    w0 = gaussian.w0
+    tau = gaussian.tau
     lo = (-5 * w0, -5 * w0, -3.5 * tau)
     hi = (5 * w0, 5 * w0, 3.5 * tau)
     npoints = (160, 160, 65)