diff --git a/dyson/__init__.py b/dyson/__init__.py
index e6b7aa5..82010e7 100644
--- a/dyson/__init__.py
+++ b/dyson/__init__.py
@@ -71,6 +71,8 @@
        self-energy.
    * - :class:`~dyson.solvers.dynamic.cpgf.CPGF`
      - Chebyshev polynomial moments of the dynamic Green's function.
+   * - :class:`~dyson.solvers.dynamic.direct.Direct`
+     - Dynamic self-energy and initial Green's function, in frequency space.
 
 For a full accounting of the inputs and their types, please see the documentation for each solver.
 
@@ -135,5 +137,6 @@
     DensityRelaxation,
     CorrectionVector,
     CPGF,
+    Direct,
 )
 from dyson.expressions import HF, CCSD, FCI, ADC2, ADC2x, Hamiltonian
diff --git a/dyson/grids/__init__.py b/dyson/grids/__init__.py
index 887cf98..4d0622e 100644
--- a/dyson/grids/__init__.py
+++ b/dyson/grids/__init__.py
@@ -11,7 +11,73 @@
 
     grid
     frequency
+    time
+    fourier
+    pade
+    util
 """
 
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, TypeVar
+
+from dyson import numpy as np
+from dyson.grids.grid import BaseGrid
 from dyson.grids.frequency import RealFrequencyGrid, GridRF
 from dyson.grids.frequency import ImaginaryFrequencyGrid, GridIF
+from dyson.grids.time import RealTimeGrid, GridRT
+from dyson.grids.time import ImaginaryTimeGrid, GridIT
+from dyson.grids.fourier import fourier_transform_imag
+from dyson.grids.pade import analytic_continuation_freq_pade
+from dyson.grids.util import are_dual, are_equal
+
+if TYPE_CHECKING:
+    from typing import Any
+
+    from dyson.representations import Dynamic
+
+GridSrcT = TypeVar("GridSrcT", bound=BaseGrid)
+GridDstT = TypeVar("GridDstT", bound=BaseGrid)
+
+
+def transform(dynamic: Dynamic[GridSrcT], grid: GridDstT, **kwargs: Any) -> Dynamic[GridDstT]:
+    """Transform a dynamic quantity to a new grid using either FFT or AC.
+
+    Currently available transformations are:
+
+    .. code-block:: bash
+
+                   AC
+               ─────────>
+        GridRF <───────── GridIF
+                   AC
+                           │ ^
+                           │ │
+                      IFFT │ │ FFT
+                           │ │
+                           v │
+
+                         GridIT
+
+    Args:
+        dynamic: The dynamic quantity to transform.
+        grid: The grid to transform to.
+        **kwargs: Additional keyword arguments passed to the transformation function.
+
+    Returns:
+        The transformed dynamic quantity.
+
+    Raises:
+        NotImplementedError: If the transformation is not implemented.
+    """
+    if isinstance(dynamic.grid, GridIT) and isinstance(grid, GridIF):
+        return fourier_transform_imag(dynamic, grid, **kwargs)  # type: ignore
+    if isinstance(dynamic.grid, GridIF) and isinstance(grid, GridIT):
+        return fourier_transform_imag(dynamic, grid, **kwargs)  # type: ignore
+    if isinstance(dynamic.grid, GridIF) and isinstance(grid, GridRF):
+        return analytic_continuation_freq_pade(dynamic, grid, **kwargs)  # type: ignore
+    if isinstance(dynamic.grid, GridRF) and isinstance(grid, GridIF):
+        return analytic_continuation_freq_pade(dynamic, grid, **kwargs)  # type: ignore
+    raise NotImplementedError(
+        f"transformation between {type(dynamic.grid)} and {type(grid)} not implemented"
+    )
diff --git a/dyson/grids/fourier.py b/dyson/grids/fourier.py
new file mode 100644
index 0000000..46d4087
--- /dev/null
+++ b/dyson/grids/fourier.py
@@ -0,0 +1,81 @@
+"""Fourier transformation."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from dyson import numpy as np
+from dyson import util
+from dyson.grids.util import are_dual
+from dyson.representations.enums import Component, Reduction
+
+if TYPE_CHECKING:
+    from dyson.grids.grid import BaseImaginaryGrid
+    from dyson.representations.dynamic import Dynamic
+
+
+def fourier_transform_imag(
+    greens_function: Dynamic[BaseImaginaryGrid], grid: BaseImaginaryGrid
+) -> Dynamic[BaseImaginaryGrid]:
+    """Fourier transform between imaginary frequency and imaginary time grids.
+
+    Args:
+        greens_function: Dynamic quantity in imaginary frequency or time domain.
+        grid: Target grid for the Fourier transform.
+
+    Returns:
+        Dynamic quantity in the target domain.
+    """
+    grid_in, grid_out = greens_function.grid, grid
+    if not are_dual(grid_in, grid_out):
+        raise ValueError("the two grids must be dual to each other.")
+    if grid_in.reality or grid_out.reality:
+        raise ValueError("only imaginary frequency and imaginary time grids is supported.")
+    if greens_function.component != Component.FULL:
+        raise ValueError("only full component is supported.")
+    if greens_function.reduction == Reduction.TRACE:
+        raise ValueError("traced reduction is not supported.")
+    forward = grid_in.domain == "time"
+
+    # Setup based on direction
+    beta = grid_in.beta
+    if forward:
+        sign = 1
+        norm = beta
+        fft = np.fft.ifft
+    else:
+        sign = -1
+        norm = 2 / beta
+        fft = np.fft.fft
+
+    # Get the shifts
+    shifts = (
+        np.exp(1.0j * sign * np.min(grid_out.points) * grid_in.points),
+        np.exp(1.0j * sign * np.min(grid_in.points) * (grid_out.points - np.min(grid_out.points))),
+    )
+
+    # Get the input array
+    array = greens_function.array.copy()
+
+    # Perform the Fourier transform
+    array = util.einsum("w...,w->w...", array, shifts[0])
+    array = fft(array, max(len(grid_in), len(grid_out)), axis=0)[: len(grid_out)]
+    array = util.einsum("w...,w->w...", array, shifts[1])
+
+    # Normalise
+    array *= norm
+
+    # Hermitise for imaginary frequency to imaginary time transform
+    if not forward:
+        if greens_function.reduction == Reduction.NONE:
+            array = 0.5 * (array + array.swapaxes(1, 2).conj())
+        else:
+            array = 0.5 * (array + array.conj())
+
+    return greens_function.__class__(
+        grid_out,
+        array,
+        reduction=greens_function.reduction,
+        ordering=greens_function.ordering,
+        hermitian=greens_function.hermitian,
+    )
diff --git a/dyson/grids/frequency.py b/dyson/grids/frequency.py
index fc13f28..dafcecb 100644
--- a/dyson/grids/frequency.py
+++ b/dyson/grids/frequency.py
@@ -6,89 +6,22 @@
 from typing import TYPE_CHECKING
 
 import scipy.special
+from typing_extensions import Self
 
 from dyson import numpy as np
-from dyson import util
-from dyson.grids.grid import BaseGrid
-from dyson.representations.enums import Component, Ordering, Reduction
+from dyson.grids.grid import BaseGrid, BaseImaginaryGrid, BaseRealGrid
+from dyson.representations.enums import Ordering
 
 if TYPE_CHECKING:
     from typing import Any
 
-    from dyson.representations.dynamic import Dynamic
-    from dyson.representations.lehmann import Lehmann
+    from dyson.grids.time import ImaginaryTimeGrid, RealTimeGrid
     from dyson.typing import Array
 
 
 class BaseFrequencyGrid(BaseGrid):
     """Base class for frequency grids."""
 
-    def evaluate_lehmann(
-        self,
-        lehmann: Lehmann,
-        reduction: Reduction = Reduction.NONE,
-        component: Component = Component.FULL,
-        **kwargs: Any,
-    ) -> Dynamic[BaseFrequencyGrid]:
-        r"""Evaluate a Lehmann representation on the grid.
-
-        The imaginary frequency representation is defined as
-
-        .. math::
-            \sum_{k} \frac{v_{pk} u_{qk}^*}{i \omega - \epsilon_k},
-
-        and the real frequency representation is defined as
-
-        .. math::
-            \sum_{k} \frac{v_{pk} u_{qk}^*}{\omega - \epsilon_k \pm i \eta},
-
-        where :math:`\omega` is the frequency grid, :math:`\epsilon_k` are the poles, and the sign
-        of the broadening factor is determined by the time ordering.
-
-        Args:
-            lehmann: Lehmann representation to evaluate.
-            reduction: The reduction of the dynamic representation.
-            component: The component of the dynamic representation.
-            kwargs: Additional keyword arguments for the resolvent.
-
-        Returns:
-            Lehmann representation, realised on the grid.
-        """
-        from dyson.representations.dynamic import Dynamic  # noqa: PLC0415
-
-        left, right = lehmann.unpack_couplings()
-        resolvent = self.resolvent(lehmann.energies, lehmann.chempot, **kwargs)
-        reduction = Reduction(reduction)
-        component = Component(component)
-
-        # Get the input and output indices based on the reduction type
-        inp = "qk"
-        out = "wpq"
-        if reduction == reduction.NONE:
-            pass
-        elif reduction == reduction.DIAG:
-            inp = "pk"
-            out = "wp"
-        elif reduction == reduction.TRACE:
-            inp = "pk"
-            out = "w"
-        else:
-            reduction.raise_invalid_representation()
-
-        # Perform the downfolding operation
-        array = util.einsum(f"pk,{inp},wk->{out}", right, left.conj(), resolvent)
-
-        # Get the required component
-        # TODO: Save time by not evaluating the full array when not needed
-        if component == Component.REAL:
-            array = array.real
-        elif component == Component.IMAG:
-            array = array.imag
-
-        return Dynamic(
-            self, array, reduction=reduction, component=component, hermitian=lehmann.hermitian
-        )
-
     @property
     def domain(self) -> str:
         """Get the domain of the grid.
@@ -100,26 +33,50 @@ def domain(self) -> str:
 
     @abstractmethod
     def resolvent(  # noqa: D417
-        self, energies: Array, chempot: float | Array, **kwargs: Any
+        self,
+        energies: Array,
+        chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
+        invert: bool = True,
     ) -> Array:
-        """Get the resolvent of the grid.
+        """Get the resolvent of a Lehmann representation on the grid.
 
         Args:
             energies: Energies of the poles.
             chempot: Chemical potential.
+            ordering: Time ordering of the resolvent.
+            invert: Whether to apply the inversion in the resolvent formula.
 
         Returns:
-            Resolvent of the grid.
+            Resolvent on the grid.
         """
         pass
 
+    def _lehmann_kernel(
+        self,
+        energies: Array,
+        chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
+    ) -> Array:
+        """Get the kernel of a Lehmann representation on the grid.
 
-class RealFrequencyGrid(BaseFrequencyGrid):
-    """Real frequency grid."""
+        Args:
+            energies: Energies of the poles.
+            chempot: Chemical potential.
+            ordering: Time ordering of the resolvent.
 
-    eta: float = 1e-2
+        Returns:
+            Kernel of a Lehmann representation on the grid.
+
+        Note:
+            The kernel is a hook to generalise the resolvent or propagator, depending on whether
+            the grid is in the frequency or time domain.
+        """
+        return self.resolvent(energies, chempot, ordering=ordering)
 
-    _options = {"eta"}
+
+class RealFrequencyGrid(BaseFrequencyGrid, BaseRealGrid):
+    """Real frequency grid."""
 
     def __init__(  # noqa: D417
         self, points: Array, weights: Array | None = None, **kwargs: Any
@@ -135,15 +92,6 @@ def __init__(  # noqa: D417
         self._weights = np.asarray(weights) if weights is not None else None
         self.set_options(**kwargs)
 
-    @property
-    def reality(self) -> bool:
-        """Get the reality of the grid.
-
-        Returns:
-            Reality of the grid.
-        """
-        return True
-
     @staticmethod
     def _resolvent_signs(energies: Array, ordering: Ordering) -> Array:
         """Get the signs for the resolvent based on the time ordering."""
@@ -165,9 +113,8 @@ def resolvent(  # noqa: D417
         chempot: float | Array,
         ordering: Ordering = Ordering.ORDERED,
         invert: bool = True,
-        **kwargs: Any,
     ) -> Array:
-        r"""Get the resolvent of the grid.
+        r"""Get the resolvent of a Lehmann representation on the grid.
 
         For real frequency grids, the resolvent is given by
 
@@ -184,20 +131,25 @@ def resolvent(  # noqa: D417
             invert: Whether to apply the inversion in the resolvent formula.
 
         Returns:
-            Resolvent of the grid.
+            Resolvent on the grid.
         """
-        if kwargs:
-            raise TypeError(f"resolvent() got unexpected keyword argument: {next(iter(kwargs))}")
-        signs = self._resolvent_signs(energies - chempot, ordering)
-        grid = np.expand_dims(self.points, axis=tuple(range(1, energies.ndim + 1)))
+        ndim = energies.ndim
+        if ndim > 2:
+            raise ValueError("energies array must be at most 2D")
+        if ndim < 2:
+            signs = self._resolvent_signs(energies - chempot, ordering)
+        else:
+            signs = self._resolvent_signs(np.diag(energies) - chempot, ordering)
+            signs = np.diag(signs)
+        grid = np.expand_dims(self.points, axis=tuple(range(1, ndim + 1)))
         energies = np.expand_dims(energies, axis=0)
         denominator = grid + (signs * 1.0j * self.eta - energies)
-        return 1.0 / denominator if invert else denominator
+        if invert:
+            return 1.0 / denominator if ndim < 2 else np.linalg.inv(denominator)
+        return denominator
 
     @classmethod
-    def from_uniform(
-        cls, start: float, stop: float, num: int, eta: float | None = None
-    ) -> RealFrequencyGrid:
+    def from_uniform(cls, start: float, stop: float, num: int, eta: float | None = None) -> Self:
         """Create a uniform real frequency grid.
 
         Args:
@@ -212,16 +164,33 @@ def from_uniform(
         points = np.linspace(start, stop, num, endpoint=True)
         return cls(points, eta=eta)
 
+    @classmethod
+    def from_dual(cls, other: RealTimeGrid) -> Self:
+        """Create a grid from another grid in the dual domain (real time).
 
-GridRF = RealFrequencyGrid
+        Args:
+            other: Other (real time) grid to create from.
 
+        Returns:
+            Real frequency grid.
+        """
+        if not other.uniformly_spaced:
+            raise NotImplementedError("only uniformly spaced grids are supported.")
+        if not other.uniformly_weighted:
+            raise NotImplementedError("only uniformly weighted grids are supported.")
+        if other.domain != "time" or not other.reality:
+            raise ValueError(f"dual grid for {cls.__name__} must be of type RealTimeGrid")
+        spacing = 2.0 * np.pi / (other.separation * len(other))
+        num = len(other)
+        points = np.linspace(-spacing * (num - 1) / 2, spacing * (num + 1) / 2, num, endpoint=False)
+        return cls(points, eta=other.eta)
 
-class ImaginaryFrequencyGrid(BaseFrequencyGrid):
-    """Imaginary frequency grid."""
 
-    beta: float = 256
+GridRF = RealFrequencyGrid
 
-    _options = {"beta"}
+
+class ImaginaryFrequencyGrid(BaseFrequencyGrid, BaseImaginaryGrid):
+    """Imaginary frequency grid."""
 
     def __init__(  # noqa: D417
         self, points: Array, weights: Array | None = None, **kwargs: Any
@@ -237,23 +206,14 @@ def __init__(  # noqa: D417
         self._weights = np.asarray(weights) if weights is not None else None
         self.set_options(**kwargs)
 
-    @property
-    def reality(self) -> bool:
-        """Get the reality of the grid.
-
-        Returns:
-            Reality of the grid.
-        """
-        return False
-
     def resolvent(  # noqa: D417
         self,
         energies: Array,
         chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
         invert: bool = True,
-        **kwargs: Any,
     ) -> Array:
-        r"""Get the resolvent of the grid.
+        r"""Get the resolvent of a Lehmann representation on the grid.
 
         For imaginary frequency grids, the resolvent is given by
 
@@ -265,20 +225,24 @@ def resolvent(  # noqa: D417
         Args:
             energies: Energies of the poles.
             chempot: Chemical potential.
+            ordering: Time ordering of the resolvent.
             invert: Whether to apply the inversion in the resolvent formula.
 
         Returns:
-            Resolvent of the grid.
+            Resolvent on the grid.
         """
-        if kwargs:
-            raise TypeError(f"resolvent() got unexpected keyword argument: {next(iter(kwargs))}")
-        grid = np.expand_dims(self.points, axis=tuple(range(1, energies.ndim + 1)))
+        ndim = energies.ndim
+        if ndim > 2:
+            raise ValueError("energies array must be at most 2D")
+        grid = np.expand_dims(self.points, axis=tuple(range(1, ndim + 1)))
         energies = np.expand_dims(energies, axis=0)
         denominator = 1.0j * grid - energies
-        return 1.0 / denominator if invert else denominator
+        if invert:
+            return 1.0 / denominator if ndim < 2 else np.linalg.inv(denominator)
+        return denominator
 
     @classmethod
-    def from_uniform(cls, num: int, beta: float | None = None) -> ImaginaryFrequencyGrid:
+    def from_uniform(cls, num: int, beta: float | None = None) -> Self:
         """Create a uniform imaginary frequency grid.
 
         Args:
@@ -290,16 +254,13 @@ def from_uniform(cls, num: int, beta: float | None = None) -> ImaginaryFrequency
         """
         if beta is None:
             beta = cls.beta
-        separation = 2.0 * np.pi / beta
-        start = 0.5 * separation
-        stop = (num - 0.5) * separation
-        points = np.linspace(start, stop, num, endpoint=True)
+        points = (2 * np.arange(num) + 1) * np.pi / beta
         return cls(points, beta=beta)
 
     @classmethod
     def from_legendre(
         cls, num: int, diffuse_factor: float = 1.0, beta: float | None = None
-    ) -> ImaginaryFrequencyGrid:
+    ) -> Self:
         """Create a Legendre imaginary frequency grid.
 
         Args:
@@ -314,5 +275,23 @@ def from_legendre(
         points = (1 - points) / (diffuse_factor * (1 + points))
         return cls(points, weights=weights, beta=beta)
 
+    @classmethod
+    def from_dual(cls, other: ImaginaryTimeGrid) -> Self:
+        """Create a grid from another grid in the dual domain (imaginary time).
+
+        Args:
+            other: Other (imaginary time) grid to create from.
+
+        Returns:
+            Imaginary frequency grid.
+        """
+        if not other.uniformly_spaced:
+            raise NotImplementedError("only uniformly spaced grids are supported.")
+        if not other.uniformly_weighted:
+            raise NotImplementedError("only uniformly weighted grids are supported.")
+        if other.domain != "time" or other.reality:
+            raise ValueError(f"dual grid for {cls.__name__} must be of type ImaginaryTimeGrid")
+        return cls.from_uniform(len(other) // 2, other.beta)  # Use τ:ω ratio of 2:1
+
 
 GridIF = ImaginaryFrequencyGrid
diff --git a/dyson/grids/grid.py b/dyson/grids/grid.py
index f22cba2..d29515d 100644
--- a/dyson/grids/grid.py
+++ b/dyson/grids/grid.py
@@ -5,8 +5,11 @@
 from abc import ABC, abstractmethod
 from typing import TYPE_CHECKING
 
+from typing_extensions import Self
+
 from dyson import numpy as np
-from dyson.representations.enums import Component, Reduction, RepresentationEnum
+from dyson import util
+from dyson.representations.enums import Component, Ordering, Reduction, RepresentationEnum
 
 if TYPE_CHECKING:
     from typing import Any
@@ -52,22 +55,111 @@ def set_options(self, **kwargs: Any) -> None:
             setattr(self, key, val)
 
     @abstractmethod
+    def _lehmann_kernel(
+        self,
+        energies: Array,
+        chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
+    ) -> Array:
+        """Get the kernel of a Lehmann representation on the grid.
+
+        Args:
+            energies: Energies of the poles.
+            chempot: Chemical potential.
+            ordering: Time ordering of the resolvent or propagator.
+
+        Returns:
+            Kernel of a Lehmann representation on the grid.
+
+        Note:
+            The kernel is a hook to generalise the resolvent or propagator, depending on whether
+            the grid is in the frequency or time domain.
+        """
+        pass
+
     def evaluate_lehmann(
         self,
         lehmann: Lehmann,
         reduction: Reduction = Reduction.NONE,
         component: Component = Component.FULL,
+        ordering: Ordering = Ordering.ORDERED,
     ) -> Dynamic[Any]:
-        """Evaluate a Lehmann representation on the grid.
+        r"""Evaluate a Lehmann representation on the grid.
+
+        The imaginary frequency representation is defined as
+
+        .. math::
+            \sum_{k} \frac{v_{pk} u_{qk}^*}{i \omega - \epsilon_k},
+
+        and the real frequency representation is defined as
+
+        .. math::
+            \sum_{k} \frac{v_{pk} u_{qk}^*}{\omega - \epsilon_k \pm i \eta},
+
+        where :math:`\omega` is the frequency grid, :math:`\epsilon_k` are the poles, and the sign
+        of the broadening factor is determined by the time ordering.
 
         Args:
             lehmann: Lehmann representation to evaluate.
             reduction: The reduction of the dynamic representation.
             component: The component of the dynamic representation.
+            ordering: The time ordering of the Lehmann representation.
 
         Returns:
             Lehmann representation, realised on the grid.
         """
+        from dyson.representations.dynamic import Dynamic  # noqa: PLC0415
+
+        left, right = lehmann.unpack_couplings()
+        kernel = self._lehmann_kernel(lehmann.energies, lehmann.chempot, ordering=ordering)
+        reduction = Reduction(reduction)
+        component = Component(component)
+        ordering = Ordering(ordering)
+
+        # Get the input and output indices based on the reduction type
+        inp = "qk"
+        out = "wpq"
+        if reduction == reduction.NONE:
+            pass
+        elif reduction == reduction.DIAG:
+            inp = "pk"
+            out = "wp"
+        elif reduction == reduction.TRACE:
+            inp = "pk"
+            out = "w"
+        else:
+            reduction.raise_invalid_representation()
+
+        # Perform the downfolding operation
+        array = util.einsum(f"pk,{inp},wk->{out}", right, left.conj(), kernel)
+
+        # Get the required component
+        # TODO: Save time by not evaluating the full array when not needed
+        if component == Component.REAL:
+            array = array.real
+        elif component == Component.IMAG:
+            array = array.imag
+
+        return Dynamic(
+            self,
+            array,
+            reduction=reduction,
+            component=component,
+            ordering=ordering,
+            hermitian=lehmann.hermitian,
+        )
+
+    @classmethod
+    @abstractmethod
+    def from_dual(cls, other: Any) -> Self:
+        """Create a grid from another grid in the dual domain.
+
+        Args:
+            other: Other grid to create from.
+
+        Returns:
+            Grid.
+        """
         pass
 
     @property
@@ -154,3 +246,37 @@ def domain(self) -> str:
     def reality(self) -> bool:
         """Get the reality of the grid."""
         pass
+
+
+class BaseRealGrid(BaseGrid):
+    """Base class for real grids."""
+
+    eta: float = 1e-2
+
+    _options = {"eta"}
+
+    @property
+    def reality(self) -> bool:
+        """Get the reality of the grid.
+
+        Returns:
+            Reality of the grid.
+        """
+        return True
+
+
+class BaseImaginaryGrid(BaseGrid):
+    """Base class for imaginary grids."""
+
+    beta: float = 256
+
+    _options = {"beta"}
+
+    @property
+    def reality(self) -> bool:
+        """Get the reality of the grid.
+
+        Returns:
+            Reality of the grid.
+        """
+        return False
diff --git a/dyson/grids/pade.py b/dyson/grids/pade.py
new file mode 100644
index 0000000..4a791bc
--- /dev/null
+++ b/dyson/grids/pade.py
@@ -0,0 +1,117 @@
+"""Pade method for analytic continuation."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from dyson import numpy as np
+from dyson import util
+from dyson.representations.enums import Component, Ordering, Reduction
+
+if TYPE_CHECKING:
+    from dyson.grids.frequency import BaseFrequencyGrid
+    from dyson.representations.dynamic import Dynamic
+    from dyson.typing import Array
+
+
+def _pade_coefficients(greens_function: Dynamic[BaseFrequencyGrid]) -> Array:
+    """Get the coefficient for the Pade approximation for a frequency domain function.
+
+    Args:
+        greens_function: Dynamic quantity in frequency domain.
+
+    Returns:
+        Coefficients for the Pade approximation.
+    """
+    grid = greens_function.grid
+    if not grid.uniformly_spaced:
+        raise ValueError("only uniform grids are supported.")
+    if not grid.uniformly_weighted:
+        raise ValueError("only uniform weights are supported.")
+
+    # Initialise the coefficients
+    resolvent = grid.resolvent(np.array(0.0), 0.0, invert=False, ordering=greens_function.ordering)
+    coefficients = greens_function.array.astype(np.result_type(greens_function.array, resolvent))
+
+    # Recursively compute the coefficients
+    # FIXME
+    for i in range(len(grid) - 1):
+        with np.errstate(over="ignore"):
+            factor = coefficients[i] * util.reciprocal(coefficients[i + 1 :]) - 1.0
+        factor[~np.isfinite(factor)] = 0.0
+        difference = resolvent[i + 1 :] - resolvent[i]
+        coefficients[i + 1 :] = util.einsum("w...,w->w...", factor, util.reciprocal(difference))
+
+    return coefficients
+
+
+def evaluate_pade(
+    coefficients: Array,
+    grid_old: BaseFrequencyGrid,
+    grid_new: BaseFrequencyGrid,
+    ordering: Ordering = Ordering.RETARDED,
+) -> Array:
+    """Evaluate the Pade approximation on a new frequency grid.
+
+    Args:
+        coefficients: Coefficients for the Pade approximation.
+        grid_old: Original frequency grid.
+        grid_new: New frequency grid.
+        ordering: Ordering of the Green's function.
+
+    Returns:
+        Dynamic quantity on the new frequency grid.
+    """
+    if not grid_old.uniformly_spaced or not grid_new.uniformly_spaced:
+        raise ValueError("only uniform grids are supported.")
+    if not grid_old.uniformly_weighted or not grid_new.uniformly_weighted:
+        raise ValueError("only uniform weights are supported.")
+
+    # Initialise the array
+    resolvent_old = grid_old.resolvent(np.array(0.0), 0.0, invert=False, ordering=ordering)
+    resolvent_new = grid_new.resolvent(np.array(0.0), 0.0, invert=False, ordering=ordering)
+    array = coefficients[[-1]].copy()
+
+    # Recursively evaluate the Pade approximation
+    for i in range(len(grid_old) - 1, -1, -1):
+        term = 1.0 + util.einsum("w,w...->w...", resolvent_new - resolvent_old[i], array)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            array = coefficients[i] / term
+        array[~np.isfinite(array)] = 0.0
+
+    return array
+
+
+def analytic_continuation_freq_pade(
+    greens_function: Dynamic[BaseFrequencyGrid],
+    grid: BaseFrequencyGrid,
+) -> Dynamic[BaseFrequencyGrid]:
+    """Perform analytic continuation in the frequency domain using the Pade approximation.
+
+    Args:
+        greens_function: Green's function in a frequency domain.
+        grid: Frequency grid to continue to.
+
+    Returns:
+        Green's function in the conjugate frequency domain.
+    """
+    if greens_function.ordering == Ordering.ORDERED:
+        raise ValueError("time-ordered ordering is not supported.")
+    if greens_function.component != Component.FULL:
+        raise ValueError("only full component is supported.")
+    if greens_function.reduction == Reduction.TRACE:
+        raise ValueError("traced reduction is not supported.")
+    coefficients = _pade_coefficients(greens_function)
+    array = evaluate_pade(
+        coefficients,
+        greens_function.grid,
+        grid,
+        ordering=greens_function.ordering,
+    )
+    return greens_function.__class__(
+        grid,
+        array,
+        reduction=greens_function.reduction,
+        hermitian=greens_function.hermitian,
+        ordering=greens_function.ordering,
+    )
diff --git a/dyson/grids/time.py b/dyson/grids/time.py
new file mode 100644
index 0000000..08171be
--- /dev/null
+++ b/dyson/grids/time.py
@@ -0,0 +1,250 @@
+"""Time grids."""
+
+from __future__ import annotations
+
+from abc import abstractmethod
+from typing import TYPE_CHECKING
+
+from typing_extensions import Self
+
+from dyson import numpy as np
+from dyson import util
+from dyson.grids.grid import BaseGrid, BaseImaginaryGrid, BaseRealGrid
+from dyson.representations.enums import Ordering
+
+if TYPE_CHECKING:
+    from typing import Any
+
+    from dyson.grids.frequency import ImaginaryFrequencyGrid, RealFrequencyGrid
+    from dyson.typing import Array
+
+
+class BaseTimeGrid(BaseGrid):
+    """Base class for time grids."""
+
+    @property
+    def domain(self) -> str:
+        """Return the domain of the grid."""
+        return "time"
+
+    @abstractmethod
+    def propagator(  # noqa: D417
+        self,
+        energies: Array,
+        chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
+    ) -> Array:
+        """Get the propagator of a Lehmann representation on the grid.
+
+        Args:
+            energies: Energies of the poles.
+            chempot: Chemical potential.
+            ordering: Time ordering of the resolvent.
+
+        Returns:
+            Propagator of the grid.
+        """
+        pass
+
+    def _lehmann_kernel(
+        self,
+        energies: Array,
+        chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
+    ) -> Array:
+        """Get the kernel of a Lehmann representation on the grid.
+
+        Args:
+            energies: Energies of the poles.
+            chempot: Chemical potential.
+            ordering: Time ordering of the propagator.
+
+        Returns:
+            Kernel of a Lehmann representation on the grid.
+
+        Note:
+            The kernel is a hook to generalise the resolvent or propagator, depending on whether
+            the grid is in the frequency or time domain.
+        """
+        return self.propagator(energies, chempot, ordering=ordering)
+
+
+class RealTimeGrid(BaseTimeGrid, BaseRealGrid):
+    """Real time grid."""
+
+    def __init__(  # noqa: D417
+        self, points: Array, weights: Array | None = None, **kwargs: Any
+    ) -> None:
+        """Initialise the grid.
+
+        Args:
+            points: Points of the grid.
+            weights: Weights of the grid.
+            eta: Broadening factor.
+        """
+        self._points = np.asarray(points)
+        self._weights = np.asarray(weights) if weights is not None else None
+        self.set_options(**kwargs)
+
+    @staticmethod
+    def _heaviside(points: Array, energies: Array, ordering: Ordering) -> Array:
+        """Get the Heaviside term."""
+        ordering = Ordering(ordering)
+        theta: Array
+        if ordering == ordering.ORDERED:
+            occ = (energies < 0).astype(np.float64)
+            vir = 1.0 - occ
+            theta = occ * np.heaviside(points, 0.5) - vir * np.heaviside(-points, 0.5)
+        elif ordering == ordering.RETARDED:
+            theta = -np.heaviside(-points, 0.5)
+        elif ordering == ordering.ADVANCED:
+            theta = np.heaviside(points, 0.5)
+        else:
+            ordering.raise_invalid_representation()
+        return theta
+
+    def propagator(  # noqa: D417
+        self,
+        energies: Array,
+        chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
+    ) -> Array:
+        """Get the propagator of a Lehmann representation on the grid.
+
+        Args:
+            energies: Energies of the poles.
+            chempot: Chemical potential.
+            ordering: Time ordering of the resolvent.
+
+        Returns:
+            Propagator on the grid.
+        """
+        ordering = Ordering(ordering)
+        grid = np.expand_dims(self.points, axis=tuple(range(1, energies.ndim + 1)))
+        energies = np.expand_dims(energies, axis=0)
+        if ordering == Ordering.RETARDED:
+            phase = np.exp(grid * self.eta)
+        elif ordering == Ordering.ADVANCED:
+            phase = np.exp(-grid * self.eta)
+        elif ordering == Ordering.ORDERED:
+            phase = np.exp(-np.abs(grid) * self.eta)
+        else:
+            ordering.raise_invalid_representation()
+        theta = self._heaviside(grid, energies - chempot, ordering)
+        propagator = 1.0j * phase * np.exp(1.0j * grid * energies) * theta
+        return propagator
+
+    @classmethod
+    def from_uniform(cls, start: float, stop: float, num: int, eta: float | None = None) -> Self:
+        """Create a uniform real time grid.
+
+        Args:
+            start: Start of the grid.
+            stop: End of the grid.
+            num: Number of points in the grid.
+            eta: Broadening factor.
+
+        Returns:
+            Uniform real time grid.
+        """
+        points = np.linspace(start, stop, num, endpoint=True)
+        return cls(points, eta=eta)
+
+    @classmethod
+    def from_dual(cls, other: RealFrequencyGrid) -> Self:
+        """Create a grid from another grid in the dual domain (real frequency).
+
+        Args:
+            other: Other (real frequency) grid to create from.
+
+        Returns:
+            Real time grid.
+        """
+        if not other.uniformly_spaced:
+            raise NotImplementedError("only uniformly spaced grids are supported.")
+        if not other.uniformly_weighted:
+            raise NotImplementedError("only uniformly weighted grids are supported.")
+        if other.domain != "frequency" or not other.reality:
+            raise ValueError(f"dual grid for {cls.__name__} must be of type RealFrequencyGrid")
+        spacing = 2.0 * np.pi / (other.separation * len(other))
+        num = len(other)
+        points = np.linspace(-spacing * num / 2, spacing * num / 2, num, endpoint=False)
+        return cls(points, eta=other.eta)
+
+
+GridRT = RealTimeGrid
+
+
+class ImaginaryTimeGrid(BaseTimeGrid, BaseImaginaryGrid):
+    """Imaginary time grid."""
+
+    def __init__(  # noqa: D417
+        self, points: Array, weights: Array | None = None, **kwargs: Any
+    ) -> None:
+        """Initialise the grid.
+
+        Args:
+            points: Points of the grid.
+            weights: Weights of the grid.
+        """
+        self._points = np.asarray(points)
+        self._weights = np.asarray(weights) if weights is not None else None
+        self.set_options(**kwargs)
+
+    def propagator(  # noqa: D417
+        self,
+        energies: Array,
+        chempot: float | Array,
+        ordering: Ordering = Ordering.ORDERED,
+    ) -> Array:
+        """Get the propagator of a Lehmann representation on the grid.
+
+        Args:
+            energies: Energies of the poles.
+            chempot: Chemical potential.
+            ordering: Time ordering of the resolvent.
+
+        Returns:
+            Propagator on the grid.
+        """
+        grid = np.expand_dims(self.points, axis=tuple(range(1, energies.ndim + 1)))
+        energies = np.expand_dims(energies - chempot, axis=0)
+        fermi = util.reciprocal(1.0 + util.exp(-self.beta * energies))  # overflow protected
+        propagator = -util.exp(-energies * grid) * fermi  # overflow protected
+        return propagator.astype(np.complex128)
+
+    @classmethod
+    def from_uniform(cls, num: int, beta: float) -> Self:
+        """Create a uniform real time grid.
+
+        Args:
+            num: Number of points in the grid.
+            beta: Inverse temperature.
+
+        Returns:
+            Uniform real time grid.
+        """
+        shift = 0.5 * beta / num
+        points = np.linspace(shift, beta + shift, num, endpoint=True)
+        return cls(points, beta=beta)
+
+    @classmethod
+    def from_dual(cls, other: ImaginaryFrequencyGrid) -> Self:
+        """Create a grid from another grid in the dual domain (imaginary frequency).
+
+        Args:
+            other: Other (imaginary frequency) grid to create from.
+
+        Returns:
+            Imaginary time grid.
+        """
+        if not other.uniformly_spaced:
+            raise NotImplementedError("only uniformly spaced grids are supported.")
+        if not other.uniformly_weighted:
+            raise NotImplementedError("only uniformly weighted grids are supported.")
+        if other.domain != "frequency" or other.reality:
+            raise ValueError(f"dual grid for {cls.__name__} must be of type ImaginaryFrequencyGrid")
+        return cls.from_uniform(len(other) * 2, other.beta)  # Use τ:ω ratio of 2:1
+
+
+GridIT = ImaginaryTimeGrid
diff --git a/dyson/grids/util.py b/dyson/grids/util.py
new file mode 100644
index 0000000..614791e
--- /dev/null
+++ b/dyson/grids/util.py
@@ -0,0 +1,57 @@
+"""Utility functions for grids."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from dyson import numpy as np
+
+if TYPE_CHECKING:
+    from dyson.grids.grid import BaseGrid
+
+
+def are_dual(grid1: BaseGrid, grid2: BaseGrid) -> bool:
+    """Check if a pair of grids are dual to each other.
+
+    Dual grids are related by a Fourier transform.
+
+    Args:
+        grid1: The first grid.
+        grid2: The second grid.
+
+    Returns:
+        Whether the grids are dual to each other.
+    """
+    if not {grid1.domain, grid2.domain} == {"time", "frequency"}:
+        raise ValueError("one grid must be in time domain and the other in frequency domain")
+    if grid1.reality != grid2.reality:
+        raise ValueError("both grids must be either real or imaginary")
+    if not grid1.uniformly_spaced or not grid2.uniformly_spaced:
+        raise ValueError("duality only supported for uniformly spaced grids")
+    if not grid1.uniformly_weighted or not grid2.uniformly_weighted:
+        raise ValueError("duality only supported for uniformly weighted grids")
+    time, freq = (grid1, grid2) if grid1.domain == "time" else (grid2, grid1)
+    time_recov = time.from_dual(freq)
+    freq_recov = freq.from_dual(time)
+    same_points = len(time) == len(time_recov) and np.allclose(time.points, time_recov.points)
+    same_points &= len(freq) == len(freq_recov) and np.allclose(freq.points, freq_recov.points)
+    same_eta = np.isclose(getattr(time, "eta", 0.0), getattr(freq, "eta", 0.0))
+    return bool(same_points and same_eta)
+
+
+def are_equal(grid1: BaseGrid, grid2: BaseGrid) -> bool:
+    """Check if a pair of grids are equal.
+
+    Args:
+        grid1: The first grid.
+        grid2: The second grid.
+
+    Returns:
+        Whether the grids are equal.
+    """
+    same_domain = grid1.domain == grid2.domain
+    same_reality = grid1.reality == grid2.reality
+    same_points = len(grid1) == len(grid2) and np.allclose(grid1.points, grid2.points)
+    same_weights = np.allclose(grid1.weights, grid2.weights)
+    same_eta = np.isclose(getattr(grid1, "eta", 0.0), getattr(grid2, "eta", 0.0))
+    return bool(same_domain and same_reality and same_points and same_weights and same_eta)
diff --git a/dyson/plotting.py b/dyson/plotting.py
index 4b3135e..ed5656e 100644
--- a/dyson/plotting.py
+++ b/dyson/plotting.py
@@ -81,7 +81,7 @@
 plt.rcParams.update(theme)
 
 
-def _unit_name(unit: str) -> str:
+def _unit_name(unit: Literal["Ha", "eV"]) -> str:
     """Return the name of the unit for SciPy."""
     if unit == "Ha":
         return "hartree"
@@ -91,13 +91,21 @@ def _unit_name(unit: str) -> str:
         raise ValueError(f"Unknown energy unit: {unit}. Use 'Ha' or 'eV'.")
 
 
-def _convert(energy: float, unit_from: str, unit_to: str) -> float:
+def _convert(
+    energy: float,
+    unit_from: Literal["Ha", "eV"],
+    unit_to: Literal["Ha", "eV"],
+    domain: Literal["frequency", "time"] = "frequency",
+) -> float:
     """Convert energies between Hartree and eV."""
     if unit_from == unit_to:
         return energy
-    unit_from = _unit_name(unit_from)
-    unit_to = _unit_name(unit_to)
-    return energy * scipy.constants.physical_constants[f"{unit_from}-{unit_to} relationship"][0]
+    if domain == "time":
+        unit_from, unit_to = unit_to, unit_from
+    convert = scipy.constants.physical_constants[
+        f"{_unit_name(unit_from)}-{_unit_name(unit_to)} relationship"
+    ][0]
+    return energy * convert
 
 
 def plot_lehmann(
@@ -172,7 +180,7 @@ def plot_dynamic(
             'real or imaginary part, use dynamic.copy(component="real") or '
             'dynamic.copy(component="imag") to create a copy with the desired component.'
         )
-    grid = _convert(dynamic.grid.points, "Ha", energy_unit)
+    grid = _convert(dynamic.grid.points, "Ha", energy_unit, domain=dynamic.grid.domain)
     array = dynamic.array
     if normalise:
         array = array / np.max(np.abs(array))
diff --git a/dyson/representations/dynamic.py b/dyson/representations/dynamic.py
index 3bf1328..2bf60cc 100644
--- a/dyson/representations/dynamic.py
+++ b/dyson/representations/dynamic.py
@@ -7,7 +7,7 @@
 from dyson import numpy as np
 from dyson import util
 from dyson.grids.grid import BaseGrid
-from dyson.representations.enums import Component, Reduction
+from dyson.representations.enums import Component, Ordering, Reduction
 from dyson.representations.representation import BaseRepresentation
 
 if TYPE_CHECKING:
@@ -32,18 +32,31 @@ def _cast_component(first: Component, second: Component) -> Component:
     return Component.FULL
 
 
+def _cast_ordering(first: Ordering, second: Ordering) -> Ordering:
+    """Find the ordering that is compatible with both orderings."""
+    if first == second:
+        return first
+    raise ValueError(f"Cannot cast orderings {first} and {second}.")
+
+
 def _cast_arrays(first: Dynamic[_TGrid], second: Dynamic[_TGrid]) -> tuple[Array, Array]:
     """Cast the arrays of two dynamic representations to the same component and reduction."""
     component = _cast_component(first.component, second.component)
     reduction = _cast_reduction(first.reduction, second.reduction)
-    array_first = first.as_dynamic(component=component, reduction=reduction).array
-    array_second = second.as_dynamic(component=component, reduction=reduction).array
+    ordering = _cast_ordering(first.ordering, second.ordering)
+    array_first = first.as_dynamic(
+        component=component, reduction=reduction, ordering=ordering
+    ).array
+    array_second = second.as_dynamic(
+        component=component, reduction=reduction, ordering=ordering
+    ).array
     return array_first, array_second
 
 
 def _same_grid(first: Dynamic[_TGrid], second: Dynamic[_TGrid]) -> bool:
     """Check if two dynamic representations have the same grid."""
     # TODO: Move to BaseGrid
+    # TODO: Cache by ID? This can be O(ngrid)
     if not isinstance(second.grid, type(first.grid)):
         return False
     if len(first.grid) != len(second.grid):
@@ -70,6 +83,7 @@ def __init__(
         array: Array,
         reduction: Reduction = Reduction.NONE,
         component: Component = Component.FULL,
+        ordering: Ordering = Ordering.ORDERED,
         hermitian: bool = False,
     ):
         """Initialise the object.
@@ -79,6 +93,7 @@ def __init__(
             array: The array of values at each point in the grid.
             reduction: The reduction of the dynamic representation.
             component: The component of the dynamic representation.
+            ordering: The time ordering of the dynamic representation.
             hermitian: Whether the array is Hermitian.
         """
         self._grid = grid
@@ -86,6 +101,7 @@ def __init__(
         self._hermitian = hermitian
         self._reduction = Reduction(reduction)
         self._component = Component(component)
+        self._ordering = Ordering(ordering)
         if array.shape[0] != len(grid):
             raise ValueError(
                 f"Array must have the same size as the grid in the first dimension, but got "
@@ -109,6 +125,7 @@ def from_lehmann(
         grid: _TGrid,
         reduction: Reduction = Reduction.NONE,
         component: Component = Component.FULL,
+        ordering: Ordering = Ordering.ORDERED,
     ) -> Dynamic[_TGrid]:
         """Construct a dynamic representation from a Lehmann representation.
 
@@ -117,16 +134,19 @@ def from_lehmann(
             grid: The grid on which the dynamic representation is defined.
             reduction: The reduction of the dynamic representation.
             component: The component of the dynamic representation.
+            ordering: The time ordering of the dynamic representation.
 
         Returns:
             A dynamic representation.
         """
-        return grid.evaluate_lehmann(lehmann, reduction=reduction, component=component)
+        return grid.evaluate_lehmann(
+            lehmann, reduction=reduction, component=component, ordering=ordering
+        )
 
     @property
     def nphys(self) -> int:
         """Get the number of physical degrees of freedom."""
-        return self.array.shape[-1]
+        return self.array.shape[-1] if self.reduction != Reduction.TRACE else 1
 
     @property
     def grid(self) -> _TGrid:
@@ -148,6 +168,11 @@ def component(self) -> Component:
         """Get the component of the dynamic representation."""
         return self._component
 
+    @property
+    def ordering(self) -> Ordering:
+        """Get the time ordering of the dynamic representation."""
+        return self._ordering
+
     @property
     def hermitian(self) -> bool:
         """Get a boolean indicating if the system is Hermitian."""
@@ -167,6 +192,7 @@ def copy(
         deep: bool = True,
         reduction: Reduction | None = None,
         component: Component | None = None,
+        ordering: Ordering | None = None,
     ) -> Dynamic[_TGrid]:
         """Return a copy of the dynamic representation.
 
@@ -174,6 +200,7 @@ def copy(
             deep: Whether to return a deep copy of the energies and couplings.
             component: The component of the dynamic representation.
             reduction: The reduction of the dynamic representation.
+            ordering: The time ordering of the dynamic representation.
 
         Returns:
             A new dynamic representation.
@@ -186,6 +213,9 @@ def copy(
         if component is None:
             component = self.component
         component = Component(component)
+        if ordering is None:
+            ordering = self.ordering
+        ordering = Ordering(ordering)
 
         # Copy the array if requested
         if deep:
@@ -225,21 +255,45 @@ def copy(
                     "representation."
                 )
 
-        return self.__class__(grid, array, hermitian=self.hermitian)
+        # Ordering must be the same
+        if ordering != self.ordering:
+            raise ValueError(
+                f"Cannot convert from {self.ordering} to {ordering} for dynamic representation. "
+                "This may be possible for a given dynamic representation, but cannot be "
+                "generalised for all grid types."
+            )
+
+        return self.__class__(
+            grid,
+            array,
+            hermitian=self.hermitian,
+            reduction=reduction,
+            component=component,
+            ordering=ordering,
+        )
 
     def as_dynamic(
-        self, component: Component | None = None, reduction: Reduction | None = None
+        self,
+        component: Component | None = None,
+        reduction: Reduction | None = None,
+        ordering: Ordering | None = None,
     ) -> Dynamic[_TGrid]:
         """Return the dynamic representation with the specified component and reduction.
 
         Args:
             component: The component of the dynamic representation.
             reduction: The reduction of the dynamic representation.
+            ordering: The time ordering of the dynamic representation.
 
         Returns:
             A new dynamic representation with the specified component and reduction.
         """
-        return self.copy(deep=False, component=component, reduction=reduction)
+        return self.copy(
+            deep=False,
+            component=component,
+            reduction=reduction,
+            ordering=ordering,
+        )
 
     def rotate(self, rotation: Array | tuple[Array, Array]) -> Dynamic[_TGrid]:
         """Rotate the dynamic representation.
@@ -273,6 +327,54 @@ def rotate(self, rotation: Array | tuple[Array, Array]) -> Dynamic[_TGrid]:
             array,
             component=component,
             reduction=Reduction.NONE,
+            ordering=self.ordering,
+            hermitian=self.hermitian,
+        )
+
+    def inverse(self) -> Dynamic[_TGrid]:
+        """Return the inverse of the dynamic representation."""
+        if self.reduction == Reduction.NONE:
+            array = np.linalg.inv(self.array)
+        elif self.reduction == Reduction.DIAG:
+            array = 1.0 / self.array
+        elif self.reduction == Reduction.TRACE:
+            raise ValueError("Cannot invert a dynamic representation with trace reduction.")
+        else:
+            self.reduction.raise_invalid_representation()
+        return self.__class__(
+            self.grid,
+            array,
+            component=self.component,
+            reduction=self.reduction,
+            ordering=self.ordering,
+            hermitian=self.hermitian,
+        )
+
+    def conjugate(self) -> Dynamic[_TGrid]:
+        """Return the complex conjugate of the dynamic representation."""
+        return self.__class__(
+            self.grid,
+            np.conjugate(self.array),
+            component=self.component,
+            reduction=self.reduction,
+            ordering=self.ordering,
+            hermitian=self.hermitian,
+        )
+
+    def transpose(self) -> Dynamic[_TGrid]:
+        """Return the transpose of the dynamic representation."""
+        if self.reduction == Reduction.NONE:
+            array = np.transpose(self.array, axes=(0, 2, 1))
+        elif self.reduction in {Reduction.DIAG, Reduction.TRACE}:
+            array = self.array
+        else:
+            self.reduction.raise_invalid_representation()
+        return self.__class__(
+            self.grid,
+            array,
+            component=self.component,
+            reduction=self.reduction,
+            ordering=self.ordering,
             hermitian=self.hermitian,
         )
 
@@ -287,6 +389,7 @@ def __add__(self, other: Dynamic[_TGrid]) -> Dynamic[_TGrid]:
             np.add(*_cast_arrays(self, other)),
             component=_cast_component(self.component, other.component),
             reduction=_cast_reduction(self.reduction, other.reduction),
+            ordering=_cast_ordering(self.ordering, other.ordering),
             hermitian=self.hermitian or other.hermitian,
         )
 
@@ -301,6 +404,7 @@ def __sub__(self, other: Dynamic[_TGrid]) -> Dynamic[_TGrid]:
             np.subtract(*_cast_arrays(self, other)),
             component=_cast_component(self.component, other.component),
             reduction=_cast_reduction(self.reduction, other.reduction),
+            ordering=_cast_ordering(self.ordering, other.ordering),
             hermitian=self.hermitian or other.hermitian,
         )
 
@@ -313,6 +417,7 @@ def __mul__(self, other: float | int) -> Dynamic[_TGrid]:
             self.array * other,
             component=self.component,
             reduction=self.reduction,
+            ordering=self.ordering,
             hermitian=self.hermitian,
         )
 
diff --git a/dyson/representations/enums.py b/dyson/representations/enums.py
index b9dba3a..602dd95 100644
--- a/dyson/representations/enums.py
+++ b/dyson/representations/enums.py
@@ -5,8 +5,10 @@
 from enum import Enum
 from typing import TYPE_CHECKING
 
+from dyson import numpy as np
+
 if TYPE_CHECKING:
-    pass
+    from dyson.typing import Array
 
 
 class RepresentationEnum(Enum):
@@ -38,6 +40,17 @@ def ndim(self) -> int:
         """Get the number of dimensions of the array for this reduction."""
         return {Reduction.NONE: 2, Reduction.DIAG: 1, Reduction.TRACE: 0}[self]
 
+    def identity(self, size: int) -> Array:
+        """Get the identity array for this reduction."""
+        if self == Reduction.NONE:
+            return np.eye(size)
+        elif self == Reduction.DIAG:
+            return np.ones(size)
+        elif self == Reduction.TRACE:
+            return np.array(size)  # \sum_i 1 = size
+        else:
+            self.raise_invalid_representation()
+
 
 class Component(RepresentationEnum):
     """Enumeration for the component of the dynamic representation.
diff --git a/dyson/solvers/__init__.py b/dyson/solvers/__init__.py
index 3b9efb6..f9c5bf9 100644
--- a/dyson/solvers/__init__.py
+++ b/dyson/solvers/__init__.py
@@ -61,6 +61,7 @@
     solver
     static
     dynamic
+    recipes
 
 """
 
@@ -73,3 +74,5 @@
 from dyson.solvers.static.density import DensityRelaxation
 from dyson.solvers.dynamic.corrvec import CorrectionVector
 from dyson.solvers.dynamic.cpgf import CPGF
+from dyson.solvers.dynamic.direct import Direct
+from dyson.solvers.recipes import greens_function_from_hamiltonian
diff --git a/dyson/solvers/dynamic/__init__.py b/dyson/solvers/dynamic/__init__.py
index 68ad18f..d253ba8 100644
--- a/dyson/solvers/dynamic/__init__.py
+++ b/dyson/solvers/dynamic/__init__.py
@@ -8,5 +8,6 @@
 
     corrvec
     cpgf
+    direct
 
 """
diff --git a/dyson/solvers/dynamic/corrvec.py b/dyson/solvers/dynamic/corrvec.py
index 19bd6cb..f5d9179 100644
--- a/dyson/solvers/dynamic/corrvec.py
+++ b/dyson/solvers/dynamic/corrvec.py
@@ -303,6 +303,7 @@ def kernel(self) -> Dynamic[RealFrequencyGrid]:
             greens_function,
             reduction=self.reduction,
             component=self.component,
+            ordering=self.ordering,
             hermitian=self.get_state_ket is None,
         )
 
diff --git a/dyson/solvers/dynamic/cpgf.py b/dyson/solvers/dynamic/cpgf.py
index 321bf29..654a8bd 100644
--- a/dyson/solvers/dynamic/cpgf.py
+++ b/dyson/solvers/dynamic/cpgf.py
@@ -222,6 +222,7 @@ def kernel(self, iteration: int | None = None) -> Dynamic[RealFrequencyGrid]:
             greens_function,
             reduction=self.reduction,
             component=self.component,
+            ordering=self.ordering,
             hermitian=np.allclose(self.moments, self.moments.transpose(0, 2, 1).conj()),
         )
 
diff --git a/dyson/solvers/dynamic/direct.py b/dyson/solvers/dynamic/direct.py
new file mode 100644
index 0000000..d223945
--- /dev/null
+++ b/dyson/solvers/dynamic/direct.py
@@ -0,0 +1,163 @@
+"""Real frequency Dyson equation solver."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from typing_extensions import Self
+
+from dyson import console, printing
+from dyson import numpy as np
+from dyson.grids.util import are_equal
+from dyson.representations.dynamic import Dynamic
+from dyson.representations.enums import Component, Ordering, Reduction
+from dyson.solvers.recipes import greens_function_from_hamiltonian
+from dyson.solvers.solver import DynamicSolver
+
+if TYPE_CHECKING:
+    from typing import Any
+
+    from dyson.expressions.expression import BaseExpression
+    from dyson.grids.frequency import BaseFrequencyGrid
+    from dyson.representations.lehmann import Lehmann
+    from dyson.typing import Array
+
+
+class Direct(DynamicSolver):
+    """Direct frequency-space Dyson equation solver."""
+
+    def __init__(  # noqa: D417
+        self,
+        greens_function_init: Dynamic[BaseFrequencyGrid],
+        self_energy: Dynamic[BaseFrequencyGrid],
+        overlap: Array | None = None,
+        **kwargs: Any,
+    ):
+        """Initialise the solver.
+
+        Args:
+            greens_function_init: The initial Green's function (i.e. :math:`G_0`).
+            self_energy: The self-energy to solve.
+        """
+        self._greens_function_init = greens_function_init
+        self._self_energy = self_energy
+        self._overlap = overlap if overlap is not None else np.eye(greens_function_init.nphys)
+        self.set_options(**kwargs)
+
+    def __post_init__(self) -> None:
+        """Hook called after :meth:`__init__`."""
+        # Check the input
+        if self.greens_function_init.nphys != self.self_energy.nphys:
+            raise ValueError("input functions must have the same number of physical indices")
+        if not are_equal(self.greens_function_init.grid, self.self_energy.grid):
+            raise ValueError("input functions must be defined on the same grid")
+        if self.greens_function_init.grid.domain != "frequency":
+            raise ValueError("input functions must be defined on a frequency grid")
+        if self.greens_function_init.reduction != self.self_energy.reduction:
+            raise ValueError("input functions must have the same reduction")
+        if self.greens_function_init.component != self.self_energy.component:
+            raise ValueError("input functions must have the same component")
+        if self.greens_function_init.ordering != self.self_energy.ordering:
+            raise ValueError("input functions must have the same ordering")
+
+        # Warn if ordering is time-ordered, which may be numerically unstable
+        if self.greens_function_init.ordering == Ordering.ORDERED:
+            console.print(
+                "[bad]Input functions are time-ordered, which may lead to numerical instability "
+                "near poles. Consider using retarded or advanced ordering if possible.[/bad]"
+            )
+
+        # Print the input informationn
+        cond = printing.format_float(
+            np.linalg.cond(self.overlap), threshold=1e10, scientific=True, precision=4
+        )
+        console.print(f"Number of physical states: [input]{self.nphys}[/input]")
+        console.print(f"Overlap condition number: {cond}")
+
+    @classmethod
+    def from_self_energy(
+        cls,
+        static: Array,
+        self_energy: Lehmann,
+        overlap: Array | None = None,
+        **kwargs: Any,
+    ) -> Self:
+        """Create a solver from a self-energy.
+
+        Args:
+            static: Static part of the self-energy.
+            self_energy: Self-energy.
+            overlap: Overlap matrix for the physical space.
+            kwargs: Additional keyword arguments for the solver.
+
+        Returns:
+            Solver instance.
+        """
+        if "grid" not in kwargs:
+            raise ValueError("Missing required argument grid.")
+        grid: BaseFrequencyGrid = kwargs.pop("grid")
+        representation = dict(
+            ordering=kwargs.pop("ordering", Ordering.ORDERED),
+            reduction=kwargs.pop("reduction", Reduction.NONE),
+            component=kwargs.pop("component", Component.FULL),
+        )
+        greens_function_init = grid.evaluate_lehmann(
+            greens_function_from_hamiltonian(static, overlap=overlap),
+            **representation,
+        )
+        self_energy_grid = grid.evaluate_lehmann(self_energy, **representation)
+        return cls(
+            greens_function_init=greens_function_init,
+            self_energy=self_energy_grid,
+            overlap=overlap,
+            **kwargs,
+        )
+
+    @classmethod
+    def from_expression(cls, expression: BaseExpression, **kwargs: Any) -> Self:
+        """Create a solver from an expression.
+
+        Args:
+            expression: Expression to be solved.
+            kwargs: Additional keyword arguments for the solver.
+
+        Returns:
+            Solver instance.
+        """
+        raise NotImplementedError(
+            f"Cannot instantiate {cls.__name__} from an expression. "
+            f"Please use {cls.__name__}.from_self_energy instead."
+        )
+
+    def kernel(self) -> Dynamic[BaseFrequencyGrid]:
+        """Run the solver.
+
+        Returns:
+            The Green's function on the frequency grid.
+        """
+        return (self.greens_function_init.inverse() - self.self_energy).inverse()
+
+    @property
+    def greens_function_init(self) -> Dynamic[BaseFrequencyGrid]:
+        """Get the initial Green's function."""
+        return self._greens_function_init
+
+    @property
+    def self_energy(self) -> Dynamic[BaseFrequencyGrid]:
+        """Get the self-energy."""
+        return self._self_energy
+
+    @property
+    def overlap(self) -> Array:
+        """Get the overlap matrix."""
+        return self._overlap
+
+    @property
+    def grid(self) -> BaseFrequencyGrid:
+        """Get the frequency grid."""
+        return self.greens_function_init.grid
+
+    @property
+    def nphys(self) -> int:
+        """Get the number of physical degrees of freedom."""
+        return self.greens_function_init.nphys
diff --git a/dyson/solvers/recipes.py b/dyson/solvers/recipes.py
new file mode 100644
index 0000000..5305920
--- /dev/null
+++ b/dyson/solvers/recipes.py
@@ -0,0 +1,48 @@
+"""Recipes."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from dyson import numpy as np
+from dyson import util
+from dyson.expressions.hamiltonian import Hamiltonian
+from dyson.solvers.static.exact import Exact
+
+if TYPE_CHECKING:
+    from dyson.representations.lehmann import Lehmann
+    from dyson.typing import Array
+
+
+def greens_function_from_hamiltonian(
+    hamiltonian: Array,
+    overlap: Array | None = None,
+    hermitian: bool | None = None,
+) -> Lehmann:
+    """Construct a Green's function from a Hamiltonian, solving the Hamiltonian exactly.
+
+    This is principally used to find the non-interacting Green's function for a given Hamiltonian,
+    which when expressed on a grid, is given by
+
+    .. math:: G(z) = (z - H)^{-1}.
+
+    Args:
+        hamiltonian: Hamiltonian matrix.
+        overlap: Overlap matrix. If ``None``, the identity matrix is used.
+        hermitian: Whether the Hamiltonian is Hermitian. If ``None``, this is inferred from the
+            Hamiltonian.
+
+    Returns:
+        Green's function.
+    """
+    if hermitian is None:
+        hermitian = np.allclose(hamiltonian, hamiltonian.T.conj())
+    if overlap is None:
+        overlap = bra = ket = np.eye(hamiltonian.shape[-1], dtype=hamiltonian.dtype)
+    else:
+        bra, ket = util.unpack_vectors(util.matrix_power(overlap, 0.5)[0])
+    exp = Hamiltonian(hamiltonian, bra=bra, ket=ket if not hermitian else None)
+    solver = Exact.from_expression(exp)
+    solver.kernel()
+    assert solver.result is not None
+    return solver.result.get_greens_function()
diff --git a/dyson/util/__init__.py b/dyson/util/__init__.py
index 13e96ac..9627239 100644
--- a/dyson/util/__init__.py
+++ b/dyson/util/__init__.py
@@ -34,6 +34,10 @@
     block_diag,
     set_subspace,
     rotate_subspace,
+    upper_bound_for_exp,
+    lower_bound_for_reciprocal,
+    exp,
+    reciprocal,
 )
 from dyson.util.moments import (
     se_moments_to_gf_moments,
diff --git a/dyson/util/linalg.py b/dyson/util/linalg.py
index 24d80f7..1469cec 100644
--- a/dyson/util/linalg.py
+++ b/dyson/util/linalg.py
@@ -626,3 +626,54 @@ def rotate_subspace(vectors: Array, rotation: Array) -> Array:
     size = rotation.shape[0]
     subspace = rotation @ vectors[:size]
     return set_subspace(vectors, subspace)
+
+
+def upper_bound_for_exp(x: Array, factor: float = 1.0) -> Array:
+    r"""Upper bound for the exponential function to avoid overflow.
+
+    For 64-bit floats, this is equivalent to ensuring that
+
+    .. math
+        \exp(x \times \mathrm{factor}) \leq \exp(709.78).
+
+    Args:
+        x: Input array.
+        factor: Scaling factor for the input array.
+
+    Returns:
+        The upper-bounded array.
+    """
+    dtype_max = np.finfo(x.dtype).max
+    threshold = np.log(dtype_max) / factor
+    return np.minimum(x, threshold)
+
+
+def exp(x: Array) -> Array:
+    """Exponentiate an array with overflow protection."""
+    return np.exp(upper_bound_for_exp(x))
+
+
+def lower_bound_for_reciprocal(x: Array, factor: float = 1.0) -> Array:
+    r"""Lower bound for the reciprocal function to avoid overflow.
+
+    For 64-bit floats, this is equivalent to ensuring that
+
+    .. math
+        \vert x \times \mathrm{factor} \vert \geq 5.56 \times 10^{-309}.
+
+    Args:
+        x: Input array.
+        factor: Scaling factor for the input array.
+
+    Returns:
+        The lower-bounded array.
+    """
+    dtype_min = np.finfo(x.dtype).tiny
+    threshold = dtype_min / factor
+    signs = np.where(x >= 0, 1, -1)
+    return np.where(np.abs(x) < threshold, signs * threshold, x)
+
+
+def reciprocal(x: Array) -> Array:
+    """Return the reciprocal of an array with overflow protection."""
+    return 1.0 / lower_bound_for_reciprocal(x)
diff --git a/examples/solver-direct.py b/examples/solver-direct.py
new file mode 100644
index 0000000..4045b6d
--- /dev/null
+++ b/examples/solver-direct.py
@@ -0,0 +1,64 @@
+"""Example of the direct frequency-space solver.
+
+This solver is the traditional form of the frequency-dependent Dyson equation, inverting the
+downfolded Green's function matrix at each frequency point.
+"""
+
+import numpy
+from pyscf import gto, scf
+
+from dyson import FCI, Direct, Exact
+from dyson.grids import RealFrequencyGrid
+from dyson.representations.spectral import Spectral
+from dyson.solvers.recipes import greens_function_from_hamiltonian
+
+# Get a molecule and mean-field from PySCF
+mol = gto.M(atom="Li 0 0 0; H 0 0 1.64", basis="sto-3g", verbose=0)
+mf = scf.RHF(mol)
+mf.kernel()
+
+# Use an FCI expression for the Hamiltonian
+exp_h = FCI.hole.from_mf(mf)
+exp_p = FCI.particle.from_mf(mf)
+
+# Initialise a real frequency grid for the direct solver
+grid = RealFrequencyGrid.from_uniform(-5.0, 5.0, 512, eta=1e-2)
+
+# Use the exact solver to get the central self-energy for demonstration purposes
+exact_h = Exact.from_expression(exp_h)
+exact_h.kernel()
+exact_p = Exact.from_expression(exp_p)
+exact_p.kernel()
+assert exact_h.result is not None and exact_p.result is not None
+result = Spectral.combine_for_self_energy(exact_h.result, exact_p.result)
+static = result.get_static_self_energy()
+self_energy = result.get_self_energy()
+overlap = result.get_overlap()
+
+# Solve the Hamiltonian using the Direct solver, initialisation via either:
+
+# 1) Create the solver from a self-energy.
+#
+#    Note that the default time-ordering (ordered) will display a warning as it is numerically
+#    unstable, but the default is kept consistent with other functionality in the package).
+solver = Direct.from_self_energy(
+    static, self_energy, overlap=overlap, grid=grid, ordering="advanced"
+)
+solver.kernel()
+
+# 2) Create the solver directly from the initial Green's function and self-energy
+gf_0 = grid.evaluate_lehmann(
+    greens_function_from_hamiltonian(static, overlap=overlap),
+    ordering="advanced",
+)
+se = grid.evaluate_lehmann(self_energy, ordering="advanced")
+solver = Direct(gf_0, se, overlap=overlap)
+solver.kernel()
+
+# 3) Don't bother using the solver at all because this is trivial!
+gf = (gf_0.inverse() - se).inverse()
+
+# Compare to that of the Exact solver, by downfolding the Green's function corresponding to the
+# exact result onto the same grid
+gf_exact = grid.evaluate_lehmann(result.get_greens_function(), ordering="advanced")
+print("Direct solver error:", numpy.max(numpy.abs(gf - gf_exact)))
diff --git a/tests/test_cpgf.py b/tests/test_cpgf.py
index d42674c..8cc7c23 100644
--- a/tests/test_cpgf.py
+++ b/tests/test_cpgf.py
@@ -9,6 +9,7 @@
 
 from dyson.expressions.hf import BaseHF
 from dyson.grids import RealFrequencyGrid
+from dyson.representations.enums import Ordering
 from dyson.solvers import CPGF
 
 if TYPE_CHECKING:
@@ -36,7 +37,7 @@ def test_vs_exact_solver(
     # Solve the Hamiltonian exactly
     exact = exact_cache(mf, expression_cls)
     assert exact.result is not None
-    gf_exact = grid.evaluate_lehmann(exact.result.get_greens_function(), ordering="advanced")
+    gf_exact = grid.evaluate_lehmann(exact.result.get_greens_function(), ordering=Ordering.ADVANCED)
 
     # Solve the Hamiltonian with CorrectionVector
     cpgf = CPGF.from_self_energy(
diff --git a/tests/test_direct.py b/tests/test_direct.py
new file mode 100644
index 0000000..aae68ee
--- /dev/null
+++ b/tests/test_direct.py
@@ -0,0 +1,76 @@
+"""Tests for :module:`~dyson.solvers.dynamic.direct`."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pytest
+
+from dyson.grids import ImaginaryFrequencyGrid, RealFrequencyGrid
+from dyson.representations.enums import Ordering
+from dyson.representations.spectral import Spectral
+from dyson.solvers import Direct, Exact
+from dyson.solvers.recipes import greens_function_from_hamiltonian
+
+if TYPE_CHECKING:
+    from pyscf import scf
+
+    from dyson.grids.grid import BaseGrid
+
+    from .conftest import ExactGetter, ExpressionCollection, Helper
+
+
+@pytest.mark.parametrize(
+    "grid",
+    [
+        RealFrequencyGrid.from_uniform(-5, 5, 1001, eta=0.2),
+        ImaginaryFrequencyGrid.from_uniform(501, beta=32),
+    ],
+)
+@pytest.mark.parametrize("ordering", [Ordering.ADVANCED, Ordering.RETARDED])
+def test_vs_exact_solver(
+    helper: Helper,
+    mf: scf.hf.RHF,
+    expression_method: ExpressionCollection,
+    exact_cache: ExactGetter,
+    grid: BaseGrid,
+    ordering: Ordering,
+) -> None:
+    """Test frequency-space solver compared to the exact solver."""
+    # Get the quantities required from the expressions
+    expression_h = expression_method.h.from_mf(mf)
+    expression_p = expression_method.p.from_mf(mf)
+    if expression_h.nconfig > 1024 or expression_p.nconfig > 1024:
+        pytest.skip("Skipping test for large Hamiltonian")
+
+    # Solve the Hamiltonian exactly
+    exact_h = exact_cache(mf, expression_method.h)
+    exact_p = exact_cache(mf, expression_method.p)
+    assert exact_h.result is not None
+    assert exact_p.result is not None
+    result = Spectral.combine_for_self_energy(exact_h.result, exact_p.result)
+
+    # Get the central self-energy
+    overlap = result.get_overlap()
+    static = result.get_static_self_energy()
+    self_energy = result.get_self_energy()
+
+    # Solve the self-energy exactly to get the exact Green's function
+    exact = Exact.from_self_energy(static, self_energy, overlap=overlap)
+    exact.kernel()
+    assert exact.result is not None
+    se_exact = grid.evaluate_lehmann(exact.result.get_self_energy(), ordering=ordering)
+    gf_exact = grid.evaluate_lehmann(exact.result.get_greens_function(), ordering=ordering)
+
+    # Get G_0
+    gf_0 = grid.evaluate_lehmann(
+        greens_function_from_hamiltonian(static, overlap=overlap, hermitian=result.hermitian),
+        ordering=ordering,
+    )
+
+    # Solve the Hamiltonian with RealFrequency
+    direct = Direct(gf_0, se_exact, overlap=overlap)
+    gf = direct.kernel()
+
+    assert np.allclose(gf.array, gf_exact.array)
diff --git a/tests/test_grids.py b/tests/test_grids.py
new file mode 100644
index 0000000..6cccf37
--- /dev/null
+++ b/tests/test_grids.py
@@ -0,0 +1,152 @@
+"""Tests for :module:`~dyson.grids`."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import numpy as np
+import pytest
+
+from dyson.grids import GridIF, GridIT, GridRF, GridRT, transform
+from dyson.grids.grid import BaseGrid
+from dyson.representations.enums import Component, Ordering, Reduction
+from dyson.representations.spectral import Spectral
+
+if TYPE_CHECKING:
+    from pyscf import scf
+
+    from dyson.expressions.expression import ExpressionCollection
+    from dyson.representations.dynamic import Dynamic
+
+    from .conftest import ExactGetter, Helper
+
+
+def with_dual(grid: BaseGrid) -> tuple[BaseGrid, BaseGrid]:
+    """Return a grid and its dual."""
+    dual_class = {
+        GridIF: GridIT,
+        GridIT: GridIF,
+        GridRF: GridRT,
+        GridRT: GridRF,
+    }[type(grid)]
+    return grid, dual_class.from_dual(grid)
+
+
+def get_dynamic_pair(
+    helper: Helper,
+    mf: scf.hf.RHF,
+    expression_method: ExpressionCollection,
+    exact_cache: ExactGetter,
+    ordering: Ordering,
+    reduction: Reduction,
+    component: Component,
+    grid1: BaseGrid,
+    grid2: BaseGrid,
+) -> tuple[Dynamic, Dynamic]:
+    """Get the pair of dynamic Green's functions for the given parameters."""
+    expression_h = expression_method.h.from_mf(mf)
+    expression_p = expression_method.p.from_mf(mf)
+    if expression_h.nconfig > 1024 or expression_p.nconfig > 1024:
+        pytest.skip("Skipping test for large Hamiltonian")
+
+    # Solve the Hamiltonian exactly
+    exact_h = exact_cache(mf, expression_method.h)
+    exact_p = exact_cache(mf, expression_method.p)
+    assert exact_h.result is not None
+    assert exact_p.result is not None
+    result = Spectral.combine_for_greens_function(exact_h.result, exact_p.result)
+    gf1 = grid1.evaluate_lehmann(
+        result.get_greens_function(),
+        ordering=ordering,
+        reduction=reduction,
+        component=component,
+    )
+    gf2 = grid2.evaluate_lehmann(
+        result.get_greens_function(),
+        ordering=ordering,
+        reduction=reduction,
+        component=component,
+    )
+
+    return gf1, gf2
+
+
+@pytest.mark.parametrize("ordering", [Ordering.RETARDED, Ordering.ADVANCED])
+@pytest.mark.parametrize("reduction", [Reduction.NONE, Reduction.DIAG])
+@pytest.mark.parametrize("component", [Component.FULL])
+@pytest.mark.parametrize("grid_rf", [GridRF.from_uniform(-8, 8, 501, eta=0.1)])
+@pytest.mark.parametrize("grid_if", [GridIF.from_uniform(501, beta=32)])
+def test_transform_rf_if(
+    helper: Helper,
+    mf: scf.hf.RHF,
+    expression_method: ExpressionCollection,
+    exact_cache: ExactGetter,
+    ordering: Ordering,
+    reduction: Reduction,
+    component: Component,
+    grid_rf: BaseGrid,
+    grid_if: BaseGrid,
+) -> None:
+    """Test Fourier transform between imaginary time and frequency."""
+    gf_rf, gf_if = get_dynamic_pair(
+        helper,
+        mf,
+        expression_method,
+        exact_cache,
+        ordering,
+        reduction,
+        component,
+        grid_rf,
+        grid_if,
+    )
+
+    # Transform the Green's functions
+    gf_if_recov = transform(gf_rf, grid_if)
+    gf_rf_recov = transform(gf_if, grid_rf)
+
+    # Check that the recovered Green's functions match the original ones
+    assert np.mean(np.abs(gf_rf_recov - gf_rf)) < 1e-1  # approximate
+    assert np.allclose(gf_if_recov, gf_if)  # exact
+
+
+@pytest.mark.parametrize("ordering", [Ordering.RETARDED, Ordering.ADVANCED, Ordering.ORDERED])
+@pytest.mark.parametrize("reduction", [Reduction.NONE, Reduction.DIAG])
+@pytest.mark.parametrize("component", [Component.FULL])
+@pytest.mark.parametrize("grid_if, grid_it", [with_dual(GridIF.from_uniform(501, beta=32))])
+def test_transform_if_it(
+    helper: Helper,
+    mf: scf.hf.RHF,
+    expression_method: ExpressionCollection,
+    exact_cache: ExactGetter,
+    ordering: Ordering,
+    reduction: Reduction,
+    component: Component,
+    grid_if: BaseGrid,
+    grid_it: BaseGrid,
+) -> None:
+    """Test Fourier transform between imaginary time and frequency."""
+    gf_if, gf_it = get_dynamic_pair(
+        helper,
+        mf,
+        expression_method,
+        exact_cache,
+        ordering,
+        reduction,
+        component,
+        grid_if,
+        grid_it,
+    )
+
+    # Transform the Green's functions
+    gf_it_recov = transform(gf_if, grid_it)
+    gf_if_recov = transform(gf_it, grid_if)
+
+    # Check that the recovered Green's functions match the original ones
+    # TODO Fix edges and precision with high frequency tail
+    if "O" in mf.mol.atom and expression_method._name in ["FCI", "ADC(2)"]:
+        pytest.skip("Skipping test for h2o-sto3g due to precision issues with high frequency tail")
+    edges = len(grid_it) // 16
+    mask = np.ones(len(grid_it), dtype=bool)
+    mask[:edges] = mask[-edges:] = False
+    assert np.allclose(gf_if_recov.array, gf_if.array, atol=0.05)
+    assert np.allclose(gf_it_recov.array[mask], gf_it.array[mask], atol=0.05)