The question I believe correctly states some key constraints, particularly that the average of an ellipse with itself rotated $\pi/2$ should be a circle.
Another important constraint on a method for computing the average of a set of ellipses is that the computed average must be independent of the ordering of the ellipses.
Developing a more concrete answer to this question requires more concise definition of the properties of the mean of a set of ellipses. One way to approach this is to define an error metric for the difference between one ellipse and another.
For this, I will assume that all the ellipses are centered on the origin, with semi-major axis $a$ oriented along the cartesian $x$-axis and semi-minor axis $b$ oriented along the cartesian $y$-axis, and each ellipse may have a rotation of angle $\alpha$ (positive counter-clockwise) from the $+x$ axis.
Let the eccentricity of the ellipse be $e = \sqrt{1 - \frac{b^2}{a^2}}$. Using the formulation $\left[ \frac{\cos^2(\alpha)}{a^2} + \frac{\sin^2(\alpha)}{b^2} \right] x^2 + 2\cos(\alpha)\sin(\alpha)\left[\frac{1}{a^2} - \frac{1}{b^2}\right]xy + \left[ \frac{\sin^2(\alpha)}{a^2} + \frac{\cos^2(\alpha)}{b^2} \right] y^2 = 1$ from https://www.maa.org/external_archive/joma/Volume8/Kalman/General.html, if we then substitute \begin{align} x \rightarrow r\cos(\theta) \\ y \rightarrow r\sin(\theta) \end{align} then working through all the tedious algebra yields an explicit formula in polar coordinates for the radius $r$ as a function of angle $\theta$: $r = \frac{b}{\sqrt{1 + e^2\left(\sin^2(\alpha) + \cos(2\alpha)\sin^2(\theta) - \frac{1}{2}\sin(2\alpha)\sin(2\theta) - 1 \right)}} .$
This generalized polar form for an arbitrary ellipse is a naturally suited coordinate system for evaluating the similarity between one ellipse and another. You could use the absolute area between two ellipses, computed in polar coordinates, as your chosen error metric for example: $A_{diff} = \int_0^{2\pi}\frac{1}{2}\|r^2_1 - r^2_0 \| d\theta$
Probably a more practical error metric, since this ellipse mean in general will be computed for $N$ ellipses and most likely to be done numerically, might be (this is what I would use): $E = \sum_{i=1}^N E_i \\ E_i = \int_0^{2\pi} \left(r^2 - r_i^2 \right)^2 d\theta$ and numerically find a discrete form of $r$ that minimizes $E$, then solve for the $e$ and $b$ that fit $r$.
Such a formulation in terms of functional minimization will ensure the mean is independent of the ordering of the set of input ellipses, and should ensure that the mean of ellipse with itself rotated by $\pi/2$ is a circle (though I haven't proven or tested this).