API

fcscumulants_recursive(L, mJ, nC, rho_ss, nu; cumulant_type="")
fcscumulants_recursive(H, J, mJ, nC, rho_ss, nu; cumulant_type="")

Calculate n-th zero-frequency cumulant of full counting statistics using a recursive scheme.

Arguments

• L: Vectorized Liouvillian matrix (sparse or dense, ComplexF64)

Alternatively, one can provide the Hamiltonian and jump operators instead of L

• H: Hamiltonian operator (sparse or dense, or a backend operator — QuantumOptics.Operator / QuantumToolbox.QuantumObject)
• J: Vector of jump operators (same accepted types as H)
• mJ: Vector of the monitored jump operators, as ordinary n×n matrices (or backend operators) — not super-operators and not vectorized. They select which jump channels are counted.
• nC: Number of cumulants to be calculated.
• rho_ss: Steady-state density matrix (sparse or dense, ComplexF64)
• nu: Vector of length length(mJ) with weights for each jump. The weights set both the sign convention and the units of the current (±1 for counts, charges for electric currents, energy quanta for heat currents).

Keyword arguments

• method: Drazin-solve backend, :lu (default) or :iterative.
• σ, τ, rtol, itmax, memory: options for the :iterative backend (see prepare_drazin_solver); ignored by :lu.
• Pl: an externally built preconditioner reused by the :iterative backend instead of building an ILU internally (see prepare_drazin_solver). nothing (default) builds the internal preconditioner; ignored by :lu and by the dense-Liouvillian method.
• cumulant_type: Optional keyword. Use "" (default) for ordinary cumulants, or "factorial" to convert the result with factorial_cumulants.

fcscumulants_recursive(problem::FCSProblem)

Solve a full-counting-statistics problem, returning its first problem.nC zero-frequency cumulants. Equivalent to calling the positional method with the problem's fields, after normalizing any backend operators to their underlying matrices.

fcscumulants_recursive(ctx::PreparedLindbladFCS; mJ, nu, nC=2, cumulant_type="")

Evaluate a full-counting-statistics observable on a prepared context, reusing its Drazin solver. Returns the first nC zero-frequency cumulants.

mJ are the monitored jump operators/matrices (backend operators are accepted and normalized) and nu their weights (length(nu) == length(mJ)). Each jump must be n×n, where n is the side of the context's steady state. See prepare_fcs_context for the "prepare once, evaluate many observables" workflow and fcscumulants_recursive for cumulant_type.

factorial_cumulants(cumulants::AbstractVector{<:Number})

Convert ordinary cumulants [c1, c2, ..., cn] to factorial cumulants [f1, f2, ..., fn].

The conversion uses the signed Stirling numbers of the first kind,

f_m = sum(s(m, j) * c_j for j in 1:m)

with s(1, 1) = 1 and the recurrence s(m, j) = s(m - 1, j - 1) - (m - 1) * s(m - 1, j).

The physical factorial-cumulant interpretation assumes dimensionless count cumulants. Dimensionful heat-current cumulants should be rescaled first.

FCSProblem

Abstract supertype for full-counting-statistics problems. A concrete subtype encapsulates everything fcscumulants_recursive needs — the Liouvillian (or the Hamiltonian/jump operators it is built from), the monitored jumps, the steady state, the weights, and the number of cumulants — so that a simulation can be organised as a single object and solved with fcscumulants_recursive(problem).

LindbladFCS(; H=nothing, J=nothing, L=nothing, mJ, rho_ss, nu, nC=2,
            method=:lu, σ=nothing, τ=0.05, Pl=nothing,
            rtol=1e-8, itmax=200, memory=30)
LindbladFCS(H, J; mJ, rho_ss, nu, nC=2, kwargs...)

A full-counting-statistics problem for a Lindblad master equation.

The second form is a convenience constructor provided by the QuantumOptics and QuantumToolbox extensions: it takes the model operators positionally and accepts the same keywords as the first.

Construct it by keyword. Either supply a (vectorized) Liouvillian L, or both a Hamiltonian H and a vector of jump operators J (in which case L is built from them at solve time by the active backend). Fields may be plain dense/sparse ComplexF64 arrays, or backend operators (QuantumOptics.Operator, QuantumToolbox.QuantumObject); the relevant package extension takes care of extracting the underlying matrices.

Fields

• H: optional Hamiltonian operator (backend type) used to build L.
• J: optional vector of jump operators (backend type) used to build L.
• L: optional (vectorized) Liouvillian. Takes precedence over H/J.
• mJ: vector of monitored jump operators/matrices.
• rho_ss: steady-state density matrix.
• nu: weights, one per monitored jump (length(nu) == length(mJ)).
• nC: number of cumulants to compute (default 2).
• method: Drazin-solve backend, :lu (default) or :iterative.
• σ, τ, rtol, itmax, memory: options for the :iterative backend (diagonal shift, ILU drop tolerance, Krylov tolerance, iteration cap, and GMRES restart memory). Ignored by :lu. σ=nothing auto-scales the shift from L.
• Pl: an externally built preconditioner for the :iterative backend to reuse instead of building its own ILU (e.g. the ILU from the steady-state solve). See prepare_drazin_solver for the contract it must satisfy. nothing (default) builds the internal preconditioner; ignored by :lu.

Solve with fcscumulants_recursive:

p = LindbladFCS(; H=H, J=J, mJ=[Jc], rho_ss=ρss, nu=[1], nC=2)
c1, c2 = fcscumulants_recursive(p)

# large sparse Liouvillian (needs `using Krylov, IncompleteLU`):
p = LindbladFCS(; L=L, mJ=[Jc], rho_ss=ρss, nu=[1], nC=3, method=:iterative)

PreparedLindbladFCS

A reusable full-counting-statistics context: the observable-independent data for a fixed Liouvillian and steady state, with a Drazin solver prepared once.

Build it with prepare_fcs_context and evaluate any number of observables that share the same L and rho_ss with

fcscumulants_recursive(ctx; mJ, nu, nC=2)

Because the monitored jumps mJ and weights nu enter only the counting-field super-operators and the recursion's right-hand sides — never the Drazin linear algebra — a single prepared solver serves every observable, so the expensive preparation (sparse LU factorization, or shifted-ILU build for the iterative backend) is done only once.

This complements LindbladFCS: use LindbladFCS for a single observable, and a PreparedLindbladFCS context when many observables share one Liouvillian and steady state.

Fields

• L: the (sparse ComplexF64) vectorized Liouvillian.
• rho_ss: the (sparse ComplexF64) steady-state density matrix, kept for introspection and dimension checks.
• vrho_ss: the normalized, vectorized steady state (right null vector).
• vId: the vectorized identity / trace functional (left null vector).
• solver: the prepared DrazinSolver reused across observables.

prepare_fcs_context(; H=nothing, J=nothing, L=nothing, rho_ss,
                    method=:lu, σ=nothing, τ=0.05, Pl=nothing,
                    rtol=1e-8, itmax=200, memory=30) -> PreparedLindbladFCS

Prepare a reusable PreparedLindbladFCS context for computing many full-counting-statistics observables that share the same Liouvillian and steady state.

Supply either a (vectorized) Liouvillian L, or both a Hamiltonian H and jump operators J (built into L by the active backend). L/H/J/rho_ss may be plain dense/sparse ComplexF64 arrays or backend operators (QuantumOptics.Operator, QuantumToolbox.QuantumObject); the relevant package extension extracts the underlying matrices. The Drazin solver is prepared once here and then reused by every call to fcscumulants_recursive(ctx; mJ, nu, nC).

The solver-backend keywords are identical to those of LindbladFCS / fcscumulants_recursive and configure this one-time preparation:

• method: :lu (default) or :iterative (needs the QuantumFCSIterativeExt extension, using Krylov, IncompleteLU).
• σ, τ, Pl, rtol, itmax, memory: options for the :iterative backend, forwarded to prepare_drazin_solver; ignored by :lu. Pl supplies an externally built preconditioner for the iterative backend to reuse.

Example

ctx  = prepare_fcs_context(; L = L, rho_ss = ρss, method = :lu)   # prepares LU once
hot  = fcscumulants_recursive(ctx; mJ = mJ_hot,  nu = nu_hot,  nC = 2)
cold = fcscumulants_recursive(ctx; mJ = mJ_cold, nu = nu_cold, nC = 2)

prepare_fcs_context(ss::TraceConstrainedSteadyState; method = :iterative,
                    σ = nothing, τ = 0.05, rtol = 1e-8, itmax = 200, memory = 30)

Bridge a solved TraceConstrainedSteadyState into a reusable PreparedLindbladFCS context, forwarding the steady-state preconditioner so the iterative FCS backend does not rebuild an ILU.

With method = :iterative the steady-state ss.Pl is passed through as the FCS preconditioner; with method = :lu the direct backend is used and Pl is ignored. The steady state is not recomputed. The σ/τ/rtol/itmax/memory keywords tune the FCS Drazin solve exactly as in the keyword prepare_fcs_context.

trace_constrained_steadystate(L; method = :iterative, kwargs...)
trace_constrained_steadystate(H, J; method = :iterative, kwargs...)
trace_constrained_steadystate(sys::TraceConstrainedSystem; method = :iterative, kwargs...)

Solve for the trace-constrained steady state, returning a lean TraceConstrainedSteadyState that carries rho_ss and the reusable preconditioner Pl.

method = :lu solves the sparse system directly (A \ b) and returns Pl = nothing — a reliable baseline for small/medium systems. method = :iterative builds a shifted-ILU preconditioner and solves with GMRES (requires the QuantumFCSIterativeExt extension); the resulting Pl is meant to be reused by the iterative FCS backend.

Iterative keyword arguments (ignored by :lu):

• Pl — reuse an externally built preconditioner instead of building one.
• u0 — initial guess for GMRES (warm start, e.g. a neighbouring solution).
• τ, shift_factor, shift — shifted_ilu_preconditioner options.
• rtol (1e-10), atol (1e-14), itmax (200), memory (60) — GMRES tolerances, iteration cap, and Krylov basis size (defaults match the applications).

Pass the low-level sys::TraceConstrainedSystem form to reuse a prebuilt system and preconditioner across a continuation sweep without rebuilding either.

TraceConstrainedSteadyState

The lean, user-facing result of trace_constrained_steadystate: a steady state together with the preconditioner built to find it, ready to feed into the iterative FCS machinery via prepare_fcs_context.

It intentionally does not retain the large trace-constrained matrix A or the right-hand side b — the target problems are memory critical.

Fields

• L: the sparse ComplexF64 Liouvillian used for FCS.
• rho_ss: the steady state as a sparse ComplexF64 matrix, hermitianized and trace-normalized.
• Pl: the preconditioner built for the steady-state solve (nothing for the direct :lu method). Reused by the :iterative FCS backend — this reuse is the whole point of the API.
• stats: a NamedTuple of scalar diagnostics only (convergence, iterations, residuals, trace/hermiticity errors, ILU/GMRES timings) — no large arrays.

trace_constrained_system(L; weight = nothing, dimensions = nothing) -> TraceConstrainedSystem
trace_constrained_system(H, J; weight = nothing, dimensions = nothing) -> TraceConstrainedSystem

Build the trace-constrained steady-state linear system for a Liouvillian.

Supply a (vectorized) Liouvillian L, or a Hamiltonian H and jump operators J (built into L by the active backend). Inputs may be plain dense/sparse ComplexF64 arrays or backend operators (QuantumOptics, QuantumToolbox); they are normalized to SparseMatrixCSC{ComplexF64, Int}.

The default constraint weight is norm(L, 1) / length(L) (the entrywise 1-norm of L divided by its total number of entries), matching the application helpers; pass weight to override it.

dimensions is carried through to the returned system's dimensions field for your own bookkeeping (e.g. the Hilbert-space factor sizes needed to rewrap rho_ss as a backend operator). It is not inferred from the inputs — see TraceConstrainedSystem.

TraceConstrainedSystem

The trace-constrained steady-state linear system A ρ⃗ = b for a Liouvillian L.

L is singular (the physical steady state is its null vector), so the redundant first equation is replaced by a trace constraint: A = L + |row 1⟩⟨diag|·w and b = w·e₁, where the constraint hits the vectorized-identity (diagonal) indices so that solving imposes tr(ρ) = 1.

This is a low-level object for constructing, testing, and warm-started continuation solves. Prefer the lean TraceConstrainedSteadyState as the user-facing result. Build one with trace_constrained_system.

Fields

• L: the sparse ComplexF64 Liouvillian (unmodified; this is what FCS uses).
• A: the trace-constrained system matrix L + trace_constraint.
• b: the right-hand side (w in the first entry, zeros elsewhere).
• vId: the vectorized identity / trace functional (diagonal indices).
• dimensions: Hilbert-space dimensions, for rewrapping rho_ss as a backend operator. Reserved: no backend currently populates it, so it is nothing unless you pass dimensions to trace_constrained_system yourself. (A backend override is not portable — a QuantumToolbox super-operator reports Liouville-space dimensions, not the Hilbert-space ones needed here.)

shifted_ilu_preconditioner(A; τ = 1e-3, shift_factor = 1e-6, shift = nothing)

Build a shifted incomplete-LU preconditioner for the trace-constrained matrix A.

A small diagonal shift keeps the incomplete factorization well-behaved: ilu(A + shift·I; τ). If shift is nothing it is scaled from the operator as shift_factor · max(norm(A, 1) / size(A, 1), eps(Float64)), matching the application helpers. The shift affects only the preconditioner, never the solved system.

Requires the QuantumFCSIterativeExt extension (using Krylov, IncompleteLU).

drazin(L, vrho_ss, vId, IdL)

Calculate the (projected) Drazin inverse of the Liouvillian L.

When a backend extension is loaded, a convenience overload drazin(H, J, vrho_ss, vId, IdL) builds L from a Hamiltonian and jump operators first.

Arguments

• L : Liouvillian matrix
• vrho_ss: vectorised density matrix specifying the steady-state of the Liouvillian.
• vId: vectorised identity matrix (1×N row or vector)
• IdL: Identity matrix in Liouville space (N×N)

Returns

A dense matrix holding the (projected) Drazin inverse Lᴰ of L.

drazin_apply(L, α, ρ, vId; F=nothing, rtol=1e-12, atol=0.0)

Apply the (projected) Drazin inverse of the Liouvillean L to the vector α by solving a linear system.

Arguments

• L: Liouvillean operator (matrix).
• α: Right-hand side vector.
• ρ: Steady-state vector.
• vId: Vectorized identity vector.
• F: Optional factorization of L to reuse (default: nothing).
• rtol: Relative tolerance for the solver (default: 1e-12).
• atol: Absolute tolerance for the solver (default: 0.0).

Returns

A (sparse) vector representing the result of applying the projected Drazin inverse.

m_jumps(mJ; n=1, nu=vcat(fill(+1, length(mJ)÷2), fill(-1, length(mJ)÷2)))

Calculate the vectorized super-operator ℒ(n) = ∑ₖ (νₖ)ⁿ (Lₖ*)⊗Lₖ.

Arguments

• mJ: List of monitored jumps
• n : Power of the weights νₖ. By default set to 1, since this case appears more often.
• nu: vector of length length(mJ) with weights for each jump operator.