diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7bfe0b8b..8c2fc264 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -27,6 +27,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   types of convolution functions, so we simply convert the grid to use only
   proton PDFs
 
+### Changed
+
+- the function `Grid::evolve` now makes use of parallelization to take advantage
+  of the number of CPU cores available using the Rayon crate; the number of CPU
+  cores to be used can be controlled via the `RAYON_NUM_THREADS` environment
+  variable
+
 ## [1.1.0] - 08/07/2025
 
 ### Added
diff --git a/Cargo.lock b/Cargo.lock
index c5559f48..7dd0fe3d 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -968,6 +968,7 @@ dependencies = [
  "pineappl 0.8.3",
  "rand",
  "rand_pcg",
+ "rayon",
  "rustc-hash 1.1.0",
  "serde",
  "serde_yaml",
diff --git a/pineappl/Cargo.toml b/pineappl/Cargo.toml
index 9cc1254c..b07573fd 100644
--- a/pineappl/Cargo.toml
+++ b/pineappl/Cargo.toml
@@ -27,6 +27,7 @@ itertools = "0.10.1"
 lz4_flex = "0.9.2"
 ndarray = { features = ["serde"], version = "0.15.4" }
 pineappl-v0 = { package = "pineappl", version = "0.8.2" }
+rayon = "1.5.1"
 rustc-hash = "1.1.0"
 serde = { features = ["derive"], version = "1.0.130" }
 thiserror = "1.0.30"
diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs
index a9225893..f7d03dbe 100644
--- a/pineappl/src/evolution.rs
+++ b/pineappl/src/evolution.rs
@@ -15,6 +15,7 @@ use ndarray::{
     s, Array1, Array2, Array3, ArrayD, ArrayView1, ArrayView4, ArrayViewD, ArrayViewMutD, Axis,
     Ix1, Ix2,
 };
+use rayon::iter::{IndexedParallelIterator, IntoParallelRefMutIterator, ParallelIterator};
 use std::iter;
 
 /// This structure captures the information needed to create an evolution kernel operator (EKO) for
@@ -400,8 +401,6 @@ pub(crate) fn evolve_slice(
     let channels0 = channels0;
 
     let mut sub_fk_tables = Vec::with_capacity(grid.bwfl().len() * channels0.len());
-
-    // TODO: generalize to `n`
     let mut last_x1 = vec![Vec::new(); infos.len()];
     let mut eko_slices = vec![Vec::new(); infos.len()];
     let dim: Vec<_> = infos.iter().map(|info| info.x0.len()).collect();
@@ -445,34 +444,35 @@ pub(crate) fn evolve_slice(
                 }
             }
 
-            for (pids1, factor) in channel1.entry() {
-                // find the tuple of EKOs that evolve the current channel entry of the grid into
-                // every channel of the FK-table
-                let tmp = channels0.iter().map(|pids0| {
-                    izip!(pids0, pids1, &pids01, &eko_slices)
-                        .map(|(&pid0, &pid1, pids, slices)| {
-                            // for each convolution ...
-                            pids.iter().zip(slices).find_map(|(&(p0, p1), op)| {
-                                // find the EKO that matches both the FK-table and the grid PID
-                                ((p0 == pid0) && (p1 == pid1)).then_some(op)
+            // for each channel in the FK-table ...
+            tables
+                .par_iter_mut()
+                .zip(&channels0)
+                .for_each(|(fk_table, pids0)| {
+                    // and each sub-channel in the grid ...
+                    for (pids1, factor) in channel1.entry() {
+                        // find the tuple of EKOs that evolve this combination
+                        let ops: Option<Box<[_]>> = izip!(pids0, pids1, &pids01, &eko_slices)
+                            .map(|(&pid0, &pid1, pids, slices)| {
+                                // for each convolution ...
+                                pids.iter().zip(slices).find_map(|(pid01, op)| {
+                                    // find the EKO that matches both the FK-table and the grid PID
+                                    (pid01 == &(pid0, pid1)).then_some(op)
+                                })
                             })
-                        })
-                        // if an EKO isn't found, it's zero and therefore the whole FK-table
-                        // channel contribution will be zero
-                        .collect::<Option<Box<[_]>>>()
-                });
-
-                for (fk_table, ops) in tables.iter_mut().zip(tmp) {
-                    // if there's one zero EKO, the entire tuple is `None`
-                    if let Some(ops) = ops {
-                        general_tensor_mul(*factor, array.view(), &ops, fk_table.view_mut());
+                            // if an EKO isn't found, it's zero and therefore the whole FK-table
+                            // channel contribution will be zero
+                            .collect();
+
+                        // if there's one zero EKO, the entire tuple is `None`
+                        if let Some(ops) = ops {
+                            general_tensor_mul(*factor, array.view(), &ops, fk_table.view_mut());
+                        }
                     }
-                }
-            }
+                });
         }
 
         let mut node_values = vec![scale_values.to_vec()];
-
         for info in infos {
             node_values.push(info.x0.clone());
         }