Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion TODO.md
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@
- [x] feedback reduction
- [x] interaction picture
- [x] general expression for A(t)
- [x] two_time_corr_matrix for Ls const (see example 02-2)
- [x] two_time_corr_matrix for Js const (see example 02-2)

- [x] Formatter/SpellCheck
- [x] JET
Expand Down
12 changes: 6 additions & 6 deletions benchmarks/correlations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ function benchmark_correlations!(SUITE)
G_cas = ▷(G_u, G_c, G_v)

H_sym = hamiltonian(G_cas)
L_sym = jump_operator(G_cas)[1]
J_sym = jump_operator(G_cas)[1]

γ_ = 1.0
σ_pulse = 1 / γ_
Expand All @@ -40,24 +40,24 @@ function benchmark_correlations!(SUITE)
dict_p_t = Dict(g_u => gu_t)

H_QO = to_numeric(H_sym, b; parameter = dict_p, time_parameter = dict_p_t)
L_QO = to_numeric(L_sym, b; parameter = dict_p, time_parameter = dict_p_t)
J_QO = to_numeric(J_sym, b; parameter = dict_p, time_parameter = dict_p_t)

ψ0 = fockstate(bu1, 1) ⊗ fockstate(bc1, 0) ⊗ fockstate(bv1, 0)
_, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [L_QO])
_, ρt = timeevolution.master_dynamic(T, ψ0, H_QO, [J_QO])

au_qo = to_numeric(au, b)
c_qo = to_numeric(c, b)
Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo
Js(t) = gu_t(t) * au_qo + √(γ_) * c_qo

## --- Two-time correlation ---
# Pass the time-dependent operators (`H_QO`, `[L_QO]`) directly to the solver rather
# Pass the time-dependent operators (`H_QO`, `[J_QO]`) directly to the solver rather
# than wrapping them in a `(t, ρ)` closure: since v0.10 `to_numeric` returns a lazy
# `TimeDependentSum`, the direct path lets the solver build its integrator once and is
# markedly faster than re-reading the operator from a function at every step.

SUITE["Correlations"]["two-time"] = BenchmarkGroup()

SUITE["Correlations"]["two-time"]["single photon cavity"] =
@benchmarkable correlation_matrix($T, $ρt, $H_QO, [$L_QO], $Ls)
@benchmarkable correlation_matrix($T, $ρt, $H_QO, [$J_QO], $Js)
return nothing
end
4 changes: 2 additions & 2 deletions benchmarks/interaction_picture.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ function benchmark_interaction_picture!(SUITE)
G_cas = ▷(G_u, G_s, G_v)

H = hamiltonian(G_cas)
L = jump_operator(G_cas)[1]
J = jump_operator(G_cas)[1]
H_uv = hamiltonian(▷(G_u, G_v))
H_int_ = simplify(H - H_uv)

Expand All @@ -82,7 +82,7 @@ function benchmark_interaction_picture!(SUITE)
SUITE["Interaction Picture"]["operator substitution"]["TLS cascade"] =
@benchmarkable begin
simplify(substitute($H_int_, $int_dict))
simplify(substitute($L, $int_dict))
simplify(substitute($J, $int_dict))
end

return nothing
Expand Down
10 changes: 5 additions & 5 deletions benchmarks/slh_algebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -99,20 +99,20 @@ function benchmark_slh_algebra!(SUITE)
G_t = G_R_t ⊞ G_L_t

H_f = hamiltonian(G_t)
L_f = jump_operator(G_t)
J_f = jump_operator(G_t)

t_mid = T[length(T)÷2]

_callable(x) = x isa Union{Function,FunctionWrappers.FunctionWrapper}
Hf = _callable(H_f) ? H_f : (t -> H_f)
L_callables = [_callable(Li) ? Li : (t -> Li) for Li in L_f]
J_callables = [_callable(Ji) ? Ji : (t -> Ji) for Ji in J_f]

SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide H(t)"] =
@benchmarkable $Hf($t_mid)

SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide L(t)"] = @benchmarkable begin
for Li in $L_callables
Li($t_mid)
SUITE["SLH Algebra"]["closure evaluation"]["2-QD waveguide J(t)"] = @benchmarkable begin
for Ji in $J_callables
Ji($t_mid)
end
end

Expand Down
16 changes: 8 additions & 8 deletions benchmarks/translation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ function benchmark_translation!(SUITE)
a = Destroy(h, :a, 1)
σ(i, j) = Transition(h, Symbol("σ"), i, j, 2)

## Derive H and L from cascade
## Derive H and J from cascade
@variables gu_sym::Real γ_sym::Real gv_sym::Real
hu_ = FockSpace(:u)
hc_ = FockSpace(:c)
Expand All @@ -30,7 +30,7 @@ function benchmark_translation!(SUITE)
G_v = SLH(1, gv_sym * av_, 0)
G_cas = ▷(G_u, G_c, G_v)
H_sym = hamiltonian(G_cas)
L_sym = jump_operator(G_cas)[1]
J_sym = jump_operator(G_cas)[1]

## --- Static translation (no time dependence) ---

Expand Down Expand Up @@ -63,7 +63,7 @@ function benchmark_translation!(SUITE)
time_parameter = $dict_p_t,
)

# Full cascade H and L translation
# Full cascade H and J translation
γ_ = 1.0
σ_pulse = 1 / γ_
T = [0:0.002:1;] * 12σ_pulse
Expand All @@ -79,25 +79,25 @@ function benchmark_translation!(SUITE)
dict_p_cav = Dict([γ_sym, Δ_sym, gv_sym] .=> [γ_, 0.0, 0])
dict_p_t_cav = Dict([gu_sym, gv_sym] .=> [gu_t, gv_t])

SUITE["Translation"]["time-dependent"]["3-cavity H+L"] = @benchmarkable begin
SUITE["Translation"]["time-dependent"]["3-cavity H+J"] = @benchmarkable begin
to_numeric($H_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav)
to_numeric($L_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav)
to_numeric($J_sym, $b_cav; parameter = $dict_p_cav, time_parameter = $dict_p_t_cav)
end

## --- Closure evaluation (the ODE hot loop) ---

SUITE["Translation"]["closure evaluation"] = BenchmarkGroup()

H_QO = to_numeric(H_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav)
L_QO = to_numeric(L_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav)
J_QO = to_numeric(J_sym, b_cav; parameter = dict_p_cav, time_parameter = dict_p_t_cav)

t_mid = T[length(T)÷2]

SUITE["Translation"]["closure evaluation"]["3-cavity H(t)"] =
@benchmarkable $H_QO($t_mid)

SUITE["Translation"]["closure evaluation"]["3-cavity L(t)"] =
@benchmarkable $L_QO($t_mid)
SUITE["Translation"]["closure evaluation"]["3-cavity J(t)"] =
@benchmarkable $J_QO($t_mid)

return nothing
end
16 changes: 8 additions & 8 deletions docs/src/implementation.md
Original file line number Diff line number Diff line change
Expand Up @@ -19,20 +19,20 @@ av = Destroy(h, :a_v, 3)
@variables γ::Real g_u::Complex g_v::Complex
```

An SLH component is represented as `(S, L, H)` by the [`SLH`](@ref) type. The cascade [`▷`](@ref), concatenation [`⊞`](@ref), and feedback reduction [`feedback`](@ref) rules implement the standard network composition from the SLH framework.
An SLH component is represented as `(S, J, H)` by the [`SLH`](@ref) type. The cascade [`▷`](@ref), concatenation [`⊞`](@ref), and feedback reduction [`feedback`](@ref) rules implement the standard network composition from the SLH framework.

The resulting effective operators are accessed by [`hamiltonian`](@ref) and [`jump_operator`](@ref) and remain symbolic until translation. This is especially useful when you want to further manipulate the expressions, e.g. to transform into the interaction picture.

```julia
G_cas = ▷(G_u, G_s, G_v)

H = hamiltonian(G_cas)
L = jump_operator(G_cas)
J = jump_operator(G_cas)
```

For networks with internal loops, the symbolic model can be reduced directly with [`feedback`](@ref), which applies the SLH feedback reduction rule before translation. This keeps the symbolic workflow consistent: build a network from cascades and concatenations, eliminate internal connections symbolically, and only then translate the reduced Hamiltonian and Lindblad operators to numerics.

If you directly want to use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions, you can pass numeric operators and callables directly to [`SLH`](@ref), which supports both symbolic and numeric operator types as well as time-dependent `L` or `H` while still using the same cascade, concatenate, and feedback rules. This can be much faster.
If you directly want to use [QuantumOptics.jl](https://github.com/qojulia/QuantumOptics.jl) operators and functions, you can pass numeric operators and callables directly to [`SLH`](@ref), which supports both symbolic and numeric operator types as well as time-dependent `J` or `H` while still using the same cascade, concatenate, and feedback rules. This can be much faster.

## Translate to numerics

Expand All @@ -55,7 +55,7 @@ gu_t = coupling_input(t -> exp(-t^2), 0:0.01:5)
dict_p_t = Dict(g_u => gu_t)

H_QO = to_numeric(H, b; parameter=dict_p, time_parameter=dict_p_t)
L_QO = to_numeric(L, b; parameter=dict_p, time_parameter=dict_p_t)
J_QO = to_numeric(J, b; parameter=dict_p, time_parameter=dict_p_t)
```

In some cases it can be useful to define your own set of numeric operators which should replace the symbolic expressions, e.g. to reduce the Hilbert space if the output cavities are not analyzed but they are already included in the symbolic derivation. Such a list of operators can be provide with the dictionary `operators`.
Expand Down Expand Up @@ -105,16 +105,16 @@ gv2_eff = coupling_output(v2_eff, T)
The dominant output modes are extracted by computing the two-time correlation matrix

```math
g^{(1)}(t_1, t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle
g^{(1)}(t_1, t_2) = \langle J_s^\dagger(t_1) J_s(t_2) \rangle
```

and diagonalizing it. In this package, [`correlation_matrix`](@ref) builds that matrix from a previously computed trajectory $\rho(t)$ and a chosen output operator $L_s(t)$, using the quantum regression theorem. This means, for each time point $t_1$ we calculate $L_s(t_1) \rho(t_1)$ and use this as the initial "state" for the propagation of $t_2$, with the same Hamiltonian and Lindblad terms.
and diagonalizing it. In this package, [`correlation_matrix`](@ref) builds that matrix from a previously computed trajectory $\rho(t)$ and a chosen output operator $J_s(t)$, using the quantum regression theorem. This means, for each time point $t_1$ we calculate $J_s(t_1) \rho(t_1)$ and use this as the initial "state" for the propagation of $t_2$, with the same Hamiltonian and Lindblad terms.

The eigenvectors of the matrix $g^{(1)}(t_1, t_2)$ correspond to temporal modes and the eigenvalues to their mean photon-number weights. The full procedure is illustrated in the [Tutorial](@ref).

```julia
Ls(t) = gu_t(t) * au_qo + √(1.0) * c_qo
g1 = correlation_matrix(T, ρt, input_output_1, Ls)
Js(t) = gu_t(t) * au_qo + √(1.0) * c_qo
g1 = correlation_matrix(T, ρt, input_output_1, Js)
F = eigen(g1)
v_mode = F.vectors[:, end] / sqrt(T[2] - T[1])
```
Expand Down
24 changes: 12 additions & 12 deletions docs/src/theory.md
Original file line number Diff line number Diff line change
Expand Up @@ -86,10 +86,10 @@ This is the multi-mode extension used by [`effective_input_mode`](@ref) and

## SLH Networks for Cascaded Pulses

In the SLH formalism, each component is specified by a triple `(S, L, H)`:
In the SLH formalism, each component is specified by a triple `(S, J, H)`:

- `S`cattering matrix
- `L`indblad operators
- `J`ump (Lindblad) operators
- `H`amiltonian

Networks are built using composition rules:
Expand All @@ -99,8 +99,8 @@ Networks are built using composition rules:
```math
G_1 \triangleright G_2 = \left(
S_2 S_1,\;
L_2 + S_2 L_1,\;
H_1 + H_2 + \frac{1}{2i}\left(L_2^\dagger S_2 L_1 - L_1^\dagger S_2^\dagger L_2\right)
J_2 + S_2 J_1,\;
H_1 + H_2 + \frac{1}{2i}\left(J_2^\dagger S_2 J_1 - J_1^\dagger S_2^\dagger J_2\right)
\right).
```

Expand All @@ -113,27 +113,27 @@ S_1 & 0 \\
0 & S_2
\end{bmatrix},\;
\begin{bmatrix}
L_1 \\
L_2
J_1 \\
J_2
\end{bmatrix},\;
H_1 + H_2
\right).
```

3. **Feedback reduction**: connecting output port $x$ back into input port $y$ of an $n$-port component
$G = (S, L, H)$ yields a reduced $(n-1)$-port model $[(G)]_{x \to y} = (\tilde S, \tilde L, \tilde H)$ with
$G = (S, J, H)$ yields a reduced $(n-1)$-port model $[(G)]_{x \to y} = (\tilde S, \tilde J, \tilde H)$ with

```math
\tilde S = S_{\bar x,\bar y} + S_{\bar x,y}(1 - S_{x,y})^{-1}S_{x,\bar y},
```

```math
\tilde L = L_{\bar x} + S_{\bar x,y}(1 - S_{x,y})^{-1}L_x,
\tilde J = J_{\bar x} + S_{\bar x,y}(1 - S_{x,y})^{-1}J_x,
```

```math
\tilde H = H + \frac{1}{2i}\left[
\left(\sum_{j=1}^n L_j^\dagger S_{j,y}\right)(1 - S_{x,y})^{-1}L_x - \mathrm{h.c.}
\left(\sum_{j=1}^n J_j^\dagger S_{j,y}\right)(1 - S_{x,y})^{-1}J_x - \mathrm{h.c.}
\right].
```

Expand All @@ -144,7 +144,7 @@ obtain an effective SLH triple for the full problem. This describes a master equ

## Common SLH Elements

Appendix 1 of [Combes et al. (2017)](https://doi.org/10.1080/23746149.2017.1343097) lists SLH triples for commonly used components. The most frequently used elements in this package are summarized below using the standard ``(S, L, H)`` ordering.
Appendix 1 of [Combes et al. (2017)](https://doi.org/10.1080/23746149.2017.1343097) lists SLH triples for commonly used components. The most frequently used elements in this package are summarized below using the standard ``(S, J, H)`` ordering.

Phase shifter (single input/output):

Expand Down Expand Up @@ -210,10 +210,10 @@ The output field is generally multimode. To determine a proper **basis of tempor
occupations, we compute the first-order correlation function

```math
g^{(1)}(t_1, t_2) = \langle \hat L_s^\dagger(t_1)\, \hat L_s(t_2) \rangle,
g^{(1)}(t_1, t_2) = \langle \hat J_s^\dagger(t_1)\, \hat J_s(t_2) \rangle,
```

where $\hat L_s$ is the output operator (e.g., $\hat L_s = g_u(t) \hat{a}_u + \sqrt{\kappa} \hat c$ for a single sided cavity interacting with a pulse$u(t)$).
where $\hat J_s$ is the output operator (e.g., $\hat J_s = g_u(t) \hat{a}_u + \sqrt{\kappa} \hat c$ for a single sided cavity interacting with a pulse$u(t)$).
The eigen-decomposition

```math
Expand Down
16 changes: 8 additions & 8 deletions docs/src/tutorial.md
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ H = hamiltonian(G_cas)
```

```@example tutorial
L = jump_operator(G_cas)[1]
J = jump_operator(G_cas)[1]
```

## 2. Numerical parameters and input pulse
Expand Down Expand Up @@ -88,7 +88,7 @@ bv1 = FockBasis(2)
b = bu1 ⊗ bc1 ⊗ bv1

H_QO = to_numeric(H, b; parameter=dict_p, time_parameter=dict_p_t)
L_QO = to_numeric(L, b; parameter=dict_p, time_parameter=dict_p_t)
J_QO = to_numeric(J, b; parameter=dict_p, time_parameter=dict_p_t)
nothing # hide
```

Expand All @@ -99,7 +99,7 @@ We now solve the master equation. The required callback for `timeevolution.maste
```@example tutorial
function input_output_1(t, ρ)
Ht = H_QO(t)
J = [L_QO(t)]
J = [J_QO(t)]
return Ht, J, dagger.(J)
end

Expand All @@ -111,15 +111,15 @@ nothing # hide
## 4. Two-time correlation function

To extract the dominant output mode, we compute the two-time correlation matrix
``g^{(1)}(t_1,t_2) = \langle L_s^\dagger(t_1) L_s(t_2) \rangle`` and diagonalize it. To this end, we first define the desired numerical operators.
``g^{(1)}(t_1,t_2) = \langle J_s^\dagger(t_1) J_s(t_2) \rangle`` and diagonalize it. To this end, we first define the desired numerical operators.

```@example tutorial
au_qo = to_numeric(au, b)
c_qo = to_numeric(c, b)
av_qo = to_numeric(av, b)

Ls(t) = gu_t(t) * au_qo + √(γ_) * c_qo
g1_m = correlation_matrix(T, ρt, input_output_1, Ls)
Js(t) = gu_t(t) * au_qo + √(γ_) * c_qo
g1_m = correlation_matrix(T, ρt, input_output_1, Js)
nothing # hide
```

Expand Down Expand Up @@ -161,11 +161,11 @@ dict_p_t_2 = Dict([g_u, g_v] .=> [gu_t, gv_t])
dict_p_2 = Dict([γ, Δ] .=> [γ_, Δ_])

H_QO_2 = to_numeric(H, b; parameter=dict_p_2, time_parameter=dict_p_t_2)
L_QO_2 = to_numeric(L, b; parameter=dict_p_2, time_parameter=dict_p_t_2)
J_QO_2 = to_numeric(J, b; parameter=dict_p_2, time_parameter=dict_p_t_2)

function input_output_2(t, ρ)
Ht = H_QO_2(t)
J = [L_QO_2(t)]
J = [J_QO_2(t)]
return Ht, J, dagger.(J)
end
nothing # hide
Expand Down
Loading
Loading