diff --git a/Project.toml b/Project.toml
index ea1bc9aa5..4475267c9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -43,6 +43,6 @@ StaticArrays = "1"
 StrideArrays = "0.1"
 TimerOutputs = "0.5"
 TrixiBase = "0.1"
-PointNeighbors = "0.2"
+PointNeighbors = "0.2.4"
 WriteVTK = "1"
 julia = "1.9"
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index f170e08d8..192870d9c 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -22,7 +22,8 @@ using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
 using TimerOutputs: TimerOutput, TimerOutputs, print_timer, reset_timer!
 using TrixiBase: trixi_include
-@reexport using PointNeighbors: TrivialNeighborhoodSearch, GridNeighborhoodSearch
+@reexport using PointNeighbors: TrivialNeighborhoodSearch, GridNeighborhoodSearch,
+                                NeighborListsNeighborhoodSearch
 using PointNeighbors: PointNeighbors, for_particle_neighbor
 using WriteVTK: vtk_grid, MeshCell, VTKCellTypes, paraview_collection, vtk_save
 
diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl
index 6a1d703b2..0e11b8c3f 100644
--- a/src/general/semidiscretization.jl
+++ b/src/general/semidiscretization.jl
@@ -148,6 +148,21 @@ function create_neighborhood_search(system, neighbor, ::Val{GridNeighborhoodSear
     return search
 end
 
+function create_neighborhood_search(system, neighbor,
+                                    ::Val{NeighborListsNeighborhoodSearch},
+                                    periodic_box_min_corner, periodic_box_max_corner,
+                                    threaded_nhs_update)
+    radius = compact_support(system, neighbor)
+    search = NeighborListsNeighborhoodSearch{ndims(system)}(radius, nparticles(neighbor),
+                                                            periodic_box_min_corner=periodic_box_min_corner,
+                                                            periodic_box_max_corner=periodic_box_max_corner)
+    # Initialize neighborhood search
+    PointNeighbors.initialize!(search, initial_coordinates(system),
+                               initial_coordinates(neighbor))
+
+    return search
+end
+
 @inline function compact_support(system, neighbor)
     (; smoothing_kernel, smoothing_length) = system
     return compact_support(smoothing_kernel, smoothing_length)