diff --git a/CHANGELOG.md b/CHANGELOG.md
index edd8b42f8..5cb34b902 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,8 +6,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
-- Add error message for bad input in state/operator generating functions ([#603])
+
+- Add error message for bad input in state/operator generating functions. ([#603])
 - Support time-dependent operators on `dsf_mesolve` ([#610])
+- Fix `thermal_dm` and `enr_thermal_dm` for extreme cases. ([#614], [#613])
 
 ## [v0.39.1]
 Release date: 2025-11-19
@@ -385,3 +387,5 @@ Release date: 2024-11-13
 [#596]: https://github.com/qutip/QuantumToolbox.jl/issues/596
 [#603]: https://github.com/qutip/QuantumToolbox.jl/issues/603
 [#610]: https://github.com/qutip/QuantumToolbox.jl/issues/610
+[#613]: https://github.com/qutip/QuantumToolbox.jl/issues/613
+[#614]: https://github.com/qutip/QuantumToolbox.jl/issues/614
diff --git a/src/qobj/energy_restricted.jl b/src/qobj/energy_restricted.jl
index f65ea84ca..31f7c19e2 100644
--- a/src/qobj/energy_restricted.jl
+++ b/src/qobj/energy_restricted.jl
@@ -1,5 +1,5 @@
 #=
-This file defines the energy restricted space structure.
+This file defines the excitation number restricted space structure.
 =#
 
 export EnrSpace, enr_state_dictionaries
@@ -194,13 +194,15 @@ function enr_thermal_dm(
     D = s_enr.size
     idx2state = s_enr.idx2state
 
-    diags = ComplexF64[prod((nvec ./ (1 .+ nvec)) .^ idx2state[idx]) for idx in 1:D]
-    diags /= sum(diags)
-
+    β = @. log(1 + 1 / nvec)
+    P = _complex_float_type(T)[
+        prod(_Boltzmann_weight(β[k], n_excite) for (k, n_excite) in pairs(idx2state[idx])) for idx in 1:D
+    ]
+    P /= sum(P)
     if getVal(sparse)
-        return QuantumObject(spdiagm(0 => diags), Operator(), s_enr)
+        return QuantumObject(spdiagm(0 => P), Operator(), s_enr)
     else
-        return QuantumObject(diagm(0 => diags), Operator(), s_enr)
+        return QuantumObject(diagm(0 => P), Operator(), s_enr)
     end
 end
 
diff --git a/src/qobj/states.jl b/src/qobj/states.jl
index 38cddd6c9..aad4ff7de 100644
--- a/src/qobj/states.jl
+++ b/src/qobj/states.jl
@@ -118,13 +118,14 @@ Density matrix for a thermal state (generating thermal state probabilities) with
 !!! warning "Beware of type-stability!"
     If you want to keep type stability, it is recommended to use `thermal_dm(N, n, sparse=Val(sparse))` instead of `thermal_dm(N, n, sparse=sparse)`. See [this link](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-value-type) and the [related Section](@ref doc:Type-Stability) about type stability for more details.
 """
-function thermal_dm(N::Int, n::Real; sparse::Union{Bool,Val} = Val(false))
-    β = log(1.0 / n + 1.0)
-    data = [exp(-β * ComplexF64(j)) for j in 0:(N-1)]
+function thermal_dm(N::Int, n::T; sparse::Union{Bool,Val} = Val(false)) where {T<:Real}
+    β = log(1 + 1 / n)
+    P = _complex_float_type(T)[_Boltzmann_weight(β, j) for j in 0:(N-1)]
+    P /= sum(P)
     if getVal(sparse)
-        return QuantumObject(spdiagm(0 => data ./ sum(data)), Operator(), N)
+        return QuantumObject(spdiagm(0 => P), Operator(), N)
     else
-        return QuantumObject(diagm(0 => data ./ sum(data)), Operator(), N)
+        return QuantumObject(diagm(0 => P), Operator(), N)
     end
 end
 
diff --git a/src/utilities.jl b/src/utilities.jl
index 6d647c08b..d50496b2b 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -137,6 +137,8 @@ _sparse_similar(A::AbstractArray, args...) = sparse(args...)
 
 _Ginibre_ensemble(n::Int, rank::Int = n) = randn(ComplexF64, n, rank) / sqrt(n)
 
+_Boltzmann_weight(β::T, E::Int) where {T<:Real} = (E != 0 || isfinite(β)) ? exp(-β * E) : one(T)
+
 makeVal(x::Val{T}) where {T} = x
 makeVal(x) = Val(x)
 
@@ -194,7 +196,8 @@ _complex_float_type(::Type{Complex{Int32}}) = ComplexF32
 _complex_float_type(::Type{Complex{Int64}}) = ComplexF64
 _complex_float_type(::Type{Complex{Float32}}) = ComplexF32
 _complex_float_type(::Type{Complex{Float64}}) = ComplexF64
-_complex_float_type(T::Type{<:Complex}) = T # Allow other untracked Complex types, like ForwardDiff.Dual
+_complex_float_type(T::Type{<:Real}) = Complex{T} # Allow other untracked Complex types, like ForwardDiff.Dual
+_complex_float_type(T::Type{<:Complex}) = T       # Allow other untracked Complex types, like ForwardDiff.Dual
 
 _convert_eltype_wordsize(::Type{T}, ::Val{64}) where {T<:Int} = Int64
 _convert_eltype_wordsize(::Type{T}, ::Val{32}) where {T<:Int} = Int32
diff --git a/test/core-test/enr_state_operator.jl b/test/core-test/enr_state_operator.jl
index 50bed88c8..bf63e362e 100644
--- a/test/core-test/enr_state_operator.jl
+++ b/test/core-test/enr_state_operator.jl
@@ -1,5 +1,6 @@
 @testitem "Excitation number restricted state space" begin
     using StaticArraysCore
+    using SparseArrays
 
     @testset "EnrSpace" begin
         s_enr = EnrSpace((2, 2, 3), 3)
@@ -32,14 +33,46 @@
         space_s = (Space(D1), Space(D2))
 
         # EnrSpace
-        dims_enr = (2, 2, 3)
-        excitations = 3
+        dims_enr = (2, 3, 2)
+        excitations = 4
         space_enr = EnrSpace(dims_enr, excitations)
-        ρ_enr = enr_thermal_dm(space_enr, rand(3))
         I_enr = enr_identity(space_enr)
         size_enr = space_enr.size
 
+        # enr_thermal_dm (extreme cases)
+        ρTd0 = enr_thermal_dm(space_enr, 0.0)
+        ρTs0 = enr_thermal_dm(space_enr, 0.0; sparse = Val(true))
+        ρTd∞ = enr_thermal_dm(space_enr, Inf)
+        ρTs∞ = enr_thermal_dm(space_enr, Inf; sparse = Val(true))
+        @test tr(ρTd0) ≈ tr(ρTs0) ≈ tr(ρTd∞) ≈ tr(ρTs∞) ≈ 1.0
+        @test ρTd0.data ≈ ρTs0.data ≈ fock_dm(size_enr, 0).data
+        @test ρTd∞.data ≈ ρTs∞.data ≈ maximally_mixed_dm(size_enr).data
+
+        # general case (also test BigFloat)
+        nvec = BigFloat[0.123, 0.456, 0.789]
+        ρTd = enr_thermal_dm(space_enr, nvec)
+        ρTs = enr_thermal_dm(space_enr, nvec; sparse = Val(true))
+        @test isoper(ρTd)
+        @test tr(ρTd) ≈ tr(ρTs) ≈ 1.0
+        @test diag(ρTd) ≈ Float64[
+            0.44317797863426783,
+            0.19545412249437527,
+            0.13879749880303996,
+            0.06121365374823841,
+            0.0434695463284246,
+            0.019171309140931812,
+            0.04854041974355738,
+            0.021407708875163092,
+            0.015202219370235007,
+            0.006704612120243388,
+            0.004761134637930744,
+            0.0020997961035927096,
+        ]
+        @test ρTs.data isa AbstractSparseMatrix
+        @test ρTd ≈ ρTs
+
         # tensor between normal and ENR space
+        ρ_enr = enr_thermal_dm(space_enr, rand(3))
         ρ_tot = tensor(ρ_s, ρ_enr)
         opstring = sprint((t, s) -> show(t, "text/plain", s), ρ_tot)
         datastring = sprint((t, s) -> show(t, "text/plain", s), ρ_tot.data)
diff --git a/test/core-test/states_and_operators.jl b/test/core-test/states_and_operators.jl
index 677d511b4..483d32d92 100644
--- a/test/core-test/states_and_operators.jl
+++ b/test/core-test/states_and_operators.jl
@@ -36,21 +36,31 @@
     end
 
     @testset "thermal state" begin
-        ρTd = thermal_dm(5, 0.123)
-        ρTs = thermal_dm(5, 0.123; sparse = Val(true))
+        N = 5
+
+        # extreme cases
+        ρTd0 = thermal_dm(N, 0.0)
+        ρTs0 = thermal_dm(N, 0.0; sparse = Val(true))
+        ρTd∞ = thermal_dm(N, Inf)
+        ρTs∞ = thermal_dm(N, Inf; sparse = Val(true))
+        @test tr(ρTd0) ≈ tr(ρTs0) ≈ tr(ρTd∞) ≈ tr(ρTs∞) ≈ 1.0
+        @test ρTd0 ≈ ρTs0 ≈ fock_dm(N, 0)
+        @test ρTd∞ ≈ ρTs∞ ≈ maximally_mixed_dm(N)
+
+        # general case (also test BigFloat)
+        ρTd = thermal_dm(N, big(0.123))
+        ρTs = thermal_dm(N, big(0.123); sparse = Val(true))
         @test isoper(ρTd)
-        @test ρTd.dims == [5]
-        @test tr(ρTd) ≈ 1.0
-        @test ρTd.data ≈ spdiagm(
-            0 => Float64[
-                0.8904859864731106,
-                0.09753319353178326,
-                0.010682620484781245,
-                0.0011700465891612583,
-                0.00012815292116369966,
-            ],
-        )
-        @test typeof(ρTs.data) <: AbstractSparseMatrix
+        @test ρTd.dims == [N]
+        @test tr(ρTd) ≈ tr(ρTs) ≈ 1.0
+        @test diag(ρTd) ≈ Float64[
+            0.8904859864731106,
+            0.09753319353178326,
+            0.010682620484781245,
+            0.0011700465891612583,
+            0.00012815292116369966,
+        ]
+        @test ρTs.data isa AbstractSparseMatrix
         @test ρTd ≈ ρTs
     end
 
diff --git a/test/ext-test/cpu/autodiff/autodiff.jl b/test/ext-test/cpu/autodiff/autodiff.jl
index 5141e8854..c95205e77 100644
--- a/test/ext-test/cpu/autodiff/autodiff.jl
+++ b/test/ext-test/cpu/autodiff/autodiff.jl
@@ -85,6 +85,26 @@ end
 n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2))
 
 @testset "Autodiff" verbose=true begin
+    @testset "ForwardDiff for thermal_dm" begin
+        N = 100  # dimension of the system
+        N_op = num(N) # number operator
+
+        function N_expect(p) # p = [n]
+            ρT = thermal_dm(N, p[1]; sparse = Val(true))
+            return real(expect(N_op, ρT))
+        end
+
+        # Average photon number for thermal state
+        n = rand(Float64)
+        p = [n]
+
+        # Use ForwardDiff.gradient to compute the gradient
+        grad = ForwardDiff.gradient(N_expect, p)[1]
+
+        # Compare the result
+        @test isapprox(grad, 1.0; atol = 1e-6)
+    end
+
     @testset "sesolve" verbose=true begin
         Ω = 1.0
         params = [Ω]