Facilitate export to XArray

Describe the problem

XArray is very convenient when doing data exploration as it often allows for a more transparent code by using dimension names rather than indexes. It would be very useful if MNE was providing XArray export functions where it already provides exports to NumPy array. For example, whereas someone would currently use

ica.get_sources(epochs).get_data().mean(2).mean(0)

by adding a to_xarray() method, we could use the following alternative code

ica.get_sources(epochs).to_xarray().mean(["time", "epochs"])
#  or ica.get_sources(epochs).to_xarray().mean("time").mean("epochs")

Personally, I found the second version much easier to understand. With the first version, I would typically have to print the shape of the data to infer what axis corresponds to what to understand what operation is being performed. The first form (with dimension indexes) is also much more error-prone because all sorts of things can affect the axes ordering whereas dimension names and coordinates generally remain consistent and generate much more transparent errors when they don't.

Start by tracing the existing ica.get_sources(epochs).get_data() export path and the other NumPy export entry points referred to in the issue. Define the scope of the proposed to_xarray() methods and identify the relevant tests; done should provide equivalent data with stable named dimensions and coordinates.