I wonder why the displacement field $\mathbf{D}$ at a normal mode are parallel to the major and minor axis of the ellipse which is the intersection of index ellipsoid and plane orthogonal to the $\mathbf{k}$.
The Maxwell equation in anisotropic medium is $$\nabla \times \mathbf{E}=-\mu_0\frac{\partial \mathbf{H}}{\partial t}$$ $$\nabla \times \mathbf{H}= \frac{\partial \mathbf{D}}{\partial t}$$ Therefore, $$\nabla \times \nabla \times \frac{\epsilon_r ^{-1}}{\epsilon_0}\mathbf{D}=-\mu_0\frac{\partial^2\mathbf{D}}{\partial t^2}$$ In the form of plane wave, we can write above equation as below: $$-\mathbf{k}\times\mathbf{k}\times\epsilon_r ^{-1}\mathbf{D}=\mu_0\epsilon_0\omega^2\mathbf{D}=k_0^2 \mathbf{D}$$ where $\mathbf{k}=k\mathbf{\hat{u}}$ and $k=k_0n$. So we get $$-\mathbf{\hat{u}}\times\mathbf{\hat{u}}\times \epsilon_r^{-1}\mathbf{D}=\frac{1}{n^2} \mathbf{D}$$ Here the operator $(-\mathbf{\hat{u}}\times\mathbf{\hat{u}}\times)$ means the projection to the normal plane of $\mathbf{\hat{u}}$, denoted by $P_\mathbf{u}$. So this is an eigenvalue problem: $$P_\mathbf{u}\epsilon_r^{-1} \mathbf{D}=\frac{1}{n^2}\mathbf{D}$$ where $\mathbf{D}$ is an eigenvector with an eigenvalue $1/n^2$ of the linear map $P_\mathbf{u}\epsilon_r^{-1}$