From 59f784055e1318102e10d2c84c28d8b4629b1083 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Tue, 11 Aug 2026 15:21:32 +0900
Subject: [PATCH 1/6] Add quantics tensor train solver

---
 Project.toml     |   6 +-
 README.md        |   7 +-
 docs/make.jl     |   1 +
 docs/src/QTT.md  | 145 +++++++++++++++++++
 src/QTT.jl       | 360 +++++++++++++++++++++++++++++++++++++++++++++++
 src/TwoBody.jl   |   1 +
 test/QTT.jl      |  36 +++++
 test/runtests.jl |   1 +
 8 files changed, 553 insertions(+), 4 deletions(-)
 create mode 100644 docs/src/QTT.md
 create mode 100644 src/QTT.jl
 create mode 100644 test/QTT.jl

diff --git a/Project.toml b/Project.toml
index 74c6852..29a79c4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,16 +10,20 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e"
+TensorCrossInterpolation = "b261b2ec-6378-4871-b32e-9173bb050604"
 
 [compat]
 ArnoldiMethod = "0.4.0"
 FiniteDifferenceMatrices = "0.1.0"
 ForwardDiff = "0.10, 1"
 Optim = "1.9.4"
+QuanticsTCI = "0.7"
 SpecialFunctions = "2.3.1"
 Subscripts = "0.1.3"
-julia = "1.7"
+TensorCrossInterpolation = "0.9.16"
+julia = "1.9"
diff --git a/README.md b/README.md
index a52e217..a182b2b 100644
--- a/README.md
+++ b/README.md
@@ -18,13 +18,14 @@ flowchart TD
   A["Hamiltonian.jl"]
   C["Rayleigh-Ritz.jl"]
   F["FDM.jl"]
+  Q["QTT.jl"]
   G["VMC.jl"]
   H["DB.jl"]
   Z["TwoBody.jl"]
   A --> H
-  A --> C & F & G
-  H --> C & F & G
-  C & F & G --> Z
+  A --> C & F & Q & G
+  H --> C & F & Q & G
+  C & F & Q & G --> Z
 ```
 
 ## Developer's Guide
diff --git a/docs/make.jl b/docs/make.jl
index 47bffd1..61543f1 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -21,6 +21,7 @@ makedocs(;
         "Database" => "DB.md",
         "Rayleigh-Ritz Method" => "Rayleigh-Ritz.md",
         "Finite Difference Method" => "FDM.md",
+        "Quantics Tensor Train" => "QTT.md",
         "Variational Monte Carlo" => "VMC.md",
         "API reference" => "API.md",
     ],
diff --git a/docs/src/QTT.md b/docs/src/QTT.md
new file mode 100644
index 0000000..0c5fa99
--- /dev/null
+++ b/docs/src/QTT.md
@@ -0,0 +1,145 @@
+```@meta
+CurrentModule = TwoBody
+```
+
+# Quantics Tensor Train
+
+The quantics tensor train (QTT) method represents a radial grid of ``2^q`` points as
+``q`` binary tensor sites. A sampled radial function is then approximated by a chain
+of small tensor cores rather than stored as a vector with ``2^q`` entries,
+
+```math
+f(r_i)=F_{b_1\ldots b_q}
+\simeq G^{(1)}(b_1)G^{(2)}(b_2)\cdots G^{(q)}(b_q).
+```
+
+The Hamiltonian is represented in the same form as a matrix product operator (MPO),
+so the energy can be evaluated by contracting tensor networks without expanding the
+wave function or Hamiltonian,
+
+```math
+E_{\mathrm{QTT}}
+=\frac{\langle u_{\mathrm{QTT}}|H_{\mathrm{QTT}}|u_{\mathrm{QTT}}\rangle}
+{\langle u_{\mathrm{QTT}}|u_{\mathrm{QTT}}\rangle}.
+```
+
+When the tensor-train ranks remain moderate, both storage and contraction costs grow
+with ``q=\log_2 N`` rather than with the full grid size ``N``.
+
+## Theory
+
+For ``N=2^q`` interior grid points, write the zero-based grid index in binary,
+
+```math
+i-1 = \sum_{k=1}^{q} b_k 2^{q-k},
+\qquad b_k \in \{0,1\}.
+```
+
+Each vector core has dimensions
+``G^{(k)}\in\mathbb{R}^{\chi_{k-1}\times2\times\chi_k}``, with
+``\chi_0=\chi_q=1``. A matrix is quantized in both its row and column indices and
+represented as
+
+```math
+A_{b_1\ldots b_q,\,c_1\ldots c_q}
+\simeq W^{(1)}(b_1,c_1)W^{(2)}(b_2,c_2)\cdots W^{(q)}(b_q,c_q).
+```
+
+The central second-difference operator can be written using the one-point shift
+matrices ``S_-`` and ``S_+`` as
+
+```math
+D^{(2)} = \frac{S_- - 2I + S_+}{\Delta r^2}.
+```
+
+This tridiagonal operator has an exact QTT/MPO representation with maximum bond
+dimension three. Potential functions and trial wave functions are compressed by
+tensor cross interpolation. If ``\chi`` and ``\rho`` bound the vector and MPO bond
+dimensions, their storage is bounded by
+
+```math
+\operatorname{storage}(F_{\mathrm{QTT}}) \leq 2q\chi^2,
+\qquad
+\operatorname{storage}(A_{\mathrm{QTT}}) \leq 4q\rho^2.
+```
+
+Thus, when the QTT ranks remain moderate, storage grows as ``O(\log N)`` instead of
+``O(N)`` for a dense vector or ``O(N^2)`` for a dense matrix.
+
+## Usage
+
+For the three-dimensional spherical oscillator in atomic units, the normalized
+ground-state radial dependence is proportional to
+``\exp(-r^2/2)``. Configure the QTT grid with `quantics`; the number of grid points is
+`2^quantics`.
+
+```@example qtt
+using TwoBody
+
+H = Hamiltonian(
+  Kinetic(hbar=1, m=1),
+  PowerLaw(coefficient=1 / 2, exponent=2),
+)
+
+method = QuanticsTensorTrainMethod(
+  quantics=8,
+  r₀=0.0,
+  rₘₐₓ=8.0,
+  l=0,
+  tolerance=1e-9,
+  maxbonddim=32,
+)
+
+result = solve(H, r -> exp(-r^2 / 2), method; info=0)
+result.E
+```
+
+The result approaches the exact ground-state energy ``E=3/2`` as the grid is
+refined. `result.ψ`, `result.u`, and `result.H` contain the compressed radial wave
+function, reduced radial wave function, and reduced Hamiltonian. Their tensor-train
+ranks can be inspected without expanding any object:
+
+```@example qtt
+ranks(result.ψ)
+```
+
+Individual entries can also be contracted directly:
+
+```@example qtt
+qttvalue(result.ψ, 1)
+```
+
+!!! note "Current scope"
+    `QuanticsTensorTrainMethod` currently supports `Kinetic`, `RestEnergy`,
+    `Constant`, `Linear`, `Coulomb`, `PowerLaw`, `Gaussian`, `Exponential`, `Yukawa`,
+    and `Custom` operators. It evaluates a user-supplied trial wave function; a
+    tensor-train eigensolver is not yet included. `r₀` and `rₘₐₓ` are zero-value
+    (Dirichlet) boundaries and the ``2^q`` grid points lie strictly between them, so
+    choose boundaries where the reduced radial wave function is zero or sufficiently
+    small and perform a grid-convergence study.
+
+## Acknowledgments
+
+This work was developed in collaboration with Lucas Arenstein at
+[CompPhysHack 2026](https://qc-hybrid.github.io/CompPhysHack2026/). The authors
+gratefully acknowledge all those involved in organizing the hackathon for their
+contributions and support.
+
+## Bibliography
+
+The implementation follows the finite-difference QTT construction explored in the
+[CompPhysHack 2026 proof of concept](https://github.com/ohno/CompPhysHack2026Ohno/blob/main/julia/qtt.jl)
+and the method described by Arenstein, Mikkelsen, and Kastoryano in
+[Fast and Flexible Quantum-Inspired Differential Equation Solvers with Data Integration](https://doi.org/10.48550/arXiv.2505.17046).
+
+## API reference
+
+```@docs; canonical=false
+QuanticsTensorTrainMethod
+solve(hamiltonian::Hamiltonian, wavefunction::Function, method::QuanticsTensorTrainMethod; info=1)
+QTTVector
+QTTMatrix
+order
+ranks
+qttvalue
+```
diff --git a/src/QTT.jl b/src/QTT.jl
new file mode 100644
index 0000000..2f5bd76
--- /dev/null
+++ b/src/QTT.jl
@@ -0,0 +1,360 @@
+export QuanticsTensorTrainMethod, QTTVector, QTTMatrix, order, ranks, qttvalue
+
+import LinearAlgebra
+import QuanticsTCI
+import TensorCrossInterpolation
+
+"A vector stored as a quantics tensor train with cores `(left rank, 2, right rank)`."
+struct QTTVector{T}
+  cores::Vector{Array{T,3}}
+end
+
+"A matrix stored as a quantics tensor train with cores `(left rank, 2, 2, right rank)`."
+struct QTTMatrix{T}
+  cores::Vector{Array{T,4}}
+end
+
+"Return the number of binary sites in a QTT."
+order(tt::Union{QTTVector,QTTMatrix}) = length(tt.cores)
+
+"Return the left boundary, internal, and right boundary ranks of a QTT."
+ranks(tt::Union{QTTVector,QTTMatrix}) = [size(tt.cores[1], 1); [size(core, ndims(core)) for core in tt.cores]]
+
+struct QuanticsTensorTrainMethod
+  quantics::Int
+  r₀::Float64
+  rₘₐₓ::Float64
+  R::LinRange{Float64,Int}
+  Δr::Float64
+  l::Int
+  tolerance::Float64
+  maxbonddim::Int
+  maxiter::Int
+
+  function QuanticsTensorTrainMethod(;
+    quantics=10,
+    r₀=0.0,
+    rₘₐₓ=50.0,
+    l=0,
+    tolerance=1e-10,
+    maxbonddim=64,
+    maxiter=100,
+  )
+    quantics >= 2 || throw(ArgumentError("quantics must be at least 2"))
+    0 <= r₀ < rₘₐₓ || throw(ArgumentError("the radial interval must satisfy 0 <= r₀ < rₘₐₓ"))
+    l >= 0 || throw(ArgumentError("l must be nonnegative"))
+    tolerance > 0 || throw(ArgumentError("tolerance must be positive"))
+    maxbonddim >= 1 || throw(ArgumentError("maxbonddim must be positive"))
+    maxiter >= 1 || throw(ArgumentError("maxiter must be positive"))
+
+    npoints = 1 << quantics
+    Δr = Float64(rₘₐₓ - r₀) / (npoints + 1)
+    R = range(Float64(r₀) + Δr, Float64(rₘₐₓ) - Δr; length=npoints)
+    new(quantics, Float64(r₀), Float64(rₘₐₓ), R, Δr, l,
+        Float64(tolerance), maxbonddim, maxiter)
+  end
+end
+
+Base.string(method::QuanticsTensorTrainMethod) =
+  "QuanticsTensorTrainMethod(" *
+  join(["$(symbol)=$(getproperty(method, symbol))" for symbol in
+        (:quantics, :r₀, :rₘₐₓ, :l, :tolerance, :maxbonddim, :maxiter)], ", ") * ")"
+Base.show(io::IO, method::QuanticsTensorTrainMethod) = print(io, string(method))
+
+function _check_cores(tt::QTTVector)
+  isempty(tt.cores) && throw(ArgumentError("a QTT must contain at least one core"))
+  size(first(tt.cores), 1) == 1 || throw(ArgumentError("the first left rank must be one"))
+  size(last(tt.cores), 3) == 1 || throw(ArgumentError("the last right rank must be one"))
+  all(size(core, 2) == 2 for core in tt.cores) || throw(ArgumentError("every physical dimension must be two"))
+  all(size(tt.cores[k], 3) == size(tt.cores[k + 1], 1) for k in 1:(order(tt) - 1)) ||
+    throw(ArgumentError("adjacent QTT ranks do not match"))
+  return tt
+end
+
+function _check_cores(tt::QTTMatrix)
+  isempty(tt.cores) && throw(ArgumentError("a QTT must contain at least one core"))
+  size(first(tt.cores), 1) == 1 || throw(ArgumentError("the first left rank must be one"))
+  size(last(tt.cores), 4) == 1 || throw(ArgumentError("the last right rank must be one"))
+  all(size(core, 2) == size(core, 3) == 2 for core in tt.cores) ||
+    throw(ArgumentError("every physical dimension must be two by two"))
+  all(size(tt.cores[k], 4) == size(tt.cores[k + 1], 1) for k in 1:(order(tt) - 1)) ||
+    throw(ArgumentError("adjacent QTT ranks do not match"))
+  return tt
+end
+
+function qttvalue(tt::QTTVector{T}, index::Integer) where {T}
+  n = 1 << order(tt)
+  1 <= index <= n || throw(BoundsError(1:n, index))
+  state = ones(T, 1)
+  index0 = index - 1
+  for site in 1:order(tt)
+    core = tt.cores[site]
+    physical = Int((index0 >> (order(tt) - site)) & 1) + 1
+    state = vec(reshape(state, 1, :) * @view(core[:, physical, :]))
+  end
+  return only(state)
+end
+
+function qttvalue(tt::QTTMatrix{T}, row::Integer, column::Integer) where {T}
+  n = 1 << order(tt)
+  1 <= row <= n || throw(BoundsError(1:n, row))
+  1 <= column <= n || throw(BoundsError(1:n, column))
+  state = ones(T, 1)
+  row0, column0 = row - 1, column - 1
+  for site in 1:order(tt)
+    rowbit = Int((row0 >> (order(tt) - site)) & 1) + 1
+    columnbit = Int((column0 >> (order(tt) - site)) & 1) + 1
+    state = vec(reshape(state, 1, :) * @view(tt.cores[site][:, rowbit, columnbit, :]))
+  end
+  return only(state)
+end
+
+function Base.:*(factor::Number, tt::QTTVector)
+  cores = copy(tt.cores)
+  cores[1] = factor .* cores[1]
+  return QTTVector(cores)
+end
+
+function Base.:*(factor::Number, tt::QTTMatrix)
+  cores = copy(tt.cores)
+  cores[1] = factor .* cores[1]
+  return QTTMatrix(cores)
+end
+
+Base.:*(tt::Union{QTTVector,QTTMatrix}, factor::Number) = factor * tt
+Base.:-(tt::QTTMatrix) = -1 * tt
+Base.:-(A::QTTMatrix, B::QTTMatrix) = A + (-B)
+
+function Base.:+(A::QTTMatrix{TA}, B::QTTMatrix{TB}) where {TA,TB}
+  order(A) == order(B) || throw(DimensionMismatch("QTT orders must match"))
+  nsites = order(A)
+  T = promote_type(TA, TB)
+  cores = Vector{Array{T,4}}(undef, nsites)
+
+  rankA, rankB = size(A.cores[1], 4), size(B.cores[1], 4)
+  cores[1] = zeros(T, 1, 2, 2, rankA + rankB)
+  cores[1][1, :, :, 1:rankA] = A.cores[1]
+  cores[1][1, :, :, (rankA + 1):end] = B.cores[1]
+
+  for site in 2:(nsites - 1)
+    a, b = A.cores[site], B.cores[site]
+    aleft, aright = size(a, 1), size(a, 4)
+    bleft, bright = size(b, 1), size(b, 4)
+    core = zeros(T, aleft + bleft, 2, 2, aright + bright)
+    core[1:aleft, :, :, 1:aright] = a
+    core[(aleft + 1):end, :, :, (aright + 1):end] = b
+    cores[site] = core
+  end
+
+  a, b = A.cores[end], B.cores[end]
+  aleft, bleft = size(a, 1), size(b, 1)
+  cores[end] = zeros(T, aleft + bleft, 2, 2, 1)
+  cores[end][1:aleft, :, :, 1] = a[:, :, :, 1]
+  cores[end][(aleft + 1):end, :, :, 1] = b[:, :, :, 1]
+  return QTTMatrix(cores)
+end
+
+function _tridiagonal_qtt(quantics::Int, lower::Real, diagonal::Real, upper::Real)
+  identity2 = Matrix{Float64}(LinearAlgebra.I, 2, 2)
+  shift = [0.0 1.0; 0.0 0.0]
+
+  firstcore = zeros(Float64, 1, 2, 2, 3)
+  firstcore[1, :, :, 1] = identity2
+  firstcore[1, :, :, 2] = transpose(shift)
+  firstcore[1, :, :, 3] = shift
+
+  middlecore = zeros(Float64, 3, 2, 2, 3)
+  middlecore[1, :, :, 1] = identity2
+  middlecore[1, :, :, 2] = transpose(shift)
+  middlecore[1, :, :, 3] = shift
+  middlecore[2, :, :, 2] = shift
+  middlecore[3, :, :, 3] = transpose(shift)
+
+  lastcore = zeros(Float64, 3, 2, 2, 1)
+  lastcore[1, :, :, 1] = diagonal * identity2 + upper * shift + lower * transpose(shift)
+  lastcore[2, :, :, 1] = lower * shift
+  lastcore[3, :, :, 1] = upper * transpose(shift)
+
+  cores = Vector{Array{Float64,4}}(undef, quantics)
+  cores[1] = firstcore
+  for site in 2:(quantics - 1)
+    cores[site] = middlecore
+  end
+  cores[end] = lastcore
+  return QTTMatrix(cores)
+end
+
+function _diagonal(tt::QTTVector{T}) where {T}
+  cores = Vector{Array{T,4}}(undef, order(tt))
+  for site in 1:order(tt)
+    vectorcore = tt.cores[site]
+    core = zeros(T, size(vectorcore, 1), 2, 2, size(vectorcore, 3))
+    core[:, 1, 1, :] = vectorcore[:, 1, :]
+    core[:, 2, 2, :] = vectorcore[:, 2, :]
+    cores[site] = core
+  end
+  return QTTMatrix(cores)
+end
+
+function _qttvector_from_tci(tci, method::QuanticsTensorTrainMethod, f::Function)
+  tensortrain = TensorCrossInterpolation.TensorTrain(tci.tci)
+  tensortrain = TensorCrossInterpolation.reverse(tensortrain)
+  sitetensors = collect(TensorCrossInterpolation.sitetensors(tensortrain))
+  length(sitetensors) == method.quantics || throw(ArgumentError("unexpected tensor-train order"))
+  qtt = _check_cores(QTTVector([Float64.(Array(core)) for core in sitetensors]))
+
+  reversed = QTTVector([permutedims(core, (3, 2, 1)) for core in reverse(qtt.cores)])
+  probes = unique(clamp.([1, 2, length(method.R) ÷ 2, length(method.R) - 1, length(method.R)], 1, length(method.R)))
+  normalerror = sum(abs(qttvalue(qtt, i) - f(method.R[i])) for i in probes)
+  reversederror = sum(abs(qttvalue(reversed, i) - f(method.R[i])) for i in probes)
+  return reversederror < normalerror ? reversed : qtt
+end
+
+function _qttvector(f::Function, method::QuanticsTensorTrainMethod)
+  tci, _, _ = QuanticsTCI.quanticscrossinterpolate(
+    Float64,
+    r -> Float64(f(r)),
+    method.R;
+    tolerance=method.tolerance,
+    maxbonddim=method.maxbonddim,
+    maxiter=method.maxiter,
+  )
+  return _qttvector_from_tci(tci, method, f)
+end
+
+function _inner(left::QTTVector{TL}, right::QTTVector{TR}) where {TL,TR}
+  order(left) == order(right) || throw(DimensionMismatch("QTT orders must match"))
+  T = promote_type(TL, TR)
+  environment = ones(T, 1, 1)
+  for site in 1:order(left)
+    lcore, rcore = left.cores[site], right.cores[site]
+    next = zeros(T, size(lcore, 3), size(rcore, 3))
+    for lleft in axes(lcore, 1), rleft in axes(rcore, 1), physical in 1:2,
+        lright in axes(lcore, 3), rright in axes(rcore, 3)
+      next[lright, rright] += environment[lleft, rleft] *
+        conj(lcore[lleft, physical, lright]) * rcore[rleft, physical, rright]
+    end
+    environment = next
+  end
+  return only(environment)
+end
+
+function _expectation(left::QTTVector{TL}, operator::QTTMatrix{TO}, right::QTTVector{TR}) where {TL,TO,TR}
+  order(left) == order(operator) == order(right) || throw(DimensionMismatch("QTT orders must match"))
+  T = promote_type(TL, TO, TR)
+  environment = ones(T, 1, 1, 1)
+  for site in 1:order(left)
+    lcore, opcore, rcore = left.cores[site], operator.cores[site], right.cores[site]
+    next = zeros(T, size(lcore, 3), size(opcore, 4), size(rcore, 3))
+    for lleft in axes(lcore, 1), oleft in axes(opcore, 1), rleft in axes(rcore, 1),
+        row in 1:2, column in 1:2, lright in axes(lcore, 3),
+        oright in axes(opcore, 4), rright in axes(rcore, 3)
+      next[lright, oright, rright] += environment[lleft, oleft, rleft] *
+        conj(lcore[lleft, row, lright]) * opcore[oleft, row, column, oright] *
+        rcore[rleft, column, rright]
+    end
+    environment = next
+  end
+  return only(environment)
+end
+
+function matrix(term::Kinetic, method::QuanticsTensorTrainMethod)
+  h = method.Δr
+  secondderivative = _tridiagonal_qtt(method.quantics, 1 / h^2, -2 / h^2, 1 / h^2)
+  radialoperator = if method.l == 0
+    secondderivative
+  else
+    centrifugal = method.l * (method.l + 1) * _diagonal(_qttvector(r -> 1 / r^2, method))
+    secondderivative - centrifugal
+  end
+  return (-term.hbar^2 / (2 * term.m)) * radialoperator
+end
+
+matrix(term::RestEnergy, method::QuanticsTensorTrainMethod) =
+  _diagonal(_qttvector(_ -> term.m * term.c^2, method))
+
+const QTTPotentialTerm = Union{Constant,Linear,Coulomb,PowerLaw,Gaussian,Exponential,Yukawa,Custom}
+
+matrix(term::QTTPotentialTerm, method::QuanticsTensorTrainMethod) =
+  _diagonal(_qttvector(r -> V(term, r), method))
+
+matrix(term::PotentialTerm, ::QuanticsTensorTrainMethod) =
+  throw(ArgumentError("$(typeof(term)) is not supported by QuanticsTensorTrainMethod"))
+
+function matrix(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod)
+  isempty(hamiltonian.terms) && throw(ArgumentError("the Hamiltonian must contain at least one term"))
+  matrices = map(hamiltonian.terms) do term
+    term isa Union{Kinetic,RestEnergy,QTTPotentialTerm} ||
+      throw(ArgumentError("$(typeof(term)) is not supported by QuanticsTensorTrainMethod"))
+    matrix(term, method)
+  end
+  return reduce(+, matrices)
+end
+
+function solve(
+  hamiltonian::Hamiltonian,
+  wavefunction::Function,
+  method::QuanticsTensorTrainMethod;
+  info=1,
+)
+  H = matrix(hamiltonian, method)
+  unnormalizedψ = _qttvector(wavefunction, method)
+  unnormalizedu = _qttvector(r -> r * wavefunction(r), method)
+  denominator = _inner(unnormalizedu, unnormalizedu)
+  numerator = _expectation(unnormalizedu, H, unnormalizedu)
+  energy = real(numerator / denominator)
+  normalization = sqrt(real(4π * method.Δr * denominator))
+  ψ = (1 / normalization) * unnormalizedψ
+  u = (1 / normalization) * unnormalizedu
+
+  if info > 0
+    println("\n# method\n")
+    println(method)
+    println("\n# energy\n")
+    println("E = $energy\n")
+  end
+
+  return info >= 0 ? (
+    hamiltonian=hamiltonian,
+    method=method,
+    H=H,
+    E=energy,
+    ψ=ψ,
+    u=u,
+  ) : (E=energy,)
+end
+
+@doc raw"""
+    QuanticsTensorTrainMethod(; quantics=10, r₀=0.0, rₘₐₓ=50.0, l=0,
+                               tolerance=1e-10, maxbonddim=64, maxiter=100)
+
+Configure a radial quantics tensor train (QTT) calculation on ``2^q`` uniformly spaced
+interior points, where `q = quantics`. `r₀` and `rₘₐₓ` are zero-value (Dirichlet)
+boundaries and are excluded from the grid, so Coulomb and centrifugal potentials are
+never evaluated at the origin. `tolerance`, `maxbonddim`, and `maxiter` are forwarded
+to tensor cross interpolation.
+""" QuanticsTensorTrainMethod
+
+@doc raw"""
+    solve(hamiltonian, wavefunction, method::QuanticsTensorTrainMethod; info=1)
+
+Compress the Hamiltonian and trial radial wave function as QTTs and evaluate the
+Rayleigh quotient for the reduced radial wave function ``u(r)=r\psi(r)``,
+
+```math
+E = \frac{\langle u | H_u | u \rangle}{\langle u | u \rangle},
+\qquad
+H_u = -\frac{\hbar^2}{2m}
+\left[\frac{\mathrm{d}^2}{\mathrm{d}r^2}-\frac{l(l+1)}{r^2}\right]+V(r).
+```
+
+The finite-difference kinetic operator is represented directly as a rank-three QTT
+matrix; functions of ``r`` are compressed with tensor cross interpolation.
+""" solve(hamiltonian::Hamiltonian, wavefunction::Function, method::QuanticsTensorTrainMethod; info=1)
+
+@doc raw"""
+    qttvalue(tt::QTTVector, index)
+    qttvalue(tt::QTTMatrix, row, column)
+
+Contract a QTT at one vector or matrix index without materializing its ``2^q`` entries.
+""" qttvalue
diff --git a/src/TwoBody.jl b/src/TwoBody.jl
index 9a74c74..c08ae79 100644
--- a/src/TwoBody.jl
+++ b/src/TwoBody.jl
@@ -12,6 +12,7 @@ include("./Basis.jl")
 # Solvers
 include("./Rayleigh-Ritz.jl")
 include("./FDM.jl")
+include("./QTT.jl")
 include("./VMC.jl")
 
 end
diff --git a/test/QTT.jl b/test/QTT.jl
new file mode 100644
index 0000000..8336902
--- /dev/null
+++ b/test/QTT.jl
@@ -0,0 +1,36 @@
+@testset "QTT.jl" begin
+  @test_throws ArgumentError QuanticsTensorTrainMethod(quantics=1)
+  @test_throws ArgumentError QuanticsTensorTrainMethod(r₀=-0.1)
+
+  method = QuanticsTensorTrainMethod(
+    quantics=7,
+    r₀=0.0,
+    rₘₐₓ=8.0,
+    tolerance=1e-9,
+    maxbonddim=32,
+  )
+  @test length(method.R) == 2^method.quantics
+  @test method.R[1] ≈ method.r₀ + method.Δr
+  @test method.R[end] ≈ method.rₘₐₓ - method.Δr
+
+  tridiagonal = TwoBody._tridiagonal_qtt(3, -1, 2, -1)
+  for row in 1:8, column in 1:8
+    expected = row == column ? 2 : abs(row - column) == 1 ? -1 : 0
+    @test qttvalue(tridiagonal, row, column) ≈ expected atol=1e-12
+  end
+
+  H = Hamiltonian(
+    Kinetic(hbar=1, m=1),
+    PowerLaw(coefficient=1 / 2, exponent=2),
+  )
+  result = solve(H, r -> exp(-r^2 / 2), method; info=0)
+  @test result.E ≈ 1.5 atol=2e-3
+  @test order(result.ψ) == method.quantics
+  @test first(ranks(result.ψ)) == last(ranks(result.ψ)) == 1
+  @test qttvalue(result.ψ, 1) isa Real
+
+  potential = TwoBody.matrix(PowerLaw(coefficient=1 / 2, exponent=2), method)
+  @test qttvalue(potential, 5, 5) ≈ method.R[5]^2 / 2 atol=1e-7
+  @test qttvalue(potential, 5, 6) ≈ 0 atol=1e-12
+  @test_throws ArgumentError TwoBody.matrix(Delta(), method)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index cc66f2c..22a0fd3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,5 +12,6 @@ using ForwardDiff
   include("Basis.jl")
   include("Rayleigh-Ritz.jl")
   include("FDM.jl")
+  include("QTT.jl")
   include("VMC.jl")
 end

From ca0334c38275d6c39d33a066fe4e8346659a753a Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Tue, 11 Aug 2026 16:25:28 +0900
Subject: [PATCH 2/6] Use TensorTrainNumerics DMRG for QTT states

---
 Project.toml     |   9 +-
 docs/src/QTT.md  |  89 ++++++----
 src/QTT.jl       | 424 +++++++++++++++++++++++------------------------
 test/QTT.jl      |  43 +++--
 test/runtests.jl |   1 +
 5 files changed, 300 insertions(+), 266 deletions(-)

diff --git a/Project.toml b/Project.toml
index 29a79c4..c8f0af7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,22 +8,21 @@ ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
 FiniteDifferenceMatrices = "a7a66f33-e7b8-47af-b618-f9b5bea05f3d"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e"
-TensorCrossInterpolation = "b261b2ec-6378-4871-b32e-9173bb050604"
+TensorTrainNumerics = "304e70ee-b558-420e-aed5-de4d05b6ca2b"
 
 [compat]
 ArnoldiMethod = "0.4.0"
 FiniteDifferenceMatrices = "0.1.0"
 ForwardDiff = "0.10, 1"
 Optim = "1.9.4"
-QuanticsTCI = "0.7"
 SpecialFunctions = "2.3.1"
 Subscripts = "0.1.3"
-TensorCrossInterpolation = "0.9.16"
-julia = "1.9"
+TensorTrainNumerics = "1.1.3"
+julia = "1.10"
diff --git a/docs/src/QTT.md b/docs/src/QTT.md
index 0c5fa99..4f2ba5a 100644
--- a/docs/src/QTT.md
+++ b/docs/src/QTT.md
@@ -13,14 +13,15 @@ f(r_i)=F_{b_1\ldots b_q}
 \simeq G^{(1)}(b_1)G^{(2)}(b_2)\cdots G^{(q)}(b_q).
 ```
 
-The Hamiltonian is represented in the same form as a matrix product operator (MPO),
-so the energy can be evaluated by contracting tensor networks without expanding the
-wave function or Hamiltonian,
+The Hamiltonian is represented in the same form as a matrix product operator (MPO).
+DMRG solves the eigenvalue problem directly in tensor-train form; previously found
+states are shifted out of the low-energy spectrum to expose successive excitations,
 
 ```math
-E_{\mathrm{QTT}}
-=\frac{\langle u_{\mathrm{QTT}}|H_{\mathrm{QTT}}|u_{\mathrm{QTT}}\rangle}
-{\langle u_{\mathrm{QTT}}|u_{\mathrm{QTT}}\rangle}.
+H^{(n)}_{\mathrm{QTT}}u_n=E_nu_n,
+\qquad
+H^{(n)}_{\mathrm{QTT}}
+=H_{\mathrm{QTT}}+\mu\sum_{k=0}^{n-1}|u_k\rangle\langle u_k|.
 ```
 
 When the tensor-train ranks remain moderate, both storage and contraction costs grow
@@ -53,9 +54,13 @@ D^{(2)} = \frac{S_- - 2I + S_+}{\Delta r^2}.
 ```
 
 This tridiagonal operator has an exact QTT/MPO representation with maximum bond
-dimension three. Potential functions and trial wave functions are compressed by
-tensor cross interpolation. If ``\chi`` and ``\rho`` bound the vector and MPO bond
-dimensions, their storage is bounded by
+dimension three. Potential functions and initial wave functions are compressed by
+tensor cross interpolation. Two-site DMRG sweeps optimize the tensor cores and adapt
+their bond dimensions. After each state is found, its projector is added to the MPO
+and the result is compressed before the next DMRG solve.
+
+If ``\chi`` and ``\rho`` bound the vector and MPO bond dimensions, their storage is
+bounded by
 
 ```math
 \operatorname{storage}(F_{\mathrm{QTT}}) \leq 2q\chi^2,
@@ -68,13 +73,15 @@ Thus, when the QTT ranks remain moderate, storage grows as ``O(\log N)`` instead
 
 ## Usage
 
-For the three-dimensional spherical oscillator in atomic units, the normalized
-ground-state radial dependence is proportional to
-``\exp(-r^2/2)``. Configure the QTT grid with `quantics`; the number of grid points is
-`2^quantics`.
+For the three-dimensional spherical oscillator in atomic units, configure the QTT
+grid with `quantics`; the number of interior grid points is `2^quantics`. Here a
+coarse grid is used so the example runs quickly.
 
 ```@example qtt
 using TwoBody
+using Random
+
+Random.seed!(1234)
 
 H = Hamiltonian(
   Kinetic(hbar=1, m=1),
@@ -82,41 +89,61 @@ H = Hamiltonian(
 )
 
 method = QuanticsTensorTrainMethod(
-  quantics=8,
+  quantics=5,
   r₀=0.0,
   rₘₐₓ=8.0,
   l=0,
-  tolerance=1e-9,
-  maxbonddim=32,
+  tolerance=1e-8,
+  maxbonddim=16,
+  maxoperatorbonddim=32,
+  sweeps=4,
 )
 
-result = solve(H, r -> exp(-r^2 / 2), method; info=0)
+result = solve(
+  H,
+  method;
+  initial=r -> exp(-r^2 / 2),
+  nₘₐₓ=2,
+  info=0,
+)
 result.E
 ```
 
-The result approaches the exact ground-state energy ``E=3/2`` as the grid is
-refined. `result.ψ`, `result.u`, and `result.H` contain the compressed radial wave
-function, reduced radial wave function, and reduced Hamiltonian. Their tensor-train
-ranks can be inspected without expanding any object:
+The values approach the exact ``l=0`` oscillator energies ``E_0=3/2`` and
+``E_1=7/2`` as the grid is refined. `result.ψ`, `result.u`, and `result.C` contain
+the compressed radial wave functions, physically normalized reduced radial wave
+functions, and unit-norm DMRG eigenvectors. Their tensor-train ranks can be inspected
+without expanding any object:
+
+```@example qtt
+ranks.(result.C)
+```
+
+The overlap matrix and residual norms provide convergence diagnostics:
+
+```@example qtt
+result.overlaps
+```
 
 ```@example qtt
-ranks(result.ψ)
+result.residuals
 ```
 
 Individual entries can also be contracted directly:
 
 ```@example qtt
-qttvalue(result.ψ, 1)
+qttvalue(result.ψ[1], 1)
 ```
 
 !!! note "Current scope"
     `QuanticsTensorTrainMethod` currently supports `Kinetic`, `RestEnergy`,
     `Constant`, `Linear`, `Coulomb`, `PowerLaw`, `Gaussian`, `Exponential`, `Yukawa`,
-    and `Custom` operators. It evaluates a user-supplied trial wave function; a
-    tensor-train eigensolver is not yet included. `r₀` and `rₘₐₓ` are zero-value
-    (Dirichlet) boundaries and the ``2^q`` grid points lie strictly between them, so
-    choose boundaries where the reduced radial wave function is zero or sufficiently
-    small and perform a grid-convergence study.
+    and `Custom` operators. Excited states use projector deflation, so
+    `deflation_shift` must exceed the relevant spectral gaps; also monitor MPO ranks,
+    overlaps, and residuals. `r₀` and `rₘₐₓ` are zero-value (Dirichlet) boundaries and
+    the ``2^q`` grid points lie strictly between them, so choose boundaries where the
+    reduced radial wave function is zero or sufficiently small and perform a
+    grid-convergence study.
 
 ## Acknowledgments
 
@@ -129,14 +156,16 @@ contributions and support.
 
 The implementation follows the finite-difference QTT construction explored in the
 [CompPhysHack 2026 proof of concept](https://github.com/ohno/CompPhysHack2026Ohno/blob/main/julia/qtt.jl)
-and the method described by Arenstein, Mikkelsen, and Kastoryano in
+and uses the DMRG and tensor-train operations provided by
+[TensorTrainNumerics.jl](https://github.com/MartinMikkelsen/TensorTrainNumerics.jl).
+The QTT construction is described by Arenstein, Mikkelsen, and Kastoryano in
 [Fast and Flexible Quantum-Inspired Differential Equation Solvers with Data Integration](https://doi.org/10.48550/arXiv.2505.17046).
 
 ## API reference
 
 ```@docs; canonical=false
 QuanticsTensorTrainMethod
-solve(hamiltonian::Hamiltonian, wavefunction::Function, method::QuanticsTensorTrainMethod; info=1)
+solve(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod)
 QTTVector
 QTTMatrix
 order
diff --git a/src/QTT.jl b/src/QTT.jl
index 2f5bd76..41283a7 100644
--- a/src/QTT.jl
+++ b/src/QTT.jl
@@ -1,24 +1,21 @@
 export QuanticsTensorTrainMethod, QTTVector, QTTMatrix, order, ranks, qttvalue
 
 import LinearAlgebra
-import QuanticsTCI
-import TensorCrossInterpolation
+import Logging
+import TensorTrainNumerics
 
-"A vector stored as a quantics tensor train with cores `(left rank, 2, right rank)`."
-struct QTTVector{T}
-  cores::Vector{Array{T,3}}
-end
+"A vector stored as a quantics tensor train."
+const QTTVector = TensorTrainNumerics.TTvector
 
-"A matrix stored as a quantics tensor train with cores `(left rank, 2, 2, right rank)`."
-struct QTTMatrix{T}
-  cores::Vector{Array{T,4}}
-end
+"A matrix stored as a quantics tensor train / matrix product operator."
+const QTTMatrix = TensorTrainNumerics.TToperator
 
 "Return the number of binary sites in a QTT."
-order(tt::Union{QTTVector,QTTMatrix}) = length(tt.cores)
+order(tt::Union{QTTVector,QTTMatrix}) = tt.N
 
 "Return the left boundary, internal, and right boundary ranks of a QTT."
-ranks(tt::Union{QTTVector,QTTMatrix}) = [size(tt.cores[1], 1); [size(core, ndims(core)) for core in tt.cores]]
+ranks(tt::QTTVector) = copy(tt.ttv_rks)
+ranks(tt::QTTMatrix) = copy(tt.tto_rks)
 
 struct QuanticsTensorTrainMethod
   quantics::Int
@@ -29,7 +26,10 @@ struct QuanticsTensorTrainMethod
   l::Int
   tolerance::Float64
   maxbonddim::Int
+  maxoperatorbonddim::Int
   maxiter::Int
+  sweeps::Int
+  deflation_shift::Float64
 
   function QuanticsTensorTrainMethod(;
     quantics=10,
@@ -38,59 +38,56 @@ struct QuanticsTensorTrainMethod
     l=0,
     tolerance=1e-10,
     maxbonddim=64,
+    maxoperatorbonddim=2 * maxbonddim,
     maxiter=100,
+    sweeps=4,
+    deflation_shift=10.0,
   )
     quantics >= 2 || throw(ArgumentError("quantics must be at least 2"))
     0 <= r₀ < rₘₐₓ || throw(ArgumentError("the radial interval must satisfy 0 <= r₀ < rₘₐₓ"))
     l >= 0 || throw(ArgumentError("l must be nonnegative"))
     tolerance > 0 || throw(ArgumentError("tolerance must be positive"))
     maxbonddim >= 1 || throw(ArgumentError("maxbonddim must be positive"))
+    maxoperatorbonddim >= 1 || throw(ArgumentError("maxoperatorbonddim must be positive"))
     maxiter >= 1 || throw(ArgumentError("maxiter must be positive"))
+    sweeps >= 1 || throw(ArgumentError("sweeps must be positive"))
+    deflation_shift > 0 || throw(ArgumentError("deflation_shift must be positive"))
 
     npoints = 1 << quantics
     Δr = Float64(rₘₐₓ - r₀) / (npoints + 1)
     R = range(Float64(r₀) + Δr, Float64(rₘₐₓ) - Δr; length=npoints)
-    new(quantics, Float64(r₀), Float64(rₘₐₓ), R, Δr, l,
-        Float64(tolerance), maxbonddim, maxiter)
+    new(
+      quantics,
+      Float64(r₀),
+      Float64(rₘₐₓ),
+      R,
+      Δr,
+      l,
+      Float64(tolerance),
+      maxbonddim,
+      maxoperatorbonddim,
+      maxiter,
+      sweeps,
+      Float64(deflation_shift),
+    )
   end
 end
 
 Base.string(method::QuanticsTensorTrainMethod) =
   "QuanticsTensorTrainMethod(" *
   join(["$(symbol)=$(getproperty(method, symbol))" for symbol in
-        (:quantics, :r₀, :rₘₐₓ, :l, :tolerance, :maxbonddim, :maxiter)], ", ") * ")"
+        (:quantics, :r₀, :rₘₐₓ, :l, :tolerance, :maxbonddim,
+         :maxoperatorbonddim, :maxiter, :sweeps, :deflation_shift)], ", ") * ")"
 Base.show(io::IO, method::QuanticsTensorTrainMethod) = print(io, string(method))
 
-function _check_cores(tt::QTTVector)
-  isempty(tt.cores) && throw(ArgumentError("a QTT must contain at least one core"))
-  size(first(tt.cores), 1) == 1 || throw(ArgumentError("the first left rank must be one"))
-  size(last(tt.cores), 3) == 1 || throw(ArgumentError("the last right rank must be one"))
-  all(size(core, 2) == 2 for core in tt.cores) || throw(ArgumentError("every physical dimension must be two"))
-  all(size(tt.cores[k], 3) == size(tt.cores[k + 1], 1) for k in 1:(order(tt) - 1)) ||
-    throw(ArgumentError("adjacent QTT ranks do not match"))
-  return tt
-end
-
-function _check_cores(tt::QTTMatrix)
-  isempty(tt.cores) && throw(ArgumentError("a QTT must contain at least one core"))
-  size(first(tt.cores), 1) == 1 || throw(ArgumentError("the first left rank must be one"))
-  size(last(tt.cores), 4) == 1 || throw(ArgumentError("the last right rank must be one"))
-  all(size(core, 2) == size(core, 3) == 2 for core in tt.cores) ||
-    throw(ArgumentError("every physical dimension must be two by two"))
-  all(size(tt.cores[k], 4) == size(tt.cores[k + 1], 1) for k in 1:(order(tt) - 1)) ||
-    throw(ArgumentError("adjacent QTT ranks do not match"))
-  return tt
-end
-
 function qttvalue(tt::QTTVector{T}, index::Integer) where {T}
   n = 1 << order(tt)
   1 <= index <= n || throw(BoundsError(1:n, index))
   state = ones(T, 1)
   index0 = index - 1
   for site in 1:order(tt)
-    core = tt.cores[site]
     physical = Int((index0 >> (order(tt) - site)) & 1) + 1
-    state = vec(reshape(state, 1, :) * @view(core[:, physical, :]))
+    state = vec(reshape(state, 1, :) * @view(tt.ttv_vec[site][physical, :, :]))
   end
   return only(state)
 end
@@ -104,179 +101,78 @@ function qttvalue(tt::QTTMatrix{T}, row::Integer, column::Integer) where {T}
   for site in 1:order(tt)
     rowbit = Int((row0 >> (order(tt) - site)) & 1) + 1
     columnbit = Int((column0 >> (order(tt) - site)) & 1) + 1
-    state = vec(reshape(state, 1, :) * @view(tt.cores[site][:, rowbit, columnbit, :]))
+    state = vec(reshape(state, 1, :) * @view(tt.tto_vec[site][rowbit, columnbit, :, :]))
   end
   return only(state)
 end
 
-function Base.:*(factor::Number, tt::QTTVector)
-  cores = copy(tt.cores)
-  cores[1] = factor .* cores[1]
-  return QTTVector(cores)
-end
-
-function Base.:*(factor::Number, tt::QTTMatrix)
-  cores = copy(tt.cores)
-  cores[1] = factor .* cores[1]
-  return QTTMatrix(cores)
-end
-
-Base.:*(tt::Union{QTTVector,QTTMatrix}, factor::Number) = factor * tt
-Base.:-(tt::QTTMatrix) = -1 * tt
-Base.:-(A::QTTMatrix, B::QTTMatrix) = A + (-B)
-
-function Base.:+(A::QTTMatrix{TA}, B::QTTMatrix{TB}) where {TA,TB}
-  order(A) == order(B) || throw(DimensionMismatch("QTT orders must match"))
-  nsites = order(A)
-  T = promote_type(TA, TB)
-  cores = Vector{Array{T,4}}(undef, nsites)
-
-  rankA, rankB = size(A.cores[1], 4), size(B.cores[1], 4)
-  cores[1] = zeros(T, 1, 2, 2, rankA + rankB)
-  cores[1][1, :, :, 1:rankA] = A.cores[1]
-  cores[1][1, :, :, (rankA + 1):end] = B.cores[1]
-
-  for site in 2:(nsites - 1)
-    a, b = A.cores[site], B.cores[site]
-    aleft, aright = size(a, 1), size(a, 4)
-    bleft, bright = size(b, 1), size(b, 4)
-    core = zeros(T, aleft + bleft, 2, 2, aright + bright)
-    core[1:aleft, :, :, 1:aright] = a
-    core[(aleft + 1):end, :, :, (aright + 1):end] = b
-    cores[site] = core
-  end
-
-  a, b = A.cores[end], B.cores[end]
-  aleft, bleft = size(a, 1), size(b, 1)
-  cores[end] = zeros(T, aleft + bleft, 2, 2, 1)
-  cores[end][1:aleft, :, :, 1] = a[:, :, :, 1]
-  cores[end][(aleft + 1):end, :, :, 1] = b[:, :, :, 1]
-  return QTTMatrix(cores)
-end
-
-function _tridiagonal_qtt(quantics::Int, lower::Real, diagonal::Real, upper::Real)
-  identity2 = Matrix{Float64}(LinearAlgebra.I, 2, 2)
-  shift = [0.0 1.0; 0.0 0.0]
-
-  firstcore = zeros(Float64, 1, 2, 2, 3)
-  firstcore[1, :, :, 1] = identity2
-  firstcore[1, :, :, 2] = transpose(shift)
-  firstcore[1, :, :, 3] = shift
-
-  middlecore = zeros(Float64, 3, 2, 2, 3)
-  middlecore[1, :, :, 1] = identity2
-  middlecore[1, :, :, 2] = transpose(shift)
-  middlecore[1, :, :, 3] = shift
-  middlecore[2, :, :, 2] = shift
-  middlecore[3, :, :, 3] = transpose(shift)
-
-  lastcore = zeros(Float64, 3, 2, 2, 1)
-  lastcore[1, :, :, 1] = diagonal * identity2 + upper * shift + lower * transpose(shift)
-  lastcore[2, :, :, 1] = lower * shift
-  lastcore[3, :, :, 1] = upper * transpose(shift)
-
-  cores = Vector{Array{Float64,4}}(undef, quantics)
-  cores[1] = firstcore
-  for site in 2:(quantics - 1)
-    cores[site] = middlecore
-  end
-  cores[end] = lastcore
-  return QTTMatrix(cores)
-end
+function _qttvector(f::Function, method::QuanticsTensorTrainMethod)
+  q = method.quantics
+  domain = [Float64[0, 1] for _ in 1:q]
 
-function _diagonal(tt::QTTVector{T}) where {T}
-  cores = Vector{Array{T,4}}(undef, order(tt))
-  for site in 1:order(tt)
-    vectorcore = tt.cores[site]
-    core = zeros(T, size(vectorcore, 1), 2, 2, size(vectorcore, 3))
-    core[:, 1, 1, :] = vectorcore[:, 1, :]
-    core[:, 2, 2, :] = vectorcore[:, 2, :]
-    cores[site] = core
+  function evaluate(bits::AbstractMatrix)
+    indices = zeros(Int, size(bits, 1))
+    for site in 1:q
+      indices .+= round.(Int, @view(bits[:, site])) .* (1 << (q - site))
+    end
+    return [Float64(f(method.R[index + 1])) for index in indices]
   end
-  return QTTMatrix(cores)
-end
-
-function _qttvector_from_tci(tci, method::QuanticsTensorTrainMethod, f::Function)
-  tensortrain = TensorCrossInterpolation.TensorTrain(tci.tci)
-  tensortrain = TensorCrossInterpolation.reverse(tensortrain)
-  sitetensors = collect(TensorCrossInterpolation.sitetensors(tensortrain))
-  length(sitetensors) == method.quantics || throw(ArgumentError("unexpected tensor-train order"))
-  qtt = _check_cores(QTTVector([Float64.(Array(core)) for core in sitetensors]))
-
-  reversed = QTTVector([permutedims(core, (3, 2, 1)) for core in reverse(qtt.cores)])
-  probes = unique(clamp.([1, 2, length(method.R) ÷ 2, length(method.R) - 1, length(method.R)], 1, length(method.R)))
-  normalerror = sum(abs(qttvalue(qtt, i) - f(method.R[i])) for i in probes)
-  reversederror = sum(abs(qttvalue(reversed, i) - f(method.R[i])) for i in probes)
-  return reversederror < normalerror ? reversed : qtt
-end
 
-function _qttvector(f::Function, method::QuanticsTensorTrainMethod)
-  tci, _, _ = QuanticsTCI.quanticscrossinterpolate(
-    Float64,
-    r -> Float64(f(r)),
-    method.R;
-    tolerance=method.tolerance,
-    maxbonddim=method.maxbonddim,
+  algorithm = TensorTrainNumerics.MaxVol(
+    tol=method.tolerance,
+    rmax=method.maxbonddim,
     maxiter=method.maxiter,
+    verbose=false,
+  )
+  return TensorTrainNumerics.tt_cross(
+    evaluate,
+    domain,
+    algorithm;
+    ranks=min(2, method.maxbonddim),
+    val_size=min(1000, 1 << q),
   )
-  return _qttvector_from_tci(tci, method, f)
 end
 
-function _inner(left::QTTVector{TL}, right::QTTVector{TR}) where {TL,TR}
-  order(left) == order(right) || throw(DimensionMismatch("QTT orders must match"))
-  T = promote_type(TL, TR)
-  environment = ones(T, 1, 1)
-  for site in 1:order(left)
-    lcore, rcore = left.cores[site], right.cores[site]
-    next = zeros(T, size(lcore, 3), size(rcore, 3))
-    for lleft in axes(lcore, 1), rleft in axes(rcore, 1), physical in 1:2,
-        lright in axes(lcore, 3), rright in axes(rcore, 3)
-      next[lright, rright] += environment[lleft, rleft] *
-        conj(lcore[lleft, physical, lright]) * rcore[rleft, physical, rright]
-    end
-    environment = next
-  end
-  return only(environment)
+function _compress_vector(vector::QTTVector, method::QuanticsTensorTrainMethod)
+  compressed = copy(vector)
+  TensorTrainNumerics.tt_compress!(
+    compressed,
+    method.maxbonddim;
+    truncerr=method.tolerance,
+  )
+  return compressed
 end
 
-function _expectation(left::QTTVector{TL}, operator::QTTMatrix{TO}, right::QTTVector{TR}) where {TL,TO,TR}
-  order(left) == order(operator) == order(right) || throw(DimensionMismatch("QTT orders must match"))
-  T = promote_type(TL, TO, TR)
-  environment = ones(T, 1, 1, 1)
-  for site in 1:order(left)
-    lcore, opcore, rcore = left.cores[site], operator.cores[site], right.cores[site]
-    next = zeros(T, size(lcore, 3), size(opcore, 4), size(rcore, 3))
-    for lleft in axes(lcore, 1), oleft in axes(opcore, 1), rleft in axes(rcore, 1),
-        row in 1:2, column in 1:2, lright in axes(lcore, 3),
-        oright in axes(opcore, 4), rright in axes(rcore, 3)
-      next[lright, oright, rright] += environment[lleft, oleft, rleft] *
-        conj(lcore[lleft, row, lright]) * opcore[oleft, row, column, oright] *
-        rcore[rleft, column, rright]
-    end
-    environment = next
-  end
-  return only(environment)
+function _compress_operator(operator::QTTMatrix, method::QuanticsTensorTrainMethod)
+  fused = TensorTrainNumerics.tto_to_ttv(operator)
+  TensorTrainNumerics.tt_compress!(
+    fused,
+    method.maxoperatorbonddim;
+    truncerr=method.tolerance,
+  )
+  return TensorTrainNumerics.ttv_to_tto(fused)
 end
 
 function matrix(term::Kinetic, method::QuanticsTensorTrainMethod)
-  h = method.Δr
-  secondderivative = _tridiagonal_qtt(method.quantics, 1 / h^2, -2 / h^2, 1 / h^2)
-  radialoperator = if method.l == 0
-    secondderivative
-  else
-    centrifugal = method.l * (method.l + 1) * _diagonal(_qttvector(r -> 1 / r^2, method))
-    secondderivative - centrifugal
+  kinetic = (term.hbar^2 / (2 * term.m * method.Δr^2)) *
+    TensorTrainNumerics.Δ(method.quantics)
+  if method.l == 0
+    return kinetic
   end
-  return (-term.hbar^2 / (2 * term.m)) * radialoperator
+  centrifugal = _qttvector(
+    r -> term.hbar^2 * method.l * (method.l + 1) / (2 * term.m * r^2),
+    method,
+  )
+  return kinetic + TensorTrainNumerics.ttv_to_diag_tto(centrifugal)
 end
 
 matrix(term::RestEnergy, method::QuanticsTensorTrainMethod) =
-  _diagonal(_qttvector(_ -> term.m * term.c^2, method))
+  TensorTrainNumerics.ttv_to_diag_tto(_qttvector(_ -> term.m * term.c^2, method))
 
 const QTTPotentialTerm = Union{Constant,Linear,Coulomb,PowerLaw,Gaussian,Exponential,Yukawa,Custom}
 
 matrix(term::QTTPotentialTerm, method::QuanticsTensorTrainMethod) =
-  _diagonal(_qttvector(r -> V(term, r), method))
+  TensorTrainNumerics.ttv_to_diag_tto(_qttvector(r -> V(term, r), method))
 
 matrix(term::PotentialTerm, ::QuanticsTensorTrainMethod) =
   throw(ArgumentError("$(typeof(term)) is not supported by QuanticsTensorTrainMethod"))
@@ -288,69 +184,157 @@ function matrix(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod)
       throw(ArgumentError("$(typeof(term)) is not supported by QuanticsTensorTrainMethod"))
     matrix(term, method)
   end
-  return reduce(+, matrices)
+  return _compress_operator(reduce(+, matrices), method)
 end
 
-function solve(
+_ttnorm(vector::QTTVector) = TensorTrainNumerics.norm(vector)
+_ttdot(left::QTTVector, right::QTTVector) = TensorTrainNumerics.dot(left, right)
+
+function _dmrg_groundstate(operator::QTTMatrix, guess::QTTVector, method::QuanticsTensorTrainMethod)
+  normalized_guess = guess / _ttnorm(guess)
+  history, state, rank_history = Logging.with_logger(Logging.NullLogger()) do
+    TensorTrainNumerics.dmrg_eigsolve(
+      operator,
+      normalized_guess;
+      N=2,
+      tol=method.tolerance,
+      sweep_schedule=[method.sweeps],
+      rmax_schedule=[method.maxbonddim],
+      it_solver=true,
+      linsolv_maxiter=method.maxiter,
+      linsolv_tol=max(sqrt(method.tolerance), 1e-8),
+    )
+  end
+  state = state / _ttnorm(state)
+  return last(history), state, history, rank_history
+end
+
+function _state_guess(initial::Function, level::Int, method::QuanticsTensorTrainMethod)
+  width = method.rₘₐₓ - method.r₀
+  modulation(r) = level == 0 ? 1.0 : cos(level * π * (r - method.r₀) / width)
+  return _qttvector(r -> r * initial(r) * modulation(r), method)
+end
+
+function _solve_qtt(
   hamiltonian::Hamiltonian,
-  wavefunction::Function,
+  initial::Function,
   method::QuanticsTensorTrainMethod;
-  info=1,
+  info=4,
+  nₘₐₓ=4,
 )
+  1 <= nₘₐₓ <= (1 << method.quantics) ||
+    throw(ArgumentError("nₘₐₓ must lie between 1 and the number of grid points"))
+
   H = matrix(hamiltonian, method)
-  unnormalizedψ = _qttvector(wavefunction, method)
-  unnormalizedu = _qttvector(r -> r * wavefunction(r), method)
-  denominator = _inner(unnormalizedu, unnormalizedu)
-  numerator = _expectation(unnormalizedu, H, unnormalizedu)
-  energy = real(numerator / denominator)
-  normalization = sqrt(real(4π * method.Δr * denominator))
-  ψ = (1 / normalization) * unnormalizedψ
-  u = (1 / normalization) * unnormalizedu
+  deflated = H
+  energies = Float64[]
+  states = QTTVector[]
+  histories = Vector{Float64}[]
+  rank_histories = Vector{Int}[]
+  residuals = Float64[]
+
+  for level in 0:(nₘₐₓ - 1)
+    guess = _state_guess(initial, level, method)
+    _, state, history, rank_history = _dmrg_groundstate(deflated, guess, method)
+
+    Hstate = H * state
+    energy = real(_ttdot(state, Hstate) / _ttdot(state, state))
+    residual = _ttnorm(Hstate - energy * state)
+    push!(energies, energy)
+    push!(states, state)
+    push!(histories, history)
+    push!(rank_histories, rank_history)
+    push!(residuals, residual)
+
+    if level < nₘₐₓ - 1
+      projector = TensorTrainNumerics.outer_product(state, state)
+      deflated = _compress_operator(deflated + method.deflation_shift * projector, method)
+    end
+  end
+
+  overlaps = [_ttdot(states[i], states[j]) for i in eachindex(states), j in eachindex(states)]
+  scale = inv(sqrt(4π * method.Δr))
+  u = [scale * state for state in states]
+  inverse_radius = _qttvector(r -> inv(r), method)
+  ψ = [_compress_vector(TensorTrainNumerics.hadamard(state, inverse_radius), method) for state in u]
 
   if info > 0
     println("\n# method\n")
     println(method)
-    println("\n# energy\n")
-    println("E = $energy\n")
+    println("\n# eigenvalues\n")
+    for level in 1:min(nₘₐₓ, info)
+      println("E[$level] = $(energies[level])")
+    end
+    println()
   end
 
   return info >= 0 ? (
     hamiltonian=hamiltonian,
     method=method,
+    nₘₐₓ=nₘₐₓ,
     H=H,
-    E=energy,
+    E=energies,
+    C=states,
     ψ=ψ,
     u=u,
-  ) : (E=energy,)
+    overlaps=overlaps,
+    residuals=residuals,
+    histories=histories,
+    rank_histories=rank_histories,
+  ) : (E=energies,)
+end
+
+function solve(
+  hamiltonian::Hamiltonian,
+  method::QuanticsTensorTrainMethod;
+  initial=r -> exp(-r),
+  info=4,
+  nₘₐₓ=4,
+)
+  return _solve_qtt(hamiltonian, initial, method; info=info, nₘₐₓ=nₘₐₓ)
+end
+
+function solve(
+  hamiltonian::Hamiltonian,
+  wavefunction::Function,
+  method::QuanticsTensorTrainMethod;
+  info=4,
+  nₘₐₓ=4,
+)
+  return _solve_qtt(hamiltonian, wavefunction, method; info=info, nₘₐₓ=nₘₐₓ)
 end
 
 @doc raw"""
     QuanticsTensorTrainMethod(; quantics=10, r₀=0.0, rₘₐₓ=50.0, l=0,
-                               tolerance=1e-10, maxbonddim=64, maxiter=100)
-
-Configure a radial quantics tensor train (QTT) calculation on ``2^q`` uniformly spaced
-interior points, where `q = quantics`. `r₀` and `rₘₐₓ` are zero-value (Dirichlet)
-boundaries and are excluded from the grid, so Coulomb and centrifugal potentials are
-never evaluated at the origin. `tolerance`, `maxbonddim`, and `maxiter` are forwarded
-to tensor cross interpolation.
+                               tolerance=1e-10, maxbonddim=64,
+                               maxoperatorbonddim=128, maxiter=100, sweeps=4,
+                               deflation_shift=10.0)
+
+Configure a radial QTT calculation on ``2^q`` uniformly spaced interior points,
+where `q = quantics`. `r₀` and `rₘₐₓ` are zero-value (Dirichlet) boundaries and are
+excluded from the grid. `maxbonddim` bounds DMRG state ranks;
+`maxoperatorbonddim` bounds the compressed Hamiltonian and deflation-projector ranks.
+`deflation_shift` controls the penalty applied to states already found.
 """ QuanticsTensorTrainMethod
 
 @doc raw"""
-    solve(hamiltonian, wavefunction, method::QuanticsTensorTrainMethod; info=1)
+    solve(hamiltonian, method::QuanticsTensorTrainMethod;
+          initial=r -> exp(-r), info=4, nₘₐₓ=4)
+    solve(hamiltonian, initial::Function, method::QuanticsTensorTrainMethod;
+          info=4, nₘₐₓ=4)
 
-Compress the Hamiltonian and trial radial wave function as QTTs and evaluate the
-Rayleigh quotient for the reduced radial wave function ``u(r)=r\psi(r)``,
+Find the lowest `nₘₐₓ` eigenstates of the reduced radial Hamiltonian with the
+TensorTrainNumerics.jl two-site DMRG eigensolver. Excited states are obtained by
+successive projector deflation,
 
 ```math
-E = \frac{\langle u | H_u | u \rangle}{\langle u | u \rangle},
-\qquad
-H_u = -\frac{\hbar^2}{2m}
-\left[\frac{\mathrm{d}^2}{\mathrm{d}r^2}-\frac{l(l+1)}{r^2}\right]+V(r).
+H^{(n)} = H + \mu \sum_{k=0}^{n-1} |u_k\rangle\langle u_k|.
 ```
 
-The finite-difference kinetic operator is represented directly as a rank-three QTT
-matrix; functions of ``r`` are compressed with tensor cross interpolation.
-""" solve(hamiltonian::Hamiltonian, wavefunction::Function, method::QuanticsTensorTrainMethod; info=1)
+The supplied `initial` function seeds the first state; node-modulated versions seed
+subsequent states. The returned energies are Rayleigh quotients of the original,
+undeflated Hamiltonian. Inspect `overlaps` and `residuals` to assess convergence.
+""" solve(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod; initial=r -> exp(-r), info=4, nₘₐₓ=4)
 
 @doc raw"""
     qttvalue(tt::QTTVector, index)
diff --git a/test/QTT.jl b/test/QTT.jl
index 8336902..f62102e 100644
--- a/test/QTT.jl
+++ b/test/QTT.jl
@@ -1,33 +1,54 @@
 @testset "QTT.jl" begin
   @test_throws ArgumentError QuanticsTensorTrainMethod(quantics=1)
   @test_throws ArgumentError QuanticsTensorTrainMethod(r₀=-0.1)
+  @test_throws ArgumentError QuanticsTensorTrainMethod(deflation_shift=0)
 
   method = QuanticsTensorTrainMethod(
-    quantics=7,
+    quantics=5,
     r₀=0.0,
     rₘₐₓ=8.0,
-    tolerance=1e-9,
-    maxbonddim=32,
+    tolerance=1e-8,
+    maxbonddim=16,
+    maxoperatorbonddim=32,
+    maxiter=100,
+    sweeps=4,
+    deflation_shift=10.0,
   )
   @test length(method.R) == 2^method.quantics
   @test method.R[1] ≈ method.r₀ + method.Δr
   @test method.R[end] ≈ method.rₘₐₓ - method.Δr
 
-  tridiagonal = TwoBody._tridiagonal_qtt(3, -1, 2, -1)
+  kinetic = TwoBody.matrix(Kinetic(hbar=1, m=1), method)
   for row in 1:8, column in 1:8
-    expected = row == column ? 2 : abs(row - column) == 1 ? -1 : 0
-    @test qttvalue(tridiagonal, row, column) ≈ expected atol=1e-12
+    expected = row == column ? 1 / method.Δr^2 :
+      abs(row - column) == 1 ? -1 / (2 * method.Δr^2) : 0
+    @test qttvalue(kinetic, row, column) ≈ expected atol=1e-10
   end
 
   H = Hamiltonian(
     Kinetic(hbar=1, m=1),
     PowerLaw(coefficient=1 / 2, exponent=2),
   )
-  result = solve(H, r -> exp(-r^2 / 2), method; info=0)
-  @test result.E ≈ 1.5 atol=2e-3
-  @test order(result.ψ) == method.quantics
-  @test first(ranks(result.ψ)) == last(ranks(result.ψ)) == 1
-  @test qttvalue(result.ψ, 1) isa Real
+  Random.seed!(1234)
+  result = solve(
+    H,
+    method;
+    initial=r -> exp(-r^2 / 2),
+    nₘₐₓ=2,
+    info=0,
+  )
+  @test result.E ≈ [1.5, 3.5] atol=6e-2
+  @test abs(result.overlaps[1, 2]) < 1e-6
+  @test result.residuals[1] < 1e-5
+  @test result.residuals[2] < 1e-5
+  @test length(result.ψ) == length(result.u) == length(result.C) == 2
+  @test all(order(state) == method.quantics for state in result.C)
+  @test all(first(ranks(state)) == last(ranks(state)) == 1 for state in result.C)
+  @test qttvalue(result.ψ[1], 1) isa Real
+
+  discrete_norm = 4π * method.Δr *
+    sum(abs2(qttvalue(result.u[1], index)) for index in eachindex(method.R))
+  @test discrete_norm ≈ 1 atol=1e-8
 
   potential = TwoBody.matrix(PowerLaw(coefficient=1 / 2, exponent=2), method)
   @test qttvalue(potential, 5, 5) ≈ method.R[5]^2 / 2 atol=1e-7
diff --git a/test/runtests.jl b/test/runtests.jl
index 22a0fd3..138354f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ using Printf
 using Antique
 using SpecialFunctions
 using ForwardDiff
+using Random
 
 @testset verbose = true "TwoBody.jl" begin
   include("Hamiltonian.jl")

From c931c2b91077e9176b8aed11ae8ac31ed27889be Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Tue, 11 Aug 2026 17:58:26 +0900
Subject: [PATCH 3/6] Remove duplicate QTT documentation

---
 docs/src/QTT.md | 42 +++++++++++++-----------------------------
 src/QTT.jl      | 25 +++++++++----------------
 2 files changed, 22 insertions(+), 45 deletions(-)

diff --git a/docs/src/QTT.md b/docs/src/QTT.md
index 4f2ba5a..3dd04bb 100644
--- a/docs/src/QTT.md
+++ b/docs/src/QTT.md
@@ -4,18 +4,19 @@ CurrentModule = TwoBody
 
 # Quantics Tensor Train
 
-The quantics tensor train (QTT) method represents a radial grid of ``2^q`` points as
-``q`` binary tensor sites. A sampled radial function is then approximated by a chain
-of small tensor cores rather than stored as a vector with ``2^q`` entries,
+The quantics tensor train (QTT) solver rewrites a radial grid of ``N=2^q`` points as
+``q`` binary sites and approximates the sampled wave function by a chain of tensor
+cores,
 
 ```math
 f(r_i)=F_{b_1\ldots b_q}
-\simeq G^{(1)}(b_1)G^{(2)}(b_2)\cdots G^{(q)}(b_q).
+\simeq G^{(1)}(b_1)\cdots G^{(q)}(b_q),
+\qquad
+i-1=\sum_{k=1}^{q}b_k2^{q-k}.
 ```
 
-The Hamiltonian is represented in the same form as a matrix product operator (MPO).
-DMRG solves the eigenvalue problem directly in tensor-train form; previously found
-states are shifted out of the low-energy spectrum to expose successive excitations,
+The Hamiltonian is stored as a matrix product operator (MPO). Two-site DMRG finds
+its low-energy states, with projector penalties exposing successive excitations,
 
 ```math
 H^{(n)}_{\mathrm{QTT}}u_n=E_nu_n,
@@ -24,18 +25,8 @@ H^{(n)}_{\mathrm{QTT}}
 =H_{\mathrm{QTT}}+\mu\sum_{k=0}^{n-1}|u_k\rangle\langle u_k|.
 ```
 
-When the tensor-train ranks remain moderate, both storage and contraction costs grow
-with ``q=\log_2 N`` rather than with the full grid size ``N``.
-
 ## Theory
 
-For ``N=2^q`` interior grid points, write the zero-based grid index in binary,
-
-```math
-i-1 = \sum_{k=1}^{q} b_k 2^{q-k},
-\qquad b_k \in \{0,1\}.
-```
-
 Each vector core has dimensions
 ``G^{(k)}\in\mathbb{R}^{\chi_{k-1}\times2\times\chi_k}``, with
 ``\chi_0=\chi_q=1``. A matrix is quantized in both its row and column indices and
@@ -56,8 +47,7 @@ D^{(2)} = \frac{S_- - 2I + S_+}{\Delta r^2}.
 This tridiagonal operator has an exact QTT/MPO representation with maximum bond
 dimension three. Potential functions and initial wave functions are compressed by
 tensor cross interpolation. Two-site DMRG sweeps optimize the tensor cores and adapt
-their bond dimensions. After each state is found, its projector is added to the MPO
-and the result is compressed before the next DMRG solve.
+their bond dimensions; intermediate MPOs are compressed between solves.
 
 If ``\chi`` and ``\rho`` bound the vector and MPO bond dimensions, their storage is
 bounded by
@@ -110,10 +100,8 @@ result.E
 ```
 
 The values approach the exact ``l=0`` oscillator energies ``E_0=3/2`` and
-``E_1=7/2`` as the grid is refined. `result.ψ`, `result.u`, and `result.C` contain
-the compressed radial wave functions, physically normalized reduced radial wave
-functions, and unit-norm DMRG eigenvectors. Their tensor-train ranks can be inspected
-without expanding any object:
+``E_1=7/2`` as the grid is refined. Inspect the eigenvector ranks without expanding
+them:
 
 ```@example qtt
 ranks.(result.C)
@@ -138,12 +126,8 @@ qttvalue(result.ψ[1], 1)
 !!! note "Current scope"
     `QuanticsTensorTrainMethod` currently supports `Kinetic`, `RestEnergy`,
     `Constant`, `Linear`, `Coulomb`, `PowerLaw`, `Gaussian`, `Exponential`, `Yukawa`,
-    and `Custom` operators. Excited states use projector deflation, so
-    `deflation_shift` must exceed the relevant spectral gaps; also monitor MPO ranks,
-    overlaps, and residuals. `r₀` and `rₘₐₓ` are zero-value (Dirichlet) boundaries and
-    the ``2^q`` grid points lie strictly between them, so choose boundaries where the
-    reduced radial wave function is zero or sufficiently small and perform a
-    grid-convergence study.
+    and `Custom` operators. Choose `deflation_shift` above the relevant spectral gaps,
+    monitor MPO ranks, overlaps, and residuals, and perform a grid-convergence study.
 
 ## Acknowledgments
 
diff --git a/src/QTT.jl b/src/QTT.jl
index 41283a7..c15f1f2 100644
--- a/src/QTT.jl
+++ b/src/QTT.jl
@@ -310,11 +310,10 @@ end
                                maxoperatorbonddim=128, maxiter=100, sweeps=4,
                                deflation_shift=10.0)
 
-Configure a radial QTT calculation on ``2^q`` uniformly spaced interior points,
-where `q = quantics`. `r₀` and `rₘₐₓ` are zero-value (Dirichlet) boundaries and are
-excluded from the grid. `maxbonddim` bounds DMRG state ranks;
-`maxoperatorbonddim` bounds the compressed Hamiltonian and deflation-projector ranks.
-`deflation_shift` controls the penalty applied to states already found.
+Configure a radial QTT grid with ``2^{\mathtt{quantics}}`` uniformly spaced interior
+points. The excluded endpoints `r₀` and `rₘₐₓ` impose zero-value (Dirichlet)
+boundaries. `maxbonddim` and `maxoperatorbonddim` bound state and MPO ranks;
+`deflation_shift` is the penalty for each previously computed state.
 """ QuanticsTensorTrainMethod
 
 @doc raw"""
@@ -323,17 +322,11 @@ excluded from the grid. `maxbonddim` bounds DMRG state ranks;
     solve(hamiltonian, initial::Function, method::QuanticsTensorTrainMethod;
           info=4, nₘₐₓ=4)
 
-Find the lowest `nₘₐₓ` eigenstates of the reduced radial Hamiltonian with the
-TensorTrainNumerics.jl two-site DMRG eigensolver. Excited states are obtained by
-successive projector deflation,
-
-```math
-H^{(n)} = H + \mu \sum_{k=0}^{n-1} |u_k\rangle\langle u_k|.
-```
-
-The supplied `initial` function seeds the first state; node-modulated versions seed
-subsequent states. The returned energies are Rayleigh quotients of the original,
-undeflated Hamiltonian. Inspect `overlaps` and `residuals` to assess convergence.
+Compute the lowest `nₘₐₓ` QTT eigenstates. `initial` seeds the first DMRG solve;
+`info` controls progress output. The result contains energies `E`, unit-norm tensor
+trains `C`, radial functions `ψ` and `u`, and the diagnostics `overlaps`, `residuals`,
+`histories`, and `rank_histories`. Each energy is evaluated with the original
+Hamiltonian rather than its deflated counterpart.
 """ solve(hamiltonian::Hamiltonian, method::QuanticsTensorTrainMethod; initial=r -> exp(-r), info=4, nₘₐₓ=4)
 
 @doc raw"""

From 5205a820823884092004627669b6f896b1a82939 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Tue, 11 Aug 2026 18:32:07 +0900
Subject: [PATCH 4/6] Clarify QTT acknowledgments and references

---
 docs/src/QTT.md | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/docs/src/QTT.md b/docs/src/QTT.md
index 3dd04bb..55e118f 100644
--- a/docs/src/QTT.md
+++ b/docs/src/QTT.md
@@ -131,18 +131,19 @@ qttvalue(result.ψ[1], 1)
 
 ## Acknowledgments
 
-This work was developed in collaboration with Lucas Arenstein at
-[CompPhysHack 2026](https://qc-hybrid.github.io/CompPhysHack2026/). The authors
-gratefully acknowledge all those involved in organizing the hackathon for their
+The proof of concept for this QTT solver was developed at
+[CompPhysHack 2026](https://qc-hybrid.github.io/CompPhysHack2026/) in collaboration
+with Lucas Arenstein. We thank him and the hackathon organizers for their
 contributions and support.
 
 ## Bibliography
 
-The implementation follows the finite-difference QTT construction explored in the
-[CompPhysHack 2026 proof of concept](https://github.com/ohno/CompPhysHack2026Ohno/blob/main/julia/qtt.jl)
-and uses the DMRG and tensor-train operations provided by
+The corresponding
+[proof-of-concept implementation](https://github.com/ohno/CompPhysHack2026Ohno/blob/main/julia/qtt.jl)
+is archived in the CompPhysHack repository. The current solver is based on the
+tensor-train operations and two-site DMRG eigensolver provided by
 [TensorTrainNumerics.jl](https://github.com/MartinMikkelsen/TensorTrainNumerics.jl).
-The QTT construction is described by Arenstein, Mikkelsen, and Kastoryano in
+Background on the QTT construction is given by Arenstein, Mikkelsen, and Kastoryano in
 [Fast and Flexible Quantum-Inspired Differential Equation Solvers with Data Integration](https://doi.org/10.48550/arXiv.2505.17046).
 
 ## API reference

From eef047bd7f371c40d737447a49a36162ec41a807 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Tue, 11 Aug 2026 19:14:05 +0900
Subject: [PATCH 5/6] Consolidate QTT acknowledgments and bibliography

---
 docs/src/QTT.md | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/docs/src/QTT.md b/docs/src/QTT.md
index 55e118f..153cc84 100644
--- a/docs/src/QTT.md
+++ b/docs/src/QTT.md
@@ -133,18 +133,19 @@ qttvalue(result.ψ[1], 1)
 
 The proof of concept for this QTT solver was developed at
 [CompPhysHack 2026](https://qc-hybrid.github.io/CompPhysHack2026/) in collaboration
-with Lucas Arenstein. We thank him and the hackathon organizers for their
-contributions and support.
+with Lucas Arenstein; its source is available in the
+[CompPhysHack repository](https://github.com/ohno/CompPhysHack2026Ohno/blob/main/julia/qtt.jl).
+The current solver builds on the tensor-train operations and two-site DMRG
+eigensolver in
+[TensorTrainNumerics.jl](https://github.com/MartinMikkelsen/TensorTrainNumerics.jl).
+We thank Lucas Arenstein, the TensorTrainNumerics.jl developers, and the hackathon
+organizers for their contributions and support.
 
 ## Bibliography
 
-The corresponding
-[proof-of-concept implementation](https://github.com/ohno/CompPhysHack2026Ohno/blob/main/julia/qtt.jl)
-is archived in the CompPhysHack repository. The current solver is based on the
-tensor-train operations and two-site DMRG eigensolver provided by
-[TensorTrainNumerics.jl](https://github.com/MartinMikkelsen/TensorTrainNumerics.jl).
-Background on the QTT construction is given by Arenstein, Mikkelsen, and Kastoryano in
-[Fast and Flexible Quantum-Inspired Differential Equation Solvers with Data Integration](https://doi.org/10.48550/arXiv.2505.17046).
+- L. Arenstein, M. Mikkelsen, and M. Kastoryano,
+  [*Fast and Flexible Quantum-Inspired Differential Equation Solvers with Data Integration*](https://doi.org/10.48550/arXiv.2505.17046),
+  arXiv:2505.17046 (2025).
 
 ## API reference
 

From 29808b5063bdc6ef52a8f6c8bb938d1f8659fcdf Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Wed, 12 Aug 2026 07:08:44 +0900
Subject: [PATCH 6/6] Require TensorTrainNumerics 1.1.4

---
 .github/workflows/CI.yml | 2 +-
 Project.toml             | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 77e2b0b..33c160d 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -24,7 +24,7 @@ jobs:
       matrix:
         version:
           - '1.10'
-          - '1.7'
+          - '1.12'
         os:
           - ubuntu-latest
         arch:
diff --git a/Project.toml b/Project.toml
index f66e0f1..8ad94a0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -31,5 +31,6 @@ Optimisers = "0.2, 0.3, 0.4"
 QuadGK = "2.11"
 SpecialFunctions = "2.3.1"
 Subscripts = "0.1.3"
+TensorTrainNumerics = "1.1.4"
 Zygote = "0.6, 0.7"
-julia = "1.7"
+julia = "1.10"