What is the vector form of an ellipsis' major and minor axis, given lissajous-figur notation for the ellipsis? A lissajous figure is described by 
$$ t \rightarrow \begin{pmatrix} A_x \sin(ω_1t + φ_1)\\A_y \sin(ω_2t + φ_2) \end{pmatrix}, t\in \left[0,\infty\right] $$
for $ω_1 = ω_2$ you get a simplified form of 
$$ t \rightarrow \begin{pmatrix} A_x \sin(ωt + φ)\\A_y \sin(ωt) \end{pmatrix}, t\in \left[0,\infty\right] $$
and the resulting shapes are ellipses. 
I am looking for the vector form of the (major and) minor axis of those ellipses, expressed through $$A_x, A_y, φ$$. I am interested in the general form, so also for tilted ellipses, where the ellipse's axis' are not parallel to the cartesian coordinat axis.
Update: 
With the matrix C as provided by @amd I was able to calcualte the eigenvalues and eigenvektors, and with those i could calculate the direction vector (and the angle) of the minor axis of the ellipse. I did not manage to do the same with the formula for $$tan(2\theta)$$ that he provided. The calculation that I ended up using was:
def calc_vector2(x_abs, y_abs, phase_shift):

square_tan = (np.square(x_abs) + np.square(y_abs)) * np.tan(phase_shift)
sqrt = np.sqrt( -(np.square(np.sin(phase_shift)) * (np.power(x_abs, 4) + 2 * np.cos(2 * phase_shift) * np.square(x_abs * y_abs) + np.power(y_abs,4))) / (np.square(np.sin(phase_shift)) - 1))
xysin = x_abs * y_abs * np.sin(phase_shift)
frac = 2 * np.square(xysin) * np.tan(phase_shift)

eigenvalue_1 = (square_tan + sqrt) / frac
eigenvalue_2 = (square_tan - sqrt) / frac

if eigenvalue_1 >= eigenvalue_2:
    x_dir = -(np.square(y_abs) * np.tan(phase_shift) - np.square(x_abs) * np.tan(phase_shift) + np.sqrt(
        -(np.square(np.sin(phase_shift)) * (
                np.power(x_abs, 4) + 2 * np.cos(2 * phase_shift) * np.square(x_abs * y_abs) + np.power(y_abs,
                                                                                                       4))) / (
            np.square(np.sin(phase_shift)) - 1))) / (2 * x_abs * y_abs * np.sin(phase_shift))
else:
    x_dir = (-np.square(y_abs) * np.tan(phase_shift) + np.square(x_abs) * np.tan(phase_shift) + np.sqrt(
        -(np.square(np.sin(phase_shift)) * (
                np.power(x_abs, 4) + 2 * np.cos(2 * phase_shift) * np.square(x_abs * y_abs) + np.power(y_abs,
                                                                                                       4))) / (
            np.square(np.sin(phase_shift)) - 1))) / (2 * x_abs * y_abs * np.sin(phase_shift))
y_dir = 1
angle_off_minor_axis = np.arctan(y_dir / x_dir)
print("Winkel gegen die X-Achse: " + str(angle_off_minor_axis / np.pi * 180))
return angle_off_minor_axis, [x_dir, y_dir]

 A: The squared distance from the origin is $ A_x^2 \sin^2 (\omega t  + \varphi) + A_y^2 \sin^2 (\omega t)$, which can be rewritten as
\begin{align}
&\frac 12 A_x^2(1-\cos (2(\omega t+\varphi)))+\frac 12A_y^2(1-\cos (2\omega t))=\\=&\frac12(A_x^2+A_y^2)-\frac 12(A_x^2\cos(2\varphi)+A_y^2)\cos(2\omega t)+\frac 12A_x^2\sin(2\varphi)\sin(2\omega t)=\\=&C_0+C_1\cos(2\omega t)+C_2\sin(2\omega t)
\end{align}
If $C_1=C_2=0$ then the curve is a circle. Otherwise there is an angle $\psi$ such that 
$\cos \psi=\frac {C_1}{\sqrt {C_1^2+C_2^2}} \text{ and }\sin \psi=\frac {C_2}{\sqrt {C_1^2+C_2^2}{}{}}$ , which leads to a squared distance $C_0+\sqrt {C_1^2+C_2^2}\cos(2\omega t-\psi)$ with clear minima and maxima.
A: Calling your expression $h(t)$, it's clear that the shape traced out by $h$ is symmetric about the origin, so the ellipse-center is $(0,0)$. The squared distance from $h(t)$ to the origin will be maximized/minimized at each $u$-value for which $h(u)$ lies along the major/minor axis. Let's look at that:
$$
k(u) = \|h(u)\|^2 = A_x^2 \sin^2 (\omega t  + \phi) + A_y^2 \sin^2 (\omega t),
$$
Computing the derivative and setting to zero gives:
\begin{align}
0 &= k'(u) \\
& = 2\omega A_x^2 \sin (\omega t  + \phi)\cos (\omega t  + \phi)+ 2\omega A_y^2 \sin (\omega t)\cos (\omega t),
\end{align}
and dividing through by $\omega$ (your whole formula makes no sense if $\omega = 0$) we get
\begin{align}
0 
&= A_x^2 \cdot 2\sin (\omega t  + \phi)\cos (\omega t  + \phi)+ A_y^2 \cdot 2\sin (\omega t)\cos (\omega t)\\
&= A_x^2 \cdot \sin (2(\omega t  + \phi))+ A_y^2 \cdot \sin (2\omega t)\\
\end{align}
so that 
\begin{align}
-A_x^2 \cdot \sin (2(\omega t  + \phi)) &=  A_y^2 \cdot \sin (2\omega t)\\
-\frac{A_y^2}{A_x^2} = \frac{\sin (2(\omega t  + \phi))}{\sin (2\omega t)}
\end{align}
Defining $r(t) = \frac{\sin (2(\omega t  + \phi))}{\sin (2\omega t)}$, it's easy to see that $r$ is undefined when $t$ is a multiple of $\pi$, but also that $r$ is surjective onto $\Bbb R$, hence that there's some $t_0$ for which $r(t_0)$ is any particular value, in particular, when it equals the value prescribed by the left-hand side. 
Your point $h(t_0)$ will lie along either the major or the minor axis. 
Sadly, I don't know an easy expression for $t_0$, so this may be as far as you can get. 
A: Here’s a somewhat roundabout approach, but it involves far less trigonometric chicanery than others that I can think of, so for my part I’m less apt to make a silly error somewhere along the way.   
One way to interpret a Lissajous figure is as the projection of a sinusoid drawn on the surface of a cylinder. So, we can view this ellipse as the projection of the parametric curve $$\mathbf r(t) = \left(A_x\sin(\omega t+\phi),A_x\cos(\omega t+\phi),A_y\sin\omega t\right)$$ onto the $x$-$z$ plane. W.l.o.g. we can reparametrize this as $\left(A_x\sin(t+\delta),A_x\cos(t+\delta),A_y\sin t\right)$, since $\omega\ne0$ just controls how quickly the curve is traced out and $\delta = \phi/\omega$ determines where on the cylinder the sinusoid starts.  
The vectors $\mathbf r(0)$ and $\mathbf r(\pi/2)$ are obviously linearly independent, and $$\begin{vmatrix}\mathbf r(t)^T\\\mathbf r(0)^T\\\mathbf r(\pi/2)^T\end{vmatrix} = \begin{vmatrix}A_x\sin(t+\delta)&A_x\cos(t+\delta)&A_y\sin t \\ A_x\sin\delta& A_x\cos\delta&0 \\ A_x\cos\delta&-A_x\sin\delta&A_y\end{vmatrix}=0,$$ so $\mathbf r(t)$ is in their span $P$ for all $t$: this space curve is planar. Thus we can obtain it as the intersection of $P$ with the cylinder $x^2+y^2=A_x^2$. The vectors $\mathbf r(0)$ and $\mathbf r(\pi/2)$ happen to be semiaxes of this ellipse, which simplifies the following computations, but that particular choice of basis isn’t essential.   
Let $$M = \begin{bmatrix} \sin\delta & \cos\delta & 0 \\ \cos\delta & -\sin\delta & 0 \\ 0 & A_y/A_x & 0 \\ 0&0&1 \end{bmatrix}.$$ This matrix represents a homogeneous coordinate system on $P$: the spatial homogeneous coordinates $\mathbf X$ of any point on the plane are a linear combination of the columns of $M$, i.e., $\mathbf X=M\mathbf x$. If you prefer, you can think of that last formula as a parameterization of $P$.  
The cylinder is represented by the matrix $Q=\operatorname{diag}(1,1,0,-A_x^2)$, and in $M$-coordinates, its intersection with $P$ has the matrix $$C = M^TQM = \begin{bmatrix}1&0&0\\0&1&0\\0&0&-A_x^2\end{bmatrix}.$$ If you’re familiar with matrix representations of conics, this might look like the matrix of a circle to you, but remember that the basis vectors that we’re using for $P$ have unequal norms.  
We now obtain a homography between $P$ and the $x$-$z$ plane by multiplying $M$ by the appropriate orthogonal projection matrix, which just amounts to deleting the second row of $M$: $$H = \begin{bmatrix}1&0&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix}M = \begin{bmatrix}\sin\delta&\cos\delta&0\\0&A_y/A_x&0\\0&0&1\end{bmatrix}.$$ Finally we apply this homography to $C$ to obtain the matrix of the ellipse represented by the original parametric equations: $$H^{-T}CH^{-1} = \begin{bmatrix}\csc^2\delta & -{A_x\over A_y}\cot\delta\csc\delta & 0 \\ -{A_x\over A_y}\cot\delta\csc\delta & {A_x^2\over A_y^2}\csc^2\delta & 0 \\ 0&0&-A_x^2\end{bmatrix}$$ which we can normalize to $$C' = \begin{bmatrix}\frac1{(A_x\sin\delta)^2} & -\frac1{A_xA_y\tan\delta\sin\delta} & 0 \\ -\frac1{A_xA_y\tan\delta\sin\delta} & \frac1{(A_y\sin\delta)^2} & 0 \\ 0&0&-1\end{bmatrix}$$ corresponding to the Cartesian equation $${x^2\over(A_x\sin\delta)^2}-{2xy\over A_xA_y\tan\delta\sin\delta}+{y^2\over(A_y\sin\delta)^2}=1.$$ This is undefined for $\delta=n\pi$, but in that case we’re looking at $\mathbf r(t)$ edge on and the resulting Lissajous figure is a line segment rather than an ellipse.  
From here you can use standard formulas and techniques to analyze the ellipse: the eigenvalues of $C'$ are inverse squares of the semiaxis lengths; the ellipse is inclined at an angle $\theta$ from the horizontal, where $$\tan{2\theta} = {2A_xA_y\cos\delta\over A_x^2-A_y^2},$$ and so on.  

Eliminating $t$ directly from the original parametric equations yields the curious equation $\arcsin{\frac x{A_x}}=\arcsin{\frac y{A_y}}+\delta$. I haven’t pursued it any farther, but it does appear to have an interesting interpretation in terms of counterrotating points.  

Update: It struck me that eliminating $t$ from the parameterization isn’t all that tricky after all. We immediately get $\sin t=y/A_y$. Expand $x=A_x\sin(t+\delta)$ using the sum-of-angles formula, substitute for $\sin t$ and solve for $\cos t$, then use the usual $\cos^2t+\sin^2t=1$ to get $$\left({x\over A_x\sin\delta}-{y\over A_y\tan\delta}\right)^2+\left({y\over A_y}\right)^2=1.$$ Expanding and collecting terms results in the same equation as previously. Much less work than the original approach, but perhaps not as geometrically satisfying.
