How can I get the axes of the polarization ellipse from the Jones vector of the light?

Let me try a second time. I use https://math.stackexchange.com/questions/1204131/converting-a-rotated-ellipse-in-parametric-form-to-cartesian-form as a resource.

Prior Manipulations

The physical electric field is

$$\mathbf{E}_{phys} = Re\left[\mathbf{E} e^{i\omega t}\right] = Re\left[\begin{pmatrix}|E_x|e^{i\varphi_x}\\|E_y|e^{i\varphi_y}\end{pmatrix}e^{i\omega t}\right] = \begin{pmatrix}|E_x|\cos(\omega t + \varphi_x)\\|E_y| \cos(\omega t + \varphi_y)\end{pmatrix} = \begin{pmatrix}|E_x|\left[\cos(\omega t)\cos(\varphi_x) - \sin(\omega t) \sin(\varphi_x) \right]\\|E_y| \left[\cos(\omega t)\cos(\varphi_y) - \sin(\omega t) \sin(\varphi_y) \right]\end{pmatrix}$$

This is a parametric equation for an ellipse, which is traced by the electric field.

Major axes angle

Let

$$a=|E_x|\cos(\varphi_x)$$

$$b=|E_x|\sin(\varphi_x)$$

$$c=|E_y|\cos(\varphi_y)$$

$$d=|E_y|\sin(\varphi_y)$$

Then by comparison with the linked math.SE question the accepted answer states that the major and minor axes point on the ellipse (which is centred on the origin) fulfill

$$\omega t={1\over2}\arctan{2(ab+cd)\over(a^2+c^2)-(b^2+d^2)}+{k\pi\over2}\qquad(0\leq k\leq3)\ .$$

Expansion required by the OP

Indeed this quantity is the angle $\theta$ in the expansion appearing in the OPs question. However depending on which value for $k$ is chosen in may be $\pi/2 - \theta$, a case distinction depending on the sector of the arctan is necessary.

This of course also yields $\hat{u}$ and $\hat{v}$, so the expansion can now be easily obtained by projecting the Jones vector this basis.

The formulae may be along, but constitute a close form solution to the problem up to the case distinction of choosing $k$. I do not see how a simpler solution can exist, since the formula give for the angle given does not seem algebraically reducible.