An ellipsoid has the equation
$$ (\mathbf{r} - \mathbf{r_0} )^T Q_e (\mathbf{r} - \mathbf{r_0} ) = 1 $$
It is positioned such that all of the ellipsoid lies above a projection plane whose equation is
$$ \mathbf{n}^T (\mathbf{r} - \mathbf{r_1}) = 0 $$
where $\mathbf{n}$ is a unit vector. The ellipsoid is subject to uniform direction light (like sun rays) that has a unit direction vector $\mathbf{d_0}$, producing a shadow of the ellipsoid on the given plane.
I want to derive a closed form expression of the area of the shadow.
My approach: is detailed in my answer below.
Your comments, or alternative answers are highly appreciated.
The shadow of the ellipsoid on the given plane is generated by the intersection of the plane and the elliptical cylinder that is tangent to the ellipsoid, and which has its axis parallel to the fixed direction of light rays $\mathbf{d_0}$.
Since the normal vector to the surface of the ellipsoid is given by $2 Q_e (\mathbf{r} - \mathbf{r_0})$, it follows that tangency points on the ellipsoid are characterized by the vector equation,
$ \mathbf{d_0}^T Q_e (\mathbf{r} - \mathbf{r_0} ) = 0 $
This is an equation of a plane passing through the center of the ellipsoid $\mathbf{r_0}$ and have a normal vector of $Q_e \mathbf{d_0}$. The intersection between this plane and the ellipsoid is of course an ellipse. To find its area, consider a family of cutting planes all having the unit normal $\mathbf{n_1} = \dfrac{ Q_e \mathbf{d_e} }{\sqrt{ \mathbf{d_0}^T Q_e^2 \mathbf{d_0} }} $ and passing through the point $\mathbf{p_1}(t) = \mathbf{r_0} + t \ \mathbf{ n_1 } $. Let $V(t)$ be the volume of the portion of the ellipsoid that is on the side of the plane from which unit vector $\mathbf{n_1} $ is originating, then the area of the cut $A(t)$ is equal to the derivative of $V(t)$ with respect to $t$.
Now it is rather simple to find $V(t)$. To this end, scale down the ellipsoid to a unit sphere, as follows:
Let the symmetric matrix $Q_e$ be decomposed as follows
$ Q_e = R D R^T = R \sqrt{D} \sqrt{D} R^T $
where $D$ is a diagonal matrix of $Q_e$'s eigenvalues, and $R$ is a rotation (orthogonal) matrix. Define $\mathbf{r'} = \sqrt{D} R^T (\mathbf{r} - \mathbf{r_0}) $, then it follows that $\mathbf{r'}^T\mathbf{r'} = 1$, that is, $\mathbf{r'}$ lies on the unit sphere with center at the origin.
The original cutting plane is the plane $\mathbf{n_1}^T ( \mathbf{r} - (\mathbf{r_0} + t \mathbf{n_1} ) ) = 0 $. The transformation of this plane into the $\mathbf{r'}$ space is obtained by substituting
$ \mathbf{r} = \mathbf{r_0} + R \sqrt{D^{-1}} \mathbf{r'} $
into the equation of the plane, this results in
$\mathbf{n_1}^T ( R \sqrt{D^{-1}} \ \mathbf{r'} - t \mathbf{n_1} ) = 0 $
The distance between this plane (now in the $\mathbf{r'}$ space), and the origin is
$$ a = \dfrac{ t }{\sqrt{ \mathbf{n_1}^T R D^{-1} R^T \mathbf{n_1} }} $$
But $ \mathbf{n_1} = \dfrac{Q_e \mathbf{d_0} } {\sqrt{ \mathbf{d_0}^T Q_e^2 \mathbf{d_0} }} $. It follows that
$ a = t \ \dfrac{ \sqrt{ \mathbf{d_0}^T Q_e^2 \mathbf{d_0} }}{ \sqrt{ \mathbf{d_0}^T Q_e \mathbf{d_0} }} =: t a_0 $
Now the volume of the portion of the unit sphere lying on the side of the plane from which the normal $\mathbf{n'} = \sqrt{D^{-1}} R^T \mathbf{n_1} $ is originating is given by (this is very simple to prove)
$V_1(t) = \pi (\dfrac{2}{3} + a - \dfrac{1}{3} a^3 ) = \pi ( \dfrac{2}{3} + a_0 t - \dfrac{1}{3} a_0^3 t^3 ) $
The volume of the corresponding portion of the original ellipsoid is
$ V(t) = \det(\sqrt{D^{-1}}) V_1(t) = \sqrt{ \det(D^{-1}) } \ V_1(t) = \dfrac{ V_1(t) }{ \sqrt{ \det(D) } } = \dfrac{ V_1(t)}{\sqrt{ \det(Q_e) }} $
Note that the cutting plane we're interested in is at $t=0$, therefore, differentiate $V(t)$ and evaluate the derivative at $t = 0$. This gives
$ A_0 = \dfrac{ \pi a_0 } { \sqrt{ \det(Q_e) } } = \dfrac{ \pi } { \sqrt{ \det(Q_e) } } \ \dfrac{ \sqrt{ \mathbf{d_0}^T Q_e^2 \mathbf{d_0} } }{ \sqrt{ \mathbf{d_0}^T Q_e \mathbf{d_0} } } $
We're almost there. Note the normal to the cut plane is along $\mathbf{n_1}$ which is parallel to $Q_e \mathbf{d_0} $, while the axis of the cylinder which is orthogonal to its cross section is along $\mathbf{d_0} $. Thus if $\theta_1$ is the angle betwen these two vectors then the cross sectional area of the cylinder is
$A_c = A_0 | \cos \theta_1 | $
Using inner product, we know that
$ \cos \theta_1 = \dfrac{ \mathbf{d_0}^T Q_e \mathbf{d_0}}{\sqrt{ \mathbf{d_0}^T Q_e^2 \mathbf{d_e} }}$ (remember that $\mathbf{d_0}$ is a unit vector.
Finally, the shoadow is related to the cross sectional area of the shadow cylinder by
$ A_{\text{shadow}} = \dfrac{A_c}{ |\cos \theta_2| } $
where $\theta_2 $ is the angle between the axis of the cylinder $\mathbf{d_0}$ and the normal vector to the projection plane (vector $\mathbf{n}). So $\cos \theta_2 = \mathbf{n}^T \mathbf{d_0} $
Putting it all together, we obtain,
$ A_{\text{shadow}} = \dfrac{\pi}{\sqrt{ \det(Q_e) }} \ \dfrac{ \sqrt{ \mathbf{d_0}^T Q_e^2 \mathbf{d_0} }}{ \sqrt{ \mathbf{d_0}^T Q_e \mathbf{d_0} } } \ \dfrac{ \mathbf{d_0}^T Q_e \mathbf{d_0} }{\sqrt{ \mathbf{d_0}^T Q_e^2 \mathbf{d_e} } } \dfrac{1}{| \mathbf{n}^T \mathbf{d_0} | } $
Simplifying we obtain the final form
$ A_{\text{shadow} } = \pi \dfrac{ \sqrt{ \mathbf{d_0}^T Q_e \mathbf{d_0} } }{ \sqrt{ \det(Q_e) } \ | \mathbf{n}^T \mathbf{d_0} | }$