From f30a3a2d4997a9f4a13c406ff21dc4e882d65463 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Thu, 25 Jun 2026 19:06:08 +0200 Subject: [PATCH] feat: susceptibility --- docs/src/api.md | 22 ++ ...lated-emission__PRL2019_123-123604_fig4.jl | 2 +- ...-1_kerr-parametric-oscillator__response.jl | 142 +++++++++ src/QuantumInputOutput.jl | 7 + src/response.jl | 291 ++++++++++++++++++ test/runtests.jl | 1 + test/test_response.jl | 110 +++++++ 7 files changed, 574 insertions(+), 1 deletion(-) create mode 100644 examples/10-1_kerr-parametric-oscillator__response.jl create mode 100644 src/response.jl create mode 100644 test/test_response.jl diff --git a/docs/src/api.md b/docs/src/api.md index 1b8eca1..392486b 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -94,6 +94,28 @@ solve_mode_evolution_symmetric correlation_matrix ``` +## [Response: S-parameters, susceptibility, spectra](@id API: Response) + +```@docs +output_field +``` + +```@docs +scattering_parameter +``` + +```@docs +susceptibility +``` + +```@docs +power_spectrum +``` + +```@docs +liouvillian_resolvent +``` + ## [Utilities](@id API: Utilities) ```@docs diff --git a/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl b/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl index e786f37..b0978cc 100644 --- a/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl +++ b/examples/01-2_stimulated-emission__PRL2019_123-123604_fig4.jl @@ -113,7 +113,7 @@ p = plot(t_, popu_u_n1; color = :red, label = L"| 1 \rangle_u") plot!(p, t_, popu_e; color = :pink, ls = :dash, label = L"| e \rangle") plot!(p, t_, popu_v_n1; color = :blue, ls = :dash, label = L"| 1 \rangle_v") plot!(p, t_, popu_v_n2; color = :green, label = L"| 2 \rangle_v") -plot!(p, t_, I_out_int; color = :black, ls = :dot, label = L"\\int I_{out} dt") +plot!(p, t_, I_out_int; color = :black, ls = :dot, label = L"\int I_{out} dt") plot!( p; xlims = (0, 4), diff --git a/examples/10-1_kerr-parametric-oscillator__response.jl b/examples/10-1_kerr-parametric-oscillator__response.jl new file mode 100644 index 0000000..c8c9b2c --- /dev/null +++ b/examples/10-1_kerr-parametric-oscillator__response.jl @@ -0,0 +1,142 @@ +# # Driven response of a Kerr parametric oscillator + +# A Kerr parametric oscillator (KPO) is a single bosonic mode with a Kerr nonlinearity +# that is pumped by a two-photon (parametric) drive. Written in the frame rotating at half +# the pump frequency, its Hamiltonian is time independent, +# +# ```math +# H = \Delta\, a^\dagger a + K\, a^\dagger a^\dagger a a + \frac{G}{2}\left(a^\dagger a^\dagger + a a\right), +# ``` +# +# with detuning ``\Delta``, Kerr coefficient ``K`` and parametric pump amplitude ``G``. We +# couple the mode to a transmission line through two ports with rates ``\kappa_1`` (input) +# and ``\kappa_2`` (output), so the device is a **two-port** scatterer. +# +# This example shows how to measure such a device the way an experimentalist would, with a +# weak probe tone, using three tools from `QuantumInputOutput`: +# +# * [`scattering_parameter`](@ref) — the transmission ``S_{21}(\omega)`` and reflection ``S_{11}(\omega)``, +# * [`power_spectrum`](@ref) — the output emission and squeezing spectra, +# * [`susceptibility`](@ref) — the underlying Kubo linear-response engine. +# +# The probe is treated only to linear order (the weak-tone regime), but the system itself is +# never linearised: the full nonlinear master equation enters through the Liouvillian. No +# steady-state Jacobian or fluctuation expansion is needed. + +using QuantumInputOutput +using SecondQuantizedAlgebra +using QuantumOptics +using Plots +using LaTeXStrings + +# ## Building the two-port KPO as an SLH component + +# The mode lives in a single Fock space. A two-port cavity has a two-component Lindblad +# vector ``L = (\sqrt{\kappa_1}\,a,\ \sqrt{\kappa_2}\,a)`` and an identity scattering matrix. + +h = FockSpace(:kpo) +a = Destroy(h, :a) + +@variables Δ::Real K::Real G::Real κ1::Real κ2::Real + +H = Δ*a'*a + K*a'*a'*a*a + (G/2)*(a'*a' + a*a) +Gkpo = SLH([1 0; 0 1], [√(κ1)*a, √(κ2)*a], H) +nothing # hide + +# The numeric basis and a working point below the parametric oscillation threshold +# ``G_\mathrm{th} = (\kappa_1+\kappa_2)/2``. + +b = FockBasis(20) +κ1_, κ2_ = 0.5, 0.5 +params(G_; Δ_ = 0.0, K_ = 0.1) = Dict(Δ => Δ_, K => K_, G => G_, κ1 => κ1_, κ2 => κ2_) +nothing # hide + +# ## (A) Steady state and output field + +# The pumped steady state is the working point we linearise the probe around, and it is the +# one expensive, method-sensitive ingredient shared by all three tools below: you compute it +# once with whichever solver suits your system and pass it in. Here we translate the network to +# numeric operators with [`translate_qo`](@ref) and take the exact null vector of the +# Liouvillian with `QuantumOptics.steadystate`. The output field operator of a port follows the +# input–output relation ``b_{\mathrm{out},k} = (S\,b_{\mathrm{in}})_k - L_k``. + +function steady(p) + Hn, Jn = translate_qo(Gkpo, b; parameter = p) + return steadystate.eigenvector(Hn, collect(Jn)) +end +ρ_passive = steady(params(0.0)) +ρ_mid = steady(params(0.3)) +ρ_near = steady(params(0.45)) +b_out = output_field(Gkpo, 2) # symbolic output operator of port 2 +nothing # hide + +# ## (B) Transmission ``S_{21}(\omega)`` and parametric gain + +# We sweep the probe detuning ``\omega``. With the pump off (``G=0``) the device is a passive +# two-port cavity and ``|S_{21}| \le 1``. As ``G`` approaches threshold the parametric drive +# adds gain and ``|S_{21}|`` exceeds one: the KPO acts as a phase-insensitive amplifier. + +ω = collect(range(-1.5, 1.5; length = 401)) +S21_passive = scattering_parameter(Gkpo, b, ρ_passive; in_port = 1, out_port = 2, omega = ω, parameter = params(0.0)) +S21_mid = scattering_parameter(Gkpo, b, ρ_mid; in_port = 1, out_port = 2, omega = ω, parameter = params(0.3)) +S21_near = scattering_parameter(Gkpo, b, ρ_near; in_port = 1, out_port = 2, omega = ω, parameter = params(0.45)) + +p1 = plot(ω, abs.(S21_passive); label = "G = 0 (passive)", color = :black) +plot!(p1, ω, abs.(S21_mid); label = "G = 0.3", color = :blue) +plot!(p1, ω, abs.(S21_near); label = "G = 0.45 (near threshold 0.5)", color = :red) +hline!(p1, [1.0]; ls = :dot, color = :gray, label = "") +plot!(p1; xlabel = L"probe detuning $\omega$", ylabel = L"|S_{21}(\omega)|", + title = "KPO transmission: parametric gain", legend = :topright, size = (650, 350)) +p1 + +# For a lossless passive two-port, transmission and reflection conserve energy, +# ``|S_{11}|^2 + |S_{21}|^2 = 1``. This is a useful sanity check on the working point. + +S11_passive = scattering_parameter(Gkpo, b, ρ_passive; in_port = 1, out_port = 1, omega = ω, parameter = params(0.0)) +@info "max energy deviation (passive): " maximum(abs.(abs2.(S11_passive) .+ abs2.(S21_passive) .- 1)) + +# ## (C) Output emission and squeezing spectra + +# With the pump on, the KPO emits parametric fluorescence even with no probe. The emission +# spectrum of the output port is the Fourier transform of ``\langle b_\mathrm{out}^\dagger(\tau)\,b_\mathrm{out}(0)\rangle`` +# and has a zero vacuum floor. Passing a `quadrature` angle instead returns the vacuum-normalised +# homodyne (squeezing) spectrum, whose shot-noise floor is one: the squeezed quadrature drops +# below it while the orthogonal one is anti-squeezed above it. + +S_emit = power_spectrum(Gkpo, b, ρ_mid; port = 2, omega = ω, parameter = params(0.3)) +S_sqz = power_spectrum(Gkpo, b, ρ_mid; port = 2, omega = ω, parameter = params(0.3), quadrature = π/4) +S_anti = power_spectrum(Gkpo, b, ρ_mid; port = 2, omega = ω, parameter = params(0.3), quadrature = 3π/4) + +p2 = plot(ω, S_emit; label = "emission", color = :black) +plot!(p2, ω, S_sqz; label = L"squeezed quadrature $\theta=\pi/4$", color = :blue) +plot!(p2, ω, S_anti; label = L"anti-squeezed $\theta=3\pi/4$", color = :red) +hline!(p2, [1.0]; ls = :dot, color = :gray, label = "vacuum / shot-noise floor") +plot!(p2; xlabel = L"frequency $\omega$", ylabel = "output spectrum", + title = "KPO output spectra", legend = :topright, size = (650, 380)) +p2 + +# ## (B, general) The susceptibility engine + +# [`scattering_parameter`](@ref) is a thin wrapper over [`susceptibility`](@ref), the Kubo response +# ``\chi_{AB}(\omega) = i\,\mathrm{Tr}[A\,(\mathcal{L}+i\omega)^{-1}[B,\rho_{ss}]]``. The +# anomalous (idler) response of a parametric system, the response of ``a`` to a drive that +# couples through ``a`` rather than ``a^\dagger``, is nonzero and is what makes the gain +# phase sensitive. It is one call away: + +χ_normal = susceptibility(Gkpo, b, ρ_mid, a, a', ω; parameter = params(0.3)) # ⟨a⟩ response to a† drive +χ_anomalous = susceptibility(Gkpo, b, ρ_mid, a, a, ω; parameter = params(0.3)) # ⟨a⟩ response to a drive (idler) + +p3 = plot(ω, abs.(χ_normal); label = L"\chi_{a,a^\dagger} (normal)", color = :blue) +plot!(p3, ω, abs.(χ_anomalous); label = L"\chi_{a,a} (anomalous/idler)", color = :red) +plot!(p3; xlabel = L"\omega", ylabel = L"|\chi(\omega)|", + title = "Normal vs anomalous response", legend = :topright, size = (650, 350)) +p3 + +# ## Package versions + +using InteractiveUtils +versioninfo() + +using Pkg +Pkg.status(["QuantumInputOutput", "SecondQuantizedAlgebra", "QuantumOptics", "Plots"]; + mode = PKGMODE_MANIFEST) diff --git a/src/QuantumInputOutput.jl b/src/QuantumInputOutput.jl index 9ee3a8f..2314432 100644 --- a/src/QuantumInputOutput.jl +++ b/src/QuantumInputOutput.jl @@ -42,6 +42,12 @@ export SLH, solve_mode_evolution_symmetric, # Correlations correlation_matrix, + # Response: S-parameters, susceptibility, spectra + output_field, + liouvillian_resolvent, + susceptibility, + scattering_parameter, + power_spectrum, # Operators substitute_operators @@ -51,5 +57,6 @@ include("utils.jl") include("pulses.jl") include("correlations.jl") include("interaction_picture.jl") +include("response.jl") end diff --git a/src/response.jl b/src/response.jl new file mode 100644 index 0000000..db05c1a --- /dev/null +++ b/src/response.jl @@ -0,0 +1,291 @@ +################################################################ +### linear & nonlinear response: steady state, S-parameters, ### +### susceptibility and output spectra ### +################################################################ +# +# These tools compute the *driven response* of an SLH network around a +# (pumped) steady state, the way a network/spectrum analyser measures it: +# a weak probe tone is sent into a port and the coherent field leaving +# another port is read out. The probe response is obtained exactly from the +# Liouvillian, without linearising the system Hamiltonian, so it stays valid +# for genuinely nonlinear systems such as a Kerr parametric oscillator. Only +# the probe is treated to linear order (the experimental weak-tone regime). +# +# All three observables reduce to the same primitive: the resolvent of the +# Liouvillian, `(ℒ + μ)⁻¹`, applied to a traceless operator. They share one +# Hessenberg factorisation (`LiouvillianResolvent`) so a frequency sweep costs +# one cheap shifted solve per point. + +using QuantumOpticsBase: AbstractOperator, Operator, identityoperator, expect, liouvillian +using LinearAlgebra: hessenberg, tr, I + +# `dagger` is imported from QuantumOpticsBase by the top-level module. + +# ────────────────────────────────────────────── +# A. Output field +# ────────────────────────────────────────────── + +""" + output_field(G::SLH, k; input=0) + +Return the **output-field operator** of port `k` of an SLH network, i.e. the +operator whose expectation value is the travelling field leaving that port. + +Following the input–output boundary relation ``b_{\\mathrm{out},k} = +(S\\,b_{\\mathrm{in}})_k - L_k``, the operator part is `-lindblad(G)[k]`. With no +field injected (`input=0`) this is just `-L_k`; passing a vector `input` of +coherent input amplitudes per port adds the directly-scattered classical part +``(S\\,\\alpha)_k``. + +This is the canonical observable behind `scattering_parameter` and `power_spectrum`. +""" +function output_field(G::SLH, k::Integer; input = 0) + L = lindblad(G)[k] + if input == 0 || (input isa AbstractVector && all(iszero, input)) + return -L + end + α = input isa AbstractVector ? input : error("`input` must be a vector of port amplitudes") + S = scattering(G) + offset = sum(S[k, j] * α[j] for j in axes(S, 2)) + return offset - L +end + +# ────────────────────────────────────────────── +# Liouvillian resolvent (shared engine) +# ────────────────────────────────────────────── + +""" + LiouvillianResolvent + +Precomputed factorisation for evaluating ``(\\mathcal{L} + \\mu)^{-1} X`` at many +shifts `μ`, where ``\\mathcal{L}`` is the Liouvillian and `X` is a traceless +operator. The zero mode of ``\\mathcal{L}`` (the steady state) is deflated so the +resolvent is regular at every real frequency, including ``\\omega = 0``; this is +exact because all response right-hand sides used here are traceless. + +Build with [`liouvillian_resolvent`](@ref); apply with +[`susceptibility`](@ref), [`scattering_parameter`](@ref) or [`power_spectrum`](@ref). +""" +struct LiouvillianResolvent{BT,OT,FT} + basis::BT + ρ_ss::OT + F::FT + d::Int +end + +""" + liouvillian_resolvent(H, J, ρ_ss; deflation=1.0) + +Construct a [`LiouvillianResolvent`](@ref) from a numeric Hamiltonian `H`, +collapse operators `J` and a precomputed steady state `ρ_ss`. Builds a Hessenberg +factorisation of the deflated Liouvillian for fast shifted solves; the zero mode +is deflated using `ρ_ss` as the right null vector. + +Compute `ρ_ss` yourself with whichever solver suits your system, e.g. +`QuantumOptics.steadystate.eigenvector(H, collect(J))`, and reuse the resulting +resolvent across [`susceptibility`](@ref), [`scattering_parameter`](@ref) and +[`power_spectrum`](@ref). +""" +function liouvillian_resolvent( + H::AbstractOperator, + J, + ρ_ss::AbstractOperator; + deflation = 1.0, +) + b = H.basis_l + d = length(b) + Jvec = collect(J) # liouvillian needs a plain Vector + ℒ = liouvillian(H, Jvec) + ρss = ρ_ss + vecI = reshape(Matrix(identityoperator(b).data), d * d) # ⟨𝟙| : the trace functional + vecρ = reshape(Matrix(ρss.data), d * d) # |ρ_ss⟩ : the right null vector + F = hessenberg(Matrix(ℒ.data) + deflation * (vecρ * vecI')) + return LiouvillianResolvent(b, ρss, F, d) +end + +# Solve (ℒ + μ) X = rhs for the operator X. `rhs` must be traceless. +function _resolvent_solve(R::LiouvillianResolvent, μ::Number, rhs::AbstractOperator) + x = (R.F + μ * I) \ reshape(Matrix(rhs.data), R.d * R.d) + return Operator(R.basis, R.basis, reshape(x, R.d, R.d)) +end + +_commutator(A, B) = A * B - B * A + +# ────────────────────────────────────────────── +# B. Susceptibility (Kubo) and S-parameters +# ────────────────────────────────────────────── + +""" + susceptibility(R::LiouvillianResolvent, A, B, ω) + susceptibility(G::SLH, b, ρ_ss, A, B, ω; parameter=Dict(), operators=Dict(), kwargs...) + +Kubo linear-response susceptibility + +```math +\\chi_{AB}(\\omega) = i\\,\\mathrm{Tr}\\!\\big[A\\,(\\mathcal{L} + i\\omega)^{-1}\\,[B, \\rho_{ss}]\\big], +``` + +the response of ``\\langle A \\rangle`` at frequency `ω` to a weak perturbation +``H_1(t) = f(t)\\,B + \\text{h.c.}`` with ``f(t) = e^{-i\\omega t}``. The system is +**not** linearised: ``\\mathcal{L}`` is the full (nonlinear) Liouvillian and only +the perturbation is taken to first order. `ω` may be a scalar or a vector. + +Pass numeric operators `A`, `B` and a prebuilt resolvent `R`, or symbolic `A`, +`B` together with an `SLH` network `G`, basis `b` and a precomputed steady state +`ρ_ss` (from `QuantumOptics.steadystate` or any solver) to translate and build the +resolvent in one call. This is the general engine; for the standard transmission +or reflection coefficient use [`scattering_parameter`](@ref). +""" +function susceptibility(R::LiouvillianResolvent, A::AbstractOperator, B::AbstractOperator, ω::Real) + X = _resolvent_solve(R, im * ω, _commutator(B, R.ρ_ss)) + return im * tr(A.data * X.data) +end + +susceptibility(R::LiouvillianResolvent, A, B, ω::AbstractVector) = + [susceptibility(R, A, B, ωk) for ωk in ω] + +function susceptibility( + G::SLH, + b, + ρ_ss::AbstractOperator, + A, + B, + ω; + parameter = Dict(), + operators = Dict(), + kwargs..., +) + H, J = translate_qo(G, b; parameter, operators) + A_ = translate_qo(A, b; parameter, operators) + B_ = translate_qo(B, b; parameter, operators) + R = liouvillian_resolvent(H, J, ρ_ss; kwargs...) + return susceptibility(R, A_, B_, ω) +end + +""" + scattering_parameter(G::SLH, b, ρ_ss; in_port=1, out_port=2, omega, parameter=Dict(), operators=Dict(), kwargs...) + +Scattering parameter ``S_{\\text{out\\_port},\\text{in\\_port}}(\\omega)`` of an SLH +network: the coherent field leaving `out_port` per unit weak probe sent into +`in_port`, as a function of probe detuning `omega` (in the rotating frame of the +pump). With `out_port ≠ in_port` this is a transmission coefficient (``S_{21}``); +with `out_port == in_port`, a reflection coefficient (``S_{11}``). + +`ρ_ss` is the precomputed pumped steady state (from `QuantumOptics.steadystate` or any +solver) at the same working point; compute it once and reuse it across +`scattering_parameter`/`power_spectrum`/`susceptibility`. + +```math +S_{\\text{out},\\text{in}}(\\omega) = S_{\\text{out},\\text{in}} + + \\mathrm{Tr}\\!\\big[L_{\\text{out}}\\,(\\mathcal{L} + i\\omega)^{-1}\\,[L_{\\text{in}}^\\dagger, \\rho_{ss}]\\big], +``` + +where ``L_k`` is the port-`k` coupling operator (`lindblad(G)[k]`) and the first +term is the direct scattering matrix element `scattering(G)[out_port, in_port]`. +The Liouvillian is built from **all** ports of `G`, so internal loss ports +correctly reduce the transmission. Because the probe enters only to linear order, +the result is the exact weak-tone response of the full nonlinear network and +satisfies ``|S_{11}|^2 + |S_{21}|^2 = 1`` for a lossless two-port. + +`kwargs` (e.g. `method`, `rho_ss`, `deflation`) are forwarded to +[`liouvillian_resolvent`](@ref). +""" +function scattering_parameter( + G::SLH, + b, + ρ_ss::AbstractOperator; + in_port::Integer = 1, + out_port::Integer = 2, + omega, + parameter = Dict(), + operators = Dict(), + kwargs..., +) + H, J = translate_qo(G, b; parameter, operators) + L_in = J[in_port] + L_out = J[out_port] + R = liouvillian_resolvent(H, J, ρ_ss; kwargs...) + direct = scattering(G)[out_port, in_port] + rhs = _commutator(dagger(L_in), R.ρ_ss) # [L_in†, ρ_ss], traceless + return [direct + tr(L_out.data * _resolvent_solve(R, im * ωk, rhs).data) for ωk in omega] +end + +# ────────────────────────────────────────────── +# C. Output power / squeezing spectrum +# ────────────────────────────────────────────── + +""" + power_spectrum(G::SLH, b, ρ_ss; port=1, omega, quadrature=nothing, + parameter=Dict(), operators=Dict(), subtract_mean=true, kwargs...) + +Steady-state output spectrum at `port` of an SLH network around the precomputed +steady state `ρ_ss` X. Two physically +distinct spectra are available. + +**Emission spectrum** (default, `quadrature=nothing`). With ``b_{\\mathrm{out}} = +-L`` the port output field and ``\\bar b_{\\mathrm{out}} = b_{\\mathrm{out}} - +\\langle b_{\\mathrm{out}}\\rangle``, + +```math +S(\\omega) = -2\\,\\mathrm{Re}\\,\\mathrm{Tr}\\!\\big[\\bar b_{\\mathrm{out}}^\\dagger\\,(\\mathcal{L} - i\\omega)^{-1}\\,(\\bar b_{\\mathrm{out}}\\,\\rho_{ss})\\big], +``` + +the Fourier transform of ``\\langle \\bar b_{\\mathrm{out}}^\\dagger(\\tau)\\,\\bar +b_{\\mathrm{out}}(0)\\rangle`` (convention ``\\int e^{-i\\omega\\tau}\\,d\\tau``, +matching `QuantumOptics.timecorrelations.spectrum`). This is the +normally-ordered parametric-fluorescence / emission spectrum; its vacuum level +is zero. + +**Squeezing spectrum** (`quadrature=θ`). The vacuum-normalised homodyne spectrum +of the output quadrature ``X_\\theta = (e^{-i\\theta} b_{\\mathrm{out}} + +e^{i\\theta} b_{\\mathrm{out}}^\\dagger)/\\sqrt{2}``, + +```math +S_\\theta(\\omega) = 1 + 2\\,\\mathrm{Re}\\Big[ + e^{-2i\\theta}\\big(C_{LL}(\\omega)+C_{LL}(-\\omega)\\big) + + C_{L^\\dagger L}(\\omega)+C_{L^\\dagger L}(-\\omega)\\Big], +\\quad +C_{PQ}(\\omega) = -\\mathrm{Tr}\\!\\big[\\delta P\\,(\\mathcal{L}+i\\omega)^{-1}(\\delta Q\\,\\rho_{ss})\\big], +``` + +with ``L`` the port coupling operator and ``\\delta O = O - \\langle O\\rangle``. +The shot-noise floor is ``1``; a squeezed quadrature dips below it. For a +degenerate parametric amplifier this reproduces the Collett–Gardiner result +``S_\\mp(\\omega) = 1 \\mp 2\\kappa G/((\\kappa/2 \\pm G)^2 + \\omega^2)`` exactly. + +Set `subtract_mean=false` to keep the coherent part. `kwargs` are forwarded to +[`liouvillian_resolvent`](@ref). +""" +function power_spectrum( + G::SLH, + b, + ρ_ss::AbstractOperator; + port::Integer = 1, + omega, + quadrature = nothing, + parameter = Dict(), + operators = Dict(), + subtract_mean = true, + kwargs..., +) + H, J = translate_qo(G, b; parameter, operators) + R = liouvillian_resolvent(H, J, ρ_ss; kwargs...) + Lp = J[port] + if quadrature === nothing + Abar = subtract_mean ? Lp - expect(Lp, R.ρ_ss) * identityoperator(b) : Lp + Adag = dagger(Abar) + rhs = Abar * R.ρ_ss # traceless (⟨Ā⟩ = 0) + return [-2 * real(tr(Adag.data * _resolvent_solve(R, -im * ωk, rhs).data)) for ωk in omega] + end + θ = quadrature + δL = subtract_mean ? Lp - expect(Lp, R.ρ_ss) * identityoperator(b) : Lp + δLd = dagger(δL) + # one-sided spectrum C_PQ(ω) = ∫₀^∞ e^{iωτ}⟨δP(τ)δQ(0)⟩dτ = -Tr[δP (ℒ+iω)⁻¹ (δQ ρ_ss)] + C(P, Q, ω) = -tr(P.data * _resolvent_solve(R, im * ω, Q * R.ρ_ss).data) + return [ + 1 + 2 * real( + exp(-2im * θ) * (C(δL, δL, ωk) + C(δL, δL, -ωk)) + + (C(δLd, δL, ωk) + C(δLd, δL, -ωk)), + ) for ωk in omega + ] +end diff --git a/test/runtests.jl b/test/runtests.jl index 2f6c82f..6a0b252 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ names = [ "test_feedback.jl", "test_translate.jl", "test_correlations.jl", + "test_response.jl", "test_example_cavity_scattering.jl", "test_compare_example_05_1_05_2.jl", "test_interaction_picture.jl", diff --git a/test/test_response.jl b/test/test_response.jl new file mode 100644 index 0000000..6438c23 --- /dev/null +++ b/test/test_response.jl @@ -0,0 +1,110 @@ +using QuantumInputOutput +using SecondQuantizedAlgebra +using QuantumOptics +using QuantumOpticsBase: dagger +using LinearAlgebra +using Test + +# Symbolic two-port Kerr parametric oscillator (KPO), rotating frame: +# H = -Δ a†a + K a†a†aa + (G/2)(a†a† + aa), ports L = (√κ1 a, √κ2 a) +h = FockSpace(:kpo) +a = Destroy(h, :a) +@variables Δ::Real K::Real G::Real κ1::Real κ2::Real +H = -Δ*a'*a + K*a'*a'*a*a + (G/2)*(a'*a' + a*a) +Gkpo = SLH([1 0; 0 1], [√(κ1)*a, √(κ2)*a], H) +b = FockBasis(16) + +# the user computes the steady state and passes it in (the package no longer wraps a solver) +function steady(p) + Hn, Jn = QuantumInputOutput.translate_qo(Gkpo, b; parameter = p) + return steadystate.eigenvector(Hn, collect(Jn)) +end + +@testset "output_field" begin + # b_out,k = -L_k with no input + Lnum = QuantumInputOutput.translate_qo(√(κ2)*a, b; parameter = Dict(κ2 => 0.6)) + of = QuantumInputOutput.translate_qo(output_field(Gkpo, 2), b; parameter = Dict(κ2 => 0.6)) + @test isapprox(Matrix(of.data), -Matrix(Lnum.data); atol = 1e-12) + # coherent offset on the same port adds (S α)_k = α_k for S = I + of2 = QuantumInputOutput.translate_qo( + output_field(Gkpo, 2; input = [0.0, 2.0]), b; parameter = Dict(κ2 => 0.6)) + @test isapprox(Matrix(of2.data), 2.0*Matrix(one(b).data) - Matrix(Lnum.data); atol = 1e-12) +end + +@testset "scattering_parameter: linear cavity vs analytic Lorentzian + energy conservation" begin + Δ_, κ1_, κ2_ = 0.7, 0.4, 0.6 + p = Dict(Δ => Δ_, K => 0.0, G => 0.0, κ1 => κ1_, κ2 => κ2_) + ω = collect(-2.0:0.1:2.0) + ρ = steady(p) + S21 = scattering_parameter(Gkpo, b, ρ; in_port = 1, out_port = 2, omega = ω, parameter = p) + S11 = scattering_parameter(Gkpo, b, ρ; in_port = 1, out_port = 1, omega = ω, parameter = p) + # Gardiner convention a_out = a_in - √κ a, with H = -Δ a†a + κ = κ1_ + κ2_ + ana(δ) = -sqrt(κ1_*κ2_) / (κ/2 + im*(-Δ_ - δ)) + @test all(isapprox.(S21, ana.(ω); atol = 1e-6)) + # lossless two-port conserves energy at every frequency + @test maximum(abs.(abs2.(S11) .+ abs2.(S21) .- 1)) < 1e-8 + # resonance peak (critical for κ1=κ2 would be 1; here 2√(κ1κ2)/κ) + @test maximum(abs.(S21)) ≈ 2*sqrt(κ1_*κ2_)/κ atol = 1e-4 +end + +@testset "scattering_parameter: parametric gain above passive bound" begin + # below threshold G < κ/2 = 0.5; expect |S21| > passive max of 1 + p = Dict(Δ => 0.0, K => 0.0, G => 0.4, κ1 => 0.5, κ2 => 0.5) + ω = collect(-1.0:0.02:1.0) + ρ = steady(p) + S21 = scattering_parameter(Gkpo, b, ρ; in_port = 1, out_port = 2, omega = ω, parameter = p) + @test maximum(abs.(S21)) > 1.0 +end + +@testset "power_spectrum: emission ≥ 0 and matches correlation integral" begin + p = Dict(Δ => 0.0, K => 0.0, G => 0.35, κ1 => 1.0, κ2 => 1e-9) + ω = collect(-2.0:0.25:2.0) + Hn, Jn = QuantumInputOutput.translate_qo(Gkpo, b; parameter = p) + Jn = collect(Jn) + ρ = steadystate.eigenvector(Hn, Jn) + S = power_spectrum(Gkpo, b, ρ; port = 1, omega = ω, parameter = p) + @test all(S .>= -1e-8) + # compare to direct numerical FT of ⟨ā†(τ)ā(0)⟩ from QuantumOptics + A = Jn[1] + amean = tr(A.data*ρ.data) + Ā = A - amean*one(b) + τ = collect(0.0:0.004:80.0) + g = timecorrelations.correlation(τ, ρ, Hn, Jn, dagger(Ā), Ā) + Squad(w) = begin + f = @. exp(-im*w*τ) * g + 2*real(sum((f[1:end-1] .+ f[2:end]) ./ 2 .* diff(τ))) + end + @test all(isapprox(S[i], Squad(ω[i]); rtol = 1e-3, atol = 1e-5) for i in eachindex(ω)) +end + +@testset "power_spectrum: squeezing matches analytic DPA + below vacuum" begin + # Degenerate parametric amplifier (K=0, single port): the vacuum-normalised output + # squeezing spectrum has the closed form (Collett–Gardiner) + # S_∓(ω) = 1 ∓ 2κG / ((κ/2 ± G)² + ω²), κ = κ1 + κ2. + κ_, G_ = 1.0, 0.35 + p = Dict(Δ => 0.0, K => 0.0, G => G_, κ1 => κ_, κ2 => 1e-9) + ω = collect(-2.5:0.25:2.5) + ρ = steady(p) + sqz = power_spectrum(Gkpo, b, ρ; port = 1, omega = ω, parameter = p, quadrature = π/4) + anti = power_spectrum(Gkpo, b, ρ; port = 1, omega = ω, parameter = p, quadrature = 3π/4) + Ssq(δ) = 1 - 2κ_*G_ / ((κ_/2 + G_)^2 + δ^2) + Santi(δ) = 1 + 2κ_*G_ / ((κ_/2 - G_)^2 + δ^2) + # the sharp anti-squeezed peak (~32) is sensitive to the Fock truncation; use a + # relative tolerance (≈0.15% at FockBasis(16), →1e-5 at FockBasis(24)) + @test all(isapprox.(sqz, Ssq.(ω); rtol = 1e-2, atol = 2e-3)) + @test all(isapprox.(anti, Santi.(ω); rtol = 1e-2, atol = 2e-3)) + @test minimum(sqz) < 1.0 # squeezed below the shot-noise floor + @test maximum(anti) > 1.0 # anti-squeezed above it +end + +@testset "susceptibility: reduces to scattering_parameter relation" begin + p = Dict(Δ => 0.2, K => 0.05, G => 0.25, κ1 => 0.5, κ2 => 0.5) + ω = collect(-1.0:0.2:1.0) + # S_{out,in} = scattering[out,in] + √(κ_out κ_in) · χ_{a,a†}/(i) is implicit; here just + # check the engine runs, returns the right length and finite complex numbers + ρ = steady(p) + χ = susceptibility(Gkpo, b, ρ, a, a', ω; parameter = p) + @test length(χ) == length(ω) + @test all(isfinite, abs.(χ)) +end