How to draw dimension reduced high dimensional gaussians in 2D for EM algorithm visualization

Question

I'm implementing the EM algorithm. The visualization works for 2D features. I'd like to visualize it for higher dimensional data using dimension reduction(PCA)

Here k= 3. Each group of elipses are the gaussians. I'm drawing them with some slightly modified code taken from matplotlib:

from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms

#n_std = 3 : 98.9% of points
def confidence_ellipse(mean_x, mean_y, covariance, ax, n_std=3.0, edgecolor= 'red', **kwargs):
    #matplotlib method of drawing confidence elipse from mean and covar


    # Calculating the standard deviation of x from
    # the squareroot of the variance and multiplying
    # with the given number of standard deviations.
    scale_x = np.sqrt(covariance[0, 0]) * n_std

    # calculating the standard deviation of y 
    scale_y = np.sqrt(covariance[1, 1]) * n_std

    # apply the transformation to make normal elipse into whatever abomination the gaussian is 
    transf = transforms.Affine2D() \
        .rotate_deg(45) \
        .scale(scale_x, scale_y) \
        .translate(mean_x, mean_y)

    # Using a special case to obtain the eigenvalues of this
    # two-dimensional dataset.
    pearson = covariance[0, 1]/np.sqrt(covariance[0, 0] * covariance[1, 1])
    a = np.sqrt(1 + pearson)
    b = np.sqrt(1 - pearson)
    ellipse = Ellipse((0, 0), width=a * 2, height=b * 2,edgecolor=edgecolor, facecolor=None,fill=None, **kwargs)

    ellipse.set_transform(transf + ax.transData)
    return ax.add_patch(ellipse)
# just draw elipses at diff std values
def multi_eclipse(mean_x, mean_y, covariance, ax, **kwargs):
    for std, color in [(2,'r'),(2.4,'g'),(2.8, 'y'),(3, 'm')]:
        confidence_ellipse(mean_x, mean_y, covariance, ax, n_std=std, edgecolor=color)

Modified from: https://matplotlib.org/stable/gallery/statistics/confidence_ellipse.html

So the elipses are drawn based on variance and mean of each gaussian.

I've implemented 2D(I take two eigen vectors) PCA and am able to make a scatter plot of 6D data from it. I've run the EM algorithm on the 6D data. But I'm unsure how to make the multidimensional means and variances of each gaussian into 2D. For means I can just project them as they are just points. But I don't know what to do for covariances.

One possible simplification would be to group the 6D points to predicted gaussians and calculate mean, covariance from the projected points. But that will obviously be inaccurate and dependent on the data to work properly.