diff --git a/Project.toml b/Project.toml
index f6c5866f0..377ddb3d8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,9 +3,6 @@ uuid = "bb1c41ca-d63c-52ed-829e-0820dda26502"
 version = "0.13.13"
 authors = "Lukas Devos, Maarten Van Damme and contributors"
 
-[workspace]
-projects = ["test", "docs", "examples"]
-
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
@@ -21,7 +18,6 @@ OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 TensorKitManifolds = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
@@ -30,9 +26,17 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [weakdeps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 
 [extensions]
 MPSKitAdaptExt = "Adapt"
+MPSKitMakieExt = ["Makie", "LaTeXStrings"]
+MPSKitPlotsExt = "RecipesBase"
+
+[workspace]
+projects = ["test", "docs", "examples"]
 
 [compat]
 Accessors = "0.1"
@@ -42,18 +46,20 @@ Compat = "3.47, 4.10"
 DocStringExtensions = "0.9.3"
 HalfIntegers = "1.6.0"
 KrylovKit = "0.8.3, 0.9.2, 0.10"
+LaTeXStrings = "1"
 LinearAlgebra = "1.6"
 LoggingExtras = "~1.0"
+Makie = "0.24, 0.25"
 MatrixAlgebraKit = "0.6.5"
 OhMyThreads = "0.7, 0.8"
 OptimKit = "0.4, 0.5"
 Preferences = "1"
 Printf = "1"
 Random = "1"
-RecipesBase = "1.1"
+RecipesBase = "1"
 TensorKit = "0.17.1"
 TensorKitManifolds = "0.7, 0.8"
 TensorOperations = "5.6.2"
 TimerOutputs = "1"
 VectorInterface = "0.2, 0.3, 0.4, 0.5, 0.6"
-julia = "1.10"
+julia = "1.10"
\ No newline at end of file
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 893a670e7..551907f5e 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -38,6 +38,12 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
   by `MPSKit.default_allocator`, instead of leaving them to the garbage collector
   (two-site DMRG: -64% allocations, -57% GC time, -23% wall time).
   Disable with `MPSKit.Defaults.set_buffering!(false)`. ([#467](https://github.com/QuantumKitHub/MPSKit.jl/pull/467))
+- Makie.jl support for `entanglementplot` and `transferplot`, through a new `MPSKitMakieExt`
+  package extension. Loading any Makie backend (e.g. `using CairoMakie`) alongside MPSKit
+  enables both plot functions, as an alternative to Plots.jl. Backend-specific styling can be
+  passed via the `plotkwargs` keyword. ([#428](https://github.com/QuantumKitHub/MPSKit.jl/pull/428))
+- `entanglementplot!` and `transferplot!`, the mutating counterparts of the plotting functions,
+  which add a spectrum to an existing figure instead of creating a new one. ([#428](https://github.com/QuantumKitHub/MPSKit.jl/pull/428))
 
 ### Changed
 
diff --git a/ext/MPSKitMakieExt.jl b/ext/MPSKitMakieExt.jl
new file mode 100644
index 000000000..34341dd0f
--- /dev/null
+++ b/ext/MPSKitMakieExt.jl
@@ -0,0 +1,288 @@
+module MPSKitMakieExt
+
+using Makie, LaTeXStrings
+using MPSKit, TensorKit
+
+#TODO?: add Colors.jl to access this, allows Plots extension to also use these colors
+const JLCOLORS = Makie.Colors.JULIA_LOGO_COLORS
+
+sector_color(i::Integer) = JLCOLORS[mod1(i, length(JLCOLORS))]
+
+convert_kwargs(kwargs::NamedTuple) = kwargs
+function convert_kwargs(kwargs) # weird convert thing
+    return NamedTuple(Symbol(k) => (v isa Observable ? v[] : v) for (k, v) in pairs(kwargs))
+end
+
+# the recipes publish the axis attributes they want as an `:axis_info` node instead of `current_axis()`
+function apply_axis_info!(ax, plot, plotkwargs = (;))
+    ax isa Makie.AbstractAxis || return ax
+    function apply!(info)
+        for (k, v) in pairs(info)
+            setproperty!(ax, k, v)
+        end
+        # user-provided attributes take precedence
+        for (k, v) in pairs(plotkwargs)
+            setproperty!(ax, k, v)
+        end
+        return nothing
+    end
+    node = plot.attributes[:axis_info]
+    apply!(node[])
+    on(apply!, node) # keep the axis in sync when the inputs change
+    return ax
+end
+
+@recipe EntanglementPlot (mps,) begin
+    site = 0
+    expand_symmetry = false
+    sortby = maximum
+    sector_margin = 1 // 10
+    sector_formatter = string
+    markersize = 12
+    marker = :circle
+end
+
+function Makie.plot!(ep::EntanglementPlot)
+    # this closure only reruns when one of the inputs changes
+    map!(ep.attributes, [:mps, :site, :expand_symmetry, :sortby], :spectrum_data) do mps, site, expand_symmetry, sortby
+        spectra = entanglement_spectrum(mps, site)
+
+        sectors = sectortype(mps)[]
+        spectrum = Vector{Vector{Float64}}()
+        for (c, b) in pairs(spectra)
+            if expand_symmetry
+                b′ = repeat(collect(b), dim(c))
+                sort!(b′; rev = true)
+            else
+                b′ = collect(b)
+            end
+            push!(spectrum, b′)
+            push!(sectors, c)
+        end
+
+        if any(v -> any(<=(0), v), spectrum)
+            @warn "Entanglement spectrum contains vanishing Schmidt values. These are omitted from the plot."
+            foreach(v -> filter!(>(0), v), spectrum)
+        end
+
+        # Sort sectors according to provided method
+        if length(spectrum) > 1
+            order = sortperm(spectrum; by = sortby, rev = true)
+            spectrum = spectrum[order]
+            sectors = sectors[order]
+        end
+        return (; sectors, spectrum)
+    end
+
+    # styling
+    # only reruns when spectrum_data or sector_margin changes
+    map!(ep.attributes, [:spectrum_data, :sector_margin], :positions) do data, margin
+        points = Point2d[]
+        for (i, partial_spectrum) in enumerate(data.spectrum)
+            n_spectrum = length(partial_spectrum)
+            xs = if n_spectrum == 1
+                range(i + 0.5, i + 0.5; length = 1)
+            else
+                range(i + float(margin), i + 1 - float(margin); length = n_spectrum)
+            end
+            for (x, λ) in zip(xs, partial_spectrum)
+                push!(points, Point2d(x, log10(λ)))
+            end
+        end
+        return points
+    end
+
+    # only reruns when spectrum_data changes
+    map!(ep.attributes, [:spectrum_data], :colors) do data
+        colors = typeof(sector_color(1))[]
+        for (i, partial_spectrum) in enumerate(data.spectrum)
+            append!(colors, fill(sector_color(i), length(partial_spectrum)))
+        end
+        return colors
+    end
+
+    map!(ep.attributes, [:mps, :site, :spectrum_data, :sector_formatter], :axis_info) do mps, site, data, sector_formatter
+        nsectors = length(data.sectors)
+        bottom = if isempty(data.spectrum)
+            -1
+        else
+            smallest = minimum(Iterators.flatten(data.spectrum); init = 1.0) # spectrum is already > 0
+            floor(Int, log10(smallest))
+        end
+        return (;
+            xticks = (1:nsectors, sector_formatter.(data.sectors)),
+            xticklabelsize = 16,
+            xticklabelrotation = 45.0,
+            xticklabelalign = (:right, :top),
+            ylabel = L"\log(\lambda)",
+            ylabelsize = 24,
+            yticks = (bottom:2:0, latexstring.(collect(bottom:2:0))),
+            yticklabelsize = 16,
+            limits = ((1, nsectors + 1), (bottom, 0 + 1.0e-1)),
+        )
+    end
+
+    scatter!(ep, ep.positions; color = ep.colors, markersize = ep.markersize, marker = ep.marker)
+    return ep
+end
+
+function MPSKit.entanglementplot(args...; plotkwargs = (;), kwargs...)
+    p = entanglementplot(args...; kwargs...)
+    apply_axis_info!(p.axis, p.plot, plotkwargs)
+    return p
+end
+
+function MPSKit.entanglementplot!(state::MPSKit.AbstractMPS; plotkwargs = (;), kwargs...)
+    p = entanglementplot!(state; kwargs...)
+    apply_axis_info!(Makie.current_axis(), p, plotkwargs)
+    return p
+end
+function MPSKit.entanglementplot!(target, state::MPSKit.AbstractMPS; plotkwargs = (;), kwargs...)
+    p = entanglementplot!(target, state; kwargs...)
+    apply_axis_info!(target, p, plotkwargs)
+    return p
+end
+
+#------------------------------------------------------------
+
+@recipe TransferPlot (above, below) begin
+    sectors = nothing
+    transferkwargs = NamedTuple()
+    thetaorigin = 0.0
+    sector_formatter = string
+    legend_position = :ct
+    markersize = 12
+    marker = :circle
+end
+
+function Makie.plot!(tp::TransferPlot)
+    #TODO: consider radial plot
+    # this only reruns when one of the inputs changes
+    map!(tp.attributes, [:above, :below, :sectors, :transferkwargs], :spectrum_data) do above, below, sectors, transferkwargs
+        kwargs = convert_kwargs(transferkwargs)
+        if sectors !== nothing && get(kwargs, :howmany, 20) isa Int
+            # restrict the computation to the requested sectors
+            howmany = Dict(c => get(kwargs, :howmany, 20) for c in sectors)
+            kwargs = (; kwargs..., howmany)
+        end
+        spectra = transfer_spectrum(above, below; kwargs...)
+
+        data = Pair{sectortype(above), Vector{complex(scalartype(above))}}[]
+        for (sector, spectrum) in pairs(spectra)
+            sectors === nothing || sector in sectors || continue
+            push!(data, sector => collect(spectrum))
+        end
+        return data
+    end
+
+    map!(tp.attributes, [:spectrum_data, :thetaorigin], :positions) do data, thetaorigin
+        points = Point2d[]
+        for (_, spectrum) in data, λ in spectrum
+            θ = mod2pi(angle(λ) + thetaorigin) - thetaorigin
+            push!(points, Point2d(θ, abs(λ)))
+        end
+        return points
+    end
+
+    map!(tp.attributes, [:spectrum_data], :colors) do data
+        colors = typeof(sector_color(1))[]
+        for (i, (_, spectrum)) in enumerate(data)
+            append!(colors, fill(sector_color(i), length(spectrum)))
+        end
+        return colors
+    end
+
+    map!(tp.attributes, [:thetaorigin], :axis_info) do thetaorigin
+        return (;
+            xlabel = L"\theta",
+            xlabelsize = 24,
+            xticks = pitick(0, 2pi, 4; mode = :latex),
+            xticklabelsize = 16,
+            ylabel = L"r",
+            ylabelsize = 24,
+            yticks = (range(0, 1.0; length = 6), latexstring.(range(0, 1.0; length = 6))),
+            yticklabelsize = 16,
+            xgridvisible = true,
+            ygridvisible = true,
+            leftspinevisible = true,
+            rightspinevisible = false,
+            bottomspinevisible = true,
+            topspinevisible = false,
+            limits = ((thetaorigin - 0.1, thetaorigin + 2π + 0.1), (nothing, 1.05)),
+        )
+    end
+
+    map!(tp.attributes, [:spectrum_data, :sector_formatter], :legend_entries) do data, sector_formatter
+        return [(sector_formatter(sector), sector_color(i)) for (i, (sector, _)) in enumerate(data)]
+    end
+
+    scatter!(tp, tp.positions; color = tp.colors, markersize = tp.markersize, marker = tp.marker)
+    return tp
+end
+
+function add_sector_legend!(ax, plot, legend_position)
+    ax isa Makie.AbstractAxis || return nothing
+    entries = plot.attributes[:legend_entries][]
+    isempty(entries) && return nothing
+    elements = [MarkerElement(; color, marker = :circle, markersize = 12) for (_, color) in entries]
+    # cannot use current_figure() when supporting in-place method
+    axislegend(ax, elements, [label for (label, _) in entries]; position = legend_position)
+    return nothing
+end
+
+function MPSKit.transferplot(above, below = above; plotkwargs = (;), kwargs...)
+    p = transferplot(above, below; kwargs...)
+    apply_axis_info!(p.axis, p.plot, plotkwargs)
+    add_sector_legend!(p.axis, p.plot, p.plot.legend_position[])
+    return p
+end
+
+function MPSKit.transferplot!(
+        above::MPSKit.AbstractMPS, below::MPSKit.AbstractMPS = above;
+        plotkwargs = (;), kwargs...
+    )
+    p = transferplot!(above, below; kwargs...)
+    ax = Makie.current_axis()
+    apply_axis_info!(ax, p, plotkwargs)
+    add_sector_legend!(ax, p, p.legend_position[])
+    return p
+end
+function MPSKit.transferplot!(
+        target, above::MPSKit.AbstractMPS, below::MPSKit.AbstractMPS = above;
+        plotkwargs = (;), kwargs...
+    )
+    p = transferplot!(target, above, below; kwargs...)
+    apply_axis_info!(target, p, plotkwargs)
+    add_sector_legend!(target, p, p.legend_position[])
+    return p
+end
+
+# utility for plotting
+
+function pitick(start, stop, denom; mode = :latex)
+    a = Int(cld(start, π / denom))
+    b = Int(fld(stop, π / denom))
+    tick = range(a * π / denom, b * π / denom; step = π / denom)
+    ticklabel = piticklabel.((a:b) .// denom, Val(mode))
+    return tick, ticklabel
+end
+
+function piticklabel(x::Rational, ::Val{:text})
+    iszero(x) && return "0"
+    S = x < 0 ? "-" : ""
+    n, d = abs(numerator(x)), denominator(x)
+    N = n == 1 ? "" : repr(n)
+    d == 1 && return S * N * "π"
+    return S * N * "π/" * repr(d)
+end
+
+function piticklabel(x::Rational, ::Val{:latex})
+    iszero(x) && return L"0"
+    S = x < 0 ? "-" : ""
+    n, d = abs(numerator(x)), denominator(x)
+    N = n == 1 ? "" : repr(n)
+    d == 1 && return L"%$S%$N\pi"
+    return L"%$S\frac{%$N\pi}{%$d}"
+end
+
+end
diff --git a/ext/MPSKitPlotsExt.jl b/ext/MPSKitPlotsExt.jl
new file mode 100644
index 000000000..464ab2efd
--- /dev/null
+++ b/ext/MPSKitPlotsExt.jl
@@ -0,0 +1,122 @@
+module MPSKitPlotsExt
+
+using RecipesBase
+using MPSKit, TensorKit
+
+@userplot EntanglementPlot
+
+@recipe function f(
+        h::EntanglementPlot; site = 0, expand_symmetry = false, sortby = maximum,
+        sector_margin = 1 // 10, sector_formatter = string
+    )
+    mps = h.args[1]
+
+    spectra = entanglement_spectrum(mps, site)
+    sectors = sectortype(mps)[]
+    spectrum = Vector{Vector{Float64}}()
+    for (c, b) in pairs(spectra)
+        if expand_symmetry # Duplicate entries according to the quantum dimension.
+            b′ = repeat(b, dim(c))
+            sort!(b′; rev = true)
+        else
+            b′ = collect(b)
+        end
+        push!(spectrum, b′)
+        push!(sectors, c)
+    end
+
+    if any(v -> any(<=(0), v), spectrum)
+        @warn "Entanglement spectrum contains vanishing Schmidt values. These are omitted from the plot."
+        foreach(v -> filter!(>(0), v), spectrum)
+    end
+
+    if length(spectrum) > 1
+        order = sortperm(spectrum; by = sortby, rev = true)
+        spectrum = spectrum[order]
+        sectors = sectors[order]
+    end
+
+    smallest = minimum(Iterators.flatten(spectrum); init = 1.0) # spectrum is already > 0
+    bottom = floor(Int, log10(smallest))
+
+    for (i, (partial_spectrum, sector)) in enumerate(zip(spectrum, sectors))
+        @series begin
+            legend --> false
+            grid --> :xy
+            widen --> true
+            bottom_margin --> (10, :mm)
+
+            xticks --> (1:length(sectors), sector_formatter.(sectors))
+            xtickfonthalign --> :center
+            xtick_direction --> :out
+            xrotation --> 45
+            xlims --> (1, length(sectors) + 1)
+
+            ylims --> (exp10(bottom), 1 + 1.0e-1)
+            yscale --> :log10
+            seriestype := :scatter
+            label := sector_formatter(sector)
+            n_spectrum = length(partial_spectrum)
+
+            # Put single dot in the middle, or a linear range with padding.
+            if n_spectrum == 1
+                x = [i + 1 // 2]
+            else
+                x = range(i + sector_margin, i + 1 - sector_margin; length = n_spectrum)
+            end
+            return x, partial_spectrum
+        end
+    end
+
+    return nothing
+end
+
+MPSKit.entanglementplot(args...; plotkwargs = (;), kwargs...) = entanglementplot(args...; kwargs..., plotkwargs...)
+MPSKit.entanglementplot!(args...; plotkwargs = (;), kwargs...) = entanglementplot!(args...; kwargs..., plotkwargs...)
+
+#-----------------------------------------------------------------------------
+
+@userplot TransferPlot
+
+@recipe function f(
+        h::TransferPlot; sectors = nothing, transferkwargs = (;), thetaorigin = 0,
+        sector_formatter = string
+    )
+    below = length(h.args) == 1 ? h.args[1] : h.args[2]
+    kwargs = (; transferkwargs...)
+    if sectors !== nothing && get(kwargs, :howmany, 20) isa Int
+        howmany = Dict(c => get(kwargs, :howmany, 20) for c in sectors)
+        kwargs = (; kwargs..., howmany)
+    end
+    spectra = transfer_spectrum(h.args[1], below; kwargs...)
+
+    for (sector, spectrum) in pairs(spectra)
+        sectors === nothing || sector in sectors || continue
+
+        @series begin
+            yguide --> "r"
+            ylims --> (-Inf, 1.05)
+
+            xguide --> "θ"
+            xlims --> (thetaorigin, thetaorigin + 2pi)
+            xticks --> range(0, 2pi; length = 7)
+            xformatter --> x -> "$(rationalize(x / π, tol = 0.05))π"
+            xwiden --> true
+            seriestype := :scatter
+            markershape --> :auto
+            label := sector_formatter(sector)
+            return mod2pi.(angle.(spectrum) .+ thetaorigin) .- thetaorigin, abs.(spectrum)
+        end
+    end
+
+    legend --> false
+    grid --> :xy
+    framestyle --> :zerolines
+
+    return nothing
+end
+
+MPSKit.transferplot(args...; plotkwargs = (;), kwargs...) = transferplot(args...; kwargs..., plotkwargs...)
+MPSKit.transferplot!(args...; plotkwargs = (;), kwargs...) = transferplot!(args...; kwargs..., plotkwargs...)
+
+end
diff --git a/src/MPSKit.jl b/src/MPSKit.jl
index f900884ae..04eef7fa5 100644
--- a/src/MPSKit.jl
+++ b/src/MPSKit.jl
@@ -50,6 +50,7 @@ export correlation_length, marek_gap, transfer_spectrum
 export entropy, entanglement_spectrum
 export open_boundary_conditions, periodic_boundary_conditions
 export entanglementplot, transferplot
+export entanglementplot!, transferplot!
 export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # TODO: rename
 
 # unexported
@@ -76,7 +77,6 @@ using KrylovKit: KrylovAlgorithm
 using OptimKit
 using Base.Threads
 using Base.Iterators
-using RecipesBase
 using VectorInterface
 using Accessors
 using HalfIntegers
diff --git a/src/utility/plotting.jl b/src/utility/plotting.jl
index b58b912cc..6fc83f879 100644
--- a/src/utility/plotting.jl
+++ b/src/utility/plotting.jl
@@ -13,83 +13,21 @@ Plot the entanglement spectrum (see [`entanglement_spectrum`](@ref)) of a given
   between `site` and `site + 1`.
 - `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_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.
+- `plotkwargs = (; )`: kwargs for the underlying plot, e.g. `plotkwargs = (; title = "custom title", xlabel =   L"\text{custom label}", xticks = (1:2, ["a", "b"]))`.
 
 !!! note
-    You will need to manually import [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to
-    be able to use this function. MPSKit.jl defines its plots based on
-    [RecipesBase.jl](https://github.com/JuliaPlots/Plots.jl/tree/v2/RecipesBase), but the
-    user still has to add `using Plots` to be able to actually produce the plots.
+    You will need to manually import any plotting backend of [Makie.jl](https://github.com/MakieOrg/Makie.jl) or
+    [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to be able to use this function.
 
+See also [`entanglementplot!`](@ref) for plotting into an existing figure.
 """
 function entanglementplot end
-@userplot EntanglementPlot
-
-@recipe function f(
-        h::EntanglementPlot; site = 0, expand_symmetry = false, sortby = maximum,
-        sector_margin = 1 // 10, sector_formatter = string
-    )
-    mps = h.args[1]
-
-    spectra = entanglement_spectrum(mps, site)
-    sectors = []
-    spectrum = []
-    for (c, b) in pairs(spectra)
-        if expand_symmetry # Duplicate entries according to the quantum dimension.
-            b′ = repeat(b, dim(c))
-            sort!(b′; rev = true)
-            push!(spectrum, b′)
-        else
-            push!(spectrum, b)
-        end
-        push!(sectors, c)
-    end
-
-    if length(spectrum) > 1
-        order = sortperm(spectrum; by = sortby, rev = true)
-        spectrum = spectrum[order]
-        sectors = sectors[order]
-    end
-
-    for (i, (partial_spectrum, sector)) in enumerate(zip(spectrum, sectors))
-        @series begin
-            seriestype := :scatter
-            label := sector_formatter(sector)
-            n_spectrum = length(partial_spectrum)
-
-            # Put single dot in the middle, or a linear range with padding.
-            if n_spectrum == 1
-                x = [i + 1 // 2]
-            else
-                x = range(i + sector_margin, i + 1 - sector_margin; length = n_spectrum)
-            end
-            return x, partial_spectrum
-        end
-    end
-
-    title --> "Entanglement Spectrum"
-    legend --> false
-    grid --> :xy
-    widen --> true
-
-    xguide --> "χ = $(dim(_firstspace(mps.C[site])))"
-    xticks --> (1:length(sectors), sector_formatter.(sectors))
-    xtickfonthalign --> :center
-    xtick_direction --> :out
-    xrotation --> 45
-    xlims --> (1, length(sectors) + 1)
-
-    ylims --> (-Inf, 1 + 1.0e-1)
-    yscale --> :log10
-    label := nothing
-
-    return []
-end
+function entanglementplot! end
 
 """
-    transferplot(above, below = above; sectors = nothing, transferkwargs = (;)[, kwargs...])
+    transferplot(above, below = above; sectors = nothing, transferkwargs = (;), plotkwargs = (;), legend_position = :ct)
 
 Plot the partial transfer matrix spectrum of two InfiniteMPS's.
 
@@ -103,56 +41,17 @@ Plot the partial transfer matrix spectrum of two InfiniteMPS's.
 - `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.
+- `plotkwargs = (; )`: kwargs for the underlying plot, e.g. `plotkwargs = (; title = "custom title", xlabel = L"latexstring", xticks = (1:2, ["a", "b"]))`.
 - `thetaorigin = 0`: origin of the angle range.
 - `sector_formatter = string`: how to convert sectors to strings.
+- `legend_position = :ct`: Makie only, the `position` passed to `axislegend`.
+  For Plots, use the standard `legend` attribute instead (e.g. `legend = :topright`) in `plotkwargs`.
 
 !!! note
-    You will need to manually import [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to
-    be able to use this function. MPSKit.jl defines its plots based on
-    [RecipesBase.jl](https://github.com/JuliaPlots/Plots.jl/tree/v2/RecipesBase), but the
-    user still has to add `using Plots` to be able to actually produce the plots.
+    You will need to manually import any plotting backend of [Makie.jl](https://github.com/MakieOrg/Makie.jl) or
+    [Plots.jl](https://github.com/JuliaPlots/Plots.jl) to be able to use this function.
 
+See also [`transferplot!`](@ref) for plotting into an existing figure.
 """
 function transferplot end
-@userplot TransferPlot
-
-@recipe function f(
-        h::TransferPlot; sectors = nothing, transferkwargs = (;), thetaorigin = 0,
-        sector_formatter = string
-    )
-    below = length(h.args) == 1 ? h.args[1] : h.args[2]
-    kwargs = (; transferkwargs...)
-    if sectors !== nothing && get(kwargs, :howmany, 20) isa Int
-        # restrict the computation to the requested sectors
-        howmany = Dict(c => get(kwargs, :howmany, 20) for c in sectors)
-        kwargs = (; kwargs..., howmany)
-    end
-    spectra = transfer_spectrum(h.args[1], below; kwargs...)
-
-    for (sector, spectrum) in pairs(spectra)
-        sectors === nothing || sector in sectors || continue
-
-        @series begin
-            yguide --> "r"
-            ylims --> (-Inf, 1.05)
-
-            xguide --> "θ"
-            xlims --> (thetaorigin, thetaorigin + 2pi)
-            xticks --> range(0, 2pi; length = 7)
-            xformatter --> x -> "$(rationalize(x / π, tol = 0.05))π"
-            xwiden --> true
-            seriestype := :scatter
-            markershape --> :auto
-            label := sector_formatter(sector)
-            return mod2pi.(angle.(spectrum) .+ thetaorigin) .- thetaorigin, abs.(spectrum)
-        end
-    end
-
-    title --> "Transfer Spectrum"
-    legend --> false
-    grid --> :xy
-    framestyle --> :zerolines
-
-    return nothing
-end
+function transferplot! end
diff --git a/test/Project.toml b/test/Project.toml
index 5597b8c84..b65a77075 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,6 +8,8 @@ BlockTensorKit = "5f87ffc2-9cf1-4a46-8172-465d160bd8cd"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502"
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
@@ -27,8 +29,10 @@ MPSKit = {path = ".."}
 [compat]
 AMDGPU = "2"
 Aqua = "0.8.9"
+CairoMakie = "0.15"
 CUDA = "6"
 Combinatorics = "1"
+Makie = "0.24, 0.25"
 ParallelTestRunner = "2"
 Plots = "1.40"
 Pkg = "1"
diff --git a/test/internals/makie.jl b/test/internals/makie.jl
new file mode 100644
index 000000000..cdd39968e
--- /dev/null
+++ b/test/internals/makie.jl
@@ -0,0 +1,63 @@
+println("
+-----------------------------------
+|     Plot tests with Makie.jl    |
+-----------------------------------
+")
+
+using .TestSetup
+using Test, TestExtras
+using MPSKit
+using TensorKit
+using TensorKit: ℙ
+using CairoMakie
+
+@testset "plot tests" begin
+    ψ = InfiniteMPS([ℙ^2], [ℙ^5])
+    @test transferplot(ψ) isa Makie.FigureAxisPlot
+    @test transferplot(ψ, ψ) isa Makie.FigureAxisPlot
+    @test transferplot(ψ; transferkwargs = (; howmany = 3)) isa Makie.FigureAxisPlot
+
+    @test entanglementplot(ψ) isa Makie.FigureAxisPlot
+
+    # mutating plots
+    fig = Figure()
+    ax = Axis(fig[1, 1])
+    @test entanglementplot!(ax, ψ) isa Makie.Plot
+    @test transferplot!(Axis(fig[1, 2]), ψ) isa Makie.Plot
+    @test transferplot!(Axis(fig[1, 3]), ψ, ψ) isa Makie.Plot
+
+    # no target -> default to current axis
+    fig1 = Figure()
+    Axis(fig1[1, 1])
+    @test entanglementplot!(ψ) isa Makie.Plot
+    @test transferplot!(ψ) isa Makie.Plot
+    @test transferplot!(ψ, ψ) isa Makie.Plot
+
+    # detect plotkwargs in targeted axis
+    fig3 = Figure()
+    ax3 = Axis(fig3[1, 1])
+    entanglementplot!(ax3, ψ; plotkwargs = (; title = "custom"))
+    @test ax3.title[] == "custom"
+
+    # plotting into a non-current axis must not style the current one
+    fig2 = Figure()
+    target = Axis(fig2[1, 1])
+    current = Axis(fig2[1, 2]) # created last, so this is the current axis
+    entanglementplot!(target, ψ; plotkwargs = (; title = "custom"))
+    @test target.title[] == "custom"
+    @test current.title[] == ""
+end
+
+@testset "graded plots" begin
+    ψ = InfiniteMPS([Z2Space(0 => 1, 1 => 1)], [Z2Space(0 => 4, 1 => 4)])
+
+    @test entanglementplot(ψ) isa Makie.FigureAxisPlot
+    @test entanglementplot(ψ; site = 1) isa Makie.FigureAxisPlot
+
+    @test transferplot(ψ) isa Makie.FigureAxisPlot
+    @test transferplot(ψ, ψ) isa Makie.FigureAxisPlot
+
+    # restrict sectors
+    triv = unit(sectortype(ψ))
+    @test transferplot(ψ; sectors = [triv]) isa Makie.FigureAxisPlot
+end
diff --git a/test/internals/plots.jl b/test/internals/plots.jl
index 412d69d93..af5c23bbb 100644
--- a/test/internals/plots.jl
+++ b/test/internals/plots.jl
@@ -1,7 +1,7 @@
 println("
----------------------
-|   Plot tests       |
----------------------
+-----------------------------------
+|     Plot tests with Plots.jl    |
+-----------------------------------
 ")
 
 using .TestSetup
diff --git a/test/runtests.jl b/test/runtests.jl
index ef29c29ad..d79e4ae20 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -41,4 +41,10 @@ const init_code = quote
     const fast_tests = $fast
 end
 
-ParallelTestRunner.runtests(MPSKit, args; testsuite, init_worker_code, init_code)
+# custom workers for different plot backends
+function test_worker(name, init_worker_code)
+    name in ("misc/plots", "misc/makie") && return addworker(; init_worker_code)
+    return nothing
+end
+
+ParallelTestRunner.runtests(MPSKit, args; testsuite, init_worker_code, init_code, test_worker)