diff --git a/cirq-core/cirq/__init__.py b/cirq-core/cirq/__init__.py
index 39490bdc8a5..f2019aafcd1 100644
--- a/cirq-core/cirq/__init__.py
+++ b/cirq-core/cirq/__init__.py
@@ -516,6 +516,7 @@
     has_kraus,
     has_mixture,
     has_stabilizer_effect,
+    has_superoperator,
     has_unitary,
     inverse,
     is_measurement,
@@ -545,6 +546,7 @@
     read_json,
     resolve_parameters,
     resolve_parameters_once,
+    superoperator,
     SerializableByKey,
     SupportsActOn,
     SupportsActOnQubits,
@@ -571,6 +573,7 @@
     SupportsQasm,
     SupportsQasmWithArgs,
     SupportsQasmWithArgsAndQubits,
+    SupportsSuperoperator,
     SupportsTraceDistanceBound,
     SupportsUnitary,
     to_json_gzip,
diff --git a/cirq-core/cirq/protocols/__init__.py b/cirq-core/cirq/protocols/__init__.py
index a9cb6e2f517..961e85064b6 100644
--- a/cirq-core/cirq/protocols/__init__.py
+++ b/cirq-core/cirq/protocols/__init__.py
@@ -153,6 +153,11 @@
     SupportsExplicitQidShape,
     SupportsExplicitNumQubits,
 )
+from cirq.protocols.superoperator_protocol import (
+    has_superoperator,
+    superoperator,
+    SupportsSuperoperator,
+)
 from cirq.protocols.unitary_protocol import (
     SupportsUnitary,
     unitary,
diff --git a/cirq-core/cirq/protocols/json_test_data/spec.py b/cirq-core/cirq/protocols/json_test_data/spec.py
index 897523a23a0..2f78ebf2398 100644
--- a/cirq-core/cirq/protocols/json_test_data/spec.py
+++ b/cirq-core/cirq/protocols/json_test_data/spec.py
@@ -157,6 +157,7 @@
         'SupportsQasm',
         'SupportsQasmWithArgs',
         'SupportsQasmWithArgsAndQubits',
+        'SupportsSuperoperator',
         'SupportsTraceDistanceBound',
         'SupportsUnitary',
         # mypy types:
diff --git a/cirq-core/cirq/protocols/superoperator_protocol.py b/cirq-core/cirq/protocols/superoperator_protocol.py
new file mode 100644
index 00000000000..dbe73dfee89
--- /dev/null
+++ b/cirq-core/cirq/protocols/superoperator_protocol.py
@@ -0,0 +1,162 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Protocol and methods for obtaining matrix representation of quantum channels."""
+
+from typing import Any, Tuple, TypeVar, Union
+from typing_extensions import Protocol
+
+import numpy as np
+
+from cirq._doc import doc_private
+from cirq.protocols import kraus_protocol
+from cirq.qis.channels import kraus_to_superoperator
+from cirq.type_workarounds import NotImplementedType
+
+# This is a special indicator value used by the superoperator protocol to determine whether
+# or not the caller provided a 'default' argument. It must be of type np.ndarray to ensure
+# the method has the correct type signature in that case. It is checked for using `is`, so
+# it won't have a false positive if the user provides a different np.array([]) value.
+RaiseTypeErrorIfNotProvided = np.array([])
+
+TDefault = TypeVar('TDefault')
+
+
+class SupportsSuperoperator(Protocol):
+    """An object that may be described using a superoperator."""
+
+    @doc_private
+    def _superoperator_(self) -> Union[np.ndarray, NotImplementedType]:
+        """Matrix of the superoperator corresponding to this quantum channel.
+
+        Every quantum channel is a linear map between the vector spaces of input and output
+        density matrices and therefore admits a representation as a matrix. The matrix depends
+        on the choice of bases in the input and output vector spaces. This function returns
+        the matrix of self in the standard basis. Assuming self acts on density matrices of
+        a d-dimensional system, the matrix of the superoperator has d**2 rows and d**2 columns
+        corresponding to the fact that a vectorized density matrix is a d**2-dimensional vector.
+
+        Note that this is different from the Choi matrix. This is easiest to notice in the
+        general case of a quantum channel mapping density matrices of a d1-dimensional system
+        to density matrices of a d2-dimensional system. The Choi matrix is then a square matrix
+        with d1*d2 rows and d1*d2 columns. On the other hand, the superoperator matrix has
+        d2*d2 rows and d1*d1 columns. Note that this implementation assumes d1=d2, so the
+        Choi matrix and the superoperator matrix have the same shape. Nevertheless, they are
+        generally different matrices in this case, too.
+
+        Returns:
+            Superoperator matrix of the channel, or NotImplemented if there is no such
+            matrix.
+        """
+
+    @doc_private
+    def _has_superoperator_(self) -> bool:
+        """Whether this value has a superoperator representation.
+
+        This method is used by the `cirq.has_superoperator` protocol. If this
+        method is not present, or returns NotImplemented, the protocol tries
+        to find the superoperator via other means and returns True if any of
+        the available fallback avenues succeed. Otherwise, `cirq.has_superoperator`
+        returns False.
+
+        Returns:
+            True if the value has a superoperator representation, False otherwise.
+        """
+
+
+def superoperator(
+    val: Any, default: Any = RaiseTypeErrorIfNotProvided
+) -> Union[np.ndarray, TDefault]:
+    """Matrix of the superoperator corresponding to this quantum channel.
+
+    Every quantum channel is a linear map between the vector spaces of input and output
+    density matrices and therefore admits a representation as a matrix. The matrix depends
+    on the choice of bases in the input and output vector spaces. This function returns
+    the matrix of self in the standard basis. Assuming self acts on density matrices of
+    a d-dimensional system, the matrix of the superoperator has d**2 rows and d**2 columns
+    corresponding to the fact that a vectorized density matrix is a d**2-dimensional vector.
+
+    Note that this is different from the Choi matrix. This is easiest to notice in the
+    general case of a quantum channel mapping density matrices of a d1-dimensional system
+    to density matrices of a d2-dimensional system. The Choi matrix is then a square matrix
+    with d1*d2 rows and d1*d2 columns. On the other hand, the superoperator matrix has
+    d2*d2 rows and d1*d1 columns. Note that this implementation assumes d1=d2, so the
+    Choi matrix and the superoperator matrix have the same shape. Nevertheless, they are
+    generally different matrices in this case, too.
+
+    Args:
+        val: The object whose superoperator representation is to be found.
+        default: Determines the fallback behavior when no superoperator representation
+            is found for `val`. If `default` is not set, a TypeError is raised. If
+            default is set to a value, that value is returned.
+
+    Returns:
+        If `val` has a `_superoperator_` method and its result is not NotImplemented,
+        that result is returned. Otherwise, the function attempts to obtain the Kraus
+        representation and compute the superoperator from it. If this is unsuccessful
+        and a default value was specified, the default value is returned.
+
+    Raises:
+        TypeError: `val` doesn't have a _superoperator_, kraus_, _mixture_ or _unitary_
+            method (or they returned NotImplemented) and also no default value was
+            specified.
+    """
+    so_getter = getattr(val, '_superoperator_', None)
+    so_result = NotImplemented if so_getter is None else so_getter()
+    if so_result is not NotImplemented:
+        return so_result
+
+    kraus_result = kraus_protocol.kraus(val, default=None)
+    if kraus_result is not None:
+        return kraus_to_superoperator(kraus_result)
+
+    if default is not RaiseTypeErrorIfNotProvided:
+        return default
+
+    if so_getter is None:
+        raise TypeError(
+            "object of type '{}' has no _superoperator_ method and we could find no Kraus "
+            "representation for it.".format(type(val))
+        )
+
+    raise TypeError(
+        "object of type '{}' does have a _superoperator_ method, but it returned "
+        "NotImplemented and we could find no Kraus representation for it.".format(type(val))
+    )
+
+
+def has_superoperator(val: Any, *, allow_decompose: bool = True) -> bool:
+    """Returns whether val has a superoperator representation.
+
+    Args:
+        val: The value to check.
+        allow_decompose: Used by internal methods to stop redundant decompositions from
+            being performed, see `cirq.has_kraus` for more details. Defaults to True.
+            When False, the decomposition strategy for determining the result is skipped.
+
+    Returns:
+        If `val` has a `_has_superoperator_` method and its result is not NotImplemented,
+        that result is returned. Otherwise, if `cirq.has_kraus` protocol returns True then
+        True is returned. Finally, if `cirq.superoperator` returns a non-default value
+        then True is returned. Otherwise, False is returned.
+    """
+    so_getter = getattr(val, '_has_superoperator_', None)
+    so_result = NotImplemented if so_getter is None else so_getter()
+    if so_result is not NotImplemented:
+        return so_result
+
+    if kraus_protocol.has_kraus(val, allow_decompose=allow_decompose):
+        return True
+
+    return superoperator(val, default=None) is not None
diff --git a/cirq-core/cirq/protocols/superoperator_protocol_test.py b/cirq-core/cirq/protocols/superoperator_protocol_test.py
new file mode 100644
index 00000000000..a0340053c95
--- /dev/null
+++ b/cirq-core/cirq/protocols/superoperator_protocol_test.py
@@ -0,0 +1,242 @@
+# Copyright 2021 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Tests for superoperator_protocol.py."""
+
+from typing import Sequence, Tuple
+
+import numpy as np
+import pytest
+
+import cirq
+
+LOCAL_DEFAULT = np.array([])
+
+
+def test_superoperator_no_methods():
+    class NoMethod:
+        pass
+
+    with pytest.raises(TypeError, match="no _superoperator_ method"):
+        _ = cirq.superoperator(NoMethod())
+
+    assert cirq.superoperator(NoMethod(), None) is None
+    assert cirq.superoperator(NoMethod(), NotImplemented) is NotImplemented
+    assert cirq.superoperator(NoMethod(), LOCAL_DEFAULT) is LOCAL_DEFAULT
+    assert np.all(cirq.superoperator(NoMethod(), np.eye(4)) == np.eye(4))
+
+    assert not cirq.has_superoperator(NoMethod())
+
+
+def test_superoperator_returns_not_implemented():
+    class ReturnsNotImplemented:
+        def _superoperator_(self):
+            return NotImplemented
+
+    with pytest.raises(TypeError, match="returned NotImplemented"):
+        _ = cirq.superoperator(ReturnsNotImplemented())
+
+    assert cirq.superoperator(ReturnsNotImplemented(), None) is None
+    assert cirq.superoperator(ReturnsNotImplemented(), NotImplemented) is NotImplemented
+    assert cirq.superoperator(ReturnsNotImplemented(), LOCAL_DEFAULT) is LOCAL_DEFAULT
+    assert np.all(cirq.superoperator(ReturnsNotImplemented(), np.eye(4)) == np.eye(4))
+
+    assert not cirq.has_superoperator(ReturnsNotImplemented())
+
+
+def test_kraus_returns_not_implemented():
+    class ReturnsNotImplemented:
+        def _kraus_(self):
+            return NotImplemented
+
+    with pytest.raises(TypeError, match="no Kraus representation"):
+        _ = cirq.superoperator(ReturnsNotImplemented())
+
+    assert cirq.superoperator(ReturnsNotImplemented(), None) is None
+    assert cirq.superoperator(ReturnsNotImplemented(), NotImplemented) is NotImplemented
+    assert cirq.superoperator(ReturnsNotImplemented(), LOCAL_DEFAULT) is LOCAL_DEFAULT
+    assert np.all(cirq.superoperator(ReturnsNotImplemented(), np.eye(4)) == np.eye(4))
+
+    assert not cirq.has_superoperator(ReturnsNotImplemented())
+
+
+def test_explicit_superoperator():
+    s = np.array([[0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]])
+
+    class ReturnsSuperoperator:
+        def _superoperator_(self) -> np.ndarray:
+            return s
+
+    assert cirq.superoperator(ReturnsSuperoperator()) is s
+    assert cirq.superoperator(ReturnsSuperoperator(), None) is s
+    assert cirq.superoperator(ReturnsSuperoperator(), NotImplemented) is s
+    assert cirq.superoperator(ReturnsSuperoperator(), np.eye(4)) is s
+    assert cirq.superoperator(ReturnsSuperoperator(), LOCAL_DEFAULT) is s
+
+    assert cirq.has_superoperator(ReturnsSuperoperator())
+
+
+def test_superoperator_fallback_to_kraus():
+    g = 0.1
+    k = (np.diag([1, np.sqrt(1 - g)]), np.array([[0, np.sqrt(g)], [0, 0]]))
+
+    class ReturnsKraus:
+        def _kraus_(self) -> Sequence[np.ndarray]:
+            return k
+
+    s = np.array(
+        [[1, 0, 0, g], [0, np.sqrt(1 - g), 0, 0], [0, 0, np.sqrt(1 - g), 0], [0, 0, 0, 1 - g]]
+    )
+
+    assert np.allclose(cirq.superoperator(ReturnsKraus()), s)
+    assert np.allclose(cirq.superoperator(ReturnsKraus(), None), s)
+    assert np.allclose(cirq.superoperator(ReturnsKraus(), NotImplemented), s)
+    assert np.allclose(cirq.superoperator(ReturnsKraus(), (1,)), s)
+    assert np.allclose(cirq.superoperator(ReturnsKraus(), LOCAL_DEFAULT), s)
+
+    assert cirq.has_superoperator(ReturnsKraus())
+
+
+def test_superoperator_fallback_to_mixture():
+    x, y, z = [cirq.unitary(g) for g in (cirq.X, cirq.Y, cirq.Z)]
+    m = ((0.2, x), (0.3, y), (0.5, z))
+
+    class ReturnsMixture:
+        def _mixture_(self) -> Sequence[Tuple[float, np.ndarray]]:
+            return m
+
+    s = 0.2 * np.kron(x, x) + 0.3 * np.kron(y, y.conj()) + 0.5 * np.kron(z, z)
+
+    assert np.allclose(cirq.superoperator(ReturnsMixture()), s)
+    assert np.allclose(cirq.superoperator(ReturnsMixture(), None), s)
+    assert np.allclose(cirq.superoperator(ReturnsMixture(), NotImplemented), s)
+    assert np.allclose(cirq.superoperator(ReturnsMixture(), (1,)), s)
+    assert np.allclose(cirq.superoperator(ReturnsMixture(), LOCAL_DEFAULT), s)
+
+    assert cirq.has_superoperator(ReturnsMixture())
+
+
+def test_superoperator_fallback_to_unitary():
+    u = cirq.unitary(cirq.Y)
+
+    class ReturnsUnitary:
+        def _unitary_(self) -> np.ndarray:
+            return u
+
+    s = np.kron(u, u.conj())
+
+    assert np.allclose(cirq.superoperator(ReturnsUnitary()), s)
+    assert np.allclose(cirq.superoperator(ReturnsUnitary(), None), s)
+    assert np.allclose(cirq.superoperator(ReturnsUnitary(), NotImplemented), s)
+    assert np.allclose(cirq.superoperator(ReturnsUnitary(), (1,)), s)
+    assert np.allclose(cirq.superoperator(ReturnsUnitary(), LOCAL_DEFAULT), s)
+
+    assert cirq.has_superoperator(ReturnsUnitary())
+
+
+class HasSuperoperator:
+    def _has_superoperator_(self) -> bool:
+        return True
+
+
+class HasKraus(cirq.SingleQubitGate):
+    def _has_kraus_(self) -> bool:
+        return True
+
+
+class HasMixture(cirq.SingleQubitGate):
+    def _has_mixture_(self) -> bool:
+        return True
+
+
+class HasUnitary(cirq.SingleQubitGate):
+    def _has_unitary_(self) -> bool:
+        return True
+
+
+class HasKrausWhenDecomposed(cirq.SingleQubitGate):
+    def __init__(self, decomposed_cls):
+        self.decomposed_cls = decomposed_cls
+
+    def _decompose_(self, qubits):
+        return [self.decomposed_cls().on(q) for q in qubits]
+
+
+@pytest.mark.parametrize('cls', [HasSuperoperator, HasKraus, HasMixture, HasUnitary])
+def test_has_superoperator(cls):
+    assert cirq.has_superoperator(cls())
+
+
+@pytest.mark.parametrize('decomposed_cls', [HasKraus, HasMixture, HasUnitary])
+def test_has_superoperator_when_decomposed(decomposed_cls):
+    op = HasKrausWhenDecomposed(decomposed_cls).on(cirq.NamedQubit('test'))
+    assert cirq.has_superoperator(op)
+    assert not cirq.has_superoperator(op, allow_decompose=False)
+
+
+def apply_superoperator(superoperator: np.ndarray, rho: np.ndarray) -> np.ndarray:
+    """Computes output density matrix from superoperator and input density matrix."""
+    vectorized_input = np.reshape(rho, np.prod(rho.shape))
+    vectorized_output = superoperator @ vectorized_input
+    return np.reshape(vectorized_output, rho.shape)
+
+
+@pytest.mark.parametrize('op', (cirq.I, cirq.X, cirq.Y, cirq.Z, cirq.H, cirq.amplitude_damp(0.1)))
+@pytest.mark.parametrize(
+    'rho',
+    (
+        np.diag([1, 0]),
+        np.diag([0, 1]),
+        np.array([[0.5, 0.5], [0.5, 0.5]]),
+        np.array([[0.5, -0.5], [-0.5, 0.5]]),
+        np.array([[0.5, -0.5j], [0.5j, 0.5]]),
+    ),
+)
+def test_compare_superoperator_action_to_simulation_one_qubit(op, rho):
+    superoperator = cirq.superoperator(op)
+    actual_state = apply_superoperator(superoperator, rho)
+
+    sim = cirq.DensityMatrixSimulator()
+    q = cirq.NamedQubit('q')
+    circuit = cirq.Circuit(op.on(q))
+    initial_state = np.array(rho, dtype=np.complex64)
+    expected_state = sim.simulate(circuit, initial_state=initial_state).final_density_matrix
+
+    assert np.allclose(actual_state, expected_state)
+
+
+@pytest.mark.parametrize('op', (cirq.CNOT, cirq.ISWAP, cirq.depolarize(0.2, n_qubits=2)))
+@pytest.mark.parametrize(
+    'rho',
+    (
+        np.diag([1, 0, 0, 0]),
+        np.diag([0, 1, 0, 0]),
+        np.diag([0, 0, 1, 0]),
+        np.diag([0, 0, 0, 1]),
+        np.array([[0.5, 0.5, 0, 0], [0.5, 0.5, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]),
+        np.array([[0.5, -0.5j, 0, 0], [0.5j, 0.5, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]),
+        np.array([[0.5, 0, 0.5j, 0], [0, 0, 0, 0], [-0.5j, 0, 0.5, 0], [0, 0, 0, 0]]),
+    ),
+)
+def test_compare_superoperator_action_to_simulation_two_qubits(op, rho):
+    superoperator = cirq.superoperator(op)
+    actual_state = apply_superoperator(superoperator, rho)
+
+    sim = cirq.DensityMatrixSimulator()
+    q0, q1 = cirq.LineQubit.range(2)
+    circuit = cirq.Circuit(op.on(q0, q1))
+    initial_state = np.array(rho, dtype=np.complex64)
+    expected_state = sim.simulate(circuit, initial_state=initial_state).final_density_matrix
+
+    assert np.allclose(actual_state, expected_state)