diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index ca3dce6b7..b4b0bfd3b 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -12,13 +12,15 @@
 
 <h3>Bug fixes</h3>
 
+* Calls `scipy.integrate.simpson` instead of `scipy.integrate.simps` to be compatible with `scipy==1.14.0` and explicitly states `dtype=object` when creating ragged `numpy` arrays to be compatible with `numpy==1.26.4`. [#748](https://github.com/XanaduAI/strawberryfields/pull/748)
+
 <h3>Documentation</h3>
 
 <h3>Contributors</h3>
 
 This release contains contributions from (in alphabetical order):
 
-Theodor Isacsson
+Luke Helt, Theodor Isacsson
 
 
 # Release 0.23.0 (current release)
diff --git a/requirements.txt b/requirements.txt
index bcd80965d..f7dd47dac 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,11 +1,11 @@
 networkx==2.7.1
 numba==0.55.1
-numpy==1.21.5
+numpy==1.26.4
 plotly==5.6.0
 python-dateutil==2.8.2
 quantum-blackbird==0.4.0
 requests==2.27.1
-scipy==1.8.0
+scipy==1.14.0
 sympy==1.10
 tensorflow==2.8.0
 tensorboard==2.8.0
diff --git a/setup.py b/setup.py
index d302dc8b8..26c28bb65 100644
--- a/setup.py
+++ b/setup.py
@@ -21,11 +21,11 @@
 requirements = [
     "networkx>=2.0",
     "numba",
-    "numpy>=1.17.4",
+    "numpy>=1.17.4,<2",
     "python-dateutil>=2.8.0",
     "quantum-blackbird>=0.4.0",
     "requests>=2.22.0",
-    "scipy>=1.0.0",
+    "scipy>=1.14.0",
     "sympy>=1.5",
     "thewalrus>=0.18.0",
     "toml",
diff --git a/strawberryfields/backends/states.py b/strawberryfields/backends/states.py
index 184532909..f5dd9f36f 100644
--- a/strawberryfields/backends/states.py
+++ b/strawberryfields/backends/states.py
@@ -24,7 +24,7 @@
 from scipy.linalg import block_diag
 from scipy.stats import multivariate_normal
 from scipy.special import factorial
-from scipy.integrate import simps
+from scipy.integrate import simpson
 
 from thewalrus.symplectic import rotation as _R
 from thewalrus.symplectic import xpxp_to_xxpp
@@ -427,7 +427,6 @@ def parity_expectation(self, modes):
         raise NotImplementedError
 
     def p_quad_values(self, mode, xvec, pvec):
-
         r"""Calculates the discretized p-quadrature probability distribution of the specified mode.
 
         Args:
@@ -443,12 +442,11 @@ def p_quad_values(self, mode, xvec, pvec):
         W = self.wigner(mode, xvec, pvec)
         y = []
         for i in range(0, len(pvec)):
-            res = simps(W[i, : len(xvec)], xvec)
+            res = simpson(W[i, : len(xvec)], x=xvec)
             y.append(res)
         return np.array(y)
 
     def x_quad_values(self, mode, xvec, pvec):
-
         r"""Calculates the discretized x-quadrature probability distribution of the specified mode.
 
         Args:
@@ -464,7 +462,7 @@ def x_quad_values(self, mode, xvec, pvec):
         W = self.wigner(mode, xvec, pvec)
         y = []
         for i in range(0, len(xvec)):
-            res = simps(W[: len(pvec), i], pvec)
+            res = simpson(W[: len(pvec), i], x=pvec)
             y.append(res)
         return np.array(y)
 
@@ -1627,7 +1625,6 @@ def wigner(self, mode, xvec, pvec):
 
         wigner = 0
         for i, weight_i in enumerate(weights):
-
             if X.shape == P.shape:
                 arr = np.array([X - means[i, 0], P - means[i, 1]])
                 arr = arr.squeeze()
diff --git a/strawberryfields/io/blackbird_io.py b/strawberryfields/io/blackbird_io.py
index eed23d6bb..777ab9dea 100644
--- a/strawberryfields/io/blackbird_io.py
+++ b/strawberryfields/io/blackbird_io.py
@@ -243,7 +243,10 @@ def to_blackbird(prog: Program, version: str = "1.0") -> blackbird.BlackbirdProg
             }
         )
         bb._var.update(
-            {f"{p.name}": np.array([prog.tdm_params[i]]) for i, p in enumerate(prog.loop_vars)}
+            {
+                f"{p.name}": np.array([prog.tdm_params[i]], dtype=object)
+                for i, p in enumerate(prog.loop_vars)
+            }
         )
 
     return bb
diff --git a/tests/apps/test_data.py b/tests/apps/test_data.py
index d51babf1c..0596784d9 100644
--- a/tests/apps/test_data.py
+++ b/tests/apps/test_data.py
@@ -259,7 +259,8 @@ def dataset(self, datasets):
                 [1.1432025, 1.10163201, 0.07714951, 0.74115172],
                 [0.63240514, 0.99179716, 0.74115172, 0.41063534],
             ],
-        ]
+        ],
+        dtype=object,
     )
 
     @pytest.fixture
diff --git a/tests/apps/train/test_embed.py b/tests/apps/train/test_embed.py
index 53e863436..537a00263 100644
--- a/tests/apps/train/test_embed.py
+++ b/tests/apps/train/test_embed.py
@@ -26,28 +26,30 @@
     [[0.1, 0.2, 0.3], [0.3, 0.1, 0.2], [0.2, 0.3, 0.1]],
     [[0.1, 0.2, 0.3, 0.4], [0.4, 0.1, 0.2, 0.3], [0.3, 0.4, 0.1, 0.2], [0.2, 0.3, 0.4, 0.1]],
 ]
-feats = np.array([np.array(f) for f in feats])
+feats = np.array([np.array(f) for f in feats], dtype=object)
 
 ps = [[1.0, 2.0], [1.0, 2.0, 3.0], [1.0, 2.0, 3.0, 4.0]]
-ps = np.array([np.array(p) for p in ps])
+ps = np.array([np.array(p) for p in ps], dtype=object)
 
 weights_f = np.array(
     [
         np.exp(-np.array([0.5, 0.4])),
         np.exp(-np.array([1.4, 1.1, 1.1])),
         np.exp(-np.array([3.0, 2.4, 2.2, 2.4])),
-    ]
+    ],
+    dtype=object,
 )
 weights = np.array(
     [
         np.exp(-np.array([1.0, 2.0])),
         np.exp(-np.array([1.0, 2.0, 3.0])),
         np.exp(-np.array([1.0, 2.0, 3.0, 4.0])),
-    ]
+    ],
+    dtype=object,
 )
 
-jacobian_f = np.array([np.zeros((d, d)) for d in range(2, 5)])
-jacobian = np.array([np.zeros((d, d)) for d in range(2, 5)])
+jacobian_f = np.array([np.zeros((d, d)) for d in range(2, 5)], dtype=object)
+jacobian = np.array([np.zeros((d, d)) for d in range(2, 5)], dtype=object)
 
 for i in range(3):
     jacobian[i] = -np.diag(weights[i])
diff --git a/tests/backend/test_squeeze_operation.py b/tests/backend/test_squeeze_operation.py
index 3b5517ca5..76973f782 100644
--- a/tests/backend/test_squeeze_operation.py
+++ b/tests/backend/test_squeeze_operation.py
@@ -36,7 +36,7 @@ def matrix_elem(n, r, m):
         return 0.0
 
     if r == 0.0:
-        return np.complex(n == m)  # delta function
+        return complex(n == m)  # delta function
 
     k = np.arange(m % 2, min([m, n]) + 1, 2)
     res = np.sum(
diff --git a/tests/backend/test_states_polyquad.py b/tests/backend/test_states_polyquad.py
index 79c8074df..ccf19a1f4 100644
--- a/tests/backend/test_states_polyquad.py
+++ b/tests/backend/test_states_polyquad.py
@@ -16,7 +16,7 @@
 
 import numpy as np
 from scipy.stats import multivariate_normal
-from scipy.integrate import simps
+from scipy.integrate import simpson
 from scipy.linalg import block_diag
 
 from thewalrus.symplectic import rotation as R
@@ -79,8 +79,8 @@ def _sample(func, correction=0, mu=None, cov=None):
             poly = func(X, P, XP)
             PDF = multivariate_normal.pdf(grid, mu, cov)
 
-            Ex = simps(simps(poly * PDF, P[0]), X.T[0])
-            ExSq = simps(simps(poly**2 * PDF, P[0]), X.T[0])
+            Ex = simpson(simpson(poly * PDF, x=P[0]), x=X.T[0])
+            ExSq = simpson(simpson(poly**2 * PDF, x=P[0]), x=X.T[0])
 
             var = ExSq - Ex**2 + correction
 
diff --git a/tests/frontend/test_post_processing.py b/tests/frontend/test_post_processing.py
index 231ebbbbb..d6a76121f 100644
--- a/tests/frontend/test_post_processing.py
+++ b/tests/frontend/test_post_processing.py
@@ -207,7 +207,7 @@ def test_invalid_samples(self, samples):
     # Arbitrary sequences that are considered to have invalid types during the
     # input checks
     invalid_type_modes_sequences = [
-        np.array([list([0]), list([1, 2])]),
+        np.array([list([0]), list([1, 2])], dtype=object),
     ]
 
     @pytest.mark.parametrize("modes", invalid_type_modes_sequences)