diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 7e5dfbf04..f7cd0741e 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -33,6 +33,15 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   shorthand: a two-argument call is always the overlap, since the second argument cannot be
   disambiguated between a ket and an operator (undecidable for density matrices, where states
   and operators share a representation). ([#436](https://github.com/QuantumKitHub/MPSKit.jl/pull/436))
+- `transfer_spectrum` now computes the spectrum for all sectors of the transfer space at once,
+  returning a `TensorKit.SectorVector` that can be indexed per sector.
+  The `sector` keyword is removed; a specific selection of sectors (with per-sector counts) can be
+  requested by passing `howmany` as an `AbstractDict`/iterable of `sector => count` pairs.
+  The `below` state and the eigensolver algorithm are optional positional arguments, and the
+  algorithm can be resolved per sector. The Krylov dimension adapts to the number of values requested
+  in each sector, controlled by the new `oversampling` and `oversampling_factor` keywords.
+  Accordingly, `marek_gap` and `correlation_length` now return a `TensorKit.SectorDict` of
+  per-sector results by default; pass `sector = ...` to obtain a single sector's result as before.
 
 ### Deprecated
 
diff --git a/docs/src/examples/classic2d/1.hard-hexagon/figure-1.png b/docs/src/examples/classic2d/1.hard-hexagon/figure-1.png
index 62ee1bbbb..78d2bc359 100644
Binary files a/docs/src/examples/classic2d/1.hard-hexagon/figure-1.png and b/docs/src/examples/classic2d/1.hard-hexagon/figure-1.png differ
diff --git a/docs/src/examples/classic2d/1.hard-hexagon/index.md b/docs/src/examples/classic2d/1.hard-hexagon/index.md
index 96d27fd21..59c72eca4 100644
--- a/docs/src/examples/classic2d/1.hard-hexagon/index.md
+++ b/docs/src/examples/classic2d/1.hard-hexagon/index.md
@@ -12,7 +12,7 @@ EditURL = "../../../../../examples/classic2d/1.hard-hexagon/main.jl"
 
 Tensor networks are a natural way to do statistical mechanics on a lattice.
 As an example of this we will extract the central charge of the hard hexagon model.
-This model is known to have central charge 0.8, and has very peculiar non-local (anyonic) symmetries.
+This model is known to have central charge ``0.8``, and has very peculiar non-local (anyonic) symmetries.
 Because TensorKit supports anyonic symmetries, so does MPSKit.
 To follow the tutorial you need the following packages.
 
@@ -23,7 +23,8 @@ using MPSKit, MPSKitModels, TensorKit, Plots, Polynomials
 The [hard hexagon model](https://en.wikipedia.org/wiki/Hard_hexagon_model) is a 2-dimensional lattice model of a gas, where particles are allowed to be on the vertices of a triangular lattice, but no two particles may be adjacent.
 This can be encoded in a transfer matrix with a local MPO tensor using anyonic symmetries, and the resulting MPO has been implemented in MPSKitModels.
 
-In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry. Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum.
+In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry.
+Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum.
 
 ````julia
 mpo = hard_hexagon()
@@ -53,21 +54,24 @@ V = virtual_space(D)
 ) # use non-hermitian eigensolver
 F = real(expectation_value(ψ, mpo))
 S = real(first(entropy(ψ)))
-ξ = correlation_length(ψ)
+ξ = correlation_length(ψ; sector = leftunit(ψ))
 println("F = $F\tS = $S\tξ = $ξ")
 ````
 
 ````
-F = 0.8839037051703849	S = 1.280782962202701	ξ = 13.849682582739709
+F = 0.8839037051703852	S = 1.2807829621826905	ξ = 13.849682581482702
 
 ````
 
 ## The scaling hypothesis
 
-The dominant eigenvector is of course only an approximation. The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system. This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.
+The dominant eigenvector is of course only an approximation.
+The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system.
+This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.
 
-First we need to know the entropy and correlation length at a bunch of different bond dimensions. Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.
-According to the scaling hypothesis we should have ``S \propto \frac{c}{6} log(ξ)``. Therefore we should find ``c`` using
+First we need to know the entropy and correlation length at a bunch of different bond dimensions.
+Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.
+According to the scaling hypothesis we should have ``S ∝ \frac{c}{6} log(ξ)``. Therefore we should find ``c`` using
 
 ````julia
 function scaling_simulations(
@@ -77,16 +81,17 @@ function scaling_simulations(
     entropies = similar(Ds, Float64)
     correlations = similar(Ds, Float64)
     alg = VUMPS(; verbosity, tol, alg_eigsolve)
+    sector = unit(sectortype(mpo)) # dominant correlation functions are in the trivial sector
 
     ψ, envs, = leading_boundary(ψ₀, mpo, alg)
     entropies[1] = real(entropy(ψ)[1])
-    correlations[1] = correlation_length(ψ)
+    correlations[1] = correlation_length(ψ; sector)
 
     for (i, d) in enumerate(diff(Ds))
         ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme = truncrank(d)), envs)
         ψ, envs, = leading_boundary(ψ, mpo, alg, envs)
         entropies[i + 1] = real(entropy(ψ)[1])
-        correlations[i + 1] = correlation_length(ψ)
+        correlations[i + 1] = correlation_length(ψ; sector)
     end
     return entropies, correlations
 end
@@ -100,7 +105,7 @@ c = f.coeffs[2]
 ````
 
 ````
-0.8025237309498792
+0.802524639544328
 ````
 
 ````julia
diff --git a/docs/src/examples/classic2d/1.hard-hexagon/main.ipynb b/docs/src/examples/classic2d/1.hard-hexagon/main.ipynb
index 0a98e0195..86a70b195 100644
--- a/docs/src/examples/classic2d/1.hard-hexagon/main.ipynb
+++ b/docs/src/examples/classic2d/1.hard-hexagon/main.ipynb
@@ -10,7 +10,7 @@
     "\n",
     "Tensor networks are a natural way to do statistical mechanics on a lattice.\n",
     "As an example of this we will extract the central charge of the hard hexagon model.\n",
-    "This model is known to have central charge 0.8, and has very peculiar non-local (anyonic) symmetries.\n",
+    "This model is known to have central charge $0.8$, and has very peculiar non-local (anyonic) symmetries.\n",
     "Because TensorKit supports anyonic symmetries, so does MPSKit.\n",
     "To follow the tutorial you need the following packages."
    ]
@@ -31,7 +31,8 @@
     "The [hard hexagon model](https://en.wikipedia.org/wiki/Hard_hexagon_model) is a 2-dimensional lattice model of a gas, where particles are allowed to be on the vertices of a triangular lattice, but no two particles may be adjacent.\n",
     "This can be encoded in a transfer matrix with a local MPO tensor using anyonic symmetries, and the resulting MPO has been implemented in MPSKitModels.\n",
     "\n",
-    "In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry. Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum."
+    "In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry.\n",
+    "Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum."
    ]
   },
   {
@@ -77,7 +78,7 @@
     ") # use non-hermitian eigensolver\n",
     "F = real(expectation_value(ψ, mpo))\n",
     "S = real(first(entropy(ψ)))\n",
-    "ξ = correlation_length(ψ)\n",
+    "ξ = correlation_length(ψ; sector = leftunit(ψ))\n",
     "println(\"F = $F\\tS = $S\\tξ = $ξ\")"
    ]
   },
@@ -87,10 +88,13 @@
    "source": [
     "## The scaling hypothesis\n",
     "\n",
-    "The dominant eigenvector is of course only an approximation. The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system. This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.\n",
+    "The dominant eigenvector is of course only an approximation.\n",
+    "The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system.\n",
+    "This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.\n",
     "\n",
-    "First we need to know the entropy and correlation length at a bunch of different bond dimensions. Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.\n",
-    "According to the scaling hypothesis we should have $S \\propto \\frac{c}{6} log(ξ)$. Therefore we should find $c$ using"
+    "First we need to know the entropy and correlation length at a bunch of different bond dimensions.\n",
+    "Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.\n",
+    "According to the scaling hypothesis we should have $S ∝ \\frac{c}{6} log(ξ)$. Therefore we should find $c$ using"
    ]
   },
   {
@@ -106,16 +110,17 @@
     "    entropies = similar(Ds, Float64)\n",
     "    correlations = similar(Ds, Float64)\n",
     "    alg = VUMPS(; verbosity, tol, alg_eigsolve)\n",
+    "    sector = unit(sectortype(mpo)) # dominant correlation functions are in the trivial sector\n",
     "\n",
     "    ψ, envs, = leading_boundary(ψ₀, mpo, alg)\n",
     "    entropies[1] = real(entropy(ψ)[1])\n",
-    "    correlations[1] = correlation_length(ψ)\n",
+    "    correlations[1] = correlation_length(ψ; sector)\n",
     "\n",
     "    for (i, d) in enumerate(diff(Ds))\n",
     "        ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme = truncrank(d)), envs)\n",
     "        ψ, envs, = leading_boundary(ψ, mpo, alg, envs)\n",
     "        entropies[i + 1] = real(entropy(ψ)[1])\n",
-    "        correlations[i + 1] = correlation_length(ψ)\n",
+    "        correlations[i + 1] = correlation_length(ψ; sector)\n",
     "    end\n",
     "    return entropies, correlations\n",
     "end\n",
diff --git a/docs/src/examples/quantum1d/5.haldane-spt/figure-1.png b/docs/src/examples/quantum1d/5.haldane-spt/figure-1.png
index f3641018e..5d8edbbce 100644
Binary files a/docs/src/examples/quantum1d/5.haldane-spt/figure-1.png and b/docs/src/examples/quantum1d/5.haldane-spt/figure-1.png differ
diff --git a/docs/src/examples/quantum1d/5.haldane-spt/figure-2.png b/docs/src/examples/quantum1d/5.haldane-spt/figure-2.png
index 6091af86d..1df134303 100644
Binary files a/docs/src/examples/quantum1d/5.haldane-spt/figure-2.png and b/docs/src/examples/quantum1d/5.haldane-spt/figure-2.png differ
diff --git a/docs/src/examples/quantum1d/5.haldane-spt/index.md b/docs/src/examples/quantum1d/5.haldane-spt/index.md
index 6f95d9769..8fb8dee4d 100644
--- a/docs/src/examples/quantum1d/5.haldane-spt/index.md
+++ b/docs/src/examples/quantum1d/5.haldane-spt/index.md
@@ -128,7 +128,7 @@ E_plus = expectation_value(ψ_plus, H)
 ````
 
 ````
--1.4014193313393009 - 5.2545708620134027e-17im
+-1.4014193313393009 - 3.851708855717825e-17im
 ````
 
 ````julia
@@ -139,7 +139,7 @@ E_minus = expectation_value(ψ_minus, H)
 ````
 
 ````
--1.401483973963084 - 2.875923408269285e-17im
+-1.4014839739630844 - 5.800167584873572e-17im
 ````
 
 ````julia
@@ -188,11 +188,12 @@ println("S_plus = $S_plus")
 ````
 
 ````
-S_minus + log(2) = 1.5486227235421324
-S_plus = 1.5450323530571919
+S_minus + log(2) = 1.548622723541372
+S_plus = 1.5450323530299226
 
 ````
 
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+
diff --git a/docs/src/examples/quantum1d/5.haldane-spt/main.ipynb b/docs/src/examples/quantum1d/5.haldane-spt/main.ipynb
index a5c803849..f85c243c3 100644
--- a/docs/src/examples/quantum1d/5.haldane-spt/main.ipynb
+++ b/docs/src/examples/quantum1d/5.haldane-spt/main.ipynb
@@ -213,4 +213,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 3
-}
+}
\ No newline at end of file
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/figure-1.png b/docs/src/examples/quantum1d/8.bose-hubbard/figure-1.png
index 358ff7c7f..2ab23add4 100644
Binary files a/docs/src/examples/quantum1d/8.bose-hubbard/figure-1.png and b/docs/src/examples/quantum1d/8.bose-hubbard/figure-1.png differ
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/figure-2.png b/docs/src/examples/quantum1d/8.bose-hubbard/figure-2.png
index 85b39fd2f..b8ea6a15c 100644
Binary files a/docs/src/examples/quantum1d/8.bose-hubbard/figure-2.png and b/docs/src/examples/quantum1d/8.bose-hubbard/figure-2.png differ
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/figure-3.png b/docs/src/examples/quantum1d/8.bose-hubbard/figure-3.png
index 85f16f318..523496c7a 100644
Binary files a/docs/src/examples/quantum1d/8.bose-hubbard/figure-3.png and b/docs/src/examples/quantum1d/8.bose-hubbard/figure-3.png differ
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/figure-4.png b/docs/src/examples/quantum1d/8.bose-hubbard/figure-4.png
index 5a88f7b67..addf62c5a 100644
Binary files a/docs/src/examples/quantum1d/8.bose-hubbard/figure-4.png and b/docs/src/examples/quantum1d/8.bose-hubbard/figure-4.png differ
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/figure-5.png b/docs/src/examples/quantum1d/8.bose-hubbard/figure-5.png
index 4048167a9..a33028f2e 100644
Binary files a/docs/src/examples/quantum1d/8.bose-hubbard/figure-5.png and b/docs/src/examples/quantum1d/8.bose-hubbard/figure-5.png differ
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/figure-6.png b/docs/src/examples/quantum1d/8.bose-hubbard/figure-6.png
index 9eec6c12d..805386dc7 100644
Binary files a/docs/src/examples/quantum1d/8.bose-hubbard/figure-6.png and b/docs/src/examples/quantum1d/8.bose-hubbard/figure-6.png differ
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/index.md b/docs/src/examples/quantum1d/8.bose-hubbard/index.md
index 9342b7f9f..0472fcdcb 100644
--- a/docs/src/examples/quantum1d/8.bose-hubbard/index.md
+++ b/docs/src/examples/quantum1d/8.bose-hubbard/index.md
@@ -154,9 +154,9 @@ println("Energy: ", expectation_value(ground_state, hamiltonian))
 ````
 
 ````
-[ Info: VUMPS init:	obj = +4.962471958690e-01	err = 5.8876e-01
-[ Info: VUMPS conv 50:	obj = -6.757777651150e-01	err = 9.9482617748e-07	time = 6.72 sec
-Energy: -0.6757777651149999 - 6.246015340721835e-17im
+[ Info: VUMPS init:	obj = +3.450895756898e-01	err = 6.6224e-01
+[ Info: VUMPS conv 72:	obj = -6.756981551609e-01	err = 9.5394492271e-07	time = 2.91 sec
+Energy: -0.6756981551608794 - 1.967768893425733e-17im
 
 ````
 
@@ -257,7 +257,7 @@ plot!(
 )
 
 scatter!(
-    p, Ds, correlation_length.(states),
+    p, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states),
     ylabel = "Correlation length", xlabel = "Bond dimension",
     xscale = :log10, yscale = :log10,
     inset = bbox(0.2, 0.51, 0.25, 0.25),
@@ -285,13 +285,13 @@ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_o
 
 ````
 7-element Vector{Float64}:
- 0.30974277207656425
- 0.28814775930068737
- 0.2702000951730164
- 0.25712728156142256
- 0.2468538617129866
- 0.2353979140629328
- 0.2279965552408022
+ 0.31098779070601257
+ 0.2881478589434881
+ 0.2702000230423913
+ 0.25712728516508654
+ 0.24685385948652017
+ 0.23539753899204166
+ 0.22799654088348645
 ````
 
 We may now also visualize the momentum distribution function, which is obtained as the
@@ -466,3 +466,4 @@ using what we have learnt in this tutorial.
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+
diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb
index d4ace3b04..615cea4a7 100644
--- a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb
+++ b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb
@@ -1,8 +1,10 @@
 {
  "cells": [
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "using Markdown\n",
     "using MPSKit, MPSKitModels, TensorKit\n",
@@ -11,12 +13,11 @@
     "\n",
     "theme(:wong)\n",
     "default(fontfamily = \"Computer Modern\", label = nothing, dpi = 100, framestyle = :box)"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "# 1D Bose-Hubbard model\n",
     "\n",
@@ -96,21 +97,21 @@
     "`ComplexSpace` (typeset as `\\bbC`). As $D$ is increased,\n",
     "one increases the amount of entanglement, i.e, quantum correlations that can be captured by\n",
     "the state."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "cutoff, D = 4, 5\n",
     "initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "This simply initializes a tensor filled with random entries (check out the documentation for\n",
     "other useful constructors). Next, we need the creation and annihilation operators. While we\n",
@@ -118,72 +119,73 @@
     "[`MPSKitModels.jl`](https://github.com/QuantumKitHub/MPSKitModels.jl) instead that has\n",
     "predefined operators and models for most well-known lattice models. In particular, we can\n",
     "use `MPSKitModels.a_min` to create the bosonic annihilation operator."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "a_op = a_min(cutoff = cutoff) # creates a bosonic annihilation operator without any symmetries\n",
     "display(a_op[])\n",
     "display((a_op' * a_op)[])"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "The [] accessor lets us see the underlying array, and indeed the operators are exactly what\n",
     "we require. Similarly, the Bose Hubbard model is also predefined in\n",
     "`MPSKitModels.bose_hubbard_model` (although we will construct our own variant\n",
     "later on)."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "hamiltonian = bose_hubbard_model(InfiniteChain(1); cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "This has created the Hamiltonian operator as a [matrix product operator](@ref\n",
     "InfiniteMPOHamiltonian) (MPO) which is a convenient form to use in conjunction with MPS.\n",
     "Finally, the ground state optimization may be performed with either `iDMRG` or\n",
     "`VUMPS`. Both should take similar arguments but it is known that VUMPS is typically\n",
     "more efficient for these systems so we proceed with that."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol = 1.0e-6, verbosity = 2, maxiter = 200))\n",
     "println(\"Energy: \", expectation_value(ground_state, hamiltonian))"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "This automatically runs the algorithm until a certain [error measure](@ref\n",
     "MPSKit.calc_galerkin) falls below the specified tolerance or the maximum iterations is\n",
     "reached. Let us wrap all this into a convenient function."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "function get_ground_state(mu, t, cutoff, D; kwargs...)\n",
     "    hamiltonian = bose_hubbard_model(InfiniteChain(); cutoff = cutoff, U = 1, mu = mu, t = t)\n",
@@ -194,53 +196,52 @@
     "end\n",
     "\n",
     "ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500)"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "Now that we have the state, we may compute observables using the `expectation_value`\n",
     "function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a\n",
     "`TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not\n",
     "necessary to specify the indices as it spans the whole lattice. We can now plot the\n",
     "correlation function $\\langle \\hat{a}^{\\dagger}_i \\hat{a}_j\\rangle$."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = \"Site index\", ylabel = \"Correlation function\", yscale = :log10)\n",
     "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "We see that the correlations drop off exponentially, indicating the existence of a gapped\n",
     "Mott insulating phase. Let us now shift our parameters to probe other phases."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500)\n",
     "\n",
     "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw = 2, xlabel = \"Site index\", ylabel = \"Correlation function\", yscale = :log10, xscale = :log10)\n",
     "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "In this case, the correlation function drops off algebraically and eventually saturates as\n",
     "$\\lim_{i \\to \\infty}\\langle\\hat{a}_i^{\\dagger} \\hat{a}_j\\rangle ≈ \\langle \\hat{a}_i^{\\dagger}\\rangle \\langle \\hat{a}_j \\rangle = |\\langle a_i \\rangle|^2 \\neq 0$.\n",
@@ -256,12 +257,13 @@
     "effects. We can see this clearly by increasing the bond dimension. We also see that the\n",
     "correlation length seems to depend algebraically on the bond dimension as expected from\n",
     "finite-entanglement scaling arguments."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "cutoff = 4\n",
     "Ds = 20:5:50\n",
@@ -293,7 +295,7 @@
     ")\n",
     "\n",
     "scatter!(\n",
-    "    p, Ds, correlation_length.(states),\n",
+    "    p, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states),\n",
     "    ylabel = \"Correlation length\", xlabel = \"Bond dimension\",\n",
     "    xscale = :log10, yscale = :log10,\n",
     "    inset = bbox(0.2, 0.51, 0.25, 0.25),\n",
@@ -305,32 +307,31 @@
     "    ylims = [20, 130],\n",
     "    xlims = [15, 60]\n",
     ")"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system,\n",
-    "forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of\n",
+    "forming a Bose-Einstein condensate which introduces erroneous long-distance behaviour of\n",
     "correlation functions. In case of finite bond dimension, it is thus reasonable to associate\n",
     "the finite expectation value of the field operator to the 'quasicondensate' density of the\n",
     "system which vanishes as $D \\to \\infty$."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_op)), states)"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "We may now also visualize the momentum distribution function, which is obtained as the\n",
     "Fourier transform of the single-particle density matrix. Starting from the definition of the\n",
@@ -366,12 +367,13 @@
     "that is not indicative of the true physics of the system. Since we know this contribution\n",
     "vanishes in the infinite bond dimension limit, we instead work with\n",
     "$\\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle_c = \\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle - |\\langle \\hat{a}\\rangle|^2$."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "ks = range(-0.05, 0.15, 500)\n",
     "momentum_distribution = map(\n",
@@ -385,12 +387,11 @@
     ")\n",
     "momentum_distribution = vcat(momentum_distribution...)'\n",
     "plot(ks, momentum_distribution, lab = \"D = \" .* string.(permutedims(Ds)), lw = 1.5, xlabel = \"Momentum k\", ylabel = L\"\\langle n_k \\rangle\", ylim = [0, 50])"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "We see that the density seems to peak around $k=0$, this time seemingly becoming more\n",
     "prominent as $D \\to \\infty$ which seems to suggest again that there is a condensate.\n",
@@ -425,12 +426,13 @@
     "of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs.\n",
     "Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of\n",
     "$e^{\\pm i\\phi}$ in front of the hopping amplitudes."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "function bose_hubbard_model_twisted_bc(\n",
     "        elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial,\n",
@@ -475,12 +477,11 @@
     "superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid\n",
     "\n",
     "superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "Now that we know what phases to expect, we can plot the phase diagram by scanning over a\n",
     "range of parameters. In general, one could do better by performing a bisection algorithm for\n",
@@ -490,12 +491,13 @@
     "using the quasi-condensate density as an order parameter since extracting the superfluid\n",
     "density accurately requires a more robust scheme to compute second derivatives which takes\n",
     "us away from the focus of this tutorial."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "cutoff, D = 4, 10\n",
     "mus = range(0, 0.75, 40)\n",
@@ -512,44 +514,42 @@
     "end\n",
     "\n",
     "heatmap(ts, mus, order_parameters, xlabel = L\"t/U\", ylabel = L\"\\mu/U\", title = L\"\\langle \\hat{a}_i \\rangle\")"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "Although the bond dimension here is quite low, we already see the deformation of the Mott\n",
     "insulator lobes to give way to the well known BKT transition that happens at commensurate\n",
     "density. One can go further and estimate the critical exponents using finite-entanglement\n",
     "scaling procedures on the correlation functions, but these may now be performed with ease\n",
     "using what we have learnt in this tutorial."
-   ],
-   "metadata": {}
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "---\n",
     "\n",
     "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
-   ],
-   "metadata": {}
+   ]
   }
  ],
- "nbformat_minor": 3,
  "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.12.6",
+   "language": "julia",
+   "name": "julia-1.12"
+  },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.12.4"
-  },
-  "kernelspec": {
-   "name": "julia-1.12",
-   "display_name": "Julia 1.12.4",
-   "language": "julia"
+   "version": "1.12.6"
   }
  },
- "nbformat": 4
+ "nbformat": 4,
+ "nbformat_minor": 3
 }
\ No newline at end of file
diff --git a/examples/Cache.toml b/examples/Cache.toml
index 76b7cebb6..1b15ac166 100644
--- a/examples/Cache.toml
+++ b/examples/Cache.toml
@@ -1,12 +1,12 @@
 [classic2d]
-"1.hard-hexagon" = "71d189b9138f98545c4239ed80a922cc4182e934db30af5ba867a0cee57305d0"
+"1.hard-hexagon" = "059f9a5162d75323c7104a8fb178458f965db7bda1481021ae51d4fb8ec9cde9"
 
 [quantum1d]
 "2.haldane" = "c5a0eb70f0930d38053535c659ab39a87121b6ebda040ceadd9deb5f30921315"
 "6.hubbard" = "cef140a6224350345735aac889ee7e33724fc0af9ce94f68be47a7e40107c09c"
 "7.xy-finiteT" = "0f330a157bea739a43a82a937791c680b2fa6e5e479171ee5ef318d1fdae7bcf"
-"8.bose-hubbard" = "a060ac3e29a9169e9420f7125cdb020c83b26a380508ac85f108bff1562bc70c"
+"8.bose-hubbard" = "cf0d9a543e784dc6053e413780d19e1b59bf956dd4598752fdf664804dd68ce6"
 "3.ising-dqpt" = "a2900eed23de7655f600943fae72d1dc35f87f33d11947f707d302ce245399fe"
-"5.haldane-spt" = "28eb0350bd7e721866ccc1c14744e78ae5c650fda81e5f0eefb91de40ee4a76a"
+"5.haldane-spt" = "c8fd3a8d406b9ff3ea9a34b596c5910d81e4eb710b665d7fa04e30f795a46086"
 "4.xxz-heisenberg" = "9033fb104c3f658ccb1b8e335b986b7abda59e0b5359be086d3b62e1cc32ebd6"
 "1.ising-cft" = "3d34757eee7b95120b746c5e30e713203a876a87353ecd9937798f29183c9916"
diff --git a/examples/classic2d/1.hard-hexagon/main.jl b/examples/classic2d/1.hard-hexagon/main.jl
index 3557caa24..3cf117fe6 100644
--- a/examples/classic2d/1.hard-hexagon/main.jl
+++ b/examples/classic2d/1.hard-hexagon/main.jl
@@ -5,7 +5,7 @@ md"""
 
 Tensor networks are a natural way to do statistical mechanics on a lattice.
 As an example of this we will extract the central charge of the hard hexagon model.
-This model is known to have central charge 0.8, and has very peculiar non-local (anyonic) symmetries.
+This model is known to have central charge ``0.8``, and has very peculiar non-local (anyonic) symmetries.
 Because TensorKit supports anyonic symmetries, so does MPSKit.
 To follow the tutorial you need the following packages.
 """
@@ -20,7 +20,8 @@ md"""
 The [hard hexagon model](https://en.wikipedia.org/wiki/Hard_hexagon_model) is a 2-dimensional lattice model of a gas, where particles are allowed to be on the vertices of a triangular lattice, but no two particles may be adjacent.
 This can be encoded in a transfer matrix with a local MPO tensor using anyonic symmetries, and the resulting MPO has been implemented in MPSKitModels.
 
-In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry. Thus, we implement away of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum.
+In order to use these anyonic symmetries, we need to generalise the notion of the bond dimension and define how it interacts with the symmetry.
+Thus, we implement a way of converting integers to symmetric spaces of the given dimension, which provides a crude guess for how the final MPS would distribute its Schmidt spectrum.
 """
 mpo = hard_hexagon()
 P = physicalspace(mpo, 1)
@@ -49,16 +50,19 @@ V = virtual_space(D)
 ) # use non-hermitian eigensolver
 F = real(expectation_value(ψ, mpo))
 S = real(first(entropy(ψ)))
-ξ = correlation_length(ψ)
+ξ = correlation_length(ψ; sector = leftunit(ψ))
 println("F = $F\tS = $S\tξ = $ξ")
 
 md"""
 ## The scaling hypothesis
 
-The dominant eigenvector is of course only an approximation. The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system. This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.
+The dominant eigenvector is of course only an approximation.
+The finite bond dimension enforces a finite correlation length, which effectively introduces a length scale in the system.
+This can be exploited to formulate a scaling hypothesis [pollmann2009](@cite), which in turn allows to extract the central charge.
 
-First we need to know the entropy and correlation length at a bunch of different bond dimensions. Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.
-According to the scaling hypothesis we should have ``S \propto \frac{c}{6} log(ξ)``. Therefore we should find ``c`` using
+First we need to know the entropy and correlation length at several different bond dimensions.
+Our approach will be to re-use the previous approximated dominant eigenvector, and then expanding its bond dimension and re-running VUMPS.
+According to the scaling hypothesis we should have ``S ∝ \frac{c}{6} log(ξ)``. Therefore we should find ``c`` using
 """
 
 function scaling_simulations(
@@ -68,16 +72,17 @@ function scaling_simulations(
     entropies = similar(Ds, Float64)
     correlations = similar(Ds, Float64)
     alg = VUMPS(; verbosity, tol, alg_eigsolve)
+    sector = unit(sectortype(mpo)) # in this example the dominant correlation functions are in the trivial sector
 
     ψ, envs, = leading_boundary(ψ₀, mpo, alg)
     entropies[1] = real(entropy(ψ)[1])
-    correlations[1] = correlation_length(ψ)
+    correlations[1] = correlation_length(ψ; sector)
 
     for (i, d) in enumerate(diff(Ds))
         ψ, envs = changebonds(ψ, mpo, OptimalExpand(; trscheme = truncrank(d)), envs)
         ψ, envs, = leading_boundary(ψ, mpo, alg, envs)
         entropies[i + 1] = real(entropy(ψ)[1])
-        correlations[i + 1] = correlation_length(ψ)
+        correlations[i + 1] = correlation_length(ψ; sector)
     end
     return entropies, correlations
 end
diff --git a/examples/quantum1d/8.bose-hubbard/main.jl b/examples/quantum1d/8.bose-hubbard/main.jl
index 679649303..86edb774d 100644
--- a/examples/quantum1d/8.bose-hubbard/main.jl
+++ b/examples/quantum1d/8.bose-hubbard/main.jl
@@ -213,7 +213,7 @@ plot!(
 )
 
 scatter!(
-    p, Ds, correlation_length.(states),
+    p, Ds, map(ψ -> correlation_length(ψ; sector = leftunit(ψ)), states),
     ylabel = "Correlation length", xlabel = "Bond dimension",
     xscale = :log10, yscale = :log10,
     inset = bbox(0.2, 0.51, 0.25, 0.25),
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index 694dd8024..b7387f3da 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -149,6 +149,7 @@ include("algorithms/derivatives/hamiltonian_derivatives.jl")
 include("algorithms/derivatives/projection_derivatives.jl")
 include("algorithms/expval.jl")
 include("algorithms/toolbox.jl")
+include("algorithms/transfer_spectrum.jl")
 include("algorithms/grassmann.jl")
 include("algorithms/correlators.jl")
 
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index 4b9df6e58..934fdde42 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -66,39 +66,6 @@ function calc_galerkin(below, operator, above, envs)
     return maximum(pos -> calc_galerkin(pos, below, operator, above, envs), eachindex(below))
 end
 
-"""
-    transfer_spectrum(above::InfiniteMPS; below=above, tol=Defaults.tol, num_vals=20,
-                           sector=leftunit(above))
-
-Calculate the partial spectrum of the left mixed transfer matrix corresponding to the
-overlap of a given `above` state and a `below` state. The `sector` keyword argument can be
-used to specify a non-trivial total charge for the transfer matrix eigenvectors.
-Specifically, an auxiliary space `ℂ[typeof(sector)](sector => 1)'` will be added to the
-domain of each eigenvector. The `tol` and `num_vals` keyword arguments are passed to
-`KrylovKit.eigsolve`
-"""
-function transfer_spectrum(
-        above::InfiniteMPS; below = above, tol = Defaults.tol, num_vals = 20,
-        sector = leftunit(above)
-    )
-    init = randomize!(
-        similar(
-            above.AL[1], left_virtualspace(below, 1),
-            spacetype(above)(sector => 1)' * left_virtualspace(above, 1)
-        )
-    )
-
-    transferspace = fuse(left_virtualspace(above, 1) * left_virtualspace(below, 1)')
-    num_vals = min(dim(transferspace, sector), num_vals) # we can ask at most this many values
-    eigenvals, eigenvecs, convhist = eigsolve(
-        flip(TransferMatrix(above.AL, below.AL)), init, num_vals, :LM; tol = tol
-    )
-    convhist.converged < num_vals &&
-        @warn "correlation length failed to converge: normres = $(convhist.normres)"
-
-    return eigenvals
-end
-
 """
     entanglement_spectrum(ψ, site::Int) -> SectorVector{T, sectortype(ψ), AbstractVector{T}}
 
@@ -120,69 +87,6 @@ function entanglement_spectrum(st::FiniteMPS, site::Int)
     return LinearAlgebra.svdvals(st.C[site])
 end
 
-"""
-Find the closest fractions of π, differing at most ```tol_angle```
-"""
-function approx_angles(spectrum; tol_angle = 0.1)
-    angles = angle.(spectrum) ./ π                          # ∈ ]-1, 1]
-    angles_approx = rationalize.(angles, tol = tol_angle)     # ∈ [-1, 1]
-
-    # Remove the effects of the branchcut.
-    angles_approx[findall(angles_approx .== -1)] .= 1       # ∈ ]-1, 1]
-
-    return angles_approx .* π                               # ∈ ]-π, π]
-end
-
-"""
-Given an InfiniteMPS, compute the gap ```ϵ``` for the asymptotics of the transfer matrix, as
-well as the Marek gap ```δ``` as a scaling measure of the bond dimension.
-"""
-function marek_gap(above::InfiniteMPS; tol_angle = 0.1, kwargs...)
-    spectrum = transfer_spectrum(above; kwargs...)
-    return marek_gap(spectrum; tol_angle)
-end
-
-function marek_gap(spectrum::AbstractVector{T}; tol_angle = 0.1) where {T <: Number}
-    # Remove 1s from the spectrum
-    inds = findall(abs.(spectrum) .< 1 - 1.0e-12)
-    length(spectrum) - length(inds) < 2 || @warn "Non-injective mps?"
-
-    spectrum = spectrum[inds]
-
-    angles = approx_angles(spectrum; tol_angle = tol_angle)
-    θ = first(angles)
-
-    spectrum_at_angle = spectrum[findall(angles .== θ)]
-
-    lambdas = -log.(abs.(spectrum_at_angle))
-
-    ϵ = first(lambdas)
-
-    δ = Inf
-    if length(lambdas) > 2
-        δ = lambdas[2] - lambdas[1]
-    end
-
-    return ϵ, δ, θ
-end
-
-"""
-    correlation_length(above::InfiniteMPS; kwargs...)
-
-Compute the correlation length of a given InfiniteMPS based on the next-to-leading
-eigenvalue of the transfer matrix. The `kwargs` are passed to [`transfer_spectrum`](@ref),
-and can for example be used to target the correlation length in a specific sector. 
-"""
-function correlation_length(above::InfiniteMPS; kwargs...)
-    ϵ, = marek_gap(above; kwargs...)
-    return 1 / ϵ
-end
-
-function correlation_length(spectrum::AbstractVector{T}; kwargs...) where {T <: Number}
-    ϵ, = marek_gap(spectrum; kwargs...)
-    return 1 / ϵ
-end
-
 """
     variance(state, hamiltonian, [envs=environments(state, hamiltonian, state)])
 
diff --git a/src/algorithms/transfer_spectrum.jl b/src/algorithms/transfer_spectrum.jl
new file mode 100644
index 000000000..4581ff01f
--- /dev/null
+++ b/src/algorithms/transfer_spectrum.jl
@@ -0,0 +1,179 @@
+"""
+    transfer_spectrum(above::InfiniteMPS, [below = above], [alg]; howmany = 20, kwargs...)
+        -> TensorKit.SectorVector
+
+Calculate the partial spectrum of the left transfer matrix corresponding to the overlap of a given `above` state and a `below` state.
+The result is returned as a `TensorKit.SectorVector`, whose values can be inspected per sector by indexing, i.e. `transfer_spectrum(above)[sector]`.
+
+# Arguments
+
+- `above::InfiniteMPS`: the state for the "above" leg of the mixed transfer matrix.
+- `below::InfiniteMPS=above`: the state for the "below" leg; defaults to the pure transfer matrix of `above`.
+- `alg`: the eigensolver algorithm specification, resolved per sector via [`MatrixAlgebraKit.select_algorithm`](@extref MatrixAlgebraKit.select_algorithm).
+  This can be a KrylovKit algorithm instance (used verbatim for every sector), a `MatrixAlgebraKit.DefaultAlgorithm` or `NamedTuple` bundling keyword arguments,
+  or `nothing` (the default) to construct the eigensolver from the keyword arguments below.
+
+# Keyword arguments
+
+- `howmany = 20`: the number of eigenvalues to compute. This can either be a single `Int`, which is used for every sector of the transfer space,
+  or an `AbstractDict`/iterable of `sector => count` pairs to restrict the computation to specific sectors and request a different number of values per sector.
+- `oversampling = 10`: additive margin on the Krylov dimension beyond the number of values requested in a sector (see `krylovdim` below).
+- `oversampling_factor = 1`: proportionality factor between the Krylov dimension and the number of values requested in a sector (see `krylovdim` below).
+- `krylovdim`: the Krylov dimension of the eigensolver. Unless given explicitly, this is chosen adaptively per sector as
+  `max(Defaults.krylovdim, ceil(Int, oversampling_factor * howmany) + oversampling)`, where `howmany` is the number of values requested in that sector.
+- `kwargs...`: further keyword arguments (e.g. `tol`, `maxiter`) are forwarded to `MatrixAlgebraKit.default_algorithm` to build the eigensolver.
+  Passing eigensolver keyword arguments when `alg` is an algorithm instance is not allowed, and will result in an error.
+"""
+function transfer_spectrum(
+        above::InfiniteMPS, below::InfiniteMPS = above, alg = nothing;
+        howmany = 20, kwargs...
+    )
+    S = TensorKit.check_spacetype(above, below)
+    transferspace = fuse(left_virtualspace(above, 1) * left_virtualspace(below, 1)')
+    howmany = _transfer_howmany(howmany, transferspace)
+
+    T = complex(scalartype(above))
+    Vsp = S(c => n for (c, n) in howmany)
+    spectrum = TensorKit.SectorVector{T}(undef, Vsp)
+
+    A = flip(TransferMatrix(above.AL, below.AL))
+    for (c, n) in howmany
+        algorithm = MatrixAlgebraKit.select_algorithm(
+            transfer_spectrum, above, alg; howmany = n, kwargs...
+        )
+        Va = S(c => 1)
+        Vt = left_virtualspace(above, 1)
+        Vb = left_virtualspace(below, 1)
+        init = randomize!(similar(above.AL[1], Vb ← Va' ⊗ Vt))
+        eigenvals, _, convhist = eigsolve(A, init, n, :LM, algorithm)
+        convhist.converged < n &&
+            @warn "transfer spectrum in sector $c not converged: normres = $(convhist.normres)"
+        spectrum[c] = view(eigenvals, 1:n)
+    end
+    return spectrum
+end
+
+# convenience: positional `alg` without repeating `below`
+transfer_spectrum(above::InfiniteMPS, alg; kwargs...) =
+    transfer_spectrum(above, above, alg; kwargs...)
+
+# normalize the `howmany` specification into a `SectorDict` of per-sector counts,
+# capped by the transfer-space dimension in each sector
+_transfer_howmany(howmany::Int, V) =
+    TensorKit.SectorDict(c => min(howmany, dim(V, c)) for c in sectors(V))
+function _transfer_howmany(howmany, V)
+    return TensorKit.SectorDict(
+        c => min(n, dim(V, c)) for (c, n) in howmany if dim(V, c) > 0
+    )
+end
+
+# `howmany` is intercepted here so it never reaches the algorithm constructor, and is
+# exempt from the no-kwargs-with-`alg` rule (it is injected per sector by `transfer_spectrum`)
+function MatrixAlgebraKit.select_algorithm(
+        ::typeof(transfer_spectrum), above::Union{InfiniteMPS, Type{<:InfiniteMPS}},
+        alg = nothing; howmany = 20, kwargs...
+    )
+    if isnothing(alg)
+        return MatrixAlgebraKit.default_algorithm(transfer_spectrum, above; howmany, kwargs...)
+    elseif alg isa NamedTuple || alg isa Base.Pairs
+        isempty(kwargs) ||
+            throw(ArgumentError("keyword arguments are not allowed when `alg` is a NamedTuple"))
+        return MatrixAlgebraKit.default_algorithm(transfer_spectrum, above; alg..., howmany)
+    elseif alg isa DefaultAlgorithm
+        isempty(kwargs) ||
+            throw(ArgumentError("keyword arguments are not allowed when `alg` is given"))
+        return MatrixAlgebraKit.default_algorithm(transfer_spectrum, above; alg.kwargs..., howmany)
+    elseif alg isa KrylovAlgorithm
+        isempty(kwargs) ||
+            throw(ArgumentError("keyword arguments are not allowed when `alg` is given"))
+        return alg
+    end
+    throw(ArgumentError("Unknown alg $alg"))
+end
+
+function MatrixAlgebraKit.default_algorithm(
+        ::typeof(transfer_spectrum), ::Type{<:InfiniteMPS}; howmany = 20,
+        tol = Defaults.tol, maxiter = Defaults.maxiter,
+        oversampling = 10, oversampling_factor = 1,
+        krylovdim = max(Defaults.krylovdim, ceil(Int, oversampling_factor * howmany) + oversampling),
+        verbosity = Defaults.VERBOSE_NONE, eager = true
+    )
+    return Arnoldi(; tol, maxiter, krylovdim, verbosity, eager)
+end
+
+"""
+Find the closest fractions of π, differing at most ```tol_angle```
+"""
+function approx_angles(spectrum; tol_angle = 0.1)
+    angles = angle.(spectrum) ./ π                          # ∈ ]-1, 1]
+    angles_approx = rationalize.(angles, tol = tol_angle)     # ∈ [-1, 1]
+
+    # Remove the effects of the branchcut.
+    angles_approx[findall(angles_approx .== -1)] .= 1       # ∈ ]-1, 1]
+
+    return angles_approx .* π                               # ∈ ]-π, π]
+end
+
+"""
+    marek_gap(above::InfiniteMPS; sector=nothing, kwargs...)
+
+Compute the gap `ϵ` for the asymptotics of the transfer matrix, as well as the Marek gap
+`δ` as a scaling measure of the bond dimension, along with the associated angle `θ`.
+
+By default this returns a `TensorKit.SectorDict` mapping each sector of the transfer spectrum
+to its `(ϵ, δ, θ)` triplet. Passing a specific `sector` returns only that sector's triplet. The
+remaining `kwargs` are passed on to [`transfer_spectrum`](@ref).
+"""
+function marek_gap(above::InfiniteMPS; sector = nothing, tol_angle = 0.1, howmany = 20, kwargs...)
+    if isnothing(sector)
+        spectrum = transfer_spectrum(above; howmany, kwargs...)
+        return TensorKit.SectorDict(c => marek_gap(spectrum[c]; tol_angle) for c in keys(spectrum))
+    end
+    spectrum = transfer_spectrum(above; howmany = Dict(sector => howmany), kwargs...)
+    return marek_gap(spectrum[sector]; tol_angle)
+end
+
+function marek_gap(spectrum::AbstractVector{T}; tol_angle = 0.1) where {T <: Number}
+    # Remove 1s from the spectrum
+    inds = findall(abs.(spectrum) .< 1 - 1.0e-12)
+    length(spectrum) - length(inds) < 2 || @warn "Non-injective mps?"
+
+    spectrum = spectrum[inds]
+
+    angles = approx_angles(spectrum; tol_angle = tol_angle)
+    θ = first(angles)
+
+    spectrum_at_angle = spectrum[findall(angles .== θ)]
+
+    lambdas = -log.(abs.(spectrum_at_angle))
+
+    ϵ = first(lambdas)
+
+    δ = Inf
+    if length(lambdas) > 2
+        δ = lambdas[2] - lambdas[1]
+    end
+
+    return ϵ, δ, θ
+end
+
+"""
+    correlation_length(above::InfiniteMPS; sector=nothing, kwargs...)
+
+Compute the correlation length of a given InfiniteMPS based on the next-to-leading
+eigenvalue of the transfer matrix.
+
+By default this returns a `TensorKit.SectorDict` mapping each sector of the transfer spectrum
+to its correlation length. Passing a specific `sector` returns only that sector's correlation
+length as a scalar. The remaining `kwargs` are passed on to [`transfer_spectrum`](@ref).
+"""
+function correlation_length(above::InfiniteMPS; sector = nothing, kwargs...)
+    gaps = marek_gap(above; sector, kwargs...)
+    isnothing(sector) || return 1 / first(gaps)                       # gaps :: (ϵ, δ, θ)
+    return TensorKit.SectorDict(c => 1 / first(g) for (c, g) in gaps)  # gaps :: SectorDict
+end
+
+function correlation_length(spectrum::AbstractVector{T}; kwargs...) where {T <: Number}
+    ϵ, = marek_gap(spectrum; kwargs...)
+    return 1 / ϵ
+end
diff --git a/src/utility/plotting.jl b/src/utility/plotting.jl
index d94f90c37..022fee33b 100644
--- a/src/utility/plotting.jl
+++ b/src/utility/plotting.jl
@@ -89,7 +89,7 @@ function entanglementplot end
 end
 
 """
-    transferplot(above, below=above; sectors=[], transferkwargs=(;)[, kwargs...])
+    transferplot(above, below=above; sectors=nothing, transferkwargs=(;)[, kwargs...])
 
 Plot the partial transfer matrix spectrum of two InfiniteMPS's.
 
@@ -98,7 +98,8 @@ Plot the partial transfer matrix spectrum of two InfiniteMPS's.
 - `below::InfiniteMPS=above`: below mps for [`transfer_spectrum`](@ref).
 
 # Keyword Arguments
-- `sectors=[]`: vector of sectors for which to compute the spectrum.
+- `sectors=nothing`: restrict the spectrum to the given sectors; by default all sectors of
+  the transfer space are included.
 - `transferkwargs`: kwargs for call to [`transfer_spectrum`](@ref).
 - `kwargs`: other kwargs are passed on to the plotting backend.
 - `thetaorigin=0`: origin of the angle range.
@@ -118,16 +119,17 @@ function transferplot end
         h::TransferPlot; sectors = nothing, transferkwargs = (;), thetaorigin = 0,
         sector_formatter = string
     )
-    if sectors === nothing
-        sectors = [leftunit(h.args[1])]
+    below = length(h.args) == 1 ? h.args[1] : h.args[2]
+    kwargs = (; transferkwargs...)
+    if sectors !== nothing && get(kwargs, :howmany, 20) isa Int
+        # restrict the computation to the requested sectors
+        howmany = Dict(c => get(kwargs, :howmany, 20) for c in sectors)
+        kwargs = (; kwargs..., howmany)
     end
+    spectra = transfer_spectrum(h.args[1], below; kwargs...)
 
-    for sector in sectors
-        below = length(h.args) == 1 ? h.args[1] : h.args[2]
-        spectrum = transfer_spectrum(
-            h.args[1]; below = below, sector = sector,
-            transferkwargs...
-        )
+    for (sector, spectrum) in pairs(spectra)
+        sectors === nothing || sector in sectors || continue
 
         @series begin
             yguide --> "r"
diff --git a/test/algorithms/correlators.jl b/test/algorithms/correlators.jl
index 52fb914ea..2730a2ccf 100644
--- a/test/algorithms/correlators.jl
+++ b/test/algorithms/correlators.jl
@@ -14,12 +14,12 @@ using TensorKit: ℙ
     ψ = InfiniteMPS([ℙ^2], [ℙ^10])
     H = force_planar(transverse_field_ising())
     ψ, = find_groundstate(ψ, H, VUMPS(; verbosity = 0))
-    len_crit = correlation_length(ψ)[1]
+    len_crit = correlation_length(ψ; sector = leftunit(ψ))
     entrop_crit = entropy(ψ)
 
     H = force_planar(transverse_field_ising(; g = 4))
     ψ, = find_groundstate(ψ, H, VUMPS(; verbosity = 0))
-    len_gapped = correlation_length(ψ)[1]
+    len_gapped = correlation_length(ψ; sector = leftunit(ψ))
     entrop_gapped = entropy(ψ)
 
     @test len_crit > len_gapped
diff --git a/test/algorithms/statmech.jl b/test/algorithms/statmech.jl
index 2101b0447..9a9e8f906 100644
--- a/test/algorithms/statmech.jl
+++ b/test/algorithms/statmech.jl
@@ -94,10 +94,10 @@ using TensorKit: ℙ
         E_tdvp = expectation_value(rho_mps_tdvp, H)
         @test E_tdvp ≈ E_taylor atol = 1.0e-2
 
-        num_vals = 2
-        vals_taylor = @constinferred(transfer_spectrum(convert(InfiniteMPS, rho_taylor); num_vals))
-        vals_mps = @constinferred(transfer_spectrum(rho_mps; num_vals))
-        @test vals_taylor[1:num_vals] ≈ vals_mps[1:num_vals]
+        howmany = 2
+        vals_taylor = @constinferred(transfer_spectrum(convert(InfiniteMPS, rho_taylor); howmany))
+        vals_mps = @constinferred(transfer_spectrum(rho_mps; howmany))
+        @test all(c -> vals_taylor[c] ≈ vals_mps[c], keys(vals_taylor))
     end
 
     @testset "2D infinite partition functions with boundary MPS" verbose = true begin
diff --git a/test/misc/multifusion.jl b/test/misc/multifusion.jl
index 3315a112a..66639ea48 100644
--- a/test/misc/multifusion.jl
+++ b/test/misc/multifusion.jl
@@ -131,8 +131,9 @@ module TestMultifusion
             @test var < 1.0e-8
         end
 
-        @test first(transfer_spectrum(ψ2; sector = C0)) ≈ 1
-        @test !(abs(first(transfer_spectrum(ψ2; sector = C1))) ≈ 1) # testing injectivity
+        spec = transfer_spectrum(ψ2)
+        @test first(spec[C0]) ≈ 1
+        @test !(abs(first(spec[C1])) ≈ 1) # testing injectivity
 
         @test only(keys(entanglement_spectrum(ψ2))) == M