Replace n_th with n_thermal for QuTiP compatibility (#249)

* make `n_th` to `n_thermal` and align with qutip

* fix fontsize in docs

* improve type stability

benchmarks/eigenvalues.jl

@@ -9,10 +9,10 @@ function benchmark_eigenvalues!(SUITE)
     ωb = 1
     g = 0.2
     κ = 0.01
-    n_thermal = 0.1
+    n_th = 0.1
 
     H = ωc * a_d * a + ωb * b_d * b + g * (a + a_d) * (b + b_d)
-    c_ops = [√((1 + n_thermal) * κ) * a, √κ * b, √(n_thermal * κ) * a_d]
+    c_ops = [√((1 + n_th) * κ) * a, √κ * b, √(n_th * κ) * a_d]
     L = liouvillian(H, c_ops)
 
     SUITE["Eigenvalues"]["eigenstates"]["dense"] = @benchmarkable eigenstates($L)

docs/src/api.md

@@ -238,7 +238,7 @@ AbstractLinearMap
 QuantumToolbox.versioninfo
 QuantumToolbox.about
 gaussian
-n_th
+n_thermal
 row_major_reshape
 meshgrid
 _calculate_expectation!

docs/src/users_guide/steadystate.md

@@ -103,14 +103,11 @@ sol_me = mesolve(H, ψ0, tlist, c_op_list, e_ops=e_ops, progress_bar=false)
 exp_me = real(sol_me.expect[1, :])
 
 # plot the results
-fig = Figure(size = (500, 350), fontsize = 15)
+fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], 
     title = L"Decay of Fock state $|10\rangle$ in a thermal environment with $\langle n\rangle=2$",
     xlabel = "Time", 
-    ylabel = "Number of excitations", 
-    titlesize = 24,
-    xlabelsize = 20, 
-    ylabelsize = 20
+    ylabel = "Number of excitations",
 )
 lines!(ax, tlist, exp_mc, label = "Monte-Carlo", linewidth = 2, color = :blue)
 lines!(ax, tlist, exp_me, label = "Master Equation", linewidth = 2, color = :orange, linestyle = :dash)

src/time_evolution/time_evolution.jl

@@ -322,7 +322,7 @@ function liouvillian_generalized(
         end
 
         # Ohmic reservoir
-        N_th = n_th.(Ωp, T_list[i])
+        N_th = n_thermal.(Ωp, T_list[i])
         Sp₀ = QuantumObject(triu(X_op, 1), type = Operator, dims = dims)
         Sp₁ = QuantumObject(droptol!((@. Ωp * N_th * Sp₀.data), tol), type = Operator, dims = dims)
         Sp₂ = QuantumObject(droptol!((@. Ωp * (1 + N_th) * Sp₀.data), tol), type = Operator, dims = dims)

src/utilities.jl

@@ -3,7 +3,7 @@ Utilities:
     internal (or external) functions which will be used throughout the entire package
 =#
 
-export gaussian, n_th
+export gaussian, n_thermal
 export row_major_reshape, meshgrid
 
 @doc raw"""
@@ -34,15 +34,18 @@ where ``\mu`` and ``\sigma^2`` are the mean and the variance respectively.
 gaussian(x::Number, μ::Number, σ::Number) = exp(-(x - μ)^2 / (2 * σ^2))
 
 @doc raw"""
-    n_th(ω::Number, T::Real)
+    n_thermal(ω::Real, ω_th::Real)
 
-Gives the mean number of excitations in a mode with frequency ω at temperature T:
-``n_{\rm th} (\omega, T) = \frac{1}{e^{\omega/T} - 1}``
+Return the number of photons in thermal equilibrium for an harmonic oscillator mode with frequency ``\omega``, at the temperature described by ``\omega_{\textrm{th}} \equiv k_B T / \hbar``:
+```math
+n(\omega, \omega_{\textrm{th}}) = \frac{1}{e^{\omega/\omega_{\textrm{th}}} - 1},
+```
+where ``\hbar`` is the reduced Planck constant, and ``k_B`` is the Boltzmann constant.
 """
-function n_th(ω::Real, T::Real)::Float64
-    (T == 0 || ω == 0) && return 0.0
-    abs(ω / T) > 50 && return 0.0
-    return 1 / (exp(ω / T) - 1)
+function n_thermal(ω::Real, ω_th::Real)::Float64
+    x = exp(ω / ω_th)
+    n = ((x != 1) && (ω_th > 0)) ? (1.0 / (x - 1.0)) : 0.0
+    return n
 end
 
 _get_dense_similar(A::AbstractArray, args...) = similar(A, args...)

test/core-test/eigenvalues_and_operators.jl

@@ -50,10 +50,10 @@
     ωb = 1
     g = 0.01
     κ = 0.1
-    n_thermal = 0.01
+    n_th = 0.01
 
     H = ωc * a_d * a + ωb * b_d * b + g * (a + a_d) * (b + b_d)
-    c_ops = [√((1 + n_thermal) * κ) * a, √κ * b, √(n_thermal * κ) * a_d]
+    c_ops = [√((1 + n_th) * κ) * a, √κ * b, √(n_th * κ) * a_d]
     L = liouvillian(H, c_ops)
 
     # eigen solve for general matrices
@@ -103,10 +103,10 @@
         ωb = 1
         g = 0.01
         κ = 0.1
-        n_thermal = 0.01
+        n_th = 0.01
 
         H = ωc * a_d * a + ωb * b_d * b + g * (a + a_d) * (b + b_d)
-        c_ops = [√((1 + n_thermal) * κ) * a, √κ * b, √(n_thermal * κ) * a_d]
+        c_ops = [√((1 + n_th) * κ) * a, √κ * b, √(n_th * κ) * a_d]
         L = liouvillian(H, c_ops)
 
         @inferred eigenstates(H, sparse = false)

test/core-test/generalized_master_equation.jl

@@ -40,7 +40,7 @@
     a2 = Qobj(dense_to_sparse((U'*a*U).data[1:N_trunc, 1:N_trunc], tol))
     sm2 = Qobj(dense_to_sparse((U'*sm*U).data[1:N_trunc, 1:N_trunc], tol))
 
-    @test abs(expect(Xp' * Xp, steadystate(L1)) - n_th(1, Tlist[1])) / n_th(1, Tlist[1]) < 1e-4
+    @test abs(expect(Xp' * Xp, steadystate(L1)) - n_thermal(1, Tlist[1])) / n_thermal(1, Tlist[1]) < 1e-4
 
     @testset "Type Inference (liouvillian_generalized)" begin
         N_c = 30