diff --git a/python/ffsim/linalg/givens.py b/python/ffsim/linalg/givens.py
index fc717a4bf..2c588bfce 100644
--- a/python/ffsim/linalg/givens.py
+++ b/python/ffsim/linalg/givens.py
@@ -72,7 +72,7 @@ def givens_decomposition(
 ) -> tuple[list[tuple[tuple[float, complex], tuple[int, int]]], np.ndarray]:
     """Givens rotation decomposition of a unitary matrix."""
     n, _ = mat.shape
-    current_matrix = mat.astype(complex, copy=False)
+    current_matrix = mat.astype(complex, copy=True)
     left_rotations = []
     right_rotations: list[tuple[tuple[float, complex], tuple[int, int]]] = []
 
diff --git a/tests/gates/orbital_rotation_test.py b/tests/gates/orbital_rotation_test.py
index b20d79871..b4ac71368 100644
--- a/tests/gates/orbital_rotation_test.py
+++ b/tests/gates/orbital_rotation_test.py
@@ -13,6 +13,7 @@
 from __future__ import annotations
 
 import itertools
+from pathlib import Path
 
 import numpy as np
 import pyscf
@@ -89,7 +90,7 @@ def test_apply_orbital_rotation_spin(dtype: type, atol: float):
             np.testing.assert_allclose(result, expected, atol=atol)
 
 
-def test_apply_orbital_rotation_no_side_effects():
+def test_apply_orbital_rotation_no_side_effects_vec():
     """Test applying orbital basis change doesn't modify the original vector."""
     norb = 5
     rng = np.random.default_rng()
@@ -262,9 +263,6 @@ def test_apply_orbital_rotation_compose():
         np.testing.assert_allclose(result, expected_state)
 
 
-@pytest.mark.skip(
-    reason="Flaky test. See https://github.com/qiskit-community/ffsim/issues/97"
-)
 def test_apply_orbital_rotation_nitrogen():
     """Test a special case that was found to cause issues."""
     mol = pyscf.gto.Mole()
@@ -300,3 +298,41 @@ def test_apply_orbital_rotation_nitrogen():
         )
         expected = scipy.sparse.linalg.expm_multiply(op, original_vec, traceA=1)
         np.testing.assert_allclose(result, expected, atol=1e-12)
+
+
+def test_apply_orbital_rotation_no_side_effects_mat():
+    """Test applying orbital rotation doesn't modify the original matrix."""
+    norb = 8
+    nelec = 5, 5
+    vec = ffsim.hartree_fock_state(norb, nelec)
+
+    rng = np.random.default_rng()
+    for _ in range(5):
+        mat = ffsim.random.random_unitary(norb, seed=rng)
+        original_mat = mat.copy()
+        _ = ffsim.apply_orbital_rotation(vec, mat, norb, nelec)
+
+        assert ffsim.linalg.is_unitary(original_mat)
+        assert ffsim.linalg.is_unitary(mat)
+        np.testing.assert_allclose(mat, original_mat, atol=1e-12)
+
+
+def test_apply_orbital_rotation_no_side_effects_special_case():
+    """Test applying orbital rotation doesn't modify the original matrix."""
+    datadir = Path(__file__).parent.parent / "test_data"
+    filepath = datadir / "orbital_rotation-0.npy"
+
+    with open(filepath, "rb") as f:
+        mat = np.load(f)
+    assert ffsim.linalg.is_unitary(mat, atol=1e-12)
+
+    norb = 8
+    nelec = 5, 5
+    vec = ffsim.hartree_fock_state(norb, nelec)
+
+    original_mat = mat.copy()
+    _ = ffsim.apply_orbital_rotation(vec, mat, norb, nelec)
+
+    assert ffsim.linalg.is_unitary(original_mat)
+    assert ffsim.linalg.is_unitary(mat)
+    np.testing.assert_allclose(mat, original_mat, atol=1e-12)
diff --git a/tests/linalg/givens_test.py b/tests/linalg/givens_test.py
index 14b2f6a66..00468ac94 100644
--- a/tests/linalg/givens_test.py
+++ b/tests/linalg/givens_test.py
@@ -12,6 +12,8 @@
 
 from __future__ import annotations
 
+from pathlib import Path
+
 import numpy as np
 from scipy.linalg.lapack import zrot
 
@@ -36,3 +38,34 @@ def test_givens_decomposition():
             reconstructed = reconstructed.T
 
         np.testing.assert_allclose(reconstructed, mat)
+
+
+def test_givens_decomposition_no_side_effects():
+    """Test that the Givens decomposition doesn't modify the original matrix."""
+    norb = 8
+    rng = np.random.default_rng()
+    for _ in range(5):
+        mat = ffsim.random.random_unitary(norb, seed=rng)
+        original_mat = mat.copy()
+        _ = givens_decomposition(mat)
+
+        assert ffsim.linalg.is_unitary(original_mat)
+        assert ffsim.linalg.is_unitary(mat)
+        np.testing.assert_allclose(mat, original_mat, atol=1e-12)
+
+
+def test_givens_decomposition_no_side_effects_special_case():
+    """Test that the Givens decomposition doesn't modify the original matrix."""
+    datadir = Path(__file__).parent.parent / "test_data"
+    filepath = datadir / "orbital_rotation-0.npy"
+
+    with open(filepath, "rb") as f:
+        mat = np.load(f)
+    assert ffsim.linalg.is_unitary(mat, atol=1e-12)
+
+    original_mat = mat.copy()
+    _ = givens_decomposition(mat)
+
+    assert ffsim.linalg.is_unitary(original_mat)
+    assert ffsim.linalg.is_unitary(mat)
+    np.testing.assert_allclose(mat, original_mat, atol=1e-12)
diff --git a/tests/test_data/orbital_rotation-0.npy b/tests/test_data/orbital_rotation-0.npy
new file mode 100644
index 000000000..4a570d497
Binary files /dev/null and b/tests/test_data/orbital_rotation-0.npy differ