diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0c38a4b1c..24aa649e5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Fix Dynamical Fock Dimension states saving due to wrong saving of dimensions. ([#375])
+- Support a list of observables for `expect`. ([#374], [#376])
 
 ## [v0.25.0]
 Release date: 2025-01-20
@@ -85,4 +86,6 @@ Release date: 2024-11-13
 [#360]: https://github.com/qutip/QuantumToolbox.jl/issues/360
 [#370]: https://github.com/qutip/QuantumToolbox.jl/issues/370
 [#371]: https://github.com/qutip/QuantumToolbox.jl/issues/371
+[#374]: https://github.com/qutip/QuantumToolbox.jl/issues/374
 [#375]: https://github.com/qutip/QuantumToolbox.jl/issues/375
+[#376]: https://github.com/qutip/QuantumToolbox.jl/issues/376
diff --git a/src/qobj/functions.jl b/src/qobj/functions.jl
index 00791f344..3ed38394e 100644
--- a/src/qobj/functions.jl
+++ b/src/qobj/functions.jl
@@ -17,7 +17,7 @@ ket2dm(ψ::QuantumObject{KetQuantumObject}) = ψ * ψ'
 ket2dm(ρ::QuantumObject{OperatorQuantumObject}) = ρ
 
 @doc raw"""
-    expect(O::AbstractQuantumObject, ψ::Union{QuantumObject,Vector{QuantumObject}})
+    expect(O::Union{AbstractQuantumObject,Vector{AbstractQuantumObject}}, ψ::Union{QuantumObject,Vector{QuantumObject}})
 
 Expectation value of the [`Operator`](@ref) `O` with the state `ψ`. The state can be a [`Ket`](@ref), [`Bra`](@ref) or [`Operator`](@ref).
 
@@ -27,48 +27,75 @@ If `ψ` is a density matrix ([`Operator`](@ref)), the function calculates ``\tex
 
 The function returns a real number if `O` is of `Hermitian` type or `Symmetric` type, and returns a complex number otherwise. You can make an operator `O` hermitian by using `Hermitian(O)`.
 
-Note that `ψ` can also be given as a list of [`QuantumObject`](@ref), it returns a list of expectation values.
+!!! note "List of observables and states"
+    The observable `O` and state `ψ` can be given as a list of [`QuantumObject`](@ref), it returns a list of expectation values. If both of them are given as a list, it returns a `Matrix` of expectation values.
 
 # Examples
 
 ```jldoctest
-julia> ψ = 1 / √2 * (fock(10,2) + fock(10,4));
+julia> ψ1 = 1 / √2 * (fock(10,2) + fock(10,4));
+
+julia> ψ2 = coherent(10, 0.6 + 0.8im);
 
 julia> a = destroy(10);
 
-julia> expect(a' * a, ψ) |> round
+julia> expect(a' * a, ψ1) |> round
 3.0 + 0.0im
 
-julia> expect(Hermitian(a' * a), ψ) |> round
+julia> expect(Hermitian(a' * a), ψ1) |> round
 3.0
+
+julia> round.(expect([a' * a, a' + a, a], [ψ1, ψ2]), digits = 1)
+3×2 Matrix{ComplexF64}:
+ 3.0+0.0im  1.0+0.0im
+ 0.0+0.0im  1.2-0.0im
+ 0.0+0.0im  0.6+0.8im
 ```
 """
-function expect(O::AbstractQuantumObject{OperatorQuantumObject}, ψ::QuantumObject{KetQuantumObject})
-    return dot(ψ.data, O.data, ψ.data)
-end
+expect(O::AbstractQuantumObject{OperatorQuantumObject}, ψ::QuantumObject{KetQuantumObject}) =
+    dot(ψ.data, O.data, ψ.data)
 expect(O::AbstractQuantumObject{OperatorQuantumObject}, ψ::QuantumObject{BraQuantumObject}) = expect(O, ψ')
 expect(O::QuantumObject{OperatorQuantumObject}, ρ::QuantumObject{OperatorQuantumObject}) = tr(O * ρ)
-function expect(
+expect(
     O::QuantumObject{OperatorQuantumObject,DimsType,<:Union{<:Hermitian{TF},<:Symmetric{TR}}},
     ψ::QuantumObject{KetQuantumObject},
-) where {DimsType<:AbstractDimensions,TF<:Number,TR<:Real}
-    return real(dot(ψ.data, O.data, ψ.data))
-end
-function expect(
+) where {DimsType<:AbstractDimensions,TF<:Number,TR<:Real} = real(dot(ψ.data, O.data, ψ.data))
+expect(
     O::QuantumObject{OperatorQuantumObject,DimsType,<:Union{<:Hermitian{TF},<:Symmetric{TR}}},
     ψ::QuantumObject{BraQuantumObject},
-) where {DimsType<:AbstractDimensions,TF<:Number,TR<:Real}
-    return real(expect(O, ψ'))
-end
-function expect(
+) where {DimsType<:AbstractDimensions,TF<:Number,TR<:Real} = real(expect(O, ψ'))
+expect(
     O::QuantumObject{OperatorQuantumObject,DimsType,<:Union{<:Hermitian{TF},<:Symmetric{TR}}},
     ρ::QuantumObject{OperatorQuantumObject},
+) where {DimsType<:AbstractDimensions,TF<:Number,TR<:Real} = real(tr(O * ρ))
+expect(
+    O::AbstractVector{<:AbstractQuantumObject{OperatorQuantumObject,DimsType,<:Union{<:Hermitian{TF},<:Symmetric{TR}}}},
+    ρ::QuantumObject,
+) where {DimsType<:AbstractDimensions,TF<:Number,TR<:Real} = expect.(O, Ref(ρ))
+function expect(O::AbstractVector{<:AbstractQuantumObject{OperatorQuantumObject}}, ρ::QuantumObject)
+    result = Vector{ComplexF64}(undef, length(O))
+    result .= expect.(O, Ref(ρ))
+    return result
+end
+expect(O::AbstractQuantumObject{OperatorQuantumObject}, ρ::AbstractVector{<:QuantumObject}) = expect.(Ref(O), ρ)
+function expect(
+    O::AbstractVector{<:AbstractQuantumObject{OperatorQuantumObject,DimsType,<:Union{<:Hermitian{TF},<:Symmetric{TR}}}},
+    ρ::AbstractVector{<:QuantumObject},
 ) where {DimsType<:AbstractDimensions,TF<:Number,TR<:Real}
-    return real(tr(O * ρ))
+    N_ops = length(O)
+    result = Matrix{Float64}(undef, N_ops, length(ρ))
+    for i in 1:N_ops
+        result[i, :] .= expect.(Ref(O[i]), ρ)
+    end
+    return result
 end
-function expect(O::QuantumObject{OperatorQuantumObject}, ρ::Vector{<:QuantumObject})
-    _expect = _ρ -> expect(O, _ρ)
-    return _expect.(ρ)
+function expect(O::AbstractVector{<:AbstractQuantumObject{OperatorQuantumObject}}, ρ::AbstractVector{<:QuantumObject})
+    N_ops = length(O)
+    result = Matrix{ComplexF64}(undef, N_ops, length(ρ))
+    for i in 1:N_ops
+        result[i, :] .= expect.(Ref(O[i]), ρ)
+    end
+    return result
 end
 
 @doc raw"""
diff --git a/test/core-test/quantum_objects.jl b/test/core-test/quantum_objects.jl
index 8761a57af..ad09f848e 100644
--- a/test/core-test/quantum_objects.jl
+++ b/test/core-test/quantum_objects.jl
@@ -479,6 +479,18 @@
         ψlist = [normalize!(basis(N, 4) + x * basis(N, 3)) for x in xlist]
         @test all(expect(a', ψlist) .≈ xlist)
 
+        # when input is a vector of observables
+        ρlist = Hermitian.(ket2dm.(ψlist)) # an alternative way to calculate expectation values for a list of density matrices
+        Olist1 = [a' * a, a' + a, a]
+        Olist2 = [Hermitian(a' * a), Hermitian(a' + a)]
+        exp_val_1 = expect(Olist1, ψlist)
+        exp_val_2 = expect(Olist2, ψlist)
+        @test size(exp_val_1) == (3, 4)
+        @test size(exp_val_2) == (2, 4)
+        @test all(exp_val_1[1, :] .≈ exp_val_2[1, :] .≈ expect(ρlist, a' * a))
+        @test all(exp_val_1[2, :] .≈ exp_val_2[2, :] .≈ expect(ρlist, a' + a))
+        @test all(exp_val_1[3, :] .≈ expect(a, ρlist))
+
         @testset "Type Inference (expect)" begin
             @inferred expect(a, ψ)
             @inferred expect(a, ψ')
@@ -488,6 +500,12 @@
             @inferred variance(a, ρ)
             @inferred expect(a, ψlist)
             @inferred variance(a, ψlist)
+            @inferred expect(ρlist, a)
+            @inferred expect(Olist1, ψ)
+            @inferred expect(Olist1, ψ')
+            @inferred expect(Olist1, ρ)
+            @inferred expect(Olist1, ψlist)
+            @inferred expect(Olist2, ψlist)
         end
     end