diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index fb4bb73a..64bc2079 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -18,7 +18,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.10", "3.11", "3.12"]
+        python-version: ["3.10", "3.11", "3.12", "3.13"]
         platform: [ubuntu-latest, macos-latest, windows-latest]
 
     steps:
diff --git a/Cargo.toml b/Cargo.toml
index f6f3cb42..5b3f3e6b 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -9,8 +9,8 @@ name = "_cylindra_ext"
 crate-type = ["cdylib"]
 
 [dependencies]
-pyo3 = "0.23.3"
-numpy = "0.23.0"
+pyo3 = "0.22.6"
+numpy = "0.22.1"  # NOTE: numpy=0.23.0 fails to compile in some OS
 rand = "0.8.5"
 mt19937 = "2.0.1"
 num = { version = "0.4.0", features = ["std"] }
diff --git a/pyproject.toml b/pyproject.toml
index b5377ee4..5cb54e4f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -33,9 +33,9 @@ dependencies = [
     "pyqtgraph>=0.12.4",
     "pyarrow>=11.0.0",
     "numpy>=1.23.0",
-    "scipy>=1.11.3",
+    "scipy>=1.11.3,!=1.15.0",  # See https://github.com/scipy/scipy/issues/22271#issuecomment-2574884519
     "pandas>=1.5.0",
-    "polars>=1.1.0",
+    "polars>=1.19.0",
     "scikit-image>=0.21.0",
     "napari>=0.5.1",
     "qtpy>=2.3.1",
diff --git a/src/alleviate.rs b/src/alleviate.rs
index 99e58668..61f12ad1 100644
--- a/src/alleviate.rs
+++ b/src/alleviate.rs
@@ -44,7 +44,7 @@ pub fn alleviate<'py>(
     let mut arr = arr.to_owned();
     if indices.len() == 0 {
         // nothing to do.
-        return Ok(arr.into_pyarray(py).unbind());
+        return Ok(arr.into_pyarray_bound(py).unbind());
     }
     let geometry = CylinderGeometry::new(ny, na, nrise);
     let mut processed = HashSet::new();
@@ -81,7 +81,7 @@ pub fn alleviate<'py>(
             return value_error!("Infinite recursion occurred.");
         }
     }
-    Ok(arr.into_pyarray(py).unbind())
+    Ok(arr.into_pyarray_bound(py).unbind())
 }
 
 fn local_average(neighbor: &Index, arr: &Array3<f32>, geometry: &CylinderGeometry) -> PyResult<(f32, f32, f32)> {
diff --git a/src/annealing/models/basic.rs b/src/annealing/models/basic.rs
index 4ee4250b..b6335ca2 100644
--- a/src/annealing/models/basic.rs
+++ b/src/annealing/models/basic.rs
@@ -103,25 +103,25 @@ impl CylindricAnnealingModel {
 
     /// Get all the existing distances of longitudinal connections as a numpy array.
     pub fn longitudinal_distances<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_longitudinal_distances().into_pyarray(py).into()
+        self.graph.get_longitudinal_distances().into_pyarray_bound(py).into()
     }
 
     /// Get all the existing distances of lateral connections as a numpy array.
     pub fn lateral_distances<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_lateral_distances().into_pyarray(py).into()
+        self.graph.get_lateral_distances().into_pyarray_bound(py).into()
     }
 
     pub fn longitudinal_angles<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_longitudinal_angles().into_pyarray(py).into()
+        self.graph.get_longitudinal_angles().into_pyarray_bound(py).into()
     }
 
     pub fn lateral_angles<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_lateral_angles().into_pyarray(py).into()
+        self.graph.get_lateral_angles().into_pyarray_bound(py).into()
     }
 
     pub fn get_edge_info<'py>(&self, py: Python<'py>) -> (Py<PyArray2<f32>>, Py<PyArray2<f32>>, Py<PyArray1<i32>>) {
         let (out0, out1, out2) = self.graph.get_edge_states();
-        (out0.into_pyarray(py).into(), out1.into_pyarray(py).into(), out2.into_pyarray(py).into())
+        (out0.into_pyarray_bound(py).into(), out1.into_pyarray_bound(py).into(), out2.into_pyarray_bound(py).into())
     }
 
     #[pyo3(signature = (indices, npf, nrise))]
@@ -199,7 +199,7 @@ impl CylindricAnnealingModel {
 
     /// Get integer shift in each local coordinates as a numpy array.
     pub fn shifts<'py>(&self, py: Python<'py>) -> Py<PyArray2<isize>> {
-        self.graph.get_shifts().into_pyarray(py).into()
+        self.graph.get_shifts().into_pyarray_bound(py).into()
     }
 
     pub fn set_shifts<'py>(
@@ -223,7 +223,7 @@ impl CylindricAnnealingModel {
 
     pub fn binding_energies<'py>(&self, py: Python<'py>) -> (Py<PyArray1<f32>>, Py<PyArray1<f32>>) {
         let (lon, lat) = self.graph.binding_energies();
-        (lon.into_pyarray(py).into(), lat.into_pyarray(py).into())
+        (lon.into_pyarray_bound(py).into(), lat.into_pyarray_bound(py).into())
     }
 
     /// Get current optimization state as a string.
diff --git a/src/annealing/models/defective.rs b/src/annealing/models/defective.rs
index 9319dd30..a10086ac 100644
--- a/src/annealing/models/defective.rs
+++ b/src/annealing/models/defective.rs
@@ -123,25 +123,25 @@ impl DefectiveCylindricAnnealingModel {
 
     /// Get all the existing distances of longitudinal connections as a numpy array.
     pub fn longitudinal_distances<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_longitudinal_distances().into_pyarray(py).into()
+        self.graph.get_longitudinal_distances().into_pyarray_bound(py).into()
     }
 
     /// Get all the existing distances of lateral connections as a numpy array.
     pub fn lateral_distances<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_lateral_distances().into_pyarray(py).into()
+        self.graph.get_lateral_distances().into_pyarray_bound(py).into()
     }
 
     pub fn longitudinal_angles<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_longitudinal_angles().into_pyarray(py).into()
+        self.graph.get_longitudinal_angles().into_pyarray_bound(py).into()
     }
 
     pub fn lateral_angles<'py>(&self, py: Python<'py>) -> Py<PyArray1<f32>> {
-        self.graph.get_lateral_angles().into_pyarray(py).into()
+        self.graph.get_lateral_angles().into_pyarray_bound(py).into()
     }
 
     pub fn get_edge_info<'py>(&self, py: Python<'py>) -> (Py<PyArray2<f32>>, Py<PyArray2<f32>>, Py<PyArray1<i32>>) {
         let (out0, out1, out2) = self.graph.get_edge_states();
-        (out0.into_pyarray(py).into(), out1.into_pyarray(py).into(), out2.into_pyarray(py).into())
+        (out0.into_pyarray_bound(py).into(), out1.into_pyarray_bound(py).into(), out2.into_pyarray_bound(py).into())
     }
 
     #[pyo3(signature = (indices, npf, nrise))]
@@ -219,7 +219,7 @@ impl DefectiveCylindricAnnealingModel {
 
     /// Get integer shift in each local coordinates as a numpy array.
     pub fn shifts<'py>(&self, py: Python<'py>) -> Py<PyArray2<isize>> {
-        self.graph.get_shifts().into_pyarray(py).into()
+        self.graph.get_shifts().into_pyarray_bound(py).into()
     }
 
     pub fn set_shifts<'py>(
@@ -243,7 +243,7 @@ impl DefectiveCylindricAnnealingModel {
 
     pub fn binding_energies<'py>(&self, py: Python<'py>) -> (Py<PyArray1<f32>>, Py<PyArray1<f32>>) {
         let (lon, lat) = self.graph.binding_energies();
-        (lon.into_pyarray(py).into(), lat.into_pyarray(py).into())
+        (lon.into_pyarray_bound(py).into(), lat.into_pyarray_bound(py).into())
     }
 
     /// Get current optimization state as a string.
diff --git a/src/array.rs b/src/array.rs
index cec4520c..a21c5a33 100644
--- a/src/array.rs
+++ b/src/array.rs
@@ -30,7 +30,7 @@ pub fn oblique_coordinates<'py>(
         out[[i, 0]] = (nth + npf * tan1) * d0 + c0;
         out[[i, 1]] = (nth * tan0 + npf) * d1 + c1;
     }
-    Ok(out.into_pyarray(py).unbind())
+    Ok(out.into_pyarray_bound(py).unbind())
 
 }
 
@@ -75,7 +75,7 @@ pub fn cylinder_faces<'py>(
         out_vert[[i, 2]] = v.2;
     }
 
-    Ok(out_vert.into_pyarray(py).unbind())
+    Ok(out_vert.into_pyarray_bound(py).unbind())
 }
 
 #[pyfunction]
diff --git a/src/filters.rs b/src/filters.rs
index 4f8dbfda..2bdf512d 100644
--- a/src/filters.rs
+++ b/src/filters.rs
@@ -83,7 +83,7 @@ impl CylindricArray {
 
     /// Convert the CylindricArray to a 2D numpy array.
     pub fn asarray(&self, py: Python) -> Py<PyArray2<f32>> {
-        self.array.clone().into_pyarray(py).unbind()
+        self.array.clone().into_pyarray_bound(py).unbind()
     }
 
     pub fn as1d(&self, py: Python) -> Py<PyArray1<f32>> {
@@ -91,7 +91,7 @@ impl CylindricArray {
         for i in 0..self.ycoords.len() {
             out[[i]] = self.array[[self.ycoords[[i]], self.acoords[[i]]]];
         }
-        out.into_pyarray(py).unbind()
+        out.into_pyarray_bound(py).unbind()
     }
 
     /// Convolution on the cylinder surface.
diff --git a/src/regionprops.rs b/src/regionprops.rs
index f2833d07..f9ac2f35 100644
--- a/src/regionprops.rs
+++ b/src/regionprops.rs
@@ -341,7 +341,7 @@ impl RegionProfiler {
                 "std" => self.intensity_std(),
                 _ => return value_error!(format!("Unknown property: {}", prop)),
             };
-            out.insert(prop, Array1::from(vec).into_pyarray(py).unbind());
+            out.insert(prop, Array1::from(vec).into_pyarray_bound(py).unbind());
         }
         Ok(out)
     }
diff --git a/src/viterbi/core.rs b/src/viterbi/core.rs
index 78cb4a61..eaaed108 100644
--- a/src/viterbi/core.rs
+++ b/src/viterbi/core.rs
@@ -109,7 +109,7 @@ impl ViterbiGrid {
         out[0] = pos.z;
         out[1] = pos.y;
         out[2] = pos.x;
-        Ok(out.into_pyarray(py).into())
+        Ok(out.into_pyarray_bound(py).into())
     }
 
     #[pyo3(signature = (dist_min, dist_max, angle_max = None))]
@@ -131,7 +131,7 @@ impl ViterbiGrid {
                 }
             }
         )?;
-        Ok((states.into_pyarray(py).into(), score))
+        Ok((states.into_pyarray_bound(py).into(), score))
     }
 
     #[pyo3(signature = (dist_min, dist_max, coords, origin, angle_max = None))]
@@ -160,7 +160,7 @@ impl ViterbiGrid {
                 }
             }
         )?;
-        Ok((states.into_pyarray(py).into(), score))
+        Ok((states.into_pyarray_bound(py).into(), score))
     }
 }