diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3e2823c36..ecb9d746f 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -3,3 +3,4 @@ repos:
     rev: v2.0.1
     hooks:
       - id: runic
+        args: [--docstrings]
diff --git a/docs/DOCSTRING_STYLE.md b/docs/DOCSTRING_STYLE.md
new file mode 100644
index 000000000..f3137133e
--- /dev/null
+++ b/docs/DOCSTRING_STYLE.md
@@ -0,0 +1,230 @@
+# Docstring style guide
+
+This document defines the conventions for docstrings across MPSKit.jl.
+The goal is a single, uniform presentation for all public-facing functionality, so that the API reference reads as one coherent whole.
+
+These rules apply to **every** docstring in `src/`, exported or not.
+Not every symbol needs every section — pick the template that fits (see [Templates](#templates)) — but when a docstring documents arguments, returns, fields, examples, or references, it does so in the one format described here.
+
+## Quick rules
+
+- **Section headers use a single hash** (`# Arguments`, `# Returns`), matching Julia Base.
+- **Leave a blank line after every section header**, before the bullets, prose, or signature block that follows.
+- **Bullet entries** are `` - `name`: description `` — a dash, the name in backticks, a colon, one space, then the description. Add the type (`` `name::Type` ``) only when it is helpful.
+- **Cross-references** use `[`name`](@ref)` for internal symbols and `[`name`](@extref Pkg.name)` for symbols in other packages.
+- **Type parameter lists** put a space after each comma: `Array{T, N}`, `Union{A, B, C}` — never `Array{T,N}`.
+- **Default values** put spaces around the `=`: `tol = 1e-10`, not `tol=1e-10`.
+- **Literature references** use `@cite` keys, never inline DOIs or URLs.
+- **Examples** that show output are runnable `jldoctest` blocks.
+- **Caveats** use `!!! note`; **unstable or experimental** features use `!!! warning`.
+
+## Section headers
+
+Use a single hash (`#`) for all section headers inside a docstring.
+The canonical section names, in the order they should appear, are:
+
+1. `# Constructors` — constructor signatures (container types with non-trivial constructors).
+2. `# Arguments` — positional arguments.
+3. `# Keyword Arguments` — keyword arguments.
+4. `# Returns` — the return value(s).
+5. `# Fields` — the struct fields, when the raw fields are the public API (algorithm structs; rendered by `$(TYPEDFIELDS)`).
+6. `# Properties` — the `getproperty` interface, when it differs from the raw storage (e.g. the gauge views of an MPS container). Use `# Fields` **or** `# Properties`, whichever describes the public surface — not both.
+7. `# Notes` — conventions and caveats worth a dedicated block.
+8. `# Examples` — runnable examples.
+9. `# See also` — related functions; for an algorithm struct, the driver(s) that consume it.
+10. `# References` — literature citations.
+
+Omit any section that does not apply.
+Do not use `## Arguments` (double hash), `# Keywords`, or other spellings.
+Always follow a section header with a blank line, so every section reads the same way:
+
+```
+# Keyword Arguments
+
+- `tol = 1e-10`: convergence tolerance
+```
+
+For docstrings long enough to warrant it, split off detail into a `# Extended help` section (a Julia Base convention) so the summary line and first paragraph stay terse.
+
+## Bullet format
+
+Document arguments, keyword arguments, returns, and manually-listed properties as bullet lists in this form:
+
+```
+- `name`: description
+- `name::Type`: description
+```
+
+A dash (not `*`), the name in backticks, a colon **with no leading space**, one trailing space, then the description.
+
+Include the type in the backtick span **only when it earns its place** — when it constrains what the caller may pass or disambiguates an overloaded name (e.g. `` `O::Union{AbstractMPO, Pair, AbstractTensorMap}` ``).
+Omit it when the type is obvious from the name, the surrounding prose, or the default value (e.g. `` `verbosity`: how much information is displayed ``).
+Never repeat a type that `$(TYPEDFIELDS)` already renders from the struct definition.
+
+Give keyword arguments their default in the backtick span when it is informative, with spaces around the `=`: `` - `tol = 1e-10`: convergence tolerance ``.
+Always put spaces around `=` when writing a default value in a docstring (both in bullet entries and in signature blocks), even where the underlying code omits them.
+Continuation lines of a long description are indented to align under the description text.
+
+## Cross-references and citations
+
+- Internal symbols: `` [`find_groundstate`](@ref) ``.
+- External symbols: `` [`Householder`](@extref MatrixAlgebraKit.Householder) ``.
+  Note the parentheses — `@extref` only expands the `[text](@extref target)` form, not `[text][target]`.
+- Literature: `[Zauner-Stauber et al. Phys. Rev. B 97 (2018)](@cite zauner-stauber2018)`, with the key defined in the bibliography.
+  Do not paste raw DOIs or arXiv links.
+
+## Templates
+
+Three templates cover the whole package.
+Choose by what the symbol is, not by how important it is.
+
+There are two flavours of type docstring — pick by what the type is.
+Algorithm and configuration structs (`A1`) are keyword-configured bags of settings; container/data types (`A2`) hold state and are built through hand-written constructors.
+
+### Template A1 — algorithm and configuration structs
+
+For the keyword-configured `@kwdef` structs: every `Algorithm` subtype, and any similar options struct.
+Each field carries a per-field string literal so that `$(TYPEDFIELDS)` renders the field documentation, including its type — the doc strings themselves do not repeat the type.
+
+```julia
+"""
+$(TYPEDEF)
+
+One paragraph: what the algorithm does and how.
+
+# Fields
+
+$(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`leading_boundary`](@ref).
+
+# References
+
+* [Author et al. Journal (Year)](@cite key)
+"""
+@kwdef struct VUMPS <: Algorithm
+    "tolerance for convergence criterium"
+    tol::Float64 = 1e-10
+    "maximal amount of iterations"
+    maxiter::Int = 200
+end
+```
+
+`$(TYPEDEF)` generates the type signature — do not hand-write it.
+The keyword constructor generated by `@kwdef` *is* the field list, so there is no `# Constructors` section; add one only if the type also offers a non-obvious convenience constructor.
+`# See also` names the driver function(s) that accept the struct, so the algorithm is discoverable from its own page.
+`# References` is optional and only appears when there is literature to cite.
+
+### Template A2 — container and data types
+
+For state-holding types with hand-written constructors: `FiniteMPS`, `InfiniteMPS`, `WindowMPS`, the MPO types, and similar.
+Here the raw fields are internal; document the public `getproperty` interface under `# Properties`, and the constructors explicitly.
+
+```julia
+"""
+$(TYPEDEF)
+
+Type that represents a finite Matrix Product State.
+
+# Constructors
+
+    FiniteMPS([f, eltype], physicalspaces, maxvirtualspaces; kwargs...)
+    FiniteMPS([f, eltype], N, physicalspace, maxvirtualspaces; kwargs...)
+    FiniteMPS(As::Vector{<:GenericMPSTensor}; kwargs...)
+
+Construct an MPS from physical and virtual spaces, or from a list of tensors `As`.
+
+# Arguments
+
+- `As`: vector of site tensors
+- `f = rand`: initializer for tensor data
+- `physicalspaces`: list of physical spaces
+
+# Keyword Arguments
+
+- `normalize = true`: normalize the constructed state
+- `left`: left-most virtual space
+
+# Properties
+
+- `AL`: left-gauged MPS tensors
+- `AR`: right-gauged MPS tensors
+- `AC`: center-gauged MPS tensors
+- `C`: gauge (bond) tensors
+
+# Notes
+
+By convention, `AL[i] * C[i] == AC[i] == C[i-1] * AR[i]`.
+"""
+```
+
+Use `$(TYPEDEF)` for the top line here too, so the type signature never drifts from the definition.
+The constructors are the user-facing interface and are documented separately: stack their signatures as an indented code block under `# Constructors`, then document their parameters in flat sibling `# Arguments` / `# Keyword Arguments` sections (not nested `### ` sub-headers).
+
+### Template B — full-contract functions
+
+Use for the user-facing verbs: `find_groundstate`, `leading_boundary`, `timestep`, `time_evolve`, `expectation_value`, `changebonds`, `approximate`, `correlator`, and the like.
+
+```julia
+"""
+    funcname(ψ₀, H, [environments]; kwargs...) -> (ψ, environments, ϵ)
+
+One paragraph describing the operation.
+
+# Arguments
+
+- `ψ₀::AbstractMPS`: initial guess
+- `H::AbstractMPO`: the operator
+
+# Keyword Arguments
+
+- `tol::Float64 = 1e-10`: convergence tolerance
+
+# Returns
+
+- `ψ::AbstractMPS`: the converged state
+- `ϵ::Float64`: final error estimate
+
+# Examples
+
+```jldoctest
+julia> # runnable example
+```
+
+# References
+
+* [...](@cite key)
+"""
+```
+
+The top line is an indented, four-space signature; stack multiple overloads as separate signature lines.
+Keep the `-> (...)` return annotation on the signature even when a `# Returns` section is present: the signature is the glanceable form, the section is the contract.
+Omit any section that does not apply (a function with no keywords has no `# Keyword Arguments`).
+
+### Template C — lightweight
+
+Use for simple helpers and most internal functions: a signature and a one- or two-sentence description, no sections.
+
+```julia
+"""
+    correlator(ψ, O1, O2, i, j)
+    correlator(ψ, O12, i, j)
+
+Compute the 2-point correlator `⟨ψ|O1[i]O2[j]|ψ⟩`.
+Also accepts a range for `j`.
+"""
+```
+
+## Admonitions
+
+- `!!! note` for caveats and conventions the reader must know (e.g. gauge conventions).
+- `!!! warning` for anything unstable or experimental — everything in `lib/internals`, current GPU support, and any feature that may change.
+
+## Attachment
+
+- Prefer a leading `"""..."""` block directly above the definition.
+- Use `@doc (@doc a) b` only to alias a genuinely identical docstring onto a sibling.
+- A comment between the docstring and the definition silently detaches the docstring; keep them adjacent and put any comment above the docstring.
+- Do not put an HTML comment inside a docstring: DocumenterVitepress escapes it, so `<!-- ... -->` renders as visible body text on the page.
diff --git a/docs/src/assets/mpskit.bib b/docs/src/assets/mpskit.bib
index 9edd3e20b..aa0075e72 100644
--- a/docs/src/assets/mpskit.bib
+++ b/docs/src/assets/mpskit.bib
@@ -129,6 +129,20 @@ @article{devos2022
   abstract = {We calculate the Haldane gap of the SU⁡(3) spin [300] Heisenberg model using variational uniform fully symmetric SU⁡(3) matrix product states, and find that the minimal gap {$\Delta$}/{$J$}=0.0263 is obtained in the [210] sector at momentum 2⁢{$\pi$}/3. We also discuss the symmetry protected topological order of the ground state, and determine the full dispersion relation of the elementary excitations and the correlation lengths of the system.}
 }
 
+@article{gleis2023,
+  title = {Controlled {{Bond Expansion}} for {{Density Matrix Renormalization Group Ground State Search}} at {{Single-Site Costs}}},
+  author = {Gleis, Andreas and Li, Jheng-Wei and {von Delft}, Jan},
+  year = {2023},
+  month = jun,
+  journal = {Physical Review Letters},
+  volume = {130},
+  number = {24},
+  pages = {246402},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevLett.130.246402},
+  url = {https://link.aps.org/doi/10.1103/PhysRevLett.130.246402}
+}
+
 @article{haegeman2011,
   title = {Time-{{Dependent Variational Principle}} for {{Quantum Lattices}}},
   author = {Haegeman, Jutho and Cirac, J. Ignacio and Osborne, Tobias J. and Pi{\v z}orn, Iztok and Verschelde, Henri and Verstraete, Frank},
@@ -931,7 +945,7 @@ @article{zong2026pseudogap
   archiveprefix = {arXiv}
 }
 
-@article{Hubig2015,
+@article{hubig2015,
   title = {Strictly single-site DMRG algorithm with subspace expansion},
   author = {Hubig, C. and McCulloch, I. P. and Schollw\"ock, U. and Wolf, F. A.},
   journal = {Phys. Rev. B},
diff --git a/src/algorithms/ED.jl b/src/algorithms/ED.jl
index 32c0e800a..7b42d6888 100644
--- a/src/algorithms/ED.jl
+++ b/src/algorithms/ED.jl
@@ -1,30 +1,31 @@
 """
-    exact_diagonalization(H::FiniteMPOHamiltonian;
-                          sector=rightunit(H),
-                          len::Int=length(H), num::Int=1, which::Symbol=:SR,
-                          alg=Defaults.alg_eigsolve(; dynamic_tols=false))
-                            -> vals, state_vecs, convhist
+    exact_diagonalization(
+            H::FiniteMPOHamiltonian;
+            sector = rightunit(H), num::Int = 1, which::Symbol = :SR,
+            alg = Defaults.alg_eigsolve(; dynamic_tols = false)
+        ) -> vals, state_vecs, convhist
 
 Use [`KrylovKit.eigsolve`](@extref) to perform exact diagonalization on a
 `FiniteMPOHamiltonian` to find its eigenvectors as `FiniteMPS` of maximal rank, essentially
 equivalent to dense eigenvectors.
 
-### Arguments
+# Arguments
+
 - `H::FiniteMPOHamiltonian`: the Hamiltonian to diagonalize.
 
-### Keyword arguments
-- `sector=rightunit(H)`: the total charge of the
+# Keyword Arguments
+
+- `sector = rightunit(H)`: the total charge of the
   eigenvectors, which is chosen trivial by default.
-- `len::Int=length(H)`: the length of the system.
-- `num::Int=1`: the number of eigenvectors to find.
-- `which::Symbol=:SR`: the kind eigenvalues to find, see [`KrylovKit.eigsolve`](@extref). 
-- `alg=Defaults.alg_eigsolve(; dynamic_tols=false)`: the diagonalization algorithm to use,
+- `num::Int = 1`: the number of eigenvectors to find.
+- `which::Symbol = :SR`: the kind eigenvalues to find, see [`KrylovKit.eigsolve`](@extref).
+- `alg = Defaults.alg_eigsolve(; dynamic_tols = false)`: the diagonalization algorithm to use,
   see [`KrylovKit.eigsolve`](@extref).
 
 !!! note "Valid `sector` values"
     The total charge of the eigenvectors is imposed by adding a charged auxiliary space as
     the leftmost virtualspace of each eigenvector. Specifically, this is achieved by passing
-    `left=Vect[typeof(sector)](sector => 1)` to the [`FiniteMPS`](@ref) constructor. As
+    `left = Vect[typeof(sector)](sector => 1)` to the [`FiniteMPS`](@ref) constructor. As
     such, the only valid `sector` values (i.e. `sector` values for which the corresponding
     eigenstates have valid fusion channels) are those that occur in the dual of the fusion
     of all the physical spaces in the system.
diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl
index 1c699397f..7cfb5dc12 100644
--- a/src/algorithms/approximate/approximate.jl
+++ b/src/algorithms/approximate/approximate.jl
@@ -14,14 +14,16 @@ of an MPS, using initial guess `ψ₀`. If only a state `ψ` is supplied instead
 **Not every algorithm supports every combination of arguments below** — see the per-algorithm
 notes at the end of this docstring before picking one.
 
-## Arguments
+# Arguments
+
 - `ψ₀::AbstractMPS`: initial guess of the approximated state
 - `(O::AbstractMPO, ψ::AbstractMPS)`: operator `O` and state `ψ` to be approximated
 - `ψ::AbstractMPS`: state to be approximated directly (without an operator)
 - `algorithm`: approximation algorithm. See below for a list of available algorithms.
 - `[environments]`: MPS environment manager
 
-## Keywords
+# Keyword Arguments
+
 The keyword-based call (no explicit `algorithm`) is a convenience method that picks an
 algorithm for you based on the type of `ψ₀` (`DMRG`/`DMRG2` for a finite MPS, `VOMPS`/`IDMRG`/
 `IDMRG2` for an infinite MPS) and only accepts the `(O, ψ)` tuple form of `toapprox`. Once you
@@ -33,7 +35,8 @@ struct itself instead (e.g. `DMRG(; tol, maxiter, verbosity)`).
 - `trunc`: if supplied, a truncated two-site sweep (`DMRG2`/`IDMRG2`) is prepended to
   refine the bond dimension before the single-site algorithm polishes the result.
 
-## Algorithms
+# Algorithms
+
 Each algorithm below only supports a subset of the general interface. Check this table before
 picking one — in particular, note that **only `DMRG`/`DMRG2` accept a bare state `ψ`**; the
 infinite algorithms always require an explicit `(O, ψ)` tuple, and **`VOMPS` has no in-place
diff --git a/src/algorithms/approximate/zipup.jl b/src/algorithms/approximate/zipup.jl
index 5ce3ee04d..06b6ec401 100644
--- a/src/algorithms/approximate/zipup.jl
+++ b/src/algorithms/approximate/zipup.jl
@@ -14,14 +14,10 @@ the in-place version simply uses `ψ` as the destination of the sweep, overwriti
 The out-of-place version allocates a destination with the promoted scalar type of `O` and `ϕ`.
 Both return the truncation error `ϵ` alongside the approximated state.
 
-## Fields
+# Constructors
 
-$(TYPEDFIELDS)
-
-## Constructors
-
-    Zipup(; trunc, alg_svd=Defaults.alg_svd(), left_to_right=true)
-    Zipup(alg_zipup, [alg_zipdown]; left_to_right=true)
+    Zipup(; trunc, alg_svd = Defaults.alg_svd(), left_to_right = true)
+    Zipup(alg_zipup, [alg_zipdown]; left_to_right = true)
 
 Create a `Zipup` algorithm with the given truncated gauge algorithm, or by passing a truncation scheme and singular value decomposition algorithm.
 The keyword `trunc` can be either one truncation strategy for a single zip-up sweep, or a tuple `(zipup_trunc, zipdown_trunc)` for a zip-up sweep followed by a zip-down sweep.
@@ -30,10 +26,14 @@ The keyword `left_to_right` selects the direction of the zip-up sweep, the zip-d
 
 Following Paeckel et al., if the desired final bond dimension is `D`, one can use a more permissive zip-up truncation, e.g. rank `2D` with stricter tolerances, and use `alg_zipdown` to impose the final truncation.
 
-## References
+# Fields
+
+$(TYPEDFIELDS)
+
+# References
 
-- [Stoudenmire and White New J. Phys. 12 (2010)](@cite stoudenmire2010)
-- [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019)
+* [Stoudenmire and White New J. Phys. 12 (2010)](@cite stoudenmire2010)
+* [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019)
 """
 struct Zipup{
         U <: MatrixAlgebraKit.TruncatedAlgorithm,
diff --git a/src/algorithms/changebonds/changebonds.jl b/src/algorithms/changebonds/changebonds.jl
index 9dce51c45..1869b347c 100644
--- a/src/algorithms/changebonds/changebonds.jl
+++ b/src/algorithms/changebonds/changebonds.jl
@@ -8,6 +8,35 @@ changedbonds! can modify both the provided state and environments, depending on
 For FiniteMPS, changebonds also modifies the environments.
 
 See also: [`SvdCut`](@ref), [`RandExpand`](@ref), [`VUMPSSvdCut`](@ref), [`OptimalExpand`](@ref)
+
+# Examples
+
+Growing the bond dimension of a product state with [`OptimalExpand`](@ref), which expands
+each bond with directions orthogonal to the current state (using the environments of `H`):
+
+```jldoctest
+julia> Z = TensorMap(Float64[1 0; 0 -1], ℂ^2, ℂ^2);
+
+julia> ψ = FiniteMPS(ones(Float64, (ℂ^2)^4));
+
+julia> H = FiniteMPOHamiltonian(fill(ℂ^2, 4), ((i, i + 1) => Z ⊗ Z for i in 1:3));
+
+julia> dim(left_virtualspace(ψ, 3))
+1
+
+julia> ψ′, envs = changebonds(ψ, H, OptimalExpand(; trunc = truncrank(4)));
+
+julia> dim(left_virtualspace(ψ′, 3))
+2
+```
+
+!!! note
+    A bond is only expanded if there is something to expand it with.
+    If the projection of the two-site update onto the orthogonal complement of the current state
+    vanishes — for instance when the state is already an exact eigenstate of the local terms, or
+    when the operator does not couple into a symmetry sector yet — that bond is left untouched.
+    Replacing `Z ⊗ Z` by `X ⊗ X` above illustrates this: `ones(Float64, (ℂ^2)^4)` is an eigenstate
+    of every `X ⊗ X` term, so every bond stays at dimension 1.
 """ changebonds, changebonds!
 function changebonds end
 function changebonds! end
@@ -16,9 +45,9 @@ function changebonds! end
     changebond(site, dir, ψ, [H], alg, [envs]) -> ψ
     changebond!(site, dir, ψ, [H], alg, [envs]) -> ψ
 
-Expand a single bond of `ψ` in place by adding directions orthogonal to the current state, keeping the state in mixed-canonical form around the enriched bond.
+Expand a single bond of `ψ` by adding directions orthogonal to the current state, keeping the state in mixed-canonical form around the expanded bond.
 The sweep direction `dir` is a `Val(:right)` or `Val(:left)` used for dispatch.
-For `Val(:right)` the bond `(site, site + 1)` is enriched on the right tensor (`ψ.AR[site + 1]`) with zero weight added at `ψ.AC[site]`, so that a subsequent single-site optimization of `site` sees the new directions;
+For `Val(:right)` the bond `(site, site + 1)` is expanded on the right tensor (`ψ.AR[site + 1]`) with zero weight added at `ψ.AC[site]`, so that a subsequent single-site optimization of `site` sees the new directions;
 for `Val(:left)` the mirror is applied to bond `(site - 1, site)`.
 
 See also [`changebonds`](@ref), [`changebonds!`](@ref).
diff --git a/src/algorithms/changebonds/optimalexpand.jl b/src/algorithms/changebonds/optimalexpand.jl
index 9a9c82bc6..ee23d69ea 100644
--- a/src/algorithms/changebonds/optimalexpand.jl
+++ b/src/algorithms/changebonds/optimalexpand.jl
@@ -5,21 +5,40 @@ An algorithm that expands the given mps as described in
 [Zauner-Stauber et al. Phys. Rev. B 97 (2018)](@cite zauner-stauber2018), by selecting the
 dominant contributions of a two-site updated MPS tensor, orthogonal to the original ψ.
 
-The expansion does not alter the state: the added directions are connected through a zero block,
-so that the expanded state is identical to the original one (as required for e.g. TDVP).
+The expansion is state-preserving: the added directions are connected through a zero block,
+so that the expanded state represents the same physical state as the original one (as required
+for e.g. TDVP).
+
+!!! note
+    `trunc` bounds how much is *added* to each bond, not the total bond dimension that is kept.
+    It is applied to the two-site update projected onto the orthogonal complement of the current
+    state, so `truncrank(k)` grows every bond by at most `k` — capped by the dimension of the
+    local two-site complement, which is why a bond can grow by less than `k`, or not at all.
+    The `trunc` of [`SvdCut`](@ref), and of the drivers [`DMRG`](@ref) and [`TDVP`](@ref), has
+    the other meaning: it bounds what is *kept*.
+
+!!! note
+    The projected block is normalized before the decomposition, so a value-based strategy
+    (`trunctol`, `truncerror`) selects a *fraction of the complement weight* rather than an
+    absolute error on the state, and the retained fraction does not shrink as the state
+    converges. `truncrank` and `truncspace` are the strategies with a robust meaning here.
 
 !!! note
     [`changebonds!`](@ref) is only defined for `FiniteMPS`, and modifies both the state and its environment.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`changebonds`](@ref) and [`changebonds!`](@ref).
 """
 @kwdef struct OptimalExpand{S} <: Algorithm
     "algorithm used for the singular value decomposition"
     alg_svd::S = Defaults.alg_svd()
 
-    "algorithm used for truncating the expanded space"
+    "[truncation strategy](@extref MatrixAlgebraKit.TruncationStrategy) selecting how many directions are *added* to each bond, rather than how much of the bond is kept"
     trunc::TruncationStrategy
 end
 
diff --git a/src/algorithms/changebonds/randexpand.jl b/src/algorithms/changebonds/randexpand.jl
index bbbdfd069..ac948042f 100644
--- a/src/algorithms/changebonds/randexpand.jl
+++ b/src/algorithms/changebonds/randexpand.jl
@@ -3,25 +3,34 @@ $(TYPEDEF)
 
 An algorithm that expands the bond dimension by adding random unitary vectors that are
 orthogonal to the existing state. This means that additional directions are added to
-`AL` and `AR` that are contained in the nullspace of both. Note that this is happens in
+`AL` and `AR` that are contained in the nullspace of both. Note that this happens in
 parallel, and therefore the expansion will never go beyond the local two-site subspace.
 
-The truncation strategy dictates the number of expanded states, by generating uniformly
-distributed weights for each state in the two-site space and truncating that.
+`trunc` bounds how much is *added* to each bond, not the total bond dimension that is kept, and
+it acts on a spectrum that carries no physical information: for an `InfiniteMPS` the weights are
+drawn uniformly at random, one per candidate direction, while for a `FiniteMPS` they are the
+singular values of a randomized two-site update restricted to the orthogonal complement. Only
+`truncrank` and `truncspace` therefore have a robust meaning — `trunctol(; atol = x)` keeps the
+directions whose *random* weight happens to exceed `x`. The `trunc` of [`SvdCut`](@ref), and of
+the drivers [`DMRG`](@ref) and [`TDVP`](@ref), has the other meaning: it bounds what is *kept*.
 
 !!! note
     The environments are not used here, but [`changebonds!`](@ref) modifies both the state
     and environment so they remain consistent.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`changebonds`](@ref) and [`changebonds!`](@ref).
 """
 @kwdef struct RandExpand{S} <: Algorithm
     "algorithm used for the singular value decomposition"
     alg_svd::S = Defaults.alg_svd()
 
-    "algorithm used for [truncation](@extref MatrixAlgebraKit.TruncationStrategy) of the expanded space"
+    "[truncation strategy](@extref MatrixAlgebraKit.TruncationStrategy) selecting how many directions are *added* to each bond, rather than how much of the bond is kept"
     trunc::TruncationStrategy
 end
 
diff --git a/src/algorithms/changebonds/sketchedexpand.jl b/src/algorithms/changebonds/sketchedexpand.jl
index 61501339a..b4d4680c9 100644
--- a/src/algorithms/changebonds/sketchedexpand.jl
+++ b/src/algorithms/changebonds/sketchedexpand.jl
@@ -3,27 +3,46 @@ $(TYPEDEF)
 
 An algorithm that expands the bond dimension like [`OptimalExpand`](@ref) — selecting the
 dominant directions of the projected two-site update orthogonal to the current state — but at
-single-site cost using the randomized "shrewd selection" of Controlled Bond Expansion
-(Gleis et al. Phys. Rev. Lett. 130, 246402 (2023)). A random sketch of the orthogonal complement
+single-site cost using the randomized "shrewd selection" of Controlled Bond Expansion.
+A random sketch of the orthogonal complement
 is folded into the effective environment, collapsing the large bond before the two-site update is
 ever formed, and the dominant directions are read off a small singular value decomposition.
 
 The state-preserving behaviour matches [`OptimalExpand`](@ref).
 
+!!! note
+    `trunc` bounds how much is *added* to each bond, not the total bond dimension that is kept:
+    it sizes the sketch target `Vk` within the orthogonal complement (see
+    [`sketch_space`](@ref)), so `truncrank(k)` aims to grow every bond by `k`, capped by the
+    dimension of the local complement. Because that target space is selected from uniformly
+    random weights rather than from a spectrum, its per-sector split is drawn at random rather
+    than ordered by importance — unlike [`OptimalExpand`](@ref), where the decomposition itself
+    picks out the dominant sectors — so only `truncrank` and `truncspace` have a robust meaning
+    here. The `trunc` of [`SvdCut`](@ref), and of the drivers [`DMRG`](@ref) and [`TDVP`](@ref),
+    has the other meaning: it bounds what is *kept*.
+
 !!! note
     Only defined for `FiniteMPS` (through [`changebond!`](@ref)), so it can be used standalone or
     as the `alg_expand` strategy of [`DMRG`](@ref). The reported `ϵ_2site` is a randomized
     estimate, and the folded application does not exploit `JordanMPO` sparsity.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`changebonds`](@ref) and [`changebonds!`](@ref).
+
+# References
+
+* [Gleis et al. Phys. Rev. Lett. 130, 246402 (2023)](@cite gleis2023)
 """
 @kwdef struct SketchedExpand{S} <: Algorithm
     "algorithm used to orthonormalize the sketched complement (passed as the `alg` of `left_orth!`/`right_orth!`); `nothing` selects QR without oversampling and an SVD-based decomposition otherwise"
     alg_orth::S = nothing
 
-    "algorithm used for truncating the expanded space"
+    "[truncation strategy](@extref MatrixAlgebraKit.TruncationStrategy) selecting how many directions are *added* to each bond, rather than how much of the bond is kept"
     trunc::TruncationStrategy
 
     "number of extra sketch columns drawn beyond the target rank (range-finder oversampling)"
diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl
index a1a2eee14..295e48557 100644
--- a/src/algorithms/changebonds/svdcut.jl
+++ b/src/algorithms/changebonds/svdcut.jl
@@ -6,13 +6,15 @@ This is achieved by a sweeping algorithm that locally performs (optimal) truncat
 
 changedbonds! is only defined for FiniteMPS and FiniteMPO.
 
-See also [`changebonds(!)`](@ref changebonds)
-
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`changebonds`](@ref) and [`changebonds!`](@ref).
+
+# References
 
 * [Parker et al. Phys. Rev. B 102 (2020)](@cite parker2020)
 """
diff --git a/src/algorithms/changebonds/vumpssvd.jl b/src/algorithms/changebonds/vumpssvd.jl
index f90a9d019..4064c68b7 100644
--- a/src/algorithms/changebonds/vumpssvd.jl
+++ b/src/algorithms/changebonds/vumpssvd.jl
@@ -6,9 +6,13 @@ An algorithm that uses a two-site update step to change the bond dimension of a
 !!! note
     [`changebonds!`](@ref) is not defined.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`changebonds`](@ref).
 """
 @kwdef struct VUMPSSvdCut <: Algorithm
     "algorithm used for gauging the `InfiniteMPS`"
diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl
index 438fe87ec..f965db9ac 100644
--- a/src/algorithms/derivatives/hamiltonian_derivatives.jl
+++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl
@@ -7,7 +7,7 @@ const _HAM_MPS_TYPES = Union{
 # Single site derivative
 # ----------------------
 """
-    JordanMPO_AC_Hamiltonian{O1,O2,O3}
+    JordanMPO_AC_Hamiltonian{O1, O2, O3}
 
 Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs.
 In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor.
@@ -149,7 +149,7 @@ end
 # Two site derivative
 # -------------------
 """
-    JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4}
+    JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}
 
 Efficient operator for representing the single-site derivative of a `MPOHamiltonian` sandwiched between two MPSs.
 In particular, this operator aims to make maximal use of the structure of the `MPOHamiltonian` to reduce the number of operations required to apply the operator to a tensor.
diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl
index db78fe768..5bc4bc943 100644
--- a/src/algorithms/derivatives/mpo_derivatives.jl
+++ b/src/algorithms/derivatives/mpo_derivatives.jl
@@ -1,5 +1,5 @@
 """
-    struct MPODerivativeOperator{L,O<:Tuple,R}
+    struct MPODerivativeOperator{L, O <: Tuple, R}
 
 Effective local operator obtained from taking the partial derivative of an MPS-MPO-MPS sandwich.
 """
diff --git a/src/algorithms/excitation/chepigaansatz.jl b/src/algorithms/excitation/chepigaansatz.jl
index 794032b68..8926aead0 100644
--- a/src/algorithms/excitation/chepigaansatz.jl
+++ b/src/algorithms/excitation/chepigaansatz.jl
@@ -3,11 +3,7 @@ $(TYPEDEF)
 
 Single-site optimization algorithm for excitations on top of MPS groundstates.
 
-## Fields
-
-$(TYPEDFIELDS)
-
-## Constructors
+# Constructors
 
     ChepigaAnsatz()
     ChepigaAnsatz(; kwargs...)
@@ -16,9 +12,17 @@ $(TYPEDFIELDS)
 Create a `ChepigaAnsatz` algorithm with the given eigensolver, or by passing the
 keyword arguments to [`Arnoldi`](@extref KrylovKit.Arnoldi).
 
-## References
+# Fields
+
+$(TYPEDFIELDS)
 
-- [Chepiga et al. Phys. Rev. B 96 (2017)](@cite chepiga2017) 
+# See also
+
+Used as the `algorithm` argument of [`excitations`](@ref).
+
+# References
+
+* [Chepiga et al. Phys. Rev. B 96 (2017)](@cite chepiga2017)
 """
 struct ChepigaAnsatz{A <: KrylovAlgorithm} <: Algorithm
     "algorithm used for the eigenvalue solvers"
@@ -64,29 +68,36 @@ function excitations(
 end
 
 """
-    ChepigaAnsatz2 <: Algorithm
+$(TYPEDEF)
 
 Two-site optimization algorithm for excitations on top of MPS groundstates.
 
-## Fields
-- `alg::A = Defaults.eigsolver`: algorithm to use for the eigenvalue problem.
-- `trunc = Defaults.trunc`: algorithm to use for truncation.
-
-## Constructors
+# Constructors
 
     ChepigaAnsatz2()
     ChepigaAnsatz2(; kwargs...)
     ChepigaAnsatz2(alg, trunc)
 
 Create a `ChepigaAnsatz2` algorithm with the given eigensolver and truncation, or by passing the
-keyword arguments to `Arnoldi`.
+keyword arguments to [`Arnoldi`](@extref KrylovKit.Arnoldi).
+
+# Fields
+
+$(TYPEDFIELDS)
+
+# See also
 
-## References
+Used as the `algorithm` argument of [`excitations`](@ref).
 
-- [Chepiga et al. Phys. Rev. B 96 (2017)](@cite chepiga2017) 
+# References
+
+* [Chepiga et al. Phys. Rev. B 96 (2017)](@cite chepiga2017)
 """
 struct ChepigaAnsatz2{A <: KrylovAlgorithm} <: Algorithm
+    "algorithm used for the eigenvalue solvers, defaults to `Arnoldi(; krylovdim = 30, tol = 1.0e-10, eager = true)`"
     alg::A
+
+    "[truncation strategy](@extref MatrixAlgebraKit.TruncationStrategy) used when splitting the optimized two-site tensor, defaults to `notrunc()`"
     trunc::Any
 end
 function ChepigaAnsatz2(; trunc = notrunc(), kwargs...)
diff --git a/src/algorithms/excitation/dmrgexcitation.jl b/src/algorithms/excitation/dmrgexcitation.jl
index 72c47a401..278585235 100644
--- a/src/algorithms/excitation/dmrgexcitation.jl
+++ b/src/algorithms/excitation/dmrgexcitation.jl
@@ -7,9 +7,13 @@ Variational optimization algorithm for excitations of finite MPS by minimizing t
 H + λᵢ |ψᵢ⟩⟨ψᵢ|
 ```
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`excitations`](@ref).
 """
 @kwdef struct FiniteExcited{A} <: Algorithm
     "optimization algorithm"
diff --git a/src/algorithms/excitation/excitations.jl b/src/algorithms/excitation/excitations.jl
index eb1247f90..3d0204fbd 100644
--- a/src/algorithms/excitation/excitations.jl
+++ b/src/algorithms/excitation/excitations.jl
@@ -1,18 +1,29 @@
 """
-    excitations(H, algorithm::QuasiparticleAnsatz, ψ::FiniteQP, [left_environments],
-                [right_environments]; num=1) -> (energies, states)
-    excitations(H, algorithm::QuasiparticleAnsatz, ψ::InfiniteQP, [left_environments],
-                [right_environments]; num=1, solver=Defaults.solver) -> (energies, states)
-    excitations(H, algorithm::FiniteExcited, ψs::NTuple{<:Any, <:FiniteMPS};
-                num=1, init=copy(first(ψs))) -> (energies, states)
-    excitations(H, algorithm::ChepigaAnsatz, ψ::FiniteMPS, [envs];
-                num=1, pos=length(ψ)÷2) -> (energies, states)
-    excitations(H, algorithm::ChepigaAnsatz2, ψ::FiniteMPS, [envs];
-                num=1, pos=length(ψ)÷2) -> (energies, states)
+    excitations(
+            H, algorithm::QuasiparticleAnsatz, ψ::FiniteQP, [left_environments],
+            [right_environments]; num = 1
+        ) -> (energies, states)
+    excitations(
+            H, algorithm::QuasiparticleAnsatz, ψ::InfiniteQP, [left_environments],
+            [right_environments]; num = 1
+        ) -> (energies, states)
+    excitations(
+            H, algorithm::FiniteExcited, ψs::NTuple{<:Any, <:FiniteMPS};
+            num = 1, init
+        ) -> (energies, states)
+    excitations(
+            H, algorithm::ChepigaAnsatz, ψ::FiniteMPS, [envs];
+            num = 1, pos = length(ψ) ÷ 2
+        ) -> (energies, states)
+    excitations(
+            H, algorithm::ChepigaAnsatz2, ψ::FiniteMPS, [envs];
+            num = 1, pos = length(ψ) ÷ 2
+        ) -> (energies, states)
 
 Compute the first excited states and their energy gap above a ground state.
 
 # Arguments
+
 - `H::AbstractMPO`: operator for which to find the excitations
 - `algorithm`: optimization algorithm
 - `ψ::QP`: initial quasiparticle guess
@@ -20,10 +31,11 @@ Compute the first excited states and their energy gap above a ground state.
 - `[left_environments]`: left ground state environment
 - `[right_environments]`: right ground state environment
 
-# Keywords
+# Keyword Arguments
+
 - `num::Int`: number of excited states to compute
 - `solver`: algorithm for the linear solver of the quasiparticle environments
-- `init`: initial excited state guess
+- `init`: initial excited state guess; defaults to a copy of the first state in `ψs`
 - `pos`: position of perturbation
 """
 function excitations end
diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl
index 424d97e7b..3ba679af0 100644
--- a/src/algorithms/excitation/quasiparticleexcitation.jl
+++ b/src/algorithms/excitation/quasiparticleexcitation.jl
@@ -7,22 +7,26 @@ $(TYPEDEF)
 
 Optimization algorithm for quasi-particle excitations on top of MPS groundstates.
 
-## Fields
+# Constructors
 
-$(TYPEDFIELDS)
-
-## Constructors
-    
     QuasiparticleAnsatz()
     QuasiparticleAnsatz(; kwargs...)
     QuasiparticleAnsatz(alg)
 
-Create a `QuasiparticleAnsatz` algorithm with the given algorithm, or by passing the 
-keyword arguments to `Arnoldi`.
+Create a `QuasiparticleAnsatz` algorithm with the given eigensolver, or by passing the
+keyword arguments to [`Arnoldi`](@extref KrylovKit.Arnoldi).
+
+# Fields
+
+$(TYPEDFIELDS)
+
+# See also
 
-## References
+Used as the `algorithm` argument of [`excitations`](@ref).
 
-- [Haegeman et al. Phys. Rev. Let. 111 (2013)](@cite haegeman2013)
+# References
+
+* [Haegeman et al. Phys. Rev. Let. 111 (2013)](@cite haegeman2013)
 """
 struct QuasiparticleAnsatz{A, E} <: Algorithm
     "algorithm used for the eigenvalue solvers"
@@ -66,14 +70,17 @@ function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP; num = 1, kw
 end
 
 """
-    excitations(H, algorithm::QuasiparticleAnsatz, momentum::Union{Number, Vector{<:Number}},
-                left_ψ::InfiniteMPS, [left_environment],
-                [right_ψ::InfiniteMPS], [right_environment];
-                kwargs...)
+    excitations(
+            H, algorithm::QuasiparticleAnsatz, momentum::Union{Number, Vector{<:Number}},
+            left_ψ::InfiniteMPS, [left_environment],
+            [right_ψ::InfiniteMPS], [right_environment];
+            kwargs...
+        ) -> (energies, states)
 
 Create and optimize infinite quasiparticle states.
 
 # Arguments
+
 - `H::AbstractMPO`: operator for which to find the excitations
 - `algorithm::QuasiparticleAnsatz`: optimization algorithm
 - `momentum::Union{Number, Vector{<:Number}}`: momentum or list of momenta
@@ -82,11 +89,12 @@ Create and optimize infinite quasiparticle states.
 - `[right_ψ::InfiniteMPS]`: right ground state
 - `[right_environment]`: right ground state environment
 
-# Keywords
+# Keyword Arguments
+
 - `num::Int`: number of excited states to compute
 - `solver`: algorithm for the linear solver of the quasiparticle environments
-- `sector=leftunit(left_ψ)`: charge of the quasiparticle state
-- `parallel=true`: enable multi-threading over different momenta
+- `sector = leftunit(left_ψ)`: charge of the quasiparticle state
+- `parallel = true`: enable multi-threading over different momenta
 """
 function excitations(
         H, alg::QuasiparticleAnsatz, momentum::Number, lmps::InfiniteMPS,
@@ -159,12 +167,15 @@ function excitations(
 end
 
 """
-    excitations(H, algorithm::QuasiparticleAnsatz, left_ψ::FiniteMPS, [left_environment],
-                [right_ψ::FiniteMPS], [right_environment]; kwargs...)
+    excitations(
+            H, algorithm::QuasiparticleAnsatz, left_ψ::FiniteMPS, [left_environment],
+            [right_ψ::FiniteMPS], [right_environment]; kwargs...
+        ) -> (energies, states)
 
 Create and optimize finite quasiparticle states.
 
 # Arguments
+
 - `H::AbstractMPO`: operator for which to find the excitations
 - `algorithm::QuasiparticleAnsatz`: optimization algorithm
 - `left_ψ::FiniteMPS`: left ground state
@@ -172,9 +183,10 @@ Create and optimize finite quasiparticle states.
 - `[right_ψ::FiniteMPS]`: right ground state
 - `[right_environment]`: right ground state environment
 
-# Keywords
+# Keyword Arguments
+
 - `num::Int`: number of excited states to compute
-- `sector=leftunit(lmps)`: charge of the quasiparticle state
+- `sector = leftunit(lmps)`: charge of the quasiparticle state
 """
 function excitations(
         H, alg::QuasiparticleAnsatz, lmps::FiniteMPS,
diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl
index d60213af2..b54de93d7 100644
--- a/src/algorithms/expval.jl
+++ b/src/algorithms/expval.jl
@@ -1,9 +1,9 @@
 """
-    expectation_value(ψ, O, [environments])
-    expectation_value(ψ, inds => O)
-    expectation_value(ψ, (mpo, site => O), [environments])
+    expectation_value(ψ, O, [environments]) -> val
+    expectation_value(ψ, inds => O) -> val
+    expectation_value(ψ, (mpo, site => O), [environments]) -> val
 
-Compute the expectation value of an operator `O` on a state `ψ`. 
+Compute the expectation value of an operator `O` on a state `ψ`, normalized by `⟨ψ|ψ⟩`.
 Optionally, it is possible to make the computations more efficient by also passing in
 previously calculated `environments`.
 
@@ -15,13 +15,25 @@ acts, while the operator is either a `AbstractTensorMap` or a `FiniteMPO`. In th
 the operator is a `AbstractTensorMap` that acts on the physical space of a single site.
 
 # Arguments
-* `ψ::AbstractMPS` : the state on which to compute the expectation value
-* `O::Union{AbstractMPO,Pair,AbstractTensorMap}` : the operator to compute the expectation value of. 
+
+- `ψ::AbstractMPS`: the state on which to compute the expectation value
+- `O::Union{AbstractMPO, Pair, AbstractTensorMap}`: the operator to compute the expectation value of.
     This can either be an `AbstractMPO`, a pair of indices and local operator, or a local MPO tensor
     represented as a `AbstractTensorMap`.
-* `environments::AbstractMPSEnvironments` : the environments to use for the calculation. If not given, they will be calculated.
+- `environments::AbstractMPSEnvironments`: the environments to use for the calculation. If not given, they will be calculated.
     Depending on the type of `O`, these will be the environments of the operator `O` or the MPO `mpo`.
+
+# Returns
+
+- `val::Number`: the (normalized) expectation value `⟨ψ|O|ψ⟩ / ⟨ψ|ψ⟩`.
+
+!!! note "Infinite operators"
+    For an infinite state and an infinite operator (e.g. an `InfiniteMPOHamiltonian`), the
+    return value is the total over one unit cell; divide by `length(ψ)` to obtain a
+    per-site value.
+
 # Examples
+
 ```jldoctest
 julia> ψ = FiniteMPS(ones(Float64, (ℂ^2)^4));
 
diff --git a/src/algorithms/fidelity_susceptibility.jl b/src/algorithms/fidelity_susceptibility.jl
index 63f028401..55dbac066 100644
--- a/src/algorithms/fidelity_susceptibility.jl
+++ b/src/algorithms/fidelity_susceptibility.jl
@@ -1,8 +1,10 @@
 """
-    fidelity_susceptibility(state::Union{FiniteMPS,InfiniteMPS}, H₀::T,
-                            Vs::AbstractVector{T}, [henvs=environments(state, H₀, state)];
-                            maxiter=Defaults.maxiter,
-                            tol=Defaults.tol) where {T<:MPOHamiltonian}
+    fidelity_susceptibility(
+        state::Union{FiniteMPS, InfiniteMPS}, H₀::T,
+        Vs::AbstractVector{T}, [henvs = environments(state, H₀, state)];
+        maxiter = Defaults.maxiter,
+        tol = Defaults.tol
+    ) where {T <: MPOHamiltonian}
 
 Computes the fidelity susceptibility of a the ground state `state` of a base Hamiltonian
 `H₀` with respect to a set of perturbing Hamiltonians `Vs`. Each of the perturbing
diff --git a/src/algorithms/grassmann.jl b/src/algorithms/grassmann.jl
index 8cd8b02e6..add98ffe7 100644
--- a/src/algorithms/grassmann.jl
+++ b/src/algorithms/grassmann.jl
@@ -90,7 +90,7 @@ end
 """
     retract(state, g, α) -> state′, ξ
 
-Retract a state a distance `α` along a direction `g`, obtaining a new state and the local tangent vector. 
+Retract a state a distance `α` along a direction `g`, obtaining a new state and the local tangent vector.
 """
 function retract(state::FiniteMPS, g, α::Real)
     state′ = copy(state)
@@ -137,7 +137,7 @@ function transport!(h, state, g, α::Real, state′)
 end
 
 """
-    fg(state, operator, envs=environments(state, operator, state))
+    fg(state, operator, envs = environments(state, operator, state))
 
 Compute the cost function and the tangent vector with respect to the `AL` parameters of the state.
 """
@@ -219,7 +219,7 @@ function fg(
 end
 
 """
-    rho_inv_regularized(C; rtol=eps(real(scalartype(C)))^(3/4))
+    rho_inv_regularized(C; rtol = eps(real(scalartype(C)))^(3 / 4))
 
 Compute the (regularized) inverse of the MPS fixed point `ρ = C * C'`.
 Here we use the Tikhonov regularization, i.e. `inv(ρ) = inv(C * C' + δ²1)`,
diff --git a/src/algorithms/groundstate/dmrg.jl b/src/algorithms/groundstate/dmrg.jl
index 49bcafb2f..034822208 100644
--- a/src/algorithms/groundstate/dmrg.jl
+++ b/src/algorithms/groundstate/dmrg.jl
@@ -22,9 +22,9 @@ eigensolve, and (3) a gauge step (`alg_gauge`). With the defaults (`alg_expand =
 `alg_gauge = nothing`, a non-truncating QR gauge derived from `trunc = notrunc()`) this is
 textbook single-site DMRG, which cannot change the bond dimension. Setting `alg_expand` to a
 bond-expansion algorithm (e.g. [`OptimalExpand`](@ref), [`RandExpand`](@ref), [`SketchedExpand`](@ref))
-enriches the bond with directions orthogonal to the current state ahead of each eigensolve,
+expands the bond with directions orthogonal to the current state ahead of each eigensolve,
 recovering Controlled Bond Expansion (CBE) DMRG. Setting `alg_gauge` to a bond-expanding gauge
-algorithm (e.g. [`DMRG3S`](@ref)) instead enriches the bond as part of the gauge step, after the
+algorithm (e.g. [`DMRG3S`](@ref)) instead expands the bond as part of the gauge step, after the
 eigensolve. Either way, a truncating gauge (see below) is then desirable to cut the enlarged
 bond back down.
 
@@ -51,9 +51,13 @@ If `alg_gauge` is instead given with its inner gauge already set (e.g. `DMRG3S(0
 some_gauge)`), `trunc`/`alg_svd`/`alg_orth` must be left at their defaults — passing both is an
 error, since it leaves two conflicting sources for the same setting.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref).
 """
 struct DMRG{A, F, E, G} <: Algorithm
     "tolerance for convergence criterium"
@@ -149,9 +153,13 @@ $(TYPEDEF)
 
 Two-site DMRG algorithm for finding the dominant eigenvector.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`approximate`](@ref).
 """
 struct DMRG2{A, G, F} <: Algorithm
     "tolerance for convergence criterium"
@@ -234,6 +242,26 @@ _sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)
 
 inner_alg_gauge(alg::Union{DMRG, DMRG2}) = alg_gauge(alg.alg_gauge)
 
+"""
+    find_groundstate!(ψ, H, algorithm, [environments]) -> (ψ, environments, ϵ)
+
+In-place version of [`find_groundstate`](@ref): optimize the finite MPS `ψ` for the
+Hamiltonian `H`, overwriting the input state instead of working on a copy.
+Currently supported for the finite-system algorithms [`DMRG`](@ref) and [`DMRG2`](@ref).
+
+# Arguments
+
+- `ψ::AbstractFiniteMPS`: initial guess, mutated in place
+- `H`: operator for which to find the ground state
+- `algorithm`: optimization algorithm
+- `[environments]`: MPS environment manager
+
+# Returns
+
+- `ψ::AbstractFiniteMPS`: converged ground state
+- `environments`: environments corresponding to the converged state
+- `ϵ::Float64`: final convergence error upon terminating the algorithm
+"""
 function find_groundstate!(
         ψ::AbstractFiniteMPS, H, alg::Union{DMRG, DMRG2}, envs = environments(ψ, H, ψ)
     )
diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl
index 80ed7212a..29d37271f 100644
--- a/src/algorithms/groundstate/find_groundstate.jl
+++ b/src/algorithms/groundstate/find_groundstate.jl
@@ -1,25 +1,67 @@
 """
     find_groundstate(ψ₀, H, [environments]; kwargs...) -> (ψ, environments, ϵ)
-    find_groundstate(ψ₀, H, algorithm, environments) -> (ψ, environments, ϵ)
+    find_groundstate(ψ₀, H, algorithm, [environments]) -> (ψ, environments, ϵ)
 
-Compute the ground state for Hamiltonian `H` with initial guess `ψ`. If not specified, an
-optimization algorithm will be attempted based on the supplied keywords.
+Compute the ground state for Hamiltonian `H` with initial guess `ψ₀`. If no `algorithm` is
+specified, one is selected automatically from the type of `ψ₀` and the supplied keywords
+(see the automatic-selection notes below).
+
+# Arguments
 
-## Arguments
 - `ψ₀::AbstractMPS`: initial guess
 - `H::AbstractMPO`: operator for which to find the ground state
 - `[environments]`: MPS environment manager
 - `algorithm`: optimization algorithm
 
-## Keywords
-- `tol::Float64`: tolerance for convergence criterium
-- `maxiter::Int`: maximum amount of iterations
-- `verbosity::Int`: display progress information
+# Keyword Arguments
+
+- `tol::Float64 = $(Defaults.tol)`: tolerance for the convergence criterion
+- `maxiter::Int = $(Defaults.maxiter)`: maximum number of iterations
+- `verbosity::Int = $(Defaults.verbosity)`: display progress information
+- `trunc = nothing`: if supplied, a truncation strategy that enables bond-dimension growth
+  through a two-site algorithm (see below)
+
+# Automatic algorithm selection
+
+When no `algorithm` is passed, the choice depends on the type of `ψ₀`:
+- `InfiniteMPS`: [`VUMPS`](@ref) (with its tolerance floored at `1e-4`), refined by
+  [`GradientGrassmann`](@ref) when `tol < 1e-4`. If `trunc` is given, an [`IDMRG2`](@ref)
+  stage is prepended to grow the bond dimension.
+- `AbstractFiniteMPS`: [`DMRG`](@ref). If `trunc` is given, a [`DMRG2`](@ref) stage is
+  prepended to grow the bond dimension.
+
+Because single-site [`DMRG`](@ref) preserves the bond dimension of `ψ₀`, passing a
+`trunc` (or an explicit two-site `algorithm`) is the usual way to converge from a
+low-bond-dimension initial guess such as a product state.
+
+# Returns
 
-## Returns
 - `ψ::AbstractMPS`: converged ground state
 - `environments`: environments corresponding to the converged state
 - `ϵ::Float64`: final convergence error upon terminating the algorithm
+
+# Examples
+
+Ground state of a 4-site transverse-field Ising chain, `H = -∑ XₖXₖ₊₁ - ∑ Zₖ`, starting
+from a product state and letting `DMRG2` grow the bond dimension:
+
+```jldoctest
+julia> X = TensorMap(Float64[0 1; 1 0], ℂ^2, ℂ^2);
+
+julia> Z = TensorMap(Float64[1 0; 0 -1], ℂ^2, ℂ^2);
+
+julia> L = 4; lattice = fill(ℂ^2, L);
+
+julia> H = FiniteMPOHamiltonian(lattice, ((i, i + 1) => -(X ⊗ X) for i in 1:(L - 1))) +
+           FiniteMPOHamiltonian(lattice, ((i,) => -Z for i in 1:L));
+
+julia> ψ₀ = FiniteMPS(ones(Float64, (ℂ^2)^L));
+
+julia> ψ, envs, ϵ = find_groundstate(ψ₀, H; verbosity = 0, trunc = truncrank(16));
+
+julia> round(real(expectation_value(ψ, H)); digits = 4)
+-4.7588
+```
 """
 function find_groundstate(
         ψ::AbstractMPS, H, envs::AbstractMPSEnvironments = environments(ψ, H, ψ);
diff --git a/src/algorithms/groundstate/gradient_grassmann.jl b/src/algorithms/groundstate/gradient_grassmann.jl
index 049cfaae5..c5be94d82 100644
--- a/src/algorithms/groundstate/gradient_grassmann.jl
+++ b/src/algorithms/groundstate/gradient_grassmann.jl
@@ -2,29 +2,33 @@
 $(TYPEDEF)
 
 Variational gradient-based optimization algorithm that keeps the MPS in left-canonical form,
-as points on a Grassmann manifold. The optimization is then a Riemannian gradient descent 
+as points on a Grassmann manifold. The optimization is then a Riemannian gradient descent
 with a preconditioner to induce the metric from the Hilbert space inner product.
 
-## Fields
+# Constructors
 
-$(TYPEDFIELDS)
+    GradientGrassmann(; kwargs...)
 
-## References
+# Keyword Arguments
 
-* [Hauru et al. SciPost Phys. 10 (2021)](@cite hauru2021)
+- `method = ConjugateGradient`: instance of optimization algorithm, or type of optimization
+    algorithm to construct
+- `finalize!`: finalizer algorithm
+- `tol = Defaults.tol`: tolerance for convergence criterium
+- `maxiter = Defaults.maxiter`: maximum amount of iterations
+- `verbosity = Defaults.verbosity - 1`: level of information display
 
----
+# Fields
 
-## Constructors
-    GradientGrassmann(; kwargs...)
+$(TYPEDFIELDS)
 
-### Keywords
-- `method=ConjugateGradient`: instance of optimization algorithm, or type of optimization
-    algorithm to construct
-- `finalize!`: finalizer algorithm
-- `tol::Float64`: tolerance for convergence criterium
-- `maxiter::Int`: maximum amount of iterations
-- `verbosity::Int`: level of information display
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`leading_boundary`](@ref).
+
+# References
+
+* [Hauru et al. SciPost Phys. 10 (2021)](@cite hauru2021)
 """
 struct GradientGrassmann{O <: OptimKit.OptimizationAlgorithm, F} <: Algorithm
     "optimization algorithm"
diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl
index c12e07a7a..a7f8d326a 100644
--- a/src/algorithms/groundstate/idmrg.jl
+++ b/src/algorithms/groundstate/idmrg.jl
@@ -3,9 +3,13 @@ $(TYPEDEF)
 
 Single site infinite DMRG algorithm for finding the dominant eigenvector.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref).
 """
 @kwdef struct IDMRG{A} <: Algorithm
     "tolerance for convergence criterium"
@@ -29,9 +33,13 @@ $(TYPEDEF)
 
 Two-site infinite DMRG algorithm for finding the dominant eigenvector.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref), [`leading_boundary`](@ref), and [`approximate`](@ref).
 """
 @kwdef struct IDMRG2{A, S} <: Algorithm
     "tolerance for convergence criterium"
diff --git a/src/algorithms/groundstate/vumps.jl b/src/algorithms/groundstate/vumps.jl
index 1f1cf51f7..149194c06 100644
--- a/src/algorithms/groundstate/vumps.jl
+++ b/src/algorithms/groundstate/vumps.jl
@@ -3,11 +3,15 @@ $(TYPEDEF)
 
 Variational optimization algorithm for uniform matrix product states, based on the combination of DMRG with matrix product state tangent space concepts.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`leading_boundary`](@ref).
+
+# References
 
 * [Zauner-Stauber et al. Phys. Rev. B 97 (2018)](@cite zauner-stauber2018)
 * [Vanderstraeten et al. SciPost Phys. Lect. Notes 7 (2019)](@cite vanderstraeten2019)
diff --git a/src/algorithms/post_expand/dmrg3s.jl b/src/algorithms/post_expand/dmrg3s.jl
index 1b8cf419a..065c8f16e 100644
--- a/src/algorithms/post_expand/dmrg3s.jl
+++ b/src/algorithms/post_expand/dmrg3s.jl
@@ -36,7 +36,7 @@ Base.:∘(s1::NoiseSchedule, s2::NoiseSchedule) =
 
 Noise schedule that shrinks geometrically: `noise -> noise * decay_rate^iter`, snapped
 to exactly zero once it falls below `threshold`. Use `decay_rate < 1` for a standalone
-[`DMRG3S`](@ref) run that gradually turns enrichment off as the state converges; a
+[`DMRG3S`](@ref) run that gradually turns the expansion off as the state converges; a
 nonzero `threshold` avoids running the (cheap, but non-free) expansion step
 indefinitely on a noise amplitude too small to matter.
 """
@@ -67,7 +67,7 @@ end
 (s::Warmup)(noise, iter, ϵ) = iter ≤ s.iters ? noise : zero(noise)
 
 """
-    DMRG3S(noise, schedule::NoiseSchedule)
+$(TYPEDEF)
 
 Gauge algorithm wrapper that, at every site update, injects a Hamiltonian-derived
 perturbation of the just-optimized tensor before gauging — the "strictly single-site DMRG with
@@ -75,10 +75,14 @@ subspace expansion" scheme. This lets
 single-site DMRG introduce basis states/quantum-number sectors absent from the initial
 state, helping it escape local minima that plain single-site DMRG can get stuck in.
 
+# Constructors
+
+    DMRG3S(noise, schedule::NoiseSchedule)
+
 `noise` is the initial perturbation amplitude; `schedule` (see [`ExponentialDecay`](@ref),
 [`Warmup`](@ref)) controls how it evolves across outer iterations, and once it decays to
 exactly zero the gauge step reverts to a plain gauge shift for the remainder of the
-run. The actual factorization used to gauge the enriched tensor
+run. The actual factorization used to gauge the expanded tensor
 is filled in by `DMRG`'s constructor, not supplied here directly — see `DMRG`'s docstring
 for the calling convention:
 
@@ -89,13 +93,24 @@ DMRG(; alg_gauge = DMRG3S(0.1, ExponentialDecay(0.7)), trunc = truncdim(50))
 A truncating `trunc` is strongly recommended alongside `DMRG3S`, to cut the perturbed
 bond back down each sweep — `DMRG`'s constructor warns if none is given.
 
-## References
+# Fields
+
+$(TYPEDFIELDS)
+
+# See also
+
+Used as the `alg_gauge` argument of [`DMRG`](@ref).
+
+# References
 
-* [Hubig et al. Phys. Rev. B 91, 155115 (2015)](@cite Hubig2015)
+* [Hubig et al. Phys. Rev. B 91, 155115 (2015)](@cite hubig2015)
 """
 struct DMRG3S{N, S <: NoiseSchedule, A} <: Algorithm
+    "initial perturbation amplitude, before `schedule` is applied"
     noise::N
+    "[`NoiseSchedule`](@ref) controlling how the amplitude evolves across outer iterations"
     schedule::S
+    "factorization used to gauge the expanded tensor; `nothing` until [`DMRG`](@ref)'s constructor fills it in"
     alg_gauge::A
 end
 
diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl
index 93bfd6637..89466e7e0 100644
--- a/src/algorithms/propagator/corvector.jl
+++ b/src/algorithms/propagator/corvector.jl
@@ -11,11 +11,15 @@ $(TYPEDEF)
 A dynamical DMRG method for calculating dynamical properties and excited states, based on a
 variational principle for dynamical correlation functions.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`propagator`](@ref).
+
+# References
 
 * [Jeckelmann. Phys. Rev. B 66 (2002)](@cite jeckelmann2002)
 """
@@ -33,13 +37,15 @@ $(TYPEDFIELDS)
 end
 
 """
-    propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init=copy(ψ₀))
+    propagator(ψ₀::AbstractFiniteMPS, z::Number, H::MPOHamiltonian, alg::DynamicalDMRG; init = copy(ψ₀)) -> (g, ψ)
 
 Calculate the action of the propagator ``\\frac{1}{z - H}|ψ₀⟩`` using the dynamical DMRG
 algorithm.
 
-Returns a tuple `(g, ψ)` where `g` is the approximation of the propagator matrix element
-``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and `ψ` is the MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩``.
+# Returns
+
+- `g`: approximation of the propagator matrix element ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩``
+- `ψ`: MPS approximation of ``\\frac{1}{z - H}|ψ₀⟩``
 """
 function propagator end
 
@@ -55,7 +61,9 @@ This algorithm minimizes the following cost function
 
 Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``.
 
-See also [`Jeckelmann`](@ref) for the original approach.
+# See also
+
+[`Jeckelmann`](@ref) for the original approach.
 """
 struct NaiveInvert <: DDMRG_Flavour end
 
@@ -121,9 +129,11 @@ Together with equation (11) from that same paper we can determine the full propa
 
 Returns the approximation of ``⟨ψ₀|\\frac{1}{z - H}|ψ₀⟩`` and ``\\frac{1}{z - H}|ψ₀⟩``.
 
-See also [`NaiveInvert`](@ref) for a less costly but less accurate alternative.
+# See also
+
+[`NaiveInvert`](@ref) for a less costly but less accurate alternative.
 
-## References
+# References
 
 * [Jeckelmann. Phys. Rev. B 66 (2002)](@cite jeckelmann2002)
 """
diff --git a/src/algorithms/statmech/leading_boundary.jl b/src/algorithms/statmech/leading_boundary.jl
index d66088eea..218a3f8f7 100644
--- a/src/algorithms/statmech/leading_boundary.jl
+++ b/src/algorithms/statmech/leading_boundary.jl
@@ -5,18 +5,21 @@
 Compute the leading boundary MPS for operator `O` with initial guess `ψ`. If not specified, an
 optimization algorithm will be attempted based on the supplied keywords.
 
-## Arguments
+# Arguments
+
 - `ψ₀::AbstractMPS`: initial guess
 - `O::AbstractMPO`: operator for which to find the leading_boundary
 - `[environments]`: MPS environment manager
 - `algorithm`: optimization algorithm
 
-## Keywords
+# Keyword Arguments
+
 - `tol::Float64`: tolerance for convergence criterium
 - `maxiter::Int`: maximum amount of iterations
 - `verbosity::Int`: display progress information
 
-## Returns
+# Returns
+
 - `ψ::AbstractMPS`: converged leading boundary MPS
 - `environments`: environments corresponding to the converged boundary
 - `ϵ::Float64`: final convergence error upon terminating the algorithm
diff --git a/src/algorithms/statmech/vomps.jl b/src/algorithms/statmech/vomps.jl
index 935d9111f..faa2406cb 100644
--- a/src/algorithms/statmech/vomps.jl
+++ b/src/algorithms/statmech/vomps.jl
@@ -1,15 +1,19 @@
 """
 $(TYPEDEF)
-    
+
 Power method algorithm for finding dominant eigenvectors of infinite MPOs.
 This method works by iteratively approximating the product of an operator and a state
 with a new state of the same bond dimension.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`leading_boundary`](@ref) and [`approximate`](@ref).
+
+# References
 
 * [Vanhecke et al. SciPost Phys. Core 4 (2021)](@cite vanhecke2021)
 """
diff --git a/src/algorithms/timestep/integrators.jl b/src/algorithms/timestep/integrators.jl
index 7b2983c34..7b75554bf 100644
--- a/src/algorithms/timestep/integrators.jl
+++ b/src/algorithms/timestep/integrators.jl
@@ -1,10 +1,11 @@
 """
-    integrate(f, y₀, t, dt, alg)
+    integrate(f, y₀, t, dt, alg) -> y
 
-Integrate the differential equation ``i dy/dt = f(y, t)`` over a time step 'dt' starting from
+Integrate the differential equation ``i dy/dt = f(y, t)`` over a time step `dt` starting from
 ``y(t₀)=y₀``, using the provided algorithm.
 
 # Arguments
+
 - `f`: driving function
 - `y₀`: object to integrate
 - `t::Number`: starting time of time-step
diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl
index 9e5943288..b129e9d0c 100644
--- a/src/algorithms/timestep/taylorcluster.jl
+++ b/src/algorithms/timestep/taylorcluster.jl
@@ -3,11 +3,15 @@ $(TYPEDEF)
 
 Algorithm for constructing the `N`th order time evolution MPO using the Taylor cluster expansion.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`make_time_mpo`](@ref).
+
+# References
 
 * [Van Damme et al. SciPost Phys. 17 (2024)](@cite vandamme2024)
 """
@@ -21,7 +25,7 @@ $(TYPEDFIELDS)
 end
 
 """
-    const WI = TaylorCluster(; N=1, extension=false, compression=false)
+    const WI = TaylorCluster(; N = 1, extension = false, compression = false)
 
 First order Taylor expansion for a time-evolution MPO.
 """
diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl
index c14d05d57..06f2ef603 100644
--- a/src/algorithms/timestep/tdvp.jl
+++ b/src/algorithms/timestep/tdvp.jl
@@ -4,7 +4,7 @@ $(TYPEDEF)
 Single site MPS time-evolution algorithm based on the Time-Dependent Variational Principle.
 
 For finite MPS, setting `alg_expand` to a bond-expansion algorithm (e.g. [`OptimalExpand`](@ref),
-[`SketchedExpand`](@ref)) enriches the bond with directions orthogonal to the current state
+[`SketchedExpand`](@ref)) expands the bond with directions orthogonal to the current state
 ahead of each local integration, recovering Controlled Bond Expansion (CBE) TDVP and lifting the
 fixed-bond limitation of plain single-site TDVP. A truncating `trunc` is then required to cut
 the enlarged bond back down (selecting the truncated-SVD gauge). The expansion is
@@ -16,11 +16,15 @@ state-preserving, as required for a consistent time evolution.
     evolution instead renormalizes at every step, like a ground-state search. CBE is only
     available for finite MPS.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`timestep`](@ref), [`timestep!`](@ref) and [`time_evolve`](@ref).
+
+# References
 
 * [Haegeman et al. Phys. Rev. Lett. 107 (2011)](@cite haegeman2011)
 """
@@ -190,11 +194,15 @@ $(TYPEDEF)
 
 Two-site MPS time-evolution algorithm based on the Time-Dependent Variational Principle.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`timestep`](@ref), [`timestep!`](@ref) and [`time_evolve`](@ref).
+
+# References
 
 * [Haegeman et al. Phys. Rev. Lett. 107 (2011)](@cite haegeman2011)
 """
diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl
index 51288abd7..5410011d3 100644
--- a/src/algorithms/timestep/time_evolve.jl
+++ b/src/algorithms/timestep/time_evolve.jl
@@ -1,25 +1,30 @@
 """
-    time_evolve(ψ₀, H, t_span, [alg], [envs]; kwargs...) -> (ψ, envs)
-    time_evolve!(ψ₀, H, t_span, [alg], [envs]; kwargs...) -> (ψ₀, envs)
+    time_evolve(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ, envs)
+    time_evolve!(ψ₀, H, t_span, alg, [envs]; kwargs...) -> (ψ₀, envs)
 
 Time-evolve the initial state `ψ₀` with Hamiltonian `H` over a given time span by stepping
 through each of the time points obtained by iterating t_span.
 
-## Arguments
+# Arguments
 
 - `ψ₀::AbstractMPS`: initial state
 - `H::AbstractMPO`: operator that generates the time evolution (can be time-dependent).
 - `t_span::AbstractVector{<:Number}`: time points over which the time evolution is stepped
-- `[alg]`: algorithm to use for the time evolution. Defaults to [`TDVP`](@ref).
-- `[envs]`: MPS environment manager
+- `alg`: algorithm to use for the time evolution, e.g. [`TDVP`](@ref) or [`TDVP2`](@ref).
+- `envs`: MPS environment manager
 
-## Keyword Arguments
+# Keyword Arguments
 
-- `verbosity::Int=0`: verbosity level for logging
-- `imaginary_evolution::Bool=false`: if true, the time evolution is done with an imaginary time step
+- `verbosity::Int = 0`: verbosity level for logging
+- `imaginary_evolution::Bool = false`: if true, the time evolution is done with an imaginary time step
     instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this
     function to compute the ground state of a Hamiltonian, or to compute finite-temperature
     properties of a system.
+
+# Returns
+
+- `ψ`: the time-evolved state
+- `envs`: the updated environment manager
 """
 function time_evolve end, function time_evolve! end
 
@@ -48,27 +53,54 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time
 end
 
 """
-    timestep(ψ₀, H, t, dt, [alg], [envs]; kwargs...) -> (ψ, envs)
-    timestep!(ψ₀, H, t, dt, [alg], [envs]; kwargs...) -> (ψ₀, envs)
+    timestep(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ, envs)
+    timestep!(ψ₀, H, t, dt, alg, [envs]; kwargs...) -> (ψ₀, envs)
 
 Time-step the state `ψ₀` with Hamiltonian `H` over a given time step `dt` at time `t`,
 solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``.
 
-## Arguments
+# Arguments
 
 - `ψ₀::AbstractMPS`: initial state
 - `H::AbstractMPO`: operator that generates the time evolution (can be time-dependent).
 - `t::Number`: starting time of time-step
 - `dt::Number`: time-step magnitude
-- `[alg]`: algorithm to use for the time evolution. Defaults to [`TDVP`](@ref).
-- `[envs]`: MPS environment manager
+- `alg`: algorithm to use for the time evolution, e.g. [`TDVP`](@ref) or [`TDVP2`](@ref).
+- `envs`: MPS environment manager
 
-## Keyword Arguments
+# Keyword Arguments
 
-- `imaginary_evolution::Bool=false`: if true, the time evolution is done with an imaginary time step
+- `imaginary_evolution::Bool = false`: if true, the time evolution is done with an imaginary time step
     instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this
     function to compute the ground state of a Hamiltonian, or to compute finite-temperature
     properties of a system.
+
+# Returns
+
+- `ψ`: the time-stepped state
+- `envs`: the updated environment manager
+
+# Examples
+
+Real-time evolution of the `|+···+⟩` product state under a transverse field `H = ∑ Zₖ`.
+Each spin precesses independently, so `⟨Xₖ(t)⟩ = cos(2t)`; after a step `dt = 0.1` this is
+`cos(0.2) ≈ 0.980067`. The initial state must be complex, since real-time evolution
+multiplies by `-i`:
+
+```jldoctest
+julia> X = TensorMap(ComplexF64[0 1; 1 0], ℂ^2, ℂ^2);
+
+julia> Z = TensorMap(ComplexF64[1 0; 0 -1], ℂ^2, ℂ^2);
+
+julia> ψ₀ = FiniteMPS(ones(ComplexF64, (ℂ^2)^4));
+
+julia> H = FiniteMPOHamiltonian(fill(ℂ^2, 4), ((i,) => Z for i in 1:4));
+
+julia> ψ, envs = timestep(ψ₀, H, 0.0, 0.1, TDVP());
+
+julia> round(real(expectation_value(ψ, 2 => X)); digits = 6)
+0.980067
+```
 """
 function timestep end, function timestep! end
 
@@ -77,9 +109,9 @@ function timestep end, function timestep! end
 
 Construct an `MPO` that approximates ``\\exp(-iHdt)``.
 
-## Keyword Arguments
+# Keyword Arguments
 
-- `imaginary_evolution::Bool=false`: if true, the time evolution operator is constructed
+- `imaginary_evolution::Bool = false`: if true, the time evolution operator is constructed
     with an imaginary time step instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``).
     This can be useful for using this function to compute the ground state of a Hamiltonian,
     or to compute finite-temperature properties of a system.
diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl
index 5c39500fd..e41405fd2 100644
--- a/src/algorithms/timestep/wii.jl
+++ b/src/algorithms/timestep/wii.jl
@@ -3,11 +3,15 @@ $(TYPEDEF)
 
 Generalization of the Euler approximation of the operator exponential for MPOs.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-## References
+# See also
+
+Used as the `algorithm` argument of [`make_time_mpo`](@ref).
+
+# References
 
 * [Zaletel et al. Phys. Rev. B 91 (2015)](@cite zaletel2015)
 * [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019)
diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl
index d1a9e07c5..de11eabcd 100644
--- a/src/algorithms/toolbox.jl
+++ b/src/algorithms/toolbox.jl
@@ -88,7 +88,7 @@ function entanglement_spectrum(st::FiniteMPS, site::Int)
 end
 
 """
-    variance(state, hamiltonian, [envs=environments(state, hamiltonian, state)])
+    variance(state, hamiltonian, [envs = environments(state, hamiltonian, state)])
 
 Compute the variance of the energy of the state with respect to the Hamiltonian.
 """
diff --git a/src/algorithms/transfer_spectrum.jl b/src/algorithms/transfer_spectrum.jl
index 4581ff01f..9b498ba8b 100644
--- a/src/algorithms/transfer_spectrum.jl
+++ b/src/algorithms/transfer_spectrum.jl
@@ -8,12 +8,12 @@ The result is returned as a `TensorKit.SectorVector`, whose values can be inspec
 # 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`.
+- `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
+# 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.
@@ -115,7 +115,7 @@ function approx_angles(spectrum; tol_angle = 0.1)
 end
 
 """
-    marek_gap(above::InfiniteMPS; sector=nothing, kwargs...)
+    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 `θ`.
@@ -158,7 +158,7 @@ function marek_gap(spectrum::AbstractVector{T}; tol_angle = 0.1) where {T <: Num
 end
 
 """
-    correlation_length(above::InfiniteMPS; sector=nothing, kwargs...)
+    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.
diff --git a/src/algorithms/unionalg.jl b/src/algorithms/unionalg.jl
index 9ef7f722a..801ffa606 100644
--- a/src/algorithms/unionalg.jl
+++ b/src/algorithms/unionalg.jl
@@ -1,11 +1,16 @@
 """
 $(TYPEDEF)
 
-Algorithm wrapper representing the sequential application of two algorithms.
+Algorithm wrapper representing the sequential application of two algorithms, as produced by
+`alg1 & alg2`.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `algorithm` argument of [`find_groundstate`](@ref) and [`changebonds`](@ref).
 """
 struct UnionAlg{A, B} <: Algorithm
     "first algorithm"
diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl
index cb11cd53f..6605ccc74 100644
--- a/src/environments/finite_envs.jl
+++ b/src/environments/finite_envs.jl
@@ -97,6 +97,16 @@ function poison!(ca::FiniteEnvironments, ind)
 end
 
 #rightenv[ind] will be contracteable with the tensor on site [ind]
+"""
+    rightenv(envs, site, state)
+
+Return the right environment stored in `envs` at `site` for the given `state`: the contraction
+of everything to the right of `site` in the network the environments were built for.
+The result is gauge-compatible with the tensor of `state` at `site` and can be contracted onto
+it directly.
+
+See also [`leftenv`](@ref) and [`environments`](@ref).
+"""
 function rightenv(ca::FiniteEnvironments, ind, state)
     a = findfirst(i -> !(state.AR[i] === ca.rdependencies[i]), length(state):-1:(ind + 1))
     a = isnothing(a) ? nothing : length(state) - a + 1
@@ -114,6 +124,16 @@ function rightenv(ca::FiniteEnvironments, ind, state)
     return ca.GRs[ind + 1]
 end
 
+"""
+    leftenv(envs, site, state)
+
+Return the left environment stored in `envs` at `site` for the given `state`: the contraction
+of everything to the left of `site` in the network the environments were built for.
+The result is gauge-compatible with the tensor of `state` at `site` and can be contracted onto
+it directly.
+
+See also [`rightenv`](@ref) and [`environments`](@ref).
+"""
 function leftenv(ca::FiniteEnvironments, ind, state)
     a = findfirst(i -> !(state.AL[i] === ca.ldependencies[i]), 1:(ind - 1))
 
diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl
index d9a763402..d7133f8b6 100644
--- a/src/operators/jordanmpotensor.jl
+++ b/src/operators/jordanmpotensor.jl
@@ -1,5 +1,5 @@
 """
-    JordanMPOTensor{T,S,A} <: AbstractBlockTensorMap{T,S,2,2}
+$(TYPEDEF)
 
 A single tensor of a matrix product operator (MPO) in upper triangular (Jordan) block form,
 as used to represent the local tensors of an [`MPOHamiltonian`](@ref).
@@ -16,28 +16,30 @@ The virtual (row, column) structure is
 where `A` is the bulk of interacting operators, `C`/`B` are the operators that start/finish
 an interaction, `D` is the on-site term, and the diagonal `1`s are identities.
 
-## Representation
+# Type parameters
 
-Rather than storing the dense block matrix, the genuine operators and the identities are kept separately:
+- `T <: Number`: the `scalartype` of the tensors.
+- `S`: the `spacetype` of the tensors.
+- `A <: DenseVector{T}`: the storage type of the underlying tensors.
 
-- `tensors::SparseBlockTensorMap` holds the non-identity operators over the *full* virtual
-  space (so `A`, `B`, `C` and `D` all live at their `(row, 1, 1, col)` position).
-- `scalars::Dict{CartesianIndex{4},T}` holds the scalar multiples of the identity, keyed by
-  their `(row, 1, 1, col)` virtual position; the diagonal corner `1`s are stored here as well.
+# Properties
 
-An index is never present in both `tensors` and `scalars`.
-This keeps the ubiquitous identity blocks free of dense storage and lets identities be materialized lazily only when needed.
+The reduced-leg `A`, `B`, `C` and `D` blocks are exposed as properties (`W.A`, `W.B`, `W.C`,
+`W.D`), reconstructed on demand from the stored `tensors` and `scalars`.
 
-## Type parameters
+# Notes
 
-- `T <: Number`: the `scalartype` of the tensors.
-- `S`: the `spacetype` of the tensors.
-- `A <: DenseVector{T}`: the storage type of the underlying tensors.
+Rather than storing the dense block matrix, the genuine operators and the identities are kept
+separately:
 
-## Block accessors
+- `tensors::SparseBlockTensorMap` holds the non-identity operators over the *full* virtual
+  space (so `A`, `B`, `C` and `D` all live at their `(row, 1, 1, col)` position).
+- `scalars::Dict{CartesianIndex{4}, T}` holds the scalar multiples of the identity, keyed by
+  their `(row, 1, 1, col)` virtual position; the diagonal corner `1`s are stored here as well.
 
-The reduced-leg `A`, `B`, `C` and `D` blocks are exposed as properties (`W.A`, `W.B`, `W.C`, `W.D`),
-reconstructed on demand from `tensors` and `scalars`.
+An index is never present in both `tensors` and `scalars`.
+This keeps the ubiquitous identity blocks free of dense storage and lets identities be
+materialized lazily only when needed.
 """
 struct JordanMPOTensor{
         T <: Number, S, A <: DenseVector{T},
diff --git a/src/operators/lazysum.jl b/src/operators/lazysum.jl
index 49af40ba4..5cdd903ba 100644
--- a/src/operators/lazysum.jl
+++ b/src/operators/lazysum.jl
@@ -1,19 +1,21 @@
 """
-    LazySum{O} <: AbstractVector{O}
+$(TYPEDEF)
 
-Type that represents a lazy sum i.e explicit summation is only done when needed.
-This type is basically an `AbstractVector` with some extra functionality to calculate things efficiently.
+Type that represents a lazy sum, i.e. explicit summation is only done when needed.
+This type is basically an `AbstractVector` with some extra functionality to calculate things
+efficiently.
 
-## Fields
-- ops -- Vector of summable objects
+# Constructors
 
----
-
-## Constructors
     LazySum(x::Vector)
+    LazySum(ops::AbstractVector, fs::AbstractVector)
+
+# Fields
 
+$(TYPEDFIELDS)
 """
 struct LazySum{O} <: AbstractVector{O}
+    "vector of summable objects"
     ops::Vector{O}
 end
 
diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl
index c5c033fe4..cec0de8f6 100644
--- a/src/operators/mpo.jl
+++ b/src/operators/mpo.jl
@@ -1,5 +1,5 @@
 """
-    struct MPO{O,V<:AbstractVector{O}} <: AbstractMPO{O}
+    struct MPO{TO, V <: AbstractVector{TO}} <: AbstractMPO{TO}
 
 Matrix Product Operator (MPO) acting on a tensor product space with a linear order.
 
@@ -11,7 +11,7 @@ end
 
 """
     FiniteMPO(Os::Vector{O}) -> FiniteMPO{O}
-    FiniteMPO(O::AbstractTensorMap{S,N,N}) where {S,N} -> FiniteMPO{O<:MPOTensor}
+    FiniteMPO(O::AbstractTensorMap{S, N, N}) where {S, N} -> FiniteMPO{O <: MPOTensor}
 
 Matrix Product Operator (MPO) acting on a finite tensor product space with a linear order.
 """
@@ -447,8 +447,8 @@ function Base.isapprox(
 end
 
 @doc """
-    swap(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=Defaults.alg_svd(), trunc)
-    swap!(mpo::FiniteMPO, i::Integer; inv::Bool=false, alg=Defaults.alg_svd(), trunc)
+    swap(mpo::FiniteMPO, i::Integer; inv::Bool = false, alg = Defaults.alg_svd(), trunc)
+    swap!(mpo::FiniteMPO, i::Integer; inv::Bool = false, alg = Defaults.alg_svd(), trunc)
 
 Compose the mpo with a swap gate applied to indices `i` and `i + 1`, effectively creating an
 operator that acts on the Hilbert spaces with those factors swapped.
diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl
index cc7e42b15..1ca4bacb6 100644
--- a/src/operators/mpohamiltonian.jl
+++ b/src/operators/mpohamiltonian.jl
@@ -1,10 +1,9 @@
 """
-    MPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, local_operators...)
-    MPOHamiltonian(lattice::AbstractArray{<:VectorSpace})
-    MPOHamiltonian(x::AbstractArray{<:Any,3})
+$(TYPEDEF)
 
-MPO representation of a Hamiltonian. This is a specific form of an [`AbstractMPO`](@ref), where
-all the sites are represented by an upper triangular block matrix of the following form:
+MPO representation of a Hamiltonian.
+This is a specific form of an [`AbstractMPO`](@ref), where all the sites are represented by an
+upper triangular block matrix of the following form:
 
 ```math
 \\begin{pmatrix}
@@ -16,17 +15,46 @@ all the sites are represented by an upper triangular block matrix of the followi
 
 where `A`, `B`, `C`, and `D` are `MPOTensor`s, or (sparse) blocks thereof.
 
-## Examples
+# Constructors
 
-For example, constructing a nearest-neighbour Hamiltonian would look like this:
+The finite and infinite variants, [`FiniteMPOHamiltonian`](@ref) and
+[`InfiniteMPOHamiltonian`](@ref), are constructed from a lattice of physical spaces together
+with a set of `inds => operator` pairs describing the local terms:
 
-```julia
-lattice = fill(ℂ^2, 10)
-H = MPOHamiltonian(lattice, (i, i+1) => O for i in 1:length(lattice)-1)
+    FiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, local_operators...)
+    InfiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, local_operators...)
+
+# Properties
+
+- `A`: bulk block of interacting operators at each site
+- `B`: operators that finish an interaction
+- `C`: operators that start an interaction
+- `D`: on-site terms
+
+# Examples
+
+A nearest-neighbour term is a two-element index tuple `(i, i + 1) => O₁₂`; an on-site term
+is a one-element tuple `(i,) => O`. For the finite variant the lattice lists every site; for
+the infinite variant it is a single unit cell and indices wrap around it periodically.
+
+```jldoctest
+julia> X = TensorMap(Float64[0 1; 1 0], ℂ^2, ℂ^2);
+
+julia> Hf = FiniteMPOHamiltonian(fill(ℂ^2, 3), ((i, i + 1) => X ⊗ X for i in 1:2));
+
+julia> Hf isa FiniteMPOHamiltonian, length(Hf)
+(true, 3)
+
+julia> Hi = InfiniteMPOHamiltonian(fill(ℂ^2, 1), (1, 2) => X ⊗ X, (1,) => X);
+
+julia> Hi isa InfiniteMPOHamiltonian, length(Hi)
+(true, 1)
 ```
 
-See also [`instantiate_operator`](@ref), which is responsible for instantiating the local
-operators in a form that is compatible with this constructor.
+# See also
+
+[`instantiate_operator`](@ref) is responsible for instantiating the local operators in a form
+that is compatible with this constructor.
 """
 struct MPOHamiltonian{TO <: JordanMPOTensor, V <: AbstractVector{TO}} <: AbstractMPO{TO}
     W::V
@@ -62,8 +90,8 @@ end
     FiniteMPOHamiltonian(Ws::Vector{<:AbstractMatrix})
 
 Create a `FiniteMPOHamiltonian` from a vector of matrices, such that `Ws[i][j, k]` represents
-the operator at site `i`, left level `j` and right level `k`. Here, the entries can be
-either `MPOTensor`, `Missing` or `Number`.
+the operator at site `i`, left level `j` and right level `k`.
+Here, the entries can be either `MPOTensor`, `Missing` or `Number`.
 """
 function FiniteMPOHamiltonian(Ws::Vector{<:AbstractMatrix})
     T = promote_type(_split_mpoham_types.(Ws)...)
@@ -151,11 +179,11 @@ function FiniteMPOHamiltonian{O}(W_mats::Vector{<:AbstractMatrix}) where {O <: J
 end
 
 """
-    InfiniteMPOHamiltonian(Ws::Vector{<:Matrix})
+    InfiniteMPOHamiltonian(Ws::Vector{<:AbstractMatrix})
 
-Create a `InfiniteMPOHamiltonian` from a vector of matrices, such that `Ws[i][j, k]` represents
-the the operator at site `i`, left level `j` and right level `k`. Here, the entries can be
-either `MPOTensor`, `Missing` or `Number`.
+Create an `InfiniteMPOHamiltonian` from a vector of matrices, such that `Ws[i][j, k]`
+represents the operator at site `i`, left level `j` and right level `k`.
+Here, the entries can be either `MPOTensor`, `Missing` or `Number`.
 """
 function InfiniteMPOHamiltonian(Ws::Vector{<:AbstractMatrix})
     T = promote_type(_split_mpoham_types.(Ws)...)
diff --git a/src/operators/multilinempo.jl b/src/operators/multilinempo.jl
index d49413e63..7432cf91a 100644
--- a/src/operators/multilinempo.jl
+++ b/src/operators/multilinempo.jl
@@ -6,10 +6,13 @@
 Type that represents multiple lines of `MPO` objects.
 
 # Constructors
-    MultilineMPO(mpos::AbstractVector{<:Union{SparseMPO,DenseMPO}})
+
+    MultilineMPO(mpos::AbstractVector{<:Union{SparseMPO, DenseMPO}})
     MultilineMPO(Os::AbstractMatrix{<:MPOTensor})
 
-See also: [`Multiline`](@ref), [`AbstractMPO`](@ref)
+# See also
+
+[`Multiline`](@ref), [`AbstractMPO`](@ref)
 """
 const MultilineMPO = Multiline{<:AbstractMPO}
 
diff --git a/src/operators/multipliedoperator.jl b/src/operators/multipliedoperator.jl
index 1cb38d912..722c9175b 100644
--- a/src/operators/multipliedoperator.jl
+++ b/src/operators/multipliedoperator.jl
@@ -1,7 +1,7 @@
 """
     Structure representing a multiplied operator. Consists of
         - An operator op (MPO, Hamiltonian, ...)
-        - An object f that gets multiplied with the operator (Number, function, ...) 
+        - An object f that gets multiplied with the operator (Number, function, ...)
 """
 struct MultipliedOperator{O, F}
     op::O
diff --git a/src/operators/windowhamiltonian.jl b/src/operators/windowhamiltonian.jl
index e95e10208..0fdd83d72 100644
--- a/src/operators/windowhamiltonian.jl
+++ b/src/operators/windowhamiltonian.jl
@@ -10,11 +10,7 @@ an infinite Hamiltonian to the right.
 
 Acts similar to just a finite Hamiltonian, but we "remember" the boundary Hamiltonians.
 
-## Fields
-
-$(TYPEDFIELDS)
-
-## Constructors
+# Constructors
 
     WindowMPOHamiltonian(ham::InfiniteMPOHamiltonian, interval::UnitRange)
 
@@ -23,6 +19,10 @@ Hamiltonian `ham`.
 The finite window consists of the sites in `interval`, while the left and right environments
 are copies of `ham` whose unit cells are circshifted so that they line up with the window
 boundaries.
+
+# Fields
+
+$(TYPEDFIELDS)
 """
 struct WindowMPOHamiltonian{O} <: AbstractMPO{O}
     "Hamiltonian acting on the infinite environment to the left of the window"
diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl
index 0b03e85b4..498daa86a 100644
--- a/src/states/abstractmps.jl
+++ b/src/states/abstractmps.jl
@@ -14,7 +14,7 @@ const GenericMPSTensor{S, N} = AbstractTensorMap{T, S, N, 1} where {T} # some fu
 const MPSTensor{S} = GenericMPSTensor{S, 2} # the usual mps tensors on which we work
 
 """
-    isfullrank(A::GenericMPSTensor; side=:both)
+    isfullrank(A::GenericMPSTensor; side = :both)
 
 Determine whether the given tensor is full rank, i.e. whether both the map from the left
 virtual space and the physical space to the right virtual space, and the map from the right
@@ -37,7 +37,7 @@ function isfullrank(V::TensorKit.TensorMapSpace; side = :both)
 end
 
 """
-    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg=Defaults.alg_orth())
+    makefullrank!(A::PeriodicVector{<:GenericMPSTensor}; alg = Defaults.alg_orth())
 
 Make the set of MPS tensors full rank by performing a series of orthogonalizations.
 """
@@ -83,12 +83,12 @@ end
 # Tensor accessors
 # ----------------
 @doc """
-    AC2(ψ::AbstractMPS, i; kind=:ACAR)
+    AC2(ψ::AbstractMPS, i; kind = :ACAR)
 
 Obtain the two-site (center) gauge tensor at site `i` of the MPS `ψ`.
 If this hasn't been computed before, this can be computed as:
-- `kind=:ACAR` : AC[i] * AR[i+1]
-- `kind=:ALAC` : AL[i] * AC[i+1]
+- `kind = :ACAR`: AC[i] * AR[i+1]
+- `kind = :ALAC`: AL[i] * AC[i+1]
 """ AC2
 
 #===========================================================================================
@@ -130,7 +130,7 @@ TensorKit.sectortype(ψtype::Type{<:AbstractMPS}) = sectortype(site_type(ψtype)
 TensorKit.storagetype(ψtype::Type{<:AbstractMPS}) = storagetype(site_type(ψtype))
 
 """
-    left_virtualspace(ψ::AbstractMPS, [pos=1:length(ψ)])
+    left_virtualspace(ψ::AbstractMPS, [pos = 1:length(ψ)])
     
 Return the virtual space of the bond to the left of sites `pos`.
 
@@ -144,7 +144,7 @@ left_virtualspace(O::MPOTensor) = space(O, 1)
 left_virtualspace(ψ::AbstractMPS) = map(Base.Fix1(left_virtualspace, ψ), eachsite(ψ))
 
 """
-    right_virtualspace(ψ::AbstractMPS, [pos=1:length(ψ)])
+    right_virtualspace(ψ::AbstractMPS, [pos = 1:length(ψ)])
 
 Return the virtual space of the bond to the right of site(s) `pos`.
 
@@ -158,7 +158,7 @@ right_virtualspace(O::MPOTensor) = space(O, 4)'
 right_virtualspace(ψ::AbstractMPS) = map(Base.Fix1(right_virtualspace, ψ), eachsite(ψ))
 
 """
-    physicalspace(ψ::AbstractMPS, [pos=1:length(ψ)])
+    physicalspace(ψ::AbstractMPS, [pos = 1:length(ψ)])
 
 Return the physical space of the site tensor at site `i`.
 """
diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl
index 1b7fa1aa9..c19af5ff0 100644
--- a/src/states/finitemps.jl
+++ b/src/states/finitemps.jl
@@ -1,61 +1,80 @@
 """
-    FiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractFiniteMPS
+$(TYPEDEF)
 
 Type that represents a finite Matrix Product State.
 
-## Properties
-- `AL` -- left-gauged MPS tensors
-- `AR` -- right-gauged MPS tensors
-- `AC` -- center-gauged MPS tensors
-- `C` -- gauge tensors
-- `center` -- location of the gauge center
+# Constructors
 
-The center property returns `center::HalfInt` that indicates the location of the MPS center:
+    FiniteMPS(
+        [f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}}},
+        maxvirtualspaces::Union{S, Vector{S}};
+        normalize = true, left = unitspace(S), right = unitspace(S)
+    ) where {S <: ElementarySpace}
+    FiniteMPS(
+        [f, eltype], N::Int, physicalspace::Union{S, CompositeSpace{S}},
+        maxvirtualspaces::Union{S, Vector{S}};
+        normalize = true, left = unitspace(S), right = unitspace(S)
+    ) where {S <: ElementarySpace}
+    FiniteMPS(As::Vector{<:GenericMPSTensor}; normalize = false, overwrite = false)
+
+Construct an MPS via a specification of physical and virtual spaces, or from a list of
+tensors `As`. All cases reduce to the latter. In particular, a state with a non-trivial
+total charge can be constructed by passing a non-trivially charged vector space as the
+`left` or `right` virtual spaces.
+
+# Arguments
+
+- `As`: vector of site tensors
+- `f = rand`: initializer function for the tensor data
+- `eltype = ComplexF64`: scalar type of the tensors
+- `physicalspaces`: list of physical spaces
+- `N`: number of sites
+- `physicalspace`: local physical space, repeated for every site
+- `maxvirtualspaces`: maximal virtual space(s), truncated to what symmetry allows
+
+# Keyword Arguments
+
+- `normalize`: normalize the constructed state
+- `overwrite = false`: overwrite the given input tensors
+- `left = unitspace(S)`: left-most virtual space
+- `right = unitspace(S)`: right-most virtual space
+
+# Properties
+
+- `AL`: left-gauged MPS tensors
+- `AR`: right-gauged MPS tensors
+- `AC`: center-gauged MPS tensors
+- `C`: gauge (bond) tensors
+- `center`: location of the gauge center
+
+The `center` property returns a `center::HalfInt` that indicates the location of the MPS center:
 - `isinteger(center)` → `center` is a whole number and indicates the location of the first `AC` tensor present in the underlying `ψ.ACs` field.
 - `ishalfodd(center)` → `center` is a half-odd-integer, meaning that there are no `AC` tensors, and indicating between which sites the bond tensor lives.
 
-e.g `mps.center = 7/2` means that the bond tensor is to the right of the 3rd site and can be accessed via `mps.C[3]`.
+For example, `mps.center = 7/2` means that the bond tensor is to the right of the 3rd site and can be accessed via `mps.C[3]`.
+
+# Notes
 
-## Notes
 By convention, we have that:
 - `AL[i] * C[i]` = `AC[i]` = `C[i-1] * AR[i]`
 - `AL[i]' * AL[i] = 1`
 - `AR[i] * AR[i]' = 1`
 
----
-
-## Constructors
-    FiniteMPS([f, eltype], physicalspaces::Vector{<:Union{S,CompositeSpace{S}}},
-              maxvirtualspaces::Union{S,Vector{S}};
-              normalize=true, left=unitspace(S), right=unitspace(S)) where {S<:ElementarySpace}
-    FiniteMPS([f, eltype], N::Int, physicalspace::Union{S,CompositeSpace{S}},
-              maxvirtualspaces::Union{S,Vector{S}};
-              normalize=true, left=unitspace(S), right=unitspace(S)) where {S<:ElementarySpace}
-    FiniteMPS(As::Vector{<:GenericMPSTensor}; normalize=false, overwrite=false)
-
-Construct an MPS via a specification of physical and virtual spaces, or from a list of
-tensors `As`. All cases reduce to the latter. In particular, a state with a non-trivial
-total charge can be constructed by passing a non-trivially charged vector space as the
-`left` or `right` virtual spaces.
-
-### Arguments
-- `As::Vector{<:GenericMPSTensor}`: vector of site tensors
+# Examples
 
-- `f::Function=rand`: initializer function for tensor data
-- `eltype::Type{<:Number}=ComplexF64`: scalar type of tensors
+Building a 3-site spin-1/2 MPS from a dense array and checking that its left-gauged tensors
+are isometries (the state is kept in canonical form even though the raw data is not
+normalized):
 
-- `physicalspaces::Vector{<:Union{S, CompositeSpace{S}}`: list of physical spaces
-- `N::Int`: number of sites
-- `physicalspace::Union{S,CompositeSpace{S}}`: local physical space
+```jldoctest
+julia> ψ = FiniteMPS(ones(Float64, (ℂ^2)^3));
 
-- `virtualspaces::Vector{<:Union{S, CompositeSpace{S}}`: list of virtual spaces
-- `maxvirtualspace::S`: maximum virtual space
+julia> length(ψ)
+3
 
-### Keywords
-- `normalize=true`: normalize the constructed state
-- `overwrite=false`: overwrite the given input tensors
-- `left=unitspace(S)`: left-most virtual space
-- `right=unitspace(S)`: right-most virtual space
+julia> ψ.AL[1]' * ψ.AL[1] ≈ id(left_virtualspace(ψ, 2))
+true
+```
 """
 struct FiniteMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractFiniteMPS
     ALs::Vector{Union{Missing, A}}
@@ -182,7 +201,8 @@ Return the location of the MPS center.
 - `isinteger(center)` → `center` is a whole number and indicates the location of the first `AC` tensor present in `ψ.ACs`
 - `ishalfodd(center)` → `center` is a half-odd-integer, meaning that there are no `AC` tensors, and indicating between which sites the bond tensor lives.
 
-## Example
+# Examples
+
 ```julia
 ψ = FiniteMPS(3, ℂ^2, ℂ^16)
 ψ.center # returns 7/2, bond tensor is to the right of the 3rd site
@@ -441,7 +461,7 @@ end
 
 """
     max_virtualspaces(ψ::FiniteMPS)
-    max_virtualspaces(Ps::Vector{<:Union{S,CompositeSpace{S}}}; left=unitspace(S), right=unitspace(S))
+    max_virtualspaces(Ps::Vector{<:Union{S, CompositeSpace{S}}}; left = unitspace(S), right = unitspace(S))
 
 Compute the maximal virtual spaces of a given finite MPS or its physical spaces.
 """
diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl
index 0d7ae935e..5c791add5 100644
--- a/src/states/infinitemps.jl
+++ b/src/states/infinitemps.jl
@@ -1,47 +1,68 @@
 """
-    InfiniteMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbtractMPS
+$(TYPEDEF)
 
 Type that represents an infinite Matrix Product State.
 
-## Fields
-- `AL` -- left-gauged MPS tensors
-- `AR` -- right-gauged MPS tensors
-- `AC` -- center-gauged MPS tensors
-- `C` -- gauge tensors
+# Constructors
+
+    InfiniteMPS(
+        [f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}}},
+        virtualspaces::Vector{<:Union{S, CompositeSpace{S}}};
+        kwargs...
+    ) where {S <: ElementarySpace}
+    InfiniteMPS(As::AbstractVector{<:GenericMPSTensor}; kwargs...)
+    InfiniteMPS(ALs::AbstractVector{<:GenericMPSTensor}, C₀::MPSBondTensor; kwargs...)
+
+Construct an MPS via a specification of physical and virtual spaces, or from a list of
+tensors `As`, or a list of left-gauged tensors `ALs`.
+
+# Arguments
+
+- `As`: vector of site tensors
+- `ALs`: vector of left-gauged site tensors
+- `C₀`: initial gauge tensor
+- `f = rand`: initializer function for the tensor data
+- `eltype = ComplexF64`: scalar type of the tensors
+- `physicalspaces`: list of physical spaces
+- `virtualspaces`: list of virtual spaces
+
+# Keyword Arguments
+
+- `tol`: gauge fixing tolerance
+- `maxiter`: gauge fixing maximum iterations
+
+# Properties
+
+- `AL`: left-gauged MPS tensors
+- `AR`: right-gauged MPS tensors
+- `AC`: center-gauged MPS tensors
+- `C`: gauge (bond) tensors
+
+# Notes
 
-## Notes
 By convention, we have that:
 - `AL[i] * C[i]` = `AC[i]` = `C[i-1] * AR[i]`
 - `AL[i]' * AL[i] = 1`
 - `AR[i] * AR[i]' = 1`
 
----
+# Examples
 
-## Constructors
-    InfiniteMPS([f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}},
-                virtualspaces::Vector{<:Union{S, CompositeSpace{S}};
-                kwargs...) where {S<:ElementarySpace}
-    InfiniteMPS(As::AbstractVector{<:GenericMPSTensor}; kwargs...)
-    InfiniteMPS(ALs::AbstractVector{<:GenericMPSTensor}, C₀::MPSBondTensor;
-                kwargs...)
+A one-site unit cell built from an explicit `(V_left ⊗ P ← V_right)` tensor. Here the bond
+dimension is one, so this is the `|+⟩` product state, for which `⟨X⟩ = 1`:
 
-Construct an MPS via a specification of physical and virtual spaces, or from a list of
-tensors `As`, or a list of left-gauged tensors `ALs`.
+```jldoctest
+julia> A = TensorMap(ones(Float64, 2, 1), ℂ^1 ⊗ ℂ^2, ℂ^1);
 
-### Arguments
-- `As::AbstractVector{<:GenericMPSTensor}`: vector of site tensors
-- `ALs::AbstractVector{<:GenericMPSTensor}`: vector of left-gauged site tensors
-- `C₀::MPSBondTensor`: initial gauge tensor
+julia> ψ = InfiniteMPS([A]);
 
-- `f::Function=rand`: initializer function for tensor data
-- `eltype::Type{<:Number}=ComplexF64`: scalar type of tensors
+julia> length(ψ)
+1
 
-- `physicalspaces::AbstractVector{<:Union{S, CompositeSpace{S}}`: list of physical spaces
-- `virtualspaces::AbstractVector{<:Union{S, CompositeSpace{S}}`: list of virtual spaces
+julia> X = TensorMap(Float64[0 1; 1 0], ℂ^2, ℂ^2);
 
-### Keywords
-- `tol`: gauge fixing tolerance
-- `maxiter`: gauge fixing maximum iterations
+julia> round(real(expectation_value(ψ, 1 => X)); digits = 6)
+1.0
+```
 """
 struct InfiniteMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractMPS
     AL::PeriodicVector{A}
diff --git a/src/states/multilinemps.jl b/src/states/multilinemps.jl
index eda3c5c35..524280ffc 100644
--- a/src/states/multilinemps.jl
+++ b/src/states/multilinemps.jl
@@ -5,18 +5,31 @@ const MultilineMPS = Multiline{<:InfiniteMPS}
 @doc """
     const MultilineMPS = Multiline{<:InfiniteMPS}
 
-Type that represents multiple lines of `InfiniteMPS` objects.
+Type that represents multiple lines of [`InfiniteMPS`](@ref) objects.
 
 # Constructors
+
     MultilineMPS(mpss::AbstractVector{<:InfiniteMPS})
-    MultilineMPS([f, eltype], physicalspaces::Matrix{<:Union{S, CompositeSpace{S}},
-                 virtualspaces::Matrix{<:Union{S, CompositeSpace{S}}) where
-                 {S<:ElementarySpace}
+    MultilineMPS(
+        [f, eltype], physicalspaces::Matrix{<:Union{S, CompositeSpace{S}}},
+        virtualspaces::Matrix{<:Union{S, CompositeSpace{S}}}
+    ) where {S <: ElementarySpace}
     MultilineMPS(As::AbstractMatrix{<:GenericMPSTensor}; kwargs...)
-    MultilineMPS(ALs::AbstractMatrix{<:GenericMPSTensor},
-                 C₀::AbstractVector{<:MPSBondTensor}; kwargs...)
+    MultilineMPS(
+        ALs::AbstractMatrix{<:GenericMPSTensor},
+        C₀::AbstractVector{<:MPSBondTensor}; kwargs...
+    )
+
+# Properties
+
+- `AL`: left-gauged MPS tensors
+- `AR`: right-gauged MPS tensors
+- `AC`: center-gauged MPS tensors
+- `C`: gauge (bond) tensors
+
+# See also
 
-See also: [`Multiline`](@ref)
+[`Multiline`](@ref)
 """
 function MultilineMPS end
 
diff --git a/src/states/ortho.jl b/src/states/ortho.jl
index d7901de08..9dfc803e4 100644
--- a/src/states/ortho.jl
+++ b/src/states/ortho.jl
@@ -8,9 +8,13 @@ $(TYPEDEF)
 
 Algorithm for bringing an `InfiniteMPS` into the left-canonical form.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `alg` argument of [`gaugefix!`](@ref).
 """
 @kwdef struct LeftCanonical <: Algorithm
     "tolerance for convergence criterium"
@@ -33,9 +37,13 @@ $(TYPEDEF)
 
 Algorithm for bringing an `InfiniteMPS` into the right-canonical form.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `alg` argument of [`gaugefix!`](@ref).
 """
 @kwdef struct RightCanonical <: Algorithm
     "tolerance for convergence criterium"
@@ -58,9 +66,13 @@ $(TYPEDEF)
 
 Algorithm for bringing an `InfiniteMPS` into the mixed-canonical form.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
+
+# See also
+
+Used as the `alg` argument of [`gaugefix!`](@ref).
 """
 struct MixedCanonical <: Algorithm
     "algorithm for bringing an `InfiniteMPS` into left-canonical form."
@@ -153,7 +165,7 @@ end
 Bring updated `AC` and `C` tensors back into a consistent set of left or right canonical
 tensors. This minimizes `∥AC_i - AL_i * C_i∥` or `∥AC_i - C_{i-1} * AR_i∥`.
 
-The `alg` is passed on to `left_orth!` and `right_orth!`, and can be used to control the kind of 
+The `alg` is passed on to `left_orth!` and `right_orth!`, and can be used to control the kind of
 factorization used. By default, this is set to a (positive) QR/LQ, even though the
 optimal algorithm would use a polar decompositions instead, sacrificing a bit of
 performance for accuracy.
diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl
index 83ef7553c..6934e2b22 100644
--- a/src/states/quasiparticle_state.jl
+++ b/src/states/quasiparticle_state.jl
@@ -4,6 +4,38 @@ I think it makes sense to see these things as an actual state instead of return
 This will allow us to plot energy density (finite qp) and measure observables.
 =#
 
+"""
+$(TYPEDEF)
+
+Left-gauged quasiparticle excitation ansatz on top of a matrix product state ground state.
+The excitation is parametrized through the left-gauge nullspace of the ground-state tensors,
+and the object behaves as a vector so it can be handed directly to the iterative eigensolvers
+used by [`excitations`](@ref).
+
+For a `FiniteMPS` ground state this represents a finite (localized) quasiparticle; for an
+`InfiniteMPS` ground state it represents a momentum eigenstate with the given `momentum`.
+When `left_gs !== right_gs` the ansatz describes a domain wall between the two ground states.
+
+# Constructors
+
+    LeftGaugedQP(datfun, left_gs, right_gs = left_gs; sector, momentum = 0.0)
+
+These states are normally produced by [`excitations`](@ref) with a
+[`QuasiparticleAnsatz`](@ref) rather than constructed directly. When constructing manually,
+`datfun` initializes the variational tensors (e.g. `rand`/`randn`), `sector` selects the
+charge sector of the excitation, and `momentum` sets the momentum for infinite ground states.
+
+# Fields
+
+- `left_gs`, `right_gs`: the ground state(s) the excitation lives on; distinct values yield a domain wall.
+- `VLs`: left-nullspace tensors of the ground-state `AL` (satisfying `AL' * VL == 0`).
+- `Xs`: the variational parameters of the ansatz.
+- `momentum`: the excitation momentum (used for infinite ground states).
+
+# See also
+
+[`RightGaugedQP`](@ref), [`QP`](@ref)
+"""
 struct LeftGaugedQP{S, T1, T2, E <: Number}
     # !(left_gs === right_gs) => domain wall excitation
     left_gs::S
@@ -15,6 +47,29 @@ struct LeftGaugedQP{S, T1, T2, E <: Number}
     momentum::E
 end
 
+"""
+$(TYPEDEF)
+
+Right-gauged counterpart of [`LeftGaugedQP`](@ref): the same quasiparticle excitation ansatz,
+but parametrized through the right-gauge nullspace of the ground-state tensors. It is most
+often obtained via `convert(RightGaugedQP, ϕ)` from a `LeftGaugedQP` rather than constructed
+directly.
+
+# Constructors
+
+    RightGaugedQP(datfun, left_gs, right_gs = left_gs; sector, momentum = 0.0)
+
+# Fields
+
+- `left_gs`, `right_gs`: the ground state(s) the excitation lives on; distinct values yield a domain wall.
+- `Xs`: the variational parameters of the ansatz.
+- `VRs`: right-nullspace tensors of the ground-state `AR`.
+- `momentum`: the excitation momentum (used for infinite ground states).
+
+# See also
+
+[`LeftGaugedQP`](@ref), [`QP`](@ref)
+"""
 struct RightGaugedQP{S, T1, T2, E <: Number}
     # !(left_gs === right_gs) => domain wall excitation
     left_gs::S
@@ -207,6 +262,15 @@ function Base.convert(
 end
 
 # gauge independent code
+"""
+    QP{S, T1, T2}
+
+Union of the quasiparticle excitation ansätze [`LeftGaugedQP`](@ref) and
+[`RightGaugedQP`](@ref). It is used for dispatch and to share their gauge-independent
+interface; it is not a concrete type and cannot be constructed on its own. The internal
+aliases `FiniteQP` and `InfiniteQP` further restrict the ground-state type to `FiniteMPS`
+or `InfiniteMPS` respectively.
+"""
 const QP{S, T1, T2} = Union{LeftGaugedQP{S, T1, T2}, RightGaugedQP{S, T1, T2}}
 const FiniteQP{S <: FiniteMPS, T1, T2} = QP{S, T1, T2}
 const InfiniteQP{S <: InfiniteMPS, T1, T2} = QP{S, T1, T2}
diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl
index e32ac35d2..887109440 100644
--- a/src/states/windowmps.jl
+++ b/src/states/windowmps.jl
@@ -1,39 +1,42 @@
 """
-    WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractFiniteMPS
+$(TYPEDEF)
 
 Type that represents a finite Matrix Product State embedded in an infinite Matrix Product State.
 
-## Fields
-
-- `left_gs::InfiniteMPS` -- left infinite environment
-- `window::FiniteMPS` -- finite window Matrix Product State
-- `right_gs::InfiniteMPS` -- right infinite environment
-
----
-
-## Constructors
+# Constructors
 
     WindowMPS(left_gs::InfiniteMPS, window_state::FiniteMPS, [right_gs::InfiniteMPS])
     WindowMPS(left_gs::InfiniteMPS, window_tensors::AbstractVector, [right_gs::InfiniteMPS])
-    WindowMPS([f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}},
-              virtualspaces::Vector{<:Union{S, CompositeSpace{S}}, left_gs::InfiniteMPS,
-              [right_gs::InfiniteMPS])
-    WindowMPS([f, eltype], physicalspaces::Vector{<:Union{S,CompositeSpace{S}}},
-              maxvirtualspace::S, left_gs::InfiniteMPS, [right_gs::InfiniteMPS])
+    WindowMPS(
+        [f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}}},
+        virtualspaces::Vector{<:Union{S, CompositeSpace{S}}}, left_gs::InfiniteMPS,
+        [right_gs::InfiniteMPS]
+    )
+    WindowMPS(
+        [f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}}},
+        maxvirtualspace::S, left_gs::InfiniteMPS, [right_gs::InfiniteMPS]
+    )
+    WindowMPS(ψ::InfiniteMPS, L::Int)
 
 Construct a WindowMPS via a specification of left and right infinite environment, and either
 a window state or a vector of tensors to construct the window. Alternatively, it is possible
 to supply the same arguments as for the constructor of [`FiniteMPS`](@ref), followed by a
-left (and right) environment to construct the WindowMPS in one step.
+left (and right) environment to construct the WindowMPS in one step. Finally, a WindowMPS can
+be constructed from an `InfiniteMPS` by promoting a region of length `L` to a `FiniteMPS`.
 
 !!! note
     By default, the right environment is chosen to be equal to the left, however no copy is
     made. In this case, changing the left state will also affect the right state.
 
-    WindowMPS(state::InfiniteMPS, L::Int)
+# Properties
 
-Construct a WindowMPS from an InfiniteMPS, by promoting a region of length `L` to a
-`FiniteMPS`.
+- `left_gs::InfiniteMPS`: left infinite environment
+- `window::FiniteMPS`: finite window Matrix Product State
+- `right_gs::InfiniteMPS`: right infinite environment
+- `AL`: left-gauged MPS tensors
+- `AR`: right-gauged MPS tensors
+- `AC`: center-gauged MPS tensors
+- `C`: gauge (bond) tensors
 """
 struct WindowMPS{A <: GenericMPSTensor, B <: MPSBondTensor} <: AbstractFiniteMPS
     left_gs::InfiniteMPS{A, B}
diff --git a/src/utility/dynamictols.jl b/src/utility/dynamictols.jl
index fe52d3188..5ef4cc9f1 100644
--- a/src/utility/dynamictols.jl
+++ b/src/utility/dynamictols.jl
@@ -44,11 +44,13 @@ Algorithm wrapper with dynamically adjusted tolerances. Only the wrapped solver'
 retuned; its Krylov budget (if any) is left fixed — this is the simpler counterpart to
 [`AdaptiveKrylov`](@ref).
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-See also [`adapt_solver`](@ref).
+# See also
+
+[`adapt_solver`](@ref)
 """
 struct DynamicTol{A} <: Algorithm
     "parent algorithm"
@@ -79,7 +81,7 @@ end
 Tighten only the wrapped solver's tolerance (its Krylov budget, if any, is left fixed), from the
 global gradient / convergence-error scalar `g_global`, damped by the iteration count:
 
-    tol = clamp(tol_factor·g_global / √iter, tol_min, tol_max)
+    tol = clamp(tol_factor · g_global / √iter, tol_min, tol_max)
 """
 function adapt_solver(alg::DynamicTol; iter::Integer = 1, g_global::Real = 0.0, kwargs...)
     tol = clamp(alg.tol_factor * g_global / sqrt(max(iter, 1)), alg.tol_min, alg.tol_max)
@@ -117,11 +119,13 @@ This is driven by the local and global gradient norm, the truncation error and t
 decay rate of previous iterations in an attempt to obtain fast convergence for gapped systems
 while avoiding stagnation for gapless ones.
 
-## Fields
+# Fields
 
 $(TYPEDFIELDS)
 
-See also [`adapt_solver`](@ref).
+# See also
+
+[`adapt_solver`](@ref)
 """
 struct AdaptiveKrylov{T, O <: KrylovKit.Orthogonalizer} <: Algorithm
     "orthogonalizer passed to the instantiated `Lanczos`/`Arnoldi`"
diff --git a/src/utility/multiline.jl b/src/utility/multiline.jl
index 53b02cfc6..ed4107cf2 100644
--- a/src/utility/multiline.jl
+++ b/src/utility/multiline.jl
@@ -1,13 +1,16 @@
 """
-    struct Multiline{T}
+$(TYPEDEF)
 
 Object that represents multiple lines of objects of type `T`. Typically used to represent
 multiple lines of `InfiniteMPS` (`MultilineMPS`) or MPO (`Multiline{<:AbstractMPO}`).
 
 # Fields
-- `data::PeriodicArray{T,1}`: the data of the multiline object
 
-See also: [`MultilineMPS`](@ref) and [`MultilineMPO`](@ref)
+- `data::PeriodicArray{T, 1}`: the data of the multiline object
+
+# See also
+
+[`MultilineMPS`](@ref) and [`MultilineMPO`](@ref)
 """
 struct Multiline{T}
     data::PeriodicArray{T, 1}
diff --git a/src/utility/periodicarray.jl b/src/utility/periodicarray.jl
index 14ce75b65..0131d10e6 100644
--- a/src/utility/periodicarray.jl
+++ b/src/utility/periodicarray.jl
@@ -1,12 +1,14 @@
 """
-    PeriodicArray{T,N} <: AbstractArray{T,N}
+$(TYPEDEF)
 
 Array wrapper with periodic boundary conditions.
 
 # Fields
-- `data::Array{T,N}`: the data of the array
+
+- `data::Array{T, N}`: the data of the array
 
 # Examples
+
 ```jldoctest
 A = PeriodicArray([1, 2, 3])
 A[0], A[2], A[4]
@@ -24,7 +26,9 @@ A[-1, 1], A[1, 1], A[4, 5]
 (1, 1, 3)
 ```
 
-See also [`PeriodicVector`](@ref), [`PeriodicMatrix`](@ref)
+# See also
+
+[`PeriodicVector`](@ref), [`PeriodicMatrix`](@ref)
 """
 struct PeriodicArray{T, N} <: AbstractArray{T, N}
     data::Array{T, N}
@@ -42,7 +46,7 @@ end
     PeriodicVector{T}
 
 One-dimensional dense array with elements of type `T` and periodic boundary conditions.
-Alias for [`PeriodicArray{T,1}`](@ref).
+Alias for [`PeriodicArray{T, 1}`](@ref).
 """
 const PeriodicVector{T} = PeriodicArray{T, 1}
 PeriodicVector(data::AbstractVector{T}) where {T} = PeriodicVector{T}(data)
@@ -51,7 +55,7 @@ PeriodicVector(data::AbstractVector{T}) where {T} = PeriodicVector{T}(data)
     PeriodicMatrix{T}
 
 Two-dimensional dense array with elements of type `T` and periodic boundary conditions.
-Alias for [`PeriodicArray{T,2}`](@ref).
+Alias for [`PeriodicArray{T, 2}`](@ref).
 """
 const PeriodicMatrix{T} = PeriodicArray{T, 2}
 PeriodicMatrix(data::AbstractMatrix{T}) where {T} = PeriodicMatrix{T}(data)
diff --git a/src/utility/plotting.jl b/src/utility/plotting.jl
index 022fee33b..ec491693e 100644
--- a/src/utility/plotting.jl
+++ b/src/utility/plotting.jl
@@ -1,18 +1,20 @@
 """
-    entanglementplot(state; site=0[, kwargs...])
+    entanglementplot(state; site = 0[, kwargs...])
 
-Plot the [entanglement spectrum](@ref entanglement_spectrum) of a given MPS `state`. 
+Plot the entanglement spectrum (see [`entanglement_spectrum`](@ref)) of a given MPS `state`.
 
 # Arguments
+
 - `state`: the MPS for which to compute the entanglement spectrum.
 
 # Keyword Arguments
-- `site::Int=0`: MPS index for multisite unit cells. The spectrum is computed for the bond
+
+- `site::Int = 0`: MPS index for multisite unit cells. The spectrum is computed for the bond
   between `site` and `site + 1`.
-- `expand_symmetry::Logical=false`: add quantum dimension degeneracies.
-- `sortby=maximum`: the method of sorting the sectors.
-- `sector_margin=1//10`: the amount of whitespace between sectors.
-- `sector_formatter=string`: how to convert sectors to strings.
+- `expand_symmetry = false`: add quantum dimension degeneracies.
+- `sortby = maximum`: the method of sorting the sectors.
+- `sector_margin = 1 // 10`: the amount of whitespace between sectors.
+- `sector_formatter = string`: how to convert sectors to strings.
 - `kwargs...`: other kwargs are passed on to the plotting backend.
 
 !!! note
@@ -89,21 +91,23 @@ function entanglementplot end
 end
 
 """
-    transferplot(above, below=above; sectors=nothing, transferkwargs=(;)[, kwargs...])
+    transferplot(above, below = above; sectors = nothing, transferkwargs = (;)[, kwargs...])
 
 Plot the partial transfer matrix spectrum of two InfiniteMPS's.
 
 # Arguments
+
 - `above::InfiniteMPS`: above mps for [`transfer_spectrum`](@ref).
-- `below::InfiniteMPS=above`: below mps for [`transfer_spectrum`](@ref).
+- `below::InfiniteMPS = above`: below mps for [`transfer_spectrum`](@ref).
 
 # Keyword Arguments
-- `sectors=nothing`: restrict the spectrum to the given sectors; by default all sectors of
+
+- `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.
-- `sector_formatter=string`: how to convert sectors to strings.
+- `thetaorigin = 0`: origin of the angle range.
+- `sector_formatter = string`: how to convert sectors to strings.
 
 !!! note
     You will need to manually import [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to
diff --git a/src/utility/show.jl b/src/utility/show.jl
index bd132dd6c..2709b7d43 100644
--- a/src/utility/show.jl
+++ b/src/utility/show.jl
@@ -196,6 +196,7 @@ Each site of the MPO is represented as a block of Unicode braille characters, wi
 This visualization is useful for quickly inspecting the structure and sparsity pattern of MPOs.
 
 # Arguments
+
 - `io::IO`: The output stream to print to (e.g., `stdout`).
 - `H::Union{SparseMPO, MPOHamiltonian}`: The `SparseMPO` or `MPOHamiltonian` to visualize.
 
diff --git a/src/utility/styles.jl b/src/utility/styles.jl
index 962859019..df83d8e55 100644
--- a/src/utility/styles.jl
+++ b/src/utility/styles.jl
@@ -4,8 +4,8 @@
     OperatorStyle(::Type{T})
 
 Trait to describe the operator behavior of the input `x` or type `T`, which can be either
-*   `MPOStyle()`: product of local factors;
-*   `HamiltonianStyle()`: sum of local terms.
+- `MPOStyle()`: product of local factors;
+- `HamiltonianStyle()`: sum of local terms.
 """
 abstract type OperatorStyle end
 OperatorStyle(x) = OperatorStyle(typeof(x))
@@ -29,8 +29,8 @@ struct HamiltonianStyle <: OperatorStyle end
     GeometryStyle(::Type{T})
 
 Trait to describe the geometry of the input `x` or type `T`, which can be either
-*   `FiniteChainStyle()`: object is defined on a finite chain;
-*   `InfiniteChainStyle()`: object is defined on an infinite chain.
+- `FiniteChainStyle()`: object is defined on a finite chain;
+- `InfiniteChainStyle()`: object is defined on an infinite chain.
 """
 abstract type GeometryStyle end
 GeometryStyle(x) = GeometryStyle(typeof(x))