diff --git a/docs/src/hamiltonians.md b/docs/src/hamiltonians.md
index 556197e6e..358fcd60d 100644
--- a/docs/src/hamiltonians.md
+++ b/docs/src/hamiltonians.md
@@ -19,17 +19,15 @@ of the Hubbard model in real and momentum space, as well as some other models.
 
 ### Real space Hubbard models
 ```@docs
+HubbardRealSpace
 HubbardReal1D
-BoseHubbardReal1D2C
 HubbardReal1DEP
-HubbardRealSpace
 ExtendedHubbardReal1D
 ```
 
 ### Momentum space Hubbard models
 ```@docs
 HubbardMom1D
-BoseHubbardMom1D2C
 HubbardMom1DEP
 ExtendedHubbardMom1D
 ```
@@ -174,3 +172,9 @@ Hamiltonians.log_abs_oscillator_zero
 ```@index
 Pages   = ["hamiltonians.md"]
 ```
+
+## Deprecated Hamiltonians
+```@docs
+BoseHubbardReal1D2C
+BoseHubbardMom1D2C
+```
diff --git a/src/Hamiltonians/BoseHubbardMom1D2C.jl b/src/Hamiltonians/BoseHubbardMom1D2C.jl
index 76c12c91e..85d6d9fe1 100644
--- a/src/Hamiltonians/BoseHubbardMom1D2C.jl
+++ b/src/Hamiltonians/BoseHubbardMom1D2C.jl
@@ -1,35 +1,48 @@
-"""
-    BoseHubbardMom1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, kwargs...)
+# """
+#     BoseHubbardMom1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, kwargs...)
+
+# Implements a one-dimensional Bose Hubbard chain in momentum space with a two-component
+# Bose gas.
 
-Implements a one-dimensional Bose Hubbard chain in momentum space with a two-component
-Bose gas.
+# ```math
+# \\hat{H} = \\hat{H}_a + \\hat{H}_b + \\frac{V}{M}\\sum_{kpqr} b^†_{r} a^†_{q} b_p a_k δ_{r+q,p+k}
+# ```
 
-```math
-\\hat{H} = \\hat{H}_a + \\hat{H}_b + \\frac{V}{M}\\sum_{kpqr} b^†_{r} a^†_{q} b_p a_k δ_{r+q,p+k}
-```
+# # Arguments
 
-# Arguments
+# * `address`: the starting address.
+# * `ua`: the `u` parameter for Hamiltonian a.
+# * `ub`: the `u` parameter for Hamiltonian b.
+# * `ta`: the `t` parameter for Hamiltonian a.
+# * `tb`: the `t` parameter for Hamiltonian b.
+# * `v`: the inter-species interaction parameter V.
+# Further keyword arguments are passed on to the constructor of [`HubbardMom1D`](@ref).
 
-* `address`: the starting address.
-* `ua`: the `u` parameter for Hamiltonian a.
-* `ub`: the `u` parameter for Hamiltonian b.
-* `ta`: the `t` parameter for Hamiltonian a.
-* `tb`: the `t` parameter for Hamiltonian b.
-* `v`: the inter-species interaction parameter V.
-Further keyword arguments are passed on to the constructor of [`HubbardMom1D`](@ref).
+# # See also
 
-# See also
+# * [`BoseFS2C`](@ref)
+# * [`BoseHubbardReal1D2C`](@ref)
 
-* [`BoseFS2C`](@ref)
-* [`BoseHubbardReal1D2C`](@ref)
+# """
+"""
+    BoseHubbardMom1D2C{T} <: TwoComponentHamiltonian{T}
+Special case `AbstractHamiltonian` for a 2 component Bose gas. Use
+[`HubbardMom1D(::BoseFS2C; kwargs...)`](@ref) to construct this type.
 
+!!! warning
+    This type is not part of the user interface and may be removed in the future.
 """
 struct BoseHubbardMom1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T}
     ha::HA
     hb::HB
 end
 
-function BoseHubbardMom1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, args...)
+function BoseHubbardMom1D2C(address; kwargs...)
+    @warn "Calling a the contructor to `BoseHubbardMom1D2C` is deprecated. Use `HubbardMom1D` instead."
+    return _BoseHubbardMom1D2C(address; kwargs...)
+end
+
+function _BoseHubbardMom1D2C(address::BoseFS2C; ua=1.0, ub=1.0, ta=1.0, tb=1.0, v=1.0, args...)
     ha = HubbardMom1D(address.bsa; u=ua, t=ta, args...)
     hb = HubbardMom1D(address.bsb; u=ub, t=tb, args...)
     T = promote_type(eltype(ha), eltype(hb))
@@ -43,7 +56,7 @@ function Base.show(io::IO, h::BoseHubbardMom1D2C)
     ta = h.ha.t
     tb = h.hb.t
     v = h.v
-    print(io, "BoseHubbardMom1D2C($addr; ua=$ua, ub=$ub, ta=$ta, tb=$tb, v=$v)")
+    print(io, "HubbardMom1D($addr; t=[$ta,$tb], u=[$ua $v; $v $ub])")
 end
 
 dimension(::BoseHubbardMom1D2C, address) = number_conserving_dimension(address)
diff --git a/src/Hamiltonians/BoseHubbardReal1D2C.jl b/src/Hamiltonians/BoseHubbardReal1D2C.jl
index 74411e6d4..63c0795b7 100644
--- a/src/Hamiltonians/BoseHubbardReal1D2C.jl
+++ b/src/Hamiltonians/BoseHubbardReal1D2C.jl
@@ -21,6 +21,9 @@ Implements a two-component one-dimensional Bose Hubbard chain in real space.
 * [`HubbardReal1D`](@ref)
 * [`BoseHubbardMom1D2C`](@ref)
 
+!!! warning
+    This type is deprecated and will be removed in a future release. Use
+    `HubbardRealSpace` instead.
 """
 struct BoseHubbardReal1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T}
     ha::HA
@@ -28,6 +31,8 @@ struct BoseHubbardReal1D2C{T,HA,HB,V} <: TwoComponentHamiltonian{T}
 end
 
 function BoseHubbardReal1D2C(address::BoseFS2C; ua=1.0,ub=1.0,ta=1.0,tb=1.0,v=1.0)
+    @warn "`BoseHubbardReal1D2C` is deprecated. Use `HubbardRealSpace` instead."
+
     ha = HubbardReal1D(address.bsa; u=ua, t=ta)
     hb = HubbardReal1D(address.bsb; u=ub, t=tb)
     T = promote_type(eltype(ha), eltype(hb))
diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl
index 3ec4ce128..49fbd1341 100644
--- a/src/Hamiltonians/Hamiltonians.jl
+++ b/src/Hamiltonians/Hamiltonians.jl
@@ -5,15 +5,13 @@ Hamiltonians.
 ## Exported concrete Hamiltonian types
 
 Real space Hubbard models
+ - [`HubbardRealSpace`](@ref)
  - [`HubbardReal1D`](@ref)
- - [`BoseHubbardReal1D2C`](@ref)
  - [`HubbardReal1DEP`](@ref)
- - [`HubbardRealSpace`](@ref)
  - [`ExtendedHubbardReal1D`](@ref)
 
 Momentum space Hubbard models
 - [`HubbardMom1D`](@ref)
-- [`BoseHubbardMom1D2C`](@ref)
 - [`HubbardMom1DEP`](@ref)
 
 Harmonic oscillator models
@@ -74,7 +72,6 @@ export MatrixHamiltonian
 export HubbardReal1D, HubbardMom1D, ExtendedHubbardReal1D, ExtendedHubbardMom1D, HubbardRealSpace
 export HubbardReal1DEP, shift_lattice, shift_lattice_inv
 export HubbardMom1DEP
-export BoseHubbardMom1D2C, BoseHubbardReal1D2C
 export GutzwillerSampling, GuidingVectorSampling
 export ParitySymmetry
 export TimeReversalSymmetry
@@ -94,6 +91,10 @@ export HOCartesianContactInteractions, HOCartesianEnergyConservedPerDim, HOCarte
 export AxialAngularMomentumHO
 export get_all_blocks, fock_to_cart
 
+# deprecated and will be removed
+export BoseHubbardMom1D2C
+export BoseHubbardReal1D2C
+
 include("abstract.jl")
 include("offdiagonals.jl")
 include("geometry.jl")
diff --git a/src/Hamiltonians/HubbardMom1D.jl b/src/Hamiltonians/HubbardMom1D.jl
index 56a66c1c9..6f78cc733 100644
--- a/src/Hamiltonians/HubbardMom1D.jl
+++ b/src/Hamiltonians/HubbardMom1D.jl
@@ -64,6 +64,52 @@ function HubbardMom1D(
     return HubbardMom1D{typeof(U),M,typeof(address),U,T}(address, ks, kes)
 end
 
+"""
+    HubbardMom1D(address::BoseFS2C; t=ones(C), u=ones(C, C), dispersion=hubbard_dispersion)
+
+Implements a one-dimensional Bose Hubbard chain in momentum space with a two-component
+Bose gas.
+
+```math
+\\hat{H} = \\hat{H}_a + \\hat{H}_b + \\frac{V}{M}\\sum_{kpqr} b^†_{r} a^†_{q} b_p a_k δ_{r+q,p+k}
+```
+
+# Arguments
+
+* `address`: the starting address.
+* `t`: the hopping strengths. Must be a vector of length `2`. The `i`-th element of the
+  vector corresponds to the hopping strength of the `i`-th component.
+* `u`: the on-site interaction parameters. Must be a symmetric 2×2 matrix. `u[i, j]`
+  corresponds to the interaction between the `i`-th and `j`-th component. `u[i, i]`
+  corresponds to the interaction of a component with itself.
+* `dispersion`: defines ``ϵ_k =``` t*dispersion(k)`
+    - [`hubbard_dispersion`](@ref): ``ϵ_k = -2t \\cos(k)``
+    - [`continuum_dispersion`](@ref): ``ϵ_k = tk^2``
+
+# See also
+
+* [`BoseFS2C`](@ref)
+"""
+function HubbardMom1D(address::BoseFS2C; t=ones(2), u=ones(2,2), args...)
+    # Sanity checks
+    if length(u) ≠ 1 && !issymmetric(u)
+        throw(ArgumentError("`u` must be symmetric"))
+    elseif length(u) ≠ 2 * 2
+        throw(ArgumentError("`u` must be a 2 × 2 matrix"))
+    elseif size(t) ≠ (2,)
+        throw(ArgumentError("`t` must be a vector of length 2"))
+    end
+
+    ta = t[1]; tb = t[2]
+    ua = u[1,1]; ub = u[2,2]; v = u[1,2]
+    ha = HubbardMom1D(address.bsa; u=ua, t=ta, args...)
+    hb = HubbardMom1D(address.bsb; u=ub, t=tb, args...)
+    T = promote_type(eltype(ha), eltype(hb), typeof(v))
+    V = T(v)
+    return BoseHubbardMom1D2C{T,typeof(ha),typeof(hb),V}(ha, hb)
+end
+
+
 function Base.show(io::IO, h::HubbardMom1D)
     compact_addr = repr(h.address, context=:compact => true) # compact print address
     print(io, "HubbardMom1D($(compact_addr); u=$(h.u), t=$(h.t))")
diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl
index 5c1dc195d..4a91ffd8d 100644
--- a/test/Hamiltonians.jl
+++ b/test/Hamiltonians.jl
@@ -202,10 +202,11 @@ end
                 FermiFS((1, 1, 1, 1, 1, 0, 0, 0)),
                 FermiFS((1, 1, 1, 1, 0, 0, 0, 0)),
             ); t=[1, 2], u=[0 3; 3 0]
-        ), BoseHubbardReal1D2C(BoseFS2C((1, 2, 3), (1, 0, 0))),
-        BoseHubbardMom1D2C(BoseFS2C((1, 2, 3), (1, 0, 0))),
+        ),
+        HubbardRealSpace(CompositeFS(BoseFS(1, 2, 3), BoseFS(1, 0, 0))),
+        HubbardMom1D(BoseFS2C((1, 2, 3), (1, 0, 0))),
         GutzwillerSampling(HubbardReal1D(BoseFS((1, 2, 3)); u=6); g=0.3),
-        GutzwillerSampling(BoseHubbardMom1D2C(BoseFS2C((3, 2, 1), (1, 2, 3)); ua=6); g=0.3),
+        GutzwillerSampling(HubbardMom1D(BoseFS2C((3, 2, 1), (1, 2, 3)); u=[6 1; 1 1]); g=0.3),
         GutzwillerSampling(HubbardReal1D(BoseFS((1, 2, 3)); u=6 + 2im); g=0.3),
         MatrixHamiltonian(Float64[1 2; 2 0]),
         GutzwillerSampling(MatrixHamiltonian([1.0 2.0; 2.0 0.0]); g=0.3),
@@ -218,7 +219,7 @@ end
         HubbardMom1DEP(CompositeFS(FermiFS((0, 1, 1, 0, 0)), FermiFS((0, 0, 1, 0, 0))), v_ho=5),
         ParitySymmetry(HubbardRealSpace(CompositeFS(BoseFS((1, 2, 0)), FermiFS((0, 1, 0))))),
         TimeReversalSymmetry(HubbardMom1D(FermiFS2C((1, 0, 1), (0, 1, 1)))),
-        TimeReversalSymmetry(BoseHubbardMom1D2C(BoseFS2C((0, 1, 1), (1, 0, 1)))),
+        TimeReversalSymmetry(HubbardMom1D(BoseFS2C((0, 1, 1), (1, 0, 1)))),
         Stoquastic(HubbardMom1D(BoseFS((0, 5, 0)))),
         momentum(HubbardMom1D(BoseFS((0, 5, 0)))),
         HOCartesianContactInteractions(BoseFS((2, 0, 0, 0))),
@@ -427,7 +428,7 @@ end
     addr1 = near_uniform(BoseFS2C{1,100,20})
     addr2 = near_uniform(BoseFS2C{100,1,20})
 
-    for Hamiltonian in (BoseHubbardReal1D2C, BoseHubbardMom1D2C)
+    for Hamiltonian in (BoseHubbardReal1D2C,)
         @testset "$Hamiltonian" begin
             H1 = BoseHubbardReal1D2C(addr1; ta=1.0, tb=2.0, ua=0.5, ub=0.7, v=0.2)
             H2 = BoseHubbardReal1D2C(addr2; ta=2.0, tb=1.0, ua=0.7, ub=0.5, v=0.2)
@@ -694,8 +695,8 @@ end
             for H in (
                 HubbardMom1D(BoseFS((2,2,2)), u=6),
                 ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
-                ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
-                BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0),
+                ExtendedHubbardMom1D(BoseFS((1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), u=6, t=2.0),
+                HubbardMom1D(BoseFS2C((1,2,3), (1,0,0)), u=[1 1; 1 2]),
             )
                 # GutzwillerSampling with parameter zero is exactly equal to the original H
                 G = GutzwillerSampling(H, 0.0)
@@ -726,7 +727,7 @@ end
                 HubbardMom1D(BoseFS((2,2,2)), u=6),
                 ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
                 ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
-                # BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), # multicomponent not implemented for G2RealCorrelator
+                # multicomponent not implemented for G2RealCorrelator
             )
                 # energy
                 g = rand()
@@ -810,7 +811,7 @@ end
                 HubbardMom1D(BoseFS((2,2,2)), u=6),
                 ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
                 ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
-                # BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0), # multicomponent not implemented for G2RealCorrelator
+                # multicomponent not implemented for G2RealCorrelator
             )
                 # energy
                 x = rand()
@@ -879,8 +880,8 @@ end
         for H in (
             HubbardMom1D(BoseFS((2,2,2)), u=6),
             ExtendedHubbardReal1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
-            ExtendedHubbardMom1D(BoseFS((1,1,1,1,1,1,1,1,1,1,1,1)), u=6, t=2.0),
-            BoseHubbardMom1D2C(BoseFS2C((1,2,3), (1,0,0)), ub=2.0),
+            ExtendedHubbardMom1D(BoseFS((1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), u=6, t=2.0),
+            HubbardMom1D(BoseFS2C((1,2,3), (1,0,0))),
         )
             @test_throws ArgumentError Rimu.Hamiltonians.TransformUndoer(H)
             @test_throws ArgumentError Rimu.Hamiltonians.TransformUndoer(H, H)
@@ -1094,7 +1095,7 @@ using Rimu.Hamiltonians: circshift_dot
     end
 
     @testset "G2MomCorrelator" begin
-        # v0 is the exact ground state from BoseHubbardMom1D2C(aIni;ua=0,ub=0,v=0.1)
+        # v0 is the exact ground state from HubbardMom1D(aIni; u=[0 0.1; 0.1 0])
         bfs1 = BoseFS([0, 2, 0])
         bfs2 = BoseFS([0, 1, 0])
         aIni = BoseFS2C(bfs1,bfs2)
@@ -1571,7 +1572,7 @@ end
 
 @testset "TimeReversalSymmetry" begin
     @test_throws ArgumentError TimeReversalSymmetry(HubbardMom1D(BoseFS((1, 1))))
-    @test_throws ArgumentError TimeReversalSymmetry(BoseHubbardMom1D2C(BoseFS2C((1, 1),(2,1))))
+    @test_throws ArgumentError TimeReversalSymmetry(HubbardMom1D(BoseFS2C((1, 1),(2,1))))
     @test_throws ArgumentError begin
         TimeReversalSymmetry(HubbardRealSpace(CompositeFS(FermiFS((1, 1)),BoseFS((2,1)))))
     end
@@ -1591,8 +1592,8 @@ end
         @test issymmetric(even_m)
         @test issymmetric(odd_m)
     end
-    @testset "2-particle BoseHubbardMom1D2C" begin
-        ham = BoseHubbardMom1D2C(BoseFS2C((0,1,1),(1,0,1)))
+    @testset "2-particle 2-component HubbardMom1D" begin
+        ham = HubbardMom1D(BoseFS2C((0,1,1),(1,0,1)))
         even = TimeReversalSymmetry(ham)
         odd = TimeReversalSymmetry(ham; even=false)