diff --git a/NEWS.md b/NEWS.md
index ad3e44c715..81fc27af62 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -18,6 +18,11 @@ used in the Julia ecosystem. Notable changes will be documented in this file for
 - Added the number of split integration time steps to the `InfoCallback` output
   when a `SplitIntegrationCallback` is used (#1194).
 
+### Performance
+
+- Improved multithreaded TLSPH performance on NUMA systems by initializing runtime
+  arrays in parallel (#1294).
+
 ### Important Bugfixes
 
 - Fixed signed-distance constraints in `ParticlePackingSystem` when using a separate
diff --git a/src/schemes/structure/total_lagrangian_sph/system.jl b/src/schemes/structure/total_lagrangian_sph/system.jl
index a1cd48fb54..339306a1e5 100644
--- a/src/schemes/structure/total_lagrangian_sph/system.jl
+++ b/src/schemes/structure/total_lagrangian_sph/system.jl
@@ -145,10 +145,26 @@ function TotalLagrangianSPHSystem(initial_condition; smoothing_kernel, smoothing
         poisson_ratio_sorted = poisson_ratio
     end
 
-    initial_coordinates = copy(initial_condition_sorted.coordinates)
-    current_coordinates = copy(initial_condition_sorted.coordinates)
-    mass = copy(initial_condition_sorted.mass)
-    material_density = copy(initial_condition_sorted.density)
+    initial_coordinates = similar(initial_condition_sorted.coordinates)
+    current_coordinates = similar(initial_condition_sorted.coordinates)
+    mass = similar(initial_condition_sorted.mass)
+    material_density = similar(initial_condition_sorted.density)
+
+    # Initialize the runtime arrays in parallel so that, on NUMA systems, the first-touch
+    # allocation policy places their memory close to the threads that will access it
+    # during the simulation.
+    # The semidiscretization backend is not available in the constructor yet, and using
+    # Polyester here is harmless for serial and GPU simulations.
+    # This makes the RHS significantly faster on large data center CPUs with multiple
+    # NUMA domains (see http://github.com/trixi-framework/TrixiParticles.jl/pull/1294).
+    parallelization_backend = PolyesterBackend()
+    copyto_threaded!(initial_coordinates, initial_condition_sorted.coordinates,
+                     parallelization_backend)
+    copyto_threaded!(current_coordinates, initial_condition_sorted.coordinates,
+                     parallelization_backend)
+    copyto_threaded!(mass, initial_condition_sorted.mass, parallelization_backend)
+    copyto_threaded!(material_density, initial_condition_sorted.density,
+                     parallelization_backend)
     correction_matrix = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles)
     pk1_rho2 = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles)
     deformation_grad = Array{ELTYPE, 3}(undef, NDIMS, NDIMS, n_particles)
diff --git a/src/util.jl b/src/util.jl
index 650a2d0e88..6ca7c9ec6a 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -301,3 +301,27 @@ function Base.similar(::Broadcast.Broadcasted{ThreadedBroadcastStyle{P}},
     # TODO we only have the type `P` here and just assume that we can do `P()`
     return ThreadedBroadcastArray(similar(Array{T}, dims), parallelization_backend=P())
 end
+
+# Like `copyto!`, but parallelize over particles to use the same memory access pattern
+# as the particle loops in the right-hand side.
+@inline function copyto_threaded!(dest::Matrix, src, parallelization_backend)
+    if axes(dest) != axes(src)
+        throw(DimensionMismatch("source and destination must have the same axes"))
+    end
+
+    @threaded parallelization_backend for particle in axes(dest, 2)
+        for dimension in axes(dest, 1)
+            @inbounds dest[dimension, particle] = src[dimension, particle]
+        end
+    end
+
+    return dest
+end
+
+@inline function copyto_threaded!(dest::Vector, src, parallelization_backend)
+    @threaded parallelization_backend for particle in eachindex(dest, src)
+        @inbounds dest[particle] = src[particle]
+    end
+
+    return dest
+end
diff --git a/test/general/util.jl b/test/general/util.jl
index 6f96607012..701b1314ea 100644
--- a/test/general/util.jl
+++ b/test/general/util.jl
@@ -14,6 +14,22 @@
     @test visited == [(1, 3), (2, 4)]
 end
 
+@testset verbose=true "copyto_threaded!" begin
+    backend = SerialBackend()
+
+    matrix = zeros(2, 3)
+    TrixiParticles.copyto_threaded!(matrix, reshape(1:6, 2, 3), backend)
+    @test matrix == reshape(1:6, 2, 3)
+    @test_throws DimensionMismatch TrixiParticles.copyto_threaded!(matrix, zeros(2, 2),
+                                                                   backend)
+
+    vector = zeros(3)
+    TrixiParticles.copyto_threaded!(vector, 1:3, backend)
+    @test vector == 1:3
+    @test_throws DimensionMismatch TrixiParticles.copyto_threaded!(vector, zeros(2),
+                                                                   backend)
+end
+
 @testset verbose=true "ThreadedBroadcastArray" begin
     A = TrixiParticles.ThreadedBroadcastArray(ones(3, 3))
     B = ones(3, 3)