The most direct way to get a small time expansion is using a gradient expansion. This can be worked out in the following way. Rescale your time variable as $t\rightarrow\lambda t$ being $\lambda$ a parameter taken to be arbitrary large and introduced to fix expansion order. Then, takes $\phi$ proportional to $\lambda$. You will get
$$\lambda u_\tau=\alpha(x,\frac{\tau}{\lambda},u)u_{xx}+\beta(x,\frac{\tau}{\lambda},u)(u_x)^2+\gamma(x,\frac{\tau}{\lambda},u)u_x+\lambda\phi(x,\frac{\tau}{\lambda},u).$$
Then, expand $u$ as
$$u(x,\tau)=\sum_{n=0}^\infty\frac{1}{\lambda^n}u_n(x,\tau).$$
At the end of the computation just put $\lambda=1$. Your solution will be given as polynomials in $t$.
Now, let us consider your equation
$$tu_t=Tr[\frac{u}{2}AD(\frac{x}{u})\otimes D(\frac{x}{u})]-t^2Tr[\frac{u^3}{8}ADu\otimes Du]+t<\alpha,Du>$$
$$+t Tr[\frac{1}{2}AD^2u]-\frac{u}{2}.$$
We firstly note that $D(\frac{x}{u})=\frac{e_x}{u}-x\frac{Du}{u^2}$ and so we can write down this equation in the form
$$tu_t=\frac{1}{2u}[A_{xx}+x^2Tr[A Du\otimes Du]]-t^2Tr[\frac{u^3}{8}ADu\otimes Du]+t<\alpha,Du>$$
$$+t Tr[\frac{1}{2}AD^2u]-\frac{u}{2}.$$
I hope I have interpreted correctly your notation. Now, multiply by $u$ both sides and you will get
$$tuu_t=\frac{1}{2}[A_{xx}+x^2Tr[A Du\otimes Du]]-t^2uTr[\frac{u^3}{8}ADu\otimes Du]+t<\alpha,Du>$$
$$+t uTr[\frac{1}{2}AD^2u]-\frac{u^2}{2}.$$
Now, change the variable as $\tau=\lambda t$. You will get the new equation
$$\tau uu_\tau=\frac{1}{2}[A_{xx}+x^2Tr[A Du\otimes Du]]-\frac{1}{\lambda^2}\tau^2uTr[\frac{u^3}{8}ADu\otimes Du]+\frac{1}{\lambda}\tau <\alpha,Du>$$
$$+\frac{1}{\lambda}\tau uTr[\frac{1}{2}AD^2u]-\frac{u^2}{2}.$$
From this you can read off immediately the leading order being
$$\tau u_0u_{0\tau}=\frac{1}{2}[A_{xx}+x^2Tr[A Du_0\otimes Du_0]]-\frac{u_0^2}{2}.$$
This is a Hamilton-Jacobi equation that can be solved by the characteristic method. Higher orders can be obtained straightforwardly as an expansion in $\frac{1}{\lambda}$.
We can spend a few words about the next-to-leading order by substituting into the scaled equation
$$u=u_0+\frac{1}{\lambda}u_1+O\left(\frac{1}{\lambda^2}\right)$$
Then, the equation becomes
$$\tau (u_0+\frac{1}{\lambda}u_1)(u_{0\tau}+\frac{1}{\lambda}u_{1\tau})=\frac{1}{2}[A_{xx}+x^2Tr[A D(u_0+\frac{1}{\lambda}u_1)\otimes D(u_0+\frac{1}{\lambda}u_1)]]-$$
$$\frac{1}{\lambda^2}\tau^2(u_0+\frac{1}{\lambda}u_1)Tr[\frac{(u_0+\frac{1}{\lambda}u_1)^3}{8}AD(u_0+\frac{1}{\lambda}u_1)\otimes D(u_0+\frac{1}{\lambda}u_1)]$$
$$+\frac{1}{\lambda}\tau <\alpha,D(u_0+\frac{1}{\lambda}u_1)>$$
$$+\frac{1}{\lambda}\tau (u_0+\frac{1}{\lambda}u_1)Tr[\frac{1}{2}AD^2(u_0+\frac{1}{\lambda}u_1)]-\frac{1}{2}\left(u_0+\frac{1}{\lambda}u_1\right)^2+O\left(\frac{1}{\lambda^2}\right).$$
So, the equation to compute $u_1$ is
$$\tau (u_0u_{1\tau}+u_1u_{0\tau})=x^2Tr[A Du_0\otimes Du_1]$$
$$-\tau <\alpha,Du_0>+\tau u_0Tr[\frac{1}{2}AD^2u_0]-u_0u_1.$$
One can repeat this procedure at any order and, finally, gets the solution into the form given at the beginning of this post.
Let me change your notations. You deal with a 1D quasilinear system with size $N=4$: the standard Cauchy problem is
$$
\frac{\partial u}{\partial t}+A(t,x,u)\frac{\partial u}{\partial x}= f(t,x),\quad u(t=0,x)=u_0(x),
$$
where $t\in \mathbb R$ (time variable) as well as $x$ (this is a 1D problem), $u$ is valued in $\mathbb R^N$, $A$ is a real-valued $N\times N$ matrix. You may also assume that $A$ depends smoothly (or even analytically) of its arguments.
(1) Your first case: all eigenvalues of $A_0:=A(0,x_0, u_0(x_0))$ are real and distinct, this is indeed the strictly hyperbolic case. In this case, you can guarantee local existence, uniqueness and continuous dependence on the data, i.e. local well-posedness. Of course you cannot expect global existence in general because of the nonlinearity (think about the scalar Burgers).
(2) Let me skip some of your cases and go directly to the case where $A_0$ has a non-real eigenvalue (and thus a pair of non-real eigenvalues). Big trouble ahead: even if the matrix $A$ is analytic, in which case, Cauchy-Kovalevskaya theorem is providing a local (unique) analytic solution, that solution is very unstable in the Hadamard sense. It means that even though $v_0-u_0$ is very small in a very strong topology, such as the $C^\infty$ topology or the $H^s$ topology for a very large $s$, you will not be able to control $u(t)-v(t)$ in a quite weak topology such as $L^2$ (all this is local of course). You will find precise statements in a paper by Métivier, Remarks on the well-posedness of the nonlinear Cauchy
problem, with MR number MR2127041.
(3) When all eigenvalues are real, some with multiplicity larger than 2, then instability could or could not occur, depending on other structural factors such as semi-simplicity of $A_0$. Generally speaking, multiple roots will trigger difficulties.
(4) A very important class of theorems, with the name of Lax-Mizohata theorems is establishing a weak converse to (1): if the problem is well-posed (e.g. meaning that you have some Sobolev norm control of $u(t)$ by some Sobolev norm of $u(0)$), then it implies that the system is
weakly hyperbolic (case (3)$\cup$(1)). So if you expect your system to be well-behaved in the sense of Hadamard, no choice, the roots must be real-valued, possibly with multiplicity.
Do not think that all physically relevant problem of that type are hyperbolic: to quote just one example, Van der Waals classical system is
$$
\partial_t u+\partial_x v=0\quad \partial_t v+\partial_x q(u)=0.
$$
When $q'(u)>0$, you are in a hyperbolic region, but when $q'(u)<0$, you have a non-real eigenvalue.
Best Answer
I highly doubt the result you actually asked for is true.
Consider the linear wave equation on $(1+3)$ Minkowski space. The analytic domain of influence of a point $x$ as Lax defined it, which morally says that $y$ is in the analytic domain only if one can find perturbations in arbitrary small neighborhoods of $x$ that change $y$ (if I interpret your question statement correctly), actually consists of only the null cone emanating from $x$ and nothing more, since strong Huygen's principle holds.
The same is true for the linear wave equation on $(1+(2k+1))$ Minkowski spaces.
The opposite conclusion can be drawn on $(1+2k)$ dimensional Minkowski spaces, where the Green's function have support inside the cone.
For more general situations, you may want to consult the classical result of Atiyah-Bott-Garding on the existence of Petrowsky lacunae. For any linear hyperbolic equation that admits a lacuna, the analytic domain of influence cannot cover the entirety of the geometric one.
But for the result that you seem to actually want, where you should replace the analytic domain of dependence by a suitable "convex" envelope of it, I don't know if such a result is proven anywhere, but my guess is that, at least for the "local" version one can approach it using some sort of geometric optics construction.
For possible references (I haven't actually finished reading either, so they may not contain what you want), maybe you want to look at Michael Beals' book on propagation of singularities (sorry, the title escapes me at the moment) or Rauch's notes on Hyperbolic PDEs and Geometric Optics which I think you can find floating around on the internet.