tfr_array_multitaper with output="complex" is misleading

Describe the bug

the output of tfr_array_multitaper when output='complex' is misleading. Not sure if one should change something in the design, or just update the documentation.

Steps to reproduce

import mne.time_frequency as tfr
import numpy as np

signal = np.zeros(1000)
signal[500:] = 1

# morlet wavelets
tfr_out1 = tfr.tfr_array_morlet([[signal]], sfreq=100, freqs=[20], output='complex')
tfr_out2 = tfr.tfr_array_morlet([[signal]], sfreq=100, freqs=[20], output='power')
print(np.allclose(np.abs(tfr_out1)**2, tfr_out2))
print(tfr_out1.shape, tfr_out2.shape)

# using only 1 taper
tfr_out1 = tfr.tfr_array_multitaper([[signal]], sfreq=100, freqs=[20], time_bandwidth=2, output='complex')
tfr_out2 = tfr.tfr_array_multitaper([[signal]], sfreq=100, freqs=[20], time_bandwidth=2, output='power')
print(np.allclose(np.abs(tfr_out1)**2, tfr_out2))
print(tfr_out1.shape, tfr_out2.shape)

# using only 2 tapers
tfr_out1 = tfr.tfr_array_multitaper([[signal]], sfreq=100, freqs=[20], time_bandwidth=3, output='complex')
tfr_out2 = tfr.tfr_array_multitaper([[signal]], sfreq=100, freqs=[20], time_bandwidth=3, output='power')
print(np.allclose(np.abs(tfr_out1)**2, tfr_out2))
print(tfr_out1.shape, tfr_out2.shape)

Expected vs Actual results

Everything is as expected, apart from the last output: when using multitaper with more than 1 taper, one should get multiple complex time series outputs, one for each taper, so that one can reconstruct the power from the complex time series. With the current design, it is not possible, and it is not clear what is the meaning of the output returned by tfr.tfr_array_multitaper([[signal]], sfreq=100, freqs=[20], time_bandwidth=3, output='complex'). Is it the convolution with one of the tapers? Is it a sum of the complex time series from different tapers?

Actual results

One gets:

Morlet: True (1, 1, 1, 1000) (1, 1, 1, 1000)

multitaper with single taper: True (1, 1, 1, 1000) (1, 1, 1, 1000)

multitaper with 2 tapers: False (1, 1, 1, 1000) (1, 1, 1, 1000) (we need more information than that to reconstruct the "power" from "complex")

Additional information

Platform: macOS-10.15.7-x86_64-i386-64bit Python: 3.8.5 (default, Aug 5 2020, 03:39:04) [Clang 10.0.0 ] Executable: /Users/.../opt/anaconda3/envs/mne/bin/python CPU: i386: 12 cores Memory: 16.0 GB mne: 0.22.0 numpy: 1.19.4 {blas=NO_ATLAS_INFO, lapack=lapack} scipy: 1.5.3 matplotlib: 3.3.1 {backend=MacOSX} sklearn: 0.24.0 numba: Not found nibabel: 3.2.1 nilearn: 0.7.0 dipy: 1.3.0 cupy: Not found pandas: 1.1.5 mayavi: 4.7.2 pyvista: 0.27.4 {pyvistaqt=0.2.0, OpenGL 4.1 ATI-3.10.19 via AMD Radeon Pro 560X OpenGL Engine} vtk: 9.0.1 PyQt5: 5.12.3

mne-tools / mne-python