diff --git a/Project.toml b/Project.toml
index 0d812871..ac47b867 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.1.0"
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
@@ -23,6 +24,7 @@ CheckConcreteStructs = "0.1"
 DataInterpolations = "8"
 DiffEqNoiseProcess = "5.27.0"
 ExplicitImports = "1.11"
+ForwardDiff = "1"
 FunctionWrappers = "1"
 JET = "0.9.18, 0.10, 0.11"
 LinearAlgebra = "1.10"
diff --git a/docs/Project.toml b/docs/Project.toml
index 99b50570..468c038a 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -19,6 +20,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 [compat]
 DataInterpolations = "8"
 Documenter = "1"
+ForwardDiff = "1"
 LaTeXStrings = "1"
 MacroTools = "0.5"
 ModelingToolkitBase = "1"
diff --git a/docs/make.jl b/docs/make.jl
index df2e65e1..c1d9c2a7 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -30,6 +30,7 @@ pages = [
         "examples/07-2_hong-ou-mandel__quantum-pulse.md",
         "examples/08-1_pulse-delay__simple.md",
         "examples/09-1_coherent-feedback-squeezing__Gough-Wildfeuer-2009.md",
+        "examples/10-1_quantum-fisher-information__automatic-differentiation.md",
     ],
 ]
 
diff --git a/docs/src/api.md b/docs/src/api.md
index 1b8eca16..41730cba 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -94,6 +94,28 @@ solve_mode_evolution_symmetric
 correlation_matrix
 ```
 
+## [Information](@id API: Information)
+
+```@docs
+quantum_fisher_information
+```
+
+```@docs
+classical_fisher_information
+```
+
+```@docs
+povm_probabilities
+```
+
+```@docs
+projective_measurement
+```
+
+```@docs
+parameter_derivative
+```
+
 ## [Utilities](@id API: Utilities)
 
 ```@docs
diff --git a/examples/10-1_quantum-fisher-information__automatic-differentiation.jl b/examples/10-1_quantum-fisher-information__automatic-differentiation.jl
new file mode 100644
index 00000000..26615504
--- /dev/null
+++ b/examples/10-1_quantum-fisher-information__automatic-differentiation.jl
@@ -0,0 +1,205 @@
+# # Quantum Fisher Information with Automatic Differentiation
+#
+# Quantum Fisher information (QFI) answers an optimistic question: if a quantum
+# state ``\rho(\theta)`` depends on a parameter ``\theta``, how much information
+# about ``\theta`` is available in the state before we commit to a particular
+# measurement? For one copy of the state it gives the quantum Cramér-Rao bound
+#
+# ```math
+# \mathrm{Var}(\hat\theta) \geq \frac{1}{F_Q}.
+# ```
+#
+# The bound is useful only if we can also ask a more practical question: how much
+# information is obtained by the measurements we can actually perform? For a POVM
+# or projective measurement with outcomes ``M_i``,
+#
+# ```math
+# p_i(\theta) = \mathrm{Tr}[M_i \rho(\theta)],
+# \qquad
+# F_C = \sum_i \frac{(\partial_\theta p_i)^2}{p_i}.
+# ```
+#
+# The classical Fisher information (CFI) satisfies ``F_C \leq F_Q``. Comparing
+# the two tells us whether a measurement is close to optimal or whether useful
+# information remains hidden in coherences that the measurement does not access.
+#
+# In this example we estimate the detuning of a driven Kerr parametric
+# oscillator. The derivative ``\partial_\Delta \rho(t)`` is obtained by
+# differentiating the SLH time evolution with `ForwardDiff.jl`. We then compare
+# the QFI with the CFI of three different projective readouts.
+
+using QuantumInputOutput
+using SecondQuantizedAlgebra
+using QuantumOptics
+using LinearAlgebra
+using Plots
+
+#
+
+# ## SLH Model
+#
+# The Kerr parametric oscillator has Hamiltonian
+#
+# ```math
+# H = -\Delta\,a^\dagger a
+#     + K\,a^\dagger a^\dagger a a
+#     - G\,(a^\dagger a^\dagger + a a),
+# ```
+#
+# and one decay channel ``L = \sqrt{\gamma}\,a``. We keep the parameters
+# symbolic in the SLH model and provide numerical values only when we translate
+# and evolve it.
+
+h = FockSpace(:kpo)
+a = Destroy(h, :a)
+
+@variables Δ::Real K::Real G::Real γ::Real
+
+H = -Δ * a' * a + K * a' * a' * a * a - G * (a' * a' + a * a)
+Gkpo = SLH(1, √(γ) * a, H)
+nothing # hide
+
+#
+
+const K_ = 0.001
+const G_ = 0.002
+const γ_ = 0.01
+const N = 20
+
+b = FockBasis(N - 1)
+ψ0 = fockstate(b, 0)
+ρ0 = dm(ψ0)
+
+params(Δ_) = Dict(Δ => Δ_, K => K_, G => G_, γ => γ_)
+nothing # hide
+
+# ## Candidate Measurements
+#
+# `classical_fisher_information` works with any list of POVM elements. For
+# projective readout, [`projective_measurement`](@ref) builds the projectors from
+# a Hermitian observable. This keeps the package API independent of the physical
+# system: the readout below happens to use a truncated oscillator basis, but the
+# Fisher-information code only sees projectors.
+#
+# We compare:
+#
+# - occupation readout, resolving each basis state;
+# - parity readout, which intentionally coarse-grains occupation into even and
+#   odd outcomes;
+# - quadrature readout, here represented by the eigenprojectors of ``X=a+a^\dagger``
+#   in the same truncated basis.
+
+aqo = destroy(b)
+
+occupation = Operator(b, b, Diagonal(collect(0:(N - 1))))
+parity = Operator(b, b, Diagonal((-1.0) .^ collect(0:(N - 1))))
+quadrature = aqo + aqo'
+
+measurements = (
+    occupation = projective_measurement(occupation),
+    parity = projective_measurement(parity),
+    quadrature = projective_measurement(quadrature),
+)
+
+function cfi_for(measurement, ρ, dρ)
+    return classical_fisher_information(
+        ρ,
+        dρ,
+        measurement;
+        probability_floor = 1e-12,
+        derivative_floor = 1e-6,
+    )
+end
+nothing # hide
+
+# ## Differentiating the SLH Trajectory
+#
+# [`parameter_derivative`](@ref) evolves the density matrix to a chosen time and
+# returns both ``\rho(t)`` and ``\partial_\Delta\rho(t)``. The `estimate` keyword
+# names the symbolic parameter to differentiate, while `parameter` fixes the
+# numerical working point.
+
+ρ, dρ = parameter_derivative(Gkpo, b, ρ0, 2000; estimate = Δ, parameter = params(0.0))
+
+final_information = (
+    QFI = quantum_fisher_information(ρ, dρ),
+    occupation = cfi_for(measurements.occupation, ρ, dρ),
+    parity = cfi_for(measurements.parity, ρ, dρ),
+    quadrature = cfi_for(measurements.quadrature, ρ, dρ),
+)
+@info "Fisher information at final time" final_information
+nothing # hide
+
+# ## Readout Comparison
+#
+# Repeating the differentiated solve over a grid of readout times gives a direct
+# comparison between the measurement-independent QFI bound and concrete
+# measurement choices. The full occupation measurement keeps more information
+# than parity because it resolves the individual outcomes that parity combines.
+# Quadrature probes a different basis and can reveal parameter sensitivity that
+# is not visible in occupation probabilities alone.
+
+ts = range(0, 2000; length = 100)
+
+F_t = map(ts) do t
+    ρ, dρ = parameter_derivative(Gkpo, b, ρ0, t; estimate = Δ, parameter = params(0.0))
+    (
+        quantum = quantum_fisher_information(ρ, dρ),
+        occupation = cfi_for(measurements.occupation, ρ, dρ),
+        parity = cfi_for(measurements.parity, ρ, dρ),
+        quadrature = cfi_for(measurements.quadrature, ρ, dρ),
+    )
+end
+
+QFI_t = getproperty.(F_t, :quantum)
+CFI_occupation_t = getproperty.(F_t, :occupation)
+CFI_parity_t = getproperty.(F_t, :parity)
+CFI_quadrature_t = getproperty.(F_t, :quadrature)
+nothing # hide
+
+# ## Visualization
+
+p_info = plot(
+    ts,
+    QFI_t;
+    lw = 2.5,
+    label = "QFI bound",
+    xlabel = "time",
+    ylabel = "Fisher information",
+    grid = true,
+    size = (700, 400),
+)
+plot!(p_info, ts, CFI_occupation_t; lw = 2, ls = :dash, label = "occupation CFI")
+plot!(p_info, ts, CFI_parity_t; lw = 2, ls = :dot, label = "parity CFI")
+plot!(p_info, ts, CFI_quadrature_t; lw = 2, ls = :dashdot, label = "quadrature CFI")
+p_info
+
+#
+
+# At the final time we can also inspect the measurement probabilities. A readout
+# with high CFI is not merely one with a broad distribution: it must have outcome
+# probabilities that change strongly when the estimated parameter changes.
+
+occupation_probabilities = povm_probabilities(ρ, measurements.occupation)
+
+p_probs = bar(
+    0:(N - 1),
+    occupation_probabilities;
+    xlabel = "occupation outcome",
+    ylabel = "probability",
+    label = "",
+    grid = false,
+    size = (700, 300),
+)
+p_probs
+
+# ## Package versions
+
+using InteractiveUtils
+versioninfo()
+
+using Pkg
+Pkg.status(
+    ["QuantumInputOutput", "SecondQuantizedAlgebra", "QuantumOptics", "ForwardDiff", "Plots"],
+    mode = PKGMODE_MANIFEST,
+)
diff --git a/examples/Project.toml b/examples/Project.toml
index 30a92a5e..250cb67d 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -18,6 +19,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Documenter = "1"
+ForwardDiff = "1"
 LaTeXStrings = "1"
 MacroTools = "0.5"
 ModelingToolkitBase = "1"
@@ -26,6 +28,6 @@ Plots = "1"
 QuantumCumulants = "0.5"
 QuantumOptics = "1"
 QuantumOpticsBase = "0.5"
-SecondQuantizedAlgebra = "0.6"
+SecondQuantizedAlgebra = "0.8.2"
 SymbolicUtils = "4"
 Symbolics = "7"
diff --git a/src/QuantumInputOutput.jl b/src/QuantumInputOutput.jl
index 9ee3a8f2..0048ebea 100644
--- a/src/QuantumInputOutput.jl
+++ b/src/QuantumInputOutput.jl
@@ -12,6 +12,7 @@ using LinearAlgebra: LinearAlgebra, I, mul!
 using OrdinaryDiffEq: OrdinaryDiffEq, ODEProblem, Tsit5, solve
 using StaticArrays: StaticArrays, SMatrix, SVector
 using FunctionWrappers: FunctionWrappers, FunctionWrapper
+using ForwardDiff: ForwardDiff
 
 const SQA = SecondQuantizedAlgebra
 
@@ -42,6 +43,12 @@ export SLH,
     solve_mode_evolution_symmetric,
     # Correlations
     correlation_matrix,
+    # Information
+    quantum_fisher_information,
+    classical_fisher_information,
+    povm_probabilities,
+    projective_measurement,
+    parameter_derivative,
     # Operators
     substitute_operators
 
@@ -51,5 +58,6 @@ include("utils.jl")
 include("pulses.jl")
 include("correlations.jl")
 include("interaction_picture.jl")
+include("information.jl")
 
 end
diff --git a/src/information.jl b/src/information.jl
new file mode 100644
index 00000000..26efa737
--- /dev/null
+++ b/src/information.jl
@@ -0,0 +1,195 @@
+"""
+    quantum_fisher_information(ρ, dρ; regularization=1e-12)
+
+Compute the quantum Fisher information from a density matrix `ρ` and its
+parameter derivative `dρ`.
+
+The implementation uses the symmetric logarithmic derivative `L`, defined by
+`dρ = (ρL + Lρ) / 2`, and returns `real(tr(dρ * L))`. A small diagonal
+`regularization` is added to `ρ` before solving the Sylvester equation so that
+rank-deficient states are numerically well behaved.
+"""
+function quantum_fisher_information(
+    ρ::AbstractMatrix,
+    dρ::AbstractMatrix;
+    regularization = 1e-12,
+)
+    size(ρ, 1) == size(ρ, 2) || throw(ArgumentError("`ρ` must be square."))
+    size(dρ) == size(ρ) || throw(ArgumentError("`dρ` must have the same size as `ρ`."))
+
+    ρ_reg = ρ + regularization * LinearAlgebra.I
+    L = LinearAlgebra.sylvester(ρ_reg, ρ_reg, -2 * dρ)
+    return real(LinearAlgebra.tr(dρ * L))
+end
+
+function quantum_fisher_information(
+    ρ::QuantumOpticsBase.AbstractOperator,
+    dρ::QuantumOpticsBase.AbstractOperator;
+    kwargs...,
+)
+    return quantum_fisher_information(Matrix(ρ.data), Matrix(dρ.data); kwargs...)
+end
+
+_matrix_data(A::AbstractMatrix) = Matrix(A)
+_matrix_data(A::QuantumOpticsBase.AbstractOperator) = Matrix(A.data)
+
+function _assert_square_same_size(A, B; names = ("first argument", "second argument"))
+    size(A, 1) == size(A, 2) || throw(ArgumentError("`$(names[1])` must be square."))
+    size(B) == size(A) || throw(ArgumentError("`$(names[2])` must have the same size as `$(names[1])`."))
+    return nothing
+end
+
+"""
+    povm_probabilities(ρ, measurement)
+
+Return the probabilities ``p_i = Tr[M_i ρ]`` for a POVM or projective
+measurement `measurement = [M₁, M₂, ...]`.
+
+`ρ` and each measurement element may be matrices or `QuantumOpticsBase`
+operators.
+"""
+function povm_probabilities(ρ, measurement)
+    ρm = _matrix_data(ρ)
+    return [real(LinearAlgebra.tr(_matrix_data(M) * ρm)) for M in measurement]
+end
+
+"""
+    classical_fisher_information(ρ, dρ, measurement; probability_floor=1e-12,
+                                 derivative_floor=1e-12)
+
+Compute the classical Fisher information for a concrete POVM or projective
+measurement `measurement`.
+
+For ``p_i = Tr[M_i ρ]`` and ``∂p_i = Tr[M_i ∂ρ]``, this returns
+
+```math
+F_C = \\sum_i \\frac{(\\partial_\\theta p_i)^2}{p_i}.
+```
+
+Outcomes with probability below `probability_floor` and derivative below
+`derivative_floor` are ignored. If an impossible outcome has a nonzero
+derivative, the Fisher information is infinite.
+"""
+function classical_fisher_information(
+    ρ,
+    dρ,
+    measurement;
+    probability_floor = 1e-12,
+    derivative_floor = 1e-12,
+)
+    ρm = _matrix_data(ρ)
+    dρm = _matrix_data(dρ)
+    _assert_square_same_size(ρm, dρm; names = ("ρ", "dρ"))
+
+    F = 0.0
+    for M in measurement
+        Mm = _matrix_data(M)
+        _assert_square_same_size(ρm, Mm; names = ("ρ", "measurement element"))
+        p = real(LinearAlgebra.tr(Mm * ρm))
+        dp = real(LinearAlgebra.tr(Mm * dρm))
+
+        if p < -probability_floor
+            throw(DomainError(p, "measurement probabilities must be nonnegative"))
+        elseif p <= probability_floor
+            abs(dp) <= derivative_floor && continue
+            return Inf
+        end
+        F += abs2(dp) / p
+    end
+    return F
+end
+
+"""
+    projective_measurement(A; atol=1e-10)
+
+Return the projectors onto the eigenspaces of a Hermitian observable `A`.
+Degenerate eigenvalues within `atol` are grouped into one projector.
+"""
+function projective_measurement(A::AbstractMatrix; atol = 1e-10)
+    isapprox(A, A'; atol) || throw(ArgumentError("observable must be Hermitian."))
+    F = LinearAlgebra.eigen(LinearAlgebra.Hermitian(Matrix(A)))
+    projectors = Matrix{eltype(F.vectors)}[]
+    used = falses(length(F.values))
+
+    for i in eachindex(F.values)
+        used[i] && continue
+        inds = findall(j -> !used[j] && abs(F.values[j] - F.values[i]) <= atol, eachindex(F.values))
+        P = zeros(eltype(F.vectors), size(A, 1), size(A, 2))
+        for j in inds
+            v = F.vectors[:, j]
+            P += v * v'
+            used[j] = true
+        end
+        push!(projectors, P)
+    end
+    return projectors
+end
+
+function projective_measurement(A::QuantumOpticsBase.AbstractOperator; kwargs...)
+    return [
+        QuantumOpticsBase.Operator(A.basis_l, A.basis_r, P)
+        for P in projective_measurement(Matrix(A.data); kwargs...)
+    ]
+end
+
+function _split_complex_matrix(v, d)
+    half = length(v) ÷ 2
+    return reshape(v[1:half] + im * v[(half + 1):end], d, d)
+end
+
+function _density_state_vector(ρ)
+    ρvec = vec(Matrix(ρ.data))
+    return vcat(real.(ρvec), imag.(ρvec))
+end
+
+function _evolve_density(H, J, ρ0, t; kwargs...)
+    iszero(t) && return ρ0
+    _, ρt = timeevolution.master((zero(float(real(t))), t), ρ0, H, J; saveat = [t], save_everystep = false, kwargs...)
+    return last(ρt)
+end
+
+"""
+    parameter_derivative(G::SLH, b, ρ0, t; estimate, parameter=Dict(), operators=Dict(), kwargs...)
+
+Evolve the SLH model `G` to time `t` and return `(ρ, dρ)`, where `dρ` is the
+derivative of the density matrix with respect to the scalar symbolic parameter
+`estimate`.
+
+The keyword `parameter` gives the numerical working point and must include
+`estimate`. The Hamiltonian and Lindblad operators are obtained with
+[`translate_qo`](@ref), so `estimate` may appear in `H` or `L`. The derivative is
+computed with `ForwardDiff` through the `QuantumOptics.timeevolution.master`
+solve. Extra keyword arguments are forwarded to the time-evolution solver.
+"""
+function parameter_derivative(
+    G::SLH,
+    b::QuantumOpticsBase.Basis,
+    ρ0::QuantumOpticsBase.AbstractOperator,
+    t;
+    estimate,
+    parameter = Dict(),
+    operators = Dict(),
+    kwargs...,
+)
+    haskey(parameter, estimate) || throw(ArgumentError("`parameter` must contain the estimated parameter `$estimate`."))
+    θ0 = parameter[estimate]
+    θ0 isa Real || throw(ArgumentError("`estimate` must have a real scalar value, got `$θ0`."))
+    d = size(ρ0.data, 1)
+
+    function state_vector(θ)
+        p = Dict{Any,Any}(parameter)
+        p[estimate] = θ[1]
+        H, J = translate_qo(G, b; parameter = p, operators)
+        return _density_state_vector(_evolve_density(H, J, ρ0, t; kwargs...))
+    end
+
+    θ = [θ0]
+    ρvec = state_vector(θ)
+    dρvec = vec(ForwardDiff.jacobian(state_vector, θ))
+    ρ = _split_complex_matrix(ρvec, d)
+    dρ = _split_complex_matrix(dρvec, d)
+    return (
+        QuantumOpticsBase.Operator(ρ0.basis_l, ρ0.basis_r, ρ),
+        QuantumOpticsBase.Operator(ρ0.basis_l, ρ0.basis_r, dρ),
+    )
+end
diff --git a/src/translate.jl b/src/translate.jl
index 12c1de4f..33a95989 100644
--- a/src/translate.jl
+++ b/src/translate.jl
@@ -24,21 +24,14 @@ function _translate_numeric(
     operators = Dict{QSym,Any}(),
     op_type = sparse,
 )
-    op_ = substitute(op, parameter)
-    iszero(op_) && return op_type(0 * _one_op(b, operators))
     result = nothing
-    for (term, c_) in op_.arguments
-        c = _coeff_num(c_)
-        _coeff_is_const(c) || throw(
-            ArgumentError(
-                "cannot translate `$op` to a static operator: coefficient `$c` still " *
-                "depends on a free variable; supply it via `parameter` or `time_parameter`.",
-            ),
-        )
-        contrib = _const_coeff(c) * _numeric_product(term.ops, b, operators, op_type)
+    for (term, c_) in op.arguments
+        coeff = _coeff_value(_coeff_num(c_), parameter)
+        iszero(coeff) && continue
+        contrib = coeff * _numeric_product(term.ops, b, operators, op_type)
         result = result === nothing ? contrib : result + contrib
     end
-    return op_type(result)
+    return result === nothing ? op_type(0 * _one_op(b, operators)) : op_type(result)
 end
 
 _translate_numeric_raw(op, b; operators = Dict{QSym,Any}(), op_type = sparse) =
@@ -147,6 +140,32 @@ _coeff_num(c::Complex{Num}) = c
 _coeff_is_const(c::Complex{Num}) =
     isempty(Symbolics.get_variables(real(c))) && isempty(Symbolics.get_variables(imag(c)))
 
+function _coeff_variables(c::Complex{Num})
+    vars = Any[]
+    for v in Symbolics.get_variables(real(c))
+        push!(vars, Symbolics.wrap(v))
+    end
+    for v in Symbolics.get_variables(imag(c))
+        push!(vars, Symbolics.wrap(v))
+    end
+    return unique(vars)
+end
+
+function _parameter_value(parameter, var)
+    haskey(parameter, var) && return parameter[var]
+    uvar = SymbolicUtils.unwrap(var)
+    for (k, v) in parameter
+        isequal(SymbolicUtils.unwrap(k), uvar) && return v
+    end
+    throw(
+        ArgumentError(
+            "cannot translate coefficient depending on `$var`: supply it via `parameter`.",
+        ),
+    )
+end
+
+_same_symbol(a, b) = isequal(SymbolicUtils.unwrap(a), SymbolicUtils.unwrap(b))
+
 # Reduce a concrete `Complex{Num}` to a plain Julia number. Coefficients that are
 # already numeric take a fast path; constant symbolic expressions (e.g. `exp(0.5im)`
 # produced by substituting a numeric value into `exp(im*ϕ)`) are compiled and evaluated.
@@ -182,19 +201,38 @@ function _compile_coeff(c::Complex{Num}, vars...)
     return (vals...) -> g_re(vals...) + im * g_im(vals...)
 end
 
+function _coeff_value(c::Complex{Num}, parameter)
+    vars = _coeff_variables(c)
+    if isempty(vars)
+        return _const_coeff(c)
+    end
+    values = map(var -> _parameter_value(parameter, var), vars)
+    pref = _compile_coeff(c, vars...)
+    return pref(values...)
+end
+
+function _time_coeff_function(c::Complex{Num}, time_parameter, parameter)
+    basevars, valuefuncs = _time_basis(time_parameter)
+    staticvars = filter(_coeff_variables(c)) do var
+        !any(base -> _same_symbol(var, base), basevars)
+    end
+    staticvals = map(var -> _parameter_value(parameter, var), staticvars)
+    pref = _compile_coeff(c, basevars..., staticvars...)
+    return t -> pref(map(g -> g(t), valuefuncs)..., staticvals...)
+end
+
 # ── Per-term translation ──
 # Translate a single `(ops, coefficient)` term into a time-dependent function
 # `t -> op`. A concrete coefficient yields a constant function; a symbolic
 # coefficient is compiled against the time-parameter base variables.
-function _translate_term(ops, c::Complex{Num}, b, time_parameter, operators, op_type)
+function _translate_term(ops, c::Complex{Num}, b, time_parameter, parameter, operators, op_type)
     prodop = _numeric_product(ops, b, operators, op_type)
     if _coeff_is_const(c)
         op = _const_coeff(c) * prodop
         return t -> op
     end
-    basevars, valuefuncs = _time_basis(time_parameter)
-    pref = _compile_coeff(c, basevars...)
-    return t -> pref(map(g -> g(t), valuefuncs)...) * prodop
+    pref = _time_coeff_function(c, time_parameter, parameter)
+    return t -> pref(t) * prodop
 end
 
 """
@@ -258,6 +296,7 @@ function _translate_qo(
             _coeff_num(c),
             b,
             time_parameter,
+            parameter,
             operators,
             op_type,
         )
@@ -270,6 +309,7 @@ function _translate_qo(
         _coeff_num(pairs[1].second),
         b,
         time_parameter,
+        parameter,
         operators,
         sparse,
     )
@@ -284,6 +324,7 @@ function _translate_qo(
                 _coeff_num(pairs[k].second),
                 b,
                 time_parameter,
+                parameter,
                 operators,
                 sparse,
             )
@@ -322,27 +363,19 @@ function _translate_qo(
     operators = Dict{QSym,Any}(),
     op_type = sparse,
 )
-    arg_c = substitute(arg_c_, parameter)
     one_b = _translate_one(b, operators, op_type)
-    c = _as_cnum(arg_c)
+    c = _as_cnum(arg_c_)
 
     if isempty(time_parameter)
-        _coeff_is_const(c) || throw(
-            ArgumentError(
-                "cannot translate symbolic scalar `$arg_c` without a value: supply it via " *
-                "`parameter` (numeric) or `time_parameter` (time-dependent).",
-            ),
-        )
-        return _const_coeff(c) * one_b
+        return _coeff_value(c, parameter) * one_b
     end
 
     if _coeff_is_const(c)
         val = _const_coeff(c)
         return t -> val * one_b
     end
-    basevars, valuefuncs = _time_basis(time_parameter)
-    pref = _compile_coeff(c, basevars...)
-    return t -> pref(map(g -> g(t), valuefuncs)...) * one_b
+    pref = _time_coeff_function(c, time_parameter, parameter)
+    return t -> pref(t) * one_b
 end
 
 function translate_qo(ops::Vector, b::QuantumOpticsBase.Basis; kwargs...)
diff --git a/test/runtests.jl b/test/runtests.jl
index 2f6c82f5..89feeece 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,7 @@ names = [
     "test_example_cavity_scattering.jl",
     "test_compare_example_05_1_05_2.jl",
     "test_interaction_picture.jl",
+    "test_information.jl",
     "test_utils.jl",
 ]
 
diff --git a/test/test_information.jl b/test/test_information.jl
new file mode 100644
index 00000000..531d6a7f
--- /dev/null
+++ b/test/test_information.jl
@@ -0,0 +1,70 @@
+using Test
+using LinearAlgebra
+using QuantumInputOutput
+using QuantumOptics
+using SecondQuantizedAlgebra
+
+@testset "quantum_fisher_information" begin
+    ρ = Diagonal([0.3, 0.7])
+    dρ = Diagonal([1.0, -1.0])
+    @test quantum_fisher_information(Matrix(ρ), Matrix(dρ); regularization = 0.0) ≈
+          1 / 0.3 + 1 / 0.7
+
+    θ = 0.37
+    c, s = cos(θ), sin(θ)
+    ρ_pure = [c^2 c*s; c*s s^2]
+    dρ_pure = [-2c*s c^2-s^2; c^2-s^2 2c*s]
+    @test quantum_fisher_information(ρ_pure, dρ_pure) ≈ 4 atol = 1e-8
+
+    b = NLevelBasis(2)
+    ρ_op = Operator(b, b, Matrix(ρ))
+    dρ_op = Operator(b, b, Matrix(dρ))
+    @test quantum_fisher_information(ρ_op, dρ_op; regularization = 0.0) ≈
+          quantum_fisher_information(Matrix(ρ), Matrix(dρ); regularization = 0.0)
+
+    @test_throws ArgumentError quantum_fisher_information(ones(2, 3), ones(2, 3))
+    @test_throws ArgumentError quantum_fisher_information(ones(2, 2), ones(3, 3))
+end
+
+@testset "classical_fisher_information" begin
+    ρ = [0.25 0; 0 0.75]
+    dρ = [0.5 0; 0 -0.5]
+    M = projective_measurement([0.0 0; 0 1.0])
+
+    @test povm_probabilities(ρ, M) ≈ [0.25, 0.75]
+    @test classical_fisher_information(ρ, dρ, M) ≈ 0.5^2 / 0.25 + 0.5^2 / 0.75
+
+    σx = [0 1; 1 0]
+    P = projective_measurement(σx)
+    @test length(P) == 2
+    @test sum(P) ≈ I(2)
+    @test all(Pi -> Pi * Pi ≈ Pi, P)
+
+    @test classical_fisher_information([1.0 0; 0 0], [0.0 0; 0 1.0], M) == Inf
+    @test_throws DomainError classical_fisher_information([-0.1 0; 0 1.1], dρ, M)
+end
+
+@testset "SLH parameter_derivative" begin
+    h = FockSpace(:c)
+    a = Destroy(h, :a)
+    @variables Δ::Real
+
+    b = FockBasis(1)
+    ψ0 = (fockstate(b, 0) + fockstate(b, 1)) / sqrt(2)
+    ρ0 = dm(ψ0)
+    G = SLH(1, 0 * a, Δ * a' * a)
+
+    t = 0.7
+    θ = 0.2
+    ρ, dρ = parameter_derivative(G, b, ρ0, t; estimate = Δ, parameter = Dict(Δ => θ))
+
+    phase = exp(1im * θ * t)
+    expected_ρ = [0.5 0.5 * phase; 0.5 * conj(phase) 0.5]
+    expected_dρ = [0 0.5im * t * phase; -0.5im * t * conj(phase) 0]
+    @test Matrix(ρ.data) ≈ expected_ρ atol = 1e-8
+    @test Matrix(dρ.data) ≈ expected_dρ atol = 1e-8
+    @test quantum_fisher_information(ρ, dρ) ≈ t^2 atol = 1e-6
+    @test classical_fisher_information(ρ, dρ, projective_measurement([0.0 0; 0 1.0])) ≈ 0 atol = 1e-10
+
+    @test_throws ArgumentError parameter_derivative(G, b, ρ0, t; estimate = Δ, parameter = Dict())
+end