diff --git a/Project.toml b/Project.toml
index 155e965..c69a6de 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,7 +11,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
-ITensors = "0.3.3"
+ITensors = "0.3.16"
 KrylovKit = "0.5"
 Observers = "0.0.8"
 TimerOutputs = "0.5"
diff --git a/src/ITensorTDVP.jl b/src/ITensorTDVP.jl
index 55498ad..e4d2f68 100644
--- a/src/ITensorTDVP.jl
+++ b/src/ITensorTDVP.jl
@@ -1,14 +1,24 @@
 module ITensorTDVP
 
 using ITensors
-using KrylovKit: exponentiate, eigsolve
+using KrylovKit: exponentiate, eigsolve, svdsolve
 using Printf
+
 using TimerOutputs
 using Observers
 
 using ITensors:
-  AbstractMPS, @debug_check, @timeit_debug, check_hascommoninds, orthocenter, set_nsite!
-
+  AbstractMPS,
+  AbstractProjMPO,
+  @debug_check,
+  @timeit_debug,
+  check_hascommoninds,
+  orthocenter,
+  nullspace,
+  set_nsite!
+using ITensors.NDTensors
+using ITensors.NDTensors: eachdiagblock, blockview
+using ITensors.ITensorNetworkMaps
 # Compatibility of ITensor observer and Observers
 include("update_observer.jl")
 
@@ -25,6 +35,8 @@ include("tdvp_generic.jl")
 include("tdvp.jl")
 include("dmrg.jl")
 include("dmrg_x.jl")
+#include("nullspace.jl")
+include("subspace_expansion.jl")
 
 export tdvp, dmrg_x, to_vec, TimeDependentSum
 
diff --git a/src/subspace_expansion.jl b/src/subspace_expansion.jl
new file mode 100644
index 0000000..8019069
--- /dev/null
+++ b/src/subspace_expansion.jl
@@ -0,0 +1,221 @@
+function replaceind_indval(IV::Tuple, iĩ::Pair)
+  i, ĩ = iĩ
+  return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV))
+end
+
+# atol controls the tolerance cutoff for determining which eigenvectors are in the null
+# space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep
+# the eigenvectors corresponding to eigenvalues of approximately 1.
+
+###the most general implementation may be passing MPS, lims, and bondtensor (defaulting to the Id-matrix)
+###then we can tensorsum both left and right tensor, and just apply the bondtensor to restore input gauge
+###ideally we would also be able to apply tensorsum to the bond tensor instead of that loop below
+"""
+expand subspace (2-site like) in a sweep 
+"""
+function subspace_expansion_sweep!(
+  ψ::MPS, PH::Union{ProjMPO,ProjMPOSum}; maxdim, cutoff, atol=1e-10, kwargs...
+)
+  N = length(ψ)
+
+  if !isortho(ψ) || orthocenter(ψ) != 1
+    orthogonalize!(ψ, 1)
+  end
+
+  PH.nsite = 2
+  nsite = 2
+  position!(PH, ψ, 1)
+
+  for (b, ha) in sweepnext(N; ncenter=2)
+
+    ##TODO: figure out whether these calls should be here or inside subspace expansion, currently we do both?
+    orthogonalize!(ψ, b)
+    position!(PH, ψ, b)
+
+    if (ha == 1 && (b + nsite - 1 != N)) || (ha == 2 && b != 1)
+      b1 = (ha == 1 ? b + 1 : b)
+      Δ = (ha == 1 ? +1 : -1)
+      inds = (ha == 1 ? (b, b + Δ) : (b + Δ, b))
+      subspace_expansion!(
+        ψ, PH, (ψ.llim, ψ.rlim), inds; maxdim, cutoff=cutoff, atol=atol, kwargs...
+      )
+    end
+  end
+  return nothing
+end
+
+function subspace_expansion_sweep(
+  ψ::MPS, PH::Union{ProjMPO,ProjMPOSum}; maxdim, cutoff, atol=1e-10, kwargs...
+)
+  ψc = copy(ψ)
+  PHc = copy(PH)
+  subspace_expansion_sweep!(ψc, PHc; maxdim, cutoff, atol=1e-10, kwargs...)
+  return ψc, PHc
+end
+
+function subspace_expansion!(
+  ψ::MPS,
+  PH,
+  lims::Tuple{Int,Int},
+  b::Tuple{Int,Int};
+  bondtensor=nothing,
+  maxdim,
+  cutoff,
+  atol=1e-10,
+  kwargs...,
+)
+  """
+  expands in nullspace of site-tensors b
+  """
+  # atol refers to the tolerance in nullspace determination (for finite MPS can probably be set rather small)
+  # cutoff refers to largest singular value of gradient (acceleration of population gain of expansion vectors) to keep
+  # this should be related to the cutoff in two-site TDVP: \Delta_rho_i = 0.5 * lambda_y * tau **2 
+  # note that in the initial SVD there is another cutoff parameter `cutoff_compress`, that we set to roughly machine precision for the time being
+  # (allows to move bond-dimension between different partial QN sectors, if no such truncation occurs distribution of bond-dimensions
+  # between different QNs locally is static once bond dimension saturates maxdim.)
+
+  llim, rlim = lims
+  n1, n2 = b
+  PH.nsite = 2
+  old_linkdim = dim(commonind(ψ[n1], ψ[n2]))
+  cutoff_compress = get(kwargs, :cutoff_compress, 1e-12)
+
+  # move orthogonality center to bond
+  # and check whether there are vanishing contributions to the wavefunctions
+  # truncate accordingly
+  # at the cost of unitarity but introducing flexibility of redistributing bond dimensions among QN sectors)
+  if llim == n1
+    @assert rlim == n2 + 1
+    U, S, V = svd(
+      ψ[n2], uniqueinds(ψ[n2], ψ[n1]); maxdim=maxdim, cutoff=cutoff_compress, kwargs...
+    )
+    ϕ_1 = ψ[n1] * V
+    ϕ_2 = U
+    old_linkdim = dim(commonind(U, S))
+    bondtensor = S
+  elseif rlim == n2
+    @assert llim == n1 - 1
+    U, S, V = svd(
+      ψ[n1], uniqueinds(ψ[n1], ψ[n2]); maxdim=maxdim, cutoff=cutoff_compress, kwargs...
+    )
+    ϕ_1 = U
+    ϕ_2 = ψ[n2] * V
+    old_linkdim = dim(commonind(U, S))
+    bondtensor = S
+  end
+
+  # don't expand if we are already at maxdim
+  if old_linkdim >= maxdim
+    println("not expanding")
+    return nothing
+  end
+  position!(PH, ψ, min(n1, n2))
+
+  # orthogonalize(ψ,n1)
+  linkind_l = commonind(ϕ_1, bondtensor)
+  linkind_r = commonind(ϕ_2, bondtensor)
+
+  # compute nullspace to the left and right 
+  NL = nullspace(ϕ_1, linkind_l; atol=atol)
+  NR = nullspace(ϕ_2, linkind_r; atol=atol)
+
+  # if nullspace is empty (happen's for product states with QNs)
+  if norm(NL) == 0.0 || norm(NR) == 0.0
+    return bondtensor
+  end
+
+  # form 2-site wavefunction
+  ϕ = ϕ_1 * bondtensor * ϕ_2
+
+  # get subspace expansion
+  newL, S, newR, success = _subspace_expand_core(
+    ϕ, PH, NL, NR, ; maxdim=maxdim - old_linkdim, cutoff, kwargs...
+  )
+
+  # ToDo: Handle failed subspace expansion properly. Maybe pass an argument whether to raise an error if subspace expansion fails
+  # @assert success
+
+  if success == false
+    println(
+      "Subspace expansion not successful. This may indicate that 2-site TDVP also fails for the given state and Hamiltonian.",
+    )
+    return nothing
+  end
+
+  # expand current site tensors
+  ALⁿ¹, newl = ITensors.directsum(
+    ϕ_1 => uniqueinds(ϕ_1, newL), dag(newL) => uniqueinds(newL, ϕ_1); tags=("Left",)
+  )
+  ARⁿ², newr = ITensors.directsum(
+    ϕ_2 => uniqueinds(ϕ_2, newR), dag(newR) => uniqueinds(newR, ϕ_2); tags=("Right",)
+  )
+
+  # Some checks
+  # ToDo: clean up
+  # ToDo: remove assertions regarding expansion not exceeding maxdim ?
+  @assert (dim(commonind(newL, S)) + old_linkdim) <= maxdim
+  @assert dim(commonind(newL, S)) == dim(commonind(newR, S))
+  @assert(dim(uniqueind(ϕ_1, newL)) + dim(uniqueind(newL, ϕ_1)) == dim(newl))
+  @assert(dim(uniqueind(ϕ_2, newR)) + dim(uniqueind(newR, ϕ_2)) == dim(newr))
+  @assert (old_linkdim + dim(commonind(newL, S))) <= maxdim
+  @assert (old_linkdim + dim(commonind(newR, S))) <= maxdim
+  @assert dim(newl) <= maxdim
+  @assert dim(newr) <= maxdim
+
+  # zero-pad bond-tensor (the orthogonality center)
+  C = ITensor(dag(newl)..., dag(newr)...)
+  ψC = bondtensor
+  for I in eachindex(ψC)
+    v = ψC[I]
+    if !iszero(v)
+      C[I] = ψC[I]
+    end
+  end
+
+  # move orthogonality center back to site (should restore input orthogonality limits)
+  if rlim == n2
+    ψ[n2] = ARⁿ²
+  elseif rlim > n2
+    ψ[n2] = noprime(ARⁿ² * C)
+  end
+  if llim == n1
+    ψ[n1] = ALⁿ¹
+  elseif llim < n1
+    ψ[n1] = noprime(ALⁿ¹ * C)
+  end
+
+  return nothing
+end
+
+function _subspace_expand_core(
+  centerwf::Vector{ITensor}, env, NL, NR; maxdim, cutoff, kwargs...
+)
+  ϕ = ITensor(1.0)
+  for atensor in centerwf
+    ϕ *= atensor
+  end
+  return _subspace_expand_core(ϕ, env, NL, NR; maxdim, cutoff, kwargs...)
+end
+
+function _subspace_expand_core(ϕ::ITensor, env, NL, NR; maxdim, cutoff, kwargs...)
+  ϕH = noprime(env(ϕ))   #add noprime?
+  ϕH = NL * ϕH * NR
+  if norm(ϕH) == 0.0
+    return false, false, false, false
+  end
+
+  U, S, V = svd(
+    ϕH,
+    commoninds(ϕH, NL);
+    maxdim,
+    cutoff,
+    use_relative_cutoff=false,
+    use_absolute_cutoff=true,
+    kwargs...,
+  )
+  @assert dim(commonind(U, S)) <= maxdim
+
+  NL *= dag(U)
+  NR *= dag(V)
+  return NL, S, NR, true
+end
diff --git a/src/tdvp_step.jl b/src/tdvp_step.jl
index 8ce5fb2..5af5dc3 100644
--- a/src/tdvp_step.jl
+++ b/src/tdvp_step.jl
@@ -5,12 +5,27 @@ function tdvp_step(
   sub_time_steps = ITensorTDVP.sub_time_steps(order)
   sub_time_steps *= time_step
   global info
+  nsite::Int = get(kwargs, :nsite, 2)
+  expand::Bool = get(kwargs, :expand, nsite == 2 ? false : (hasqns(psi) ? false : false))
+
+  if expand
+    # do subspace expansion
+    maxdim = get(kwargs, :maxdim, 100)
+    cutoff = get(kwargs, :cutoff, 1e-11)
+    atol = get(kwargs, :atol, 1e-12)
+    # cutoff_trunc is should scale with acceleration of population gain of expansion vectors, thus rescaling with timestep
+    cutoff_trunc = 0.5 * abs(time_step)^2 * cutoff
+    psi, PH = ITensorTDVP.subspace_expansion_sweep(
+      psi, PH; maxdim, cutoff=cutoff_trunc, atol=atol
+    )
+  end
   for substep in 1:length(sub_time_steps)
     psi, PH, info = tdvp_sweep(
       orderings[substep], solver, PH, sub_time_steps[substep], psi; current_time, kwargs...
     )
     current_time += sub_time_steps[substep]
   end
+
   return psi, PH, info
 end
 
diff --git a/test/test_tdvp.jl b/test/test_tdvp.jl
index 1d34fb4..e4c7590 100644
--- a/test/test_tdvp.jl
+++ b/test/test_tdvp.jl
@@ -8,7 +8,7 @@ using Test
 @testset "Basic TDVP" begin
   N = 10
   cutoff = 1e-12
-
+  #maxdim = 100
   s = siteinds("S=1/2", N)
 
   os = OpSum()
@@ -266,6 +266,7 @@ end
     H,
     -im * ttotal,
     phi;
+    nsite=2,
     time_step=-im * tau,
     cutoff,
     normalize=false,
@@ -276,6 +277,102 @@ end
   @test norm(En1 - En2) < 1e-3
 end
 
+@testset "Subspace expansions vs 2-site evolution" begin
+  N = 16
+  cutoff = 1e-12
+  tau = 0.05
+  ttotal = 2.0
+
+  s = siteinds("S=1/2", N; conserve_qns=true)
+
+  os = OpSum()
+  for j in 1:(N - 1)
+    os += 0.5, "S+", j, "S-", j + 1
+    os += 0.5, "S-", j, "S+", j + 1
+    os += "Sz", j, "Sz", j + 1
+  end
+
+  H = MPO(os, s)
+
+  psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
+  phi = copy(psi)
+  c = div(N, 2)
+
+  #
+  # Evolve using TDVP
+  # 
+  struct TDVPObserver <: AbstractObserver end
+
+  Nsteps = convert(Int, ceil(abs(ttotal / tau)))
+  maxdims1 = zeros(Nsteps)
+  maxdims2 = zeros(Nsteps)
+
+  Sz1 = zeros(Nsteps)
+  En1 = zeros(Nsteps)
+  Sz2 = zeros(Nsteps)
+  En2 = zeros(Nsteps)
+
+  function ITensors.measure!(obs::TDVPObserver; sweep, bond, half_sweep, psi, kwargs...)
+    if bond == 1 && half_sweep == 2
+      Sz1[sweep] = expect(psi, "Sz"; sites=c)
+      En1[sweep] = real(inner(psi', H, psi))
+      maxdims1[sweep] = maxlinkdim(psi)
+    end
+  end
+
+  phi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
+
+  phi = tdvp(H, -im * tau, phi; nsite=2, time_step=-im * tau, cutoff, normalize=false)
+
+  psi = deepcopy(phi)
+
+  phi = tdvp(
+    H,
+    -im * ttotal,
+    phi;
+    nsite=2,
+    time_step=-im * tau,
+    cutoff=cutoff,
+    normalize=false,
+    (observer!)=TDVPObserver(),
+  )
+
+  function ITensors.measure!(obs::TDVPObserver; sweep, bond, half_sweep, psi, kwargs...)
+    if bond == 1 && half_sweep == 2
+      Sz2[sweep] = expect(psi, "Sz"; sites=c)
+      En2[sweep] = real(inner(psi', H, psi))
+      maxdims2[sweep] = maxlinkdim(psi)
+    end
+  end
+  psi = tdvp(
+    H,
+    -im * ttotal,
+    psi;
+    time_step=-im * tau,
+    cutoff=cutoff * 1e8,      ##not clear why this is so much larger than the 2-site TDVP cutoff for comparable bond dimensions?
+    cutoff_compress=cutoff,
+    normalize=false,
+    nsite=1,
+    maxdim=100,
+    expand=true,
+    atol=1e-10,
+    #cutoff=5e-2,
+    #atol=1e-11,
+    (observer!)=TDVPObserver(),
+  )
+  #display(En1)
+  #display(En2)
+  #display(Sz1)
+  #display(Sz2)
+  #@show norm(Sz1 - Sz2)
+  #@show norm(En1 - En2)
+  #@show(Sz1)
+  #@show(Sz2)
+
+  @test norm(Sz1 - Sz2) < 1e-3
+  @test norm(En1 - En2) < 1e-3
+end
+
 @testset "Imaginary Time Evolution" for reverse_step in [true, false]
   N = 10
   cutoff = 1e-12