diff --git a/src/AbstractSD.jl b/src/AbstractSD.jl
index ffa0887..a6b62e2 100644
--- a/src/AbstractSD.jl
+++ b/src/AbstractSD.jl
@@ -82,4 +82,40 @@ function correlations(J::AbstractSD, τ, β)
     IRe = quadgk(IntRe, 0.0, Inf)[1]
     IIm = quadgk(IntIm, 0.0, Inf)[1]
     return IRe + 1im*IIm
-end
\ No newline at end of file
+end
+
+"""
+    memory_kernel(J::AbstractSD, τ)
+
+Calculate the memory kernel for a spectral density `J` at a given time delay `τ`,
+that is
+```math
+\\mathcal{K}(\\tau) = 2\\Theta(\\tau)\\int_0^\\infty J(\\omega)\\sin(\\omega)\\mathrm{d}\\omega,
+```
+where ``\\Theta`` is the Heavisde theta function.
+
+# Arguments
+- `J::AbstractSD`: The spectral density.
+- `τ`: The time delay at which the memory kernel is evaluated.
+
+# Returns
+- The memory kernel for the spectral density `J` at the given time delay `τ`.
+
+"""
+memory_kernel(J::AbstractSD, τ) = τ <= zero(τ) ? zero(τ) : quadgk(ω -> 2*J(ω)*sin(ω*τ), zero(τ), Inf)[1]
+
+"""
+    imag_memory_kernel_ft(J::AbstractSD, ω)
+
+Calculate the imaginary part of the Fourier-transform of the  memory kernel for a
+spectral density `J` at a given frequency `ω`.
+
+# Arguments
+- `J::AbstractSD`: The spectral density.
+- `ω`: The frequency at which the imaginary part of the Fourier-transform of the memory kernel is evaluated.
+
+# Returns
+- The imaginary part of the Fourier-transform of the memory kernel for the spectral density `J` at the given frequency `ω`.
+
+"""
+imag_memory_kernel_ft(J::AbstractSD, ω) = π*J(ω)
\ No newline at end of file
diff --git a/src/LorentzianSD.jl b/src/LorentzianSD.jl
index 9cc3742..bba8f01 100644
--- a/src/LorentzianSD.jl
+++ b/src/LorentzianSD.jl
@@ -41,6 +41,16 @@ sdoverω(J::LorentzianSD, ω) = (J.α*J.Γ/π)/((ω^2 - J.ω0^2)^2 + (J.Γ*ω)^2
 
 reorganisation_energy(J::LorentzianSD) = J.α/J.ω0^2/2
 
+function memory_kernel(J::LorentzianSD, τ)
+    if τ <= zero(τ)
+        return zero(τ)
+    end
+    ω1 = sqrt(J.ω0^2 - J.Γ^2/4)
+    return J.α*exp(-J.Γ*τ/2)*sin(ω1*τ)/ω1
+end
+
+imag_memory_kernel_ft(J::LorentzianSD, ω) = J.α/(J.ω0^2 - ω^2 - 1im*ω*J.Γ)
+
 """
     struct UnderdampedSD <: AbstractSD