diff --git a/Project.toml b/Project.toml
index 9076f4fd6b8..1698206412d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -70,14 +70,14 @@ TrixiSparseConnectivityTracerExt = "SparseConnectivityTracer"
 
 [compat]
 Accessors = "0.1.36"
-Adapt = "4"
-CUDA = "5.8"
+Adapt = "4.1"
+CUDA = "5.8.2"
 CodeTracking = "1.0.5, 2, 3"
 ConstructionBase = "1.5"
 Convex = "0.16"
 DataStructures = "0.18.15, 0.19"
 DelimitedFiles = "1"
-DiffEqBase = "6.155.2"
+DiffEqBase = "6.174"
 DiffEqCallbacks = "2.35, 3, 4"
 Downloads = "1.6"
 ECOS = "1.1.2"
@@ -89,7 +89,7 @@ KernelAbstractions = "0.9.36"
 LinearAlgebra = "1"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.171"
-MPI = "0.20.6"
+MPI = "0.20.22"
 Makie = "0.21, 0.22, 0.23, 0.24"
 MuladdMacro = "0.2.4"
 NLsolve = "4.5.1"
@@ -104,7 +104,7 @@ RecipesBase = "1.3.4"
 RecursiveArrayTools = "3.31.1"
 Reexport = "1.2"
 Requires = "1.3"
-SciMLBase = "2.67.0"
+SciMLBase = "2.92.0"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 SparseConnectivityTracer = "1.0.1"
diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_imex_operator.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_imex_operator.jl
new file mode 100644
index 00000000000..bfa1e0073a2
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_imex_operator.jl
@@ -0,0 +1,74 @@
+using OrdinaryDiffEqSDIRK
+using SciMLBase, SciMLOperators, SparseArrays
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+advection_velocity = (1.5, 1.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+diffusivity() = 0.5f0
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                periodicity = true,
+                n_cells_max = 30_000)
+
+function initial_condition_diffusive_convergence_test(x, t,
+                                                      equation::LinearScalarAdvectionEquation2D)
+    # Store translated coordinate for easy use of exact solution
+    RealT = eltype(x)
+    x_trans = x - equation.advection_velocity * t
+
+    nu = diffusivity()
+    c = 1
+    A = 0.5f0
+    L = 2
+    f = 1.0f0 / L
+    omega = 2 * convert(RealT, pi) * f
+    scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t)
+    return SVector(scalar)
+end
+initial_condition = initial_condition_diffusive_convergence_test
+
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE setup
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan) # For accessing the initial condition u0 only
+
+D_map, _ = linear_structure_parabolic(semi)
+# Cannot directly construct `MatrixOperator` from `LinearMap`, need detour via sparse matrix
+D_op = MatrixOperator(sparse(D_map))
+
+split_func = SplitFunction(D_op, Trixi.rhs!)
+ode_operator = SplitODEProblem{true}(split_func, ode.u0, tspan, semi)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode_operator, KenCarp4();
+            adaptive = false, dt = 1e-2,
+            ode_default_options()..., callback = callbacks)
diff --git a/src/Trixi.jl b/src/Trixi.jl
index f22e4051a33..2c267a745ff 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -253,6 +253,7 @@ export entropy, energy_total, energy_kinetic, energy_internal,
        enstrophy, vorticity
 export lake_at_rest_error
 export ncomponents, eachcomponent
+export have_constant_speed
 
 export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh,
        P4estMeshView, P4estMeshCubedSphere, T8codeMesh
@@ -321,7 +322,7 @@ export ode_norm, ode_unstable_check
 
 export convergence_test,
        jacobian_fd, jacobian_ad_forward, jacobian_ad_forward_parabolic,
-       linear_structure
+       linear_structure, linear_structure_parabolic
 
 export DGMulti, DGMultiBasis, estimate_dt, DGMultiMesh, GaussSBP
 
diff --git a/src/equations/acoustic_perturbation_2d.jl b/src/equations/acoustic_perturbation_2d.jl
index 00763a11904..26986a813ed 100644
--- a/src/equations/acoustic_perturbation_2d.jl
+++ b/src/equations/acoustic_perturbation_2d.jl
@@ -399,6 +399,24 @@ end
     return SVector(diss[1], diss[2], diss[3], z, z, z, z)
 end
 
+"""
+    have_constant_speed(::AcousticPerturbationEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref).
+
+The acoustic perturbation equations are in principle linear for **constant** mean flow fields.
+However, the mean flow variables are part of the solution vector in
+[`AcousticPerturbationEquations2D`](@ref) and only the **global** mean flow variables are constant,
+similar to the [`LinearizedEulerEquations2D`](@ref).
+
+Moreover, when coupling to the [`CompressibleEulerEquations2D`](@ref) equations via 
+[`SemidiscretizationEulerAcoustics`](@ref), the mean field variables are updated
+on the fly, see [`EulerAcousticsCouplingCallback`](@ref).
+
+# Returns
+- `False()`
+"""
 @inline have_constant_speed(::AcousticPerturbationEquations2D) = False()
 
 @inline function max_abs_speeds(u, equations::AcousticPerturbationEquations2D)
diff --git a/src/equations/compressible_navier_stokes.jl b/src/equations/compressible_navier_stokes.jl
index 7320c5493b5..3cfeaaf5018 100644
--- a/src/equations/compressible_navier_stokes.jl
+++ b/src/equations/compressible_navier_stokes.jl
@@ -110,7 +110,7 @@ dynamic_viscosity(u, mu::Real, equations) = mu
 dynamic_viscosity(u, mu::T, equations) where {T} = mu(u, equations)
 
 """
-    max_diffusivity(::AbstractCompressibleNavierStokesDiffusion)
+    have_constant_diffusivity(::AbstractCompressibleNavierStokesDiffusion)
 
 # Returns
 - `False()`
@@ -118,5 +118,8 @@ dynamic_viscosity(u, mu::T, equations) where {T} = mu(u, equations)
 Used in diffusive CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the
 diffusivity is not constant in space and that [`max_diffusivity`](@ref) needs to be computed
 at every node in every element.
+
+Also employed in [`linear_structure`](@ref) and [`linear_structure_parabolic`](@ref) to check
+if the diffusion term is linear in the variables/constant.
 """
 @inline have_constant_diffusivity(::AbstractCompressibleNavierStokesDiffusion) = False()
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 9075022dd13..87e9744dd7d 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -362,6 +362,17 @@ See also the test section P4estMesh2D with combine_conservative_and_nonconservat
 """
 combine_conservative_and_nonconservative_fluxes(flux, ::AbstractEquations) = False()
 
+"""
+    have_constant_speed(::AbstractEquations)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref).
+
+This is the default fallback for nonlinear equations.
+
+# Returns
+- `False()`
+"""
 have_constant_speed(::AbstractEquations) = False()
 
 """
diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl
index e00e43daeed..bfd33cfae78 100644
--- a/src/equations/equations_parabolic.jl
+++ b/src/equations/equations_parabolic.jl
@@ -12,7 +12,7 @@ abstract type AbstractLaplaceDiffusion{NDIMS, NVARS} <:
               AbstractEquationsParabolic{NDIMS, NVARS, GradientVariablesConservative} end
 
 """
-    max_diffusivity(::AbstractLaplaceDiffusion)
+    have_constant_diffusivity(::AbstractLaplaceDiffusion)
 
 # Returns
 - `True()`
@@ -20,6 +20,9 @@ abstract type AbstractLaplaceDiffusion{NDIMS, NVARS} <:
 Used in diffusive CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the
 diffusivity is constant in space and that [`max_diffusivity`](@ref) needs **not** to be re-computed
 at every node in every element.
+
+Also employed in [`linear_structure`](@ref) and [`linear_structure_parabolic`](@ref) to check
+if the diffusion term is linear in the variables/constant.
 """
 @inline have_constant_diffusivity(::AbstractLaplaceDiffusion) = True()
 
diff --git a/src/equations/hyperbolic_diffusion_1d.jl b/src/equations/hyperbolic_diffusion_1d.jl
index 48601dfd675..3cf82b9c95c 100644
--- a/src/equations/hyperbolic_diffusion_1d.jl
+++ b/src/equations/hyperbolic_diffusion_1d.jl
@@ -172,6 +172,15 @@ end
     return sqrt(equations.nu * equations.inv_Tr)
 end
 
+"""
+    have_constant_speed(::HyperbolicDiffusionEquations1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::HyperbolicDiffusionEquations1D) = True()
 
 @inline function max_abs_speeds(eq::HyperbolicDiffusionEquations1D)
diff --git a/src/equations/hyperbolic_diffusion_2d.jl b/src/equations/hyperbolic_diffusion_2d.jl
index a4d52ea6265..63b1856771f 100644
--- a/src/equations/hyperbolic_diffusion_2d.jl
+++ b/src/equations/hyperbolic_diffusion_2d.jl
@@ -224,6 +224,15 @@ end
     return SVector(f1, f2, f3)
 end
 
+"""
+    have_constant_speed(::HyperbolicDiffusionEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::HyperbolicDiffusionEquations2D) = True()
 
 @inline function max_abs_speeds(eq::HyperbolicDiffusionEquations2D)
diff --git a/src/equations/hyperbolic_diffusion_3d.jl b/src/equations/hyperbolic_diffusion_3d.jl
index 12dfa393185..6b917263cd1 100644
--- a/src/equations/hyperbolic_diffusion_3d.jl
+++ b/src/equations/hyperbolic_diffusion_3d.jl
@@ -223,6 +223,15 @@ Godunov (upwind) flux for the 3D hyperbolic diffusion equations.
     return SVector(f1, f2, f3, f4)
 end
 
+"""
+    have_constant_speed(::HyperbolicDiffusionEquations3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::HyperbolicDiffusionEquations3D) = True()
 
 @inline function max_abs_speeds(eq::HyperbolicDiffusionEquations3D)
diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl
index 4e359b84e3a..9c6d22772df 100644
--- a/src/equations/lattice_boltzmann_2d.jl
+++ b/src/equations/lattice_boltzmann_2d.jl
@@ -371,6 +371,15 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
     return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
 end
 
+"""
+    have_constant_speed(::LatticeBoltzmannEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LatticeBoltzmannEquations2D) = True()
 
 @inline function max_abs_speeds(equations::LatticeBoltzmannEquations2D)
diff --git a/src/equations/lattice_boltzmann_3d.jl b/src/equations/lattice_boltzmann_3d.jl
index 5e5cb7e7b59..7dfabfeddb2 100644
--- a/src/equations/lattice_boltzmann_3d.jl
+++ b/src/equations/lattice_boltzmann_3d.jl
@@ -384,6 +384,15 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
     return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
 end
 
+"""
+    have_constant_speed(::LatticeBoltzmannEquations3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LatticeBoltzmannEquations3D) = True()
 
 @inline function max_abs_speeds(equations::LatticeBoltzmannEquations3D)
diff --git a/src/equations/linear_elasticity_1d.jl b/src/equations/linear_elasticity_1d.jl
index a7d2e782ee3..33b7ed2a71b 100644
--- a/src/equations/linear_elasticity_1d.jl
+++ b/src/equations/linear_elasticity_1d.jl
@@ -97,6 +97,15 @@ end
     return SVector(f1, f2)
 end
 
+"""
+    have_constant_speed(::LinearElasticityEquations1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearElasticityEquations1D) = True()
 
 @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl
index 4692f56af99..9297d9758d9 100644
--- a/src/equations/linear_scalar_advection_1d.jl
+++ b/src/equations/linear_scalar_advection_1d.jl
@@ -170,6 +170,15 @@ function flux_engquist_osher(u_ll, u_rr, orientation::Int,
                     abs(equation.advection_velocity[orientation]) * (u_R - u_L)))
 end
 
+"""
+    have_constant_speed(::LinearScalarAdvectionEquation1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearScalarAdvectionEquation1D) = True()
 
 @inline function max_abs_speeds(equation::LinearScalarAdvectionEquation1D)
diff --git a/src/equations/linear_scalar_advection_2d.jl b/src/equations/linear_scalar_advection_2d.jl
index b90f24c1947..07a16d67e64 100644
--- a/src/equations/linear_scalar_advection_2d.jl
+++ b/src/equations/linear_scalar_advection_2d.jl
@@ -269,6 +269,15 @@ function flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
     end
 end
 
+"""
+    have_constant_speed(::LinearScalarAdvectionEquation2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearScalarAdvectionEquation2D) = True()
 
 @inline function max_abs_speeds(equation::LinearScalarAdvectionEquation2D)
diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl
index 61f035ea703..70a713e7796 100644
--- a/src/equations/linear_scalar_advection_3d.jl
+++ b/src/equations/linear_scalar_advection_3d.jl
@@ -188,6 +188,15 @@ function flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
     end
 end
 
+"""
+    have_constant_speed(::LinearScalarAdvectionEquation3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearScalarAdvectionEquation3D) = True()
 
 @inline function max_abs_speeds(equation::LinearScalarAdvectionEquation3D)
diff --git a/src/equations/linearized_euler_1d.jl b/src/equations/linearized_euler_1d.jl
index 19a3bdcb3bd..3328cabc9f9 100644
--- a/src/equations/linearized_euler_1d.jl
+++ b/src/equations/linearized_euler_1d.jl
@@ -108,6 +108,15 @@ end
     return SVector(f1, f2, f3)
 end
 
+"""
+    have_constant_speed(::LinearizedEulerEquations1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearizedEulerEquations1D) = True()
 
 @inline function max_abs_speeds(equations::LinearizedEulerEquations1D)
diff --git a/src/equations/linearized_euler_2d.jl b/src/equations/linearized_euler_2d.jl
index 985fd04e604..9584c40f4fe 100644
--- a/src/equations/linearized_euler_2d.jl
+++ b/src/equations/linearized_euler_2d.jl
@@ -144,6 +144,15 @@ end
     return SVector(f1, f2, f3, f4)
 end
 
+"""
+    have_constant_speed(::LinearizedEulerEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearizedEulerEquations2D) = True()
 
 @inline function max_abs_speeds(equations::LinearizedEulerEquations2D)
diff --git a/src/equations/linearized_euler_3d.jl b/src/equations/linearized_euler_3d.jl
index 7059189eb15..8a48151c84a 100644
--- a/src/equations/linearized_euler_3d.jl
+++ b/src/equations/linearized_euler_3d.jl
@@ -180,6 +180,15 @@ end
     return SVector(f1, f2, f3, f4, f5)
 end
 
+"""
+    have_constant_speed(::LinearizedEulerEquations3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearizedEulerEquations3D) = True()
 
 @inline function max_abs_speeds(equations::LinearizedEulerEquations3D)
diff --git a/src/equations/maxwell_1d.jl b/src/equations/maxwell_1d.jl
index 67d601c8606..70047381ae7 100644
--- a/src/equations/maxwell_1d.jl
+++ b/src/equations/maxwell_1d.jl
@@ -79,6 +79,15 @@ end
     return equations.speed_of_light
 end
 
+"""
+    have_constant_speed(::MaxwellEquations1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::MaxwellEquations1D) = True()
 
 @inline function max_abs_speeds(equations::MaxwellEquations1D)
diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl
index c2ca8e9ffa8..61faefd9fc6 100644
--- a/src/semidiscretization/semidiscretization.jl
+++ b/src/semidiscretization/semidiscretization.jl
@@ -245,7 +245,8 @@ and a vector `b`:
 ```math
 \\partial_t u(t) = A u(t) - b.
 ```
-Works only for linear equations, i.e., equations with `have_constant_speed(equations) == True()`.
+Works only for linear equations, i.e.,
+equations which `have_constant_speed(equations) == True()`.
 
 This has the benefit of greatly reduced memory consumption compared to constructing
 the full system matrix explicitly, as done for instance in
@@ -253,6 +254,14 @@ the full system matrix explicitly, as done for instance in
 
 The returned linear operator `A` is a matrix-free operator which can be
 supplied to iterative solvers from, e.g., [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
+
+It is also possible to use this to construct a sparse matrix without the detour of constructing
+first the full Jacobian by calling
+```julia
+using SparseArrays
+A_map, b = linear_structure(semi, t0 = t0)
+A_sparse = sparse(A_map)
+```
 """
 function linear_structure(semi::AbstractSemidiscretization;
                           t0 = zero(real(semi)))
@@ -268,6 +277,7 @@ function linear_structure(semi::AbstractSemidiscretization;
     u_ode .= zero(eltype(u_ode))
     rhs!(du_ode, u_ode, semi, t0)
     b = -du_ode
+
     # Create a copy of `b` used internally to extract the linear part of `semi`.
     # This is necessary to get everything correct when the user updates the
     # returned vector `b`.
diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
index 86035ca2269..a0e67ce2ded 100644
--- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
+++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
@@ -394,7 +394,8 @@ and a vector `b`:
 ```math
 \\partial_t u(t) = A u(t) - b.
 ```
-Works only for linear equations, i.e., equations with `have_constant_speed(equations) == True()`.
+Works only for linear equations, i.e.,
+equations which `have_constant_speed(equations) == True()` and `have_constant_diffusivity(equations_parabolic) == True()`
 
 This has the benefit of greatly reduced memory consumption compared to constructing
 the full system matrix explicitly, as done for instance in
@@ -405,7 +406,8 @@ supplied to iterative solvers from, e.g., [Krylov.jl](https://github.com/JuliaSm
 """
 function linear_structure(semi::SemidiscretizationHyperbolicParabolic;
                           t0 = zero(real(semi)))
-    if have_constant_speed(semi.equations) == False()
+    if (have_constant_speed(semi.equations) == False() ||
+        have_constant_diffusivity(semi.equations_parabolic) == False())
         throw(ArgumentError("`linear_structure` expects linear equations."))
     end
 
@@ -457,6 +459,70 @@ function _jacobian_ad_forward(semi::SemidiscretizationHyperbolicParabolic, t0, u
     return J
 end
 
+"""
+    linear_structure_parabolic(semi::SemidiscretizationHyperbolicParabolic;
+                               t0 = zero(real(semi)))
+
+Wraps the **parabolic part** right-hand side operator of the hyperbolic-parabolic semidiscretization `semi`
+at time `t0` as an affine-linear operator given by a linear operator `A`
+and a vector `b`:
+```math
+\\partial_t u(t) = A u(t) - b.
+```
+Works only for linear parabolic equations, i.e.,
+equations which `have_constant_diffusivity(equations_parabolic) == True()`.
+
+This has the benefit of greatly reduced memory consumption compared to constructing
+the full system matrix explicitly, as done for instance in
+[`jacobian_ad_forward_parabolic`](@ref).
+
+The returned linear operator `A` is a matrix-free representation which can be
+supplied to iterative solvers from, e.g., [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
+
+It is also possible to use this to construct a sparse matrix without the detour of constructing
+first the full Jacobian by calling
+```julia
+using SparseArrays
+A_map, b = linear_structure_parabolic(semi, t0 = t0)
+A_sparse = sparse(A_map)
+```
+which can then be further used to construct for instance a
+[`MatrixOperator`](https://docs.sciml.ai/SciMLOperators/stable/tutorials/getting_started/#Simplest-Operator:-MatrixOperator)
+from [SciMLOperators.jl](https://docs.sciml.ai/SciMLOperators/stable/).
+This is especially useful for IMEX schemes where the parabolic part is implicitly,
+and for a `MatrixOperator` only factorized once.
+"""
+function linear_structure_parabolic(semi::SemidiscretizationHyperbolicParabolic;
+                                    t0 = zero(real(semi)))
+    if have_constant_diffusivity(semi.equations_parabolic) == False()
+        throw(ArgumentError("`linear_structure_parabolic` expects equations with constant diffusive terms."))
+    end
+
+    # allocate memory
+    u_ode = allocate_coefficients(mesh_equations_solver_cache(semi)...)
+    du_ode = similar(u_ode)
+
+    # get the parabolic right hand side from boundary conditions and optional source terms
+    u_ode .= zero(eltype(u_ode))
+    rhs_parabolic!(du_ode, u_ode, semi, t0)
+    b = -du_ode
+
+    # Create a copy of `b` used internally to extract the linear part of `semi`.
+    # This is necessary to get everything correct when the user updates the
+    # returned vector `b`.
+    b_tmp = copy(b)
+
+    # wrap the linear operator
+    A = LinearMap(length(u_ode), ismutating = true) do dest, src
+        rhs_parabolic!(dest, src, semi, t0)
+
+        @. dest += b_tmp
+        return dest
+    end
+
+    return A, b
+end
+
 """
     jacobian_ad_forward_parabolic(semi::SemidiscretizationHyperbolicParabolic;
                                   t0 = zero(real(semi)),
diff --git a/test/Project.toml b/test/Project.toml
index e2d2243c855..824aade8a6c 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -30,6 +30,8 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
@@ -38,10 +40,10 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TrixiTest = "0a316866-cbd0-4425-8bcb-08103b2c1f26"
 
 [compat]
-ADTypes = "1.11"
-Adapt = "4"
+ADTypes = "1.14"
+Adapt = "4.1"
 Aqua = "0.8"
-CUDA = "5.8"
+CUDA = "5.8.2"
 CairoMakie = "0.12, 0.13, 0.14, 0.15"
 Convex = "0.16"
 DelimitedFiles = "1"
@@ -51,13 +53,13 @@ ECOS = "1.1.2"
 ExplicitImports = "1.0.1"
 FiniteDiff = "2.27.0"
 ForwardDiff = "0.10.36, 1"
-Krylov = "0.9.10, 0.10"
+Krylov = "0.10"
 LinearAlgebra = "1"
-LinearSolve = "2.36.1, 3"
-MPI = "0.20.6"
+LinearSolve = "3.13"
+MPI = "0.20.22"
 NLsolve = "4.5.1"
 OrdinaryDiffEqBDF = "1.1"
-OrdinaryDiffEqCore = "1.6"
+OrdinaryDiffEqCore = "1.26"
 OrdinaryDiffEqFeagin = "1"
 OrdinaryDiffEqHighOrderRK = "1.1"
 OrdinaryDiffEqLowOrderRK = "1.2"
@@ -65,10 +67,12 @@ OrdinaryDiffEqLowStorageRK = "1.2"
 OrdinaryDiffEqSDIRK = "1.1"
 OrdinaryDiffEqSSPRK = "1.2"
 OrdinaryDiffEqTsit5 = "1.1"
-Plots = "1.26"
+Plots = "1.38.9"
 Printf = "1"
 Quadmath = "0.5.10"
 Random = "1"
+SciMLBase = "2.92.0"
+SciMLOperators = "1"
 SparseArrays = "1"
 SparseConnectivityTracer = "1.0.1"
 SparseMatrixColorings = "0.4.21"
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index f69dcd7df5d..a7778021724 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -253,6 +253,15 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "TreeMesh2D: elixir_advection_diffusion_imex_operator.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
+                                 "elixir_advection_diffusion_imex_operator.jl"),
+                        l2=[7.542670562162156e-8], linf=[3.972014046560446e-7])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl (LDG)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                  "elixir_advection_diffusion_nonperiodic.jl"),