Skip to content

Commit fbe80da

Browse files
TendonFFFTendonFFF
authored
Bloch Redfield master equation implementation (#473)
1 parent 6b6f674 commit fbe80da

File tree

5 files changed

+317
-0
lines changed

5 files changed

+317
-0
lines changed

‎CHANGELOG.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
88
## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
99

1010
- Introduce `Lanczos` solver for `spectrum`. ([#476])
11+
- Add Bloch-Redfield master equation solver. ([#473])
1112

1213
## [v0.31.1]
1314
Release date: 2025-05-16
@@ -230,4 +231,5 @@ Release date: 2024-11-13
230231
[#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
231232
[#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
232233
[#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
234+
[#473]: https://github.com/qutip/QuantumToolbox.jl/issues/473
233235
[#476]: https://github.com/qutip/QuantumToolbox.jl/issues/476

‎docs/src/resources/api.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -208,6 +208,9 @@ smesolve
208208
dfd_mesolve
209209
liouvillian
210210
liouvillian_generalized
211+
bloch_redfield_tensor
212+
brterm
213+
brmesolve
211214
```
212215

213216
### [Steady State Solvers](@id doc-API:Steady-State-Solvers)

‎src/QuantumToolbox.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,7 @@ include("time_evolution/callback_helpers/mcsolve_callback_helpers.jl")
106106
include("time_evolution/callback_helpers/ssesolve_callback_helpers.jl")
107107
include("time_evolution/callback_helpers/smesolve_callback_helpers.jl")
108108
include("time_evolution/mesolve.jl")
109+
include("time_evolution/brmesolve.jl")
109110
include("time_evolution/lr_mesolve.jl")
110111
include("time_evolution/sesolve.jl")
111112
include("time_evolution/mcsolve.jl")

‎src/time_evolution/brmesolve.jl‎

Lines changed: 196 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,196 @@
1+
export bloch_redfield_tensor, brterm, brmesolve
2+
3+
@doc raw"""
4+
bloch_redfield_tensor(
5+
H::QuantumObject{Operator},
6+
a_ops::Union{AbstractVector, Tuple},
7+
c_ops::Union{AbstractVector, Tuple, Nothing}=nothing;
8+
sec_cutoff::Real=0.1,
9+
fock_basis::Union{Val,Bool}=Val(false)
10+
)
11+
12+
Calculates the Bloch-Redfield tensor ([`SuperOperator`](@ref)) for a system given a set of operators and corresponding spectral functions that describes the system's coupling to its environment.
13+
14+
## Arguments
15+
16+
- `H`: The system Hamiltonian. Must be an [`Operator`](@ref)
17+
- `a_ops`: Nested list with each element is a `Tuple` of operator-function pairs `(a_op, spectra)`, and the coupling [`Operator`](@ref) `a_op` must be hermitian with corresponding `spectra` being a `Function` of transition energy
18+
- `c_ops`: List of collapse operators corresponding to Lindblad dissipators
19+
- `sec_cutoff`: Cutoff for secular approximation. Use `-1` if secular approximation is not used when evaluating bath-coupling terms.
20+
- `fock_basis`: Whether to return the tensor in the input (fock) basis or the diagonalized (eigen) basis.
21+
22+
## Return
23+
24+
The return depends on `fock_basis`.
25+
26+
- `fock_basis=Val(true)`: return the Bloch-Redfield tensor (in the fock basis) only.
27+
- `fock_basis=Val(false)`: return the Bloch-Redfield tensor (in the eigen basis) along with the transformation matrix from eigen to fock basis.
28+
"""
29+
function bloch_redfield_tensor(
30+
H::QuantumObject{Operator},
31+
a_ops::Union{AbstractVector,Tuple},
32+
c_ops::Union{AbstractVector,Tuple,Nothing} = nothing;
33+
sec_cutoff::Real = 0.1,
34+
fock_basis::Union{Val,Bool} = Val(false),
35+
)
36+
rst = eigenstates(H)
37+
U = QuantumObject(rst.vectors, Operator(), H.dimensions)
38+
sec_cutoff = float(sec_cutoff)
39+
40+
# in fock basis
41+
R0 = liouvillian(H, c_ops)
42+
43+
# set fock_basis=Val(false) and change basis together at the end
44+
R1 = 0
45+
isempty(a_ops) || (R1 += mapreduce(x -> _brterm(rst, x[1], x[2], sec_cutoff, Val(false)), +, a_ops))
46+
47+
SU = sprepost(U, U') # transformation matrix from eigen basis back to fock basis
48+
if getVal(fock_basis)
49+
return R0 + SU * R1 * SU'
50+
else
51+
return SU' * R0 * SU + R1, U
52+
end
53+
end
54+
55+
@doc raw"""
56+
brterm(
57+
H::QuantumObject{Operator},
58+
a_op::QuantumObject{Operator},
59+
spectra::Function;
60+
sec_cutoff::Real=0.1,
61+
fock_basis::Union{Bool, Val}=Val(false)
62+
)
63+
64+
Calculates the contribution of one coupling operator to the Bloch-Redfield tensor.
65+
66+
## Argument
67+
68+
- `H`: The system Hamiltonian. Must be an [`Operator`](@ref)
69+
- `a_op`: The operator coupling to the environment. Must be hermitian.
70+
- `spectra`: The corresponding environment spectra as a `Function` of transition energy.
71+
- `sec_cutoff`: Cutoff for secular approximation. Use `-1` if secular approximation is not used when evaluating bath-coupling terms.
72+
- `fock_basis`: Whether to return the tensor in the input (fock) basis or the diagonalized (eigen) basis.
73+
74+
## Return
75+
76+
The return depends on `fock_basis`.
77+
78+
- `fock_basis=Val(true)`: return the Bloch-Redfield term (in the fock basis) only.
79+
- `fock_basis=Val(false)`: return the Bloch-Redfield term (in the eigen basis) along with the transformation matrix from eigen to fock basis.
80+
"""
81+
function brterm(
82+
H::QuantumObject{Operator},
83+
a_op::QuantumObject{Operator},
84+
spectra::Function;
85+
sec_cutoff::Real = 0.1,
86+
fock_basis::Union{Bool,Val} = Val(false),
87+
)
88+
rst = eigenstates(H)
89+
term = _brterm(rst, a_op, spectra, sec_cutoff, makeVal(fock_basis))
90+
if getVal(fock_basis)
91+
return term
92+
else
93+
return term, Qobj(rst.vectors, Operator(), rst.dimensions)
94+
end
95+
end
96+
97+
function _brterm(
98+
rst::EigsolveResult,
99+
a_op::T,
100+
spectra::F,
101+
sec_cutoff::Real,
102+
fock_basis::Union{Val{true},Val{false}},
103+
) where {T<:QuantumObject{Operator},F<:Function}
104+
_check_br_spectra(spectra)
105+
106+
U = rst.vectors
107+
Id = I(prod(rst.dimensions))
108+
109+
skew = @. rst.values - rst.values' |> real
110+
spectrum = spectra.(skew)
111+
112+
A_mat = U' * a_op.data * U
113+
114+
ac_term = (A_mat .* spectrum) * A_mat
115+
bd_term = A_mat * (A_mat .* trans(spectrum))
116+
117+
if sec_cutoff != -1
118+
m_cut = similar(skew)
119+
map!(x -> abs(x) < sec_cutoff, m_cut, skew)
120+
ac_term .*= m_cut
121+
bd_term .*= m_cut
122+
123+
vec_skew = vec(skew)
124+
M_cut = @. abs(vec_skew - vec_skew') < sec_cutoff
125+
end
126+
127+
out =
128+
0.5 * (
129+
+ _sprepost(A_mat .* trans(spectrum), A_mat) + _sprepost(A_mat, A_mat .* spectrum) - _spost(ac_term, Id) -
130+
_spre(bd_term, Id)
131+
)
132+
133+
(sec_cutoff != -1) && (out .*= M_cut)
134+
135+
if getVal(fock_basis)
136+
SU = _sprepost(U, U')
137+
return QuantumObject(SU * out * SU', SuperOperator(), rst.dimensions)
138+
else
139+
return QuantumObject(out, SuperOperator(), rst.dimensions)
140+
end
141+
end
142+
143+
@doc raw"""
144+
brmesolve(
145+
H::QuantumObject{Operator},
146+
ψ0::QuantumObject,
147+
tlist::AbstractVector,
148+
a_ops::Union{Nothing, AbstractVector, Tuple},
149+
c_ops::Union{Nothing, AbstractVector, Tuple}=nothing;
150+
sec_cutoff::Real=0.1,
151+
e_ops::Union{Nothing, AbstractVector}=nothing,
152+
kwargs...,
153+
)
154+
155+
Solves for the dynamics of a system using the Bloch-Redfield master equation, given an input Hamiltonian, Hermitian bath-coupling terms and their associated spectral functions, as well as possible Lindblad collapse operators.
156+
157+
## Arguments
158+
159+
- `H`: The system Hamiltonian. Must be an [`Operator`](@ref)
160+
- `ψ0`: Initial state of the system $|\psi(0)\rangle$. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
161+
- `tlist`: List of times at which to save either the state or the expectation values of the system.
162+
- `a_ops`: Nested list with each element is a `Tuple` of operator-function pairs `(a_op, spectra)`, and the coupling [`Operator`](@ref) `a_op` must be hermitian with corresponding `spectra` being a `Function` of transition energy
163+
- `c_ops`: List of collapse operators corresponding to Lindblad dissipators
164+
- `sec_cutoff`: Cutoff for secular approximation. Use `-1` if secular approximation is not used when evaluating bath-coupling terms.
165+
- `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
166+
- `kwargs`: Keyword arguments for [`mesolve`](@ref).
167+
168+
## Notes
169+
170+
- This function will automatically generate the [`bloch_redfield_tensor`](@ref) and solve the time evolution with [`mesolve`](@ref).
171+
172+
# Returns
173+
174+
- `sol::TimeEvolutionSol`: The solution of the time evolution. See also [`TimeEvolutionSol`](@ref)
175+
"""
176+
function brmesolve(
177+
H::QuantumObject{Operator},
178+
ψ0::QuantumObject,
179+
tlist::AbstractVector,
180+
a_ops::Union{Nothing,AbstractVector,Tuple},
181+
c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
182+
sec_cutoff::Real = 0.1,
183+
e_ops::Union{Nothing,AbstractVector} = nothing,
184+
kwargs...,
185+
)
186+
R = bloch_redfield_tensor(H, a_ops, c_ops; sec_cutoff = sec_cutoff, fock_basis = Val(true))
187+
188+
return mesolve(R, ψ0, tlist, nothing; e_ops = e_ops, kwargs...)
189+
end
190+
191+
function _check_br_spectra(f::Function)
192+
meths = methods(f, [Real])
193+
length(meths.ms) == 0 &&
194+
throw(ArgumentError("The following function must accept one argument: `$(meths.mt.name)(ω)` with ω<:Real"))
195+
return nothing
196+
end

‎test/core-test/brmesolve.jl‎

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
@testitem "Bloch-Redfield tensor sec_cutoff" begin
2+
N = 5
3+
H = num(N)
4+
a = destroy(N)
5+
A_op = a+a'
6+
spectra(x) = (x>0) * 0.5
7+
for sec_cutoff in [0, 0.1, 1, 3, -1]
8+
R = bloch_redfield_tensor(H, [(A_op, spectra)], [a^2], sec_cutoff = sec_cutoff, fock_basis = true)
9+
R_eig, evecs = bloch_redfield_tensor(H, [(A_op, spectra)], [a^2], sec_cutoff = sec_cutoff, fock_basis = false)
10+
@test isa(R, QuantumObject)
11+
@test isa(R_eig, QuantumObject)
12+
@test isa(evecs, QuantumObject)
13+
14+
state = rand_dm(N) |> mat2vec
15+
fock_computed = R * state
16+
eig_computed = (sprepost(evecs, evecs') * R_eig * sprepost(evecs', evecs)) * state
17+
@test isapprox(fock_computed, eig_computed, atol = 1e-15)
18+
end
19+
end
20+
21+
@testitem "Compare brterm and Lindblad" begin
22+
N = 5
23+
H = num(N)
24+
a = destroy(N) + destroy(N)^2/2
25+
A_op = a+a'
26+
spectra(x) = x>0
27+
28+
# this test applies for limited cutoff
29+
lindblad = lindblad_dissipator(a)
30+
computation = brterm(H, A_op, spectra, sec_cutoff = 1.5, fock_basis = true)
31+
@test isapprox(lindblad, computation, atol = 1e-15)
32+
end
33+
34+
@testitem "brterm basis" begin
35+
N = 5
36+
H = num(N)
37+
a = destroy(N) + destroy(N)^2/2
38+
A_op = a+a'
39+
spectra(x) = x>0
40+
for sec_cutoff in [0, 0.1, 1, 3, -1]
41+
R = brterm(H, A_op, spectra, sec_cutoff = sec_cutoff, fock_basis = true)
42+
R_eig, evecs = brterm(H, A_op, spectra, sec_cutoff = sec_cutoff, fock_basis = false)
43+
@test isa(R, QuantumObject)
44+
@test isa(R_eig, QuantumObject)
45+
@test isa(evecs, QuantumObject)
46+
47+
state = rand_dm(N) |> mat2vec
48+
fock_computed = R * state
49+
eig_computed = (sprepost(evecs, evecs') * R_eig * sprepost(evecs', evecs)) * state
50+
@test isapprox(fock_computed, eig_computed, atol = 1e-15)
51+
end;
52+
end
53+
54+
@testitem "brterm sprectra function" begin
55+
f(x) = exp(x)/10
56+
function g(x)
57+
nbar = n_thermal(abs(x), 1)
58+
if x > 0
59+
return nbar
60+
elseif x < 0
61+
return 1 + nbar
62+
else
63+
return 0.0
64+
end
65+
end
66+
67+
spectra_list = [
68+
x -> (x>0), # positive frequency filter
69+
x -> one(x), # no filter
70+
f, # smooth filter
71+
g, # thermal field
72+
]
73+
74+
N = 5
75+
H = num(N)
76+
a = destroy(N) + destroy(N)^2/2
77+
A_op = a+a'
78+
for spectra in spectra_list
79+
R = brterm(H, A_op, spectra, sec_cutoff = 0.1, fock_basis = true)
80+
R_eig, evecs = brterm(H, A_op, spectra, sec_cutoff = 0.1, fock_basis = false)
81+
@test isa(R, QuantumObject)
82+
@test isa(R_eig, QuantumObject)
83+
@test isa(evecs, QuantumObject)
84+
85+
state = rand_dm(N) |> mat2vec
86+
fock_computed = R * state
87+
eig_computed = (sprepost(evecs, evecs') * R_eig * sprepost(evecs', evecs)) * state
88+
@test isapprox(fock_computed, eig_computed, atol = 1e-15)
89+
end
90+
end
91+
92+
@testitem "simple qubit system" begin
93+
pauli_vectors = [sigmax(), sigmay(), sigmaz()]
94+
γ = 0.25
95+
spectra(x) = γ * (x>=0)
96+
_m_c_op = √γ * sigmam()
97+
_z_c_op = √γ * sigmaz()
98+
_x_a_op = (sigmax(), spectra)
99+
100+
arg_sets = [[[_m_c_op], [], [_x_a_op]], [[_m_c_op], [_m_c_op], []], [[_m_c_op, _z_c_op], [_z_c_op], [_x_a_op]]]
101+
102+
δ = 0
103+
ϵ = 0.5 * 2π
104+
e_ops = pauli_vectors
105+
H = δ * 0.5 * sigmax() + ϵ * 0.5 * sigmaz()
106+
ψ0 = unit(2basis(2, 0) + basis(2, 1))
107+
times = LinRange(0, 10, 100)
108+
109+
for (me_c_ops, brme_c_ops, brme_a_ops) in arg_sets
110+
me = mesolve(H, ψ0, times, me_c_ops, e_ops = e_ops, progress_bar = Val(false))
111+
brme = brmesolve(H, ψ0, times, brme_a_ops, brme_c_ops, e_ops = e_ops, progress_bar = Val(false))
112+
113+
@test all(me.expect .== brme.expect)
114+
end
115+
end;

0 commit comments

Comments
 (0)