Quantum Fisher Information

author: Orjan Ameye

This notebook demonstrates the computation of Quantum Fisher Information (QFI) for a driven-dissipative Kerr Parametric Oscillator using automatic differentiation. The QFI quantifies the ultimate precision limit for parameter estimation in quantum systems.

We import the necessary packages for quantum simulations and automatic differentiation:

using QuantumToolbox      # Quantum optics simulations
using DifferentiationInterface  # Unified automatic differentiation interface
using SciMLSensitivity   # Allows for ODE sensitivity analysis
using FiniteDiff         # Finite difference methods
using LinearAlgebra      # Linear algebra operations
using CairoMakie              # Plotting

System Parameters and Hamiltonian

The KPO system is governed by the Hamiltonian: $H = -p_1 a^\dagger a + K (a^\dagger)^2 a^2 - G (a^\dagger a^\dagger + a a)$

where:

\[p_1\]
is the parameter we want to estimate (detuning)
\[K\]
is the Kerr nonlinearity
\[G\]
is the parametric drive strength
\[\gamma\]
is the decay rate

function final_state(p)
    G, K, γ = 0.002, 0.001, 0.01

    N = 20 # cutoff of the Hilbert space dimension
    a = destroy(N) # annihilation operator

    coef(p,t) = - p[1]
    H = QobjEvo(a' * a , coef) + K * a' * a' * a * a - G * (a' * a' + a * a)
    c_ops = [sqrt(γ)*a]
    ψ0 = fock(N, 0) # initial state

    tlist = range(0, 2000, 100)

    sol = mesolve(H, ψ0, tlist, c_ops; params = p, progress_bar = Val(false), saveat = [tlist[end]])
    return sol.states[end].data
end

final_state (generic function with 1 method)

function final_state(p, t)
    G, K, γ = 0.002, 0.001, 0.01

    N = 20 # cutoff of the Hilbert space dimension
    a = destroy(N) # annihilation operator

    coef(p,t) = - p[1]
    H = QobjEvo(a' * a , coef) + K * a' * a' * a * a - G * (a' * a' + a * a)
    c_ops = [sqrt(γ)*a]
    ψ0 = fock(N, 0) # initial state

    tlist = range(0, 2000, 100)

    sol = mesolve(H, ψ0, tlist, c_ops; params = p, progress_bar = Val(false), saveat = [t])
    return sol.states[end].data
end

final_state (generic function with 2 methods)

Quantum Fisher Information Calculation

The QFI is computed using the symmetric logarithmic derivative (SLD). For a density matrix $\rho(\theta)$ parametrized by $\theta$:

\[F_Q = \text{Tr}[\partial_\theta \rho \cdot L]\]

where $L$ is the SLD satisfying: $\partial_\theta \rho = \frac{1}{2}(\rho L + L \rho)$

function compute_fisher_information(ρ, dρ)
    # Add small regularization to avoid numerical issues with zero eigenvalues
    reg = 1e-12 * I
    ρ_reg = ρ + reg

    # Solve for the symmetric logarithmic derivative L
    # dρ = (1/2)(ρL + Lρ)
    # This is a Sylvester equation: ρL + Lρ = 2*dρ
    L = sylvester(ρ_reg, ρ_reg, -2*dρ)

    # Fisher information F = Tr(dρ * L)
    F = real(tr(dρ * L))

    return F
end

compute_fisher_information (generic function with 1 method)

Automatic Differentiation Setup

We use finite differences through DifferentiationInterface.jl to compute the derivative of the quantum state with respect to the parameter. This is a key step that enables efficient QFI computation without manual derivative calculations.

# Test the system
final_state([0], 100)

# Define state function for automatic differentiation
state(p) = final_state(p, 2000)

# Compute both the state and its derivative
ρ, dρ = DifferentiationInterface.value_and_jacobian(state, AutoFiniteDiff(), [0.0])

# Reshape the derivative back to matrix form
dρ = QuantumToolbox.vec2mat(vec(dρ))

# Compute QFI at final time
qfi_final = compute_fisher_information(ρ, dρ)
println("QFI at final time: ", qfi_final)

QFI at final time: 19795.821946511165

Time Evolution of Quantum Fisher Information

Now we compute how the QFI evolves over time to understand the optimal measurement time for parameter estimation:

ts = range(0, 2000, 100)

QFI_t = map(ts) do t
    state(p) = final_state(p, t)
    ρ, dρ = DifferentiationInterface.value_and_jacobian(state, AutoFiniteDiff(), [0.0])
    dρ = QuantumToolbox.vec2mat(vec(dρ))
    compute_fisher_information(ρ, dρ)
end

println("QFI computed for ", length(ts), " time points")

QFI computed for 100 time points

Visualization

Plot the time evolution of the Quantum Fisher Information:

fig = Figure(size=(900, 350))
ax = Axis(
 fig[1,1],
 xlabel = "Time",
 ylabel = "Quantum Fisher Information"
)
lines!(ax, ts, QFI_t)
fig

Version Information

QuantumToolbox.versioninfo()


 QuantumToolbox.jl: Quantum Toolbox in Julia
≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
Copyright © QuTiP team 2022 and later.
Current admin team:
    Alberto Mercurio and Yi-Te Huang

Package information:
====================================
Julia              Ver. 1.11.5
QuantumToolbox     Ver. 0.32.1
SciMLOperators     Ver. 0.4.0
LinearSolve        Ver. 3.18.2
OrdinaryDiffEqCore Ver. 1.26.1

System information:
====================================
OS       : Linux (x86_64-linux-gnu)
CPU      : 4 × AMD EPYC 7763 64-Core Processor
Memory   : 15.621 GB
WORD_SIZE: 64
LIBM     : libopenlibm
LLVM     : libLLVM-16.0.6 (ORCJIT, znver3)
BLAS     : libopenblas64_.so (ilp64)
Threads  : 1 (on 4 virtual cores)

using Pkg
Pkg.status()

Status `~/work/SystemLevelHamiltonian.jl/SystemLevelHamiltonian.jl/docs/Project.toml`
⌃ [13f3f980] CairoMakie v0.13.10
⌅ [a0c0ee7d] DifferentiationInterface v0.6.54
  [e30172f5] Documenter v1.13.0
  [35a29f4d] DocumenterTools v0.1.20
  [6a86dc24] FiniteDiff v2.27.0
⌃ [e9467ef8] GLMakie v0.11.11
  [16fef848] LiveServer v1.5.0
⌅ [961ee093] ModelingToolkit v9.80.6
  [1dea7af3] OrdinaryDiffEq v6.98.0
⌃ [35bcea6d] QuantumCumulants v0.3.7
  [6c2fb7c5] QuantumToolbox v0.32.1
  [1ed8b502] SciMLSensitivity v7.86.1
  [a0b4ee4f] SystemLevelHamiltonian v1.0.0-DEV `~/work/SystemLevelHamiltonian.jl/SystemLevelHamiltonian.jl`
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`