diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml
index 5f0fe4967..329ff810d 100644
--- a/.github/workflows/actions.yml
+++ b/.github/workflows/actions.yml
@@ -21,7 +21,7 @@ jobs:
           version: ${{ matrix.julia-version }}
           arch: ${{ matrix.julia-arch }}
       - name: "Load cache"
-        uses: julia-actions/cache@v1
+        uses: julia-actions/cache@v2
       - name: "Build"
         uses: julia-actions/julia-buildpkg@v1
       - name: "Run tests"
diff --git a/.github/workflows/allocations.yml b/.github/workflows/allocations.yml
index a74ea8bd2..4b576dab9 100644
--- a/.github/workflows/allocations.yml
+++ b/.github/workflows/allocations.yml
@@ -21,7 +21,7 @@ jobs:
         with:
           version: ${{ matrix.julia-version }}
           arch: ${{ matrix.julia-arch }}
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v2
       - name: "Run test"
         run: |
           julia --project=@. -e "using Pkg; Pkg.instantiate(); Pkg.build()"
diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
index b315723a4..10756b75e 100644
--- a/.github/workflows/benchmark.yml
+++ b/.github/workflows/benchmark.yml
@@ -13,7 +13,7 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: 1
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - name: Install dependencies
         run: julia -e 'using Pkg; pkg"add PkgBenchmark BenchmarkTools BenchmarkCI@0.1"'
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index f3a4442c6..fc943f72e 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -15,8 +15,8 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: '1.7'
-      - uses: julia-actions/cache@v1
+          version: '1'
+      - uses: julia-actions/cache@v2
       - name: Install dependencies
         run: |
           export JULIA_PROJECT=@.
diff --git a/docs/src/hamiltonians.md b/docs/src/hamiltonians.md
index e463c492d..850a8e61b 100644
--- a/docs/src/hamiltonians.md
+++ b/docs/src/hamiltonians.md
@@ -102,6 +102,7 @@ Observables are [`AbstractHamiltonian`](@ref)s that represent a physical
 observable. Their ground state expectation values can be sampled by passing
 them into [`AllOverlaps`](@ref).
 ```@docs
+ParticleNumberOperator
 G2RealCorrelator
 G2RealSpace
 G2MomCorrelator
diff --git a/src/Hamiltonians/HOCartesianContactInteractions.jl b/src/Hamiltonians/HOCartesianContactInteractions.jl
index 1025f742f..bb7d3d47e 100644
--- a/src/Hamiltonians/HOCartesianContactInteractions.jl
+++ b/src/Hamiltonians/HOCartesianContactInteractions.jl
@@ -31,7 +31,16 @@ function four_oscillator_integral_general(i, j, k, l; max_level = typemax(Int))
     p = sqrt(2 * gamma(i + 1) * gamma(j + 1) * gamma(k + 1) * gamma(l + 1)) * pi^2
     q = gamma(a/2) * gamma(b/2) * gamma(c/2) * gamma(k + 1) / gamma(k - l + 1)
     
-    f = _₃F₂(float(-l), b/2, c/2, float(1 + k - l), d/2, 1.0)
+    f1 = _₃F₂(float(-l), b/2, c/2, float(1 + k - l), d/2, 1.0)
+
+    if isnan(f1)
+        # workaround for some issues with large arguments
+        fp = _₃F₂(float(-l), b/2, c/2, float(1 + k - l), d/2, 1.0 + eps())
+        fm = _₃F₂(float(-l), b/2, c/2, float(1 + k - l), d/2, 1.0 - eps())
+        f = (fp + fm)/2
+    else
+        f = f1
+    end
     
     return f * q / p
 end
diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl
index 1182183ec..eff6b9d06 100644
--- a/src/Hamiltonians/Hamiltonians.jl
+++ b/src/Hamiltonians/Hamiltonians.jl
@@ -36,6 +36,7 @@ Other
 - [`Stoquastic`](@ref)
 
 ## [Observables](#Observables)
+- [`ParticleNumberOperator`](@ref)
 - [`G2MomCorrelator`](@ref)
 - [`G2RealCorrelator`](@ref)
 - [`DensityMatrixDiagonal`](@ref)
@@ -83,6 +84,7 @@ export Stoquastic
 export Transcorrelated1D
 export hubbard_dispersion, continuum_dispersion
 export FroehlichPolaron
+export ParticleNumberOperator
 
 export G2MomCorrelator, G2RealCorrelator, G2RealSpace, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum
 export StringCorrelator
@@ -124,6 +126,7 @@ include("Transcorrelated1D.jl")
 include("correlation_functions.jl")
 include("DensityMatrixDiagonal.jl")
 include("Momentum.jl")
+include("particle_number.jl")
 
 include("HOCartesianContactInteractions.jl")
 include("HOCartesianEnergyConservedPerDim.jl")
diff --git a/src/Hamiltonians/Momentum.jl b/src/Hamiltonians/Momentum.jl
index 63335a95f..c5f4fe414 100644
--- a/src/Hamiltonians/Momentum.jl
+++ b/src/Hamiltonians/Momentum.jl
@@ -27,7 +27,6 @@ Momentum(component=0; fold=true) = Momentum{component}(fold)
 Base.show(io::IO, mom::Momentum{C}) where {C} = print(io, "Momentum($C; fold=$(mom.fold))")
 
 LOStructure(::Type{<:Momentum}) = IsDiagonal()
-num_offdiagonals(ham::Momentum, add) = 0
 
 @inline function _momentum(add::SingleComponentFockAddress, fold)
     M = num_modes(add)
diff --git a/src/Hamiltonians/particle_number.jl b/src/Hamiltonians/particle_number.jl
new file mode 100644
index 000000000..c41788a7a
--- /dev/null
+++ b/src/Hamiltonians/particle_number.jl
@@ -0,0 +1,36 @@
+"""
+    ParticleNumberOperator([address]) <: AbstractHamiltonian
+
+The number operator in Fock space. This operator is diagonal in the Fock basis and
+returns the number of particles in the Fock state. Passing an address is optional.
+
+```jldoctest
+julia> h = FroehlichPolaron(fs"|0 0⟩{}"; mode_cutoff=5, v=3); bsr = BasisSetRep(h);
+
+julia> gs = DVec(zip(bsr.basis, eigen(Matrix(bsr)).vectors[:,1])); # ground state
+
+julia> dot(gs, ParticleNumberOperator(), gs) # particle number expectation value
+2.8823297252925917
+```
+
+See also [`AbstractHamiltonian`](@ref).
+"""
+struct ParticleNumberOperator{A} <: AbstractHamiltonian{Float64}
+    address::A
+end
+ParticleNumberOperator() = ParticleNumberOperator(BoseFS(1,))
+
+function Base.show(io::IO, n::ParticleNumberOperator)
+    io = IOContext(io, :compact => true)
+    print(io, "ParticleNumberOperator(")
+    n.address === BoseFS(1,) || show(io, n.address) # suppress if default
+    print(io, ")")
+end
+
+LOStructure(::Type{<:ParticleNumberOperator}) = IsDiagonal()
+starting_address(n::ParticleNumberOperator) = n.address
+
+function diagonal_element(::ParticleNumberOperator, addr::AbstractFockAddress)
+    return float(num_particles(addr))
+end
+allowed_address_type(::ParticleNumberOperator) = AbstractFockAddress
diff --git a/src/Interfaces/hamiltonians.jl b/src/Interfaces/hamiltonians.jl
index 553e179e4..25eb52db5 100644
--- a/src/Interfaces/hamiltonians.jl
+++ b/src/Interfaces/hamiltonians.jl
@@ -252,6 +252,12 @@ LOStructure(op) = LOStructure(typeof(op))
 LOStructure(::Type) = AdjointUnknown()
 LOStructure(::AbstractMatrix) = AdjointKnown()
 
+# diagonal matrices have zero offdiagonal elements
+function num_offdiagonals(h::H, addr) where {H<:AbstractHamiltonian}
+    return num_offdiagonals(LOStructure(H), h, addr)
+end
+num_offdiagonals(::IsDiagonal, _, _) = 0
+
 """
     has_adjoint(op)
 
diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl
index c24c18f94..837bcad22 100644
--- a/test/Hamiltonians.jl
+++ b/test/Hamiltonians.jl
@@ -196,7 +196,9 @@ end
         HOCartesianCentralImpurity(BoseFS((1,0,0,0,0))),
 
         FroehlichPolaron(OccupationNumberFS(1,1,1)),
-        FroehlichPolaron(OccupationNumberFS(1,1,1); momentum_cutoff = 10.0)
+        FroehlichPolaron(OccupationNumberFS(1,1,1); momentum_cutoff = 10.0),
+
+        ParticleNumberOperator(OccupationNumberFS(1, 1, 1))
     )
         test_hamiltonian_interface(H)
     end