diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
index d5f0cb91f0..6ee95201ad 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
@@ -690,6 +690,20 @@ class Persistent_cohomology {
     return result;
   }
 
+  /** @brief Returns persistence intervals for each dimension.
+   * @return A vector of diagrams, one per dimension starting from 0, where each diagram is a vector of persistence intervals (birth and death).
+   */
+  std::vector<std::vector<std::pair<Filtration_value, Filtration_value>>>
+  intervals_by_dimension() {
+    std::vector<std::vector<std::pair<Filtration_value, Filtration_value>>> result;
+    result.resize(dim_max_);
+    for (auto && pair : persistent_pairs_) {
+      auto b = get<0>(pair);
+      result[cpx_->dimension(b)].emplace_back(cpx_->filtration(b), cpx_->filtration(get<1>(pair)));
+    }
+    return result;
+  }
+
  private:
   /*
    * Structure representing a cocycle.
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 547c86c11c..e6f02b756c 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -95,6 +95,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
         vector[int] betti_numbers() nogil
         vector[int] persistent_betti_numbers(double from_value, double to_value) nogil
         vector[pair[double,double]] intervals_in_dimension(int dimension) nogil
+        vector[vector[pair[double,double]]] intervals_by_dimension() nogil
         void write_output_diagram(string diagram_file_name) nogil except +
         vector[pair[vector[int], vector[int]]] persistence_pairs() nogil
         pair[vector[vector[int]], vector[vector[int]]] lower_star_generators() nogil
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 7cf45ec61d..8a4aa32ad3 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -612,7 +612,7 @@ cdef class SimplexTree:
         """
         self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func)
 
-    def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False):
+    def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False, output_type = 'old'):
         """This function computes and returns the persistence of the simplicial complex.
 
         :param homology_coeff_field: The homology coefficient field. Must be a
@@ -627,11 +627,18 @@ cdef class SimplexTree:
             maximal dimension in the complex is computed. If false, it is
             ignored. Default is false.
         :type persistence_dim_max: bool
+        :param output_type: Format of the output. 'old' for the legacy list of (dim,(birth,death)),
+            'array by dimension' for a list of nx2 numpy arrays (one per dimension).
+        :type output_type: str
         :returns: The persistence of the simplicial complex.
         :rtype:  list of pairs(dimension, pair(birth, death))
         """
         self.compute_persistence(homology_coeff_field, min_persistence, persistence_dim_max)
-        return self.pcohptr.get_persistence()
+        if output_type == 'array by dimension':
+            v = self.pcohptr.intervals_by_dimension()
+            return [ np.asarray(dgm) for dgm in v ] # std::move(dgm)?
+        else:
+            return self.pcohptr.get_persistence()
 
     def compute_persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False):
         """This function computes the persistence of the simplicial complex, so it can be accessed through