diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index 722b9fe3de4..3df4a6752a9 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -1980,7 +1980,7 @@ namespace aspect
       /**
        * The world holding the particles
        */
-      std::unique_ptr<Particle::World<dim>> particle_world;
+      std::vector<std::unique_ptr<Particle::World<dim>>> particle_worlds;
 
       /**
        * @}
diff --git a/include/aspect/simulator_access.h b/include/aspect/simulator_access.h
index af167715830..a65c0cf4879 100644
--- a/include/aspect/simulator_access.h
+++ b/include/aspect/simulator_access.h
@@ -937,20 +937,20 @@ namespace aspect
       n_particle_worlds() const;
 
       /**
-       * Returns a const reference to the particle world, in case anyone
+       * Returns a const reference to a single particle world, in case anyone
        * wants to query something about particles.
        */
       const Particle::World<dim> &
-      get_particle_world() const;
+      get_particle_world(const unsigned int particle_world_index) const;
 
       /**
-       * Returns a reference to the particle world, in case anyone wants to
+       * Returns a reference to a single particle world, in case anyone wants to
        * change something within the particle world. Use with care, usually
        * you want to only let the functions within the particle subsystem
        * change member variables of the particle world.
        */
       Particle::World<dim> &
-      get_particle_world();
+      get_particle_world(const unsigned int particle_world_index);
 
       /**
        * Return true if using the block GMG Stokes solver.
diff --git a/source/material_model/rheology/strain_dependent.cc b/source/material_model/rheology/strain_dependent.cc
index 70eae8ceb39..d50c760d20c 100644
--- a/source/material_model/rheology/strain_dependent.cc
+++ b/source/material_model/rheology/strain_dependent.cc
@@ -412,7 +412,7 @@ namespace aspect
         // Currently this functionality only works in field composition
         if (healing_mechanism != no_healing && this->n_particle_worlds() > 0)
           {
-            const Particle::Property::Manager<dim> &particle_property_manager = this->get_particle_world().get_property_manager();
+            const Particle::Property::Manager<dim> &particle_property_manager = this->get_particle_world(0).get_property_manager();
             AssertThrow(particle_property_manager.plugin_name_exists("viscoplastic strain invariants") == false, ExcMessage("This healing mechanism currently does not work if the strain is tracked on particles."));
           }
 
diff --git a/source/mesh_refinement/particle_density.cc b/source/mesh_refinement/particle_density.cc
index ddb125d1f80..9113dcd71e4 100644
--- a/source/mesh_refinement/particle_density.cc
+++ b/source/mesh_refinement/particle_density.cc
@@ -36,7 +36,7 @@ namespace aspect
                              "postprocessor plugin `particles' to be selected. Please activate the "
                              "particles or deactivate this mesh refinement plugin."));
 
-      const Particle::ParticleHandler<dim> &particle_handler = this->get_particle_world().get_particle_handler();
+      const Particle::ParticleHandler<dim> &particle_handler = this->get_particle_world(0).get_particle_handler();
 
       for (const auto &cell : this->get_dof_handler().active_cell_iterators())
         if (cell->is_locally_owned())
diff --git a/source/particle/integrator/rk_2.cc b/source/particle/integrator/rk_2.cc
index adefbe3daf2..eec91bde25c 100644
--- a/source/particle/integrator/rk_2.cc
+++ b/source/particle/integrator/rk_2.cc
@@ -41,7 +41,7 @@ namespace aspect
       void
       RK2<dim>::initialize ()
       {
-        const auto &property_information = this->get_particle_world().get_property_manager().get_data_info();
+        const auto &property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
         property_index_old_location = property_information.get_position_by_field_name("internal: integrator properties");
       }
 
diff --git a/source/particle/integrator/rk_4.cc b/source/particle/integrator/rk_4.cc
index 59faa2cad7d..e94b454783f 100644
--- a/source/particle/integrator/rk_4.cc
+++ b/source/particle/integrator/rk_4.cc
@@ -41,7 +41,7 @@ namespace aspect
       void
       RK4<dim>::initialize ()
       {
-        const auto &property_information = this->get_particle_world().get_property_manager().get_data_info();
+        const auto &property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
 
         property_indices[0] = property_information.get_position_by_field_name("internal: integrator properties");
         property_indices[1] = property_indices[0] + dim;
diff --git a/source/particle/interpolator/bilinear_least_squares.cc b/source/particle/interpolator/bilinear_least_squares.cc
index 2b1ccf75e0c..da2fc320bce 100644
--- a/source/particle/interpolator/bilinear_least_squares.cc
+++ b/source/particle/interpolator/bilinear_least_squares.cc
@@ -281,7 +281,7 @@ namespace aspect
         {
           prm.enter_subsection("Bilinear least squares");
           {
-            const auto &particle_property_information = this->get_particle_world().get_property_manager().get_data_info();
+            const auto &particle_property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
             const unsigned int n_property_components = particle_property_information.n_components();
             const unsigned int n_internal_components = particle_property_information.get_components_by_field_name("internal: integrator properties");
 
diff --git a/source/particle/interpolator/quadratic_least_squares.cc b/source/particle/interpolator/quadratic_least_squares.cc
index 80a14824ec9..4dfdfe715f6 100644
--- a/source/particle/interpolator/quadratic_least_squares.cc
+++ b/source/particle/interpolator/quadratic_least_squares.cc
@@ -558,7 +558,7 @@ namespace aspect
         {
           prm.enter_subsection("Quadratic least squares");
           {
-            const auto &particle_property_information = this->get_particle_world().get_property_manager().get_data_info();
+            const auto &particle_property_information = this->get_particle_world(this->get_particle_world_index()).get_property_manager().get_data_info();
             const unsigned int n_property_components = particle_property_information.n_components();
             const unsigned int n_internal_components = particle_property_information.get_components_by_field_name("internal: integrator properties");
 
diff --git a/source/particle/property/cpo_bingham_average.cc b/source/particle/property/cpo_bingham_average.cc
index 3e85ce069b9..8ccba63ae27 100644
--- a/source/particle/property/cpo_bingham_average.cc
+++ b/source/particle/property/cpo_bingham_average.cc
@@ -35,7 +35,7 @@ namespace aspect
       {
         const unsigned int my_rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
         this->random_number_generator.seed(random_number_seed+my_rank);
-        const auto &manager = this->get_particle_world().get_property_manager();
+        const auto &manager = this->get_particle_world(this->get_particle_world_index()).get_property_manager();
         AssertThrow(manager.plugin_name_exists("crystal preferred orientation"),
                     ExcMessage("No crystal preferred orientation property plugin found."));
 
@@ -252,7 +252,7 @@ namespace aspect
         {
           // Get a pointer to the CPO particle property.
           cpo_particle_property = std::make_unique<const Particle::Property::CrystalPreferredOrientation<dim>> (
-                                    this->get_particle_world().get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>());
+                                    this->get_particle_world(this->get_particle_world_index()).get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>());
 
           random_number_seed = prm.get_integer ("Random number seed");
           n_grains = cpo_particle_property->get_number_of_grains();
diff --git a/source/particle/property/cpo_elastic_tensor.cc b/source/particle/property/cpo_elastic_tensor.cc
index ddc2080fc14..e15c4436b70 100644
--- a/source/particle/property/cpo_elastic_tensor.cc
+++ b/source/particle/property/cpo_elastic_tensor.cc
@@ -66,7 +66,7 @@ namespace aspect
       void
       CpoElasticTensor<dim>::initialize ()
       {
-        const auto &manager = this->get_particle_world().get_property_manager();
+        const auto &manager = this->get_particle_world(this->get_particle_world_index()).get_property_manager();
         AssertThrow(manager.plugin_name_exists("crystal preferred orientation"),
                     ExcMessage("No crystal preferred orientation property plugin found."));
 
@@ -147,7 +147,7 @@ namespace aspect
       {
         // Get a reference to the CPO particle property.
         const Particle::Property::CrystalPreferredOrientation<dim> &cpo_particle_property =
-          this->get_particle_world().get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>();
+          this->get_particle_world(this->get_particle_world_index()).get_property_manager().template get_matching_property<Particle::Property::CrystalPreferredOrientation<dim>>();
 
 
         const SymmetricTensor<2,6> C_average = voigt_average_elastic_tensor(cpo_particle_property,
diff --git a/source/particle/property/elastic_tensor_decomposition.cc b/source/particle/property/elastic_tensor_decomposition.cc
index 4fd5a7f91ac..25308ff2623 100644
--- a/source/particle/property/elastic_tensor_decomposition.cc
+++ b/source/particle/property/elastic_tensor_decomposition.cc
@@ -270,7 +270,7 @@ namespace aspect
       void
       ElasticTensorDecomposition<dim>::initialize ()
       {
-        const Particle::Property::Manager<dim> &manager = this->get_particle_world().get_property_manager();
+        const Particle::Property::Manager<dim> &manager = this->get_particle_world(this->get_particle_world_index()).get_property_manager();
         AssertThrow(manager.plugin_name_exists("crystal preferred orientation"),
                     ExcMessage("No cpo property plugin found."));
         AssertThrow(manager.plugin_name_exists("cpo elastic tensor"),
diff --git a/source/postprocess/crystal_preferred_orientation.cc b/source/postprocess/crystal_preferred_orientation.cc
index 50e38562de3..974ffa07e79 100644
--- a/source/postprocess/crystal_preferred_orientation.cc
+++ b/source/postprocess/crystal_preferred_orientation.cc
@@ -173,8 +173,8 @@ namespace aspect
     CrystalPreferredOrientation<dim>::execute (TableHandler &statistics)
     {
 
-      const Particle::Property::Manager<dim> &manager = this->get_particle_world().get_property_manager();
-      const Particle::Property::ParticleHandler<dim> &particle_handler = this->get_particle_world().get_particle_handler();
+      const Particle::Property::Manager<dim> &manager = this->get_particle_world(0).get_property_manager();
+      const Particle::Property::ParticleHandler<dim> &particle_handler = this->get_particle_world(0).get_particle_handler();
 
       const bool cpo_elastic_decomposition_plugin_exists = manager.plugin_name_exists("elastic tensor decomposition");
 
@@ -220,7 +220,7 @@ namespace aspect
                                                                         + "isotropic_norm_square") : "") << std::endl;
 
       // get particle data
-      const Particle::Property::ParticlePropertyInformation &property_information = this->get_particle_world().get_property_manager().get_data_info();
+      const Particle::Property::ParticlePropertyInformation &property_information = this->get_particle_world(0).get_property_manager().get_data_info();
 
       AssertThrow(property_information.fieldname_exists("cpo mineral 0 type") ,
                   ExcMessage("No CPO particle properties found. Make sure that the CPO particle property plugin is selected."));
diff --git a/source/postprocess/load_balance_statistics.cc b/source/postprocess/load_balance_statistics.cc
index adecef1cc8b..2749d974c50 100644
--- a/source/postprocess/load_balance_statistics.cc
+++ b/source/postprocess/load_balance_statistics.cc
@@ -45,7 +45,7 @@ namespace aspect
 
       if (this->n_particle_worlds() > 0)
         {
-          const unsigned int locally_owned_particles = this->get_particle_world().
+          const unsigned int locally_owned_particles = this->get_particle_world(0).
                                                        get_particle_handler().n_locally_owned_particles();
           const dealii::Utilities::MPI::MinMaxAvg particles_per_process =
             dealii::Utilities::MPI::min_max_avg(locally_owned_particles,this->get_mpi_communicator());
diff --git a/source/postprocess/particle_count_statistics.cc b/source/postprocess/particle_count_statistics.cc
index 3004f37ce9b..8b2ca0c952b 100644
--- a/source/postprocess/particle_count_statistics.cc
+++ b/source/postprocess/particle_count_statistics.cc
@@ -35,11 +35,11 @@ namespace aspect
     ParticleCountStatistics<dim>::execute (TableHandler &statistics)
     {
       const Particle::ParticleHandler<dim> &particle_handler =
-        this->get_particle_world().get_particle_handler();
+        this->get_particle_world(0).get_particle_handler();
 
       unsigned int local_min_particles = std::numeric_limits<unsigned int>::max();
       unsigned int local_max_particles = 0;
-      const Particle::types::particle_index global_particles = this->get_particle_world().n_global_particles();
+      const Particle::types::particle_index global_particles = this->get_particle_world(0).n_global_particles();
 
       // compute local min/max
       for (const auto &cell : this->get_dof_handler().active_cell_iterators())
diff --git a/source/postprocess/particles.cc b/source/postprocess/particles.cc
index 348b4b9da7a..126f6a1f607 100644
--- a/source/postprocess/particles.cc
+++ b/source/postprocess/particles.cc
@@ -323,7 +323,7 @@ namespace aspect
       if (std::isnan(last_output_time))
         last_output_time = this->get_time() - output_interval;
 
-      const Particle::World<dim> &world = this->get_particle_world();
+      const Particle::World<dim> &world = this->get_particle_world(0);
 
       statistics.add_value("Number of advected particles",world.n_global_particles());
 
@@ -562,7 +562,7 @@ namespace aspect
       std::ostringstream os;
       aspect::oarchive oa (os);
 
-      this->get_particle_world().save(os);
+      this->get_particle_world(0).save(os);
       oa << (*this);
 
       status_strings["Particles"] = os.str();
@@ -580,7 +580,7 @@ namespace aspect
           aspect::iarchive ia (is);
 
           // Load the particle world
-          this->get_particle_world().load(is);
+          this->get_particle_world(0).load(is);
 
           ia >> (*this);
         }
diff --git a/source/postprocess/visualization/particle_count.cc b/source/postprocess/visualization/particle_count.cc
index 715a44f72c0..f84d67b7386 100644
--- a/source/postprocess/visualization/particle_count.cc
+++ b/source/postprocess/visualization/particle_count.cc
@@ -43,7 +43,7 @@ namespace aspect
       ParticleCount<dim>::execute() const
       {
         const Particle::ParticleHandler<dim> &particle_handler =
-          this->get_particle_world().get_particle_handler();
+          this->get_particle_world(0).get_particle_handler();
 
         std::pair<std::string, std::unique_ptr<Vector<float>>>
         return_value ("particles_per_cell",
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index aaa31fbd920..d74c564d1ea 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -189,7 +189,7 @@ namespace aspect
                    nullptr),
 #endif
     boundary_heat_flux (BoundaryHeatFlux::create_boundary_heat_flux<dim>(prm)),
-    particle_world(nullptr),
+    particle_worlds(),
     time (numbers::signaling_nan<double>()),
     time_step (numbers::signaling_nan<double>()),
     old_time_step (numbers::signaling_nan<double>()),
@@ -448,12 +448,16 @@ namespace aspect
 
     if (postprocess_manager.template has_matching_postprocessor<Postprocess::Particles<dim>>())
       {
-        particle_world = std::make_unique<Particle::World<dim>>();
-        if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(particle_world.get()))
-          sim->initialize_simulator (*this);
+        particle_worlds.emplace_back(std::move(std::make_unique<Particle::World<dim>>()));
+        for (unsigned int particle_world_index = 0 ; particle_world_index < particle_worlds.size(); ++particle_world_index)
+          {
+            if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(particle_worlds[particle_world_index].get()))
+              sim->initialize_simulator (*this);
+
+            particle_worlds.back()->parse_parameters(prm,particle_world_index);
+            particle_worlds.back()->initialize();
 
-        particle_world->parse_parameters(prm,0);
-        particle_world->initialize();
+          }
       }
 
     mesh_refinement_manager.initialize_simulator (*this);
@@ -568,7 +572,9 @@ namespace aspect
     // is destroyed after the latter. But it stores a pointer to the
     // triangulation and uses it during destruction. This results in
     // trouble. So destroy it first.
-    particle_world.reset();
+
+    for (auto &particle_world : particle_worlds)
+      particle_world.reset();
 
     // wait if there is a thread that's still writing the statistics
     // object (set from the output_statistics() function)
@@ -608,9 +614,10 @@ namespace aspect
 
     // Copy particle handler to restore particle location and properties
     // before repeating a timestep
-    if (particle_world.get() != nullptr)
+    for (auto &particle_world : particle_worlds)
       particle_world->backup_particles();
 
+
     // then interpolate the current boundary velocities. copy constraints
     // into current_constraints and then add to current_constraints
     compute_current_constraints ();
@@ -642,7 +649,7 @@ namespace aspect
     if (prescribed_stokes_solution.get())
       prescribed_stokes_solution->update();
 
-    if (particle_world.get() != nullptr)
+    for (auto &particle_world : particle_worlds)
       particle_world->update();
 
     // do the same for the traction boundary conditions and other things
@@ -2133,7 +2140,7 @@ namespace aspect
             // Restore particles through stored copy of particle handler,
             // created in start_timestep(),
             // but only if this timestep is to be repeated.
-            if (particle_world.get() != nullptr)
+            for (auto &particle_world : particle_worlds)
               particle_world->restore_particles();
 
             continue; // repeat time step loop
diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc
index 6f21ded96ab..836a95c0f66 100644
--- a/source/simulator/initial_conditions.cc
+++ b/source/simulator/initial_conditions.cc
@@ -284,8 +284,8 @@ namespace aspect
     // need to write into it and we can not
     // write into vectors with ghost elements
 
-    const Particle::Interpolator::Interface<dim> &particle_interpolator = particle_world->get_interpolator();
-    const Particle::Property::Manager<dim> &particle_property_manager = particle_world->get_property_manager();
+    const Particle::Interpolator::Interface<dim> &particle_interpolator = particle_worlds[0]->get_interpolator();
+    const Particle::Property::Manager<dim> &particle_property_manager = particle_worlds[0]->get_property_manager();
 
     std::vector<unsigned int> particle_property_indices;
     ComponentMask property_mask  (particle_property_manager.get_data_info().n_components(),false);
@@ -356,7 +356,7 @@ namespace aspect
           try
             {
               particle_properties =
-                particle_interpolator.properties_at_points(particle_world->get_particle_handler(),
+                particle_interpolator.properties_at_points(particle_worlds[0]->get_particle_handler(),
                                                            quadrature_points,
                                                            property_mask,
                                                            cell);
diff --git a/source/simulator/simulator_access.cc b/source/simulator/simulator_access.cc
index e538cc06cc2..990c55a8564 100644
--- a/source/simulator/simulator_access.cc
+++ b/source/simulator/simulator_access.cc
@@ -821,30 +821,27 @@ namespace aspect
   unsigned int
   SimulatorAccess<dim>::n_particle_worlds() const
   {
-    // operator () returns whether an object is managed by a unique ptr
-    return (simulator->particle_world ? 1 : 0);
+    return simulator->particle_worlds.size();
   }
 
 
 
   template <int dim>
   const Particle::World<dim> &
-  SimulatorAccess<dim>::get_particle_world() const
+  SimulatorAccess<dim>::get_particle_world(unsigned int particle_world_index) const
   {
-    Assert (simulator->particle_world.get() != nullptr,
-            ExcMessage("You can not call this function if there is no particle world."));
-    return *simulator->particle_world.get();
+    AssertThrow (particle_world_index < simulator->particle_worlds.size(), ExcInternalError());
+    return *simulator->particle_worlds[particle_world_index].get();
   }
 
 
 
   template <int dim>
   Particle::World<dim> &
-  SimulatorAccess<dim>::get_particle_world()
+  SimulatorAccess<dim>::get_particle_world(unsigned int particle_world_index)
   {
-    Assert (simulator->particle_world.get() != nullptr,
-            ExcMessage("You can not call this function if there is no particle world."));
-    return *simulator->particle_world.get();
+    AssertThrow (particle_world_index < simulator->particle_worlds.size(), ExcInternalError());
+    return *simulator->particle_worlds[particle_world_index].get();
   }
 
 
diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc
index de5efc2a230..ec59003c77c 100644
--- a/source/simulator/solver_schemes.cc
+++ b/source/simulator/solver_schemes.cc
@@ -226,7 +226,7 @@ namespace aspect
   {
     // Advect the particles before they are potentially used to
     // set up the compositional fields.
-    if (particle_world.get() != nullptr)
+    for (auto &particle_world : particle_worlds)
       {
         // Do not advect the particles in the initial refinement stage
         const bool in_initial_refinement = (timestep_number == 0)
@@ -958,8 +958,10 @@ namespace aspect
         // Restore particles through stored copy of particle handler,
         // but only if they have already been displaced in a nonlinear
         // iteration (in the assemble_and_solve_composition call).
-        if ((particle_world.get() != nullptr) && (nonlinear_iteration > 0))
-          particle_world->restore_particles();
+
+        if (nonlinear_iteration > 0)
+          for (auto &particle_world : particle_worlds)
+            particle_world->restore_particles();
 
         const double relative_temperature_residual =
           assemble_and_solve_temperature(initial_temperature_residual,
@@ -1115,8 +1117,10 @@ namespace aspect
         // Restore particles through stored copy of particle handler,
         // but only if they have already been displaced in a nonlinear
         // iteration (in the assemble_and_solve_composition call).
-        if ((particle_world.get() != nullptr) && (nonlinear_iteration > 0))
-          particle_world->restore_particles();
+
+        if (nonlinear_iteration > 0)
+          for (auto &particle_world : particle_worlds)
+            particle_world->restore_particles();
 
         const double relative_temperature_residual =
           assemble_and_solve_temperature(initial_temperature_residual,
diff --git a/tests/mesh_deformation_particles.cc b/tests/mesh_deformation_particles.cc
index dbcb9cf516c..9e7806dd0cd 100644
--- a/tests/mesh_deformation_particles.cc
+++ b/tests/mesh_deformation_particles.cc
@@ -34,7 +34,7 @@ void post_mesh_deformation (const SimulatorAccess<dim> &simulator_access)
   // Get the reference location of the one particle,
   // and check that it has the correct value.
   const Point<dim> predicted_particle_reference_location =
-    simulator_access.get_particle_world().get_particle_handler().begin()->get_reference_location();
+    simulator_access.get_particle_world(0).get_particle_handler().begin()->get_reference_location();
   // At timestep 0, the mesh is not deformed
   // other than through the initial topography/
   // initial mesh deformation plugins, which
diff --git a/tests/particle_handler_copy.cc b/tests/particle_handler_copy.cc
index 5e6766b2e12..86abc4a9c1b 100644
--- a/tests/particle_handler_copy.cc
+++ b/tests/particle_handler_copy.cc
@@ -33,7 +33,7 @@ void duplicate_particle_handler(const SimulatorAccess<dim> &simulator_access,
                                 const unsigned int,
                                 const SolverControl &)
 {
-  const auto &particle_world = simulator_access.get_particle_world();
+  const auto &particle_world = simulator_access.get_particle_world(0);
   dealii::Particles::ParticleHandler<dim> particle_handler;
   std::cout << "duplicating particle handler" << std::endl;
 
diff --git a/tests/particle_property_post_initialize_function.cc b/tests/particle_property_post_initialize_function.cc
index 14ce6fe7adc..d5c7d46c01e 100644
--- a/tests/particle_property_post_initialize_function.cc
+++ b/tests/particle_property_post_initialize_function.cc
@@ -44,7 +44,7 @@ namespace aspect
           initialize ()
           {
             std::cout << "initialize" << std::endl;
-            const Particle::Property::Manager<dim> &manager = this->get_particle_world().get_property_manager();
+            const Particle::Property::Manager<dim> &manager = this->get_particle_world(this->get_particle_world_index()).get_property_manager();
             post_initialized_info = manager.get_data_info().get_field_index_by_name("initial position");
             std::cout << "initial position: post_initialized_info = " << post_initialized_info << std::endl;