diff --git a/docs/tutorials/dynamics_backend.rst b/docs/tutorials/dynamics_backend.rst
index aca14ba86..371b365ac 100644
--- a/docs/tutorials/dynamics_backend.rst
+++ b/docs/tutorials/dynamics_backend.rst
@@ -318,8 +318,9 @@ To enable running of the single qubit experiments, we add the following to the `
   backend. 
 - Add definitions of ``RZ`` gates as phase shifts. These instructions control the phase of the drive
   channels, as well as any control channels acting on a given qubit.
-- Add a ``CX`` gate which applies to all qubits. While this tutorial will not be utilizing it, this
-  ensures that validation steps checking that the device is fully connected will pass.
+- Add a ``CX`` gate between qubits :math:`(0, 1)` and :math:`(1, 0)`. While this tutorial will not 
+  be utilizing it, this ensures that validation steps checking that the device is fully connected 
+  will pass.
 
 .. jupyter-execute::
 
@@ -336,7 +337,7 @@ To enable running of the single qubit experiments, we add the following to the `
     target.add_instruction(XGate(), properties={(0,): None, (1,): None})
     target.add_instruction(SXGate(), properties={(0,): None, (1,): None})
 
-    target.add_instruction(CXGate())
+    target.add_instruction(CXGate(), properties={(0, 1): None, (1, 0): None})
     
     # Add RZ instruction as phase shift for drag cal
     phi = Parameter("phi")
@@ -472,11 +473,6 @@ values for the single qubit gates calibrated above.
 
     from qiskit_experiments.library import CrossResonanceHamiltonian
 
-    backend.target.add_instruction(
-        instruction=CrossResonanceHamiltonian.CRPulseGate(width=Parameter("width")), 
-        properties={(0, 1): None, (1, 0): None}
-    )
-
     cr_ham_experiment = CrossResonanceHamiltonian(
         qubits=(0, 1), 
         flat_top_widths=np.linspace(0, 5000, 17), 
diff --git a/qiskit_dynamics/arraylias/arraylias_state.py b/qiskit_dynamics/arraylias/arraylias_state.py
new file mode 100644
index 000000000..92bb29912
--- /dev/null
+++ b/qiskit_dynamics/arraylias/arraylias_state.py
@@ -0,0 +1,52 @@
+from typing import Union
+from collections.abc import Iterable
+from arraylias import numpy_alias
+from scipy.sparse import spmatrix
+from qiskit.quantum_info.operators import Operator
+from .asarray import register_to_asarray
+from .to_dense import register_to_dense
+from .to_numeric_matrix_type import register_to_numeric_matrix_type
+from .to_sparse import register_to_sparse
+from .matmul import register_matmul
+from .rmatmul import register_rmatmul
+from .multiply import register_multiply
+
+from ..array import Array
+
+
+DYNAMICS_ALIAS = numpy_alias()
+
+# Set qiskit_dynamics.array.Array to be dispatched to numpy
+DYNAMICS_ALIAS.register_type(Array, "numpy")
+
+# register required custom versions of functions for sparse type here
+DYNAMICS_ALIAS.register_type(spmatrix, lib="scipy_sparse")
+
+try:
+    from jax.experimental.sparse import BCOO
+
+    # register required custom versions of functions for BCOO type here
+    DYNAMICS_ALIAS.register_type(BCOO, lib="jax_sparse")
+except ImportError:
+    pass
+
+# register required custom versions of functions for Operator type here
+DYNAMICS_ALIAS.register_type(Operator, lib="operator")
+
+# register required custom versions of functions for List type here
+# need to discuss registering List type because the coverage of List is too broad.
+DYNAMICS_ALIAS.register_type(list, lib="list")
+
+
+register_to_asarray(alias=DYNAMICS_ALIAS)
+register_to_dense(alias=DYNAMICS_ALIAS)
+register_to_numeric_matrix_type(alias=DYNAMICS_ALIAS)
+register_to_sparse(alias=DYNAMICS_ALIAS)
+register_matmul(alias=DYNAMICS_ALIAS)
+register_multiply(alias=DYNAMICS_ALIAS)
+register_rmatmul(alias=DYNAMICS_ALIAS)
+
+DYNAMICS_NUMPY = DYNAMICS_ALIAS()
+
+
+ArrayLike = Union[DYNAMICS_ALIAS.registered_types()]
diff --git a/qiskit_dynamics/arraylias/asarray.py b/qiskit_dynamics/arraylias/asarray.py
new file mode 100644
index 000000000..838e755be
--- /dev/null
+++ b/qiskit_dynamics/arraylias/asarray.py
@@ -0,0 +1,32 @@
+import numpy as np
+from scipy.sparse import csr_matrix, issparse
+
+
+def register_to_asarray(alias):
+    @alias.register_function(lib="scipy_sparse", path="asarray")
+    def _(arr):
+        if issparse(arr):
+            return arr
+        return csr_matrix(arr)
+
+    @alias.register_function(lib="list", path="asarray")
+    def _(arr):
+        if len(arr) == 0 or isinstance(arr[0], (list, tuple)):
+            return np.asarray(arr)
+        return alias(like=arr[0]).asarray([alias().asarray(sub_arr) for sub_arr in arr])
+
+    @alias.register_fallback(path="asarray")
+    def _(arr):
+        return np.asarray(arr)
+
+    try:
+        from jax.experimental.sparse import BCOO
+
+        @alias.register_function(lib="jax_sparse", path="asarray")
+        def _(arr):
+            if type(arr).__name__ == "BCOO":
+                return arr
+            return BCOO.fromdense(arr)
+
+    except ImportError:
+        pass
diff --git a/qiskit_dynamics/arraylias/matmul.py b/qiskit_dynamics/arraylias/matmul.py
new file mode 100644
index 000000000..bd2d87f3e
--- /dev/null
+++ b/qiskit_dynamics/arraylias/matmul.py
@@ -0,0 +1,17 @@
+def register_matmul(alias):
+    @alias.register_function(lib="scipy_sparse", path="matmul")
+    def _(x, y):
+        return x * y
+
+    try:
+        from jax.experimental import sparse as jsparse
+        import jax.numpy as jnp
+
+        jsparse_matmul = jsparse.sparsify(jnp.matmul)
+
+        @alias.register_function(lib="jax_sparse", path="matmul")
+        def _(x, y):
+            return jsparse_matmul(x, y)
+
+    except:
+        pass
diff --git a/qiskit_dynamics/arraylias/multiply.py b/qiskit_dynamics/arraylias/multiply.py
new file mode 100644
index 000000000..bc671bcb1
--- /dev/null
+++ b/qiskit_dynamics/arraylias/multiply.py
@@ -0,0 +1,17 @@
+def register_multiply(alias):
+    @alias.register_function(lib="scipy_sparse", path="multiply")
+    def _(x, y):
+        return x.multiply(y)
+
+    try:
+        from jax.experimental import sparse as jsparse
+        import jax.numpy as jnp
+
+        jsparse_multiply = jsparse.sparsify(jnp.multiply)
+
+        @alias.register_function(lib="jax_sparse", path="multiply")
+        def _(x, y):
+            return jsparse_multiply(x, y)
+
+    except:
+        pass
diff --git a/qiskit_dynamics/arraylias/rmatmul.py b/qiskit_dynamics/arraylias/rmatmul.py
new file mode 100644
index 000000000..cc90e4ada
--- /dev/null
+++ b/qiskit_dynamics/arraylias/rmatmul.py
@@ -0,0 +1,28 @@
+import numpy as np
+
+
+def register_rmatmul(alias):
+    @alias.register_function(lib="numpy", path="rmatmul")
+    def _(x, y):
+        return np.matmul(y, x)
+
+    @alias.register_function(lib="scipy_sparse", path="rmatmul")
+    def _(x, y):
+        return y * x
+
+    try:
+        from jax.experimental import sparse as jsparse
+        import jax.numpy as jnp
+
+        jsparse_matmul = jsparse.sparsify(jnp.matmul)
+
+        @alias.register_function(lib="jax", path="rmatmul")
+        def _(x, y):
+            return jnp.matmul(y, x)
+
+        @alias.register_function(lib="jax_sparse", path="rmatmul")
+        def _(x, y):
+            return jsparse_matmul(y, x)
+
+    except:
+        pass
diff --git a/qiskit_dynamics/arraylias/to_dense.py b/qiskit_dynamics/arraylias/to_dense.py
new file mode 100644
index 000000000..9953f8cc3
--- /dev/null
+++ b/qiskit_dynamics/arraylias/to_dense.py
@@ -0,0 +1,38 @@
+import numpy as np
+
+
+def register_to_dense(alias):
+    @alias.register_default(path="to_dense")
+    def _(op):
+        if op is None:
+            return None
+        return op
+
+    @alias.register_function(lib="numpy", path="to_dense")
+    def _(op):
+        return op
+
+    try:
+
+        @alias.register_function(lib="jax", path="to_dense")
+        def _(op):
+            return op
+
+        @alias.register_function(lib="jax_sparse", path="to_dense")
+        def _(op):
+            return op.todense()
+
+    except:
+        pass
+
+    @alias.register_function(lib="scipy_sparse", path="to_dense")
+    def _(op):
+        return op.toarray()
+
+    @alias.register_fallback(path="to_dense")
+    def _(op):
+        return np.asarray(op)
+
+    @alias.register_function(lib="list", path="to_dense")
+    def _(op):
+        return alias().asarray([alias().to_dense(sub_op) for sub_op in op])
diff --git a/qiskit_dynamics/arraylias/to_numeric_matrix_type.py b/qiskit_dynamics/arraylias/to_numeric_matrix_type.py
new file mode 100644
index 000000000..816bfc189
--- /dev/null
+++ b/qiskit_dynamics/arraylias/to_numeric_matrix_type.py
@@ -0,0 +1,47 @@
+from scipy.sparse import spmatrix
+from qiskit_dynamics.type_utils import isinstance_qutip_qobj
+
+
+def register_to_numeric_matrix_type(alias):
+    @alias.register_default(path="to_numeric_matrix_type")
+    def _(op):
+        if op is None:
+            return None
+        if isinstance_qutip_qobj(op):
+            return alias().to_sparse(op.data)
+        return op
+
+    @alias.register_function(lib="numpy", path="to_numeric_matrix_type")
+    def _(op):
+        return op
+
+    @alias.register_function(lib="operator", path="to_numeric_matrix_type")
+    def _(op):
+        return alias().to_dense(op)
+
+    try:
+
+        @alias.register_function(lib="jax", path="to_numeric_matrix_type")
+        def _(op):
+            return op
+
+        @alias.register_function(lib="jax_sparse", path="to_numeric_matrix_type")
+        def _(op):
+            return op
+
+    except:
+        pass
+
+    @alias.register_function(lib="scipy_sparse", path="to_numeric_matrix_type")
+    def _(op):
+        return op
+
+    @alias.register_function(lib="list", path="to_numeric_matrix_type")
+    def _(op):
+        if isinstance(op[0], spmatrix) or isinstance_qutip_qobj(op[0]):
+            return [alias().to_sparse(sub_op) for sub_op in op]
+        return alias().asarray([alias().to_dense(sub_op) for sub_op in op])
+
+    @alias.register_fallback(path="to_numeric_matrix_type")
+    def _(op):
+        return op
diff --git a/qiskit_dynamics/arraylias/to_sparse.py b/qiskit_dynamics/arraylias/to_sparse.py
new file mode 100644
index 000000000..cee9f5098
--- /dev/null
+++ b/qiskit_dynamics/arraylias/to_sparse.py
@@ -0,0 +1,52 @@
+import numpy as np
+from scipy.sparse import csr_matrix
+from qiskit_dynamics.type_utils import isinstance_qutip_qobj
+
+
+def register_to_sparse(alias):
+    @alias.register_default(path="to_sparse")
+    def _(op):
+        if op is None:
+            return None
+        if isinstance_qutip_qobj(op):
+            return op.data
+        return op
+
+    @alias.register_function(lib="numpy", path="to_sparse")
+    def _(op):
+        if op.ndim < 3:
+            return csr_matrix(op)
+        return np.array([csr_matrix(sub_op) for sub_op in op])
+
+    try:
+        from jax.experimental.sparse import BCOO
+
+        @alias.register_function(lib="jax", path="to_sparse")
+        def _(op):
+            return BCOO.fromdense(op)
+
+        @alias.register_function(lib="jax_sparse", path="to_sparse")
+        def _(op):
+            return op
+
+    except ImportError:
+        pass
+
+    @alias.register_function(lib="scipy_sparse", path="to_sparse")
+    def _(op):
+        return op
+
+    @alias.register_fallback(path="to_sparse")
+    def _(op):
+        return csr_matrix(op)
+
+    @alias.register_function(lib="list", path="to_sparse")
+    def _(op):
+        try:
+            import jax.numpy as jnp
+
+            if isinstance(op[0], jnp.ndarray):
+                return BCOO.fromdense(op)
+        except ImportError:
+            return np.array([csr_matrix(sub_op) for sub_op in op])
+        return np.array([csr_matrix(sub_op) for sub_op in op])
diff --git a/qiskit_dynamics/arraylias_state.py b/qiskit_dynamics/arraylias_state.py
deleted file mode 100644
index 86ea48b4c..000000000
--- a/qiskit_dynamics/arraylias_state.py
+++ /dev/null
@@ -1,29 +0,0 @@
-# -*- coding: utf-8 -*-
-
-# This code is part of Qiskit.
-#
-# (C) Copyright IBM 2023.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-
-"""Configure custom instance of numpy alias for Dynamics."""
-
-from typing import Union
-from arraylias import numpy_alias
-from .array import Array
-
-
-DYNAMICS_ALIAS = numpy_alias()
-
-# Set qiskit_dynamics.array.Array to be dispatched to numpy
-DYNAMICS_ALIAS.register_type(Array, "numpy")
-
-DYNAMICS_NUMPY = DYNAMICS_ALIAS()
-
-ArrayLike = Union[DYNAMICS_ALIAS.registered_types()]
diff --git a/qiskit_dynamics/models/rotating_frame.py b/qiskit_dynamics/models/rotating_frame.py
index 4d4a4a097..71f63c920 100644
--- a/qiskit_dynamics/models/rotating_frame.py
+++ b/qiskit_dynamics/models/rotating_frame.py
@@ -17,13 +17,14 @@
 
 from typing import Union, List, Optional
 import numpy as np
-from scipy.sparse import issparse, csr_matrix
+from scipy.sparse import issparse
 
 from qiskit import QiskitError
 from qiskit.quantum_info.operators import Operator
 from qiskit.quantum_info.operators.predicates import is_hermitian_matrix
-from qiskit_dynamics.array import Array
-from qiskit_dynamics.type_utils import to_array, to_BCOO, to_numeric_matrix_type
+from qiskit_dynamics.arraylias.arraylias_state import ArrayLike
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_ALIAS
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_NUMPY as unp
 
 try:
     import jax.numpy as jnp
@@ -58,7 +59,7 @@ class RotatingFrame:
     """
 
     def __init__(
-        self, frame_operator: Union[Array, Operator], atol: float = 1e-10, rtol: float = 1e-10
+        self, frame_operator: Union[ArrayLike, Operator], atol: float = 1e-10, rtol: float = 1e-10
     ):
         """Initialize with a frame operator.
 
@@ -73,7 +74,7 @@ def __init__(
             frame_operator = frame_operator.frame_operator
 
         self._frame_operator = frame_operator
-        frame_operator = to_array(frame_operator)
+        frame_operator = unp.to_dense(frame_operator)
 
         if frame_operator is None:
             self._dim = None
@@ -86,7 +87,7 @@ def __init__(
             # if Hermitian convert to anti-Hermitian
             frame_operator = _enforce_anti_herm(frame_operator, atol=atol, rtol=rtol)
 
-            self._frame_diag = Array(frame_operator)
+            self._frame_diag = unp.asarray(frame_operator)
             self._frame_basis = None
             self._frame_basis_adjoint = None
             self._dim = len(self._frame_diag)
@@ -97,10 +98,12 @@ def __init__(
             frame_operator = _enforce_anti_herm(frame_operator, atol=atol, rtol=rtol)
 
             # diagonalize with eigh, utilizing assumption of anti-hermiticity
-            frame_diag, frame_basis = np.linalg.eigh(1j * frame_operator)
 
-            self._frame_diag = Array(-1j * frame_diag)
-            self._frame_basis = Array(frame_basis)
+            # assume a type of frame_operator is np.array or jnp.array
+            frame_diag, frame_basis = unp.linalg.eigh(1j * frame_operator)
+
+            self._frame_diag = -1j * frame_diag
+            self._frame_basis = frame_basis
             self._frame_basis_adjoint = frame_basis.conj().transpose()
             self._dim = len(self._frame_diag)
 
@@ -114,50 +117,50 @@ def dim(self) -> int:
         return self._dim
 
     @property
-    def frame_operator(self) -> Array:
+    def frame_operator(self) -> ArrayLike:
         """The original frame operator."""
         return self._frame_operator
 
     @property
-    def frame_diag(self) -> Array:
+    def frame_diag(self) -> ArrayLike:
         """The diagonal of the frame operator."""
         return self._frame_diag
 
     @property
-    def frame_basis(self) -> Array:
+    def frame_basis(self) -> ArrayLike:
         """The array containing diagonalizing unitary."""
         return self._frame_basis
 
     @property
-    def frame_basis_adjoint(self) -> Array:
+    def frame_basis_adjoint(self) -> ArrayLike:
         """The adjoint of the diagonalizing unitary."""
         return self._frame_basis_adjoint
 
-    def state_into_frame_basis(self, y: Array) -> Array:
+    def state_into_frame_basis(self, y: ArrayLike) -> ArrayLike:
         r"""Take a state into the frame basis, i.e. return ``self.frame_basis_adjoint @ y``.
 
         Args:
             y: The state.
 
         Returns:
-            Array: The state in the frame basis.
+            ArrayLike: The state in the frame basis.
         """
-        y = to_numeric_matrix_type(y)
+        y = unp.to_numeric_matrix_type(y)
         if self.frame_basis_adjoint is None:
             return y
 
         return self.frame_basis_adjoint @ y
 
-    def state_out_of_frame_basis(self, y: Array) -> Array:
+    def state_out_of_frame_basis(self, y: ArrayLike) -> ArrayLike:
         r"""Take a state out of the frame basis, i.e. ``return self.frame_basis @ y``.
 
         Args:
             y: The state.
 
         Returns:
-            Array: The state in the frame basis.
+            ArrayLike: The state in the frame basis.
         """
-        y = to_numeric_matrix_type(y)
+        y = unp.to_numeric_matrix_type(y)
         if self.frame_basis is None:
             return y
 
@@ -165,9 +168,9 @@ def state_out_of_frame_basis(self, y: Array) -> Array:
 
     def operator_into_frame_basis(
         self,
-        op: Union[Operator, List[Operator], Array, csr_matrix, None],
+        op: Union[Operator, List[Operator], ArrayLike, None],
         convert_type: bool = True,
-    ) -> Array:
+    ) -> ArrayLike:
         r"""Take an operator into the frame basis, i.e. return
         ``self.frame_basis_adjoint @ A @ self.frame_basis``
 
@@ -178,27 +181,24 @@ def operator_into_frame_basis(
                           ``op`` is a handled input type.
 
         Returns:
-            Array: The operator in the frame basis.
+            ArrayLike: The operator in the frame basis.
         """
         if convert_type:
-            op = to_numeric_matrix_type(op)
+            op = unp.to_numeric_matrix_type(op)
 
         if self.frame_basis is None or op is None:
             return op
 
         if isinstance(op, list):
             return [self.operator_into_frame_basis(x, convert_type=False) for x in op]
-        elif type(op).__name__ == "BCOO":
-            return self.frame_basis_adjoint @ jsparse_matmul(op, self.frame_basis.data)
-        else:
-            # parentheses are necessary for sparse op evaluation
-            return self.frame_basis_adjoint @ (op @ self.frame_basis)
+
+        return unp.rmatmul(unp.matmul(op, self.frame_basis), self.frame_basis_adjoint)
 
     def operator_out_of_frame_basis(
         self,
-        op: Union[Operator, List[Operator], Array, csr_matrix, None],
+        op: Union[Operator, List[Operator], ArrayLike, None],
         convert_type: bool = True,
-    ) -> Array:
+    ) -> ArrayLike:
         r"""Take an operator out of the frame basis, i.e. return
         ``self.frame_basis @ to_array(op) @ self.frame_basis_adjoint``.
 
@@ -209,26 +209,22 @@ def operator_out_of_frame_basis(
                           that ``op`` is a handled input type.
 
         Returns:
-            Array: The operator in the frame basis.
+            ArrayLike: The operator in the frame basis.
         """
         if convert_type:
-            op = to_numeric_matrix_type(op)
-
+            op = unp.to_numeric_matrix_type(op)
         if self.frame_basis is None or op is None:
             return op
 
         if isinstance(op, list):
             return [self.operator_out_of_frame_basis(x, convert_type=False) for x in op]
-        elif type(op).__name__ == "BCOO":
-            return self.frame_basis @ jsparse_matmul(op, self.frame_basis_adjoint.data)
-        else:
-            # parentheses are necessary for sparse op evaluation
-            return self.frame_basis @ (op @ self.frame_basis_adjoint)
+
+        return unp.rmatmul(unp.matmul(op, self.frame_basis_adjoint), self.frame_basis)
 
     def state_into_frame(
         self,
         t: float,
-        y: Array,
+        y: ArrayLike,
         y_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
     ):
@@ -242,9 +238,9 @@ def state_into_frame(
             return_in_frame_basis: Whether or not to return the result in the frame basis.
 
         Returns:
-            Array: The state in the rotating frame.
+            ArrayLike: The state in the rotating frame.
         """
-        y = to_numeric_matrix_type(y)
+        y = unp.to_numeric_matrix_type(y)
         if self._frame_operator is None:
             return y
 
@@ -255,7 +251,7 @@ def state_into_frame(
             out = self.state_into_frame_basis(out)
 
         # go into the frame using fast diagonal matrix multiplication
-        out = (np.exp(self.frame_diag * (-t)) * out.transpose()).transpose()  # = e^{tF}out
+        out = (unp.exp(self.frame_diag * (-t)) * out.transpose()).transpose()  # = e^{tF}out
 
         # if output is requested to not be in the frame basis, convert it
         if not return_in_frame_basis:
@@ -266,10 +262,10 @@ def state_into_frame(
     def state_out_of_frame(
         self,
         t: float,
-        y: Array,
+        y: ArrayLike,
         y_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
-    ) -> Array:
+    ) -> ArrayLike:
         r"""Take a state out of the rotating frame, i.e. return ``exp(tF) @ y``.
 
         Calls ``self.state_into_frame`` with time reversed.
@@ -282,19 +278,19 @@ def state_out_of_frame(
             return_in_frame_basis: Whether or not to return the result in the frame basis.
 
         Returns:
-            Array: The state out of the rotating frame.
+            ArrayLike: The state out of the rotating frame.
         """
         return self.state_into_frame(-t, y, y_in_frame_basis, return_in_frame_basis)
 
     def _conjugate_and_add(
         self,
         t: float,
-        operator: Union[Array, csr_matrix],
-        op_to_add_in_fb: Optional[Array] = None,
+        operator: ArrayLike,
+        op_to_add_in_fb: Optional[ArrayLike] = None,
         operator_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
         vectorized_operators: Optional[bool] = False,
-    ) -> Union[Array, csr_matrix]:
+    ) -> ArrayLike:
         r"""General helper function for computing :math:`\exp(-tF)G\exp(tF) + B`.
 
         Note: :math:`B` is added in the frame basis before any potential final change
@@ -314,13 +310,13 @@ def _conjugate_and_add(
             operator_in_fame_basis: Whether ``operator`` is already in the basis in which the frame
                 operator is diagonal.
             vectorized_operators: Whether ``operator`` is passed as a vectorized, ``(dim**2,)``
-                Array, rather than a ``(dim,dim)`` Array.
+                array, rather than a ``(dim,dim)`` array.
 
         Returns:
-            Array of the newly conjugated operator.
+            ArrayLike: Array of the newly conjugated operator.
         """
-        operator = to_numeric_matrix_type(operator)
-        op_to_add_in_fb = to_numeric_matrix_type(op_to_add_in_fb)
+        operator = unp.to_numeric_matrix_type(operator)
+        op_to_add_in_fb = unp.to_numeric_matrix_type(op_to_add_in_fb)
 
         if vectorized_operators:
             # If passing vectorized operator, undo vectorization temporarily
@@ -338,10 +334,8 @@ def _conjugate_and_add(
                 return operator
             else:
                 if op_to_add_in_fb is not None:
-                    if issparse(operator):
-                        op_to_add_in_fb = csr_matrix(op_to_add_in_fb)
-                    elif type(operator).__name__ == "BCOO":
-                        op_to_add_in_fb = to_BCOO(op_to_add_in_fb)
+                    if issparse(out) or type(out).__name__ == "BCOO":
+                        op_to_add_in_fb = DYNAMICS_ALIAS(like=out).asarray(op_to_add_in_fb)
 
                 return operator + op_to_add_in_fb
 
@@ -354,20 +348,18 @@ def _conjugate_and_add(
         # get frame transformation matrix in diagonal basis
         # assumption that F is anti-Hermitian implies conjugation of
         # diagonal gives inversion
-        exp_freq = np.exp(self.frame_diag * t)
+        exp_freq = unp.exp(self.frame_diag * t)
         frame_mat = exp_freq.conj().reshape(self.dim, 1) * exp_freq
+        # need to define unp.matmul
         if issparse(out):
             out = out.multiply(frame_mat)
         elif type(out).__name__ == "BCOO":
             out = out * frame_mat.data
         else:
             out = frame_mat * out
-
         if op_to_add_in_fb is not None:
-            if issparse(out):
-                op_to_add_in_fb = csr_matrix(op_to_add_in_fb)
-            elif type(out).__name__ == "BCOO":
-                op_to_add_in_fb = to_BCOO(op_to_add_in_fb)
+            if issparse(out) or type(out).__name__ == "BCOO":
+                op_to_add_in_fb = DYNAMICS_ALIAS(like=out).asarray(op_to_add_in_fb)
 
             out = out + op_to_add_in_fb
 
@@ -386,11 +378,11 @@ def _conjugate_and_add(
     def operator_into_frame(
         self,
         t: float,
-        operator: Union[Operator, Array, csr_matrix],
+        operator: Union[Operator, ArrayLike],
         operator_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
         vectorized_operators: Optional[bool] = False,
-    ) -> Array:
+    ) -> ArrayLike:
         r"""Bring an operator into the frame, i.e. return
         ``exp(-tF) @ operator @ exp(tF)``.
 
@@ -403,10 +395,10 @@ def operator_into_frame(
                 the frame is diagonal.
             return_in_frame_basis: Whether or not to return the result in the frame basis.
             vectorized_operators: Whether ``operator`` is passed as a vectorized,
-                ``(dim**2,)`` Array, rather than a ``(dim,dim)`` Array.
+                ``(dim**2,)`` array, rather than a ``(dim,dim)`` array.
 
         Returns:
-            Array: The operator in the rotating frame.
+            ArrayLike: The operator in the rotating frame.
         """
         return self._conjugate_and_add(
             t,
@@ -419,11 +411,11 @@ def operator_into_frame(
     def operator_out_of_frame(
         self,
         t: float,
-        operator: Union[Operator, Array, csr_matrix],
+        operator: Union[Operator, ArrayLike],
         operator_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
         vectorized_operators: Optional[bool] = False,
-    ):
+    ) -> ArrayLike:
         r"""Bring an operator into the rotating frame, i.e. return
         ``exp(tF) @ operator @ exp(-tF)``.
 
@@ -436,10 +428,10 @@ def operator_out_of_frame(
                 the frame is diagonal.
             return_in_frame_basis: Whether or not to return the result in the frame basis.
             vectorized_operators: Whether ``operator`` is passed as a vectorized, ``(dim**2,)``
-                Array, rather than a ``(dim,dim)`` Array.
+                array, rather than a ``(dim,dim)`` array.
 
         Returns:
-            Array: The operator out of the rotating frame.
+            ArrayLike: The operator out of the rotating frame.
         """
         return self.operator_into_frame(
             -t,
@@ -452,11 +444,11 @@ def operator_out_of_frame(
     def generator_into_frame(
         self,
         t: float,
-        operator: Union[Operator, Array, csr_matrix],
+        operator: Union[Operator, ArrayLike],
         operator_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
         vectorized_operators: Optional[bool] = False,
-    ):
+    ) -> ArrayLike:
         r"""Take an generator into the rotating frame, i.e. return
         ``exp(-tF) @ operator @ exp(tF) - F``.
 
@@ -469,19 +461,19 @@ def generator_into_frame(
                 the frame is diagonal.
             return_in_frame_basis: Whether or not to return the result in the frame basis.
             vectorized_operators: Whether ``operator`` is passed as a vectorized, ``(dim**2,)``
-                Array, rather than a ``(dim,dim)`` Array.
+                array, rather than a ``(dim,dim)`` array.
 
         Returns:
-            Array: The generator in the rotating frame.
+            ArrayLike: The generator in the rotating frame.
         """
         if self.frame_operator is None:
-            return to_numeric_matrix_type(operator)
+            return unp.to_numeric_matrix_type(operator)
         else:
             # conjugate and subtract the frame diagonal
             return self._conjugate_and_add(
                 t,
                 operator,
-                op_to_add_in_fb=-np.diag(self.frame_diag),
+                op_to_add_in_fb=-unp.diag(self.frame_diag),
                 operator_in_frame_basis=operator_in_frame_basis,
                 return_in_frame_basis=return_in_frame_basis,
                 vectorized_operators=vectorized_operators,
@@ -490,10 +482,10 @@ def generator_into_frame(
     def generator_out_of_frame(
         self,
         t: float,
-        operator: Union[Operator, Array, csr_matrix],
+        operator: Union[Operator, ArrayLike],
         operator_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
-    ) -> Array:
+    ) -> ArrayLike:
         r"""Take an operator out of the frame using the generator transformaton rule, i.e. return
         ``exp(tF) @ operator @ exp(-tF) + F``.
 
@@ -507,16 +499,16 @@ def generator_out_of_frame(
             return_in_frame_basis: Whether or not to return the result in the frame basis.
 
         Returns:
-            Array: The generator out of the rotating frame.
+            ArrayLike: The generator out of the rotating frame.
         """
         if self.frame_operator is None:
-            return to_numeric_matrix_type(operator)
+            return unp.to_numeric_matrix_type(operator)
         else:
             # conjugate and add the frame diagonal
             return self._conjugate_and_add(
                 -t,
                 operator,
-                op_to_add_in_fb=Array(np.diag(self.frame_diag)),
+                op_to_add_in_fb=unp.diag(self.frame_diag),
                 operator_in_frame_basis=operator_in_frame_basis,
                 return_in_frame_basis=return_in_frame_basis,
             )
@@ -529,7 +521,7 @@ def vectorized_frame_basis(self):
             return None
 
         if self._vectorized_frame_basis is None:
-            self._vectorized_frame_basis = np.kron(self.frame_basis.conj(), self.frame_basis)
+            self._vectorized_frame_basis = unp.kron(self.frame_basis.conj(), self.frame_basis)
             self._vectorized_frame_basis_adjoint = self._vectorized_frame_basis.conj().transpose()
 
         return self._vectorized_frame_basis
@@ -551,10 +543,10 @@ def vectorized_frame_basis_adjoint(self):
     def vectorized_map_into_frame(
         self,
         time: float,
-        op: Array,
+        op: ArrayLike,
         operator_in_frame_basis: Optional[bool] = False,
         return_in_frame_basis: Optional[bool] = False,
-    ) -> Array:
+    ) -> ArrayLike:
         r"""Given an operator ``op`` of dimension ``dim**2`` assumed to represent vectorized linear
         map in column stacking convention, returns:
 
@@ -567,12 +559,12 @@ def vectorized_map_into_frame(
 
         Args:
             time: The time :math:`t`.
-            op: The ``(dim**2,dim**2)`` Array.
+            op: The ``(dim**2,dim**2)`` array.
             operator_in_frame_basis: Whether the operator is in the frame basis.
             return_in_frame_basis: Whether the operator should be returned in the frame basis.
 
         Returns:
-            ``op`` in the frame.
+            ArrayLike: ``op`` in the frame.
         """
         if self.frame_diag is not None:
             # Put the vectorized operator into the frame basis
@@ -596,7 +588,9 @@ def vectorized_map_into_frame(
         return op
 
 
-def _enforce_anti_herm(mat: Array, atol: Optional[float] = 1e-10, rtol: Optional[float] = 1e-10):
+def _enforce_anti_herm(
+    mat: ArrayLike, atol: Optional[float] = 1e-10, rtol: Optional[float] = 1e-10
+) -> ArrayLike:
     r"""Given ``mat``, the logic of this function is:
         - if ``mat`` is hermitian, return ``-1j * mat``
         - if ``mat`` is anti-hermitian, return ``mat``
@@ -613,20 +607,17 @@ def _enforce_anti_herm(mat: Array, atol: Optional[float] = 1e-10, rtol: Optional
         rtol: The relative tolerance.
 
     Returns:
-        Array: Anti-hermitian version of ``mat``, if applicable.
+        ArrayLike: Anti-hermitian version of ``mat``, if applicable.
 
     Raises:
         ImportError: If the backend is jax and jax is not installed.
         QiskitError: If ``mat`` is not Hermitian or anti-Hermitian.
     """
-    mat = to_array(mat)
-    mat = Array(mat, dtype=complex)
+    mat = unp.to_dense(mat)
 
-    if mat.backend == "jax":
+    if not isinstance(mat, np.ndarray):
         from jax.lax import cond
 
-        mat = mat.data
-
         if mat.ndim == 1:
             # this function checks if pure imaginary. If yes it returns the
             # array, otherwise it multiplies it by jnp.nan to raise an error
@@ -637,7 +628,7 @@ def anti_herm_conditional(b):
 
             # Check if it is purely real, if not apply anti_herm_conditional
             herm_pred = jnp.allclose(mat, mat.conj(), atol=atol, rtol=rtol)
-            return Array(cond(herm_pred, lambda A: -1j * A, anti_herm_conditional, mat))
+            return unp.asarray(cond(herm_pred, lambda A: -1j * A, anti_herm_conditional, mat))
         else:
             # this function checks if anti-hermitian, if yes returns the array,
             # otherwise it multiplies it by jnp.nan
@@ -648,7 +639,7 @@ def anti_herm_conditional(b):
             # the following lines check if a is hermitian, otherwise it feeds
             # it into the anti_herm_conditional
             herm_pred = jnp.allclose(mat, mat.conj().transpose(), atol=atol, rtol=rtol)
-            return Array(cond(herm_pred, lambda A: -1j * A, anti_herm_conditional, mat))
+            return unp.asarray(cond(herm_pred, lambda A: -1j * A, anti_herm_conditional, mat))
 
     else:
         if mat.ndim == 1:
diff --git a/qiskit_dynamics/signals/signals.py b/qiskit_dynamics/signals/signals.py
index ab72fe9e3..e6f308211 100644
--- a/qiskit_dynamics/signals/signals.py
+++ b/qiskit_dynamics/signals/signals.py
@@ -27,9 +27,9 @@
 from matplotlib import pyplot as plt
 
 from qiskit import QiskitError
-from qiskit_dynamics.arraylias_state import ArrayLike
-from qiskit_dynamics.arraylias_state import DYNAMICS_ALIAS as alias
-from qiskit_dynamics.arraylias_state import DYNAMICS_NUMPY as unp
+from qiskit_dynamics.arraylias.arraylias_state import ArrayLike
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_ALIAS as alias
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_NUMPY as unp
 
 
 class Signal:
diff --git a/qiskit_dynamics/type_utils.py b/qiskit_dynamics/type_utils.py
index bd1d9c748..c7912416e 100644
--- a/qiskit_dynamics/type_utils.py
+++ b/qiskit_dynamics/type_utils.py
@@ -432,13 +432,10 @@ def to_csr(
 @requires_backend("jax")
 def to_BCOO(op: Union[Operator, Array, List[Operator], List[Array], spmatrix, "BCOO"]) -> "BCOO":
     """Convert input op or list of ops to a jax BCOO sparse array.
-
     Calls ``to_array`` to handle general conversion to a numpy or jax array, then
     builds the BCOO sparse array from the result.
-
     Args:
         op: Operator or list of operators to convert.
-
     Returns:
         BCOO: BCOO sparse version of the operator.
     """
diff --git a/setup.py b/setup.py
index 35c9f48de..0404a2bc7 100644
--- a/setup.py
+++ b/setup.py
@@ -23,8 +23,8 @@
     "multiset>=3.0.1",
 ]
 
-jax_extras = ['jax>=0.2.26',
-              'jaxlib>=0.1.75']
+jax_extras = ['jax>=0.2.26, <= 0.4.6',
+              'jaxlib>=0.1.75, <= 0.4.6']
 
 PACKAGES = setuptools.find_packages(exclude=['test*'])
 
diff --git a/test/dynamics/common.py b/test/dynamics/common.py
index 6483c0efb..496d848ad 100644
--- a/test/dynamics/common.py
+++ b/test/dynamics/common.py
@@ -28,6 +28,7 @@
 except ImportError:
     pass
 
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_NUMPY as unp
 from qiskit_dynamics.array import Array, wrap
 
 
@@ -36,8 +37,8 @@ class QiskitDynamicsTestCase(unittest.TestCase):
 
     def assertAllClose(self, A, B, rtol=1e-8, atol=1e-8):
         """Call np.allclose and assert true."""
-        A = Array(A)
-        B = Array(B)
+        A = unp.asarray(A)
+        B = unp.asarray(B)
 
         self.assertTrue(np.allclose(A, B, rtol=rtol, atol=atol))
 
diff --git a/test/dynamics/models/test_rotating_frame.py b/test/dynamics/models/test_rotating_frame.py
index e17a83ad7..0cf12dd7d 100644
--- a/test/dynamics/models/test_rotating_frame.py
+++ b/test/dynamics/models/test_rotating_frame.py
@@ -19,27 +19,27 @@
 from qiskit.quantum_info.operators import Operator
 from scipy.sparse import csr_matrix
 from qiskit_dynamics.models.rotating_frame import RotatingFrame
-from qiskit_dynamics.array import Array
-from qiskit_dynamics.type_utils import to_BCOO, to_array
-from ..common import QiskitDynamicsTestCase, TestJaxBase
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_NUMPY as unp
+from qiskit_dynamics.arraylias.arraylias_state import ArrayLike
+from ..common import QiskitDynamicsTestCase, TestJaxBase, TestNumpyBase
 
 
-class TestRotatingFrame(QiskitDynamicsTestCase):
+class TestRotatingFrame(QiskitDynamicsTestCase, TestNumpyBase):
     """Tests for RotatingFrame."""
 
     def setUp(self):
-        self.X = Array(Operator.from_label("X").data)
-        self.Y = Array(Operator.from_label("Y").data)
-        self.Z = Array(Operator.from_label("Z").data)
+        self.X = self.asarray(Operator.from_label("X").data)
+        self.Y = self.asarray(Operator.from_label("Y").data)
+        self.Z = self.asarray(Operator.from_label("Z").data)
 
     def test_instantiation_errors(self):
         """Check different modes of error raising for frame setting."""
 
         with self.assertRaisesRegex(QiskitError, "Hermitian or anti-Hermitian"):
-            RotatingFrame(Array([1.0, 1j]))
+            RotatingFrame(self.asarray([1.0, 1j]))
 
         with self.assertRaisesRegex(QiskitError, "Hermitian or anti-Hermitian"):
-            RotatingFrame(Array([[1.0, 0.0], [0.0, 1j]]))
+            RotatingFrame(self.asarray([[1.0, 0.0], [0.0, 1j]]))
 
         with self.assertRaisesRegex(QiskitError, "Hermitian or anti-Hermitian"):
             RotatingFrame(self.Z + 1j * self.X)
@@ -52,12 +52,12 @@ def test_state_out_of_frame_basis(self):
             low=-10, high=10, size=(6, 6)
         )
 
-        frame_op = Array(rand_op - rand_op.conj().transpose())
+        frame_op = unp.asarray(rand_op - rand_op.conj().transpose())
         rotating_frame = RotatingFrame(frame_op)
 
         _, U = np.linalg.eigh(1j * frame_op)
 
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -78,13 +78,13 @@ def test_operator_into_frame_basis(self):
             low=-10, high=10, size=(10, 10)
         )
 
-        frame_op = Array(rand_op - rand_op.conj().transpose())
+        frame_op = unp.asarray(rand_op - rand_op.conj().transpose())
         rotating_frame = RotatingFrame(frame_op)
 
         _, U = np.linalg.eigh(1j * frame_op)
         Uadj = U.conj().transpose()
 
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(10, 10))
             + 1j * rng.uniform(low=-10, high=10, size=(10, 10))
         )
@@ -101,7 +101,7 @@ def test_operator_into_frame_basis_sparse_list(self):
         """Test state_into_frame_basis for a list of sparse arrays."""
 
         ops = [csr_matrix([[0.0, 1.0], [1.0, 0.0]]), csr_matrix([[1.0, 0.0], [0.0, -1.0]])]
-        rotating_frame = RotatingFrame(np.array([[0.0, 1.0], [1.0, 0.0]]))
+        rotating_frame = RotatingFrame(unp.asarray([[0.0, 1.0], [1.0, 0.0]]))
 
         val = rotating_frame.operator_into_frame_basis(ops)
         U = rotating_frame.frame_basis
@@ -113,8 +113,7 @@ def test_operator_out_of_frame_basis_sparse_list(self):
         """Test state_out_of_frame_basis for a list of sparse arrays."""
 
         ops = [csr_matrix([[0.0, 1.0], [1.0, 0.0]]), csr_matrix([[1.0, 0.0], [0.0, -1.0]])]
-        rotating_frame = RotatingFrame(np.array([[0.0, 1.0], [1.0, 0.0]]))
-
+        rotating_frame = RotatingFrame(unp.asarray([[0.0, 1.0], [1.0, 0.0]]))
         val = rotating_frame.operator_out_of_frame_basis(ops)
         U = rotating_frame.frame_basis
         Uadj = rotating_frame.frame_basis_adjoint
@@ -124,17 +123,17 @@ def test_operator_out_of_frame_basis_sparse_list(self):
     def test_state_transformations_no_frame(self):
         """Test frame transformations with no frame."""
 
-        rotating_frame = RotatingFrame(Array(np.zeros(2)))
+        rotating_frame = RotatingFrame(unp.asarray(np.zeros(2)))
 
         t = 0.123
-        y = Array([1.0, 1j])
+        y = unp.asarray([1.0, 1j])
         out = rotating_frame.state_into_frame(t, y)
         self.assertAllClose(out, y)
         out = rotating_frame.state_out_of_frame(t, y)
         self.assertAllClose(out, y)
 
         t = 100.12498
-        y = Array(np.eye(2))
+        y = unp.asarray(np.eye(2))
         out = rotating_frame.state_into_frame(t, y)
         self.assertAllClose(out, y)
         out = rotating_frame.state_out_of_frame(t, y)
@@ -142,9 +141,9 @@ def test_state_transformations_no_frame(self):
 
     def test_state_into_frame_2_level(self):
         """Test state_into_frame with a non-trival frame."""
-        frame_op = -1j * np.pi * (self.X + 0.1 * self.Y + 12.0 * self.Z).data
+        frame_op = -1j * np.pi * (self.X + 0.1 * self.Y + 12.0 * self.Z)
         t = 1312.132
-        y0 = Array([[1.0, 2.0], [3.0, 4.0]])
+        y0 = unp.asarray([[1.0, 2.0], [3.0, 4.0]])
 
         self._test_state_into_frame(t, frame_op, y0)
         self._test_state_into_frame(t, frame_op, y0, y_in_frame_basis=True)
@@ -160,10 +159,10 @@ def test_state_into_frame_pseudo_random(self):
             low=-10, high=10, size=(5, 5)
         )
 
-        frame_op = Array(rand_op - rand_op.conj().transpose())
+        frame_op = unp.asarray(rand_op - rand_op.conj().transpose())
 
         t = 1312.132
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(5, 5))
             + 1j * rng.uniform(low=-10, high=10, size=(5, 5))
         )
@@ -189,7 +188,7 @@ def _test_state_into_frame(
         if not y_in_frame_basis:
             expected = U.conj().transpose() @ expected
 
-        expected = np.diag(np.exp(-t * Array(evals))) @ expected
+        expected = np.diag(np.exp(-t * unp.asarray(evals))) @ expected
 
         if not return_in_frame_basis:
             expected = U @ expected
@@ -198,9 +197,9 @@ def _test_state_into_frame(
 
     def test_state_out_of_frame_2_level(self):
         """Test state_out_of_frame with a non-trival frame."""
-        frame_op = -1j * np.pi * (3.1 * self.X + 1.1 * self.Y + 12.0 * self.Z).data
+        frame_op = -1j * np.pi * (3.1 * self.X + 1.1 * self.Y + 12.0 * self.Z)
         t = 122.132
-        y0 = Array([[1.0, 2.0], [3.0, 4.0]])
+        y0 = unp.asarray([[1.0, 2.0], [3.0, 4.0]])
         self._test_state_out_of_frame(t, frame_op, y0)
         self._test_state_out_of_frame(t, frame_op, y0, y_in_frame_basis=True)
         self._test_state_out_of_frame(t, frame_op, y0, return_in_frame_basis=True)
@@ -215,10 +214,10 @@ def test_state_out_of_frame_pseudo_random(self):
             low=-10, high=10, size=(6, 6)
         )
 
-        frame_op = Array(rand_op - rand_op.conj().transpose())
+        frame_op = unp.asarray(rand_op - rand_op.conj().transpose())
 
         t = rng.uniform(low=-100, high=100)
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -235,7 +234,7 @@ def _test_state_out_of_frame(
         self, t, frame_op, y, y_in_frame_basis=False, return_in_frame_basis=False
     ):
         evals, U = np.linalg.eigh(1j * frame_op)
-        evals = -1j * Array(evals)
+        evals = -1j * unp.asarray(evals)
 
         rotating_frame = RotatingFrame(frame_op)
 
@@ -258,10 +257,10 @@ def test_operator_into_frame(self):
             low=-10, high=10, size=(6, 6)
         )
 
-        frame_op = Array(rand_op - rand_op.conj().transpose())
+        frame_op = unp.asarray(rand_op - rand_op.conj().transpose())
 
         t = rng.uniform(low=-100, high=100)
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -278,7 +277,7 @@ def _test_operator_into_frame(
         self, t, frame_op, y, y_in_frame_basis=False, return_in_frame_basis=False
     ):
         evals, U = np.linalg.eigh(1j * frame_op)
-        evals = -1j * Array(evals)
+        evals = -1j * unp.asarray(evals)
         Uadj = U.conj().transpose()
 
         rotating_frame = RotatingFrame(frame_op)
@@ -298,7 +297,7 @@ def _test_operator_into_frame(
     def test_operator_out_of_frame(self):
         """Test operator_out_of_frame."""
         rng = np.random.default_rng(37164093)
-        rand_op = Array(
+        rand_op = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -306,7 +305,7 @@ def test_operator_out_of_frame(self):
         frame_op = rand_op - rand_op.conj().transpose()
 
         t = rng.uniform(low=-100, high=100)
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -323,7 +322,7 @@ def _test_operator_out_of_frame(
         self, t, frame_op, y, y_in_frame_basis=False, return_in_frame_basis=False
     ):
         evals, U = np.linalg.eigh(1j * frame_op)
-        evals = -1j * Array(evals)
+        evals = -1j * unp.asarray(evals)
         Uadj = U.conj().transpose()
 
         rotating_frame = RotatingFrame(frame_op)
@@ -343,7 +342,7 @@ def _test_operator_out_of_frame(
     def test_generator_into_frame(self):
         """Test operator_out_of_frame."""
         rng = np.random.default_rng(111)
-        rand_op = Array(
+        rand_op = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -351,7 +350,7 @@ def test_generator_into_frame(self):
         frame_op = rand_op - rand_op.conj().transpose()
 
         t = rng.uniform(low=-100, high=100)
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -369,7 +368,7 @@ def _test_generator_into_frame(
     ):
         """Helper function for testing generator_into_frame."""
         evals, U = np.linalg.eigh(1j * frame_op)
-        evals = -1j * Array(evals)
+        evals = -1j * unp.asarray(evals)
         Uadj = U.conj().transpose()
 
         rotating_frame = RotatingFrame(frame_op)
@@ -394,10 +393,10 @@ def test_generator_out_of_frame(self):
             low=-10, high=10, size=(6, 6)
         )
 
-        frame_op = Array(rand_op - rand_op.conj().transpose())
+        frame_op = unp.asarray(rand_op - rand_op.conj().transpose())
 
         t = rng.uniform(low=-100, high=100)
-        y0 = Array(
+        y0 = unp.asarray(
             rng.uniform(low=-10, high=10, size=(6, 6))
             + 1j * rng.uniform(low=-10, high=10, size=(6, 6))
         )
@@ -415,7 +414,7 @@ def _test_generator_out_of_frame(
     ):
         """Helper function for testing generator_into_frame."""
         evals, U = np.linalg.eigh(1j * frame_op)
-        evals = -1j * Array(evals)
+        evals = -1j * unp.asarray(evals)
         Uadj = U.conj().transpose()
 
         rotating_frame = RotatingFrame(frame_op)
@@ -437,12 +436,12 @@ def test_vectorized_conjugate_and_add_conventions(self):
         """Test whether passing a vectorized (dim**2, k) operator to _conjugate_and_add
         with vectorized_operators = True is the same as passing a (k,dim,dim) array of
         operators."""
-        vectorized_rhos = np.array(
+        vectorized_rhos = self.asarray(
             [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]
         ).transpose()
 
         nonvectord_rhos = vectorized_rhos.reshape((2, 2, 3), order="F").transpose([2, 0, 1])
-        rotating_frame = RotatingFrame(np.array([1j, 2j]))
+        rotating_frame = RotatingFrame(self.asarray([1j, 2j]))
 
         vectorized_result = rotating_frame._conjugate_and_add(
             0.1, vectorized_rhos, vectorized_operators=True
@@ -562,22 +561,22 @@ def test_state_transformations_no_frame_array_type(self):
         rotating_frame = RotatingFrame(None)
 
         t = 0.123
-        y = Array([1.0, 1j])
+        y = unp.asarray([1.0, 1j])
         out = rotating_frame.state_into_frame(t, y)
         self.assertAllClose(out, y)
-        self.assertTrue(isinstance(out, Array))
+        self.assertTrue(isinstance(out, ArrayLike))
         out = rotating_frame.state_out_of_frame(t, y)
         self.assertAllClose(out, y)
-        self.assertTrue(isinstance(out, Array))
+        self.assertTrue(isinstance(out, ArrayLike))
 
         t = 100.12498
-        y = Array(np.eye(2))
+        y = unp.asarray(np.eye(2))
         out = rotating_frame.state_into_frame(t, y)
         self.assertAllClose(out, y)
-        self.assertTrue(isinstance(out, Array))
+        self.assertTrue(isinstance(out, ArrayLike))
         out = rotating_frame.state_out_of_frame(t, y)
         self.assertAllClose(out, y)
-        self.assertTrue(isinstance(out, Array))
+        self.assertTrue(isinstance(out, ArrayLike))
 
 
 class TestRotatingJAXBCOO(QiskitDynamicsTestCase, TestJaxBase):
@@ -586,16 +585,16 @@ class TestRotatingJAXBCOO(QiskitDynamicsTestCase, TestJaxBase):
     def test_conjugate_and_add_BCOO(self):
         """Test _conjugate_and_add with operator being BCOO."""
 
-        rotating_frame = RotatingFrame(np.array([1.0, -1.0]))
+        rotating_frame = RotatingFrame(unp.asarray([1.0, -1.0]))
 
         t = 0.123
-        op = to_BCOO(np.array([[1.0, -1j], [0.0, 1.0]]))
-        op_to_add = to_BCOO(np.array([[0.0, -0.11j], [0.0, 1.0]]))
+        op = unp.to_sparse(self.asarray([[1.0, -1j], [0.0, 1.0]]))
+        op_to_add = unp.to_sparse(self.asarray([[0.0, -0.11j], [0.0, 1.0]]))
         out = rotating_frame._conjugate_and_add(t, op, op_to_add)
         self.assertTrue(type(out).__name__ == "BCOO")
-
         self.assertAllClose(
-            to_array(out), rotating_frame._conjugate_and_add(t, to_array(op), to_array(op_to_add))
+            unp.to_dense(out),
+            rotating_frame._conjugate_and_add(t, unp.to_dense(op), unp.to_dense(op_to_add)),
         )
 
     def test_operator_into_frame_basis(self):
@@ -603,11 +602,11 @@ def test_operator_into_frame_basis(self):
         frame specified as full matrix.
         """
 
-        rotating_frame = RotatingFrame(np.array([[1.0, 0.0], [0.0, -1.0]]))
+        rotating_frame = RotatingFrame(unp.asarray([[1.0, 0.0], [0.0, -1.0]]))
 
-        op = to_BCOO(np.array([[1.0, -1j], [0.0, 1.0]]))
+        op = unp.to_sparse(self.asarray([[1.0, -1j], [0.0, 1.0]]))
         output = rotating_frame.operator_into_frame_basis(op)
-        expected = rotating_frame.operator_into_frame_basis(to_array(op))
+        expected = rotating_frame.operator_into_frame_basis(unp.to_dense(op))
 
         self.assertAllClose(output, expected)
 
@@ -616,16 +615,16 @@ def test_operator_out_of_frame_basis(self):
         frame specified as full matrix.
         """
 
-        rotating_frame = RotatingFrame(np.array([[1.0, 0.0], [0.0, -1.0]]))
+        rotating_frame = RotatingFrame(unp.asarray([[1.0, 0.0], [0.0, -1.0]]))
 
-        op = to_BCOO(np.array([[1.0, -1j], [0.0, 1.0]]))
+        op = unp.to_sparse(self.asarray([[1.0, -1j], [0.0, 1.0]]))
         output = rotating_frame.operator_out_of_frame_basis(op)
-        expected = rotating_frame.operator_out_of_frame_basis(to_array(op))
+        expected = rotating_frame.operator_out_of_frame_basis(unp.to_dense(op))
 
         self.assertAllClose(output, expected)
 
 
-class TestRotatingFrameJax(TestRotatingFrame, TestJaxBase):
+class TestRotatingFrameJax(TestJaxBase, TestRotatingFrame):
     """Jax version of TestRotatingFrame tests.
 
     Note: This class has more tests due to inheritance.
@@ -639,10 +638,10 @@ def test_instantiation_errors(self):
         # pylint: disable=import-outside-toplevel
         import jax.numpy as jnp
 
-        rotating_frame = RotatingFrame(Array([1.0, 1j]))
+        rotating_frame = RotatingFrame(jnp.asarray([1.0, 1j]))
         self.assertTrue(jnp.isnan(rotating_frame.frame_diag[0]))
 
-        rotating_frame = RotatingFrame(Array([[1.0, 0.0], [0.0, 1j]]))
+        rotating_frame = RotatingFrame(jnp.asarray([[1.0, 0.0], [0.0, 1j]]))
         self.assertTrue(jnp.isnan(rotating_frame.frame_diag[0]))
 
         rotating_frame = RotatingFrame(self.Z + 1j * self.X)
@@ -651,19 +650,19 @@ def test_instantiation_errors(self):
     def test_jitting(self):
         """Test jitting of state_into_frame and _conjugate_and_add."""
 
-        rotating_frame = RotatingFrame(Array([1.0, -1.0]))
+        rotating_frame = RotatingFrame(unp.asarray([1.0, -1.0]))
 
-        self.jit_wrap(rotating_frame.state_into_frame)(t=0.1, y=np.array([0.0, 1.0]))
+        self.jit_wrap(rotating_frame.state_into_frame)(t=0.1, y=unp.asarray([0.0, 1.0]))
         self.jit_wrap(rotating_frame._conjugate_and_add)(
-            t=0.1, operator=np.array([[0.0, 1.0], [1.0, 0.0]])
+            t=0.1, operator=unp.asarray([[0.0, 1.0], [1.0, 0.0]])
         )
 
     def test_jit_and_grad(self):
         """Test jitting and gradding of state_into_frame and _conjugate_and_add."""
 
-        rotating_frame = RotatingFrame(Array([1.0, -1.0]))
+        rotating_frame = RotatingFrame(unp.asarray([1.0, -1.0]))
 
-        self.jit_grad_wrap(rotating_frame.state_into_frame)(0.1, np.array([0.0, 1.0]))
+        self.jit_grad_wrap(rotating_frame.state_into_frame)(0.1, unp.asarray([0.0, 1.0]))
         self.jit_grad_wrap(rotating_frame._conjugate_and_add)(
-            0.1, np.array([[0.0, 1.0], [1.0, 0.0]])
+            0.1, unp.asarray([[0.0, 1.0], [1.0, 0.0]])
         )
diff --git a/test/dynamics/signals/test_signals.py b/test/dynamics/signals/test_signals.py
index 07653f9f3..e823eecbc 100644
--- a/test/dynamics/signals/test_signals.py
+++ b/test/dynamics/signals/test_signals.py
@@ -19,7 +19,7 @@
 
 from qiskit_dynamics.signals import Signal, DiscreteSignal, DiscreteSignalSum, SignalList
 from qiskit_dynamics.signals.signals import to_SignalSum
-from qiskit_dynamics.arraylias_state import DYNAMICS_NUMPY as unp
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_NUMPY as unp
 
 from ..common import QiskitDynamicsTestCase, TestJaxBase, TestNumpyBase
 
diff --git a/test/dynamics/test_type_utils.py b/test/dynamics/test_type_utils.py
index 177fa77d3..c5d68a875 100644
--- a/test/dynamics/test_type_utils.py
+++ b/test/dynamics/test_type_utils.py
@@ -26,15 +26,17 @@
     vec_dissipator,
     vec_commutator,
     to_array,
-    to_csr,
     to_BCOO,
     to_numeric_matrix_type,
 )
+from qiskit_dynamics.arraylias.arraylias_state import DYNAMICS_NUMPY as unp
 
-from .common import QiskitDynamicsTestCase, TestJaxBase, TestQutipBase
+
+from .common import QiskitDynamicsTestCase, TestJaxBase, TestNumpyBase, TestQutipBase
 
 try:
     from jax.experimental import sparse as jsparse
+    import jax.numpy as jnp
 except ImportError:
     pass
 
@@ -274,177 +276,175 @@ def test_sparse_commutator_dissipator(self):
         )
 
 
-class Test_to_array(QiskitDynamicsTestCase):
-    """Tests for to_array."""
+class Test_to_dense(QiskitDynamicsTestCase, TestNumpyBase):
+    """Tests for to_dense."""
 
     def test_None_to_None(self):
         """Test that None input returns None."""
-        self.assertTrue(to_array(None) is None)
+        self.assertTrue(unp.to_dense(None) is None)
 
-    def test_to_array_Operator(self):
-        """Tests for to_array with a single operator"""
+    def test_to_dense_Operator(self):
+        """Tests for to_dense with a single operator"""
         op = Operator.from_label("X")
-        self.assertAllClose(to_array(op), Array([[0, 1], [1, 0]]))
+        self.assertAllClose(unp.to_dense(op), np.array([[0, 1], [1, 0]]))
 
-    def test_to_array_nparray(self):
-        """Tests for to_array with a single numpy array"""
+    def test_to_dense_nparray(self):
+        """Tests for to_dense with a single numpy array"""
         ndarray = np.array([[0, 1], [1, 0]])
-        self.assertAllClose(to_array(ndarray), ndarray)
+        self.assertAllClose(unp.to_dense(ndarray), ndarray)
 
-    def test_to_array_Array(self):
-        """Tests for to_array from a qiskit Array"""
+    def test_to_dense_Array(self):
+        """Tests for to_dense from a qiskit Array"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
-        arr_in = Array(list_of_ops)
-        self.assertAllClose(to_array(arr_in), arr_in)
+        arr_in = unp.asarray(list_of_ops)
+        self.assertAllClose(unp.to_dense(arr_in), arr_in)
 
-    def test_to_array_sparse_matrix(self):
-        """Tests for to_array with a single sparse matrix"""
+    def test_to_dense_sparse_matrix(self):
+        """Tests for to_dense with a single sparse matrix"""
         op = Operator.from_label("X")
         spm = csr_matrix(op)
-        ar = Array(op)
-        self.assertAllClose(to_array(spm), ar)
+        ar = unp.asarray(op)
+        self.assertAllClose(unp.to_dense(spm), ar)
 
-    def test_to_array_Operator_list(self):
-        """Tests for to_array with a list of operators"""
+    def test_to_dense_Operator_list(self):
+        """Tests for to_dense with a list of operators"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
         op_arr = [Operator.from_label(s) for s in "XYZ"]
-        normal_array = Array(np.array(list_of_ops))
-        self.assertAllClose(to_array(op_arr), normal_array)
+        normal_array = unp.asarray(np.array(list_of_ops))
+        self.assertAllClose(unp.to_dense(op_arr), normal_array)
 
-    def test_to_array_nparray_list(self):
-        """Tests for to_array with a list of numpy arrays"""
+    def test_to_dense_nparray_list(self):
+        """Tests for to_dense with a list of numpy arrays"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
         ndarray_list = [np.array(op) for op in list_of_ops]
-        list_of_arrays = [Array(op) for op in list_of_ops]
-        self.assertAllClose(to_array(ndarray_list), list_of_arrays)
+        list_of_arrays = [unp.asarray(op) for op in list_of_ops]
+        self.assertAllClose(unp.to_dense(ndarray_list), list_of_arrays)
 
-    def test_to_array_list_of_arrays(self):
-        """Tests for to_array with a list of numpy arrays"""
+    def test_to_dense_list_of_arrays(self):
+        """Tests for to_dense with a list of numpy arrays"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
         ndarray_list = [np.array(op) for op in list_of_ops]
-        list_of_arrays = [Array(op) for op in list_of_ops]
-        big_array = Array(list_of_arrays)
-        self.assertAllClose(to_array(ndarray_list), big_array)
+        list_of_arrays = [unp.asarray(op) for op in list_of_ops]
+        big_array = unp.asarray(list_of_arrays)
+        self.assertAllClose(unp.to_dense(ndarray_list), big_array)
 
-    def test_to_array_sparse_matrix_list(self):
-        """Tests for to_array with a list of sparse matrices"""
+    def test_to_dense_sparse_matrix_list(self):
+        """Tests for to_dense with a list of sparse matrices"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
-        list_of_arrays = [Array(op) for op in list_of_ops]
+        list_of_arrays = [unp.asarray(op) for op in list_of_ops]
         sparse_matrices = [csr_matrix(op) for op in list_of_ops]
-        self.assertAllClose(to_array(sparse_matrices), list_of_arrays)
+        self.assertAllClose(unp.to_dense(sparse_matrices), list_of_arrays)
 
-    def test_to_array_types(self):
-        """Type conversion tests for to_array"""
+    def test_to_dense_types(self):
+        """Type conversion tests for to_dense"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
         numpy_ops = np.array(list_of_ops)
-        normal_array = Array(np.array(list_of_ops))
+        normal_array = unp.asarray(np.array(list_of_ops))
         op_arr = [Operator.from_label(s) for s in "XYZ"]
         single_op = op_arr[0]
         list_of_arrays = [Array(op) for op in list_of_ops]
-        assert isinstance(to_array(numpy_ops), np.ndarray)
-        assert isinstance(to_array(normal_array), Array)
-        assert isinstance(to_array(op_arr), np.ndarray)
-        assert isinstance(to_array(single_op), np.ndarray)
-        assert isinstance(to_array(list_of_arrays), np.ndarray)
+        assert isinstance(unp.to_dense(numpy_ops), np.ndarray)
+        assert isinstance(unp.to_dense(normal_array), np.ndarray)
+        assert isinstance(unp.to_dense(op_arr), np.ndarray)
+        assert isinstance(unp.to_dense(single_op), np.ndarray)
+        assert isinstance(unp.to_dense(list_of_arrays), np.ndarray)
 
 
-class Test_to_array_Jax(Test_to_array, TestJaxBase):
-    """Jax version of Test_to_array tests."""
+class Test_to_dense_Jax(TestJaxBase, Test_to_dense):
+    """Jax version of Test_to_dense tests."""
 
-    def test_to_array_types(self):
-        """Type conversion tests for to_array with jax backend"""
+    def test_to_dense_types(self):
+        """Type conversion tests for to_dense with jax backend"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
         numpy_ops = np.array(list_of_ops)
-        normal_array = Array(np.array(list_of_ops))
+        normal_array = unp.asarray(np.array(list_of_ops))
         op_arr = [Operator.from_label(s) for s in "XYZ"]
         single_op = op_arr[0]
-        list_of_arrays = [Array(op) for op in list_of_ops]
-        assert isinstance(to_array(numpy_ops), Array)
-        assert isinstance(to_array(normal_array), Array)
-        assert isinstance(to_array(op_arr), Array)
-        assert isinstance(to_array(single_op), Array)
-        assert isinstance(to_array(list_of_arrays), Array)
-
-    def test_to_array_BCOO(self):
+        list_of_arrays = [unp.asarray(op) for op in list_of_ops]
+        assert isinstance(unp.to_dense(numpy_ops), np.ndarray)
+        assert isinstance(unp.to_dense(normal_array), np.ndarray)
+        assert isinstance(unp.to_dense(op_arr), np.ndarray)
+        assert isinstance(unp.to_dense(single_op), np.ndarray)
+        assert isinstance(unp.to_dense(list_of_arrays), np.ndarray)
+
+    def test_to_dense_BCOO(self):
         """Convert BCOO type to array."""
 
         bcoo = jsparse.BCOO.fromdense(np.array([[0.0, 1.0], [1.0, 0.0]]))
-        out = to_array(bcoo)
-        self.assertTrue(isinstance(out, Array))
+        out = unp.to_dense(bcoo)
+        self.assertTrue(isinstance(out, jnp.ndarray))
         self.assertAllClose(out, np.array([[0.0, 1.0], [1.0, 0.0]]))
 
 
-class Test_to_csr(QiskitDynamicsTestCase):
-    """Tests for to_csr."""
+class Test_to_sparse(QiskitDynamicsTestCase, TestNumpyBase):
+    """Tests for to_sparse."""
 
     def test_None_to_None(self):
         """Test that None input returns None."""
-        self.assertTrue(to_csr(None) is None)
+        self.assertTrue(unp.to_sparse(None) is None)
 
-    def test_to_csr_Operator(self):
-        """Tests for to_csr with a single operator"""
+    def test_to_sparse_Operator(self):
+        """Tests for to_sparse with a single operator"""
         op = Operator.from_label("X")
-        self.assertAllCloseSparse(to_csr(op), csr_matrix([[0, 1], [1, 0]]))
+        self.assertAllCloseSparse(unp.to_sparse(op), csr_matrix([[0, 1], [1, 0]]))
 
-    def test_to_csr_nparray(self):
-        """Tests for to_csr with a single numpy array"""
+    def test_to_sparse_nparray(self):
+        """Tests for to_sparse with a single numpy array"""
         nparray = np.array([[0, 1], [1, 0]])
-        self.assertAllCloseSparse(to_csr(nparray), csr_matrix(nparray))
+        self.assertAllCloseSparse(unp.to_sparse(nparray), csr_matrix(nparray))
 
-    def test_to_csr_array(self):
-        """Tests for to_csr with a single sparse matrix"""
+    def test_to_sparse_array(self):
+        """Tests for to_sparse with a single sparse matrix"""
         op = Operator.from_label("X")
         spm = csr_matrix(op)
-        ar = Array(op)
-        self.assertAllCloseSparse(to_csr(ar), spm)
+        ar = unp.asarray(op)
+        self.assertAllCloseSparse(unp.to_sparse(ar), spm)
 
-    def test_to_csr_sparse_matrix(self):
-        """Tests for to_csr with a single sparse matrix"""
+    def test_to_sparse_sparse_matrix(self):
+        """Tests for to_sparse with a single sparse matrix"""
         op = Operator.from_label("X")
         spm = csr_matrix(op)
-        self.assertAllCloseSparse(to_csr(spm), spm)
+        self.assertAllCloseSparse(unp.to_sparse(spm), spm)
 
-    def test_to_csr_Operator_list(self):
-        """Tests for to_csr with a list of operators"""
+    def test_to_sparse_Operator_list(self):
+        """Tests for to_sparse with a list of operators"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
         op_arr = [Operator.from_label(s) for s in "XYZ"]
         sparse_matrices = [csr_matrix(op) for op in list_of_ops]
-        self.assertAllCloseSparse(to_csr(op_arr), sparse_matrices)
+        self.assertAllCloseSparse(unp.to_sparse(op_arr), sparse_matrices)
 
-    def test_to_csr_nparray_list(self):
-        """Tests for to_csr with a list of numpy arrays"""
+    def test_to_sparse_nparray_list(self):
+        """Tests for to_sparse with a list of numpy arrays"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
-        nparray_list = [np.array(op) for op in list_of_ops]
+        nparray_list = [unp.asarray(op) for op in list_of_ops]
         sparse_matrices = [csr_matrix(op) for op in list_of_ops]
-        self.assertAllCloseSparse(to_csr(nparray_list), sparse_matrices)
+        self.assertAllCloseSparse(unp.to_sparse(nparray_list), sparse_matrices)
 
-    def test_to_csr_sparse_matrix_list(self):
-        """Tests for to_csr with a list of sparse matrices"""
+    def test_to_sparse_sparse_matrix_list(self):
+        """Tests for to_sparse with a list of sparse matrices"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
         sparse_matrices = [csr_matrix(op) for op in list_of_ops]
-        self.assertAllCloseSparse(to_csr(sparse_matrices), sparse_matrices)
+        self.assertAllCloseSparse(unp.to_sparse(sparse_matrices), sparse_matrices)
 
-    def test_to_csr_types(self):
-        """Type conversion tests for to_csr"""
+    def test_to_sparse_types(self):
+        """Type conversion tests for to_sparse"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
-        numpy_ops = np.array(list_of_ops)
-        normal_array = Array(np.array(list_of_ops))
+        normal_array = unp.asarray(list_of_ops)
         op_arr = [Operator.from_label(s) for s in "XYZ"]
         single_op = op_arr[0]
-        list_of_arrays = [Array(op) for op in list_of_ops]
+        list_of_arrays = [unp.asarray(op) for op in list_of_ops]
         single_array = list_of_arrays[0]
         sparse_matrices = [csr_matrix(op) for op in list_of_ops]
-        assert isinstance(to_csr(normal_array)[0], csr_matrix)
-        assert isinstance(to_csr(numpy_ops)[0], csr_matrix)
-        assert isinstance(to_csr(op_arr)[0], csr_matrix)
-        assert isinstance(to_csr(single_op), csr_matrix)
-        assert isinstance(to_csr(single_array), csr_matrix)
-        assert isinstance(to_csr(list_of_arrays[0]), csr_matrix)
-        assert isinstance(to_csr(sparse_matrices), Iterable)
-        assert isinstance(to_csr(sparse_matrices)[0], csr_matrix)
+        assert isinstance(unp.to_sparse(normal_array)[0], csr_matrix)
+        assert isinstance(unp.to_sparse(op_arr)[0], csr_matrix)
+        assert isinstance(unp.to_sparse(single_op), csr_matrix)
+        assert isinstance(unp.to_sparse(single_array), csr_matrix)
+        assert isinstance(unp.to_sparse(list_of_arrays[0]), csr_matrix)
+        assert isinstance(unp.to_sparse(sparse_matrices), Iterable)
+        assert isinstance(unp.to_sparse(sparse_matrices)[0], csr_matrix)
 
 
-class Testto_BCOO(QiskitDynamicsTestCase, TestJaxBase):
+class Test_to_BCOO(QiskitDynamicsTestCase, TestJaxBase):
     """Test the to_BCOO function."""
 
     def test_None_to_None(self):
@@ -521,16 +521,16 @@ class Test_to_numeric_matrix_type(QiskitDynamicsTestCase):
     def test_to_numeric_matrix_type(self):
         """Tests for to_numeric_matrix_type"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
-        normal_array = Array(np.array(list_of_ops))
-        list_of_arrays = [Array(op) for op in list_of_ops]
+        normal_array = unp.asarray(list_of_ops)
+        list_of_arrays = [unp.asarray(op) for op in list_of_ops]
         op_arr = [Operator.from_label(s) for s in "XYZ"]
         sparse_matrices = [csr_matrix(op) for op in list_of_ops]
-        self.assertAllClose(to_numeric_matrix_type(list_of_ops), normal_array)
-        self.assertAllClose(to_numeric_matrix_type(list_of_arrays), normal_array)
-        self.assertAllClose(to_numeric_matrix_type(op_arr), list_of_arrays)
+        self.assertAllClose(unp.to_numeric_matrix_type(list_of_ops), normal_array)
+        self.assertAllClose(unp.to_numeric_matrix_type(list_of_arrays), normal_array)
+        self.assertAllClose(unp.to_numeric_matrix_type(op_arr), list_of_arrays)
         for i in range(3):
             self.assertAllCloseSparse(
-                to_numeric_matrix_type(sparse_matrices)[i], sparse_matrices[i]
+                unp.to_numeric_matrix_type(sparse_matrices)[i], sparse_matrices[i]
             )
 
 
@@ -540,10 +540,10 @@ class Test_to_numeric_matrix_type_Jax(QiskitDynamicsTestCase, TestJaxBase):
     def test_to_numeric_matrix_type(self):
         """Tests for to_numeric_matrix_type"""
         list_of_ops = [[[0, 1], [1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]]
-        bcoo = jsparse.BCOO.fromdense(to_array(list_of_ops))
-        bcoo2 = to_numeric_matrix_type(bcoo)
+        bcoo = jsparse.BCOO.fromdense(unp.to_dense(list_of_ops))
+        bcoo2 = unp.to_numeric_matrix_type(bcoo)
         self.assertTrue(type(bcoo2).__name__ == "BCOO")
-        self.assertAllClose(bcoo.todense(), bcoo2.todense())
+        self.assertAllClose(bcoo.todense(), unp.to_dense(bcoo2))
 
 
 class TestTypeUtilsQutip(QiskitDynamicsTestCase, TestQutipBase):
@@ -558,5 +558,5 @@ def test_qutip_conversion(self):
         sparse_matrices = [csr_matrix(op) for op in list_of_ops]
         qutip_qobj = [Qobj(op) for op in list_of_ops]
         self.assertAllCloseSparse(to_numeric_matrix_type(qutip_qobj)[0], sparse_matrices[0])
-        self.assertAllCloseSparse(to_csr(qutip_qobj)[0], sparse_matrices[0])
+        self.assertAllCloseSparse(unp.to_sparse(qutip_qobj)[0], sparse_matrices[0])
         self.assertAllClose(to_array(qutip_qobj)[0], normal_array[0])
diff --git a/tox.ini b/tox.ini
index 6e65c3521..b66d0513d 100644
--- a/tox.ini
+++ b/tox.ini
@@ -15,15 +15,15 @@ commands = stestr run {posargs}
 [testenv:jax]
 deps =
     -r{toxinidir}/requirements-dev.txt
-    jax
-    jaxlib
+    jax<=0.4.6
+    jaxlib<=0.4.6
     diffrax
 
 [testenv:lint]
 deps =
     -r{toxinidir}/requirements-dev.txt
-    jax
-    jaxlib
+    jax<=0.4.6
+    jaxlib<=0.4.6
     diffrax
 commands =
   black --check {posargs} qiskit_dynamics test
@@ -37,8 +37,8 @@ commands = black {posargs} qiskit_dynamics test
 [testenv:docs]
 deps =
     -r{toxinidir}/requirements-dev.txt
-    jax
-    jaxlib
+    jax<=0.4.6
+    jaxlib<=0.4.6
     diffrax
 commands =
   sphinx-build -b html -W {posargs} docs/ docs/_build/html