diff --git a/src/ITensorNetworks.jl b/src/ITensorNetworks.jl
index 28da183a..a864d7c1 100644
--- a/src/ITensorNetworks.jl
+++ b/src/ITensorNetworks.jl
@@ -38,6 +38,7 @@ include("solvers/local_solvers/exponentiate.jl")
 include("solvers/local_solvers/dmrg_x.jl")
 include("solvers/local_solvers/contract.jl")
 include("solvers/local_solvers/linsolve.jl")
+include("solvers/local_solvers/expand.jl")
 include("treetensornetworks/abstracttreetensornetwork.jl")
 include("treetensornetworks/treetensornetwork.jl")
 include("treetensornetworks/opsum_to_ttn/matelem.jl")
@@ -51,6 +52,7 @@ include("solvers/solver_utils.jl")
 include("solvers/defaults.jl")
 include("solvers/insert/insert.jl")
 include("solvers/extract/extract.jl")
+include("solvers/subspace_expansions/linalg/standard_svd.jl")
 include("solvers/alternating_update/alternating_update.jl")
 include("solvers/alternating_update/region_update.jl")
 include("solvers/tdvp.jl")
diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl
index b5d315ff..fc3f5652 100644
--- a/src/solvers/defaults.jl
+++ b/src/solvers/defaults.jl
@@ -7,6 +7,8 @@ default_extracter() = default_extracter
 default_inserter() = default_inserter
 default_checkdone() = (; kws...) -> false
 default_transform_operator() = nothing
+default_scale_cutoff_by_timestep(cutoff;time_step, kwargs...) = isinf(time_step) ? cutoff : cutoff / abs(time_step)
+default_scale_maxdim(maxdim;kwargs...) = Int(ceil(2^(1.0/3.0) * maxdim))
 function default_region_printer(;
   cutoff,
   maxdim,
diff --git a/src/solvers/extract/extract.jl b/src/solvers/extract/extract.jl
index feb57c2f..98506ec5 100644
--- a/src/solvers/extract/extract.jl
+++ b/src/solvers/extract/extract.jl
@@ -24,3 +24,25 @@ function default_extracter(state, projected_operator, region, ortho; internal_kw
   projected_operator = position(projected_operator, state, region)
   return state, projected_operator, local_tensor
 end
+
+function extract_and_truncate(state,projected_operator, region, ortho; maxdim=nothing, cutoff=nothing,internal_kwargs)
+  svd_kwargs= (;
+  (isnothing(maxdim) ? (;) : (;maxdim))...,
+  (isnothing(cutoff) ? (;) : (;cutoff))...
+  )
+  state = orthogonalize(state, ortho)
+  if isa(region, AbstractEdge)
+    other_vertex = only(setdiff(support(region), [ortho]))
+    left_inds = uniqueinds(state[ortho], state[other_vertex])
+    #ToDo: replace with call to factorize
+    U, S, V = svd(
+      state[ortho], left_inds; lefttags=tags(state, region), righttags=tags(state, region), svd_kwargs...
+    )
+    state[ortho] = U
+    local_tensor = S * V
+  else
+    local_tensor = prod(state[v] for v in region)
+  end
+  projected_operator = position(projected_operator, state, region)
+  return state, projected_operator, local_tensor
+end
diff --git a/src/solvers/local_solvers/expand.jl b/src/solvers/local_solvers/expand.jl
new file mode 100644
index 00000000..5228768a
--- /dev/null
+++ b/src/solvers/local_solvers/expand.jl
@@ -0,0 +1,206 @@
+#ToDo: is this the correct scaling for infinite timestep? DMRG will not have infinite timestep, unless we set it explicitly
+#ToDo: implement this as a closure, to be constructed in sweep_plan (where time_step is known)
+#scale_cutoff_by_timestep(cutoff;time_step) = isinf(time_step) ? cutoff : cutoff / abs(time_step)
+
+function two_site_expansion_updater(
+    init;
+    state!,
+    projected_operator!,
+    outputlevel,
+    which_sweep,
+    sweep_plan,
+    which_region_update,
+    internal_kwargs,
+    svd_func_expand=ITensorNetworks.rsvd_iterative,
+    maxdim,
+    maxdim_func= (arg;kwargs...) -> identity(arg),
+    cutoff,
+    cutoff_func= (arg;kwargs...) -> identity(arg),
+    use_relative_cutoff=false,
+    use_absolute_cutoff=true,
+  )
+    maxdim=maxdim_func(maxdim;internal_kwargs...)
+    cutoff=cutoff_func(cutoff;internal_kwargs...)
+    region = first(sweep_plan[which_region_update])
+    typeof(region) <: NamedEdge && return init, (;)
+    region = only(region)
+    # figure out next region, since we want to expand there
+    # ToDo: account for non-integer indices into which_region_update
+    next_region = if which_region_update == length(sweep_plan)
+      nothing
+    else
+      first(sweep_plan[which_region_update + 1])
+    end
+    previous_region =
+      which_region_update == 1 ? nothing : first(sweep_plan[which_region_update - 1])
+    isnothing(next_region) && return init, (;)
+    !(typeof(next_region) <: NamedEdge) && return init, (;)
+    !(region == src(next_region) || region == dst(next_region)) && return init, (;)
+    next_vertex = src(next_region) == region ? dst(next_region) : src(next_region)
+  
+    phi, has_changed = _two_site_expand_core(
+      init, region, region => next_vertex; 
+      projected_operator!,
+      state!,
+      svd_func_expand,
+      maxdim,
+      cutoff,
+      use_relative_cutoff,
+      use_absolute_cutoff,
+    )
+    !has_changed && return init, (;)
+    return phi, (;)
+  end
+
+  """
+Returns a function which applies Nullspace projector to an ITensor with matching indices without
+explicitly constructing the projector as an ITensor or Network.
+"""
+function implicit_nullspace(A, linkind)
+  # only works when applied in the direction of the environment tensor, not the other way (due to use of ITensorNetworkMap which is not reversible)
+  outind = uniqueinds(A, linkind)
+  inind = outind  #?
+  # ToDo: implement without ITensorNetworkMap
+  P = ITensors.ITensorNetworkMaps.ITensorNetworkMap(
+    [prime(dag(A), linkind), prime(A, linkind)], inind, outind
+  )
+  return x::ITensor -> x - P(x)
+end
+
+"""
+Performs a local expansion using the two-site projected variance.
+The input state should have orthogonality center on a single vertex (region), and phi0 is that site_tensors.
+Expansion is performed on vertex that would be visited next (if next vertex is a neighbor of region). 
+"""
+function _two_site_expand_core(
+  phi0,
+  region,
+  vertexpair::Pair;
+  projected_operator!,
+  state!,
+  svd_func_expand,
+  cutoff,
+  maxdim,
+  use_relative_cutoff,
+  use_absolute_cutoff,
+)
+  # preliminaries
+  theflux = flux(phi0)
+  v1 = first(vertexpair)
+  v2 = last(vertexpair)
+  verts = [v1, v2]
+  PH = copy(projected_operator![])
+  state = copy(state![])
+  old_linkdim = dim(commonind(state[v1],state[v2]))
+  (old_linkdim >= maxdim) && return phi0, false
+  # orthogonalize on the edge
+  #update of environment not strictly necessary here
+  state, PH, phi = default_extracter(state, PH, edgetype(state)(vertexpair), v1;
+                                  internal_kwargs=(;))
+  psis = map(n -> state[n], verts)  # extract local site tensors
+  linkinds = map(psi -> commonind(psi, phi), psis)
+
+  # compute nullspace to the left and right 
+  #@timeit_debug timer "nullvector" begin
+    nullVecs = map(zip(psis, linkinds)) do (psi, linkind)
+      #return nullspace(psi, linkind; atol=atol)
+      return ITensorNetworks.implicit_nullspace(psi, linkind)
+    end
+  #end
+
+  # update the projected operator
+  #ToDo: may not be necessary if we use the extracter anyway
+  PH = set_nsite(PH, 2)
+  PH = position(PH, state, verts)
+
+  # build environments
+  g = underlying_graph(PH)
+  #@timeit_debug timer "build environments" begin
+    envs = map(zip(verts, psis)) do (n, psi)
+      return noprime(
+        mapreduce(*, [v => n for v in neighbors(g, n) if !(v ∈ verts)]; init=psi) do e
+          return PH.environments[NamedEdge(e)]
+        end * PH.operator[n],
+      )
+    end
+  #end
+
+  # apply the projection into nullspace,
+  # FIXME: think through again whether we really want to evaluate null space projector already
+  envs=[ first(nullVecs)(first(envs)),last(nullVecs)(last(envs))]
+  ininds = uniqueinds(last(psis), phi)
+  outinds = uniqueinds(first(psis), phi)
+  cin = combiner(ininds)
+  cout = combiner(outinds)
+  envMap=[cout * envs[1], phi* (cin * envs[2])]
+
+  # factorize
+  #@timeit_debug timer "svd_func" begin
+  U, S, V = svd_func_expand(
+        envMap,
+        uniqueinds(inds(cout), outinds);
+        maxdim=maxdim - old_linkdim,
+        cutoff=cutoff,
+        use_relative_cutoff,
+        use_absolute_cutoff,
+        )
+  #end
+
+  # catch cases when we decompose a map of norm==0.0
+  (isnothing(U) || iszero(norm(U))) && return phi0, false
+  #FIXME: somehow the svd funcs sometimes return empty ITensors instead of nothing, that should be caught in the SVD routines instead...
+  all(isempty.([U, S, V])) && return phi0, false #ToDo: do not rely on isempty here
+
+  # uncombine indices on the non-link-indices
+  U *= dag(cout)
+  V *= dag(cin)
+
+  # direct sum the site tensors
+  #@timeit_debug timer "direct sum" begin
+    new_psis = map(zip(psis, [U, V])) do (psi, exp_basis)
+      return ITensors.directsum(
+        psi => commonind(psi, phi),
+        exp_basis => uniqueind(exp_basis, psi);
+        tags=tags(commonind(psi, phi)),
+      )
+    end
+  #end
+  new_inds = [last(x) for x in new_psis]
+  new_psis = [first(x) for x in new_psis]
+
+  # extract the expanded linkinds from expanded site tensors and create a zero-tensor
+  phi_indices = replace(
+    inds(phi), (commonind(phi, psis[n]) => dag(new_inds[n]) for n in 1:2)...
+  )
+  if hasqns(phi)
+    new_phi = ITensor(eltype(phi), flux(phi), phi_indices...)
+    #ToDo: Remove? This used to be necessary, but may have been fixed with bugfixes in ITensor
+    #fill!(new_phi, 0.0)
+  else
+    new_phi = ITensor(eltype(phi), phi_indices...)
+  end
+
+  # set the new bond tensor elements from the old bond tensor
+  map(eachindex(phi)) do I
+    v = phi[I]
+    !iszero(v) && (return new_phi[I] = v) # I think this line errors without the fill! with zeros above
+  end
+
+  # apply combiners on linkinds #ToDo: figure out why this is strictly necessary 
+  combiners = map(
+    new_ind -> combiner(new_ind; tags=tags(new_ind), dir=dir(new_ind)), new_inds
+  )
+  for (v, new_psi, C) in zip(verts, new_psis, combiners)
+    state[v] = noprime(new_psi * C)
+  end
+  new_phi = dag(first(combiners)) * new_phi * dag(last(combiners))
+
+  # apply expanded bond tensor to site tensor and reset projected operator to region
+  state[v1] *= new_phi
+  new_phi = state[v1]
+  PH = set_nsite(PH, 1)
+  PH = position(PH, state, [region])
+  projected_operator![] = PH
+  state![] = state
+  return new_phi, true
+end
diff --git a/src/solvers/solver_utils.jl b/src/solvers/solver_utils.jl
index b898fbd9..ad42fea1 100644
--- a/src/solvers/solver_utils.jl
+++ b/src/solvers/solver_utils.jl
@@ -87,3 +87,27 @@ function cache_operator_to_disk(
   end
   return operator
 end
+
+#ToDo: Move? This belongs more into local_solvers
+function compose_updaters(;kwargs...)
+  funcs=values(kwargs)
+  kwarg_symbols=Symbol.(keys(kwargs),"_kwargs")
+  info_symbols=Symbol.(keys(kwargs),"_info")
+
+  function composed_updater(init; kwargs...)
+    info=(;)
+    kwargs_for_updaters=map(x -> kwargs[x], kwarg_symbols)
+    other_kwargs=Base.structdiff((;kwargs...),NamedTuple(kwarg_symbols .=> kwargs_for_updaters))
+  
+    for (func,kwargs_for_updater,info_symbol) in zip(funcs,kwargs_for_updaters,info_symbols)
+      init, new_info=func(init;
+      other_kwargs...,
+      kwargs_for_updater...
+      )
+      #figure out a better way to handle composing the info?
+      info=(;info...,NamedTuple((info_symbol=>new_info,))...)
+    end
+    return init, info
+  end
+  return composed_updater
+end
\ No newline at end of file
diff --git a/src/solvers/subspace_expansions/linalg/standard_svd.jl b/src/solvers/subspace_expansions/linalg/standard_svd.jl
new file mode 100644
index 00000000..3b8b66c0
--- /dev/null
+++ b/src/solvers/subspace_expansions/linalg/standard_svd.jl
@@ -0,0 +1,15 @@
+#ToDo: pass use_*_cutoff as kwargs
+function _svd_solve_normal(envMap, left_ind;
+  maxdim,
+  cutoff,
+  use_relative_cutoff,
+  use_absolute_cutoff)
+  M = ITensors.ITensorNetworkMaps.contract(envMap)
+  norm(M) ≤ eps(real(eltype(M))) && return nothing, nothing, nothing
+  U, S, V = svd(
+    M, left_ind; maxdim, cutoff, use_relative_cutoff, use_absolute_cutoff
+  )
+  vals = diag(array(S))
+  (length(vals) == 1 && vals[1]^2 ≤ cutoff) && return nothing, nothing, nothing
+  return U, S, V
+end
diff --git a/test/test_treetensornetworks/test_solvers/notatest_tdvp_subspace_minimal.jl b/test/test_treetensornetworks/test_solvers/notatest_tdvp_subspace_minimal.jl
new file mode 100644
index 00000000..a273dedf
--- /dev/null
+++ b/test/test_treetensornetworks/test_solvers/notatest_tdvp_subspace_minimal.jl
@@ -0,0 +1,133 @@
+using ITensors
+using ITensorNetworks
+using ITensorNetworks: ModelHamiltonians.tight_binding
+using ITensorNetworks: exponentiate_updater, _svd_solve_normal, ttn, tdvp
+using ITensorNetworks: local_expand_and_exponentiate_updater, vertices
+using ITensorNetworks: two_site_expansion_updater, extract_and_truncate, compose_updaters
+
+#using KrylovKit: exponentiate
+using Observers
+using Random
+using Test
+using ITensorGaussianMPS: hopping_operator
+using LinearAlgebra: exp, I
+
+"""
+Propagate a Gaussian density matrix `c` with a noninteracting Hamiltonian `h` up to time `tau`.
+"""
+function propagate_noninteracting(h::AbstractMatrix, c::AbstractMatrix, tau::Number)
+  U = exp(-im * tau * h)
+  return U * c * conj(transpose(U))
+end
+
+@testset "MPS TDVP" begin
+  @testset "2s vs 1s + local subspace" begin
+    ITensors.enable_auto_fermion()
+    Random.seed!(1234)
+    cutoff = 1e-14
+    N = 20
+    D = 8
+    t = 1.0
+    tp = 0.0
+    g = ITensorNetworks.named_path_graph(N)
+    s = siteinds("Fermion", g; conserve_qns=true)
+    os = tight_binding(g; t, tp)
+    H = ttn(os, s)
+    hmat = hopping_operator(os)
+    # get the exact result
+    tf = 1.0
+    taus = LinRange(0, tf, 2)
+    #init=I(N)[:,StatsBase.sample(1:N, div(N,2); replace=false)]
+    init = I(N)[:, 1:(div(N, 2))] #domain wall initial condition
+    states = i -> i <= div(N, 2) ? "Occ" : "Emp"
+    init = conj(init) * transpose(init)
+    #res=[]
+    res = zeros(N, length(taus))
+    for (i, tau) in enumerate(taus)
+      res[:, i] = real.(diag(propagate_noninteracting(hmat, init, tau))) #densities only
+      #plot!(1:N, real.(diag(last(res))))
+    end
+    psi = ttn(ITensorNetwork(ComplexF64,states,s;link_space=1))
+    psi0 = deepcopy(psi)
+    function my_sweep_printer(;which_sweep,state,kwargs...)
+      m = maximum(ITensorNetworks.edge_data(ITensorNetworks.linkdims(state)))
+      println("Sweep ", which_sweep,", maximum bond dimension: ",m)
+    end
+    ###compose_updater
+    dt = 0.05
+    tdvp_kwargs=( ;
+      updater=compose_updaters(;expander=two_site_expansion_updater,solver=exponentiate_updater),
+      updater_kwargs=(;
+      expander_kwargs=(;
+        cutoff=cutoff,
+        maxdim=D,
+        svd_func_expand=ITensorNetworks._svd_solve_normal,
+        maxdim_func=ITensorNetworks.default_scale_maxdim,
+        cutoff_func=ITensorNetworks.default_scale_cutoff_by_timestep,
+        ),
+      solver_kwargs=(;
+        tol=1E-8
+      )
+      ),
+      extracter=extract_and_truncate,
+      extracter_kwargs=(;maxdim=D,cutoff=cutoff),
+      sweep_printer=my_sweep_printer,
+      time_step=-im * dt,
+      reverse_step=true,
+      order=2,
+      normalize=true,
+      maxdim=D,
+      cutoff=cutoff,
+      outputlevel=1,
+      )
+
+    #=
+    tdvp_kwargs = (
+      time_step=-im * dt,
+      reverse_step=true,
+      order=2,
+      normalize=true,
+      maxdim=D,
+      cutoff=cutoff,
+      outputlevel=1,
+      updater_kwargs=(;
+        expand_kwargs=(;
+          cutoff=cutoff,
+          maxdim=D,
+          svd_func_expand=ITensorNetworks._svd_solve_normal
+        ),
+        exponentiate_kwargs=(;),
+      ),
+    )#,exponentiate_kwargs=(;tol=1e-8)))
+    =#
+    success = false
+    psife = nothing
+    #while !success
+    #  try
+    psife = tdvp(
+      H,
+      -1im * tf,
+      psi;
+      nsites=1,
+      #updater=ITensorNetworks.local_expand_and_exponentiate_updater,
+      tdvp_kwargs...,
+    )
+    #    success=true
+    #  catch
+    #    println("trying again")
+    #  end
+    #  
+    #end
+
+    @test inner(psi0, psi) ≈ 1 atol = 1e-12
+    #@show maxlinkdim(psife)
+    #mag_2s=expect("N",psif)
+    mag_exp = expect("N", psife)
+    # @test real.([mag_2s[v] for v in vertices(psif)]) ≈ res[:,end] atol=1e-3
+    @show maximum(abs.(real.([mag_exp[v] for v in vertices(psife)]) .- res[:, end]))
+    @test all(
+      i -> i < 5e-3, abs.(real.([mag_exp[v] for v in vertices(psife)]) .- res[:, end])
+    )
+  end
+end
+nothing