add phase plot to plot_psd()

Describe the new feature or enhancement

I confronted this PSD plot from an eeg recording, which I don't know where this 25 Hz artifact is coming from, any idea on it?
So, I was thinking of adding phase spectrum next to the PSD plot, but AI suggests that adding raw phase spectrum is uniformly random for stochastic signals, so it is not interpretable, instead adding:

• "Inter-segment / inter-epoch phase consistency": which is similar to compute_tfr(output="itc")
• "Cross-channel relative phase": relative to a ref channel or to the common mode. A common-mode artifact gives ≈0 rad on all channels, so it could be best viewed with plot_topomap

Describe your proposed implementation

MNE already has compute_psd(output="complex") which returns per-segment Fourier coefficients but for all plotting functions it passes through _prepare_data_for_plot and silently discards phase information. I would suggest adding 2 new methods to BaseSpectrum so the existing plotting methods can be reused.

spec = raw.compute_psd(method="welch", output="complex", average=False)
pc = spec.phase_consistency()
rp = spec.relative_phase(ref="...")

pc.plot_topomap(...) 
rp.plot_topomap(...) 

Describe possible alternatives

None.

Additional context

No response

Start with BaseSpectrum and _prepare_data_for_plot, tracing how compute_psd(output="complex", average=False) currently loses phase data. Clarify the phase-consistency and relative-phase entry points and how existing plotting methods, including plot_topomap, should render them; done when the agreed behavior no longer silently discards the relevant phase information.