mne.minimum_norm.source_band_induced_power has a label param, but in principle it should be possible to return results for multiple labels. Looking at the algorithm, our implementation appears to be:

1. Compute inverse operator K restricted to a label:

   mne-python/mne/minimum_norm/time_frequency.py

   Line 372 in 6a701d3

   K, sel, Vh, vertno, is_free_ori, noise_norm = _prepare_source_params(

2. Compute PSDs for each epoch in parallel:

   mne-python/mne/minimum_norm/time_frequency.py

   Lines 396 to 397 in 6a701d3

   	out = parallel(
   	my_compute_source_tfrs(

   Within each parallel split, take the sum across each epochs:

   mne-python/mne/minimum_norm/time_frequency.py

   Lines 292 to 294 in 9547c13

   	power += power_e
   	if with_plv:
   	plv += plv_e

3. Take the sum across parallel-splits-of-epochs:

   mne-python/mne/minimum_norm/time_frequency.py

   Lines 412 to 413 in 6a701d3

   	power = sum(o[0] for o in out)
   	power /= len(epochs_data) # average power over epochs

4. Compute PLV if requested (only used by source_induced_power, never source_band_induced_power), noise normalize (i.e., multiply each vertex by some value), baseline correct, return power and plv.

So I think we should in principle be able to add support for list-of-labels for source_band_induced_power and probably also source_induced_power when return_plv=False by doing the following:

1. Document the shape of the ndarray that is currently returned (something like (n_vertices, n_freqs, n_times))
2. Add return_plv=True to source_induced_power (the way it acts currently)
3. Fix documentation of source_induced_power to note that it returns both power and plv (when return_plv=True) and add return shapes
4. Fix documentation to use (baseline_mode)s) or whatever rather than repeat all options.
5. Add some small test that actually checks our source_band_induced_power values (e.g., assert_allclose(phase_lock[some_indices], [[some_array_values]])). This will be useful for all following changes.
6. Move noise_norm inside the parallel function.
7. Add support for list-of-label by computing a new label_op where label_op @ K averages across vertices to convert the (n_vertices, n_channels) operator K to be shape (n_labels, n_channels). This should only be allowed when with_plv/return_plv=False, which is always the case for source_band_induced_power and will now optionally be possible for source_band_power by step (3) above.
8. Add test that label=label with a .mean(axis=0) gives the same result as label=[label] when baseline_mode is None or mean, and hopefully sufficiently similar for other modes (shouldn't have a huge impact).

For point (8), it would be nice if list-of-label with one element [label] returned an identical to the label=label (single label) in all cases, but it wouldn't with these changes. The baseline correction in the list-of-label case will be done after the power is averaged within the label, whereas in the label=label case (like in our existing example) the averaging is done over already baseline-corrected vertex values within the label. I'm not certain which way is better, but 1) I doubt it makes much difference and 2) I think in principle averaging first is probably better / higher SNR, as averaging across vertices within a label first should in principle increase SNR before something like a z-transform, if used, could otherwise amplify noisy values. But either way I think if we explain the difference clearly in Notes I think we're okay.