diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 3631b3fdb..d12dbc4d7 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -123,6 +123,15 @@
   method.
   [(#564)](https://github.com/XanaduAI/strawberryfields/pull/564)
 
+* The backend utility module `shared_ops.py` has been removed, with all of its
+  functionality now provided by The Walrus.
+  [(#573)](https://github.com/XanaduAI/strawberryfields/pull/573)
+
+<h3>Breaking changes</h3>
+
+* Removes support for Python 3.6.
+  [(#573)](https://github.com/XanaduAI/strawberryfields/pull/573)
+
 <h3>Bug fixes</h3>
 
 * `Connection` objects now send requests to the platform API at version `0.2.0`
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index cf989d50e..a61d80bd0 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -41,7 +41,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v2
         with:
-          python-version: 3.6
+          python-version: 3.8
 
       - name: Install dependencies
         run: |
diff --git a/.readthedocs.yml b/.readthedocs.yml
index 6daec9bc9..2de688ba5 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -4,7 +4,7 @@ sphinx:
   configuration: doc/conf.py
 
 python:
-    version: 3.6
+    version: 3.8
     install:
         - requirements: doc/requirements.txt
         - method: pip
diff --git a/doc/code/sf_decompositions.rst b/doc/code/sf_decompositions.rst
index cdfddb729..f60182356 100644
--- a/doc/code/sf_decompositions.rst
+++ b/doc/code/sf_decompositions.rst
@@ -15,4 +15,4 @@ sf.decompositions
 .. automodapi:: strawberryfields.decompositions
     :no-heading:
     :include-all-objects:
-    :skip: block_diag, sqrtm, polar, schur, sympmat, changebasis, adj_scaling
+    :skip: block_diag, sqrtm, polar, schur, sympmat, xpxp_to_xxpp, adj_scaling
diff --git a/doc/requirements.txt b/doc/requirements.txt
index 79f8b73e0..fa84097fa 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -4,7 +4,7 @@ ipykernel
 sphinx==2.2.2
 m2r
 networkx>=2.0
-numpy>=1.16.3
+numpy>=1.20
 plotly
 quantum-blackbird
 scipy>=1.0.0
@@ -13,6 +13,5 @@ sphinx-autodoc-typehints==1.10.3
 sphinx-copybutton
 sphinx-automodapi
 sphinxcontrib-bibtex==0.4.2
-tensorflow-tensorboard==0.1.8
-tensorflow==2.0.3
+tensorflow==2.2.0
 toml
diff --git a/requirements.txt b/requirements.txt
index ebb6c641d..467451b9a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,12 +1,12 @@
-numpy>=1.17.4
-scipy==1.4.1
+numpy>=1.20
+scipy
 sympy>=1.5
 tensorflow>=2.0
 tensorboard>=2.0
 networkx>=2.0
 quantum-blackbird>=0.3.0
 python-dateutil==2.8.0
-thewalrus>=0.14.0
+thewalrus>=0.15.0
 toml
 appdirs
 numba>=0.48.0
diff --git a/setup.py b/setup.py
index 49dfb3ef1..c98c80a2b 100644
--- a/setup.py
+++ b/setup.py
@@ -28,7 +28,7 @@
     "networkx>=2.0",
     "quantum-blackbird>=0.3.0",
     "python-dateutil>=2.8.0",
-    "thewalrus>=0.14.0",
+    "thewalrus>=0.15.0",
     "numba",
     "toml",
     "appdirs",
diff --git a/strawberryfields/backends/__init__.py b/strawberryfields/backends/__init__.py
index 04c715f2a..245fe9ede 100644
--- a/strawberryfields/backends/__init__.py
+++ b/strawberryfields/backends/__init__.py
@@ -66,17 +66,6 @@
     BaseGaussian
     BaseBosonic
 
-Utility modules
----------------
-
-The following utility modules are provided for
-backend development.
-
-.. currentmodule:: strawberryfields.backends
-.. autosummary::
-    :toctree: api
-
-    shared_ops
 """
 
 from .base import BaseBackend, BaseFock, BaseGaussian, BaseBosonic, ModeMap
diff --git a/strawberryfields/backends/bosonicbackend/backend.py b/strawberryfields/backends/bosonicbackend/backend.py
index e81be0850..a78219469 100644
--- a/strawberryfields/backends/bosonicbackend/backend.py
+++ b/strawberryfields/backends/bosonicbackend/backend.py
@@ -23,8 +23,9 @@
 from scipy.special import comb
 from scipy.linalg import block_diag
 
+from thewalrus.symplectic import xxpp_to_xpxp
+
 from strawberryfields.backends import BaseBosonic
-from strawberryfields.backends.shared_ops import changebasis
 from strawberryfields.backends.states import BaseBosonicState
 
 from strawberryfields.backends.bosonicbackend.bosoniccircuit import BosonicModes
@@ -377,8 +378,7 @@ def prepare_gaussian_state(self, r, V, modes):
 
         # convert xp-ordering to symmetric ordering
         means = np.vstack([r[:N], r[N:]]).reshape(-1, order="F")
-        C = changebasis(N)
-        cov = C @ V @ C.T
+        cov = xxpp_to_xpxp(V)
 
         self.circuit.from_covmat(cov, modes)
         self.circuit.from_mean(means, modes)
diff --git a/strawberryfields/backends/gaussianbackend/backend.py b/strawberryfields/backends/gaussianbackend/backend.py
index ba9b13acc..bfef2b8b4 100644
--- a/strawberryfields/backends/gaussianbackend/backend.py
+++ b/strawberryfields/backends/gaussianbackend/backend.py
@@ -20,8 +20,6 @@
     concatenate,
     array,
     identity,
-    arctan2,
-    angle,
     sqrt,
     vstack,
     zeros_like,
@@ -29,9 +27,9 @@
     ix_,
 )
 from thewalrus.samples import hafnian_sample_state, torontonian_sample_state
+from thewalrus.symplectic import xxpp_to_xpxp
 
 from strawberryfields.backends import BaseGaussian
-from strawberryfields.backends.shared_ops import changebasis
 from strawberryfields.backends.states import BaseGaussianState
 
 from .gaussiancircuit import GaussianModes
@@ -199,8 +197,7 @@ def prepare_gaussian_state(self, r, V, modes):
 
         # convert xp-ordering to symmetric ordering
         means = vstack([r[:N], r[N:]]).reshape(-1, order="F")
-        C = changebasis(N)
-        cov = C @ V @ C.T
+        cov = xxpp_to_xpxp(V)
 
         self.circuit.fromscovmat(cov, modes)
         self.circuit.fromsmean(means, modes)
diff --git a/strawberryfields/backends/gaussianbackend/gaussiancircuit.py b/strawberryfields/backends/gaussianbackend/gaussiancircuit.py
index 948779f17..894da35ed 100644
--- a/strawberryfields/backends/gaussianbackend/gaussiancircuit.py
+++ b/strawberryfields/backends/gaussianbackend/gaussiancircuit.py
@@ -15,9 +15,9 @@
 # pylint: disable=duplicate-code,attribute-defined-outside-init
 import numpy as np
 from thewalrus.quantum import Xmat
+from thewalrus.symplectic import xxpp_to_xpxp, xpxp_to_xxpp
 
 from . import ops
-from ..shared_ops import changebasis
 
 
 class GaussianModes:
@@ -250,9 +250,6 @@ def scovmatxp(self):
         The order for the canonical operators is :math:`q_1,..,q_n, p_1,...,p_n`.
         This differs from the ordering used in [1] which is :math:`q_1,p_1,q_2,p_2,...,q_n,p_n`
         Note that one ordering can be obtained from the other by using a permutation matrix.
-
-        Said permutation matrix is implemented in the function changebasis(n) where n is
-        the number of modes.
         """
         mm11 = (
             self.nmat
@@ -288,9 +285,6 @@ def smeanxp(self):
         The order for the canonical operators is :math:`q_1, \ldots, q_n, p_1, \ldots, p_n`.
         This differs from the ordering used in [1] which is :math:`q_1, p_1, q_2, p_2, \ldots, q_n, p_n`.
         Note that one ordering can be obtained from the other by using a permutation matrix.
-
-        Said permutation matrix is implemented in the function changebasis(n) where n is
-        the number of modes.
         """
         nmodes = self.nlen
         r = np.empty(2 * nmodes)
@@ -300,8 +294,7 @@ def smeanxp(self):
 
     def scovmat(self):
         """Constructs and returns the symmetric ordered covariance matrix as defined in [1]"""
-        rotmat = changebasis(self.nlen)
-        return np.dot(np.dot(rotmat, self.scovmatxp()), np.transpose(rotmat))
+        return xxpp_to_xpxp(self.scovmatxp())
 
     def smean(self):
         r"""the symmetric mean $[q_1,p_1,q_2,p_2,...,q_n,p_n]$"""
@@ -347,8 +340,7 @@ def fromscovmat(self, V, modes=None):
                 raise ValueError("Covariance matrix is larger than the number of subsystems.")
 
         # convert to xp ordering
-        rotmat = changebasis(n)
-        VV = np.dot(np.dot(np.transpose(rotmat), V), rotmat)
+        VV = xpxp_to_xxpp(V)
 
         A = VV[0:n, 0:n]
         B = VV[0:n, n : 2 * n]
diff --git a/strawberryfields/backends/shared_ops.py b/strawberryfields/backends/shared_ops.py
deleted file mode 100644
index 959729cb4..000000000
--- a/strawberryfields/backends/shared_ops.py
+++ /dev/null
@@ -1,41 +0,0 @@
-# Copyright 2019 Xanadu Quantum Technologies Inc.
-
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-
-#     http://www.apache.org/licenses/LICENSE-2.0
-
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-"""Common shared operations that can be used by backends"""
-
-import functools
-
-import numpy as np
-
-# ================================+
-# Phase space shared operations  |
-# ================================+
-
-
-@functools.lru_cache()
-def changebasis(n):
-    r"""Change of basis matrix between the two Gaussian representation orderings.
-
-    This is the matrix necessary to transform covariances matrices written
-    in the (x_1,...,x_n,p_1,...,p_n) to the (x_1,p_1,...,x_n,p_n) ordering
-
-    Args:
-        n (int): number of modes
-    Returns:
-        array: :math:`2n\times 2n` matrix
-    """
-    m = np.zeros((2 * n, 2 * n))
-    for i in range(n):
-        m[2 * i, i] = 1
-        m[2 * i + 1, i + n] = 1
-    return m
diff --git a/strawberryfields/backends/states.py b/strawberryfields/backends/states.py
index 7367c1eaf..c344847c7 100644
--- a/strawberryfields/backends/states.py
+++ b/strawberryfields/backends/states.py
@@ -27,13 +27,12 @@
 from scipy.integrate import simps
 
 from thewalrus.symplectic import rotation as _R
+from thewalrus.symplectic import xpxp_to_xxpp
 
 import thewalrus.quantum as twq
 
 import strawberryfields as sf
 
-from .shared_ops import changebasis
-
 indices = string.ascii_lowercase
 
 
@@ -1257,8 +1256,7 @@ def poly_quad_expectation(self, A, d=None, k=0, phi=0, **kwargs):
         if phi != 0:
             # rotate all modes of the covariance matrix and vector of means
             R = _R(phi)
-            C = changebasis(self._modes)
-            rot = C.T @ block_diag(*([R] * self._modes)) @ C
+            rot = xpxp_to_xxpp(block_diag(*([R] * self._modes)))
 
             mu = rot.T @ mu
             cov = rot.T @ cov @ rot
diff --git a/strawberryfields/decompositions.py b/strawberryfields/decompositions.py
index 9f8a2897b..f591609f5 100644
--- a/strawberryfields/decompositions.py
+++ b/strawberryfields/decompositions.py
@@ -21,9 +21,7 @@
 import numpy as np
 from scipy.linalg import block_diag, sqrtm, polar, schur
 from thewalrus.quantum import adj_scaling
-from thewalrus.symplectic import sympmat
-
-from .backends.shared_ops import changebasis
+from thewalrus.symplectic import sympmat, xpxp_to_xxpp
 
 
 def takagi(N, tol=1e-13, rounding=13):
@@ -670,7 +668,6 @@ def williamson(V, tol=1e-11):
 
     n = n // 2
     omega = sympmat(n)
-    rotmat = changebasis(n)
     vals = np.linalg.eigvalsh(V)
 
     for val in vals:
@@ -698,8 +695,9 @@ def williamson(V, tol=1e-11):
     p = block_diag(*seq)
     Kt = K @ p
     s1t = p @ s1 @ p
-    dd = np.transpose(rotmat) @ s1t @ rotmat
-    Ktt = Kt @ rotmat
+    dd = xpxp_to_xxpp(s1t)
+    perm_indices = xpxp_to_xxpp(np.arange(2 * n))
+    Ktt = Kt[:, perm_indices]
     Db = np.diag([1 / dd[i, i + n] for i in range(n)] + [1 / dd[i, i + n] for i in range(n)])
     S = Mm12 @ Ktt @ sqrtm(Db)
     return Db, np.linalg.inv(S).T
diff --git a/strawberryfields/ops.py b/strawberryfields/ops.py
index bb47000a2..9a02e0f8e 100644
--- a/strawberryfields/ops.py
+++ b/strawberryfields/ops.py
@@ -25,11 +25,12 @@
 from scipy.linalg import block_diag
 import scipy.special as ssp
 
+from thewalrus.symplectic import xxpp_to_xpxp
+
 import strawberryfields as sf
 import strawberryfields.program_utils as pu
 import strawberryfields.decompositions as dec
 from .backends.states import BaseFockState, BaseGaussianState, BaseBosonicState
-from .backends.shared_ops import changebasis
 from .program_utils import Command, RegRef, MergeFailure
 from .parameters import (
     par_regref_deps,
@@ -2842,7 +2843,7 @@ def _decompose(self, reg, **kwargs):
         D = np.diag(V)
         is_diag = np.all(V == np.diag(D))
 
-        BD = changebasis(self.ns) @ V @ changebasis(self.ns).T
+        BD = xxpp_to_xpxp(V)
         BD_modes = [BD[i * 2 : (i + 1) * 2, i * 2 : (i + 1) * 2] for i in range(BD.shape[0] // 2)]
         is_block_diag = (not is_diag) and np.all(BD == block_diag(*BD_modes))
 
diff --git a/tests/backend/test_states_polyquad.py b/tests/backend/test_states_polyquad.py
index 065cbc4d3..4bec140a6 100644
--- a/tests/backend/test_states_polyquad.py
+++ b/tests/backend/test_states_polyquad.py
@@ -20,10 +20,10 @@
 from scipy.linalg import block_diag
 
 from thewalrus.symplectic import rotation as R
+from thewalrus.symplectic import xpxp_to_xxpp
 
 from strawberryfields import backends
 from strawberryfields import utils
-from strawberryfields.backends.shared_ops import changebasis
 
 # some tests require a higher cutoff for accuracy
 CUTOFF = 12
@@ -453,9 +453,8 @@ def test_three_mode_arbitrary(self, setup_backend, pure, hbar, tol):
         S1 = block_diag(BS, np.identity(2))
         S2 = block_diag(np.identity(2), BS)
 
-        C = changebasis(3)
-        mu = C.T @ S2 @ S1 @ mu
-        cov = C.T @ S2 @ S1 @ cov @ S1.T @ S2.T @ C
+        mu = xpxp_to_xxpp(S2 @ S1 @ mu)
+        cov = xpxp_to_xxpp(S2 @ S1 @ cov @ S1.T @ S2.T)
 
         modes = list(np.arange(6).reshape(2, -1).T)
 
diff --git a/tests/frontend/test_shared_ops.py b/tests/frontend/test_shared_ops.py
deleted file mode 100644
index 952a8468f..000000000
--- a/tests/frontend/test_shared_ops.py
+++ /dev/null
@@ -1,49 +0,0 @@
-# Copyright 2019 Xanadu Quantum Technologies Inc.
-
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-
-#     http://www.apache.org/licenses/LICENSE-2.0
-
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-r"""Unit tests for the Strawberry Fields shared_ops module"""
-import pytest
-
-pytestmark = pytest.mark.frontend
-
-import numpy as np
-
-import strawberryfields.backends.shared_ops as so
-
-class TestPhaseSpaceFunctions:
-    """Tests for the shared phase space operations"""
-
-    def test_means_changebasis(self):
-        """Test the change of basis function applied to vectors. This function
-        converts from xp to symmetric ordering, and vice versa."""
-        C = so.changebasis(3)
-        means_xp = [1, 2, 3, 4, 5, 6]
-        means_symmetric = [1, 4, 2, 5, 3, 6]
-
-        assert np.all(C @ means_xp == means_symmetric)
-        assert np.all(C.T @ means_symmetric == means_xp)
-
-    def test_cov_changebasis(self):
-        """Test the change of basis function applied to matrices. This function
-        converts from xp to symmetric ordering, and vice versa."""
-        C = so.changebasis(2)
-        cov_xp = np.array(
-            [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11], [12, 13, 14, 15]]
-        )
-
-        cov_symmetric = np.array(
-            [[0, 2, 1, 3], [8, 10, 9, 11], [4, 6, 5, 7], [12, 14, 13, 15]]
-        )
-
-        assert np.all(C @ cov_xp @ C.T == cov_symmetric)
-        assert np.all(C.T @ cov_symmetric @ C == cov_xp)
diff --git a/tests/integration/test_decompositions_integration.py b/tests/integration/test_decompositions_integration.py
index 7c68b7cb1..232f70e26 100644
--- a/tests/integration/test_decompositions_integration.py
+++ b/tests/integration/test_decompositions_integration.py
@@ -18,6 +18,7 @@
 from scipy.linalg import qr, block_diag
 
 from thewalrus.symplectic import rotation as rot
+from thewalrus.symplectic import xpxp_to_xxpp, xxpp_to_xpxp
 
 import strawberryfields as sf
 from strawberryfields import decompositions as dec
@@ -28,7 +29,6 @@
     squeezed_state,
 )
 from strawberryfields import ops
-from strawberryfields.backends.shared_ops import changebasis
 
 from strawberryfields.utils import random_interferometer as haar_measure
 # make the test file deterministic
@@ -292,8 +292,7 @@ def test_rotated_squeezed(self, setup_eng, hbar, tol):
         r = 0.1
         phi = 0.2312
         v1 = (hbar / 2) * np.diag([np.exp(-r), np.exp(r)])
-        A = changebasis(3)
-        cov = A.T @ block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3) @ A
+        cov = xpxp_to_xxpp(block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3))
 
         with prog.context as q:
             ops.Gaussian(cov, decomp=False) | q
@@ -301,19 +300,6 @@ def test_rotated_squeezed(self, setup_eng, hbar, tol):
         state = eng.run(prog).state
         assert np.allclose(state.cov(), cov, atol=tol)
 
-def from_xp(num_modes):
-    r"""Provides array of indices to order quadratures as (x1,p1,...,xn,pn)
-    starting from all x followed by all p i.e., (x1,...,xn,p1,..., pn).
-
-    Args:
-        num_modes (int): number of modes
-
-    Returns:
-        list: quadrature ordering for (x1,p1,...,xn,pn)
-    """
-    perm_inds_list = [(i, i + num_modes) for i in range(num_modes)]
-    perm_inds = [a for tup in perm_inds_list for a in tup]
-    return perm_inds
 
 @pytest.mark.backends("bosonic")
 class TestBosonicBackendPrepareState:
@@ -328,7 +314,7 @@ def test_vacuum(self, setup_eng, hbar, tol):
 
         state = eng.run(prog).state
         
-        indices = from_xp(3)
+        indices = xxpp_to_xpxp(np.arange(2 * 3))
         cov = cov[:,indices][indices,:]
         assert np.allclose(state.covs(), np.expand_dims(cov,axis=0), atol=tol)
         assert np.all(state.means() == np.zeros((1,6)))
@@ -345,7 +331,7 @@ def test_squeezed(self, setup_eng, hbar, tol):
 
         state = eng.run(prog).state
             
-        indices = from_xp(3)
+        indices = xxpp_to_xpxp(np.arange(2 * 3))
         cov = cov[:,indices][indices,:]
         assert np.allclose(state.covs(), np.expand_dims(cov,axis=0), atol=tol)
 
@@ -360,7 +346,7 @@ def test_displaced_squeezed(self, setup_eng, hbar, tol):
 
         state = eng.run(prog).state
             
-        indices = from_xp(3)
+        indices = xxpp_to_xpxp(np.arange(2 * 3))
         cov = cov[:,indices][indices,:]
         means = means[indices]
         assert np.allclose(state.covs(), np.expand_dims(cov,axis=0), atol=tol)
@@ -376,7 +362,7 @@ def test_thermal(self, setup_eng, hbar, tol):
 
         state = eng.run(prog).state
             
-        indices = from_xp(3)
+        indices = xxpp_to_xpxp(np.arange(2 * 3))
         cov = cov[:,indices][indices,:]
         assert np.allclose(state.covs(), np.expand_dims(cov,axis=0), atol=tol)
 
@@ -387,15 +373,14 @@ def test_rotated_squeezed(self, setup_eng, hbar, tol):
         r = 0.1
         phi = 0.2312
         v1 = (hbar / 2) * np.diag([np.exp(-r), np.exp(r)])
-        A = changebasis(3)
-        cov = A.T @ block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3) @ A
+        cov = xpxp_to_xxpp(block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3))
 
         with prog.context as q:
             ops.Gaussian(cov, decomp=False) | q
 
         state = eng.run(prog).state
             
-        indices = from_xp(3)
+        indices = xxpp_to_xpxp(np.arange(2 * 3))
         cov = cov[:,indices][indices,:]
         assert np.allclose(state.covs(), np.expand_dims(cov,axis=0), atol=tol)
 
@@ -462,8 +447,7 @@ def test_rotated_squeezed(self, setup_eng, hbar, tol):
         r = 0.1
         phi = 0.2312
         v1 = (hbar / 2) * np.diag([np.exp(-r), np.exp(r)])
-        A = changebasis(3)
-        cov = A.T @ block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3) @ A
+        cov = xpxp_to_xxpp(block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3))
 
         with prog.context as q:
             ops.Gaussian(cov) | q
@@ -534,8 +518,7 @@ def test_rotated_squeezed(self, setup_eng, cutoff, hbar, tol):
         in_state = squeezed_state(r, phi, basis="fock", fock_dim=cutoff)
 
         v1 = (hbar / 2) * np.diag([np.exp(-2 * r), np.exp(2 * r)])
-        A = changebasis(3)
-        cov = A.T @ block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3) @ A
+        cov = xpxp_to_xxpp(block_diag(*[rot(phi) @ v1 @ rot(phi).T] * 3))
 
         with prog.context as q:
             ops.Gaussian(cov) | q
@@ -615,7 +598,7 @@ def test_CXgate(self, setup_eng, pure, hbar, tol):
             assert np.allclose(state.means(), rexpected, atol=tol, rtol=0)
             
         elif eng.backend_name == "bosonic":
-            indices = from_xp(2)
+            indices = xxpp_to_xpxp(np.arange(2 * 2))
             Vexpected = Vexpected[:,indices][indices,:]
             rexpected = rexpected[indices]
             assert np.allclose(state.covs(), np.expand_dims(Vexpected,axis=0), atol=tol, rtol=0)
@@ -653,7 +636,7 @@ def test_CZgate(self, setup_eng, pure, hbar, tol):
             assert np.allclose(state.means(), rexpected, atol=tol, rtol=0)
             
         elif eng.backend_name == "bosonic":
-            indices = from_xp(2)
+            indices = xxpp_to_xpxp(np.arange(2 * 2))
             Vexpected = Vexpected[:,indices][indices,:]
             rexpected = rexpected[indices]
             assert np.allclose(state.covs(), np.expand_dims(Vexpected,axis=0), atol=tol, rtol=0)