diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 26c34ebd41d..f5468617a4f 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -228,6 +228,13 @@ class Cell {
   //! \return Density in [g/cm3]
   double density(int32_t instance = -1) const;
 
+  //! Get the importance of a cell instance
+  //! \param[in] index 0 for neutrons and 1 for photons.
+  //! \param[in] instance Instance index. If -1 is given, the importance for
+  //! the first instance is returned.
+  //! \return Importance
+  double importance(int index, int32_t instance = -1) const;
+
   //! Set the temperature of a cell instance
   //! \param[in] T Temperature in [K]
   //! \param[in] instance Instance index. If -1 is given, the temperature for
@@ -359,6 +366,11 @@ class Cell {
   //! \brief Unitless density multiplier(s) within this cell.
   vector<double> density_mult_;
 
+  //! \brief Importance for this cell per particle type
+  //!
+  //! May be multiple for distribcell.
+  vector<vector<double>> importance_ {{1.0}, {1.0}};
+
   //! \brief Neighboring cells in the same universe.
   NeighborList neighbors_;
 
diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h
index 107cc7d1f3e..e3f1d8c0cf0 100644
--- a/include/openmc/geometry.h
+++ b/include/openmc/geometry.h
@@ -6,6 +6,7 @@
 
 #include "openmc/array.h"
 #include "openmc/constants.h"
+#include "openmc/particle.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -52,6 +53,21 @@ bool check_cell_overlap(GeometryState& p, bool error = true);
 
 int cell_instance_at_level(const GeometryState& p, int level);
 
+//==============================================================================
+//! Get the cell importance for a particle at the specified
+//! type and universe level.
+//!
+//! \param p A particle for which to compute the instance using
+//!   its coordinates
+//! \param type ParticleType to get importance for
+//! \param level The level (zero indexed) of the geometry where the instance
+//! should be computed. \return The importance of the cell at the specified
+//! level.
+//==============================================================================
+
+double cell_importance_at_level(
+  const GeometryState& p, const ParticleType type, int level);
+
 //==============================================================================
 //! Locate a particle in the geometry tree and set its geometry data fields.
 //!
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 5b632bbce53..717f6da03b4 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -501,6 +501,7 @@ class ParticleData : public GeometryState {
   int g_ {0};
   int g_last_;
 
+  double importance_last_ {1.0};
   double wgt_ {1.0};
   double wgt_born_ {1.0};
   double wgt_ww_born_ {-1.0};
@@ -610,6 +611,9 @@ class ParticleData : public GeometryState {
   int& g_last() { return g_last_; }
   const int& g_last() const { return g_last_; }
 
+  double& importance_last() { return importance_last_; }
+  const double& importance_last() const { return importance_last_; }
+
   // Statistic weight of particle. Setting to zero indicates that the particle
   // is dead.
   double& wgt() { return wgt_; }
diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
index 9a6cf1b2131..fa23ddff971 100644
--- a/include/openmc/simulation.h
+++ b/include/openmc/simulation.h
@@ -22,6 +22,7 @@ constexpr int STATUS_EXIT_ON_TRIGGER {2};
 
 namespace simulation {
 
+extern bool cell_importances; //!< Use cell importances variance reduction
 extern int ct_current_file;   //!< current collision track file index
 extern "C" int current_batch; //!< current batch
 extern "C" int current_gen;   //!< current fission generation
diff --git a/include/openmc/source.h b/include/openmc/source.h
index 2f32aa2a058..0b3979c1291 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -91,7 +91,7 @@ class Source {
 
   // Methods for constraints
   void read_constraints(pugi::xml_node node);
-  bool satisfies_spatial_constraints(Position r) const;
+  bool satisfies_spatial_constraints(Position r, ParticleType type) const;
   bool satisfies_energy_constraints(double E) const;
   bool satisfies_time_constraints(double time) const;
 
diff --git a/openmc/cell.py b/openmc/cell.py
index 499cf950429..22c36270731 100644
--- a/openmc/cell.py
+++ b/openmc/cell.py
@@ -1,4 +1,4 @@
-from collections.abc import Iterable
+from collections.abc import Iterable, Mapping
 from math import cos, sin, pi
 from numbers import Real
 
@@ -69,6 +69,9 @@ class Cell(IDManagerMixin):
     rotation_matrix : numpy.ndarray
         The rotation matrix defined by the angles specified in the
         :attr:`Cell.rotation` property.
+    importance : dict
+        Dictionary containing importance of the cell for neutrons and photons. Defaults to 1.0.
+        Multiple importances can be given to give each distributed cell instance a unique importance.
     temperature : float or iterable of float
         Temperature of the cell in Kelvin.  Multiple temperatures can be given
         to give each distributed cell instance a unique temperature.
@@ -110,6 +113,7 @@ def __init__(self, cell_id=None, name='', fill=None, region=None):
         self.name = name
         self.fill = fill
         self.region = region
+        self._importance = {}
         self._rotation = None
         self._rotation_matrix = None
         self._temperature = None
@@ -238,6 +242,35 @@ def rotation(self, rotation):
     @property
     def rotation_matrix(self):
         return self._rotation_matrix
+    
+    @property
+    def importance(self) -> dict:
+        return self._importance
+
+    @importance.setter
+    def importance(self, importance):
+        # Make sure importances are positive
+        cv.check_type('cell importance', importance, (Mapping, Iterable, Real))
+        if isinstance(importance, Mapping):
+            for key, value in importance.items():
+                cv.check_type('cell importance key', key, ('neutron', 'photon'))
+                cv.check_type('cell importance value', value, (Iterable, Real))
+                if isinstance(value, Iterable):
+                    cv.check_type('cell importance value', value, Iterable, Real)
+                    for i, val in enumerate(importance):
+                        cv.check_greater_than(f'cell importance value[{i}]', val, 0.0, True)
+                else:
+                    cv.check_greater_than('cell importance value', value, 0.0, True)
+            self._importance = importance
+                
+        else:
+            if isinstance(importance, Iterable):
+                cv.check_type('cell importance', importance, Iterable, Real)
+                for i, val in enumerate(importance):
+                    cv.check_greater_than(f'cell importance[{i}]', val, 0.0, True)
+            else:
+                cv.check_greater_than('cell importance', importance, 0.0, True)
+            self._importance = {'neutron': importance, 'photon': importance}
 
     @property
     def temperature(self):
@@ -692,6 +725,14 @@ def create_surface_elements(node, element, memo=None):
             else:
                 element.set("density", str(self.density))
 
+        for key in ('neutron', 'photon'):
+            if self.importance.get(key, 1.0) != 1.0:
+                if isinstance(self.importance[key], Iterable):
+                    importance_subelement = ET.SubElement(element, f"importance_{key}")
+                    importance_subelement.text = ' '.join(str(val) for val in self.importance[key])
+                else:
+                    element.set("importance", str(self.importance[key]))
+
         if self.translation is not None:
             element.set("translation", ' '.join(map(str, self.translation)))
 
@@ -750,6 +791,16 @@ def from_xml_element(cls, elem, surfaces, materials, get_universe):
         v = get_text(elem, 'volume')
         if v is not None:
             c.volume = float(v)
+        importance = {}
+        
+        for key in ('neutron', 'photon'):
+            values = get_elem_list(elem, f"importance_{key}", float)
+            if values is not None:
+                if len(values) == 1:
+                    values = values[0]
+                importance[f"importance_{key}"] = values
+        setattr(c, "importance", importance)
+        
         for key in ('temperature', 'density', 'rotation', 'translation'):
             values = get_elem_list(elem, key, float)
             if values is not None:
diff --git a/src/cell.cpp b/src/cell.cpp
index ebe28c3d2ce..8c2d19054da 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -22,6 +22,7 @@
 #include "openmc/material.h"
 #include "openmc/nuclide.h"
 #include "openmc/settings.h"
+#include "openmc/simulation.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -101,6 +102,15 @@ double Cell::temperature(int32_t instance) const
   }
 }
 
+double Cell::importance(int index, int32_t instance) const
+{
+  if (instance >= 0) {
+    return importance_[index][importance_[index].size() > 1 ? instance : 0];
+  } else {
+    return importance_[index][0];
+  }
+}
+
 double Cell::density_mult(int32_t instance) const
 {
   if (instance >= 0) {
@@ -337,6 +347,18 @@ void Cell::to_hdf5(hid_t cell_group) const
     write_dataset(group, "lattice", model::lattices[fill_]->id_);
   }
 
+  if (importance_[0].size() == 1) {
+    if (importance_[0][0] != 1.0)
+      write_dataset(group, "importance_neutron", importance_[0][0]);
+  } else {
+    write_dataset(group, "importance_neutron", importance_[0]);
+  }
+  if (importance_[1].size() == 1) {
+    if (importance_[1][0] != 1.0)
+      write_dataset(group, "importance_photon", importance_[1][0]);
+  } else {
+    write_dataset(group, "importance_photon", importance_[1]);
+  }
   close_group(group);
 }
 
@@ -472,6 +494,38 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
     }
   }
 
+  if (check_for_node(cell_node, "importance_neutron")) {
+    importance_[0] = get_node_array<double>(cell_node, "importance_neutron");
+    importance_[0].shrink_to_fit();
+
+    // Make sure all importancs are non-negative and greater than zero.
+    for (auto importance : importance_[0]) {
+      if (importance < 0)
+        fatal_error(fmt::format(
+          "Cell {} was specified with a negative neutron importance.", id_));
+      if (!settings::weight_windows_on) {
+        if (importance != 1.0 && importance > 0.0)
+          simulation::cell_importances = true;
+      }
+    }
+  }
+
+  if (check_for_node(cell_node, "importance_photon")) {
+    importance_[1] = get_node_array<double>(cell_node, "importance_photon");
+    importance_[1].shrink_to_fit();
+
+    // Make sure all importancs are non-negative and greater than zero.
+    for (auto importance : importance_[1]) {
+      if (importance < 0)
+        fatal_error(fmt::format(
+          "Cell {} was specified with a negative photon importance.", id_));
+      if (!settings::weight_windows_on) {
+        if (importance != 1.0 && importance > 0.0)
+          simulation::cell_importances = true;
+      }
+    }
+  }
+
   // Read the region specification.
   std::string region_spec;
   if (check_for_node(cell_node, "region")) {
@@ -1114,6 +1168,7 @@ vector<int32_t> Region::surfaces() const
 
 void read_cells(pugi::xml_node node)
 {
+  simulation::cell_importances = false;
   // Count the number of cells.
   int n_cells = 0;
   for (pugi::xml_node cell_node : node.children("cell")) {
diff --git a/src/geometry.cpp b/src/geometry.cpp
index ddb61385f18..6743c4cbb1d 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -97,6 +97,47 @@ int cell_instance_at_level(const GeometryState& p, int level)
   return instance;
 }
 
+double cell_importance_at_level(
+  const GeometryState& p, ParticleType type, int level)
+{
+  // throw error if the requested level is too deep for the geometry
+  if (level > model::n_coord_levels) {
+    fatal_error(
+      fmt::format("Cell importance at level {} requested, but only {} "
+                  "levels exist in the geometry.",
+        level, p.n_coord()));
+  }
+  int j = -1;
+  if (type.is_neutron())
+    j = 0;
+  else if (type.is_photon())
+    j = 1;
+  else
+    return 1.0;
+
+  // determine the cell instance
+  Cell& c {*model::cells[p.coord(level).cell()]};
+
+  // compute the cell's instance and importance
+  int instance = 0.0;
+  double importance = 1.0;
+  for (int i = 0; i < level; i++) {
+    const auto& c_i {*model::cells[p.coord(i).cell()]};
+    if (c_i.type_ == Fill::UNIVERSE) {
+      instance += c_i.offset_[c.distribcell_index_];
+    } else if (c_i.type_ == Fill::LATTICE) {
+      instance += c_i.offset_[c.distribcell_index_];
+      auto& lat {*model::lattices[p.coord(i + 1).lattice()]};
+      const auto& i_xyz {p.coord(i + 1).lattice_index()};
+      if (lat.are_valid_indices(i_xyz)) {
+        instance += lat.offset(c.distribcell_index_, i_xyz);
+      }
+    }
+    importance *= c_i.importance(j, instance);
+  }
+  return importance;
+}
+
 //==============================================================================
 
 bool find_cell_inner(
diff --git a/src/particle.cpp b/src/particle.cpp
index 0a063559824..510a3598c55 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -298,6 +298,9 @@ void Particle::event_cross_surface()
   }
   n_coord_last() = n_coord();
 
+  if (simulation::cell_importances && material() != MATERIAL_VOID)
+    importance_last() = cell_importance_at_level(*this, type(), n_coord() - 1);
+
   // Set surface that particle is on and adjust coordinate levels
   surface() = boundary().surface();
   n_coord() = boundary().coord_level();
@@ -319,6 +322,11 @@ void Particle::event_cross_surface()
       add_surf_source_to_bank(*this, surf);
     }
     this->cross_surface(surf);
+    double importance = cell_importance_at_level(*this, type(), n_coord() - 1);
+    if (importance == 0.0) {
+      wgt() = 0.0;
+      return;
+    }
     // If no BC, add particle to surface source after crossing surface
     if (surf.surf_source_ && !surf.bc_) {
       add_surf_source_to_bank(*this, surf);
@@ -326,6 +334,35 @@ void Particle::event_cross_surface()
     if (settings::weight_window_checkpoint_surface) {
       apply_weight_windows(*this);
     }
+    if (simulation::cell_importances && material() != MATERIAL_VOID) {
+      if (importance != importance_last()) {
+        if (importance < importance_last()) {
+          if (importance_last() * prn(current_seed()) < importance) {
+            wgt() *= importance_last() / importance;
+          } else {
+            wgt() = 0.;
+            return;
+          }
+        } else {
+          // do not further split the particle if above the limit
+          if (n_split() >= settings::max_history_splits)
+            return;
+          double n_s = importance / importance_last();
+          int n = static_cast<int>(n_s);
+          if (prn(current_seed()) <= n_s - n)
+            ++n;
+          n = std::min(n, settings::max_history_splits);
+          n_split() += n;
+
+          // Create secondaries and divide weight among all particles
+          for (int l = 0; l < n - 1; l++) {
+            split(wgt() / n);
+          }
+          // remaining weight is applied to current particle
+          wgt() /= n;
+        }
+      }
+    }
     event() = TallyEvent::SURFACE;
   }
   // Score cell to cell partial currents
diff --git a/src/physics.cpp b/src/physics.cpp
index 72b4a661192..c471e2da85d 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -2,12 +2,14 @@
 
 #include "openmc/bank.h"
 #include "openmc/bremsstrahlung.h"
+#include "openmc/cell.h"
 #include "openmc/chain.h"
 #include "openmc/constants.h"
 #include "openmc/distribution_multi.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/endf.h"
 #include "openmc/error.h"
+#include "openmc/geometry.h"
 #include "openmc/ifp.h"
 #include "openmc/material.h"
 #include "openmc/math_functions.h"
@@ -139,7 +141,9 @@ void sample_neutron_reaction(Particle& p)
 
   // Create secondary photons
   if (settings::photon_transport) {
-    sample_secondary_photons(p, i_nuclide);
+    if (cell_importance_at_level(p, ParticleType::photon(), p.n_coord() - 1) >
+        0.0)
+      sample_secondary_photons(p, i_nuclide);
   }
 
   // If survival biasing is being used, the following subroutine adjusts the
diff --git a/src/physics_common.cpp b/src/physics_common.cpp
index 10760ce2da1..b1d6c62a874 100644
--- a/src/physics_common.cpp
+++ b/src/physics_common.cpp
@@ -2,6 +2,7 @@
 
 #include "openmc/random_lcg.h"
 #include "openmc/settings.h"
+#include "openmc/simulation.h"
 
 namespace openmc {
 
@@ -24,6 +25,10 @@ void apply_russian_roulette(Particle& p)
   if (!settings::survival_biasing)
     return;
 
+  // Exit if cell importances are used
+  if (simulation::cell_importances)
+    return;
+
   // if survival normalization is on, use normalized weight cutoff and
   // normalized weight survive
   if (settings::survival_normalization) {
diff --git a/src/simulation.cpp b/src/simulation.cpp
index d0d6f037f9e..0d206061f15 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -298,6 +298,7 @@ namespace openmc {
 
 namespace simulation {
 
+bool cell_importances {false};
 int ct_current_file;
 int current_batch;
 int current_gen;
diff --git a/src/source.cpp b/src/source.cpp
index 5300a685f78..bc80a5c273a 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -203,7 +203,7 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const
       accepted = true;
     } else {
       // Check whether sampled site satisfies constraints
-      accepted = satisfies_spatial_constraints(site.r) &&
+      accepted = satisfies_spatial_constraints(site.r, site.particle) &&
                  satisfies_energy_constraints(site.E) &&
                  satisfies_time_constraints(site.time);
       if (!accepted) {
@@ -237,7 +237,7 @@ bool Source::satisfies_time_constraints(double time) const
   return time > time_bounds_.first && time < time_bounds_.second;
 }
 
-bool Source::satisfies_spatial_constraints(Position r) const
+bool Source::satisfies_spatial_constraints(Position r, ParticleType type) const
 {
   GeometryState geom_state;
   geom_state.r() = r;
@@ -248,6 +248,11 @@ bool Source::satisfies_spatial_constraints(Position r) const
   if (!found)
     return false;
 
+  // Reject particle if it is in a zero importance cell
+  if (cell_importance_at_level(geom_state, type, geom_state.n_coord() - 1) ==
+      0.0)
+    return false;
+
   // Check the geometry state against specified domains
   bool accepted = true;
   if (!domain_ids_.empty()) {
@@ -379,7 +384,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
     r_wgt = r_wgt_temp;
 
     // Check if sampled position satisfies spatial constraints
-    accepted = satisfies_spatial_constraints(site.r);
+    accepted = satisfies_spatial_constraints(site.r, site.particle);
 
     // Check for rejection
     if (!accepted) {