diff --git a/src/models/hamiltonians.jl b/src/models/hamiltonians.jl
index 7da3a25..92826dc 100644
--- a/src/models/hamiltonians.jl
+++ b/src/models/hamiltonians.jl
@@ -331,11 +331,14 @@ end
 
 MPO for the hamiltonian of the Bose-Hubbard model, as defined by
 ```math
-H = -t \\sum_{\\langle i,j \\rangle} \\left( a_{i}^+ a_{j}^- + a_{i}^- a_{j}^+ \\right) - \\sum_i \\mu N_i + \\frac{U}{2} \\sum_i N_i(N_i - 1).
+H = -t \\sum_{\\langle i,j \\rangle} \\left( a_{i}^+ a_{j}^- + a_{i}^- a_{j}^+ \\right) - \\mu \\sum_i N_i + \\frac{U}{2} \\sum_i N_i(N_i - 1).
 ```
 where ``N`` is the bosonic number operator [`a_number`](@ref).
 
-By default, the model is defined on an infinite chain with unit lattice spacing, without any symmetries and with `ComplexF64` entries of the tensors. The Hilbert space is truncated such that at maximum of `cutoff` bosons can be at a single site. If the `symmetry` is not `Trivial`, a fixed particle number density `n` can be imposed.
+By default, the model is defined on an infinite chain with unit lattice spacing, without any
+symmetries and with `ComplexF64` entries of the tensors. The Hilbert space is truncated such
+that at maximum of `cutoff` bosons can be at a single site. If the `symmetry` is not
+`Trivial`, a fixed (halfinteger) particle number density `n` can be imposed.
 """
 function bose_hubbard_model end
 function bose_hubbard_model(lattice::AbstractLattice; kwargs...)
@@ -348,7 +351,7 @@ end
 function bose_hubbard_model(elt::Type{<:Number}=ComplexF64,
                             symmetry::Type{<:Sector}=Trivial,
                             lattice::AbstractLattice=InfiniteChain(1);
-                            cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, n::Integer=0)
+                            cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, n=0)
     hopping_term = a_plusmin(elt, symmetry; cutoff=cutoff) +
                    a_minplus(elt, symmetry; cutoff=cutoff)
     N = a_number(elt, symmetry; cutoff=cutoff)
@@ -368,7 +371,7 @@ function bose_hubbard_model(elt::Type{<:Number}=ComplexF64,
     elseif symmetry === U1Irrep
         isinteger(2n) ||
             throw(ArgumentError("`U₁` symmetry requires halfinteger particle number"))
-        H = MPSKit.add_physical_charge(H, fill(U1Irrep(n), length(H)))
+        H = MPSKit.add_physical_charge(H, fill(U1Irrep(-n), length(H)))
     else
         throw(ArgumentError("symmetry not implemented"))
     end
diff --git a/test/bose_hubbard.jl b/test/bose_hubbard.jl
new file mode 100644
index 0000000..6616129
--- /dev/null
+++ b/test/bose_hubbard.jl
@@ -0,0 +1,40 @@
+using Test
+using TensorKit
+using MPSKit
+using MPSKitModels
+
+## Setup
+
+symmetry = U1Irrep
+cutoff = 3
+t = 1.0
+U = 10.0
+mu = 0.0
+n = 1
+lattice = InfiniteChain()
+
+Vspace = U1Space(0 => 8, 1 => 6, -1 => 6, 2 => 4, -2 => 4, 3 => 2, -3 => 2)
+
+alg = VUMPS(; maxiter=25, verbosity=0)
+
+# compare against higher-order analytic expansion from https://arxiv.org/pdf/1507.06426
+function exact_bose_hubbard_energy(; t=1.0, U=1.0)
+    J = t / U
+
+    E = 4 * U *
+        (-J^2 + J^4 + 68 / 9 * J^6 - 1267 / 81 * J^8 + 44171 / 1458 * J^10 -
+         4902596 / 6561 * J^12 -
+         8020902135607 / 2645395200 * J^14 - 32507578587517774813 / 466647713280000 * J^16)
+
+    return E
+end
+
+## Test
+
+@testset "Bose-Hubbard ground state" begin
+    H = bose_hubbard_model(symmetry, lattice; cutoff, t, U, mu, n)
+    ψ = InfiniteMPS([physicalspace(H, 1)], [Vspace])
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test expectation_value(ψ, H, envs) ≈ exact_bose_hubbard_energy(; t, U) atol = 1e-2
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 853d174..9fc5e38 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -37,6 +37,10 @@ end
     include("heisenberg.jl")
 end
 
+@testset "bose-hubbard model" begin
+    include("bose_hubbard.jl")
+end
+
 @testset "quantum potts model" begin
     include("quantum_potts.jl")
 end