diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index b5b9c8ff..e7d279f4 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -67,12 +67,14 @@ jobs:
       - uses: julia-actions/cache@v2
       - name: Build package
         uses: julia-actions/julia-buildpkg@v1
-      - name: Run tests
+      - name: Run unit tests
         uses: julia-actions/julia-runtest@v1
         with:
           annotate: true
           # Only run coverage in one Job (Ubuntu and latest Julia version)
           coverage: ${{ matrix.os == 'ubuntu-latest' && matrix.version == '1' }}
+        env:
+          POINTNEIGHBORS_TEST: unit
       - name: Process coverage results
         # Only run coverage in one Job (Ubuntu and latest Julia version)
         if: matrix.os == 'ubuntu-latest' && matrix.version == '1'
@@ -89,3 +91,11 @@ jobs:
           flags: unit
         env:
           CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+      - name: Run benchmark tests
+        uses: julia-actions/julia-runtest@v1
+        with:
+          annotate: true
+          coverage: false
+        env:
+          POINTNEIGHBORS_TEST: benchmarks
+
diff --git a/benchmarks/benchmarks.jl b/benchmarks/benchmarks.jl
new file mode 100644
index 00000000..59e0bd2d
--- /dev/null
+++ b/benchmarks/benchmarks.jl
@@ -0,0 +1,4 @@
+include("count_neighbors.jl")
+include("n_body.jl")
+
+include("plot.jl")
diff --git a/benchmarks/count_neighbors.jl b/benchmarks/count_neighbors.jl
new file mode 100644
index 00000000..9f73add7
--- /dev/null
+++ b/benchmarks/count_neighbors.jl
@@ -0,0 +1,29 @@
+using PointNeighbors
+using BenchmarkTools
+
+"""
+    benchmark_count_neighbors(neighborhood_search, coordinates; parallel = true)
+
+A very cheap and simple neighborhood search benchmark, only counting the neighbors of each
+point. For each point-neighbor pair, only an array entry is incremented.
+
+Due to the minimal computational cost, differences between neighborhood search
+implementations are highlighted. On the other hand, this is the least realistic benchmark.
+
+For a computationally heavier benchmark, see [`benchmark_n_body`](@ref).
+"""
+function benchmark_count_neighbors(neighborhood_search, coordinates; parallel = true)
+    n_neighbors = zeros(Int, size(coordinates, 2))
+
+    function count_neighbors!(n_neighbors, coordinates, neighborhood_search, parallel)
+        n_neighbors .= 0
+
+        foreach_point_neighbor(coordinates, coordinates, neighborhood_search,
+                               parallel = parallel) do i, _, _, _
+            n_neighbors[i] += 1
+        end
+    end
+
+    return @belapsed $count_neighbors!($n_neighbors, $coordinates,
+                                       $neighborhood_search, $parallel)
+end
diff --git a/benchmarks/n_body.jl b/benchmarks/n_body.jl
new file mode 100644
index 00000000..5b04a266
--- /dev/null
+++ b/benchmarks/n_body.jl
@@ -0,0 +1,41 @@
+using PointNeighbors
+using BenchmarkTools
+
+"""
+    benchmark_n_body(neighborhood_search, coordinates; parallel = true)
+
+A simple neighborhood search benchmark, computing the right-hand side of an n-body
+simulation with a cutoff (corresponding to the search radius of `neighborhood_search`).
+
+This is a more realistic benchmark for particle-based simulations than
+[`benchmark_count_neighbors`](@ref).
+However, due to the higher computational cost, differences between neighborhood search
+implementations are less pronounced.
+"""
+function benchmark_n_body(neighborhood_search, coordinates; parallel = true)
+    mass = 1e10 * (rand(size(coordinates, 2)) .+ 1)
+    G = 6.6743e-11
+
+    dv = similar(coordinates)
+
+    function compute_acceleration!(dv, coordinates, mass, G, neighborhood_search, parallel)
+        dv .= 0.0
+
+        foreach_point_neighbor(coordinates, coordinates, neighborhood_search,
+                               parallel = parallel) do i, j, pos_diff, distance
+            # Only consider particles with a distance > 0
+            distance < sqrt(eps()) && return
+
+            dv_ = -G * mass[j] * pos_diff / distance^3
+
+            for dim in axes(dv, 1)
+                @inbounds dv[dim, i] += dv_[dim]
+            end
+        end
+
+        return dv
+    end
+
+    return @belapsed $compute_acceleration!($dv, $coordinates, $mass, $G,
+                                            $neighborhood_search, $parallel)
+end
diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl
new file mode 100644
index 00000000..28f0beb7
--- /dev/null
+++ b/benchmarks/plot.jl
@@ -0,0 +1,80 @@
+using Plots
+using BenchmarkTools
+
+# Generate a rectangular point cloud
+include("../test/point_cloud.jl")
+
+"""
+    plot_benchmarks(benchmark, n_points_per_dimension, iterations;
+                    seed = 1, perturbation_factor_position = 1.0,
+                    parallel = true, title = "")
+
+Run a benchmark for with several neighborhood searches multiple times for increasing numbers
+of points and plot the results.
+
+# Arguments
+- `benchmark`:              The benchmark function. See [`benchmark_count_neighbors`](@ref)
+                            and [`benchmark_n_body`](@ref).
+- `n_points_per_dimension`: Initial resolution as tuple. The product is the initial number
+                            of points. For example, use `(100, 100)` for a 2D benchmark or
+                            `(10, 10, 10)` for a 3D benchmark.
+- `iterations`:             Number of refinement iterations
+
+# Keywords
+- `parallel = true`:        Loop over all points in parallel
+- `title = ""`:             Title of the plot
+- `seed = 1`:               Seed to perturb the point positions. Different seeds yield
+                            slightly different point positions.
+- `perturbation_factor_position = 1.0`: Perturb point positions by this factor. A factor of
+                                        `1.0` corresponds to points being moved by
+                                        a maximum distance of `0.5` along each axis.
+"""
+function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
+                         parallel = true, title = "",
+                         seed = 1, perturbation_factor_position = 1.0)
+    neighborhood_searches_names = ["TrivialNeighborhoodSearch";;
+                                   "GridNeighborhoodSearch";;
+                                   "PrecomputedNeighborhoodSearch"]
+
+    # Multiply number of points in each iteration (roughly) by this factor
+    scaling_factor = 4
+    per_dimension_factor = scaling_factor^(1 / length(n_points_per_dimension))
+    sizes = [round.(Int, n_points_per_dimension .* per_dimension_factor^(iter - 1))
+             for iter in 1:iterations]
+
+    n_particles_vec = prod.(sizes)
+    times = zeros(iterations, length(neighborhood_searches_names))
+
+    for iter in 1:iterations
+        coordinates = point_cloud(sizes[iter], seed = seed,
+                                  perturbation_factor_position = perturbation_factor_position)
+
+        search_radius = 3.0
+        NDIMS = size(coordinates, 1)
+        n_particles = size(coordinates, 2)
+
+        neighborhood_searches = [
+            TrivialNeighborhoodSearch{NDIMS}(; search_radius, eachpoint = 1:n_particles),
+            GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles),
+            PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles),
+        ]
+
+        for i in eachindex(neighborhood_searches)
+            neighborhood_search = neighborhood_searches[i]
+            initialize!(neighborhood_search, coordinates, coordinates)
+
+            time = benchmark(neighborhood_search, coordinates, parallel = parallel)
+            times[iter, i] = time
+            time_string = BenchmarkTools.prettytime(time * 1e9)
+            println("$(neighborhood_searches_names[i])")
+            println("with $(join(sizes[iter], "x")) = $(prod(sizes[iter])) particles finished in $time_string\n")
+        end
+    end
+
+    plot(n_particles_vec, times,
+         xaxis = :log, yaxis = :log,
+         xticks = (n_particles_vec, n_particles_vec),
+         xlabel = "#particles", ylabel = "Runtime [s]",
+         legend = :outerright, size = (750, 400), dpi = 200,
+         label = neighborhood_searches_names, title = title)
+end
diff --git a/test/Project.toml b/test/Project.toml
index c3ee0727..a44e84ee 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,6 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
diff --git a/test/benchmarks.jl b/test/benchmarks.jl
new file mode 100644
index 00000000..15469a4d
--- /dev/null
+++ b/test/benchmarks.jl
@@ -0,0 +1,17 @@
+# Check that all benchmarks are running without errors.
+# Note that these are only smoke tests, not verifying the result.
+# Also note that these tests are run without coverage checks, since we want to
+# cover everything with unit tests.
+@testset verbose=true "Benchmarks" begin
+    include("../benchmarks/benchmarks.jl")
+
+    @testset verbose=true "$(length(size))D" for size in [(50,), (10, 10), (5, 5, 5)]
+        @testset verbose=true "`benchmark_count_neighbors`" begin
+            @test_nowarn_mod plot_benchmarks(benchmark_count_neighbors, size, 2)
+        end
+
+        @testset verbose=true "`benchmark_n_body`" begin
+            @test_nowarn_mod plot_benchmarks(benchmark_n_body, size, 2)
+        end
+    end
+end;
diff --git a/test/runtests.jl b/test/runtests.jl
index 0eb0f4b5..cc1f9f9b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,13 @@
 include("test_util.jl")
 
+const POINTNEIGHBORS_TEST = lowercase(get(ENV, "POINTNEIGHBORS_TEST", "all"))
+
 @testset verbose=true "PointNeighbors.jl Tests" begin
-    include("nhs_trivial.jl")
-    include("nhs_grid.jl")
-    include("neighborhood_search.jl")
+    if POINTNEIGHBORS_TEST in ("all", "unit")
+        include("unittest.jl")
+    end
+
+    if POINTNEIGHBORS_TEST in ("all", "benchmarks")
+        include("benchmarks.jl")
+    end
 end;
diff --git a/test/unittest.jl b/test/unittest.jl
new file mode 100644
index 00000000..983eafe2
--- /dev/null
+++ b/test/unittest.jl
@@ -0,0 +1,7 @@
+# Separate file that can be executed to only run unit tests.
+# Include `test_util.jl` first.
+@testset verbose=true "Unit Tests" begin
+    include("nhs_trivial.jl")
+    include("nhs_grid.jl")
+    include("neighborhood_search.jl")
+end;