diff --git a/src/qutip_qip/operations/gateclass.py b/src/qutip_qip/operations/gateclass.py
index 155ba50b..bd38e7b4 100644
--- a/src/qutip_qip/operations/gateclass.py
+++ b/src/qutip_qip/operations/gateclass.py
@@ -3,7 +3,7 @@
 from itertools import product, chain
 from functools import partial, reduce
 from operator import mul
-from typing import Optional
+from typing import Optional, Union, List, Tuple, Dict, Callable, Any
 
 import warnings
 import inspect
@@ -51,12 +51,14 @@
 
 from packaging.version import parse as parse_version
 
+QubitSpecifier = Union[int, Iterable[int]]
+
 if parse_version(qutip.__version__) >= parse_version("5.dev"):
-    is_qutip5 = True
+    is_qutip5: bool = True
 else:
-    is_qutip5 = False
+    is_qutip5: bool = False
 
-__all__ = [
+__all__: List[str] = [
     "Gate",
     "GATE_CLASS_MAP",
     "X",
@@ -96,81 +98,21 @@
     "RZX",
 ]
 
-"""
-.. testsetup::
-
-   import numpy as np
-   np.set_printoptions(5)
-"""
-
 
 class Gate:
-    r"""
-    Base class for a quantum gate,
-    concrete gate classes need to be defined as subclasses.
-
-    Parameters
-    ----------
-    targets : list or int
-        The target qubits fo the gate.
-    controls : list or int
-        The controlling qubits of the gate.
-    arg_value : object
-        Argument value of the gate. It will be saved as an attributes and
-        can be accessed when generating the `:obj:qutip.Qobj`.
-    classical_controls : int or list of int, optional
-        Indices of classical bits to control the unitary operator.
-    control_value : int, optional
-        The decimal value of controlling bits for executing
-        the unitary operator on the target qubits.
-        E.g. if the gate should be executed when the zero-th bit is 1,
-        ``controll_value=1``;
-        If the gate should be executed when the two bits are 1 and 0,
-        ``controll_value=2``.
-    classical_control_value : int, optional
-        The decimal value of controlling classical bits for executing
-        the unitary operator on the target qubits.
-        E.g. if the gate should be executed when the zero-th bit is 1,
-        ``controll_value=1``;
-        If the gate should be executed when the two bits are 1 and 0,
-        ``controll_value=2``.
-        The default is ``2**len(classical_controls)-1``
-        (i.e. all classical controls are 1).
-    arg_label : string
-        Label for the argument, it will be shown in the circuit plot,
-        representing the argument value provided to the gate, e.g,
-        if ``arg_label="\phi" the latex name for the gate in the circuit plot
-        will be ``$U(\phi)$``.
-    name : string, optional
-        The name of the gate. This is kept for backward compatibility
-        to identify different gates.
-        In most cases it is identical to the class name,
-        but that is not guaranteed.
-        It is recommended to use ``isinstance``
-        or ``issubclass`` to identify a gate rather than
-        comparing the name string.
-    style : dict, optional
-        A dictionary of style options for the gate.
-        The options are passed to the `matplotlib` plotter.
-        The default is None.
-    """
-
     def __init__(
         self,
-        name=None,
-        targets=None,
-        controls=None,
-        arg_value=None,
-        control_value=None,
-        classical_controls=None,
+        name: Optional[str] = None,
+        targets: Optional[QubitSpecifier] = None,
+        controls: Optional[QubitSpecifier] = None,
+        arg_value: Optional[Any] = None,
+        control_value: Optional[int] = None,
+        classical_controls: Optional[QubitSpecifier] = None,
         classical_control_value: Optional[int] = None,
-        arg_label=None,
-        style=None,
-        **kwargs,
-    ):
-        """
-        Create a gate with specified parameters.
-        """
+        arg_label: Optional[str] = None,
+        style: Optional[Dict[str, Any]] = None,
+        **kwargs: Any,
+    ) -> None:
         self.name = name if name is not None else self.__class__.__name__
         self.targets = None
         self.controls = None
@@ -218,18 +160,7 @@ def __init__(
             if not all_integer:
                 raise ValueError("Index of a qubit must be an integer")
 
-    def get_all_qubits(self):
-        """
-        Return a list of all qubits that the gate operator
-        acts on.
-        The list concatenates the two lists representing
-        the controls and the targets qubits while retains the order.
-
-        Returns
-        -------
-        targets_list : list of int
-            A list of all qubits, including controls and targets.
-        """
+    def get_all_qubits(self) -> List[int]:
         if self.controls is not None:
             return self.controls + self.targets
         if self.targets is not None:
@@ -238,7 +169,7 @@ def get_all_qubits(self):
             # Special case: the global phase gate
             return []
 
-    def __str__(self):
+    def __str__(self) -> str:
         str_name = (
             "Gate(%s, targets=%s, controls=%s,"
             " classical controls=%s, control_value=%s, classical_control_value=%s)"
@@ -252,22 +183,13 @@ def __str__(self):
         )
         return str_name
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         return str(self)
 
-    def _repr_latex_(self):
+    def _repr_latex_(self) -> str:
         return str(self)
 
-    def _to_qasm(self, qasm_out):
-        """
-        Pipe output of gate signature and application to QasmOutput object.
-
-        Parameters
-        ----------
-        qasm_out: QasmOutput
-            object to store QASM output.
-        """
-
+    def _to_qasm(self, qasm_out: "QasmOutput") -> None:
         qasm_gate = qasm_out.qasm_name(self.name)
 
         if not qasm_gate:
@@ -286,22 +208,7 @@ def _to_qasm(self, qasm_out):
                 )
             )
 
-    def get_compact_qobj(self):
-        """
-        Get the compact :class:`qutip.Qobj` representation of the gate
-        operator, ignoring the controls and targets.
-        In the unitary representation,
-        it always assumes that the first few qubits are controls,
-        then targets.
-
-        Returns
-        -------
-        qobj : :obj:`qutip.Qobj`
-            The compact gate operator as a unitary matrix.
-        """
-        # TODO This will be moved to each sub-class of Gate.
-        # However, one first needs to replace the direct use of Gate in
-        # other modules.
+    def get_compact_qobj(self) -> Qobj:
         if self.name == "RX":
             qobj = rx(self.arg_value)
         elif self.name == "RY":
@@ -374,27 +281,11 @@ def get_compact_qobj(self):
             raise NotImplementedError(f"{self.name} is an unknown gate.")
         return qobj
 
-    def get_qobj(self, num_qubits=None, dims=None):
-        """
-        Get the :class:`qutip.Qobj` representation of the gate operator.
-        The operator is expanded to the full Herbert space according to
-        the controls and targets qubits defined for the gate.
-
-        Parameters
-        ----------
-        num_qubits : int, optional
-            The number of qubits.
-            If not given, use the minimal number of qubits required
-            by the target and control qubits.
-        dims : list, optional
-            A list representing the dimensions of each quantum system.
-            If not given, it is assumed to be an all-qubit system.
-
-        Returns
-        -------
-        qobj : :obj:`qutip.Qobj`
-            The compact gate operator as a unitary matrix.
-        """
+    def get_qobj(
+        self,
+        num_qubits: Optional[int] = None,
+        dims: Optional[List[int]] = None,
+    ) -> Qobj:
         if self.name == "GLOBALPHASE":
             if num_qubits is not None:
                 return globalphase(self.arg_value, num_qubits)
@@ -415,7 +306,7 @@ def get_qobj(self, num_qubits=None, dims=None):
 
 
 class SingleQubitGate(Gate):
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         if self.targets is None or len(self.targets) != 1:
             raise ValueError(
@@ -441,11 +332,11 @@ class X(SingleQubitGate):
      [1. 0.]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"X"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return qutip.sigmax() if not is_qutip5 else qutip.sigmax(dtype="dense")
 
 
@@ -463,11 +354,11 @@ class Y(SingleQubitGate):
      [0.+1.j 0.+0.j]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"Y"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return qutip.sigmay() if not is_qutip5 else qutip.sigmay(dtype="dense")
 
 
@@ -485,11 +376,11 @@ class Z(SingleQubitGate):
      [ 0. -1.]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"Z"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return qutip.sigmaz() if not is_qutip5 else qutip.sigmaz(dtype="dense")
 
 
@@ -507,11 +398,13 @@ class RX(SingleQubitGate):
      [0.     -0.70711j 0.70711+0.j     ]]
     """
 
-    def __init__(self, targets, arg_value, **kwargs):
+    def __init__(
+        self, targets: QubitSpecifier, arg_value: float, **kwargs
+    ) -> None:
         super().__init__(targets=targets, arg_value=arg_value, **kwargs)
         self.latex_str = r"R_x"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return rx(self.arg_value)
 
 
@@ -529,11 +422,13 @@ class RY(SingleQubitGate):
      [ 0.70711  0.70711]]
     """
 
-    def __init__(self, targets, arg_value, **kwargs):
+    def __init__(
+        self, targets: QubitSpecifier, arg_value: float, **kwargs
+    ) -> None:
         super().__init__(targets=targets, arg_value=arg_value, **kwargs)
         self.latex_str = r"R_y"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return ry(self.arg_value)
 
 
@@ -551,11 +446,13 @@ class RZ(SingleQubitGate):
      [0.     +0.j      0.70711+0.70711j]]
     """
 
-    def __init__(self, targets, arg_value, **kwargs):
+    def __init__(
+        self, targets: QubitSpecifier, arg_value: float, **kwargs
+    ) -> None:
         super().__init__(targets=targets, arg_value=arg_value, **kwargs)
         self.latex_str = r"R_z"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return rz(self.arg_value)
 
 
@@ -573,11 +470,11 @@ class H(SingleQubitGate):
      [ 0.70711 -0.70711]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{\rm H}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return snot()
 
 
@@ -599,11 +496,11 @@ class SQRTNOT(SingleQubitGate):
      [0.5-0.5j 0.5+0.5j]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"\sqrt{\rm NOT}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return sqrtnot()
 
 
@@ -621,11 +518,11 @@ class S(SingleQubitGate):
      [0.+0.j 0.+1.j]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{\rm S}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return s_gate()
 
 
@@ -643,11 +540,11 @@ class T(SingleQubitGate):
      [0.     +0.j      0.70711+0.70711j]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{\rm T}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return t_gate()
 
 
@@ -672,11 +569,13 @@ class R(SingleQubitGate):
      [ 0.70711  0.70711]]
     """
 
-    def __init__(self, targets, arg_value=None, **kwargs):
+    def __init__(
+        self, targets: QubitSpecifier, arg_value: Tuple[float, float], **kwargs
+    ) -> None:
         super().__init__(targets=targets, arg_value=arg_value, **kwargs)
         self.latex_str = r"{\rm R}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return qrot(*self.arg_value)
 
 
@@ -697,18 +596,23 @@ class QASMU(SingleQubitGate):
      [ 0.5+0.5j -0.5+0.5j]]
     """
 
-    def __init__(self, targets, arg_value=None, **kwargs):
+    def __init__(
+        self,
+        targets: QubitSpecifier,
+        arg_value: Tuple[float, float, float],
+        **kwargs,
+    ) -> None:
         super().__init__(targets=targets, arg_value=arg_value, **kwargs)
         self.latex_str = r"{\rm QASM-U}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return qasmu_gate(self.arg_value)
 
 
 class TwoQubitGate(Gate):
     """Abstract two-qubit gate."""
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         if len(self.get_all_qubits()) != 2:
             raise ValueError(
@@ -732,11 +636,11 @@ class SWAP(TwoQubitGate):
      [0. 0. 0. 1.]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{\rm SWAP}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return swap()
 
 
@@ -756,11 +660,11 @@ class ISWAP(TwoQubitGate):
      [0.+0.j 0.+0.j 0.+0.j 1.+0.j]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{i}{\rm SWAP}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return iswap()
 
 
@@ -780,11 +684,11 @@ class SQRTSWAP(TwoQubitGate):
      [0. +0.j  0. +0.j  0. +0.j  1. +0.j ]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"\sqrt{\rm SWAP}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return sqrtswap()
 
 
@@ -804,11 +708,11 @@ class SQRTISWAP(TwoQubitGate):
      [0.     +0.j      0.     +0.j      0.     +0.j      1.     +0.j     ]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"\sqrt{{i}\rm SWAP}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return sqrtiswap()
 
 
@@ -837,11 +741,11 @@ class BERKELEY(TwoQubitGate):
      [0.     +0.38268j 0.     +0.j      0.     +0.j      0.92388+0.j     ]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{\rm BERKELEY}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return berkeley()
 
 
@@ -870,11 +774,13 @@ class SWAPALPHA(TwoQubitGate):
      [0. +0.j  0. +0.j  0. +0.j  1. +0.j ]]
     """
 
-    def __init__(self, targets, arg_value, **kwargs):
+    def __init__(
+        self, targets: QubitSpecifier, arg_value: float, **kwargs
+    ) -> None:
         super().__init__(targets=targets, arg_value=arg_value, **kwargs)
         self.latex_str = r"{\rm SWAPALPHA}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return swapalpha(self.arg_value)
 
 
@@ -903,11 +809,13 @@ class MS(TwoQubitGate):
      [0.     -0.70711j 0.     +0.j      0.     +0.j      0.70711+0.j     ]]
     """
 
-    def __init__(self, targets, arg_value, **kwargs):
+    def __init__(
+        self, targets: QubitSpecifier, arg_value: Tuple[float, float], **kwargs
+    ) -> None:
         super().__init__(targets=targets, arg_value=arg_value, **kwargs)
         self.latex_str = r"{\rm MS}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return molmer_sorensen(*self.arg_value)
 
 
@@ -931,11 +839,11 @@ class TOFFOLI(Gate):
      [0. 0. 0. 0. 0. 0. 1. 0.]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{\rm TOFFOLI}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return toffoli()
 
 
@@ -959,18 +867,23 @@ class FREDKIN(Gate):
      [0. 0. 0. 0. 0. 0. 0. 1.]]
     """
 
-    def __init__(self, targets, **kwargs):
+    def __init__(self, targets: QubitSpecifier, **kwargs: Any) -> None:
         super().__init__(targets=targets, **kwargs)
         self.latex_str = r"{\rm FREDKIN}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return fredkin()
 
 
 class ControlledGate(Gate):
     def __init__(
-        self, controls, targets, control_value, target_gate, **kwargs
-    ):
+        self,
+        controls: QubitSpecifier,
+        targets: QubitSpecifier,
+        control_value: Optional[int],
+        target_gate: Callable[..., Gate],
+        **kwargs: Any,
+    ) -> None:
         super().__init__(
             controls=controls,
             targets=targets,
@@ -990,7 +903,7 @@ def __init__(
             targets=self.targets, **self.kwargs
         ).latex_str
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return controlled_gate(
             U=self.target_gate(
                 targets=self.targets, **self.kwargs
@@ -1013,12 +926,18 @@ class _OneControlledGate(ControlledGate, TwoQubitGate):
     and raise an error if it is 0.
     """
 
-    def __init__(self, controls, targets, target_gate, **kwargs):
+    def __init__(
+        self,
+        controls: QubitSpecifier,
+        targets: QubitSpecifier,
+        target_gate: Callable[..., Gate],
+        **kwargs: Any,
+    ) -> None:
         _control_value = kwargs.get("control_value", None)
         if _control_value is not None:
             if _control_value != 1:
                 raise ValueError(
-                    f"{self.__class__.__name__} must has control_value=1"
+                    f"{self.__class__.__name__} must have control_value=1"
                 )
         else:
             kwargs["control_value"] = 1
@@ -1046,7 +965,9 @@ class CNOT(_OneControlledGate):
      [0. 0. 1. 0.]]
     """
 
-    def __init__(self, controls, targets, **kwargs):
+    def __init__(
+        self, controls: QubitSpecifier, targets: QubitSpecifier, **kwargs: Any
+    ) -> None:
         self.target_gate = X
         super().__init__(
             targets=targets,
@@ -1056,7 +977,7 @@ def __init__(self, controls, targets, **kwargs):
         )
         self.latex_str = r"{\rm CNOT}"
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return cnot()
 
 
@@ -1076,7 +997,9 @@ class CZ(_OneControlledGate):
      [ 0.  0.  0. -1.]]
     """
 
-    def __init__(self, controls, targets, **kwargs):
+    def __init__(
+        self, controls: QubitSpecifier, targets: QubitSpecifier, **kwargs: Any
+    ) -> None:
         self.target_gate = Z
         super().__init__(
             targets=targets,
@@ -1085,7 +1008,7 @@ def __init__(self, controls, targets, **kwargs):
             **kwargs,
         )
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return csign()
 
 
@@ -1105,7 +1028,9 @@ class CSIGN(_OneControlledGate):
      [ 0.  0.  0. -1.]]
     """
 
-    def __init__(self, controls, targets, **kwargs):
+    def __init__(
+        self, controls: QubitSpecifier, targets: QubitSpecifier, **kwargs: Any
+    ) -> None:
         self.target_gate = Z
         super().__init__(
             targets=targets,
@@ -1114,7 +1039,7 @@ def __init__(self, controls, targets, **kwargs):
             **kwargs,
         )
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return csign()
 
 
@@ -1144,18 +1069,24 @@ class CPHASE(_OneControlledGate):
     """
 
     def __init__(
-        self, controls, targets, arg_value, control_value=1, **kwargs
-    ):
+        self,
+        controls: QubitSpecifier,
+        targets: QubitSpecifier,
+        arg_value: float,
+        control_value: int = 1,
+        **kwargs: Any,
+    ) -> None:
         self.target_gate = RZ
         super().__init__(
             targets=targets,
             controls=controls,
             arg_value=arg_value,
             target_gate=self.target_gate,
+            control_value=control_value,
             **kwargs,
         )
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         return cphase(self.arg_value).tidyup()
 
 
@@ -1184,7 +1115,9 @@ class RZX(TwoQubitGate):
     [0.+0.j 0.+0.j 0.+1.j 0.+0.j]]
     """
 
-    def __init__(self, targets, arg_value, **kwargs):
+    def __init__(
+        self, targets: QubitSpecifier, arg_value: float, **kwargs: Any
+    ) -> None:
         self.target_gate = RZ
         super().__init__(
             targets=targets,
@@ -1193,7 +1126,7 @@ def __init__(self, targets, arg_value, **kwargs):
             **kwargs,
         )
 
-    def get_compact_qobj(self):
+    def get_compact_qobj(self) -> Qobj:
         theta = self.arg_value
         return Qobj(
             np.array(
diff --git a/src/qutip_qip/operations/gates.py b/src/qutip_qip/operations/gates.py
index 5d078713..87bd58d6 100644
--- a/src/qutip_qip/operations/gates.py
+++ b/src/qutip_qip/operations/gates.py
@@ -4,16 +4,17 @@
 from itertools import product
 from functools import partial, reduce
 from operator import mul
+from typing import Optional, List, Union, Tuple, Generator
 
 import warnings
-import inspect
-from copy import deepcopy
 
 import numpy as np
 import scipy.sparse as sp
 
 import qutip
-from qutip import Qobj, identity, qeye, sigmax, sigmay, sigmaz, tensor, fock_dm
+from qutip import Qobj, identity, qeye, sigmax, sigmay, sigmaz, tensor
+
+QubitSpecifier = Union[int, Iterable[int]]
 
 __all__ = [
     "rx",
@@ -58,18 +59,9 @@
 #
 # Single Qubit Gates
 #
-def _deprecation_warnings_gate_expansion():
-    warnings.warn(
-        "The expansion of output gate matrix is no longer included "
-        "in the gate functions. "
-        "To expand the output `Qobj` or permute the qubits, "
-        "please use expand_operator.",
-        DeprecationWarning,
-        stacklevel=2,
-    )
 
 
-def x_gate(N=None, target=0):
+def x_gate() -> Qobj:
     """Pauli-X gate or sigmax operator.
 
     Returns
@@ -79,13 +71,10 @@ def x_gate(N=None, target=0):
         a single-qubit rotation through pi radians around the x-axis.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(x_gate(), dims=[2] * N, targets=target)
     return sigmax()
 
 
-def y_gate(N=None, target=0):
+def y_gate() -> Qobj:
     """Pauli-Y gate or sigmay operator.
 
     Returns
@@ -95,13 +84,10 @@ def y_gate(N=None, target=0):
         a single-qubit rotation through pi radians around the y-axis.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(y_gate(), dims=[2] * N, targets=target)
     return sigmay()
 
 
-def cy_gate(N=None, control=0, target=1):
+def cy_gate() -> Qobj:
     """Controlled Y gate.
 
     Returns
@@ -110,21 +96,13 @@ def cy_gate(N=None, control=0, target=1):
         Quantum object for operator describing the rotation.
 
     """
-    if (control == 1 and target == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            cy_gate(), dims=[2] * N, targets=(control, target)
-        )
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1j], [0, 0, 1j, 0]],
         dims=[[2, 2], [2, 2]],
     )
 
 
-def z_gate(N=None, target=0):
+def z_gate() -> Qobj:
     """Pauli-Z gate or sigmaz operator.
 
     Returns
@@ -134,13 +112,10 @@ def z_gate(N=None, target=0):
         a single-qubit rotation through pi radians around the z-axis.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(z_gate(), dims=[2] * N, targets=target)
     return sigmaz()
 
 
-def cz_gate(N=None, control=0, target=1):
+def cz_gate() -> Qobj:
     """Controlled Z gate.
 
     Returns
@@ -149,21 +124,13 @@ def cz_gate(N=None, control=0, target=1):
         Quantum object for operator describing the rotation.
 
     """
-    if (control == 1 and target == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            cz_gate(), dims=[2] * N, targets=(control, target)
-        )
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]],
         dims=[[2, 2], [2, 2]],
     )
 
 
-def s_gate(N=None, target=0):
+def s_gate() -> Qobj:
     """Single-qubit rotation also called Phase gate or the Z90 gate.
 
     Returns
@@ -173,13 +140,10 @@ def s_gate(N=None, target=0):
         a 90 degree rotation around the z-axis.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(s_gate(), dims=[2] * N, targets=target)
     return Qobj([[1, 0], [0, 1j]])
 
 
-def cs_gate(N=None, control=0, target=1):
+def cs_gate() -> Qobj:
     """Controlled S gate.
 
     Returns
@@ -188,21 +152,13 @@ def cs_gate(N=None, control=0, target=1):
         Quantum object for operator describing the rotation.
 
     """
-    if (control == 1 and target == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            cs_gate(), dims=[2] * N, targets=(control, target)
-        )
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1j]],
         dims=[[2, 2], [2, 2]],
     )
 
 
-def t_gate(N=None, target=0):
+def t_gate() -> Qobj:
     """Single-qubit rotation related to the S gate by the relationship S=T*T.
 
     Returns
@@ -211,13 +167,10 @@ def t_gate(N=None, target=0):
         Quantum object for operator describing a phase shift of pi/4.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(t_gate(), dims=[2] * N, targets=target)
     return Qobj([[1, 0], [0, np.exp(1j * np.pi / 4)]])
 
 
-def ct_gate(N=None, control=0, target=1):
+def ct_gate() -> Qobj:
     """Controlled T gate.
 
     Returns
@@ -226,14 +179,6 @@ def ct_gate(N=None, control=0, target=1):
         Quantum object for operator describing the rotation.
 
     """
-    if (control == 1 and target == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            ct_gate(), dims=[2] * N, targets=(control, target)
-        )
     return Qobj(
         [
             [1, 0, 0, 0],
@@ -245,18 +190,15 @@ def ct_gate(N=None, control=0, target=1):
     )
 
 
-def rx(phi, N=None, target=0):
+def rx(phi: float) -> Qobj:
     """Single-qubit rotation for operator sigmax with angle phi.
 
     Returns
     -------
-    result : qobj
+    result : :class:`qutip.Qobj`
         Quantum object for operator describing the rotation.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(rx(phi), dims=[2] * N, targets=target)
     return Qobj(
         [
             [np.cos(phi / 2), -1j * np.sin(phi / 2)],
@@ -265,18 +207,15 @@ def rx(phi, N=None, target=0):
     )
 
 
-def ry(phi, N=None, target=0):
+def ry(phi: float) -> Qobj:
     """Single-qubit rotation for operator sigmay with angle phi.
 
     Returns
     -------
-    result : qobj
+    result : :class:`qutip.Qobj`
         Quantum object for operator describing the rotation.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(ry(phi), dims=[2] * N, targets=target)
     return Qobj(
         [
             [np.cos(phi / 2), -np.sin(phi / 2)],
@@ -285,42 +224,36 @@ def ry(phi, N=None, target=0):
     )
 
 
-def rz(phi, N=None, target=0):
+def rz(phi: float) -> Qobj:
     """Single-qubit rotation for operator sigmaz with angle phi.
 
     Returns
     -------
-    result : qobj
+    result : :class:`qutip.Qobj`
         Quantum object for operator describing the rotation.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(rz(phi), dims=[2] * N, targets=target)
     return Qobj([[np.exp(-1j * phi / 2), 0], [0, np.exp(1j * phi / 2)]])
 
 
-def sqrtnot(N=None, target=0):
+def sqrtnot() -> Qobj:
     """Single-qubit square root NOT gate.
 
     Returns
     -------
-    result : qobj
+    result : :class:`qutip.Qobj`
         Quantum object for operator describing the square root NOT gate.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(sqrtnot(), dims=[2] * N, targets=target)
     return Qobj([[0.5 + 0.5j, 0.5 - 0.5j], [0.5 - 0.5j, 0.5 + 0.5j]])
 
 
-def snot(N=None, target=0):
+def snot() -> Qobj:
     """Quantum object representing the SNOT (Hadamard) gate.
 
     Returns
     -------
-    snot_gate : qobj
+    snot_gate : :class:`qutip.Qobj`
         Quantum object representation of SNOT gate.
 
     Examples
@@ -333,13 +266,10 @@ def snot(N=None, target=0):
      [ 0.70710678+0.j -0.70710678+0.j]]
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(snot(), dims=[2] * N, targets=target)
     return 1 / np.sqrt(2.0) * Qobj([[1, 1], [1, -1]])
 
 
-def phasegate(theta, N=None, target=0):
+def phasegate(theta: float) -> Qobj:
     """
     Returns quantum object representing the phase shift gate.
 
@@ -350,7 +280,7 @@ def phasegate(theta, N=None, target=0):
 
     Returns
     -------
-    phase_gate : qobj
+    phase_gate : :class:`qutip.Qobj`
         Quantum object representation of phase shift gate.
 
     Examples
@@ -363,13 +293,10 @@ def phasegate(theta, N=None, target=0):
      [ 0.00000000+0.j          0.70710678+0.70710678j]]
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(phasegate(theta), dims=[2] * N, targets=target)
     return Qobj([[1, 0], [0, np.exp(1.0j * theta)]], dims=[[2], [2]])
 
 
-def qrot(theta, phi, N=None, target=0):
+def qrot(theta: float, phi: float) -> Qobj:
     """
     Single qubit rotation driving by Rabi oscillation with 0 detune.
 
@@ -390,9 +317,6 @@ def qrot(theta, phi, N=None, target=0):
         Quantum object representation of physical qubit rotation under
         a rabi pulse.
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(qrot(theta, phi), dims=[2] * N, targets=target)
     return Qobj(
         [
             [
@@ -407,23 +331,21 @@ def qrot(theta, phi, N=None, target=0):
     )
 
 
-def qasmu_gate(args, N=None, target=0):
+def qasmu_gate(args: Tuple[float, float, float]) -> Qobj:
     """
     QASM U-gate as defined in the OpenQASM standard.
 
     Parameters
     ----------
+    args : tuple
+        Three parameters:
 
-    theta : float
-        The argument supplied to the last RZ rotation.
-    phi : float
-        The argument supplied to the middle RY rotation.
-    gamma : float
-        The argument supplied to the first RZ rotation.
-    N : int
-        Number of qubits in the system.
-    target : int
-        The index of the target qubit.
+        - theta : float
+            The argument supplied to the last RZ rotation.
+        - phi : float
+            The argument supplied to the middle RY rotation.
+        - gamma : float
+            The argument supplied to the first RZ rotation.
 
     Returns
     -------
@@ -433,11 +355,6 @@ def qasmu_gate(args, N=None, target=0):
     """
 
     theta, phi, gamma = args
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            qasmu_gate([theta, phi, gamma]), dims=[2] * N, targets=target
-        )
     return Qobj(rz(phi) * ry(theta) * rz(gamma))
 
 
@@ -446,7 +363,7 @@ def qasmu_gate(args, N=None, target=0):
 #
 
 
-def cphase(theta, N=2, control=0, target=1):
+def cphase(theta: float) -> Qobj:
     """
     Returns quantum object representing the controlled phase shift gate.
 
@@ -455,48 +372,23 @@ def cphase(theta, N=2, control=0, target=1):
     theta : float
         Phase rotation angle.
 
-    N : integer
-        The number of qubits in the target space.
-
-    control : integer
-        The index of the control qubit.
-
-    target : integer
-        The index of the target qubit.
-
     Returns
     -------
-    U : qobj
+    U : :class:`qutip.Qobj`
         Quantum object representation of controlled phase gate.
     """
-    if N != 2 or control != 0 or target != 1:
-        _deprecation_warnings_gate_expansion()
-
-    if N < 1 or target < 0 or control < 0:
-        raise ValueError("Minimum value: N=1, control=0 and target=0")
-
-    if control >= N or target >= N:
-        raise ValueError("control and target need to be smaller than N")
-
-    U_list1 = [identity(2)] * N
-    U_list2 = [identity(2)] * N
-
-    U_list1[control] = fock_dm(2, 1)
-    U_list1[target] = phasegate(theta)
-
-    U_list2[control] = fock_dm(2, 0)
+    mat = np.identity(4, dtype=np.complex128)
+    mat[2:, 2:] = phasegate(theta).full()
+    return Qobj(mat, dims=[[2, 2], [2, 2]])
 
-    U = tensor(U_list1) + tensor(U_list2)
-    return U
 
-
-def cnot(N=None, control=0, target=1):
+def cnot() -> Qobj:
     """
     Quantum object representing the CNOT gate.
 
     Returns
     -------
-    cnot_gate : qobj
+    cnot_gate : :class:`qutip.Qobj`
         Quantum object representation of CNOT gate
 
     Examples
@@ -511,25 +403,19 @@ def cnot(N=None, control=0, target=1):
          [ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]]
 
     """
-    if (control == 1 and target == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(cnot(), dims=[2] * N, targets=(control, target))
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]],
         dims=[[2, 2], [2, 2]],
     )
 
 
-def csign(N=None, control=0, target=1):
+def csign() -> Qobj:
     """
     Quantum object representing the CSIGN gate.
 
     Returns
     -------
-    csign_gate : qobj
+    csign_gate : :class:`qutip.Qobj`
         Quantum object representation of CSIGN gate
 
     Examples
@@ -544,25 +430,19 @@ def csign(N=None, control=0, target=1):
          [ 0.+0.j  0.+0.j  0.+0.j  -1.+0.j]]
 
     """
-    if (control == 1 and target == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(csign(), N, (control, target))
     return Qobj(
         [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]],
         dims=[[2, 2], [2, 2]],
     )
 
 
-def berkeley(N=None, targets=[0, 1]):
+def berkeley() -> Qobj:
     """
     Quantum object representing the Berkeley gate.
 
     Returns
     -------
-    berkeley_gate : qobj
+    berkeley_gate : :class:`qutip.Qobj`
         Quantum object representation of Berkeley gate
 
     Examples
@@ -577,12 +457,6 @@ def berkeley(N=None, targets=[0, 1]):
          [ 0.+sin(pi/8).j  0.+0.j           0.+0.j           cos(pi/8).+0.j]]
 
     """
-    if (targets[0] == 1 and targets[1] == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(berkeley(), N, targets=targets)
     return Qobj(
         [
             [np.cos(np.pi / 8), 0, 0, 1.0j * np.sin(np.pi / 8)],
@@ -594,13 +468,13 @@ def berkeley(N=None, targets=[0, 1]):
     )
 
 
-def swapalpha(alpha, N=None, targets=[0, 1]):
+def swapalpha(alpha: float) -> Qobj:
     """
     Quantum object representing the SWAPalpha gate.
 
     Returns
     -------
-    swapalpha_gate : qobj
+    swapalpha_gate : :class:`qutip.Qobj`
         Quantum object representation of SWAPalpha gate
 
     Examples
@@ -615,12 +489,6 @@ def swapalpha(alpha, N=None, targets=[0, 1]):
      [ 0.+0.j  0.+0.j                    0.+0.j                    1.+0.j]]
 
     """
-    if (targets[0] == 1 and targets[1] == 0) and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(swapalpha(alpha), N, targets=targets)
     return Qobj(
         [
             [1, 0, 0, 0],
@@ -642,12 +510,12 @@ def swapalpha(alpha, N=None, targets=[0, 1]):
     )
 
 
-def swap(N=None, targets=[0, 1]):
+def swap() -> Qobj:
     """Quantum object representing the SWAP gate.
 
     Returns
     -------
-    swap_gate : qobj
+    swap_gate : :class:`qutip.Qobj`
         Quantum object representation of SWAP gate
 
     Examples
@@ -662,24 +530,18 @@ def swap(N=None, targets=[0, 1]):
      [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]
 
     """
-    if targets != [0, 1] and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(swap(), dims=[2] * N, targets=targets)
     return Qobj(
         [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]],
         dims=[[2, 2], [2, 2]],
     )
 
 
-def iswap(N=None, targets=[0, 1]):
+def iswap() -> Qobj:
     """Quantum object representing the iSWAP gate.
 
     Returns
     -------
-    iswap_gate : qobj
+    iswap_gate : :class:`qutip.Qobj`
         Quantum object representation of iSWAP gate
 
     Examples
@@ -693,33 +555,21 @@ def iswap(N=None, targets=[0, 1]):
      [ 0.+0.j  0.+1.j  0.+0.j  0.+0.j]
      [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]
     """
-    if targets != [0, 1] and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(iswap(), dims=[2] * N, targets=targets)
     return Qobj(
         [[1, 0, 0, 0], [0, 0, 1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]],
         dims=[[2, 2], [2, 2]],
     )
 
 
-def sqrtswap(N=None, targets=[0, 1]):
+def sqrtswap() -> Qobj:
     """Quantum object representing the square root SWAP gate.
 
     Returns
     -------
-    sqrtswap_gate : qobj
+    sqrtswap_gate : :class:`qutip.Qobj`
         Quantum object representation of square root SWAP gate
 
     """
-    if targets != [0, 1] and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(sqrtswap(), dims=[2] * N, targets=targets)
     return Qobj(
         np.array(
             [
@@ -733,12 +583,12 @@ def sqrtswap(N=None, targets=[0, 1]):
     )
 
 
-def sqrtiswap(N=None, targets=[0, 1]):
+def sqrtiswap() -> Qobj:
     """Quantum object representing the square root iSWAP gate.
 
     Returns
     -------
-    sqrtiswap_gate : qobj
+    sqrtiswap_gate : :class:`qutip.Qobj`
         Quantum object representation of square root iSWAP gate
 
     Examples
@@ -757,12 +607,6 @@ def sqrtiswap(N=None, targets=[0, 1]):
        0.00000000+0.j          1.00000000+0.j]]
 
     """
-    if targets != [0, 1] and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(sqrtiswap(), N, targets=targets)
     return Qobj(
         np.array(
             [
@@ -776,7 +620,7 @@ def sqrtiswap(N=None, targets=[0, 1]):
     )
 
 
-def molmer_sorensen(theta, phi=0.0, N=None, targets=[0, 1]):
+def molmer_sorensen(theta: float, phi: float = 0.0) -> Qobj:
     """
     Quantum object of a Mølmer–Sørensen gate.
 
@@ -786,25 +630,12 @@ def molmer_sorensen(theta, phi=0.0, N=None, targets=[0, 1]):
         The duration of the interaction pulse.
     phi: float
         Rotation axis. phi = 0 for XX; phi=pi for YY
-    N: int
-        Number of qubits in the system.
-    target: int
-        The indices of the target qubits.
 
     Returns
     -------
     molmer_sorensen_gate : :class:`qutip.Qobj`
         Quantum object representation of the Mølmer–Sørensen gate.
     """
-    if targets != [0, 1] and N is None:
-        N = 2
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            molmer_sorensen(theta, phi), dims=[2] * N, targets=targets
-        )
-
     return Qobj(
         [
             [
@@ -831,12 +662,12 @@ def molmer_sorensen(theta, phi=0.0, N=None, targets=[0, 1]):
 #
 
 
-def fredkin(N=None, control=0, targets=[1, 2]):
+def fredkin() -> Qobj:
     """Quantum object representing the Fredkin gate.
 
     Returns
     -------
-    fredkin_gate : qobj
+    fredkin_gate : :class:`qutip.Qobj`
         Quantum object representation of Fredkin gate.
 
     Examples
@@ -854,14 +685,6 @@ def fredkin(N=None, control=0, targets=[1, 2]):
          [ 0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  0.+0.j  1.+0.j]]
 
     """
-    if [control, targets[0], targets[1]] != [0, 1, 2] and N is None:
-        N = 3
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            fredkin(), dims=[2] * N, targets=(control,) + tuple(targets)
-        )
     return Qobj(
         [
             [1, 0, 0, 0, 0, 0, 0, 0],
@@ -877,12 +700,12 @@ def fredkin(N=None, control=0, targets=[1, 2]):
     )
 
 
-def toffoli(N=None, controls=[0, 1], target=2):
+def toffoli() -> Qobj:
     """Quantum object representing the Toffoli gate.
 
     Returns
     -------
-    toff_gate : qobj
+    toff_gate : :class:`qutip.Qobj`
         Quantum object representation of Toffoli gate.
 
     Examples
@@ -901,14 +724,6 @@ def toffoli(N=None, controls=[0, 1], target=2):
 
 
     """
-    if [controls[0], controls[1], target] != [0, 1, 2] and N is None:
-        N = 3
-
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(
-            toffoli(), dims=[2] * N, targets=tuple(controls) + (target,)
-        )
     return Qobj(
         [
             [1, 0, 0, 0, 0, 0, 0, 0],
@@ -929,28 +744,25 @@ def toffoli(N=None, controls=[0, 1], target=2):
 #
 
 
-def rotation(op, phi, N=None, target=0):
+def rotation(op: Qobj, phi: float) -> Qobj:
     """Single-qubit rotation for operator op with angle phi.
 
     Returns
     -------
-    result : qobj
+    result : ) -> Qobj:
         Quantum object for operator describing the rotation.
 
     """
-    if N is not None:
-        _deprecation_warnings_gate_expansion()
-        return expand_operator(rotation(op, phi), N, target)
     return (-1j * op * phi / 2).expm()
 
 
 def controlled_gate(
-    U,
-    controls=0,
-    targets=1,
-    N=None,
-    control_value=1,
-):
+    U: Qobj,
+    controls: QubitSpecifier = 0,
+    targets: QubitSpecifier = 1,
+    N: Optional[int] = None,
+    control_value: int = 1,
+) -> Qobj:
     """
     Create an N-qubit controlled gate from a single-qubit gate U with the given
     control and target qubits.
@@ -970,7 +782,7 @@ def controlled_gate(
 
     Returns
     -------
-    result : qobj
+    result : :class:`qutip.Qobj`
         Quantum object representing the controlled-U gate.
     """
     # Compatibility
@@ -993,13 +805,13 @@ def controlled_gate(
     result = Qobj(result, dims=[[2] * (num_controls + num_targets)] * 2)
 
     # Expand it to N qubits and permute qubits labelling
-    if controls + targets == list(range(N)):
+    if set(controls + targets) == set(range(N)):
         return result
     else:
         return expand_operator(result, N, targets=controls + targets)
 
 
-def globalphase(theta, N=1):
+def globalphase(theta: float, N: int = 1) -> Qobj:
     """
     Returns quantum object representing the global phase shift gate.
 
@@ -1010,7 +822,7 @@ def globalphase(theta, N=1):
 
     Returns
     -------
-    phase_gate : qobj
+    phase_gate : :class:`qutip.Qobj`
         Quantum object representation of global phase shift gate.
 
     Examples
@@ -1034,7 +846,7 @@ def globalphase(theta, N=1):
 #
 
 
-def _hamming_distance(x, bits=32):
+def _hamming_distance(x: int, bits: int = 32) -> int:
     """
     Calculate the bit-wise Hamming distance of x from 0: That is, the number
     1s in the integer x.
@@ -1046,12 +858,12 @@ def _hamming_distance(x, bits=32):
     return tot
 
 
-def hadamard_transform(N=1):
+def hadamard_transform(N: int = 1) -> Qobj:
     """Quantum object representing the N-qubit Hadamard gate.
 
     Returns
     -------
-    q : qobj
+    q : :class:`qutip.Qobj`
         Quantum object representation of the N-qubit Hadamard gate.
 
     """
@@ -1061,7 +873,7 @@ def hadamard_transform(N=1):
     return tensor([H] * N)
 
 
-def _powers(op, N):
+def _powers(op: Qobj, N: int) -> Generator[Qobj, None, None]:
     """
     Generator that yields powers of an operator `op`,
     through to `N`.
@@ -1074,22 +886,12 @@ def _powers(op, N):
         yield acc
 
 
-def qubit_clifford_group(N=None, target=0):
+def qubit_clifford_group() -> Generator[Qobj, None, None]:
     """
     Generates the Clifford group on a single qubit,
     using the presentation of the group given by Ross and Selinger
     (http://www.mathstat.dal.ca/~selinger/newsynth/).
 
-    Parameters
-    ----------
-
-    N : int or None
-        Number of qubits on which each operator is to be defined
-        (default: 1).
-    target : int
-        Index of the target qubit on which the single-qubit
-        Clifford operators are to act.
-
     Yields
     ------
 
@@ -1125,10 +927,7 @@ def qubit_clifford_group(N=None, target=0):
         # partial(reduce, mul) acting on the tuple yields E**i * X**j * S**k.
 
         # Finally, we optionally expand the gate.
-        if N is not None:
-            yield expand_operator(op, N, target)
-        else:
-            yield op
+        yield op
 
 
 #
@@ -1136,7 +935,11 @@ def qubit_clifford_group(N=None, target=0):
 #
 
 
-def _check_oper_dims(oper, dims=None, targets=None):
+def _check_oper_dims(
+    oper: Qobj,
+    dims: Optional[List[int]] = None,
+    targets: Optional[QubitSpecifier] = None,
+) -> None:
     """
     Check if the given operator is valid.
 
@@ -1166,7 +969,11 @@ def _check_oper_dims(oper, dims=None, targets=None):
             )
 
 
-def _targets_to_list(targets, oper=None, N=None):
+def _targets_to_list(
+    targets: Optional[QubitSpecifier],
+    oper: Optional[Qobj] = None,
+    N: Optional[int] = None,
+) -> List[int]:
     """
     transform targets to a list and check validity.
 
@@ -1206,8 +1013,13 @@ def _targets_to_list(targets, oper=None, N=None):
 
 
 def expand_operator(
-    oper, N=None, targets=None, dims=None, cyclic_permutation=False, dtype=None
-):
+    oper: Qobj,
+    N: Optional[int] = None,
+    targets: Optional[QubitSpecifier] = None,
+    dims: Optional[List[int]] = None,
+    cyclic_permutation: bool = False,
+    dtype: Optional[str] = None,
+) -> Union[Qobj, List[Qobj]]:
     """
     Expand an operator to one that acts on a system with desired dimensions.
 
@@ -1330,8 +1142,11 @@ def expand_operator(
 
 
 def gate_sequence_product(
-    U_list, left_to_right=True, inds_list=None, expand=False
-):
+    U_list: List[Qobj],
+    left_to_right: bool = True,
+    inds_list: Optional[List[List[int]]] = None,
+    expand: bool = False,
+) -> Union[Qobj, Tuple[Qobj, List[int]]]:
     """
     Calculate the overall unitary matrix for a given list of unitary operations.
 
@@ -1352,7 +1167,7 @@ def gate_sequence_product(
 
     Returns
     -------
-    U_overall : qobj
+    U_overall : :class:`qutip.Qobj`
         Unitary matrix corresponding to U_list.
 
     overall_inds : list of int, optional
diff --git a/src/qutip_qip/operations/measurement.py b/src/qutip_qip/operations/measurement.py
index 77461577..35fe7344 100644
--- a/src/qutip_qip/operations/measurement.py
+++ b/src/qutip_qip/operations/measurement.py
@@ -1,6 +1,6 @@
+from typing import Union, List, Optional
 from collections.abc import Iterable
 import numbers
-import os
 
 import numpy as np
 
@@ -9,6 +9,7 @@
 from qutip.measurement import measurement_statistics
 from .gates import expand_operator
 
+QubitSpecifier = Union[int, Iterable[int]]
 
 __all__ = ["Measurement"]
 
@@ -29,7 +30,13 @@ class Measurement:
         classical register of the circuit.
     """
 
-    def __init__(self, name, targets=None, index=None, classical_store=None):
+    def __init__(
+        self,
+        name: str,
+        targets: QubitSpecifier = None,
+        index: Optional[int] = None,
+        classical_store: Optional[int] = None,
+    ) -> None:
         """
         Create a measurement with specified parameters.
         """
@@ -52,7 +59,9 @@ def __init__(self, name, targets=None, index=None, classical_store=None):
                 if not all_integer:
                     raise ValueError("Index of a qubit must be an integer")
 
-    def measurement_comp_basis(self, state):
+    def measurement_comp_basis(
+        self, state: Qobj
+    ) -> tuple[List[Optional[Qobj]], List[float]]:
         """
         Measures a particular qubit (determined by the target)
         whose ket vector/ density matrix is specified in the
@@ -106,7 +115,7 @@ def measurement_comp_basis(self, state):
         ]
         return states, probabilities
 
-    def __str__(self):
+    def __str__(self) -> str:
         str_name = ("Measurement(%s, target=%s, classical_store=%s)") % (
             self.name,
             self.targets,
@@ -114,13 +123,13 @@ def __str__(self):
         )
         return str_name
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         return str(self)
 
-    def _repr_latex_(self):
+    def _repr_latex_(self) -> str:
         return str(self)
 
-    def _to_qasm(self, qasm_out):
+    def _to_qasm(self, qasm_out) -> None:
         """
         Pipe output of measurement to QasmOutput object.
 
diff --git a/tests/test_gates.py b/tests/test_gates.py
index 7a2ac6e0..203418e0 100644
--- a/tests/test_gates.py
+++ b/tests/test_gates.py
@@ -177,82 +177,6 @@ def test_gate_normalises_pauli_group(self, gate):
         assert len(pauli_gates) == 0
 
 
-@pytest.mark.filterwarnings('ignore::DeprecationWarning')
-class TestGateExpansion:
-    """
-    Test that gates act correctly when supplied with controls and targets, i.e.
-    that they pick out the correct qubits to act on, the action is correct, and
-    all other qubits are untouched.
-    """
-    n_qubits = 5
-
-    @pytest.mark.parametrize(["gate", "n_angles"], [
-        pytest.param(gates.rx, 1, id="Rx"),
-        pytest.param(gates.ry, 1, id="Ry"),
-        pytest.param(gates.rz, 1, id="Rz"),
-        pytest.param(gates.x_gate, 0, id="X"),
-        pytest.param(gates.y_gate, 0, id="Y"),
-        pytest.param(gates.z_gate, 0, id="Z"),
-        pytest.param(gates.s_gate, 0, id="S"),
-        pytest.param(gates.t_gate, 0, id="T"),
-        pytest.param(gates.phasegate, 1, id="phase"),
-        pytest.param(gates.qrot, 2, id="Rabi rotation"),
-    ])
-    def test_single_qubit_rotation(self, gate, n_angles):
-        base = qutip.rand_ket(2)
-        angles = np.random.rand(n_angles) * 2*np.pi
-        applied = gate(*angles) * base
-        random = [qutip.rand_ket(2) for _ in [None]*(self.n_qubits - 1)]
-        for target in range(self.n_qubits):
-            start = qutip.tensor(random[:target] + [base] + random[target:])
-            test = gate(*angles, self.n_qubits, target) * start
-            expected = qutip.tensor(random[:target] + [applied]
-                                    + random[target:])
-            assert _infidelity(test, expected) < 1e-12
-
-    @pytest.mark.parametrize(['gate', 'n_controls'], [
-        pytest.param(gates.cnot, 1, id="cnot"),
-        pytest.param(gates.cy_gate, 1, id="cY"),
-        pytest.param(gates.cz_gate, 1, id="cZ"),
-        pytest.param(gates.cs_gate, 1, id="cS"),
-        pytest.param(gates.ct_gate, 1, id="cT"),
-        pytest.param(gates.swap, 0, id="swap"),
-        pytest.param(gates.iswap, 0, id="iswap"),
-        pytest.param(gates.sqrtswap, 0, id="sqrt(swap)"),
-        pytest.param(functools.partial(gates.molmer_sorensen, 0.5*np.pi, 0.), 0,
-                     id="Molmer-Sorensen")
-    ])
-    def test_two_qubit(self, gate, n_controls):
-        targets = [qutip.rand_ket(2) for _ in [None]*2]
-        others = [qutip.rand_ket(2) for _ in [None]*self.n_qubits]
-        reference = gate() * qutip.tensor(*targets)
-        for q1, q2 in itertools.permutations(range(self.n_qubits), 2):
-            qubits = others.copy()
-            qubits[q1], qubits[q2] = targets
-            args = [[q1, q2]] if n_controls == 0 else [q1, q2]
-            test = gate(self.n_qubits, *args) * qutip.tensor(*qubits)
-            expected = _tensor_with_entanglement(qubits, reference, [q1, q2])
-            assert _infidelity(test, expected) < 1e-12
-
-    @pytest.mark.parametrize(['gate', 'n_controls'], [
-        pytest.param(gates.fredkin, 1, id="Fredkin"),
-        pytest.param(gates.toffoli, 2, id="Toffoli"),
-        pytest.param(_make_random_three_qubit_gate(), 2, id="random"),
-    ])
-    def test_three_qubit(self, gate, n_controls):
-        targets = [qutip.rand_ket(2) for _ in [None]*3]
-        others = [qutip.rand_ket(2) for _ in [None]*self.n_qubits]
-        reference = gate() * qutip.tensor(targets)
-        for q1, q2, q3 in itertools.permutations(range(self.n_qubits), 3):
-            qubits = others.copy()
-            qubits[q1], qubits[q2], qubits[q3] = targets
-            args = [q1, [q2, q3]] if n_controls == 1 else [[q1, q2], q3]
-            test = gate(self.n_qubits, *args) * qutip.tensor(*qubits)
-            expected = _tensor_with_entanglement(qubits, reference,
-                                                 [q1, q2, q3])
-            assert _infidelity(test, expected) < 1e-12
-
-
 class Test_expand_operator:
     # Conceptually, a lot of these tests are complete duplicates of
     # `TestGateExpansion`, except that they explicitly target an underlying