diff --git a/doc/conf.py b/doc/conf.py
index 720f1872c..7ea7e7909 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -133,6 +133,8 @@
     "n_events",
     "n_cons",
     "max_n_chans",
+    "n_seeds",
+    "n_targets",
     "n_unique_seeds",
     "n_unique_targets",
     "variable",
diff --git a/mne_connectivity/base.py b/mne_connectivity/base.py
index a178c6f8b..aee5ff758 100644
--- a/mne_connectivity/base.py
+++ b/mne_connectivity/base.py
@@ -24,17 +24,16 @@
 class SpectralMixin:
     """Mixin class for spectral connectivities.
 
-    Note: In mne-connectivity, we associate the word
-    "spectral" with time-frequency. Reference to
-    eigenvalue structure is not captured in this mixin.
+    Note: In mne-connectivity, we associate the word "spectral" with time-frequency.
+    Reference to eigenvalue structure is not captured in this mixin.
     """
 
     @property
     def freqs(self):
         """The frequency points of the connectivity data.
 
-        If these are computed over a frequency band, it will
-        be the median frequency of the frequency band.
+        If these are computed over a frequency band, it will be the median frequency of
+        the frequency band.
         """
         return self.xarray.coords.get("freqs").values.tolist()
 
@@ -83,12 +82,12 @@ def append(self, epoch_conn):
         Parameters
         ----------
         epoch_conn : instance of Connectivity
-            The Epoched Connectivity class to append.
+            The epoched Connectivity class to append.
 
         Returns
         -------
         self : instance of Connectivity
-            The altered Epoched Connectivity class.
+            The altered epoched Connectivity class.
         """
         if not isinstance(self, type(epoch_conn)):
             raise ValueError(
@@ -137,18 +136,18 @@ def combine(self, combine="mean"):
 
         Parameters
         ----------
-        combine : 'mean' | 'median' | callable
-            How to combine correlation estimates across epochs.
-            Default is 'mean'. If callable, it must accept one
-            positional input. For example::
+        combine : ``'mean'`` | ``'median'`` | callable
+            How to combine correlation estimates across epochs. Default is ``'mean'``.
+            If callable, it must accept one positional input. For example::
 
                 combine = lambda data: np.median(data, axis=0)
 
         Returns
         -------
         conn : instance of Connectivity
-            The combined connectivity data structure.
-        """
+            The combined connectivity data structure. Instance type reflects that of the
+            input instance, without the epoch dimension.
+        """  # noqa: E501
         from .io import _xarray_to_conn
 
         if not self.is_epoched:
@@ -216,19 +215,17 @@ def predict(self, data):
 
         Parameters
         ----------
-        data : array
-            Epoched or continuous data set. Has shape
-            (n_epochs, n_signals, n_times) or (n_signals, n_times).
+        data : array, shape ([n_epochs,] n_signals, n_times)
+            Epoched or continuous data set.
 
         Returns
         -------
-        predicted : array
-            Data as predicted by the VAR model of
-            shape same as ``data``.
+        predicted : array, shape ([n_epochs,] n_signals, n_times)
+            Data as predicted by the VAR model of shape same as ``data``.
 
         Notes
         -----
-        Residuals are obtained by r = x - var.predict(x).
+        Residuals are obtained by ``r = x - var.predict(x)``.
 
         To compute residual covariances::
 
@@ -236,8 +233,7 @@ def predict(self, data):
             # row are observations, columns are variables
             t = residuals.shape[0]
             sampled_residuals = np.concatenate(
-                np.split(residuals[:, :, lags:], t, 0),
-                axis=2
+                np.split(residuals[:, :, lags:], t, 0), axis=2
             ).squeeze(0)
             rescov = np.cov(sampled_residuals)
         """
@@ -299,10 +295,9 @@ def simulate(self, n_samples, noise_func=None, random_state=None):
         ----------
         n_samples : int
             Number of samples to generate.
-        noise_func : func, optional
-            This function is used to create the generating noise process. If
-            set to None, Gaussian white noise with zero mean and unit variance
-            is used.
+        noise_func : callable | None
+            This function is used to create the generating noise process. If ``None``,
+            Gaussian white noise with zero mean and unit variance is used.
         %(random_state)s
 
         Returns
@@ -350,14 +345,13 @@ def noise_func():
 class BaseConnectivity(DynamicMixin, EpochMixin):
     """Base class for connectivity data.
 
-    This class should not be instantiated directly, but should be used
-    to do type-checking. All connectivity classes will be returned from
-    corresponding connectivity computing functions.
+    This class should not be instantiated directly, but should be used to do
+    type-checking. All connectivity classes will be returned from corresponding
+    connectivity computing functions.
 
-    Connectivity data is anything that represents "connections"
-    between nodes as a (N, N) array. It can be symmetric, or
-    asymmetric (if it is symmetric, storage optimization will
-    occur).
+    Connectivity data is anything that represents "connections" between nodes as a
+    ``(N, N)`` array. It can be symmetric, or asymmetric (if it is symmetric, storage
+    optimization will occur).
 
     Parameters
     ----------
@@ -369,33 +363,30 @@ class BaseConnectivity(DynamicMixin, EpochMixin):
     %(events)s
     %(event_id)s
     metadata : instance of pandas.DataFrame | None
-        The metadata data frame that would come from the :class:`mne.Epochs`
-        class. See :class:`mne.Epochs` docstring for details.
+        The metadata data frame that would come from the :class:`mne.Epochs` class. See
+        :class:`mne.Epochs` docstring for details.
     %(connectivity_kwargs)s
 
     Notes
     -----
-    Connectivity data can be generally represented as a square matrix
-    with values intending the connectivity function value between two
-    nodes. We optimize storage of symmetric connectivity data
-    and allow support for computing connectivity data on a subset of nodes.
-    We store connectivity data as a raveled ``(n_estimated_nodes, ...)``
-    where ``n_estimated_nodes`` can be ``n_nodes_in * n_nodes_out`` if a
-    full connectivity structure is computed, or a subset of the nodes
-    (equal to the length of the indices passed in).
-
-    Since we store connectivity data as a raveled array, one can
-    easily optimize the storage of "symmetric" connectivity data.
-    One can use numpy to convert a full all-to-all connectivity
-    into an upper triangular portion, and set ``indices='symmetric'``.
-    This would reduce the RAM needed in half.
-
-    The underlying data structure is an ``xarray.DataArray``,
-    with a similar API to ``xarray``. We provide support for storing
-    connectivity data in a subset of nodes. Thus the underlying
-    data structure instead of a ``(n_nodes_in, n_nodes_out)`` 2D array
-    would be a ``(n_nodes_in * n_nodes_out,)`` raveled 1D array. This
-    allows us to optimize storage also for symmetric connectivity.
+    Connectivity data can be generally represented as a square matrix with values
+    intending the connectivity function value between two nodes. We optimize storage of
+    symmetric connectivity data and allow support for computing connectivity data on a
+    subset of nodes. We store connectivity data as a raveled ``(n_estimated_nodes,
+    ...)`` where ``n_estimated_nodes`` can be ``n_nodes_in * n_nodes_out`` if a full
+    connectivity structure is computed, or a subset of the nodes (equal to the length of
+    the indices passed in).
+
+    Since we store connectivity data as a raveled array, one can easily optimize the
+    storage of "symmetric" connectivity data. One can use numpy to convert a full
+    all-to-all connectivity into an upper triangular portion, and set
+    ``indices='symmetric'``. This would reduce the RAM needed in half.
+
+    The underlying data structure is an :class:`xarray.DataArray`, with a similar API to
+    ``xarray``. We provide support for storing connectivity data in a subset of nodes.
+    Thus the underlying data structure instead of a ``(n_nodes_in, n_nodes_out)`` 2D
+    array would be a ``(n_nodes_in * n_nodes_out,)`` raveled 1D array. This allows us to
+    optimize storage also for symmetric connectivity.
     """
 
     # whether or not the connectivity occurs over epochs
@@ -628,9 +619,8 @@ def shape(self):
     def n_nodes(self):
         """The number of nodes in the original dataset.
 
-        Even if ``indices`` defines a subset of nodes that
-        were computed, this should be the total number of
-        nodes in the original dataset.
+        Even if ``indices`` defines a subset of nodes that were computed, this should be
+        the total number of nodes in the original dataset.
         """
         return self.attrs["n_nodes"]
 
@@ -645,11 +635,10 @@ def indices(self):
 
         Returns
         -------
-        indices : str | tuple of lists
-            Either 'all' for all-to-all connectivity,
-            'symmetric' for symmetric all-to-all connectivity,
-            or a tuple of lists representing the node-to-nodes
-            that connectivity was computed.
+        indices : ``'all'`` | ``'symmetric'`` | tuple of list
+            Either ``'all'`` for all-to-all connectivity, ``'symmetric'`` for symmetric
+            connectivity, or a tuple of lists representing the node-to-nodes that
+            connectivity was computed for.
         """
         return self.attrs["indices"]
 
@@ -667,9 +656,8 @@ def xarray(self):
     def n_epochs_used(self):
         """Number of epochs used in computation of connectivity.
 
-        Can be 'None', if there was no epochs used. This is
-        equivalent to the number of epochs, if there is no
-        combining of epochs.
+        Can be ``None``, if there was no epochs used. This is equivalent to the number
+        of epochs, if there is no combining of epochs.
         """
         return self.attrs.get("n_epochs_used")
 
@@ -689,18 +677,21 @@ def get_data(self, output="compact"):
 
         Parameters
         ----------
-        output : str, optional
-            How to format the output, by default 'raveled', which
-            will represent each connectivity matrix as a
-            ``(n_nodes_in * n_nodes_out,)`` list. If 'dense', then
-            will return each connectivity matrix as a 2D array. If 'compact'
-            (default) then will return 'raveled' if ``indices`` were defined as
-            a list of tuples, or ``dense`` if indices is 'all'. Multivariate
-            connectivity data cannot be returned in a dense form.
+        output : ``'compact'`` | ``'raveled'`` | ``'dense'``
+            How to format the output:
+
+            - ``'raveled'`` will represent each connectivity matrix as a
+              ``(..., n_nodes_in * n_nodes_out, ...)`` array
+            - ``'dense'`` will return each connectivity matrix as a ``(..., n_nodes_in,
+              n_nodes_out, ...)`` array
+            - ``'compact'`` (default) will return ``'raveled'`` if ``indices`` were
+              defined as a tuple of arrays, or ``'dense'`` if ``indices='all'``
+
+            Multivariate connectivity data cannot be returned in a dense form.
 
         Returns
         -------
-        data : np.ndarray
+        data : array
             The output connectivity data.
         """
         _check_option("output", output, ["raveled", "dense", "compact"])
@@ -830,13 +821,13 @@ def save(self, fname):
         """Save connectivity data to disk.
 
         Can later be loaded using the function
-        :func:`mne_connectivity.read_connectivity`.
+        :func:`~mne_connectivity.read_connectivity`.
 
         Parameters
         ----------
         fname : str | pathlib.Path
-            The filepath to save the data. Data is saved
-            as netCDF files (``.nc`` extension).
+            The filepath to save the data. Data is saved as netCDF files (``.nc``
+            extension).
         """
         method = self.method
         indices = self.indices
@@ -882,9 +873,10 @@ class SpectralConnectivity(BaseConnectivity, SpectralMixin):
     """Spectral connectivity class.
 
     This class stores connectivity data that varies over frequencies. The underlying
-    data is an array of shape (n_connections, [n_components], n_freqs), or (n_nodes,
-    n_nodes, [n_components], n_freqs). ``n_components`` is an optional dimension for
-    multivariate methods where each connection has multiple components of connectivity.
+    data is an array of shape ``(n_connections, [n_components,] n_freqs)``, or
+    ``(n_nodes, n_nodes, [n_components,] n_freqs)``. ``n_components`` is an optional
+    dimension for multivariate methods where each connection has multiple components of
+    connectivity.
 
     Parameters
     ----------
@@ -936,12 +928,12 @@ def __init__(
 class TemporalConnectivity(BaseConnectivity, TimeMixin):
     """Temporal connectivity class.
 
-    This is an array of shape (n_connections, [n_components], n_times), or (n_nodes,
-    n_nodes, [n_components], n_times). This describes how connectivity varies over
-    time. It describes sample-by-sample time-varying connectivity (usually on the order
-    of milliseconds). Here time (t=0) is the same for all connectivity measures.
-    ``n_components`` is an optional dimension for multivariate methods where each
-    connection has multiple components of connectivity.
+    This is an array of shape ``(n_connections, [n_components,] n_times)``, or
+    ``(n_nodes, n_nodes, [n_components,] n_times)``. This describes how connectivity
+    varies over time. It describes sample-by-sample time-varying connectivity (usually
+    on the order of milliseconds). Here time (t=0) is the same for all connectivity
+    measures. ``n_components`` is an optional dimension for multivariate methods where
+    each connection has multiple components of connectivity.
 
     Parameters
     ----------
@@ -956,11 +948,11 @@ class TemporalConnectivity(BaseConnectivity, TimeMixin):
 
     Notes
     -----
-    `mne_connectivity.EpochConnectivity` is a similar connectivity class to this one.
-    However, that describes one connectivity snapshot for each epoch. These epochs might
-    be chunks of time that have different meaning for time ``t=0``. Epochs can mean
-    separate trials, where the beginning of the trial implies t=0. These Epochs may also
-    be discontiguous.
+    :class:`mne_connectivity.EpochConnectivity` is a similar connectivity class to this
+    one. However, that describes one connectivity snapshot for each epoch. These epochs
+    might be chunks of time that have different meaning for time ``t=0``. Epochs can
+    mean separate trials, where the beginning of the trial implies t=0. These epochs may
+    also be discontiguous.
     """
 
     expected_n_dim = 2
@@ -994,12 +986,12 @@ class SpectroTemporalConnectivity(BaseConnectivity, SpectralMixin, TimeMixin):
 
     This class stores connectivity data that varies over both frequency and time. The
     temporal part describes sample-by-sample time-varying connectivity (usually on the
-    order of milliseconds). Note the difference relative to Epochs.
+    order of milliseconds). Note the difference relative to epochs.
 
-    The underlying data is an array of shape (n_connections, [n_components], n_freqs,
-    n_times), or (n_nodes, n_nodes, [n_components], n_freqs, n_times). ``n_components``
-    is an optional dimension for multivariate methods where each connection has multiple
-    components of connectivity.
+    The underlying data is an array of shape ``(n_connections, [n_components,] n_freqs,
+    n_times)``, or ``(n_nodes, n_nodes, [n_components,] n_freqs, n_times)``.
+    ``n_components`` is an optional dimension for multivariate methods where each
+    connection has multiple components of connectivity.
 
     Parameters
     ----------
@@ -1049,10 +1041,10 @@ def __init__(
 
 @fill_doc
 class EpochSpectralConnectivity(SpectralConnectivity):
-    """Spectral connectivity class over Epochs.
+    """Spectral connectivity class over epochs.
 
-    This is an array of shape (n_epochs, n_connections, [n_components], n_freqs), or
-    (n_epochs, n_nodes, n_nodes, [n_components], n_freqs). This describes how
+    This is an array of shape ``(n_epochs, n_connections, [n_components,] n_freqs)``, or
+    ``(n_epochs, n_nodes, n_nodes, [n_components,] n_freqs)``. This describes how
     connectivity varies over frequencies for different epochs. ``n_components`` is an
     optional dimension for multivariate methods where each connection has multiple
     components of connectivity.
@@ -1101,10 +1093,10 @@ def __init__(
 
 @fill_doc
 class EpochTemporalConnectivity(TemporalConnectivity):
-    """Temporal connectivity class over Epochs.
+    """Temporal connectivity class over epochs.
 
-    This is an array of shape (n_epochs, n_connections, [n_components], n_times), or
-    (n_epochs, n_nodes, n_nodes, [n_components], n_times). This describes how
+    This is an array of shape ``(n_epochs, n_connections, [n_components,] n_times)``, or
+    ``(n_epochs, n_nodes, n_nodes, [n_components,] n_times)``. This describes how
     connectivity varies over time for different epochs. ``n_components`` is an optional
     dimension for multivariate methods where each connection has multiple components of
     connectivity.
@@ -1144,13 +1136,13 @@ def __init__(
 
 @fill_doc
 class EpochSpectroTemporalConnectivity(SpectroTemporalConnectivity):
-    """Spectrotemporal connectivity class over Epochs.
+    """Spectrotemporal connectivity class over epochs.
 
-    This is an array of shape (n_epochs, n_connections, [n_components], n_freqs,
-    n_times), or (n_epochs, n_nodes, n_nodes, [n_components], n_freqs, n_times). This
-    describes how connectivity varies over frequencies and time for different epochs.
-    ``n_components`` is an optional dimension for multivariate methods where each
-    connection has multiple components of connectivity.
+    This is an array of shape ``(n_epochs, n_connections, [n_components,] n_freqs,
+    n_times)``, or ``(n_epochs, n_nodes, n_nodes, [n_components,] n_freqs, n_times)``.
+    This describes how connectivity varies over frequencies and time for different
+    epochs. ``n_components`` is an optional dimension for multivariate methods where
+    each connection has multiple components of connectivity.
 
     Parameters
     ----------
@@ -1197,9 +1189,9 @@ def __init__(
 class Connectivity(BaseConnectivity):
     """Connectivity class without frequency or time component.
 
-    This is an array of shape (n_connections, [n_components]), or (n_nodes, n_nodes,
-    [n_components]). This describes a connectivity matrix/graph that does not vary
-    over time, frequency, or epochs. ``n_components`` is an optional dimension for
+    This is an array of shape ``(n_connections[, n_components])``, or ``(n_nodes,
+    n_nodes[, n_components])``. This describes a connectivity matrix/graph that does not
+    vary over time, frequency, or epochs. ``n_components`` is an optional dimension for
     multivariate methods where each connection has multiple components of connectivity.
 
     Parameters
@@ -1242,10 +1234,10 @@ def __init__(
 class EpochConnectivity(BaseConnectivity):
     """Epoch connectivity class.
 
-    This is an array of shape (n_epochs, n_connections, [n_components]), or (n_epochs,
-    n_nodes, n_nodes, [n_components]). This describes how connectivity varies for
-    different epochs. ``n_components`` is an optional dimension for multivariate methods
-    where each connection has multiple components of connectivity.
+    This is an array of shape ``(n_epochs, n_connections[, n_components])``, or
+    ``(n_epochs, n_nodes, n_nodes[, n_components])``. This describes how connectivity
+    varies for different epochs. ``n_components`` is an optional dimension for
+    multivariate methods where each connection has multiple components of connectivity.
 
     Parameters
     ----------
diff --git a/mne_connectivity/datasets/frequency.py b/mne_connectivity/datasets/frequency.py
index 8fd293089..fbad482e1 100644
--- a/mne_connectivity/datasets/frequency.py
+++ b/mne_connectivity/datasets/frequency.py
@@ -14,11 +14,11 @@ def make_signals_in_freq_bands(
     freq_band,
     n_epochs=10,
     n_times=200,
-    sfreq=100,
-    trans_bandwidth=1,
+    sfreq=100.0,
+    trans_bandwidth=1.0,
     snr=0.7,
     connection_delay=5,
-    tmin=0,
+    tmin=0.0,
     ch_names=None,
     ch_types="eeg",
     rng_seed=None,
@@ -31,7 +31,7 @@ def make_signals_in_freq_bands(
         Number of seed channels to simulate.
     n_targets : int
         Number of target channels to simulate.
-    freq_band : tuple of int or float
+    freq_band : tuple of float
         Frequency band where the connectivity should be simulated, where the first entry
         corresponds to the lower frequency, and the second entry to the higher
         frequency.
@@ -39,9 +39,9 @@ def make_signals_in_freq_bands(
         Number of epochs in the simulated data.
     n_times : int (default 200)
         Number of timepoints each epoch of the simulated data.
-    sfreq : int | float (default 100)
+    sfreq : float (default 100.0)
         Sampling frequency of the simulated data, in Hz.
-    trans_bandwidth : int | float (default 1)
+    trans_bandwidth : float (default 1.0)
         Transition bandwidth of the filter to apply to isolate activity in
         ``freq_band``, in Hz. These are passed to the ``l_bandwidth`` and
         ``h_bandwidth`` keyword arguments in :func:`mne.filter.create_filter`.
@@ -51,23 +51,24 @@ def make_signals_in_freq_bands(
         Number of timepoints for the delay of connectivity between the seeds and
         targets. If > 0, the target data is a delayed form of the seed data. If < 0, the
         seed data is a delayed form of the target data.
-    tmin : int | float (default 0)
+    tmin : float (default 0.0)
         Earliest time of each epoch.
     ch_names : list of str | None (default None)
-        Names of the channels in the simulated data. If `None`, the channels are named
+        Names of the channels in the simulated data. If ``None``, the channels are named
         according to their index and the frequency band of interaction. If specified,
         must be a list of ``n_seeds + n_targets`` channel names.
-    ch_types : str | list of str (default "eeg")
-        Types of the channels in the simulated data. If specified as a list, must be a
+    ch_types : str | list of str (default ``'eeg'``)
+        Types of the channels in the simulated data. Must be a recognised data channel
+        type (see :term:`mne:data channels`). If specified as a list, must be a
         list of ``n_seeds + n_targets`` channel names.
     rng_seed : int | None (default None)
-        Seed to use for the random number generator. If `None`, no seed is specified.
+        Seed to use for the random number generator. If ``None``, no seed is specified.
 
     Returns
     -------
-    epochs : mne.EpochsArray of shape (n_epochs, ``n_seeds + n_targets``, n_times)
-        The simulated data stored in an `mne.EpochsArray` object. The channels are
-        arranged according to seeds, then targets.
+    epochs : ~mne.EpochsArray, shape (n_epochs, n_seeds + n_targets, n_times)
+        The simulated data stored in an :class:`mne.EpochsArray` object. The channels
+        are arranged according to seeds, then targets.
 
     Notes
     -----
diff --git a/mne_connectivity/datasets/surrogate.py b/mne_connectivity/datasets/surrogate.py
index a482b5c60..6e0adbdeb 100644
--- a/mne_connectivity/datasets/surrogate.py
+++ b/mne_connectivity/datasets/surrogate.py
@@ -3,7 +3,7 @@
 # License: BSD (3-clause)
 
 import numpy as np
-from mne.time_frequency import EpochsSpectrum, EpochsSpectrumArray
+from mne.time_frequency import EpochsSpectrum
 from mne.utils import _validate_type
 
 
@@ -12,25 +12,25 @@ def make_surrogate_data(data, n_shuffles=1000, rng_seed=None, return_generator=T
 
     Parameters
     ----------
-    data : ~mne.time_frequency.EpochsSpectrum | ~mne.time_frequency.EpochsSpectrumArray
+    data : ~mne.time_frequency.EpochsSpectrum
         The Fourier coefficients to create the null hypothesis surrogate data for. Can
         be generated from :meth:`mne.Epochs.compute_psd` with ``output='complex'``
         (requires ``mne >= 1.8``).
     n_shuffles : int (default 1000)
         The number of surrogate datasets to create.
     rng_seed : int | None (default None)
-        The seed to use for the random number generator. If `None`, no seed is
+        The seed to use for the random number generator. If ``None``, no seed is
         specified.
     return_generator : bool (default True)
-        Whether or not to return the surrogate data as a :term:`generator` object
-        instead of a :class:`list`. This allows iterating over the surrogates without
-        having to keep them all in memory.
+        Whether or not to return the surrogate data as a generator object instead of a
+        list. This allows iterating over the surrogates without having to keep them all
+        in memory.
 
     Returns
     -------
     surrogate_data : list of ~mne.time_frequency.EpochsSpectrum
         The surrogate data for the null hypothesis with ``n_shuffles`` entries. Returned
-        as a :term:`generator` if ``return_generator=True``.
+        as a generator if ``return_generator=True``.
 
     Notes
     -----
@@ -78,7 +78,7 @@ def make_surrogate_data(data, n_shuffles=1000, rng_seed=None, return_generator=T
     # Validate inputs
     _validate_type(
         data,
-        (EpochsSpectrum, EpochsSpectrumArray),
+        (EpochsSpectrum),
         "data",
         "mne.time_frequency.EpochsSpectrum or mne.time_frequency.EpochsSpectrumArray",
     )
diff --git a/mne_connectivity/decoding/decomposition.py b/mne_connectivity/decoding/decomposition.py
index b119cceeb..99f3bbf9d 100644
--- a/mne_connectivity/decoding/decomposition.py
+++ b/mne_connectivity/decoding/decomposition.py
@@ -62,17 +62,17 @@ class CoherencyDecomposition(BaseEstimator, TransformerMixin):
     Once fit, the filters can be used to transform data into the underlying connectivity
     components. Connectivity can be computed on this transformed data using the
     bivariate coherency-based methods of the
-    `mne_connectivity.spectral_connectivity_epochs` and
-    `mne_connectivity.spectral_connectivity_time` functions. These bivariate methods
-    are:
+    :func:`~mne_connectivity.spectral_connectivity_epochs` and
+    :func:`~mne_connectivity.spectral_connectivity_time` functions. These bivariate
+    methods are:
 
     * ``"cohy"`` and ``"coh"`` for CaCoh :footcite:`VidaurreEtAl2019`
     * ``"imcoh"`` for MIC :footcite:`EwaldEtAl2012`
 
     The approach taken here is to optimise the connectivity in a given frequency band.
     Frequency bin-wise optimisation is offered in the multivariate coherency-based
-    methods of the `mne_connectivity.spectral_connectivity_epochs` and
-    `mne_connectivity.spectral_connectivity_time` functions.
+    methods of the :func:`~mne_connectivity.spectral_connectivity_epochs` and
+    :func:`~mne_connectivity.spectral_connectivity_time` functions.
 
     References
     ----------
@@ -131,7 +131,7 @@ def __init__(
         mt_adaptive=False,
         mt_low_bias=True,
         cwt_freqs=None,
-        cwt_n_cycles=7,
+        cwt_n_cycles=7.0,
         n_components=None,
         rank=None,
         n_jobs=1,
@@ -277,7 +277,7 @@ def fit(self, X, y=None):
 
         Parameters
         ----------
-        X : array, shape=(n_epochs, n_signals, n_times)
+        X : array, shape (n_epochs, n_signals, n_times)
             The input data which the connectivity decomposition filters should be fit
             to.
         y : None
@@ -397,12 +397,12 @@ def transform(self, X):
 
         Parameters
         ----------
-        X : array, shape=((n_epochs, ) n_signals, n_times)
+        X : array, shape ([n_epochs,] n_signals, n_times)
             The data to be transformed by the connectivity decomposition filters.
 
         Returns
         -------
-        X_transformed : array, shape=((n_epochs, ) n_components*2, n_times)
+        X_transformed : array, shape ([n_epochs,] n_components*2, n_times)
             The transformed data. The first ``n_components`` channels are the
             transformed seeds, and the last ``n_components`` channels are the
             transformed targets.
@@ -424,7 +424,7 @@ def fit_transform(self, X, y=None, **fit_params):
 
         Parameters
         ----------
-        X : array, shape=(n_epochs, n_signals, n_times)
+        X : array, shape (n_epochs, n_signals, n_times)
             The input data which the connectivity decomposition filters should be fit to
             and subsequently transformed.
         y : None
@@ -435,7 +435,7 @@ def fit_transform(self, X, y=None, **fit_params):
 
         Returns
         -------
-        X_transformed : array, shape=(n_epochs, n_components*2, n_times)
+        X_transformed : array, shape (n_epochs, n_components*2, n_times)
             The transformed data. The first ``n_components`` channels are the
             transformed seeds, and the last ``n_components`` channels are the
             transformed targets.
@@ -449,11 +449,11 @@ def get_transformed_indices(self):
         Returns
         -------
         indices_transformed : tuple of array
-            Indices of seeds and targets in the transformed data with the form (seeds,
-            targets) to be used when passing the data to
-            `~mne_connectivity.spectral_connectivity_epochs` and
-            `~mne_connectivity.spectral_connectivity_time`. Entries of the indices are
-            arranged such that connectivity would be computed between the first seed
+            Indices of seeds and targets in the transformed data with the form ``(seeds,
+            targets)`` to be used when passing the data to
+            :func:`~mne_connectivity.spectral_connectivity_epochs` and
+            :func:`~mne_connectivity.spectral_connectivity_time`. Entries of the indices
+            are arranged such that connectivity would be computed between the first seed
             component and first target component, second seed component and second
             target component, etc...
         """
@@ -495,8 +495,8 @@ def plot_patterns(
     ):
         """Plot topographic patterns of components.
 
-        The patterns explain how the measured data was generated from the
-        neural sources (a.k.a. the forward model) :footcite:`HaufeEtAl2014`.
+        The patterns explain how the measured data was generated from the neural sources
+        (a.k.a. the forward model) :footcite:`HaufeEtAl2014`.
 
         Seed and target patterns are plotted separately.
 
diff --git a/mne_connectivity/effective.py b/mne_connectivity/effective.py
index ab598418c..a2bcb10f4 100644
--- a/mne_connectivity/effective.py
+++ b/mne_connectivity/effective.py
@@ -35,68 +35,66 @@ def phase_slope_index(
 ):
     """Compute the Phase Slope Index (PSI) connectivity measure.
 
-    The PSI is an effective connectivity measure, i.e., a measure which can
-    give an indication of the direction of the information flow (causality).
-    For two time series, and one computes the PSI between the first and the
-    second time series as follows
+    The PSI is an effective connectivity measure, i.e., a measure which can give an
+    indication of the direction of the information flow (causality). For two time
+    series, and one computes the PSI between the first and the second time series as
+    follows::
 
-    indices = (np.array([0]), np.array([1]))
-    psi = phase_slope_index(data, indices=indices, ...)
+        indices = (np.array([0]), np.array([1]))
+        psi = phase_slope_index(data, indices=indices, ...)
 
-    A positive value means that time series 0 is ahead of time series 1 and
-    a negative value means the opposite.
+    A positive value means that time series 0 is ahead of time series 1 and a negative
+    value means the opposite.
 
     The PSI is computed from the coherency (see :func:`spectral_connectivity_epochs`),
     details can be found in :footcite:`NolteEtAl2008`.
 
     Parameters
     ----------
-    data : array-like, shape=(n_epochs, n_signals, n_times)
-        Can also be a list/generator of array, shape =(n_signals, n_times);
-        list/generator of SourceEstimate; or Epochs.
-        The data from which to compute connectivity. Note that it is also
-        possible to combine multiple signals by providing a list of tuples,
-        e.g., data = [(arr_0, stc_0), (arr_1, stc_1), (arr_2, stc_2)],
-        corresponds to 3 epochs, and arr_* could be an array with the same
-        number of time points as stc_*.
+    data : array_like, shape (n_epochs, n_signals, n_times) | ~mne.Epochs | generator
+        Can also be a list/generator of arrays, shape ``(n_signals, n_times)``;
+        list/generator of :class:`mne.SourceEstimate`; or :class:`mne.Epochs`. The
+        data from which to compute connectivity. Note that it is also possible to
+        combine multiple signals by providing a list of tuples, e.g., ``data = [(arr_0,
+        stc_0), (arr_1, stc_1), (arr_2, stc_2)]``, corresponds to 3 epochs, and
+        ``arr_*`` could be an array with the same number of time points as ``stc_*``.
     %(names)s
-    indices : tuple of array | None
-        Two arrays with indices of connections for which to compute
-        connectivity. If None, all connections are computed.
+    indices : tuple of array_like | None
+        Two array-likes with indices of connections for which to compute connectivity.
+        If ``None``, all connections are computed. See Notes of
+        :func:`~mne_connectivity.spectral_connectivity_epochs` for details.
     sfreq : float
         The sampling frequency.
-    mode : str
-        Spectrum estimation mode can be either: 'multitaper', 'fourier', or
-        'cwt_morlet'.
+    mode : ``'multitaper'`` | ``'fourier'`` | ``'cwt_morlet'``
+        Spectrum estimation mode.
     fmin : float | tuple of float
-        The lower frequency of interest. Multiple bands are defined using
-        a tuple, e.g., (8., 20.) for two bands with 8Hz and 20Hz lower freq.
-        If None the frequency corresponding to an epoch length of 5 cycles
-        is used.
+        The lower frequency of interest. Multiple bands are defined using a tuple, e.g.,
+        (8., 20.) for two bands with 8 Hz and 20 Hz lower freq. If ``None`` the
+        frequency corresponding to an epoch length of 5 cycles is used.
     fmax : float | tuple of float
-        The upper frequency of interest. Multiple bands are dedined using
-        a tuple, e.g. (13., 30.) for two band with 13Hz and 30Hz upper freq.
+        The upper frequency of interest. Multiple bands are defined using a tuple, e.g.,
+        (13., 30.) for two bands with 13 Hz and 30 Hz upper freq.
     tmin : float | None
         Time to start connectivity estimation.
     tmax : float | None
         Time to end connectivity estimation.
     mt_bandwidth : float | None
-        The bandwidth of the multitaper windowing function in Hz.
-        Only used in 'multitaper' mode.
+        The bandwidth of the multitaper windowing function in Hz. Only used in
+        ``'multitaper'`` mode.
     mt_adaptive : bool
-        Use adaptive weights to combine the tapered spectra into PSD.
-        Only used in 'multitaper' mode.
+        Use adaptive weights to combine the tapered spectra into PSD. Only used in
+        ``'multitaper'`` mode.
     mt_low_bias : bool
         Only use tapers with more than 90 percent spectral concentration within
-        bandwidth. Only used in 'multitaper' mode.
-    cwt_freqs : array
-        Array of frequencies of interest. Only used in 'cwt_morlet' mode.
-    cwt_n_cycles : float | array of float
+        bandwidth. Only used in ``'multitaper'`` mode.
+    cwt_freqs : array_like
+        Array-like of frequencies of interest. Only used in ``'cwt_morlet'`` mode.
+    cwt_n_cycles : float | array_like
         Number of cycles. Fixed number or one per frequency. Only used in
-        'cwt_morlet' mode.
+        ``'cwt_morlet'`` mode.
     block_size : int
-        How many connections to compute at once (higher numbers are faster
-        but require more memory).
+        How many connections to compute at once (higher numbers are faster but require
+        more memory).
     n_jobs : int
         How many epochs to process in parallel.
     %(verbose)s
@@ -105,11 +103,12 @@ def phase_slope_index(
     -------
     conn : instance of SpectralConnectivity or SpectroTemporalConnectivity
         Computed connectivity measure(s). Either a :class:`SpectralConnectivity`, or
-        :class:`SpectroTemporalConnectivity` container. The shape of each array is
-        either (n_signals ** 2, n_bands) mode: 'multitaper' or 'fourier' (n_signals **
-        2, n_bands, n_times) mode: 'cwt_morlet' when "indices" is None, or (n_con,
-        n_bands) mode: 'multitaper' or 'fourier' (n_con, n_bands, n_times) mode:
-        'cwt_morlet' when "indices" is specified and "n_con = len(indices[0])".
+        :class:`SpectroTemporalConnectivity` container. The shape of each array is:
+
+        - ``(n_cons, n_bands)`` for ``'multitaper'`` or ``'fourier'`` modes
+        - ``(n_cons, n_bands, n_times)`` for ``'cwt_morlet'`` mode
+        - ``n_cons = n_signals ** 2`` when ``indices=None``
+        - ``n_cons = len(indices[0])`` when ``indices`` is supplied
 
     See Also
     --------
diff --git a/mne_connectivity/envelope.py b/mne_connectivity/envelope.py
index 43c92f26b..61336b516 100644
--- a/mne_connectivity/envelope.py
+++ b/mne_connectivity/envelope.py
@@ -24,38 +24,34 @@ def envelope_correlation(
 
     Parameters
     ----------
-    data : array-like, shape=(n_epochs, n_signals, n_times) | Epochs | generator
-        The data from which to compute connectivity.
-        The array-like object can also be a list/generator of array,
-        each with shape (n_signals, n_times), or a :class:`~mne.SourceEstimate`
-        object (and ``stc.data`` will be used). If it's float data,
-        the Hilbert transform will be applied; if it's complex data,
-        it's assumed the Hilbert has already been applied.
-    names : list | array-like | None
-        A list of names associated with the signals in ``data``.
-        If None, will be a list of indices of the number of nodes.
-    orthogonalize : 'pairwise' | False
-        Whether to orthogonalize with the pairwise method or not.
-        Defaults to 'pairwise'. Note that when False,
-        the correlation matrix will not be returned with
-        absolute values.
+    data : array_like, shape (n_epochs, n_signals, n_times) | Epochs | generator
+        The data from which to compute connectivity. The array-like object can also be a
+        list/generator of arrays, each with shape ``(n_signals, n_times)``, or a
+        :class:`mne.SourceEstimate` object (and ``stc.data`` will be used). If float
+        data, the Hilbert transform will be applied; if complex data, it is assumed the
+        Hilbert has already been applied.
+    names : array_like | None
+        A list of names associated with the signals in ``data``. If ``None``, will be a
+        list of indices of the number of nodes.
+    orthogonalize : ``'pairwise'`` | False
+        Whether to orthogonalize with the pairwise method or not. Defaults to
+        ``'pairwise'``. Note that when ``False``, the correlation matrix will not be
+        returned with absolute values.
     log : bool
-        If True (default False), square and take the log before orthogonalizing
+        If ``True`` (default ``False``), square and take the log before orthogonalizing
         envelopes or computing correlations.
     absolute : bool
-        If True (default), then take the absolute value of correlation
-        coefficients before making each epoch's correlation matrix
-        symmetric (and thus before combining matrices across epochs).
-        Only used when ``orthogonalize='pairwise'``.
+        If ``True`` (default), then take the absolute value of correlation coefficients
+        before making each epoch's correlation matrix symmetric (and thus before
+        combining matrices across epochs). Only used when ``orthogonalize='pairwise'``.
     %(verbose)s
 
     Returns
     -------
     corr : instance of EpochTemporalConnectivity
-        The pairwise orthogonal envelope correlations.
-        This matrix is symmetric. The array
-        will have three dimensions, the first of which is ``n_epochs``.
-        The data shape would be ``(n_epochs, (n_nodes + 1) * n_nodes / 2)``.
+        The pairwise orthogonal envelope correlations. This matrix is symmetric. The
+        array will have three dimensions, the first of which is ``n_epochs``. The data
+        shape is ``(n_epochs, (n_nodes + 1) * n_nodes / 2)``.
 
     See Also
     --------
@@ -63,12 +59,12 @@ def envelope_correlation(
 
     Notes
     -----
-    This function computes the power envelope correlation between
-    orthogonalized signals :footcite:`HippEtAl2012,KhanEtAl2018`.
+    This function computes the power envelope correlation between orthogonalized signals
+    :footcite:`HippEtAl2012,KhanEtAl2018`.
 
-    If you would like to combine Epochs after the fact using some function over the
-    Epochs axis, see the :meth:`~EpochTemporalConnectivity.combine` method of the
-    :class:`EpochTemporalConnectivity` class.
+    If you would like to combine epochs after the fact using some function over the
+    epochs axis, see the :meth:`~mne_connectivity.EpochTemporalConnectivity.combine`
+    method of the :class:`~mne_connectivity.EpochTemporalConnectivity` class.
 
     References
     ----------
@@ -225,14 +221,13 @@ def envelope_correlation(
 def symmetric_orth(data, *, n_iter=50, tol=1e-6, verbose=None):
     """Perform symmetric orthogonalization.
 
-    Uses the method from :footcite:`ColcloughEtAl2015` to jointly
-    orthogonalize the time series.
+    Uses the method from :footcite:`ColcloughEtAl2015` to jointly orthogonalize the time
+    series.
 
     Parameters
     ----------
-    data : ndarray, shape ([n_epochs, ]n_signals, n_times) or generator
-        The data to process. If a generator, it must return 2D arrays to
-        process.
+    data : array, shape ([n_epochs,] n_signals, n_times) | generator
+        The data to process. If a generator, it must return 2D arrays to process.
     n_iter : int
         The maximum number of iterations to perform.
     tol : float
@@ -241,9 +236,8 @@ def symmetric_orth(data, *, n_iter=50, tol=1e-6, verbose=None):
 
     Returns
     -------
-    data_orth : ndarray, shape (n_epochs, n_signals, n_times) | generator
-        The orthogonalized data. If ``data`` is a generator, a generator
-        is returned.
+    data_orth : array, shape ([n_epochs,] n_signals, n_times) | generator
+        The orthogonalized data. If ``data`` is a generator, a generator is returned.
 
     References
     ----------
diff --git a/mne_connectivity/spectral/epochs.py b/mne_connectivity/spectral/epochs.py
index 0bf902e3c..eb512adfd 100644
--- a/mne_connectivity/spectral/epochs.py
+++ b/mne_connectivity/spectral/epochs.py
@@ -171,7 +171,7 @@ def _prepare_connectivity(
         tmin_idx = None
         tmax_idx = None
         warn_times = False
-    else:  # data has a time dimension (timeseries or TFR object)
+    else:  # data has a time dimension (time series or TFR object)
         if spectrum_computed:  # is a TFR object
             if mode == "cwt_morlet":
                 first_epoch = (first_epoch[0][:, 0],)  # just take first freq
@@ -191,7 +191,7 @@ def _prepare_connectivity(
         )
 
     # Sort freqs
-    if not spectrum_computed:  # is an (ordinary) timeseries
+    if not spectrum_computed:  # is an (ordinary) time series
         # check that fmin corresponds to at least 5 cycles
         fmin = _check_freqs(sfreq=sfreq, fmin=fmin, n_times=n_times)
         # compute frequencies to analyze based on number of samples, sampling rate,
@@ -772,44 +772,44 @@ def spectral_connectivity_epochs(
 ):
     r"""Compute frequency- and time-frequency-domain connectivity measures.
 
-    The connectivity method(s) are specified using the "method" parameter.
-    All methods are based on estimates of the cross- and power spectral
-    densities (CSD/PSD) Sxy and Sxx, Syy.
+    The connectivity method(s) are specified using the "method" parameter. All methods
+    are based on estimates of the cross- and power spectral densities (CSD/PSD) Sxy and
+    Sxx, Syy.
 
     Parameters
     ----------
-    data : array-like, shape=(n_epochs, n_signals, n_times) | ~mne.Epochs | ~mne.time_frequency.EpochsSpectrum | ~mne.time_frequency.EpochsTFR
-        The data from which to compute connectivity. Can be epoched timeseries data as
-        an :term:`array-like` or :class:`~mne.Epochs` object, or Fourier coefficients
-        for each epoch as an :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object. If timeseries data, the spectral
+    data : array_like, shape (n_epochs, n_signals, n_times) | ~mne.Epochs | generator | ~mne.time_frequency.EpochsSpectrum | ~mne.time_frequency.EpochsTFR
+        The data from which to compute connectivity. Can be epoched time series data as
+        an array-like or :class:`mne.Epochs` object, or Fourier coefficients for each
+        epoch as an :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object. If time series data, the spectral
         information will be computed according to the spectral estimation mode (see the
-        ``mode`` parameter). If an :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object, existing spectral information
+        ``mode`` parameter). If an :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object, existing spectral information
         will be used and the ``mode`` parameter will be ignored.
 
-        Note that it is also possible to combine multiple timeseries signals by
+        Note that it is also possible to combine multiple time series signals by
         providing a list of tuples, e.g.: ::
 
             data = [(arr_0, stc_0), (arr_1, stc_1), (arr_2, stc_2)]
 
         which corresponds to 3 epochs where ``arr_*`` is an array with the same number
-        of time points as ``stc_*``. Data can also be a :class:`list`/:term:`generator`
-        of arrays, ``shape (n_signals, n_times)``, or a :class:`list`/:term:`generator`
-        of :class:`~mne.SourceEstimate` or :class:`~mne.VolSourceEstimate` objects.
+        of time points as ``stc_*``. Data can also be a list/generator of arrays, shape
+        ``(n_signals, n_times)``, or a list/generator of :class:`mne.SourceEstimate` or
+        :class:`mne.VolSourceEstimate` objects.
 
         .. versionchanged:: 0.8
-           Fourier coefficients stored in an :class:`~mne.time_frequency.EpochsSpectrum`
-           or :class:`~mne.time_frequency.EpochsTFR` object can also be passed in as
+           Fourier coefficients stored in an :class:`mne.time_frequency.EpochsSpectrum`
+           or :class:`mne.time_frequency.EpochsTFR` object can also be passed in as
            data. Storing Fourier coefficients in
-           :class:`~mne.time_frequency.EpochsSpectrum` objects requires ``mne >= 1.8``.
-           Storing multitaper weights in :class:`~mne.time_frequency.EpochsTFR` objects
+           :class:`mne.time_frequency.EpochsSpectrum` objects requires ``mne >= 1.8``.
+           Storing multitaper weights in :class:`mne.time_frequency.EpochsTFR` objects
            requires ``mne >= 1.10``.
     %(names)s
     method : str | list of str
-        Connectivity measure(s) to compute. These can be ``['coh', 'cohy',
-        'imcoh', 'cacoh', 'mic', 'mim', 'plv', 'ciplv', 'ppc', 'pli', 'dpli',
-        'wpli', 'wpli2_debiased', 'gc', 'gc_tr']``. These are:
+        Connectivity measure(s) to compute. These can be ``['coh', 'cohy', 'imcoh',
+        'cacoh', 'mic', 'mim', 'plv', 'ciplv', 'ppc', 'pli', 'dpli', 'wpli',
+        'wpli2_debiased', 'gc', 'gc_tr']``. These are:
 
         * %(coh)s
         * %(cohy)s
@@ -828,93 +828,88 @@ def spectral_connectivity_epochs(
         * %(gc)s
         * %(gc_tr)s
 
-        Multivariate methods (``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``)
-        cannot be called with the other methods.
-    indices : tuple of array | None
-        Two arrays with indices of connections for which to compute
-        connectivity. If a bivariate method is called, each array for the seeds
-        and targets should contain the channel indices for each bivariate
-        connection. If a multivariate method is called, each array for the
-        seeds and targets should consist of nested arrays containing
-        the channel indices for each multivariate connection. If ``None``,
-        connections between all channels are computed, unless a Granger
+        Multivariate methods (``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``) cannot be
+        called with the other methods.
+    indices : tuple of array_like | None
+        Two array-likes with indices of connections for which to compute connectivity.
+        If a bivariate method is called, each array for the seeds and targets should
+        contain the channel indices for each bivariate connection. If a multivariate
+        method is called, each array for the seeds and targets should consist of nested
+        arrays containing the channel indices for each multivariate connection. If
+        ``None``, connections between all channels are computed, unless a Granger
         causality method is called, in which case an error is raised.
-    sfreq : float
-        The sampling frequency. Required if data is an :term:`array-like`.
-    mode : str
-        Spectrum estimation mode can be either: 'multitaper', 'fourier', or
-        'cwt_morlet'. Ignored if ``data`` is an
-        :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object.
+    sfreq : float | None
+        The sampling frequency. Required if ``data`` is an array-like.
+    mode : ``'multitaper'`` | ``'fourier'`` | ``'cwt_morlet'``
+        Spectrum estimation mode. Ignored if ``data`` is an
+        :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object.
     fmin : float | tuple of float | None
-        The lower frequency of interest. Multiple bands are defined using
-        a tuple, e.g., (8., 20.) for two bands with 8 Hz and 20 Hz lower freq.
-        If ``None``, the frequency corresponding to 5 cycles based on the epoch
-        length is used. For example, with an epoch length of 1 sec, the lower
-        frequency would be 5 / 1 sec = 5 Hz.
+        The lower frequency of interest. Multiple bands are defined using a tuple, e.g.,
+        (8., 20.) for two bands with 8 Hz and 20 Hz lower freq. If ``None``, the
+        frequency corresponding to 5 cycles based on the epoch length is used. For
+        example, with an epoch length of 1 sec, the lower frequency would be 5 / 1 sec =
+        5 Hz.
     fmax : float | tuple of float
-        The upper frequency of interest. Multiple bands are defined using
-        a tuple, e.g. (13., 30.) for two band with 13 Hz and 30 Hz upper freq.
+        The upper frequency of interest. Multiple bands are defined using a tuple, e.g.,
+        (13., 30.) for two bands with 13 Hz and 30 Hz upper freq.
     fskip : int
-        Omit every "(fskip + 1)-th" frequency bin to decimate in frequency
-        domain.
+        Omit every "(fskip + 1)-th" frequency bin to decimate in frequency domain.
     faverage : bool
-        Average connectivity scores for each frequency band. If True,
-        the output freqs will be a list with arrays of the frequencies
-        that were averaged.
+        Average connectivity scores for each frequency band. If ``True``, the output
+        freqs will be a list with arrays of the frequencies that were averaged.
     tmin : float | None
-        Time to start connectivity estimation. Note: when ``data`` is an
-        :term:`array-like`, the first sample is assumed to be at time 0. For
-        :class:`~mne.Epochs`, the time information contained in the object is used to
-        compute the time indices. Ignored if ``data`` is an
-        :class:`~mne.time_frequency.EpochsSpectrum` object.
+        Time to start connectivity estimation. Note: when ``data`` is an array-like, the
+        first sample is assumed to be at time 0. For :class:`mne.Epochs`, the time
+        information contained in the object is used to compute the time indices. Ignored
+        if ``data`` is an :class:`mne.time_frequency.EpochsSpectrum` object.
     tmax : float | None
-        Time to end connectivity estimation. Note: when ``data`` is an
-        :term:`array-like`, the first sample is assumed to be at time 0. For
-        :class:`~mne.Epochs`, the time information contained in the object is used to
-        compute the time indices. Ignored if ``data`` is an
-        :class:`~mne.time_frequency.EpochsSpectrum` object.
+        Time to end connectivity estimation. Note: when ``data`` is an array-like, the
+        first sample is assumed to be at time 0. For :class:`mne.Epochs`, the time
+        information contained in the object is used to compute the time indices. Ignored
+        if ``data`` is an :class:`mne.time_frequency.EpochsSpectrum` object.
     mt_bandwidth : float | None
         The bandwidth of the multitaper windowing function in Hz. Only used in
-        'multitaper' mode. Ignored if ``data`` is an
-        :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object.
+        ``'multitaper'`` mode. Ignored if ``data`` is an
+        :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object.
     mt_adaptive : bool
         Use adaptive weights to combine the tapered spectra into PSD. Only used in
-        'multitaper' mode. Ignored if ``data`` is an
-        :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object.
+        ``'multitaper'`` mode. Ignored if ``data`` is an
+        :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object.
     mt_low_bias : bool
         Only use tapers with more than 90 percent spectral concentration within
-        bandwidth. Only used in 'multitaper' mode. Ignored if ``data`` is an
-        :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object.
-    cwt_freqs : array
-        Array of frequencies of interest. Only used in 'cwt_morlet' mode. Only
-        the frequencies within the range specified by ``fmin`` and ``fmax`` are
-        used. Ignored if ``data`` is an :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object.
-    cwt_n_cycles : float | array of float
-        Number of cycles. Fixed number or one per frequency. Only used in 'cwt_morlet'
-        mode. Ignored if ``data`` is an :class:`~mne.time_frequency.EpochsSpectrum` or
-        :class:`~mne.time_frequency.EpochsTFR` object.
+        bandwidth. Only used in ``'multitaper'`` mode. Ignored if ``data`` is an
+        :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object.
+    cwt_freqs : array_like
+        Array-like of frequencies of interest. Only used in ``'cwt_morlet'`` mode. Only
+        the frequencies within the range specified by ``fmin`` and ``fmax`` are used.
+        Ignored if ``data`` is an :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object.
+    cwt_n_cycles : float | array_like
+        Number of cycles. Fixed number or one per frequency. Only used in
+        ``'cwt_morlet'`` mode. Ignored if ``data`` is an
+        :class:`mne.time_frequency.EpochsSpectrum` or
+        :class:`mne.time_frequency.EpochsTFR` object.
     gc_n_lags : int
-        Number of lags to use for the vector autoregressive model when
-        computing Granger causality. Higher values increase computational cost,
-        but reduce the degree of spectral smoothing in the results. Only used
-        if ``method`` contains any of ``['gc', 'gc_tr']``.
-    rank : tuple of array | None
-        Two arrays with the rank to project the seed and target data to,
-        respectively, using singular value decomposition. If None, the rank of
-        the data is computed and projected to. Only used if ``method`` contains
-        any of ``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``.
-    n_components : int
-        Number of connectivity components to extract from the data. If an `int`, the
-        number of components must be <= the minimum rank of the seeds and targets. E.g.
+        Number of lags to use for the vector autoregressive model when computing Granger
+        causality. Higher values increase computational cost, but reduce the degree of
+        spectral smoothing in the results. Only used if ``method`` contains any of
+        ``['gc', 'gc_tr']``.
+    rank : tuple of array_like | None
+        Two array-likes with the rank to project the seed and target data to,
+        respectively, using singular value decomposition. If ``None``, the rank of the
+        data is computed and projected to. Only used if ``method`` contains any of
+        ``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``.
+    n_components : int | None
+        Number of connectivity components to extract from the data. If an int, the
+        number of components must be <= the minimum rank of the seeds and targets. E.g.,
         if the seed channels had a rank of 5 and the target channels had a rank of 3,
-        ``n_components`` must be <= 3. If `None`, the number of components equal to the
-        minimum rank of the seeds and targets is extracted (see the ``rank`` parameter).
-        Only used if ``method`` contains any of ``['cacoh', 'mic']``.
+        ``n_components`` must be <= 3. If ``None``, the number of components equal to
+        the minimum rank of the seeds and targets is extracted (see the ``rank``
+        parameter). Only used if ``method`` contains any of ``['cacoh', 'mic']``.
 
         .. versionadded:: 0.8
     block_size : int
@@ -932,14 +927,14 @@ def spectral_connectivity_epochs(
         connectivity measures if several connectivity measures are specified. The shape
         of the connectivity result will be:
 
-        - ``(n_cons, n_freqs)`` for multitaper or fourier modes
-        - ``(n_cons, n_freqs, n_times)`` for cwt_morlet mode
-        - ``(n_cons, n_comps, n_freqs (, n_times))`` for valid multivariate methods if
+        - ``(n_cons, n_freqs)`` for ``'multitaper'`` or ``'fourier'`` modes
+        - ``(n_cons, n_freqs, n_times)`` for ``'cwt_morlet'`` mode
+        - ``(n_cons, n_comps, n_freqs[, n_times])`` for valid multivariate methods if
           ``n_components > 1``
         - ``n_cons = n_signals ** 2`` for bivariate methods with ``indices=None``
         - ``n_cons = 1`` for multivariate methods with ``indices=None``
         - ``n_cons = len(indices[0])`` for bivariate and multivariate methods when
-          indices is supplied
+          ``indices`` is supplied
 
     See Also
     --------
@@ -949,31 +944,32 @@ def spectral_connectivity_epochs(
 
     Notes
     -----
-    Please note that the interpretation of the measures in this function
-    depends on the data and underlying assumptions and does not necessarily
-    reflect a causal relationship between brain regions.
-
-    These measures are not to be interpreted over time. Each Epoch passed into
-    the dataset is interpreted as an independent sample of the same
-    connectivity structure. Within each Epoch, it is assumed that the spectral
-    measure is stationary. The spectral measures implemented in this function
-    are computed across Epochs. **Thus, spectral measures computed with only
-    one Epoch will result in errorful values and spectral measures computed
-    with few Epochs will be unreliable.** Please see :func:`spectral_connectivity_time`
-    for time-resolved connectivity estimation.
-
-    The spectral densities can be estimated using a multitaper method with
-    digital prolate spheroidal sequence (DPSS) windows, a discrete Fourier
-    transform with Hanning windows, or a continuous wavelet transform using
-    Morlet wavelets. The spectral estimation mode is specified using the
-    "mode" parameter.
-
-    By default, the connectivity between all signals is computed (only
-    connections corresponding to the lower-triangular part of the connectivity
-    matrix). If one is only interested in the connectivity between some
-    signals, the "indices" parameter can be used. For example, to compute the
-    connectivity between the signal with index 0 and signals "2, 3, 4" (a total
-    of 3 connections) one can use the following::
+    Please note that the interpretation of the measures in this function depends on the
+    data and underlying assumptions and does not necessarily reflect a causal
+    relationship between brain regions.
+
+    These measures are not to be interpreted over time. Each epoch passed into the
+    dataset is interpreted as an independent sample of the same connectivity structure.
+    Within each epoch, it is assumed that the spectral measure is stationary. The
+    spectral measures implemented in this function are computed across epochs. **Thus,
+    spectral measures computed with only one epoch will result in errorful values and
+    spectral measures computed with few Epochs will be unreliable.** Please see
+    :func:`~mne_connectivity.spectral_connectivity_time` for time-resolved connectivity
+    estimation.
+
+    The spectral densities can be estimated using a multitaper method with digital
+    prolate spheroidal sequence (DPSS) windows, a discrete Fourier transform with
+    Hanning windows, or a continuous wavelet transform using Morlet wavelets. The
+    spectral estimation mode is specified using the ``mode`` parameter. Complex Welch,
+    multitaper, or Morlet coefficients can also be passed in as data in the form of
+    :class:`mne.time_frequency.EpochsSpectrum` or :class:`mne.time_frequency.EpochsTFR`
+    objects.
+
+    By default, the connectivity between all signals is computed (only connections
+    corresponding to the lower-triangular part of the connectivity matrix). If one is
+    only interested in the connectivity between some signals, the ``indices`` parameter
+    can be used. For example, to compute the connectivity between the signal with index
+    0 and signals "2, 3, 4" (a total of 3 connections) one can use the following::
 
         indices = (np.array([0, 0, 0]),    # row indices
                    np.array([2, 3, 4]))    # col indices
@@ -981,30 +977,29 @@ def spectral_connectivity_epochs(
         con = spectral_connectivity_epochs(data, method='coh',
                                            indices=indices, ...)
 
-    In this case con.get_data().shape = (3, n_freqs). The connectivity scores
-    are in the same order as defined indices.
+    In this case ``con.get_data().shape = (3, n_freqs)``. The connectivity scores are in
+    the same order as defined indices.
 
-    For multivariate methods, this is handled differently. If "indices" is
-    None, connectivity between all signals will be computed and a single
-    connectivity spectrum will be returned (this is not possible if a Granger
-    causality method is called). If "indices" is specified, seed and target
-    indices for each connection should be specified as nested array-likes. For
-    example, to compute the connectivity between signals (0, 1) -> (2, 3) and
-    (0, 1) -> (4, 5), indices should be specified as::
+    For multivariate methods, this is handled differently. If ``indices`` is ``None``,
+    connectivity between all signals will be computed and a single connectivity spectrum
+    will be returned (this is not possible if a Granger causality method is called). If
+    ``indices`` is specified, seed and target indices for each connection should be
+    specified as nested array-likes. For example, to compute the connectivity between
+    signals (0, 1) -> (2, 3) and (0, 1) -> (4, 5), indices should be specified as::
 
         indices = (np.array([[0, 1], [0, 1]]),  # seeds
                    np.array([[2, 3], [4, 5]]))  # targets
 
-    More information on working with multivariate indices and handling
-    connections where the number of seeds and targets are not equal can be
-    found in the :doc:`../auto_examples/handling_ragged_arrays` example.
+    More information on working with multivariate indices and handling connections where
+    the number of seeds and targets are not equal can be found in the
+    :doc:`../auto_examples/handling_ragged_arrays` example.
 
     **Supported Connectivity Measures**
 
-    The connectivity method(s) is specified using the "method" parameter. The
-    following methods are supported (note: ``E[]`` denotes average over
-    epochs). Multiple measures can be computed at once by using a list/tuple,
-    e.g., ``['coh', 'pli']`` to compute coherence and PLI.
+    The connectivity method(s) is specified using the ``method`` parameter. The
+    following methods are supported (note: ``E[]`` denotes average over epochs).
+    Multiple measures can be computed at once by using a list/tuple, e.g., ``['coh',
+    'pli']`` to compute coherence and PLI.
 
         'coh' : Coherence given by::
 
@@ -1024,98 +1019,90 @@ def spectral_connectivity_epochs(
             C = ----------------------
                 sqrt(E[Sxx] * E[Syy])
 
-        'cacoh' : Canonical Coherency (CaCoh) :footcite:`VidaurreEtAl2019`
-        given by:
+        'cacoh' : Canonical Coherency (CaCoh) :footcite:`VidaurreEtAl2019` given by:
 
             :math:`\textrm{CaCoh}=\Large{\frac{\boldsymbol{a}^T\boldsymbol{D}
             (\Phi)\boldsymbol{b}}{\sqrt{\boldsymbol{a}^T\boldsymbol{a}
             \boldsymbol{b}^T\boldsymbol{b}}}}`
 
-            where: :math:`\boldsymbol{D}(\Phi)` is the cross-spectral density
-            between seeds and targets transformed for a given phase angle
-            :math:`\Phi`; and :math:`\boldsymbol{a}` and :math:`\boldsymbol{b}`
-            are eigenvectors for the seeds and targets, such that
-            :math:`\boldsymbol{a}^T\boldsymbol{D}(\Phi)\boldsymbol{b}`
-            maximises coherency between the seeds and targets. Taking the
-            absolute value of the results gives maximised coherence.
+            where: :math:`\boldsymbol{D}(\Phi)` is the cross-spectral density between
+            seeds and targets transformed for a given phase angle :math:`\Phi`; and
+            :math:`\boldsymbol{a}` and :math:`\boldsymbol{b}` are eigenvectors for the
+            seeds and targets, such that :math:`\boldsymbol{a}^T\boldsymbol{D}(\Phi)
+            \boldsymbol{b}` maximises coherency between the seeds and targets. Taking
+            the absolute value of the results gives maximised coherence.
 
-        'mic' : Maximised Imaginary part of Coherency (MIC)
-        :footcite:`EwaldEtAl2012` given by:
+        'mic' : Maximised Imaginary part of Coherency (MIC) :footcite:`EwaldEtAl2012`
+        given by:
 
             :math:`\textrm{MIC}=\Large{\frac{\boldsymbol{\alpha}^T
             \boldsymbol{E \beta}}{\parallel\boldsymbol{\alpha}\parallel
             \parallel\boldsymbol{\beta}\parallel}}`
 
-            where: :math:`\boldsymbol{E}` is the imaginary part of the
-            transformed cross-spectral density between seeds and targets; and
-            :math:`\boldsymbol{\alpha}` and :math:`\boldsymbol{\beta}` are
-            eigenvectors for the seeds and targets, such that
-            :math:`\boldsymbol{\alpha}^T \boldsymbol{E \beta}` maximises the
-            imaginary part of coherency between the seeds and targets.
+            where: :math:`\boldsymbol{E}` is the imaginary part of the transformed
+            cross-spectral density between seeds and targets; and
+            :math:`\boldsymbol{\alpha}` and :math:`\boldsymbol{\beta}` are eigenvectors
+            for the seeds and targets, such that :math:`\boldsymbol{\alpha}^T
+            \boldsymbol{E \beta}` maximises the imaginary part of coherency between the
+            seeds and targets.
 
-        'mim' : Multivariate Interaction Measure (MIM)
-        :footcite:`EwaldEtAl2012` given by:
+        'mim' : Multivariate Interaction Measure (MIM) :footcite:`EwaldEtAl2012` given
+        by:
 
             :math:`\textrm{MIM}=tr(\boldsymbol{EE}^T)`
 
-            where :math:`\boldsymbol{E}` is the imaginary part of the
-            transformed cross-spectral density between seeds and targets.
+            where :math:`\boldsymbol{E}` is the imaginary part of the transformed
+            cross-spectral density between seeds and targets.
 
-        'plv' : Phase-Locking Value (PLV) :footcite:`LachauxEtAl1999` given
-        by::
+        'plv' : Phase-Locking Value (PLV) :footcite:`LachauxEtAl1999` given by::
 
             PLV = |E[Sxy/|Sxy|]|
 
-        'ciplv' : corrected imaginary PLV (ciPLV)
-        :footcite:`BrunaEtAl2018` given by::
+        'ciplv' : corrected imaginary PLV (ciPLV) :footcite:`BrunaEtAl2018` given by::
 
                              |E[Im(Sxy/|Sxy|)]|
             ciPLV = ------------------------------------
                      sqrt(1 - |E[real(Sxy/|Sxy|)]| ** 2)
 
-        'ppc' : Pairwise Phase Consistency (PPC), an unbiased estimator
-        of squared PLV :footcite:`VinckEtAl2010`.
+        'ppc' : Pairwise Phase Consistency (PPC), an unbiased estimator of squared PLV
+        :footcite:`VinckEtAl2010`.
 
         'pli' : Phase Lag Index (PLI) :footcite:`StamEtAl2007` given by::
 
             PLI = |E[sign(Im(Sxy))]|
 
-        'pli2_unbiased' : Unbiased estimator of squared PLI
-        :footcite:`VinckEtAl2011`.
+        'pli2_unbiased' : Unbiased estimator of squared PLI :footcite:`VinckEtAl2011`.
 
-        'dpli' : Directed Phase Lag Index (DPLI) :footcite:`StamEtAl2012`
-        given by (where H is the Heaviside function)::
+        'dpli' : Directed Phase Lag Index (DPLI) :footcite:`StamEtAl2012` given by
+        (where H is the Heaviside function)::
 
             DPLI = E[H(Im(Sxy))]
 
-        'wpli' : Weighted Phase Lag Index (WPLI) :footcite:`VinckEtAl2011`
-        given by::
+        'wpli' : Weighted Phase Lag Index (WPLI) :footcite:`VinckEtAl2011` given by::
 
                       |E[Im(Sxy)]|
             WPLI = ------------------
                       E[|Im(Sxy)|]
 
-        'wpli2_debiased' : Debiased estimator of squared WPLI
-        :footcite:`VinckEtAl2011`.
+        'wpli2_debiased' : Debiased estimator of squared WPLI :footcite:`VinckEtAl2011`.
 
-        'gc' : State-space Granger Causality (GC) :footcite:`BarnettSeth2015`
-        given by:
+        'gc' : State-space Granger Causality (GC) :footcite:`BarnettSeth2015` given by:
 
             :math:`GC = ln\Large{(\frac{\lvert\boldsymbol{S}_{tt}\rvert}{\lvert
             \boldsymbol{S}_{tt}-\boldsymbol{H}_{ts}\boldsymbol{\Sigma}_{ss
             \lvert t}\boldsymbol{H}_{ts}^*\rvert}})`
 
             where: :math:`s` and :math:`t` represent the seeds and targets,
-            respectively; :math:`\boldsymbol{H}` is the spectral transfer
-            function; :math:`\boldsymbol{\Sigma}` is the residuals matrix of
-            the autoregressive model; and :math:`\boldsymbol{S}` is
-            :math:`\boldsymbol{\Sigma}` transformed by :math:`\boldsymbol{H}`.
+            respectively; :math:`\boldsymbol{H}` is the spectral transfer function;
+            :math:`\boldsymbol{\Sigma}` is the residuals matrix of the autoregressive
+            model; and :math:`\boldsymbol{S}` is :math:`\boldsymbol{\Sigma}` transformed
+            by :math:`\boldsymbol{H}`.
 
         'gc_tr' : State-space GC on time-reversed signals
-        :footcite:`BarnettSeth2015,WinklerEtAl2016` given by the same equation
-        as for 'gc', but where the autocovariance sequence from which the
-        autoregressive model is produced is transposed to mimic the reversal of
-        the original signal in time :footcite:`HaufeEtAl2012`.
+        :footcite:`BarnettSeth2015,WinklerEtAl2016` given by the same equation as for
+        ``'gc'``, but where the autocovariance sequence from which the autoregressive
+        model is produced is transposed to mimic the reversal of the original signal in
+        time :footcite:`HaufeEtAl2012`.
 
     References
     ----------
diff --git a/mne_connectivity/spectral/epochs_multivariate.py b/mne_connectivity/spectral/epochs_multivariate.py
index 4257bd402..ada7234b3 100644
--- a/mne_connectivity/spectral/epochs_multivariate.py
+++ b/mne_connectivity/spectral/epochs_multivariate.py
@@ -397,14 +397,14 @@ def _invsqrtm(self, C_r, n_seeds):
 
         Parameters
         ----------
-        C_r : np.ndarray, shape=(n_freqs, n_times, n_channels, n_channels)
+        C_r : array, shape (n_freqs, n_times, n_channels, n_channels)
             Real part of the CSD. Expected to be symmetric and non-singular.
         n_seeds : int
             Number of seed channels for the connection.
 
         Returns
         -------
-        T : np.ndarray, shape=(n_freqs, n_times, n_channels, n_channels)
+        T : array, shape (n_freqs, n_times, n_channels, n_channels)
             Inverse square root of the real-valued CSD. Name comes from Ewald
             et al. (2012).
 
@@ -667,13 +667,12 @@ def _first_optimise_phi(self, C_ab, T_aa, T_bb):
     def _final_optimise_phi(self, C_ab, T_aa, T_bb, max_coh, max_phis):
         """Fine-tune the angle at which coherence is maximised.
 
-        Uses a 2nd order Taylor expansion to approximate change in coherence
-        w.r.t. phi, and determining the next phi to evaluate coherence on (over
-        a total of 10 iterations).
+        Uses a 2nd order Taylor expansion to approximate change in coherence w.r.t. phi,
+        and determining the next phi to evaluate coherence on (over a total of 10
+        iterations).
 
-        Depending on how the new phi affects coherence, the step size for the
-        subsequent iteration is adjusted, like that in the Levenberg-Marquardt
-        algorithm.
+        Depending on how the new phi affects coherence, the step size for the subsequent
+        iteration is adjusted, like that in the Levenberg-Marquardt algorithm.
 
         Each time-freq. entry of coherence has its own corresponding phi.
         """
@@ -792,7 +791,7 @@ def _deflate_csd(self, C, W_a, W_b, n_seeds):
 
         See "Methods - Extracting further source pairs" of Dähne et al. (2014),
         NeuroImage, DOI: 10.1016/j.neuroimage.2014.03.075, for an example of applying
-        this approach to timeseries data.
+        this approach to time series data.
         """
         # get orthogonal basis space for filters
         # (streamlined version of scipy.linalg.null_space() suited for our purposes)
@@ -1065,8 +1064,8 @@ def _full_var_to_iss(self, A_f):
         """Compute innovations-form parameters for a state-space model.
 
         Parameters computed from a full VAR model using Aoki's method. For a
-        non-moving-average full VAR model, the state-space parameter C
-        (observation matrix) is identical to AF of the VAR model.
+        non-moving-average full VAR model, the state-space parameter C (observation
+        matrix) is identical to AF of the VAR model.
 
         See: Barnett, L. & Seth, A.K., 2015, Physical Review, DOI:
         10.1103/PhysRevE.91.040101.
@@ -1112,14 +1111,13 @@ def _iss_to_ugc(self, A, C, K, V, seeds, targets):
     def _iss_to_tf(self, A, C, K, z):
         """Compute transfer function for innovations-form state-space params.
 
-        In the frequency domain, the back-shift operator, z, is a vector of
-        points on a unit circle in the complex plane. z = e^-iw, where -pi < w
-        <= pi.
+        In the frequency domain, the back-shift operator, z, is a vector of points on a
+        unit circle in the complex plane. z = e^-iw, where -pi < w <= pi.
 
-        A note on efficiency: solving over the 4D time-freq. tensor is slower
-        than looping over times and freqs when n_times and n_freqs high, and
-        when n_times and n_freqs low, looping over times and freqs very fast
-        anyway (plus tensor solving doesn't allow for parallelisation).
+        A note on efficiency: solving over the 4D time-freq. tensor is slower than
+        looping over times and freqs when n_times and n_freqs high, and when n_times and
+        n_freqs low, looping over times and freqs very fast anyway (plus tensor solving
+        doesn't allow for parallelisation).
 
         See: Barnett, L. & Seth, A.K., 2015, Physical Review, DOI:
         10.1103/PhysRevE.91.040101.
@@ -1147,9 +1145,9 @@ def _iss_to_tf(self, A, C, K, z):
     def _partial_covar(self, V, seeds, targets):
         """Compute partial covariance of a matrix.
 
-        Given a covariance matrix V, the partial covariance matrix of V between
-        indices i and j, given k (V_ij|k), is equivalent to V_ij - V_ik *
-        V_kk^-1 * V_kj. In this case, i and j are seeds, and k are targets.
+        Given a covariance matrix V, the partial covariance matrix of V between indices
+        i and j, given k (V_ij|k), is equivalent to V_ij - V_ik * V_kk^-1 * V_kj. In
+        this case, i and j are seeds, and k are targets.
 
         See: Barnett, L. & Seth, A.K., 2015, Physical Review, DOI:
         10.1103/PhysRevE.91.040101.
diff --git a/mne_connectivity/spectral/smooth.py b/mne_connectivity/spectral/smooth.py
index e67f4ab28..6505f4c1c 100644
--- a/mne_connectivity/spectral/smooth.py
+++ b/mne_connectivity/spectral/smooth.py
@@ -9,9 +9,8 @@ def _create_kernel(sm_times, sm_freqs, kernel="hanning"):
     Parameters
     ----------
     sm_times : int, array_like
-        Number of points to consider for the temporal smoothing,
-        if it is an array it will be considered that the kernel
-        if frequency dependent.
+        Number of points to consider for the temporal smoothing, if it is an array it
+        will be considered that the kernel if frequency dependent.
     sm_freqs : int
         Number of points to consider for the frequency smoothing
     kernel : {'square', 'hanning'}
@@ -71,8 +70,8 @@ def __pad_kernel(s):
 def _smooth_spectra(spectra, kernel, scale=False, decim=1):
     """Smoothing spectra.
 
-    This function assumes that the frequency and time axis are respectively
-    located at positions (..., freqs, times).
+    This function assumes that the frequency and time axis are respectively located at
+    positions (..., freqs, times).
 
     Parameters
     ----------
diff --git a/mne_connectivity/spectral/tests/test_spectral.py b/mne_connectivity/spectral/tests/test_spectral.py
index bb5abc268..a5d8acc52 100644
--- a/mne_connectivity/spectral/tests/test_spectral.py
+++ b/mne_connectivity/spectral/tests/test_spectral.py
@@ -58,9 +58,9 @@ def create_test_dataset(
 
     Returns
     -------
-    data : np.ndarray of shape (n_epochs, n_signals, n_times)
+    data : array of shape (n_epochs, n_signals, n_times)
         The epoched dataset.
-    times_data : np.ndarray of shape (n_times, )
+    times_data : array of shape (n_times, )
         The times at which each sample of the ``data`` occurs at.
     """
     # Use a case known to have no spurious correlations (it would bad if
diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py
index 1876fc833..745965dba 100644
--- a/mne_connectivity/spectral/time.py
+++ b/mne_connectivity/spectral/time.py
@@ -47,13 +47,13 @@ def spectral_connectivity_time(
     fmax=None,
     fskip=0,
     faverage=False,
-    sm_times=0,
+    sm_times=0.0,
     sm_freqs=1,
     sm_kernel="hanning",
-    padding=0,
+    padding=0.0,
     mode="cwt_morlet",
     mt_bandwidth=None,
-    n_cycles=7,
+    n_cycles=7.0,
     gc_n_lags=40,
     rank=None,
     n_components=1,
@@ -63,39 +63,37 @@ def spectral_connectivity_time(
 ):
     r"""Compute time-frequency-domain connectivity measures.
 
-    This function computes spectral connectivity over time from epoched data.
-    The data may consist of a single epoch.
+    This function computes spectral connectivity over time from epoched data. The data
+    may consist of a single epoch.
 
-    The connectivity method(s) are specified using the ``method`` parameter.
-    All methods are based on time-resolved estimates of the cross- and
-    power spectral densities (CSD/PSD) Sxy and Sxx, Syy.
+    The connectivity method(s) are specified using the ``method`` parameter. All methods
+    are based on time-resolved estimates of the cross- and power spectral densities
+    (CSD/PSD) Sxy and Sxx, Syy.
 
     Parameters
     ----------
     data : array_like, shape (n_epochs, n_signals, n_times) | ~mne.Epochs | ~mne.time_frequency.EpochsTFR
-        The data from which to compute connectivity. Can be epoched timeseries data as
-        an :term:`array-like` or :class:`~mne.Epochs` object, or Fourier coefficients
-        for each epoch as an :class:`~mne.time_frequency.EpochsTFR` object. If
-        timeseries data, the spectral information will be computed according to the
-        spectral estimation mode (see the ``mode`` parameter). If an
-        :class:`~mne.time_frequency.EpochsTFR` object, existing spectral information
-        will be used and the ``mode`` parameter will be ignored.
+        The data from which to compute connectivity. Can be epoched time series data as
+        an array-like or :class:`mne.Epochs` object, or Fourier coefficients for each
+        epoch as an :class:`mne.time_frequency.EpochsTFR` object. If time series data,
+        the spectral information will be computed according to the spectral estimation
+        mode (see the ``mode`` parameter). If an :class:`mne.time_frequency.EpochsTFR`
+        object, existing spectral information will be used and the ``mode`` parameter
+        will be ignored.
 
         .. versionchanged:: 0.8
-           Fourier coefficients stored in an :class:`~mne.time_frequency.EpochsTFR`
+           Fourier coefficients stored in an :class:`mne.time_frequency.EpochsTFR`
            object can also be passed in as data. Storing multitaper weights in
-           :class:`~mne.time_frequency.EpochsTFR` objects requires ``mne >= 1.10``.
+           :class:`mne.time_frequency.EpochsTFR` objects requires ``mne >= 1.10``.
     freqs : array_like | None
-        Array of frequencies of interest for time-frequency decomposition. Only the
+        Array-like of frequencies of interest for time-frequency decomposition. Only the
         frequencies within the range specified by ``fmin`` and ``fmax`` are used. If
-        ``data`` is an :term:`array-like` or :class:`~mne.Epochs` object, the
-        frequencies must be specified. If ``data`` is an
-        :class:`~mne.time_frequency.EpochsTFR` object, ``data.freqs`` is used and this
-        parameter is ignored.
+        ``data`` is an array-like or :class:`mne.Epochs` object, the frequencies must
+        be specified. If ``data`` is an :class:`mne.time_frequency.EpochsTFR` object,
+        ``data.freqs`` is used and this parameter is ignored.
     method : str | list of str
-        Connectivity measure(s) to compute. These can be ``['coh', 'cacoh',
-        'mic', 'mim', 'plv', 'ciplv', 'pli', 'wpli', 'gc', 'gc_tr']``. These
-        are:
+        Connectivity measure(s) to compute. These can be ``['coh', 'cacoh', 'mic',
+        'mim', 'plv', 'ciplv', 'pli', 'wpli', 'gc', 'gc_tr']``. These are:
 
         * %(coh)s
         * %(cacoh)s
@@ -108,94 +106,90 @@ def spectral_connectivity_time(
         * %(gc)s
         * %(gc_tr)s
 
-        Multivariate methods (``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``)
-        cannot be called with the other methods.
+        Multivariate methods (``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``) cannot be
+        called with the other methods.
     average : bool
-        Average connectivity scores over epochs. If ``True``, output will be
-        an instance of :class:`SpectralConnectivity`, otherwise
-        :class:`EpochSpectralConnectivity`.
+        Average connectivity scores over epochs. If ``True``, output will be an instance
+        of :class:`SpectralConnectivity`, otherwise :class:`EpochSpectralConnectivity`.
     indices : tuple of array_like | None
-        Two arrays with indices of connections for which to compute
-        connectivity. If a bivariate method is called, each array for the seeds
-        and targets should contain the channel indices for the each bivariate
-        connection. If a multivariate method is called, each array for the
-        seeds and targets should consist of nested arrays containing
-        the channel indices for each multivariate connection. If None,
-        connections between all channels are computed, unless a Granger
+        Two array-likes with indices of connections for which to compute connectivity.
+        If a bivariate method is called, each array for the seeds and targets should
+        contain the channel indices for the each bivariate connection. If a multivariate
+        method is called, each array for the seeds and targets should consist of nested
+        arrays containing the channel indices for each multivariate connection. If
+        ``None``, connections between all channels are computed, unless a Granger
         causality method is called, in which case an error is raised.
-    sfreq : float
-        The sampling frequency. Required if ``data`` is not an :class:`~mne.Epochs` or
-        :class:`~mne.time_frequency.EpochsTFR` object.
+    sfreq : float | None
+        The sampling frequency. Required if ``data`` is not an :class:`mne.Epochs` or
+        :class:`mne.time_frequency.EpochsTFR` object.
     fmin : float | tuple of float | None
-        The lower frequency of interest. Multiple bands are defined using
-        a tuple, e.g., ``(8., 20.)`` for two bands with 8 Hz and 20 Hz lower
-        bounds. If `None`, the lowest frequency in ``freqs`` is used.
+        The lower frequency of interest. Multiple bands are defined using a tuple, e.g.,
+        ``(8., 20.)`` for two bands with 8 Hz and 20 Hz lower bounds. If ``None``, the
+        lowest frequency in ``freqs`` is used.
     fmax : float | tuple of float | None
-        The upper frequency of interest. Multiple bands are defined using
-        a tuple, e.g. ``(13., 30.)`` for two band with 13 Hz and 30 Hz upper
-        bounds. If `None`, the highest frequency in ``freqs`` is used.
+        The upper frequency of interest. Multiple bands are defined using a tuple, e.g.
+        ``(13., 30.)`` for two band with 13 Hz and 30 Hz upper bounds. If ``None``, the
+        highest frequency in ``freqs`` is used.
     fskip : int
-        Omit every ``(fskip + 1)``-th frequency bin to decimate in frequency
-        domain.
+        Omit every ``(fskip + 1)``-th frequency bin to decimate in frequency domain.
     faverage : bool
-        Average connectivity scores for each frequency band. If `True`,
-        the output ``freqs`` will be an array of the median frequencies of each
-        band.
+        Average connectivity scores for each frequency band. If ``True``, the output
+        ``freqs`` will be an array of the median frequencies of each band.
     sm_times : float
-        Amount of time to consider for the temporal smoothing in seconds.
-        If zero, no temporal smoothing is applied.
+        Amount of time to consider for the temporal smoothing in seconds. If 0, no
+        temporal smoothing is applied.
     sm_freqs : int
-        Number of points for frequency smoothing. By default, 1 is used which
-        is equivalent to no smoothing.
-    sm_kernel : {'square', 'hanning'}
-        Smoothing kernel type. Choose either 'square' or 'hanning'.
+        Number of points for frequency smoothing. By default, 1 is used which is
+        equivalent to no smoothing.
+    sm_kernel : ``'square'`` | ``'hanning'``
+        Smoothing kernel type. For ``'hanning'``, see :func:`numpy.hanning`.
     padding : float
-        Amount of time to consider as padding at the beginning and end of each
-        epoch in seconds. See Notes for more information.
-    mode : str
-        Time-frequency decomposition method. Can be either: ``'multitaper'``, or
-        ``'cwt_morlet'``. See :func:`mne.time_frequency.tfr_array_multitaper` and
+        Amount of time to consider as padding at the beginning and end of each epoch in
+        seconds. See Notes for more information.
+    mode : ``'multitaper'`` | ``'cwt_morlet'``
+        Time-frequency decomposition method. See
+        :func:`mne.time_frequency.tfr_array_multitaper` and
         :func:`mne.time_frequency.tfr_array_morlet` for reference. Ignored if ``data``
-        is an :class:`~mne.time_frequency.EpochsTFR` object.
+        is an :class:`mne.time_frequency.EpochsTFR` object.
     mt_bandwidth : float | None
         Product between the temporal window length (in seconds) and the full frequency
         bandwidth (in Hz). This product can be seen as the surface of the window on the
         time/frequency plane and controls the frequency bandwidth (thus the frequency
         resolution) and the number of good tapers. See
         :func:`mne.time_frequency.tfr_array_multitaper` documentation. Ignored if
-        ``data`` is an :class:`~mne.time_frequency.EpochsTFR` object.
-    n_cycles : float | array_like of float
+        ``data`` is an :class:`mne.time_frequency.EpochsTFR` object.
+    n_cycles : float | array_like
         Number of cycles in the wavelet, either a fixed number or one per frequency. The
         number of cycles ``n_cycles`` and the frequencies of interest ``freqs`` define
         the temporal window length. For details, see
         :func:`mne.time_frequency.tfr_array_multitaper` and
         :func:`mne.time_frequency.tfr_array_morlet` documentation. Ignored if ``data``
-        is an :class:`~mne.time_frequency.EpochsTFR` object.
+        is an :class:`mne.time_frequency.EpochsTFR` object.
     gc_n_lags : int
-        Number of lags to use for the vector autoregressive model when
-        computing Granger causality. Higher values increase computational cost,
-        but reduce the degree of spectral smoothing in the results. Only used
-        if ``method`` contains any of ``['gc', 'gc_tr']``.
-    rank : tuple of array | None
-        Two arrays with the rank to project the seed and target data to,
-        respectively, using singular value decomposition. If `None`, the rank
-        of the data is computed and projected to. Only used if ``method``
-        contains any of ``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``.
-    n_components : int
-        Number of connectivity components to extract from the data. If an `int`, the
-        number of components must be <= the minimum rank of the seeds and targets. E.g.
+        Number of lags to use for the vector autoregressive model when computing Granger
+        causality. Higher values increase computational cost, but reduce the degree of
+        spectral smoothing in the results. Only used if ``method`` contains any of
+        ``['gc', 'gc_tr']``.
+    rank : tuple of array_like | None
+        Two array-likes with the rank to project the seed and target data to,
+        respectively, using singular value decomposition. If ``None``, the rank of the
+        data is computed and projected to. Only used if ``method`` contains any of
+        ``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``.
+    n_components : int | None
+        Number of connectivity components to extract from the data. If an int, the
+        number of components must be <= the minimum rank of the seeds and targets. E.g.,
         if the seed channels had a rank of 5 and the target channels had a rank of 3,
-        ``n_components`` must be <= 3. If `None`, the number of components equal to the
-        minimum rank of the seeds and targets is extracted (see the ``rank`` parameter).
-        Only used if ``method`` contains any of ``['cacoh', 'mic']``.
+        ``n_components`` must be <= 3. If ``None``, the number of components equal to
+        the minimum rank of the seeds and targets is extracted (see the ``rank``
+        parameter). Only used if ``method`` contains any of ``['cacoh', 'mic']``.
 
         .. versionadded:: 0.8
     decim : int
-        To reduce memory usage, decimation factor after time-frequency
-        decomposition. Returns ``tfr[…, ::decim]``.
+        To reduce memory usage, decimation factor after time-frequency decomposition.
+        Returns ``tfr[…, ::decim]``.
     n_jobs : int
-        Number of connections to compute in parallel. Memory mapping must be
-        activated. Please see the Notes section for details.
+        Number of connections to compute in parallel. Memory mapping must be activated.
+        Please see the Notes section for details.
     %(verbose)s
 
     Returns
@@ -204,12 +198,14 @@ def spectral_connectivity_time(
         Computed connectivity measure(s). An instance of
         :class:`EpochSpectralConnectivity`, :class:`SpectralConnectivity`, or a list of
         instances corresponding to connectivity measures if several connectivity
-        measures are specified. The shape of each connectivity dataset is ([n_epochs,]
-        n_cons, [n_comps,] n_freqs). ``n_comps`` is present for valid multivariate
-        methods if ``n_components > 1``. When "indices" is None and a bivariate method
-        is called, "n_cons = n_signals ** 2", or if a multivariate method is called
-        "n_cons = 1". When "indices" is specified, "n_con = len(indices[0])" for
-        bivariate and multivariate methods.
+        measures are specified. The shape of each connectivity dataset is
+        ``([n_epochs,] n_cons, [n_comps,] n_freqs)``:
+
+        - ``n_comps`` is present for valid multivariate methods if ``n_components > 1``
+        - When ``indices`` is ``None`` and a bivariate method is called, ``n_cons =
+          n_signals ** 2``, or if a multivariate method is called ``n_cons = 1``
+        - When ``indices`` is specified, ``n_con = len(indices[0])`` for bivariate and
+          multivariate methods.
 
     See Also
     --------
@@ -219,39 +215,39 @@ def spectral_connectivity_time(
 
     Notes
     -----
-    Please note that the interpretation of the measures in this function
-    depends on the data and underlying assumptions and does not necessarily
-    reflect a causal relationship between brain regions.
-
-    The connectivity measures are computed over time within each epoch and
-    optionally averaged over epochs. High connectivity values indicate that
-    the phase coupling (interpreted as estimated connectivity) differences
-    between signals stay consistent over time.
-
-    The spectral densities can be estimated using a multitaper method with
-    digital prolate spheroidal sequence (DPSS) windows, or a continuous wavelet
-    transform using Morlet wavelets. The spectral estimation mode is specified
-    using the ``mode`` parameter.
-
-    When using the multitaper spectral estimation method, the
-    cross-spectral density is computed separately for each taper and aggregated
-    using a weighted average, where the weights correspond to the concentration
-    ratios between the DPSS windows.
-
-    Spectral estimation using multitaper or Morlet wavelets introduces edge
-    effects that depend on the length of the wavelet. To remove edge effects,
-    the parameter ``padding`` can be used to prune the edges of the signal.
-    Please see the documentation of
-    :func:`mne.time_frequency.tfr_array_multitaper` and
-    :func:`mne.time_frequency.tfr_array_morlet` for details on wavelet length
-    (i.e., time window length).
-
-    By default, the connectivity between all signals is computed (only
-    connections corresponding to the lower-triangular part of the connectivity
-    matrix). If one is only interested in the connectivity between some
-    signals, the "indices" parameter can be used. For example, to compute the
-    connectivity between the signal with index 0 and signals "2, 3, 4" (a total
-    of 3 connections) one can use the following::
+    Please note that the interpretation of the measures in this function depends on the
+    data and underlying assumptions and does not necessarily reflect a causal
+    relationship between brain regions.
+
+    The connectivity measures are computed over time within each epoch and optionally
+    averaged over epochs. High connectivity values indicate that the phase coupling
+    (interpreted as estimated connectivity) differences between signals stay consistent
+    over time.
+
+    The spectral densities can be estimated using a multitaper method with digital
+    prolate spheroidal sequence (DPSS) windows, or a continuous wavelet transform using
+    Morlet wavelets. The spectral estimation mode is specified using the ``mode``
+    parameter.
+
+    When using the multitaper spectral estimation method, the cross-spectral density is
+    computed separately for each taper and aggregated using a weighted average, where
+    the weights correspond to the concentration ratios between the DPSS windows.
+
+    Spectral estimation using multitaper or Morlet wavelets introduces edge effects that
+    depend on the length of the wavelet. To remove edge effects, the parameter
+    ``padding`` can be used to prune the edges of the signal. Please see the
+    documentation of :func:`mne.time_frequency.tfr_array_multitaper` and
+    :func:`mne.time_frequency.tfr_array_morlet` for details on wavelet length (i.e.,
+    time window length).
+
+    Complex multitaper, or Morlet coefficients can also be passed in as data in the form
+    of :class:`mne.time_frequency.EpochsTFR` objects.
+
+    By default, the connectivity between all signals is computed (only connections
+    corresponding to the lower-triangular part of the connectivity matrix). If one is
+    only interested in the connectivity between some signals, the ``indices`` parameter
+    can be used. For example, to compute the connectivity between the signal with index
+    0 and signals "2, 3, 4" (a total of 3 connections) one can use the following::
 
         indices = (np.array([0, 0, 0]),    # row indices
                    np.array([2, 3, 4]))    # col indices
@@ -259,30 +255,29 @@ def spectral_connectivity_time(
         con = spectral_connectivity_time(data, method='coh',
                                          indices=indices, ...)
 
-    In this case ``con.get_data().shape = (3, n_freqs)``. The connectivity
-    scores are in the same order as defined indices.
+    In this case ``con.get_data().shape = (3, n_freqs)``. The connectivity scores are in
+    the same order as defined indices.
 
-    For multivariate methods, this is handled differently. If "indices" is
-    None, connectivity between all signals will be computed and a single
-    connectivity spectrum will be returned (this is not possible if a Granger
-    causality method is called). If "indices" is specified, seed and target
-    indices for each connection should be specified as nested array-likes. For
-    example, to compute the connectivity between signals (0, 1) -> (2, 3) and
-    (0, 1) -> (4, 5), indices should be specified as::
+    For multivariate methods, this is handled differently. If ``indices`` is ``None``,
+    connectivity between all signals will be computed and a single connectivity spectrum
+    will be returned (this is not possible if a Granger causality method is called). If
+    ``indices`` is specified, seed and target indices for each connection should be
+    specified as nested array-likes. For example, to compute the connectivity between
+    signals (0, 1) -> (2, 3) and (0, 1) -> (4, 5), indices should be specified as::
 
         indices = (np.array([[0, 1], [0, 1]]),  # seeds
                    np.array([[2, 3], [4, 5]]))  # targets
 
-    More information on working with multivariate indices and handling
-    connections where the number of seeds and targets are not equal can be
-    found in the :doc:`../auto_examples/handling_ragged_arrays` example.
+    More information on working with multivariate indices and handling connections where
+    the number of seeds and targets are not equal can be found in the
+    :doc:`../auto_examples/handling_ragged_arrays` example.
 
     **Supported Connectivity Measures**
 
     The connectivity method(s) is specified using the ``method`` parameter. The
-    following methods are supported (note: ``E[]`` denotes average over
-    epochs). Multiple measures can be computed at once by using a list/tuple,
-    e.g., ``['coh', 'pli']`` to compute coherence and PLI.
+    following methods are supported (note: ``E[]`` denotes average over epochs).
+    Multiple measures can be computed at once by using a list/tuple, e.g., ``['coh',
+    'pli']`` to compute coherence and PLI.
 
         'coh' : Coherence given by::
 
@@ -290,50 +285,46 @@ def spectral_connectivity_time(
             C = ---------------------
                 sqrt(E[Sxx] * E[Syy])
 
-        'cacoh' : Canonical Coherency (CaCoh) :footcite:`VidaurreEtAl2019`
-        given by:
+        'cacoh' : Canonical Coherency (CaCoh) :footcite:`VidaurreEtAl2019` given by:
 
             :math:`\textrm{CaCoh}=\Large{\frac{\boldsymbol{a}^T\boldsymbol{D}
             (\Phi)\boldsymbol{b}}{\sqrt{\boldsymbol{a}^T\boldsymbol{a}
             \boldsymbol{b}^T\boldsymbol{b}}}}`
 
-            where: :math:`\boldsymbol{D}(\Phi)` is the cross-spectral density
-            between seeds and targets transformed for a given phase angle
-            :math:`\Phi`; and :math:`\boldsymbol{a}` and :math:`\boldsymbol{b}`
-            are eigenvectors for the seeds and targets, such that
-            :math:`\boldsymbol{a}^T\boldsymbol{D}(\Phi)\boldsymbol{b}`
-            maximises coherency between the seeds and targets. Taking the
-            absolute value of the results gives maximised coherence.
+            where: :math:`\boldsymbol{D}(\Phi)` is the cross-spectral density between
+            seeds and targets transformed for a given phase angle :math:`\Phi`; and
+            :math:`\boldsymbol{a}` and :math:`\boldsymbol{b}` are eigenvectors for the
+            seeds and targets, such that :math:`\boldsymbol{a}^T\boldsymbol{D}(\Phi)
+            \boldsymbol{b}` maximises coherency between the seeds and targets. Taking
+            the absolute value of the results gives maximised coherence.
 
-        'mic' : Maximised Imaginary part of Coherency (MIC)
-        :footcite:`EwaldEtAl2012` given by:
+        'mic' : Maximised Imaginary part of Coherency (MIC) :footcite:`EwaldEtAl2012`
+        given by:
 
             :math:`\textrm{MIC}=\Large{\frac{\boldsymbol{\alpha}^T
             \boldsymbol{E \beta}}{\parallel\boldsymbol{\alpha}\parallel
             \parallel\boldsymbol{\beta}\parallel}}`
 
-            where: :math:`\boldsymbol{E}` is the imaginary part of the
-            transformed cross-spectral density between seeds and targets; and
-            :math:`\boldsymbol{\alpha}` and :math:`\boldsymbol{\beta}` are
-            eigenvectors for the seeds and targets, such that
-            :math:`\boldsymbol{\alpha}^T \boldsymbol{E \beta}` maximises the
-            imaginary part of coherency between the seeds and targets.
+            where: :math:`\boldsymbol{E}` is the imaginary part of the transformed
+            cross-spectral density between seeds and targets; and
+            :math:`\boldsymbol{\alpha}` and :math:`\boldsymbol{\beta}` are eigenvectors
+            for the seeds and targets, such that :math:`\boldsymbol{\alpha}^T
+            \boldsymbol{E \beta}` maximises the imaginary part of coherency between the
+            seeds and targets.
 
-        'mim' : Multivariate Interaction Measure (MIM)
-        :footcite:`EwaldEtAl2012` given by:
+        'mim' : Multivariate Interaction Measure (MIM) :footcite:`EwaldEtAl2012` given
+        by:
 
             :math:`\textrm{MIM}=tr(\boldsymbol{EE}^T)`
 
-            where :math:`\boldsymbol{E}` is the imaginary part of the
-            transformed cross-spectral density between seeds and targets.
+            where :math:`\boldsymbol{E}` is the imaginary part of the transformed
+            cross-spectral density between seeds and targets.
 
-        'plv' : Phase-Locking Value (PLV) :footcite:`LachauxEtAl1999` given
-        by::
+        'plv' : Phase-Locking Value (PLV) :footcite:`LachauxEtAl1999` given by::
 
             PLV = |E[Sxy/|Sxy|]|
 
-        'ciplv' : Corrected imaginary PLV (ciPLV) :footcite:`BrunaEtAl2018`
-        given by::
+        'ciplv' : Corrected imaginary PLV (ciPLV) :footcite:`BrunaEtAl2018` given by::
 
                              |E[Im(Sxy/|Sxy|)]|
             ciPLV = ------------------------------------
@@ -343,49 +334,45 @@ def spectral_connectivity_time(
 
             PLI = |E[sign(Im(Sxy))]|
 
-        'wpli' : Weighted Phase Lag Index (WPLI) :footcite:`VinckEtAl2011`
-        given by::
+        'wpli' : Weighted Phase Lag Index (WPLI) :footcite:`VinckEtAl2011` given by::
 
                       |E[Im(Sxy)]|
             WPLI = ------------------
                       E[|Im(Sxy)|]
 
-        'gc' : State-space Granger Causality (GC) :footcite:`BarnettSeth2015`
-        given by:
+        'gc' : State-space Granger Causality (GC) :footcite:`BarnettSeth2015` given by:
 
             :math:`GC = ln\Large{(\frac{\lvert\boldsymbol{S}_{tt}\rvert}{\lvert
             \boldsymbol{S}_{tt}-\boldsymbol{H}_{ts}\boldsymbol{\Sigma}_{ss
             \lvert t}\boldsymbol{H}_{ts}^*\rvert}})`
 
             where: :math:`s` and :math:`t` represent the seeds and targets,
-            respectively; :math:`\boldsymbol{H}` is the spectral transfer
-            function; :math:`\boldsymbol{\Sigma}` is the residuals matrix of
-            the autoregressive model; and :math:`\boldsymbol{S}` is
-            :math:`\boldsymbol{\Sigma}` transformed by :math:`\boldsymbol{H}`.
+            respectively; :math:`\boldsymbol{H}` is the spectral transfer function;
+            :math:`\boldsymbol{\Sigma}` is the residuals matrix of the autoregressive
+            model; and :math:`\boldsymbol{S}` is :math:`\boldsymbol{\Sigma}` transformed
+            by :math:`\boldsymbol{H}`.
 
         'gc_tr' : State-space GC on time-reversed signals
-        :footcite:`BarnettSeth2015,WinklerEtAl2016` given by the same equation
-        as for 'gc', but where the autocovariance sequence from which the
-        autoregressive model is produced is transposed to mimic the reversal of
-        the original signal in time :footcite:`HaufeEtAl2012`.
-
-    Parallel computation can be activated by setting the ``n_jobs`` parameter.
-    Under the hood, this utilizes the ``joblib`` library. For effective
-    parallelization, you should activate memory mapping in MNE-Python by
-    setting ``MNE_MEMMAP_MIN_SIZE`` and ``MNE_CACHE_DIR``. Activating memory
-    mapping will make ``joblib`` store arrays greater than the minimum size on
-    disc, and forego direct RAM access for more efficient processing.
-    For example, in your code, run
+        :footcite:`BarnettSeth2015,WinklerEtAl2016` given by the same equation as for
+        ``'gc'``, but where the autocovariance sequence from which the autoregressive
+        model is produced is transposed to mimic the reversal of the original signal in
+        time :footcite:`HaufeEtAl2012`.
+
+    Parallel computation can be activated by setting the ``n_jobs`` parameter. Under the
+    hood, this utilizes the ``joblib`` library. For effective parallelization, you
+    should activate memory mapping in MNE-Python by setting ``MNE_MEMMAP_MIN_SIZE`` and
+    ``MNE_CACHE_DIR``. Activating memory mapping will make ``joblib`` store arrays
+    greater than the minimum size on disc, and forego direct RAM access for more
+    efficient processing. For example, in your code, run::
 
         mne.set_config('MNE_MEMMAP_MIN_SIZE', '10M')
         mne.set_config('MNE_CACHE_DIR', '/dev/shm')
 
-    When ``MNE_MEMMAP_MIN_SIZE=None``, the underlying joblib implementation
-    results in pickling and unpickling the whole array each time a pair of
-    indices is accessed, which is slow, compared to memory mapping the array.
+    When ``MNE_MEMMAP_MIN_SIZE=None``, the underlying ``joblib`` implementation results
+    in pickling and unpickling the whole array each time a pair of indices is accessed,
+    which is slow, compared to memory mapping the array.
 
-    This function is based on the ``frites.conn.conn_spec`` implementation in
-    Frites.
+    This function is based on the ``frites.conn.conn_spec`` implementation in Frites.
 
     .. versionadded:: 0.3
 
@@ -807,30 +794,27 @@ def _spectral_connectivity(
     sfreq : float
         Sampling frequency.
     freqs : array_like
-        Array of frequencies of interest for time-frequency decomposition.
-        Only the frequencies within the range specified by ``fmin`` and
-        ``fmax`` are used.
+        Array of frequencies of interest for time-frequency decomposition. Only the
+        frequencies within the range specified by ``fmin`` and ``fmax`` are used.
     faverage : bool
         Average over frequency bands.
     n_cycles : float | array_like of float
-        Number of cycles in the wavelet, either a fixed number or one per
-        frequency.
+        Number of cycles in the wavelet, either a fixed number or one per frequency.
     mt_bandwidth : float | None
         Multitaper time-bandwidth.
     gc_n_lags : int
-        Number of lags to use for the vector autoregressive model when
-        computing Granger causality.
+        Number of lags to use for the vector autoregressive model when computing Granger
+        causality.
     rank : tuple of array
         Ranks to project the seed and target data to.
     n_components : int
         Number of connectivity components to extract from the data. If 0, only the first
         component is extracted.
     decim : int
-        Decimation factor after time-frequency
-        decomposition.
+        Decimation factor after time-frequency decomposition.
     padding : float
-        Amount of time to consider as padding at the beginning and end of each
-        epoch in seconds.
+        Amount of time to consider as padding at the beginning and end of each epoch in
+        seconds.
     weights : array, shape (n_tapers, n_freqs) | None
         Taper weights for multitaper spectral estimation.
     multivariate_con : bool
@@ -842,15 +826,15 @@ def _spectral_connectivity(
     -------
     scores : dict
         Dictionary containing the connectivity estimates corresponding to the metrics in
-        ``method``. Each element is an array of shape (n_cons, [n_comps], n_freqs) or
-        (n_cons, [n_comps], n_fbands) if ``faverage`` is `True`. ``n_comps`` is present
-        for valid multivariate methods if ``n_components > 0``.
+        ``method``. Each element is an array of shape (n_cons, [n_comps,] n_freqs) or
+        (n_cons, [n_comps,] n_fbands) if ``faverage`` is ``True``. ``n_comps`` is
+        present for valid multivariate methods if ``n_components > 0``.
     patterns : dict
         Dictionary containing the connectivity patterns (for reconstructing the
         connectivity components in channel-space) corresponding to the metrics in
         ``method``, if multivariate methods are called, else an empty dictionary. Each
-        element is an array of shape (2, [n_comps], n_channels, n_freqs) or (2,
-        n_channels, n_fbands) if ``faverage`` is `True`, where 2 corresponds to the seed
+        element is an array of shape ``(2, [n_comps,] n_channels, n_freqs)`` or ``(2,
+        n_channels, n_fbands)`` if ``faverage=True``, where 2 corresponds to the seed
         and target signals (respectively). ``n_comps`` is present for valid multivariate
         methods if ``n_components > 0``.
     """
@@ -999,8 +983,8 @@ def _parallel_con(
     signals_use : list of int
         The unique signals on which connectivity is to be computed.
     gc_n_lags : int
-        Number of lags to use for the vector autoregressive model when
-        computing Granger causality.
+        Number of lags to use for the vector autoregressive model when computing Granger
+        causality.
     rank : tuple of array of int
         Ranks to project the seed and target data to.
     n_components : int
@@ -1020,11 +1004,11 @@ def _parallel_con(
     Returns
     -------
     out : tuple of list of array
-        Connectivity estimates for each signal pair, method, and frequency or
-        frequency band. If bivariate methods are called, the output is a tuple
-        of a list of arrays containing the connectivity scores. If multivariate
-        methods are called, the output is a tuple of lists containing arrays
-        for the connectivity scores and patterns, respectively.
+        Connectivity estimates for each signal pair, method, and frequency or frequency
+        band. If bivariate methods are called, the output is a tuple of a list of arrays
+        containing the connectivity scores. If multivariate methods are called, the
+        output is a tuple of lists containing arrays for the connectivity scores and
+        patterns, respectively.
     """
     if "coh" in method:
         # psd
@@ -1105,9 +1089,9 @@ def _pairwise_con(w, psd, x, y, method, kernel, foi_idx, faverage, weights):
     Returns
     -------
     out : list
-        List of connectivity estimates between signals ``x`` and ``y``
-        corresponding to the methods in ``method``. Each element is an array
-        with shape (n_freqs,) or (n_fbands) depending on ``faverage``.
+        List of connectivity estimates between signals ``x`` and ``y`` corresponding to
+        the methods in ``method``. Each element is an array with shape ``(n_freqs,)`` or
+        ``(n_fbands)`` depending on ``faverage``.
     """
     w_x, w_y = w[x], w[y]
     if weights is not None:
@@ -1159,12 +1143,11 @@ def _multivariate_con(
     w : array_like, shape (n_chans, n_tapers, n_freqs, n_times)
         Time-frequency data.
     seeds : array, shape of (n_cons, n_channels)
-        Seed channel indices. ``n_channels`` is the largest number of channels
-        across all connections, with missing entries padded with ``-1``.
+        Seed channel indices. ``n_channels`` is the largest number of channels across
+        all connections, with missing entries padded with ``-1``.
     targets : array, shape of (n_cons, n_channels)
-        Target channel indices. ``n_channels`` is the largest number of
-        channels across all connections, with missing entries padded with
-        ``-1``.
+        Target channel indices. ``n_channels`` is the largest number of channels across
+        all connections, with missing entries padded with ``-1``.
     signals_use : list of int
         The unique signals on which connectivity is to be computed.
     method : str
@@ -1178,8 +1161,8 @@ def _multivariate_con(
     weights : array_like, shape (n_tapers, n_freqs, n_times) | None
         Multitaper weights.
     gc_n_lags : int
-        Number of lags to use for the vector autoregressive model when
-        computing Granger causality.
+        Number of lags to use for the vector autoregressive model when computing Granger
+        causality.
     rank : tuple of array, shape of (2, n_cons)
         Ranks to project the seed and target data to.
     n_components : int
@@ -1192,17 +1175,17 @@ def _multivariate_con(
     -------
     scores : list
         List of connectivity scores between seed and target signals for each
-        connectivity method. Each element is an array with shape ([n_comps], n_freqs) or
-        ([n_comps], n_fbands) depending on ``faverage``. ``n_comps`` is present for
-        valid multivariate methods if ``n_components > 0``.
+        connectivity method. Each element is an array with shape
+        ``([n_comps,] n_freqs)`` or ``([n_comps,] n_fbands)`` depending on ``faverage``.
+        ``n_comps`` is present for valid multivariate methods if ``n_components > 0``.
     patterns : list
         List of connectivity patterns between seed and target signals for each
         connectivity method. Each element is an array of length 2 corresponding to the
-        seed and target patterns, respectively, each with shape ([n_comps], n_channels,
-        n_freqs) or ([n_comps], n_channels, n_fbands) depending on ``faverage``.
-        ``n_comps`` is present for valid multivariate methods if ``n_components > 0``.
-        ``n_channels`` is the largest number of channels across all connections, with
-        missing entries padded with ``np.nan``.
+        seed and target patterns, respectively, each with shape ``([n_comps,]
+        n_channels, n_freqs)`` or ``([n_comps], n_channels, n_fbands) depending on
+        ``faverage``. ``n_comps`` is present for valid multivariate methods if
+        ``n_components > 0``. ``n_channels`` is the largest number of channels across
+        all connections, with missing entries padded with ``np.nan``.
     """
     csd = []
     for x in signals_use:
@@ -1265,13 +1248,13 @@ def _plv(s_xy):
 
     Parameters
     ----------
-    s_xy : array-like, shape (n_freqs, n_times)
-        The cross PSD between channel 'x' and channel 'y' across
-        frequency and time points.
+    s_xy : array, shape (n_freqs, n_times)
+        The cross PSD between channel 'x' and channel 'y' across frequency and time
+        points.
 
     Returns
     -------
-    plv : array-like, shape (n_freqs, n_times)
+    plv : array, shape (n_freqs, n_times)
         The estimated PLV.
     """
     s_xy = s_xy / np.abs(s_xy)
@@ -1284,13 +1267,13 @@ def _ciplv(s_xy):
 
     Parameters
     ----------
-    s_xy : array-like, shape (n_freqs, n_times)
-        The cross PSD between channel 'x' and channel 'y' across
-        frequency and time points.
+    s_xy : array, shape (n_freqs, n_times)
+        The cross PSD between channel 'x' and channel 'y' across frequency and time
+        points.
 
     Returns
     -------
-    ciplv : array-like, shape (n_freqs, n_times)
+    ciplv : array, shape (n_freqs, n_times)
         The estimated ciPLV.
     """
     s_xy = s_xy / np.abs(s_xy)
@@ -1305,13 +1288,13 @@ def _pli(s_xy):
 
     Parameters
     ----------
-    s_xy : array-like, shape (n_freqs, n_times)
-        The cross PSD between channel 'x' and channel 'y' across
-        frequency and time points.
+    s_xy : array, shape (n_freqs, n_times)
+        The cross PSD between channel 'x' and channel 'y' across frequency and time
+        points.
 
     Returns
     -------
-    pli : array-like, shape (n_freqs, n_times)
+    pli : array, shape (n_freqs, n_times)
         The estimated PLI.
     """
     pli = np.abs(np.mean(np.sign(np.imag(s_xy)), axis=-1, keepdims=True))
@@ -1323,13 +1306,13 @@ def _wpli(s_xy):
 
     Parameters
     ----------
-    s_xy : array-like, shape (n_freqs, n_times)
-        The cross PSD between channel 'x' and channel 'y' across
-        frequency and time points.
+    s_xy : array, shape (n_freqs, n_times)
+        The cross PSD between channel 'x' and channel 'y' across frequency and time
+        points.
 
     Returns
     -------
-    wpli : array-like, shape (n_freqs, n_times)
+    wpli : array, shape (n_freqs, n_times)
         The estimated wPLI.
     """
     con_num = np.abs(s_xy.imag.mean(axis=-1, keepdims=True))
@@ -1343,17 +1326,17 @@ def _coh(s_xx, s_yy, s_xy):
 
     Parameters
     ----------
-    s_xx : array-like, shape (n_freqs, n_times)
+    s_xx : array, shape (n_freqs, n_times)
         The PSD of channel 'x'.
-    s_yy : array-like, shape (n_freqs, n_times)
+    s_yy : array, shape (n_freqs, n_times)
         The PSD of channel 'y'.
-    s_xy : array-like, shape (n_freqs, n_times)
-        The cross PSD between channel 'x' and channel 'y' across
-        frequency and time points.
+    s_xy : array, shape (n_freqs, n_times)
+        The cross PSD between channel 'x' and channel 'y' across frequency and time
+        points.
 
     Returns
     -------
-    coh : array-like, shape (n_freqs, n_times)
+    coh : array, shape (n_freqs, n_times)
         The estimated COH.
     """
     con_num = np.abs(s_xy.mean(axis=-1, keepdims=True))
@@ -1382,14 +1365,14 @@ def _foi_average(conn, foi_idx):
 
     Parameters
     ----------
-    conn : array_like, shape (..., n_freqs, n_times)
+    conn : array, shape (..., n_freqs, n_times)
         Connectivity estimate array.
-    foi_idx : array_like, shape (n_foi, 2)
+    foi_idx : array, shape (n_foi, 2)
         Upper and lower frequency bounds of each frequency band.
 
     Returns
     -------
-    conn_f : np.ndarray, shape (..., n_fbands, n_times)
+    conn_f : array, shape (..., n_fbands, n_times)
         Connectivity estimate array, averaged within frequency bands.
     """
     # get the number of foi
diff --git a/mne_connectivity/utils/docs.py b/mne_connectivity/utils/docs.py
index e8f62ab1c..47993f978 100644
--- a/mne_connectivity/utils/docs.py
+++ b/mne_connectivity/utils/docs.py
@@ -17,35 +17,33 @@
 
 # Connectivity
 docdict["data"] = """
-data : np.ndarray ([epochs], n_estimated_nodes, [components], [freqs], [times])
-    The connectivity data that is a raveled array of ``(n_estimated_nodes, ...)`` shape.
-    The ``n_estimated_nodes`` is equal to ``n_nodes_in * n_nodes_out`` if one is
-    computing the full connectivity, or a subset of nodes equal to the length of
-    ``indices`` passed in.
+data : array, shape ([epochs,] n_estimated_nodes[, components, freqs, times])
+    The connectivity data that is a raveled array of ``(..., n_estimated_nodes, ...)``
+    shape. The ``n_estimated_nodes`` is equal to ``n_nodes_in * n_nodes_out`` if one has
+    computed the full connectivity, or a subset of nodes equal to the length of the
+    arrays in ``indices`` passed in.
 """
 
 docdict["names"] = """
-names : list | np.ndarray | None
-    The names of the nodes of the dataset used to compute
-    connectivity. If 'None' (default), then names will be
-    a list of integers from 0 to ``n_nodes``. If a list
+names : array_like | None
+    The names of the nodes of the dataset used to compute connectivity. If ``None``
+    (default), then names will be a list of integers from 0 to ``n_nodes``. If a list
     of names, then it must be equal in length to ``n_nodes``.
 """
 
 docdict["indices"] = """
-indices : tuple of arrays | str | None
-    The indices of relevant connectivity data. If ``'all'`` (default),
-    then data is connectivity between all nodes. If ``'symmetric'``,
-    then data is symmetric connectivity between all nodes. If a tuple,
-    then the first list represents the "in nodes", and the second list
-    represents the "out nodes". See "Notes" for more information.
+indices : tuple of array_like | ``'all'`` | ``'symmetric'`` | None
+    The indices of relevant connectivity data. If ``'all'`` (default), then data is
+    connectivity between all nodes. If ``'symmetric'``, then data is symmetric
+    connectivity between all nodes. If a tuple, then contains two array-likes where the
+    first array represents the "in nodes" (seeds), and the second array represents the
+    "out nodes" (targets).
 """
 
 docdict["n_nodes"] = """
 n_nodes : int
-    The number of nodes in the dataset used to compute connectivity.
-    This should be equal to the number of signals in the original
-    dataset.
+    The number of nodes in the dataset used to compute connectivity. This should be
+    equal to the number of signals in the original dataset.
 """
 
 docdict["connectivity_kwargs"] = """
@@ -57,13 +55,12 @@
 """
 
 docdict["mode"] = """
-mode : str (default 'multitaper')
-    The cross-spectral density computation method. Can be ``'multitaper'``,
-    ``'fourier'``, or ``'cwt_morlet'``.
+mode : ``'multitaper'`` | ``'fourier'`` | ``'cwt_morlet'``
+    The cross-spectral density computation method.
 """
 
 docdict["mt_bandwidth"] = """
-mt_bandwidth : int | float | None (default None)
+mt_bandwidth : float | None (default None)
     The bandwidth of the multitaper windowing function in Hz to use when computing the
     cross-spectral density. Only used if ``mode='multitaper'``.
 """
@@ -82,13 +79,13 @@
 """
 
 docdict["cwt_freqs"] = """
-cwt_freqs : array of int or float | None (default None)
+cwt_freqs : array_like | None (default None)
     The frequencies of interest in Hz. Must not be `None` and only used if
     ``mode='cwt_morlet'``.
 """
 
 docdict["cwt_n_cycles"] = """
-cwt_n_cycles : int | float | array of int or float (default 7)
+cwt_n_cycles : float | array_like (default 7.0)
     The number of cycles to use when constructing the Morlet wavelets. Fixed number or
     one per frequency. Only used if ``mode='cwt_morlet'``.
 """
@@ -112,73 +109,65 @@
 
 # Downstream container variables
 docdict["freqs"] = """
-freqs : list | np.ndarray
-    The frequencies at which the connectivity data is computed over.
-    If the frequencies are "frequency bands" (i.e. gamma band), then
-    these are the median of those bands.
+freqs : list | array
+    The frequencies at which the connectivity data is computed over. If the frequencies
+    are "frequency bands" (i.e. gamma band), then these are the median of those bands.
 """
 
 docdict["times"] = """
-times : list | np.ndarray
+times : list | array
     The times at which the connectivity data is computed over.
 """
 
 docdict["method"] = """
-method : str, optional
-    The method name used to compute connectivity.
+method : str | None
+    The method name used to compute connectivity (default ``None``).
 """
 
 docdict["spec_method"] = """
-spec_method : str, optional
-    The type of method used to compute spectral analysis,
-    by default None.
+spec_method : str | None
+    The type of method used to compute spectral analysis (default ``None``).
 """
 
 docdict["n_epochs_used"] = """
-n_epochs_used : int, optional
-    The number of epochs used in the computation of connectivity,
-    by default None.
+n_epochs_used : int | None
+    The number of epochs used in the computation of connectivity (default ``None``).
 """
 
 docdict["events"] = """
 events : array of int, shape (n_events, 3)
-    The events typically returned by the read_events function.
-    If some events don't match the events of interest as specified
-    by event_id, they will be marked as 'IGNORED' in the drop log.
+    The events typically returned by the read_events function. If some events don't
+    match the events of interest as specified by ``event_id``, they will be marked as
+    'IGNORED' in the drop log.
 """
 
 docdict["event_id"] = """
 event_id : int | list of int | dict | None
-    The id of the event to consider. If dict,
-    the keys can later be used to access associated events. Example:
-    dict(auditory=1, visual=3). If int, a dict will be created with
-    the id as string. If a list, all events with the IDs specified
-    in the list are used. If None, all events will be used with
-    and a dict is created with string integer names corresponding
-    to the event id integers.
+    The ID of the event to consider. If a dict, the keys can later be used to access
+    associated events. Example: ``dict(auditory=1, visual=3)``. If an int, a dict will
+    be created with the ID as string. If a list, all events with the IDs specified in
+    the list are used. If ``None``, all events will be used and a dict is created with
+    string integer names corresponding to the event ID integers.
 """
 
 # Verbose
 docdict["verbose"] = """
-verbose : bool, str, int, or None
-    If not None, override default verbose level (see :func:`mne.verbose`
-    for more info). If used, it should be passed as a
-    keyword-argument only."""
+verbose : bool | str | int | None
+    If not ``None``, override default verbose level (see :func:`mne.verbose` for more
+    info). If used, it should be passed as a keyword-argument only."""
 
 # Parallelization
 docdict["n_jobs"] = """
 n_jobs : int
-    The number of jobs to run in parallel (default 1).
-    Requires the joblib package.
+    The number of jobs to run in parallel (default 1). Requires the joblib package.
 """
 
 # Random state
 docdict["random_state"] = """
 random_state : None | int | instance of ~numpy.random.RandomState
-    If ``random_state`` is an :class:`int`, it will be used as a seed for
-    :class:`~numpy.random.RandomState`. If ``None``, the seed will be
-    obtained from the operating system (see
-    :class:`~numpy.random.RandomState` for details). Default is
+    If ``random_state`` is an int, it will be used as a seed for
+    :class:`numpy.random.RandomState`. If ``None``, the seed will be obtained from the
+    operating system (see :class:`numpy.random.RandomState` for details). Default is
     ``None``.
 """
 
@@ -191,7 +180,7 @@
 """
 
 docdict["method_decoding"] = """
-method : str
+method : ``'cacoh'`` | ``'mic'``
     The multivariate method to use for the decomposition. Can be:
 
     * ``'cacoh'`` - Canonical Coherency (CaCoh) :footcite:`VidaurreEtAl2019`
@@ -199,38 +188,40 @@
 """
 
 docdict["fmin_decoding"] = """
-fmin : int | float | None (default None)
-    The lowest frequency of interest in Hz. Must not be `None` and only used if
+fmin : float | None (default None)
+    The lowest frequency of interest in Hz. Must not be ``None`` and only used if
     ``mode in ['multitaper', 'fourier']``.
 """
 
 docdict["fmax_decoding"] = """
-fmax : int | float | None (default None)
-    The highest frequency of interest in Hz. Must not be `None` and only used if
+fmax : float | None (default None)
+    The highest frequency of interest in Hz. Must not be ``None`` and only used if
     ``mode in ['multitaper', 'fourier']``.
 """
 
 docdict["indices_decoding"] = """
-indices : tuple of array
-    A tuple of two arrays, containing the indices of the seed and target channels in the
-    input data, respectively. The indices of only a single connection (i.e. between one
-    group of seeds and one group of targets) is supported.
+indices : tuple of array_like
+    A tuple of two array-likes, containing the indices of the seed and target channels
+    in the input data, respectively. The indices of only a single connection (i.e.
+    between one group of seeds and one group of targets) is supported, such that the
+    array-likes have shape ``(n_seed_channels,)`` and ``(n_target_channels,)``.
 """
 
 docdict["n_components"] = """
 n_components : int | None (default None)
-    The number of connectivity components (sources) to extract from the data. If `None`,
-    the number of components equal to the minimum rank of the seeds and targets is
-    extracted (see the ``rank`` parameter). If an `int`, the number of components must
-    be <= the minimum rank of the seeds and targets. E.g. if the seed channels had a
-    rank of 5 and the target channels had a rank of 3, ``n_components`` must be <= 3.
+    The number of connectivity components (sources) to extract from the data. If
+    ``None``, the number of components equal to the minimum rank of the seeds and
+    targets is extracted (see the ``rank`` parameter). If an int, the number of
+    components must be <= the minimum rank of the seeds and targets. E.g. if the seed
+    channels had a ank of 5 and the target channels had a rank of 3, ``n_components``
+    must be <= 3.
 """
 
 docdict["rank"] = """
 rank : tuple of int | None (default None)
     A tuple of two ints, containing the degree of rank subspace projection to apply to
-    the seed and target data, respectively, before filters are fit. If `None`, the rank
-    of the seed and target data is used. If a tuple of ints, the entries must be <= the
+    the seed and target data, respectively, before filters are fit. If ``None``, the
+    rank of the seed and target data is used. If a tuple, the entries must be <= the
     rank of the seed and target data. The minimum rank of the seeds and targets
     determines the maximum number of connectivity components (sources) which can be
     extracted from the data (see the ``n_components`` parameter). Specifying ranks below
@@ -239,15 +230,17 @@
 
 # Decoding attrs
 docdict["filters_"] = """
-filters_ : tuple of array, shape=(n_signals, n_components)
+filters_ : tuple of length 2
     A tuple of two arrays containing the spatial filters for transforming the seed and
-    target data, respectively.
+    target data, respectively. The arrays have shape ``(n_seeds, n_components)`` and
+    ``(n_targets, n_components)``.
 """
 
 docdict["patterns_"] = """
-patterns_ : tuple of array, shape=(n_components, n_signals)
+patterns_ : tuple of length 2
     A tuple of two arrays containing the spatial patterns corresponding to the spatial
-    filters for the seed and target data, respectively.
+    filters for the seed and target data, respectively. The arrays have shape
+    ``(n_components, n_seeds)`` and ``(n_components, n_targets)``.
 """
 
 # Decoding plotting
@@ -259,55 +252,55 @@
 
 # Topomaps
 docdict["components_topomap"] = """
-components : int | array of int | None (default None)
-    The components to plot. If `None`, all components are shown.
+components : int | array_like of int | None (default None)
+    The components to plot. If ``None``, all components are shown.
 """
 
 docdict["ch_type_topomap"] = """
-ch_type : 'mag' | 'grad' | 'planar1' | 'planar2' | 'eeg' | None (default None)
+ch_type : ``'mag'`` | ``'grad'`` | ``'planar1'`` | ``'planar2'`` | ``'eeg'`` | None (default None)
     The channel type to plot. For ``'grad'``, the gradiometers are collected in pairs
-    and the RMS for each pair is plotted. If `None`, the first available channel type
+    and the RMS for each pair is plotted. If ``None``, the first available channel type
     from the order shown above is used.
-"""
+"""  # noqa E501
 
 docdict["scalings_topomap"] = """
 scalings : dict | float | None (default None)
-    The scalings of the channel types to be applied for plotting. If `None`, uses
+    The scalings of the channel types to be applied for plotting. If ``None``, uses
     ``dict(eeg=1e6, grad=1e13, mag=1e15)``.
 """
 
 docdict["sensors_topomap"] = """
 sensors : bool | str (default True)
-    Whether to add markers for sensor locations. If `str`, should be a valid
-    matplotlib format string (e.g., ``'r+'`` for red plusses; see the Notes section of
-    :meth:`~matplotlib.axes.Axes.plot`). If `True`, black circles are used.
+    Whether to add markers for sensor locations. If a str, should be a valid matplotlib
+    format string (e.g., ``'r+'`` for red plusses; see the Notes section of
+    :meth:`matplotlib.axes.Axes.plot`). If ``True``, black circles are used.
 """
 
 docdict["show_names_topomap"] = """
 show_names : bool | callable (default False)
-    Whether to show channel names next to each sensor marker. If `callable`, channel
+    Whether to show channel names next to each sensor marker. If a callable, channel
     names will be formatted using the callable; e.g., to delete the prefix 'MEG ' from
     all channel names, pass the function ``lambda x: x.replace('MEG ', '')``. If
-    ``mask`` is not `None`, only non-masked sensor names will be shown.
+    ``mask`` is not ``None``, only non-masked sensor names will be shown.
 """
 
 docdict["mask_filters_topomap"] = """
-mask : array of bool, shape=(n_channels, n_filters) | None (default None)
+mask : array of bool, shape (n_channels, n_filters) | None (default None)
     An array specifying channel-filter combinations to highlight with a distinct
-    plotting style. Array elements set to `True` will be plotted with the parameters
-    given in ``mask_params``. If `None`, no combinations will be highlighted.
+    plotting style. Array elements set to ``True`` will be plotted with the parameters
+    given in ``mask_params``. If ``None``, no combinations will be highlighted.
 """
 docdict["mask_patterns_topomap"] = """
-mask : array of bool, shape=(n_channels, n_patterns) | None (default None)
+mask : array of bool, shape (n_channels, n_patterns) | None (default None)
     An array specifying channel-pattern combinations to highlight with a distinct
-    plotting style. Array elements set to `True` will be plotted with the parameters
-    given in ``mask_params``. If `None`, no combinations will be highlighted.
+    plotting style. Array elements set to ``True`` will be plotted with the parameters
+    given in ``mask_params``. If ``None``, no combinations will be highlighted.
 """
 
 docdict["mask_params_topomap"] = """
 mask_params : dict | None (default None)
     The plotting parameters for distinct combinations given in ``mask``.
-    Default `None` equals::
+    Default ``None`` equals::
 
         dict(marker='o', markerfacecolor='w', markeredgecolor='k',
              linewidth=0, markersize=4)
@@ -324,29 +317,29 @@
 """
 
 docdict["outlines_topomap"] = """
-outlines : 'head' | dict | None (default 'head')
-    The outlines to be drawn. If 'head', the default head scheme will be drawn. If dict,
-    each key refers to a tuple of x and y positions, the values in 'mask_pos' will serve
-    as image mask. Alternatively, a matplotlib patch object can be passed for advanced
-    masking options, either directly or as a function that returns patches (required for
-    multi-axis plots). If `None`, nothing will be drawn.
+outlines : ``'head'`` | dict | None (default ``'head'``)
+    The outlines to be drawn. If ``'head'``, the default head scheme will be drawn. If a
+    dict, each key refers to a tuple of x and y positions, the values in ``'mask_pos'``
+    will serve as an image mask. Alternatively, a matplotlib patch object can be passed
+    for advanced masking options, either directly or as a function that returns patches
+    (required for multi-axis plots). If ``None``, nothing will be drawn.
 """
 
 docdict["sphere_topomap"] = """
-sphere : float | array | mne.bem.ConductorModel | None  | 'auto' | 'eeglab' (default None)
-    The sphere parameters to use for the head outline. Can be array-like of shape (4,)
-    to give the X/Y/Z origin and radius in meters, or a single float to give just the
-    radius (origin assumed 0, 0, 0). Can also be an instance of a spherical
-    :class:`~mne.bem.ConductorModel` to use the origin and radius from that object. If
+sphere : float | array | mne.bem.ConductorModel | None  | ``'auto'`` | ``'eeglab'`` (default None)
+    The sphere parameters to use for the head outline. Can be an array-like of shape
+    (4,) to give the X/Y/Z origin and radius in meters, or a single float to give just
+    the radius (origin assumed 0, 0, 0). Can also be an instance of a spherical
+    :class:`mne.bem.ConductorModel` to use the origin and radius from that object. If
     ``'auto'`` the sphere is fit to digitization points. If ``'eeglab'`` the head circle
     is defined by EEG electrodes ``'Fpz'``, ``'Oz'``, ``'T7'``, and ``'T8'`` (if
     ``'Fpz'`` is not present, it will be approximated from the coordinates of ``'Oz'``).
-    `None` is equivalent to ``'auto'`` when enough extra digitization points are
+    ``None`` is equivalent to ``'auto'`` when enough extra digitization points are
     available, and (0, 0, 0, 0.95) otherwise.
 """  # noqa E501
 
 docdict["image_interp_topomap"] = """
-image_interp : str (default 'cubic')
+image_interp : ``'cubic'`` | ``'nearest'`` | ``'linear'`` (default ``'cubic'``)
     The image interpolation to be used. Options are ``'cubic'`` to use
     :class:`scipy.interpolate.CloughTocher2DInterpolator`, ``'nearest'`` to use
     :class:`scipy.spatial.Voronoi`, or ``'linear'`` to use
@@ -354,7 +347,7 @@
 """
 
 docdict["extrapolate_topomap"] = """
-extrapolate : str
+extrapolate : ``'box'`` | ``'local'`` | ``'head'``
     The extrapolation options. Can be one of:
 
     - ``'box'``
@@ -372,7 +365,7 @@
 """
 
 docdict["border_topomap"] = """
-border : float | 'mean' (default 'mean')
+border : float | ``'mean'`` (default ``'mean'``)
     The value to extrapolate to on the topomap borders. If ``'mean'``, each extrapolated
     point has the average value of its neighbours.
 """
@@ -383,21 +376,21 @@
 """
 
 docdict["size_topomap"] = """
-size : int | float (default 1)
+size : float (default 1.0)
     The side length of each subplot in inches.
 """
 
 docdict["cmap_topomap"] = """
-cmap : str | matplotlib.colors.Colormap | (matplotlib.colors.Colormap, bool) | 'interactive' | None (default 'RdBu_r')
-    The colormap to use. If a `str`, should be a valid matplotlib colormap. If a
-    `tuple`, the first value is `matplotlib.colors.Colormap` object to use and the
-    second value is a boolean defining interactivity. In interactive mode the colors are
-    adjustable by clicking and dragging the colorbar with left and right mouse button.
-    Left mouse button moves the scale up and down and right mouse button adjusts the
-    range. Hitting space bar resets the range. Up and down arrows can be used to change
-    the colormap. If `None`, ``'Reds'`` is used for data that is either all positive or
-    all negative, and ``'RdBu_r'`` is used otherwise. ``'interactive'`` is equivalent to
-    ``(None, True)``.
+cmap : str | matplotlib.colors.Colormap | (matplotlib.colors.Colormap, bool) | ``'interactive'`` | None (default ``'RdBu_r'``)
+    The colormap to use. If a str, should be a valid matplotlib colormap. If a tuple,
+    the first value is :class:`matplotlib.colors.Colormap` object to use and the
+    second value is a boolean defining interactivity. In interactive mode, the colors
+    are adjustable by clicking and dragging the colorbar with left and right mouse
+    button. Left mouse button moves the scale up and down and right mouse button adjusts
+    the range. Hitting space bar resets the range. Up and down arrows can be used to
+    change the colormap. If ``None``, ``'Reds'`` is used for data that is either all
+    positive or all negative, and ``'RdBu_r'`` is used otherwise. ``'interactive'`` is
+    equivalent to ``(None, True)``.
 
     .. warning::  Interactive mode works smoothly only for a small amount
         of topomaps. Interactive mode is disabled by default for more than
@@ -405,16 +398,16 @@
 """  # noqa E501
 
 docdict["vlim_topomap"] = """
-vlim : tuple of length 2 (default (None, None))
-    The lower and upper colormap bounds, respectively. If both entries are `None`, sets
-    bounds to ``(min(data), max(data))``. If one entry is `None`, the corresponding
-    boundary is set at the min/max of the data.
+vlim : tuple of float or None (default (None, None))
+    The lower and upper colormap bounds, respectively. If both entries are ``None``,
+    sets bounds to ``(min(data), max(data))``. If one entry is ``None``, the
+    corresponding boundary is set at the min/max of the data.
 """
 
 docdict["cnorm_topomap"] = """
 cnorm : matplotlib.colors.Normalize | None (default None)
-    How to normalize the colormap. If `None`, standard linear normalization is used. If
-    not `None`, ``vlim`` is ignored. See the :ref:`Matplotlib docs
+    How to normalize the colormap. If ``None``, standard linear normalization is used.
+    If not ``None``, ``vlim`` is ignored. See the :ref:`Matplotlib docs
     <matplotlib:colormapnorms>` for more details on colormap normalization.
 """
 
@@ -424,39 +417,40 @@
 """
 
 docdict["colorbar_format_topomap"] = r"""
-cbar_fmt : str (default '%.1E')
+cbar_fmt : str (default ``'%.1E'``)
     The formatting string for colorbar tick labels. See :ref:`formatspec` for details.
 """
 
 docdict["units_topomap"] = """
-units : str (default 'AU')
+units : str (default ``'AU'``)
     The units for the colorbar label. Ignored if ``colorbar=False``.
 """
 
 docdict["axes_topomap"] = """
-axes : length-2 tuple of list of matplotlib.axes.Axes | None (default None)
-    The axes to plot to. If `None`, a new figure will be created with the correct number
-    of axes. If not `None`, there must be two lists containing the axes for the seeds
-    and targets, respectively. In each of these two lists, the number of axes must match
-    ``components`` if ``colorbar=False``, or ``components * 2`` if ``colorbar=True``.
+axes : tuple of length 2 | None (default None)
+    The axes to plot to. If ``None``, a new figure will be created with the correct
+    number of axes. If not ``None``, there must be two lists containing the axes for the
+    seeds and targets, respectively. In each of these two lists, the number of axes must
+    match ``components`` if ``colorbar=False``, or ``components * 2`` if
+    ``colorbar=True``.
 """
 
 docdict["name_format_topomap"] = r"""
 name_format : str | None (default None)
-    The string format for axes titles. If `None`, uses ``f"{method}%01d_{group}"``,
+    The string format for axes titles. If ``None``, uses ``f"{method}%01d ({group})"``,
     i.e., the method name followed by the component number and the group being plotted
-    (seeds or targets). If not `None`, it must contain a formatting specifier for the
+    (seeds or targets). If not ``None``, it must contain a formatting specifier for the
     component number, and the group will be appended to the end.
 """
 
 docdict["nrows_topomap"] = """
-nrows : int | 'auto' (default 'auto')
+nrows : int | ``'auto'`` (default ``'auto'``)
     The number of rows of components to plot. If ``'auto'``, the necessary number will
     be inferred.
 """
 
 docdict["ncols_topomap"] = """
-ncols : int | 'auto' (default 'auto')
+ncols : int | ``'auto'`` (default ``'auto'``)
     The number of columns of components to plot. If ``'auto'``, the necessary number
     will be inferred. If ``nrows='auto'`` and ``ncols='auto'``, becomes ``nrows=1,
     ncols='auto'``.
diff --git a/mne_connectivity/utils/utils.py b/mne_connectivity/utils/utils.py
index 23ffdf628..875b523ab 100644
--- a/mne_connectivity/utils/utils.py
+++ b/mne_connectivity/utils/utils.py
@@ -14,14 +14,14 @@ def parallel_loop(func, n_jobs=1, verbose=1):
     func : function
         function to be executed in parallel
     n_jobs : int | None
-        Number of jobs. If set to None, do not attempt to use joblib.
+        Number of jobs. If set to ``None``, do not attempt to use joblib.
     verbose : int
         verbosity level
 
     Notes
     -----
-    Execution of the main script must be guarded with
-    `if __name__ == '__main__':` when using parallelization.
+    Execution of the main script must be guarded with ``if __name__ == '__main__':``
+    when using parallelization.
     """
     if n_jobs:
         try:
@@ -53,20 +53,21 @@ def check_indices(indices):
 
     Parameters
     ----------
-    indices : tuple of array of int, shape (2, n_cons)
-        Tuple containing index pairs.
+    indices : tuple of length 2 of array_like, shape (n_cons,)
+        A tuple of 2 arrays containing the seed and target channel indices,
+        respectively.
 
     Returns
     -------
-    indices : tuple of array of int, shape (2, n_cons)
-        The indices.
+    indices : tuple of length 2 of array_like, shape (n_cons,)
+        The indices to use for connectivity computation.
 
     Notes
     -----
-    Indices for bivariate connectivity should be a tuple of length 2,
-    containing the channel indices for the seed and target channel pairs,
-    respectively. Seed and target indices should be equal-length array-likes of
-    integers representing the indices of the individual channels in the data.
+    Indices for bivariate connectivity should be a tuple of length 2, containing the
+    channel indices for the seed and target channel pairs, respectively. Seed and target
+    indices should be equal-length array-likes of integers representing the indices of
+    the individual channels in the data.
     """
     if not isinstance(indices, tuple) or len(indices) != 2:
         raise ValueError("indices must be a tuple of length 2")
@@ -90,37 +91,37 @@ def _check_multivariate_indices(indices, n_chans):
 
     Parameters
     ----------
-    indices : tuple of array of array of int, shape (2, n_cons, variable)
-        Tuple containing index sets.
-
+    indices : tuple of length 2 of array_like
+        A tuple of two array-likes containing the seed and target indices, respectively,
+        to use for connectivity computation. Each array-like has ``n_cons`` array-like
+        entries containing the channel indices for each connection. The number of
+        channels within each connection can vary.
     n_chans : int
-        The number of channels in the data. Used when converting negative
-        indices to positive indices.
+        The number of channels in the data. Used when converting negative indices to
+        positive indices.
 
     Returns
     -------
-    indices : array of array of int, shape of (2, n_cons, max_n_chans)
+    indices : array, shape (2, n_cons, max_n_chans)
         The padded indices as a masked array.
 
     Notes
     -----
-    Indices for multivariate connectivity should be a tuple of length 2
-    containing the channel indices for the seed and target channel sets,
-    respectively. Seed and target indices should be equal-length array-likes
-    representing the indices of the channel sets in the data for each
-    connection. The indices for each connection should be an array-like of
-    integers representing the individual channels in the data. The length of
-    indices for each connection do not need to be equal. All indices within a
+    Indices for multivariate connectivity should be a tuple of length 2 containing the
+    channel indices for the seed and target channel sets, respectively. Seed and target
+    indices should be equal-length array-likes representing the indices of the channel
+    sets in the data for each connection. The indices for each connection should be an
+    array-like of integers representing the individual channels in the data. The length
+    of indices for each connection do not need to be equal. All indices within a
     connection must be unique.
 
-    If the seed and target indices are given as lists or tuples, they will be
-    converted to numpy arrays. Because the number of channels can differ across
-    connections or between the seeds and targets for a given connection (i.e.
-    ragged/jagged indices), the returned array will be padded out to a 'full'
-    array with an invalid index (``-1``) according to the maximum number of
-    channels in the seed or target of any one connection. These invalid
-    entries are then masked and returned as numpy masked arrays. E.g. the
-    ragged indices of shape ``(2, n_cons, variable)``::
+    If the seed and target indices are given as lists or tuples, they will be converted
+    to numpy arrays. Because the number of channels can differ across connections or
+    between the seeds and targets for a given connection (i.e. ragged/jagged indices),
+    the returned array will be padded out to a 'full' array with an invalid index
+    (``-1``) according to the maximum number of channels in the seed or target of any
+    one connection. These invalid entries are then masked and returned as numpy masked
+    arrays. E.g. the ragged indices of shape ``(2, n_cons, variable)``::
 
             indices = ([[0, 1], [0, 1   ]],  # seeds
                        [[2, 3], [4, 5, 6]])  # targets
@@ -130,19 +131,18 @@ def _check_multivariate_indices(indices, n_chans):
             indices = ([[0, 1, -1], [0, 1, -1]],  # seeds
                        [[2, 3, -1], [4, 5,  6]])  # targets
 
-    to have shape ``(2, n_cons, max_n_chans)``, where ``max_n_chans = 3``. The
-    invalid entries are then masked::
+    to have shape ``(2, n_cons, max_n_chans)``, where ``max_n_chans = 3``. The invalid
+    entries are then masked::
 
             indices = ([[0, 1, --], [0, 1, --]],  # seeds
                        [[2, 3, --], [4, 5,  6]])  # targets
 
     In case "indices" contains negative values to index channels, these will be
-    converted to the corresponding positive-valued index before any masking is
-    applied.
+    converted to the corresponding positive-valued index before any masking is applied.
 
-    More information on working with multivariate indices and handling
-    connections where the number of seeds and targets are not equal can be
-    found in the :doc:`../auto_examples/handling_ragged_arrays` example.
+    More information on working with multivariate indices and handling connections where
+    the number of seeds and targets are not equal can be found in the
+    :doc:`../auto_examples/handling_ragged_arrays` example.
     """
     if not isinstance(indices, tuple) or len(indices) != 2:
         raise ValueError("indices must be a tuple of length 2")
@@ -196,22 +196,23 @@ def seed_target_indices(seeds, targets):
 
     Parameters
     ----------
-    seeds : array of int | int, shape (n_unique_seeds)
-        Seed indices.
-    targets : array of int | int, shape (n_unique_targets)
-        Indices of signals for which to compute connectivity.
+    seeds : array_like, shape (n_unique_seeds) | int
+        Indices of signals for which to compute connectivity from.
+    targets : array_like, shape (n_unique_targets) | int
+        Indices of signals for which to compute connectivity to.
 
     Returns
     -------
-    indices : tuple of array of int, shape (2, n_cons)
-        The indices parameter used for connectivity computation.
+    indices : tuple of length 2 of array, shape (n_cons,)
+        A tuple of 2 arrays containing the seed and target indices, respectively, to use
+        for connectivity computation.
 
     Notes
     -----
-    ``seeds`` and ``targets`` should be array-likes or integers representing
-    the indices of the channel pairs in the data for each connection. ``seeds``
-    and ``targets`` will be expanded such that connectivity will be computed
-    between each seed and each target. E.g. the seeds and targets::
+    ``seeds`` and ``targets`` should be array-likes or integers representing the indices
+    of the channel pairs in the data for each connection. ``seeds`` and ``targets`` will
+    be expanded such that connectivity will be computed between each seed and each
+    target. E.g. the seeds and targets::
 
             seeds   = [0, 1]
             targets = [2, 3, 4]
@@ -244,28 +245,34 @@ def seed_target_multivariate_indices(seeds, targets):
 
     Parameters
     ----------
-    seeds : array of array of int, shape (n_unique_seeds, variable)
-        Seed indices.
-
-    targets : array of array of int, shape (n_unique_targets, variable)
-        Target indices.
+    seeds : array_like
+        Indices of signals for which to compute connectivity from. Has
+        ``n_unique_seeds`` array-like entries containing the channel indices for each
+        seed. The number of channels within each seed can vary.
+    targets : array_like
+        Indices of signals for which to compute connectivity to. Has
+        ``n_unique_targets`` array-like entries containing the channel indices for each
+        target. The number of channels within each target can vary.
 
     Returns
     -------
-    indices : tuple of array of array of int, shape (2, n_cons, variable)
-        The indices as a numpy object array.
+    indices : tuple of length 2 of array
+        A tuple of two numpy object arrays containing the seed and target indices,
+        respectively, to use for connectivity computation. Each array has ``n_cons``
+        array entries containing the channel indices for each connection, where
+        ``n_cons = n_unique_seeds * n_unique_targets``. The number of channels within
+        each connection can vary.
 
     Notes
     -----
-    ``seeds`` and ``targets`` should be array-likes representing the indices of
-    the channel sets in the data for each connection. The indices for each
-    connection should be an array-like of integers representing the individual
-    channels in the data. The length of indices for each connection do not need
-    to be equal. Furthermore, all indices within a connection must be unique.
+    ``seeds`` and ``targets`` should be array-likes representing the indices of the
+    channel sets in the data for each connection. The indices for each connection should
+    be an array-like of integers representing the individual channels in the data. The
+    length of indices for each connection do not need to be equal. Furthermore, all
+    indices within a connection must be unique.
 
-    Because the number of channels per connection can vary, the indices are
-    stored as numpy arrays with ``dtype=object``. E.g. ``seeds`` and
-    ``targets``::
+    Because the number of channels per connection can vary, the indices are stored as
+    numpy arrays with ``dtype=object``. E.g. ``seeds`` and ``targets``::
 
             seeds   = [[0]]
             targets = [[1, 2], [3, 4, 5]]
@@ -275,12 +282,12 @@ def seed_target_multivariate_indices(seeds, targets):
             indices = (np.array([[0   ], [0      ]], dtype=object),  # seeds
                        np.array([[1, 2], [3, 4, 5]], dtype=object))  # targets
 
-    Even if the number of channels does not vary, the indices will still be
-    stored as object arrays for compatibility.
+    Even if the number of channels does not vary, the indices will still be stored as
+    object arrays for compatibility.
 
-    More information on working with multivariate indices and handling
-    connections where the number of seeds and targets are not equal can be
-    found in the :doc:`../auto_examples/handling_ragged_arrays` example.
+    More information on working with multivariate indices and handling connections where
+    the number of seeds and targets are not equal can be found in the
+    :doc:`../auto_examples/handling_ragged_arrays` example.
     """
     array_like = (np.ndarray, list, tuple)
 
@@ -309,22 +316,21 @@ def degree(connectivity, threshold_prop=0.2):
 
     Parameters
     ----------
-    connectivity : ndarray, shape (n_nodes, n_nodes)
+    connectivity : array, shape (n_nodes, n_nodes) | Connectivity
         The connectivity matrix.
     threshold_prop : float
-        The proportion of edges to keep in the graph before
-        computing the degree. The value should be between 0
-        and 1.
+        The proportion of edges to keep in the graph before computing the degree. The
+        value should be between 0 and 1.
 
     Returns
     -------
-    degree : ndarray, shape (n_nodes,)
+    degree : array, shape (n_nodes,)
         The computed degree.
 
     Notes
     -----
-    During thresholding, the symmetry of the connectivity matrix is
-    auto-detected based on :func:`numpy.allclose` of it with its transpose.
+    During thresholding, the symmetry of the connectivity matrix is auto-detected based
+    on :func:`numpy.allclose` of it with its transpose.
     """
     from mne_connectivity.base import BaseConnectivity
 
diff --git a/mne_connectivity/vector_ar/model_selection.py b/mne_connectivity/vector_ar/model_selection.py
index b507b1235..e027c33aa 100644
--- a/mne_connectivity/vector_ar/model_selection.py
+++ b/mne_connectivity/vector_ar/model_selection.py
@@ -9,22 +9,21 @@
 def select_order(X, maxlags=None):
     """Compute lag order selections based on information criterion.
 
-    Selects a lag order based on each of the available information
-    criteria.
+    Selects a lag order based on each of the available information criteria.
 
     Parameters
     ----------
-    X : np.ndarray, shape (n_times, n_channels)
+    X : array, shape (n_times, n_channels)
         Endogenous variable, that predicts the exogenous.
     maxlags : int
-        The maximum number of lags to check. Will then check from
-        ``1`` to ``maxlags``. If None, defaults to
-        ``12 * (n_times / 100.)**(1./4)``.
+        The maximum number of lags to check. Will then check from ``1`` to ``maxlags``.
+        If ``None``, defaults to ``12 * (n_times / 100.)**(1./4)``.
 
     Returns
     -------
     selected_orders : dict
-        The selected orders based on the following information criterion.
+        The selected orders based on the following information criterion:
+
         * aic : Akaike
         * fpe : Final prediction error
         * hqic : Hannan-Quinn
@@ -75,7 +74,7 @@ def _logdet_symm(m):
 
     Parameters
     ----------
-    m : np.ndarray, shape (N, N)
+    m : array, shape (N, N)
         2d array that is positive-definite (and symmetric)
 
     Returns
@@ -96,7 +95,7 @@ def _sigma_u_mle(df_resid, nobs, sigma_u):
         Number of observations minus number of estimated parameters.
     nobs : int
         Number of observations/samples in the dataset.
-    sigma_u : np.ndarray, shape (n_channels, n_channels)
+    sigma_u : array, shape (n_channels, n_channels)
         Estimate of white noise process variance
 
     Returns
@@ -112,11 +111,11 @@ def _info_criteria(params, X, sigma_u, lags):
 
     Parameters
     ----------
-    params : np.ndarray, shape (lags, n_channels, n_channels)
+    params : array, shape (lags, n_channels, n_channels)
         The coefficient state matrix that governs the linear system (VAR).
-    X : np.ndarray (n_times, n_channels)
+    X : array (n_times, n_channels)
         Endogenous variable, that predicts the exogenous.
-    sigma_u : np.ndarray, shape (n_channels, n_channels)
+    sigma_u : array, shape (n_channels, n_channels)
         Estimate of white noise process variance
     lags : int
         Lags of the endogenous variable.
diff --git a/mne_connectivity/vector_ar/tests/test_var.py b/mne_connectivity/vector_ar/tests/test_var.py
index 6848df9a9..d14c5afc9 100644
--- a/mne_connectivity/vector_ar/tests/test_var.py
+++ b/mne_connectivity/vector_ar/tests/test_var.py
@@ -55,11 +55,11 @@ def create_noisy_data(
 
     Returns
     -------
-    sample_data : ndarray, shape (n_channels, n_samples)
+    sample_data : array, shape (n_channels, n_samples)
         Observed sample data. Possibly with noise.
-    sample_eigs : np.ndarray
+    sample_eigs : array
         The true eigenvalues of the system.
-    sample_A : np.ndarray
+    sample_A : array
         (Optional) if ``return_A`` is True, then returns the
         true linear system matrix.
     """
diff --git a/mne_connectivity/vector_ar/var.py b/mne_connectivity/vector_ar/var.py
index cab43ba15..afa765191 100644
--- a/mne_connectivity/vector_ar/var.py
+++ b/mne_connectivity/vector_ar/var.py
@@ -28,27 +28,25 @@ def vector_auto_regression(
 
     Parameters
     ----------
-    data : array-like, shape=(n_epochs, n_signals, n_times) | Epochs | generator
-        The data from which to compute connectivity. The epochs dimension
-        is interpreted differently, depending on ``'output'`` argument.
-    times : array-like
-        (Optional) The time points used to construct the epoched ``data``. If
-        ``None``, then ``times_used`` in the Connectivity will not be
-        available.
+    data : array_like, shape (n_epochs, n_signals, n_times) | Epochs | generator
+        The data from which to compute connectivity. The epochs dimension is interpreted
+        differently, depending on ``'output'`` argument.
+    times : array_like | None
+        The time points used to construct the epoched ``data``. If ``None``, then
+        ``times_used`` in the returned ``conn`` will not be available.
     %(names)s
-    lags : int, optional
-        Autoregressive model order, by default 1.
-    l2_reg : float, optional
-        Ridge penalty (l2-regularization) parameter, by default 0.0.
+    lags : int | None
+        Autoregressive model order (default 1).
+    l2_reg : float | None
+        Ridge penalty (l2-regularization) parameter (default 0.0).
     compute_fb_operator : bool
-        Whether to compute the backwards operator and average with
-        the forward operator. Addresses bias in the least-square
-        estimation :footcite:`Dawson_2016`.
-    model : str
-        Whether to compute one VAR model using all epochs as multiple
-        samples of the same VAR model ('avg-epochs'), or to compute
-        a separate VAR model for each epoch ('dynamic'), which results
-        in a time-varying VAR model. See Notes.
+        Whether to compute the backwards operator and average with the forward operator.
+        Addresses bias in the least-square estimation :footcite:`Dawson_2016`.
+    model : ``'dynamic'`` | ``'avg-epochs'``
+        Whether to compute one VAR model using all epochs as multiple samples of the
+        same VAR model (``'avg-epochs'``), or to compute a separate VAR model for each
+        epoch (``'dynamic'``, default), which results in a time-varying VAR model. See
+        Notes.
     %(n_jobs)s
     %(verbose)s
 
@@ -65,44 +63,41 @@ def vector_auto_regression(
 
     Notes
     -----
-    Names can be passed in, which are then used to instantiate the nodes
-    of the connectivity class. For example, they can be the electrode names
-    of EEG.
-
-    For higher-order VAR models, there are n_order ``A`` matrices,
-    representing the linear dynamics with respect to that lag. These
-    are represented by vertically concatenated matrices. For example, if
-    the input is data where n_signals is 3, then an order-1 VAR model will
-    result in a 3x3 connectivity matrix. An order-2 VAR model will result in a
-    6x3 connectivity matrix, with two 3x3 matrices representing the dynamics
-    at lag 1 and lag 2, respectively.
-
-    When computing a VAR model (i.e. linear dynamical system), we require
-    the input to be a ``(n_epochs, n_signals, n_times)`` 3D array. There
-    are two ways one can interpret the data in the model.
-
-    First, epochs can be treated as multiple samples observed for a single
-    VAR model. That is, we have $X_1, X_2, ..., X_n$, where each $X_i$
-    is a ``(n_signals, n_times)`` data array, with n epochs. We are
-    interested in estimating the parameters, $(A_1, A_2, ..., A_{order})$
-    from the following model over **all** epochs:
+    Names can be passed in, which are then used to instantiate the nodes of the
+    connectivity class. For example, they can be the electrode names of EEG.
+
+    For higher-order VAR models, there are ``n_order`` ``A`` matrices, representing the
+    linear dynamics with respect to that lag. These are represented by vertically
+    concatenated matrices. For example, if the input is data where ``n_signals`` is 3,
+    then an order-1 VAR model will result in a 3x3 connectivity matrix. An order-2 VAR
+    model will result in a 6x3 connectivity matrix, with two 3x3 matrices representing
+    the dynamics at lag 1 and lag 2, respectively.
+
+    When computing a VAR model (i.e. linear dynamical system), we require the input to
+    be a ``(n_epochs, n_signals, n_times)`` 3D array. There are two ways one can
+    interpret the data in the model.
+
+    First, epochs can be treated as multiple samples observed for a single VAR model.
+    That is, we have :math:`X_1, X_2, ..., X_n`, where each :math:`X_i` is a
+    ``(n_signals, n_times)`` data array, with ``n_epochs``. We are interested in
+    estimating the parameters, :math:`(A_1, A_2, ..., A_{order})`, from the
+    following model over **all** epochs:
 
     .. math::
         X(t+1) = \sum_{i=0}^{order} A_i X(t-i)
 
     This results in one VAR model over all the epochs.
 
-    The second approach treats each epoch as a different VAR model,
-    estimating a time-varying VAR model. Using the same
-    data as above,  we now are interested in estimating the
-    parameters, $(A_1, A_2, ..., A_{order})$ for **each** epoch. The model
-    would be the following for **each** epoch:
+    The second approach treats each epoch as a different VAR model, estimating a
+    time-varying VAR model. Using the same data as above,  we now are interested in
+    estimating the parameters, :math:`(A_1, A_2, ..., A_{order})` for **each** epoch.
+    The model would be the following for **each** epoch:
 
     .. math::
         X(t+1) = \sum_{i=0}^{order} A_i X(t-i)
 
-    This results in one VAR model for each epoch. This is done according
-    to the model in :footcite:`li_linear_2017`.
+    This results in one VAR model for each epoch. This is done according to the model in
+    :footcite:`li_linear_2017`.
 
     *b* is of shape [m, m*p], with sub matrices arranged as follows:
 
@@ -116,13 +111,13 @@ def vector_auto_regression(
     | b_m0 | b_m1 | ...  | b_mm |
     +------+------+------+------+
 
-    Each sub matrix b_ij is a column vector of length p that contains the
-    filter coefficients from channel j (source) to channel i (sink).
+    Each sub matrix b_ij is a column vector of length p that contains the filter
+    coefficients from channel j (source) to channel i (sink).
 
-    In order to optimize RAM usage, the estimating equations are set up
-    by iterating over sample points. This assumes that there are in general
-    more sample points then channels. You should not estimate a VAR model
-    using less sample points then channels, unless you have good reason.
+    In order to optimize RAM usage, the estimating equations are set up by iterating
+    over sample points. This assumes that there are in general more sample points than
+    channels. You should not estimate a VAR model using less sample points than
+    channels, unless you have good reason.
 
     References
     ----------
@@ -130,7 +125,7 @@ def vector_auto_regression(
     """
     if model not in ["avg-epochs", "dynamic"]:
         raise ValueError(
-            f'"model" parameter must be one of ' f"(avg-epochs, dynamic), not {model}."
+            f'"model" parameter must be one of (avg-epochs, dynamic), not {model}.'
         )
 
     events = None
@@ -254,24 +249,22 @@ def _construct_var_eqns(data, lags, l2_reg=None):
 
     Parameters
     ----------
-    data : np.ndarray (n_epochs, n_signals, n_times)
+    data : array (n_epochs, n_signals, n_times)
         The multivariate data.
     lags : int
         The order of the VAR model.
     l2_reg : float, optional
-        The l2 penalty term for ridge regression, by default None, which
-        will result in ordinary VAR equation.
+        The l2 penalty term for ridge regression (default ``None``) which will result in
+        an ordinary VAR equation.
 
     Returns
     -------
-    X : np.ndarray
+    X : array
         The predictor multivariate time-series. This will have shape
-        ``(model_order * (n_times - model_order),
-        n_signals * model_order)``. See Notes.
-    Y : np.ndarray
+        ``(model_order * (n_times - model_order), n_signals * model_order)``. See Notes.
+    Y : array
         The predicted multivariate time-series. This will have shape
-        ``(model_order * (n_times - model_order),
-        n_signals * model_order)``. See Notes.
+        ``(model_order * (n_times - model_order), n_signals * model_order)``. See Notes.
 
     Notes
     -----
@@ -279,8 +272,7 @@ def _construct_var_eqns(data, lags, l2_reg=None):
 
         Y = A X
 
-    where Y is time-shifted data copy of X and ``A`` defines
-    how X linearly maps to Y.
+    where Y is time-shifted data copy of X and ``A`` defines how X linearly maps to Y.
     """
     # n_epochs, n_signals, n_times
     n_epochs, n_signals, n_times = np.shape(data)
@@ -314,8 +306,8 @@ def _system_identification(data, lags, l2_reg=0, n_jobs=-1, compute_fb_operator=
     .. math::
         X(t+1) = \sum_{i=0}^{order} A_i X(t - i)
 
-    where ``data`` comprises of ``(n_signals, n_times)`` and ``X(t)`` are
-    the data snapshots.
+    where ``data`` comprises of ``(n_signals, n_times)`` and ``X(t)`` are the data
+    snapshots.
     """
     # 1. determine shape of the window of data
     n_epochs, n_nodes, n_times = data.shape
@@ -383,8 +375,8 @@ def _compute_lds_func(data, lags, l2_reg, compute_fb_operator):
 
     Note
     ----
-    The ``_estimate_var`` function returns a set of A matrices that represent
-    the system:
+    The ``_estimate_var`` function returns a set of A matrices that represent the
+    system:
 
         X(t+1) = X(t) A
 
@@ -392,8 +384,8 @@ def _compute_lds_func(data, lags, l2_reg, compute_fb_operator):
 
         X(t+1) = A X(t)
 
-    Therefore, a transpose is needed. If there are additional lags, then each
-    of these matrices need to be transposed.
+    Therefore, a transpose is needed. If there are additional lags, then each of these
+    matrices need to be transposed.
     """
     # make sure data is T x K (samples, coefficients) to make use of underlying
     # functions
@@ -420,29 +412,29 @@ def _estimate_var(X, lags, offset=0, l2_reg=0):
 
     Parameters
     ----------
-    X : np.ndarray (n_times, n_channels)
+    X : array (n_times, n_channels)
         Endogenous variable, that predicts the exogenous.
     lags : int
         Lags of the endogenous variable.
-    offset : int, optional
-        Periods to drop from the beginning of the time-series, by default 0.
-        Used for order selection, so it's an apples-to-apples comparison
+    offset : int
+        Periods to drop from the beginning of the time-series (default 0). Used for
+        order selection, so it's an apples-to-apples comparison
     l2_reg : int
-        The amount of l2-regularization to use. Default of 0.
+        The amount of l2-regularization to use (default 0).
 
     Returns
     -------
-    params : np.ndarray (lags, n_channels, n_channels)
+    params : array (lags, n_channels, n_channels)
         The coefficient state matrix that governs the linear system (VAR).
-    resid : np.ndarray
+    resid : array
         The residuals.
-    omega : np.ndarray (n_channels, n_channels)
+    omega : array (n_channels, n_channels)
         Estimate of white noise process variance
 
     Notes
     -----
-    This function was originally copied from statsmodels VAR model computation
-    and modified for MNE-connectivity usage.
+    This function was originally copied from statsmodels VAR model computation and
+    modified for MNE-connectivity usage.
     """
     # get the number of equations we want
     n_equations = X.shape[1]
@@ -571,16 +563,16 @@ def _get_var_predictor_matrix(y, lags):
 
     Parameters
     ----------
-    y : np.ndarray (n_samples, n_channels)
+    y : array (n_samples, n_channels)
         The passed in data array.
     lags : int
         The number of lags.
 
     Returns
     -------
-    Z : np.ndarray (n_samples, n_channels * lag_order)
-        Z is a (T x Kp) matrix, with K the number of channels,
-        p the lag order, and T the number of samples.
+    Z : array (n_samples, n_channels * lag_order)
+        Z is a (T x Kp) matrix, with K the number of channels, p the lag order, and T
+        the number of samples.
         Z := (Z_0, ..., Z_T).T (T x Kp)
         Z_t = [1 y_t y_{t-1} ... y_{t - p + 1}] (Kp x 1)
 
diff --git a/mne_connectivity/viz/_3d.py b/mne_connectivity/viz/_3d.py
index 6444b120b..3e51b4169 100644
--- a/mne_connectivity/viz/_3d.py
+++ b/mne_connectivity/viz/_3d.py
@@ -33,18 +33,18 @@ def plot_sensors_connectivity(
 
     Parameters
     ----------
-    info : dict | None
+    info : mne.Info
         The measurement info.
     con : array, shape (n_channels, n_channels) | Connectivity
-        The computed connectivity measure(s).
+        The connectivity data to plot.
     %(picks_good_data)s
         Indices of selected channels.
     cbar_label : str
         Label for the colorbar.
     n_con : int
-        Number of strongest connections shown. By default 20.
+        Number of strongest connections shown (default 20).
     cmap : str | instance of matplotlib.colors.Colormap
-        Colormap for coloring connections by strength. If :class:`str`, must be a valid
+        Colormap for coloring connections by strength. If a str, must be a valid
         Matplotlib colormap (i.e. a valid key of `matplotlib.colormaps`). Default is
         ``"RdBu"``.
     min_distance : float
diff --git a/mne_connectivity/viz/circle.py b/mne_connectivity/viz/circle.py
index 30aec8906..9b466be7d 100644
--- a/mne_connectivity/viz/circle.py
+++ b/mne_connectivity/viz/circle.py
@@ -42,44 +42,43 @@ def plot_connectivity_circle(
     Parameters
     ----------
     con : array | Connectivity
-        Connectivity scores. Can be a square matrix, or a 1D array. If a 1D
-        array is provided, "indices" has to be used to define the connection
-        indices.
+        Connectivity scores. Can be a square matrix, or a 1D array. If a 1D array is
+        provided, ``indices`` has to be used to define the connection indices.
     node_names : list of str
-        Node names. The order corresponds to the order in con.
-    indices : tuple of array | None
-        Two arrays with indices of connections for which the connections
-        strengths are defined in con. Only needed if con is a 1D array.
+        Node names. The order corresponds to the order in ``con``.
+    indices : tuple of array_like | None
+        Two arrays with indices of connections for which the connections strengths are
+        defined in ``con``. Only needed if ``con`` is a 1D array.
     n_lines : int | None
-        If not None, only the n_lines strongest connections (strength=abs(con))
-        are drawn.
+        If not ``None``, only the ``n_lines`` strongest connections
+        (``strength=abs(con)``) are drawn.
     node_angles : array, shape (n_node_names,) | None
-        Array with node positions in degrees. If None, the nodes are equally
-        spaced on the circle. See mne.viz.circular_layout.
+        Array with node positions in degrees. If ``None``, the nodes are equally spaced
+        on the circle. See :func:`mne.viz.circular_layout`.
     node_width : float | None
-        Width of each node in degrees. If None, the minimum angle between any
-        two nodes is used as the width.
+        Width of each node in degrees. If ``None``, the minimum angle between any two
+        nodes is used as the width.
     node_height : float
-        The relative height of the colored bar labeling each node. Default 1.0
-        is the standard height.
+        The relative height of the colored bar labeling each node. Default 1.0 is the
+        standard height.
     node_colors : list of tuple | list of str
-        List with the color to use for each node. If fewer colors than nodes
-        are provided, the colors will be repeated. Any color supported by
-        matplotlib can be used, e.g., RGBA tuples, named colors.
+        List with the color to use for each node. If fewer colors than nodes are
+        provided, the colors will be repeated. Any color supported by matplotlib can be
+        used, e.g., RGBA tuples, named colors.
     facecolor : str
-        Color to use for background. See matplotlib.colors.
+        Color to use for background. See :mod:`matplotlib.colors`.
     textcolor : str
-        Color to use for text. See matplotlib.colors.
+        Color to use for text. See :mod:`matplotlib.colors`.
     node_edgecolor : str
-        Color to use for lines around nodes. See matplotlib.colors.
+        Color to use for lines around nodes. See :mod:`matplotlib.colors`.
     linewidth : float
         Line width to use for connections.
     colormap : str | instance of matplotlib.colors.LinearSegmentedColormap
         Colormap to use for coloring the connections.
     vmin : float | None
-        Minimum value for colormap. If None, it is determined automatically.
+        Minimum value for colormap. If ``None``, it is determined automatically.
     vmax : float | None
-        Maximum value for colormap. If None, it is determined automatically.
+        Maximum value for colormap. If ``None``, it is determined automatically.
     colorbar : bool
         Display a colorbar or not.
     title : str
@@ -96,15 +95,15 @@ def plot_connectivity_circle(
         Font size to use for colorbar.
     padding : float
         Space to add around figure to accommodate long labels.
-    ax : instance of matplotlib PolarAxes | None
+    ax : instance of matplotlib.projections.polar.PolarAxes | None
         The axes to use to plot the connectivity circle.
     interactive : bool
-        When enabled, left-click on a node to show only connections to that
-        node. Right-click shows all connections.
+        When enabled, left-click on a node to show only connections to that node.
+        Right-click shows all connections.
     node_linewidth : float
-        Line with for nodes.
+        Line width for nodes.
     show : bool
-        Show figure if True.
+        Show figure if ``True``.
 
     Returns
     -------
@@ -115,16 +114,16 @@ def plot_connectivity_circle(
 
     Notes
     -----
-    This code is based on a circle graph example by Nicolas P. Rougier
+    This code is based on a circle graph example by Nicolas P. Rougier.
 
-    By default, :func:`matplotlib.pyplot.savefig` does not take ``facecolor``
-    into account when saving, even if set when a figure is generated. This
-    can be addressed via, e.g.::
+    By default, :func:`matplotlib.pyplot.savefig` does not take ``facecolor`` into
+    account when saving, even if set when a figure is generated. This can be addressed
+    via, e.g.::
 
     >>> fig.savefig(fname_fig, facecolor='black') # doctest:+SKIP
 
-    If ``facecolor`` is not set via :func:`matplotlib.pyplot.savefig`, the
-    figure labels, title, and legend may be cut off in the output figure.
+    If ``facecolor`` is not set via :func:`matplotlib.pyplot.savefig`, the figure
+    labels, title, and legend may be cut off in the output figure.
     """
     from mne_connectivity.base import BaseConnectivity