diff --git a/CHANGELOG.md b/CHANGELOG.md
index 64145e2db..f7dee851a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Check for orthogonality breakdown in `Lanczos` solver for `spectrum`. ([#501])
+
 ## [v0.32.1]
 Release date: 2025-06-24
 
@@ -253,3 +255,4 @@ Release date: 2024-11-13
 [#487]: https://github.com/qutip/QuantumToolbox.jl/issues/487
 [#489]: https://github.com/qutip/QuantumToolbox.jl/issues/489
 [#494]: https://github.com/qutip/QuantumToolbox.jl/issues/494
+[#501]: https://github.com/qutip/QuantumToolbox.jl/issues/501
diff --git a/src/spectrum.jl b/src/spectrum.jl
index 367470605..b97488575 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -186,10 +186,14 @@ function _spectrum(
     Ivec = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(cT, D)) # same as vec(system_identity_matrix)
     wₖ = spre(A).data' * vT(Ivec)
 
-    # Store the norm of the Green's function before renormalizing |v₁> and <w₁|
-    gfNorm = abs(dot(wₖ, vₖ))
-    vₖ ./= sqrt(gfNorm)
-    wₖ ./= sqrt(gfNorm)
+    # Store the normalization factor for the Green's function before renormalizing |v₁> and <w₁|
+    gfNorm = dot(wₖ, vₖ)
+    if gfNorm ≈ 0.0
+        throw(AssertionError("⟨w₀|v₀⟩ = 0, please check your A and B operators."))
+    end
+    scalingF = sqrt(abs(gfNorm))
+    vₖ ./= scalingF
+    wₖ ./= conj(gfNorm/scalingF)
 
     # Handle input frequency range
     ωList = vT(convert(Vector{fT}, ωlist))  # Make sure they're real frequencies and potentially on GPU
@@ -203,6 +207,7 @@ function _spectrum(
     αₖ = cT(0)
     βₖ = cT(-1)
     δₖ = fT(+1)
+    βₖδₖ = βₖ * δₖ
 
     # Current and up to second-to-last A and B Euler sequences
     A₋₂ = vT(ones(cT, Length))
@@ -232,8 +237,8 @@ function _spectrum(
         αₖ = dot(w₊₁, vₖ)
 
         # Update A(k), B(k) and continuous fraction; normalization avoids overflow
-        Aₖ .= (-1im .* ωList .+ αₖ) .* A₋₁ .- (βₖ * δₖ) .* A₋₂
-        Bₖ .= (-1im .* ωList .+ αₖ) .* B₋₁ .- (βₖ * δₖ) .* B₋₂
+        Aₖ .= (-1im .* ωList .+ αₖ) .* A₋₁ .- βₖδₖ .* A₋₂
+        Bₖ .= (-1im .* ωList .+ αₖ) .* B₋₁ .- βₖδₖ .* B₋₂
         lanczosFactor₋₁ .= lanczosFactor
         lanczosFactor .= Aₖ ./ Bₖ
 
@@ -268,9 +273,16 @@ function _spectrum(
         wₖ, w₊₁ = w₊₁, wₖ
 
         # k-th off-diagonal elements
-        buf = dot(wₖ, vₖ)
-        δₖ = sqrt(abs(buf))
-        βₖ = buf / δₖ
+        βₖδₖ = dot(wₖ, vₖ)
+        if βₖδₖ ≈ 0.0
+            if solver.verbose > 0
+                @warn "spectrum(): solver::Lanczos experienced orthogonality breakdown after $(k) iterations"
+                @warn "spectrum(): βₖδₖ = $(βₖδₖ)"
+            end
+            break
+        end
+        δₖ = sqrt(abs(βₖδₖ))
+        βₖ = βₖδₖ / δₖ
 
         # Normalize (k+1)-th left/right vectors
         vₖ ./= δₖ
diff --git a/test/core-test/correlations_and_spectrum.jl b/test/core-test/correlations_and_spectrum.jl
index bd9a2a918..71c875034 100644
--- a/test/core-test/correlations_and_spectrum.jl
+++ b/test/core-test/correlations_and_spectrum.jl
@@ -50,6 +50,10 @@
         @test last(outlines) == "spectrum(): Consider increasing maxiter and/or tol"
     end
 
+    @testset "Orthogonal input vectors Lanczos" begin
+        @test_throws AssertionError spectrum(H, ω_l2, [c_ops[1]], a', a; solver = Lanczos())
+    end
+
     # tlist and τlist checks
     t_fft_wrong = [0, 1, 10]
     t_wrong1 = [1, 2, 3]