Add quantum fisher information for a driven-dissipative Kerr parametric oscillator using auto-diff

Project.toml

@@ -3,3 +3,6 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 HierarchicalEOM = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
 QuartoNotebookRunner = "4c0109c6-14e9-4c88-93f0-2b974d3468f4"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"

QuantumToolbox.jl/time_evolution/Quantum_Fisher_Information_with_automatic_differentiation.qmd

@@ -0,0 +1,157 @@
+---
+title: Quantum Fisher Information
+author: Orjan Ameye
+date: last-modified
+---
+
+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:
+
+```{julia}
+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
+
+```{julia}
+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
+```
+
+```{julia}
+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
+```
+
+## 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)$
+
+```{julia}
+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
+```
+
+## 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.
+
+```{julia}
+# 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)
+```
+
+## Time Evolution of Quantum Fisher Information
+
+Now we compute how the QFI evolves over time to understand the optimal measurement time for parameter estimation:
+
+```{julia}
+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")
+```
+
+## Visualization
+
+Plot the time evolution of the Quantum Fisher Information:
+
+```{julia}
+fig = Figure(size=(900, 350))
+ax = Axis(
+ fig[1,1],
+ xlabel = "Time",
+ ylabel = "Quantum Fisher Information"
+)
+lines!(ax, ts, QFI_t)
+fig
+```
+
+## Version Information
+```{julia}
+QuantumToolbox.versioninfo()
+```
+
+```{julia}
+using Pkg
+Pkg.status()
+```