diff --git a/examples/basics/getting_started.ipynb b/examples/basics/getting_started.ipynb
index 59290c2..6e3b31f 100644
--- a/examples/basics/getting_started.ipynb
+++ b/examples/basics/getting_started.ipynb
@@ -25,7 +25,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 1,
    "id": "gs_imports",
    "metadata": {},
    "outputs": [],
diff --git a/examples/misc/nonhermitian_loss_chain.ipynb b/examples/misc/nonhermitian_loss_chain.ipynb
index 512fae5..50fd23f 100644
--- a/examples/misc/nonhermitian_loss_chain.ipynb
+++ b/examples/misc/nonhermitian_loss_chain.ipynb
@@ -16,7 +16,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "id": "794ec812",
    "metadata": {},
    "outputs": [],
diff --git a/src/gpu/GPU_tk.jl b/src/gpu/GPU_tk.jl
index 561c6a8..83a25bd 100644
--- a/src/gpu/GPU_tk.jl
+++ b/src/gpu/GPU_tk.jl
@@ -2497,6 +2497,11 @@ Use `type=ComplexF32` (default, faster) or `type=ComplexF64` (safer at tight cut
 or on large systems where F32 eigendecomposition can produce NaN). `dtype` is accepted
 as an alias for `type` for consistency with other GPU entry points.
 
+`return_maxlinkdim=true` returns `(result, linkdims)` instead of just `result`, where
+`linkdims::Vector{Int}` is the reached MPS bond dimension per output column (the χ the
+Chebyshev recursion hit under the given `maxdim`/`cutoff`). Useful for cutoff/tolerance
+studies where χ is the observable.
+
 !!! note "Block averaging not supported"
     `reduce=:block` is **not available** for the exciton LDOS. In the MPO-based LDOS
     functions (`get_ldos_spatial_gpu`), block averaging is a cheap O(1) partial trace
@@ -2526,7 +2531,8 @@ function get_exciton_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals
                                        type::Type{<:Number}                  = ComplexF32,
                                        dtype::Union{Nothing,Type{<:Number}}  = nothing,
                                        verbose::Bool    = false,
-                                       printinfo::Bool  = false)
+                                       printinfo::Bool  = false,
+                                       return_maxlinkdim::Bool = false)
 
     _check_gpu("get_exciton_ldos_spatial_gpu")
     gpu_type = dtype === nothing ? type : dtype
@@ -2595,6 +2601,8 @@ function get_exciton_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals
 
     printinfo && println("  [gpu] exciton ldos dtype=$gpu_type")
 
+    linkdims = zeros(Int, nX)   # reached MPS bond dim per output column (see return_maxlinkdim)
+
     for (j, group) in enumerate(groups)
         last_linkdim = 0
 
@@ -2634,11 +2642,12 @@ function get_exciton_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals
         for iω in 1:Nω
             result[iω, j] /= length(group)
         end
+        linkdims[j] = last_linkdim
         (verbose || printinfo) && (j % 5 == 0 || j == nX) &&
             println("  [gpu] exciton ldos $j/$nX (X=$(Xs[j]), n_avg=$(length(group)))  maxlinkdim=$last_linkdim")
     end
 
-    return result
+    return return_maxlinkdim ? (result, linkdims) : result
 end
 
 
diff --git a/src/solvers/KPM_tk.jl b/src/solvers/KPM_tk.jl
index bb577ca..b91fd23 100644
--- a/src/solvers/KPM_tk.jl
+++ b/src/solvers/KPM_tk.jl
@@ -1633,6 +1633,11 @@ are generated over `x_start:x_end`, with `num_avg` subpositions per group.
 `kernel=:hodc` uses the HODC reconstruction (`eta`, `m_order`); otherwise the
 standard KPM kernels are available (`:jackson`, `:lorentz`, `:fejer`,
 `:dirichlet`).
+
+`return_maxlinkdim=true` returns `(result, linkdims)` instead of just `result`,
+where `linkdims::Vector{Int}` is the reached MPS bond dimension per output column
+(the χ the Chebyshev recursion hit under the given `maxdim`/`cutoff`). Mirrors the
+GPU entry point; useful for cutoff/tolerance studies where χ is the observable.
 """
 function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
                                   X_list           = nothing,
@@ -1649,7 +1654,8 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
                                   maxdim::Int      = 100,
                                   cutoff::Real     = 1e-8,
                                   verbose::Bool    = false,
-                                  printinfo::Bool  = false)
+                                  printinfo::Bool  = false,
+                                  return_maxlinkdim::Bool = false)
     _ensure_scale!(H)
     length(H.sites) == 2 * H.L ||
         error("get_exciton_ldos_spatial: H is not an exciton Hamiltonian (expected length(H.sites) == 2*H.L).")
@@ -1702,6 +1708,7 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
     nX     = length(groups)
     Xs     = first.(groups)
     result = zeros(Float64, Nomega, nX)
+    linkdims = zeros(Int, nX)   # reached MPS bond dim per output column (see return_maxlinkdim)
 
     for (j, group) in enumerate(groups)
         last_linkdim = 0
@@ -1718,12 +1725,13 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
             valid[iomega] || continue
             result[iomega, j] = accum_group[iomega] / denom[iomega]
         end
+        linkdims[j] = last_linkdim
 
         (verbose || printinfo) && (j % 5 == 0 || j == nX) &&
             println("  exciton ldos $j/$nX (X=$(Xs[j]), n_avg=$(length(group)))  maxlinkdim=$last_linkdim")
     end
 
-    return result
+    return return_maxlinkdim ? (result, linkdims) : result
 end
 
 function get_exciton_ldos(H::TBHamiltonian, X::Int, omega_phys::Real;