def apply_op(self, op: Any, axes: Sequence[int], prng: np.random.RandomState):
        """Applies a unitary operation, mutating the object to represent the new state.

        op:
            The operation that mutates the object. Note that currently, only 1-
            and 2- qubit operations are currently supported.
        """

        old_inds = tuple(map(self.i_str, axes))
        new_inds = tuple(['new_' + old_ind for old_ind in old_inds])

        if protocols.has_unitary(op):
            U = protocols.unitary(op)
        else:
            mixtures = protocols.mixture(op)
            mixture_idx = int(prng.choice(len(mixtures), p=[mixture[0] for mixture in mixtures]))
            U = mixtures[mixture_idx][1]
        U = qtn.Tensor(
            U.reshape([self._qid_shape[axis] for axis in axes] * 2), inds=(new_inds + old_inds)
        )

        # TODO(tonybruguier): Explore using the Quimb's tensor network natively.

        if len(axes) == 1:
            n = self._grouping[axes[0]]

            self._M[n] = (U @ self._M[n]).reindex({new_inds[0]: old_inds[0]})
        elif len(axes) == 2:
            n, p = [self._grouping[axis] for axis in axes]

            if n == p:
                self._M[n] = (U @ self._M[n]).reindex(
                    {new_inds[0]: old_inds[0], new_inds[1]: old_inds[1]}
                )
            else:
                # This is the index on which we do the contraction. We need to add it iff it's
                # the first time that we do the joining for that specific pair.
                mu_ind = self.mu_str(n, p)
                if mu_ind not in self._M[n].inds:
                    self._M[n].new_ind(mu_ind)
                if mu_ind not in self._M[p].inds:
                    self._M[p].new_ind(mu_ind)

                T = U @ self._M[n] @ self._M[p]

                left_inds = tuple(set(T.inds) & set(self._M[n].inds)) + (new_inds[0],)
                X, Y = T.split(
                    left_inds,
                    method=self._simulation_options.method,
                    max_bond=self._simulation_options.max_bond,
                    cutoff=self._simulation_options.cutoff,
                    cutoff_mode=self._simulation_options.cutoff_mode,
                    get='tensors',
                    absorb='both',
                    bond_ind=mu_ind,
                )

                # Equations (13), (14), and (15):
                # TODO(tonybruguier): When Quimb 2.0.0 is released, the split()
                # function should have a 'renorm' that, when set to None, will
                # allow to compute e_n exactly as:
                # np.sum(abs((X @ Y).data) ** 2).real / np.sum(abs(T) ** 2).real
                #
                # The renormalization would then have to be done manually.
                #
                # However, for now, e_n are just the estimated value.
                e_n = self._simulation_options.cutoff
                self._estimated_gate_error_list.append(e_n)

                self._M[n] = X.reindex({new_inds[0]: old_inds[0]})
                self._M[p] = Y.reindex({new_inds[1]: old_inds[1]})
        else:
            # NOTE(tonybruguier): There could be a way to handle higher orders. I think this could
            # involve HOSVDs:
            # https://en.wikipedia.org/wiki/Higher-order_singular_value_decomposition
            #
            # TODO(tonybruguier): Evaluate whether it's even useful to implement and learn more
            # about HOSVDs.
            raise ValueError('Can only handle 1 and 2 qubit operations')
        return True
Applies a unitary operation, mutating the object to represent the new state.

        op:
            The operation that mutates the object. Note that currently, only 1-
            and 2- qubit operations are currently supported.
        
cirq-core/cirq/contrib/quimb/mps_simulator.py