Minor improvements of spectrum ES

src/correlations.jl

@@ -168,18 +168,24 @@ function _spectrum(H::QuantumObject{<:AbstractArray{T1},HOpType},
 
     ω_l = ω_list
 
-    vals, vecs = eigen(L)
-    ρss = steadystate(L)
-    ρ0 = B.data * ρss.data
+    rates, vecs = eigen(L)
+
+    # Get the steady state and update the corresponding vector
+    ss_idx = argmin(abs2.(rates))
+    ρss = vec2mat(@view(vecs[:,ss_idx]))
+    ρss2 = (ρss + ρss') / 2
+    ρss2 ./= tr(ρss2)
+    ρss .= ρss2
+
+    ρ0 = B.data * ρss
     v = vecs \ mat2vec(ρ0)
 
-    amps = map(i->tr(A.data * reshape(v[i] * vecs[:,i], prod(Hdims), prod(Hdims))), eachindex(vals))
-    rates = vals
-    idxs = sortperm(rates, by=abs)
-    # For stability reasons, we take the modulus of the amplitudes
-    amps, rates = abs.(amps[idxs]), rates[idxs]
-    idxs = @. abs(amps) > solver.tol
+    amps = map(i->v[i] * tr(A.data * vec2mat(@view(vecs[:,i]))), eachindex(rates))
+    idxs = findall(x -> abs(x) > solver.tol, amps)
     amps, rates = amps[idxs], rates[idxs]
+    @. amps = abs(amps)
+    # idxs = findall(x -> real(x) < 0, amps)
+    # @. amps[idxs] -= 2*real(amps[idxs])
 
     spec = map(ω->2*real(sum(@. amps * (1 / (1im * ω - rates)))), ω_l)
 

src/general_functions.jl

@@ -107,7 +107,7 @@ end
 Returns the gaussian function ``\exp \left[- \frac{(x - \mu)^2}{2 \sigma^2} \right]``,
 where ``\mu`` and ``\sigma^2`` are the mean and the variance respectively.
 """
-gaussian(x::Number, μ::Number, σ::Number) = exp(-0.5 * (x - μ)^2 / σ^2)
+gaussian(x::Number, μ::Number, σ::Number) = exp(-(x - μ)^2 / (2 * σ^2))
 
 @doc raw"""
     ptrace(QO::QuantumObject, sel::Vector{Int})

src/quantum_object.jl

@@ -157,6 +157,14 @@ for op in (:(+), :(-), :(*), :(/), :(^))
         function Base.Broadcast.broadcasted(::typeof($op), x::Number, y::QuantumObject)
             return QuantumObject(broadcast($op, x, y.data), y.type, y.dims)
         end
+
+        function Base.Broadcast.broadcasted(::typeof($op), x::QuantumObject, y::AbstractArray)
+            return QuantumObject(broadcast($op, x.data, y), x.type, x.dims)
+        end
+
+        function Base.Broadcast.broadcasted(::typeof($op), x::AbstractArray, y::QuantumObject)
+            return QuantumObject(broadcast($op, x, y.data), y.type, y.dims)
+        end
     end
 end