The identity
$$E(\mathbf{r},t)=A(\mathbf{r},t)e^{i(\langle\omega\rangle t-\phi(\mathbf{r},t))}\tag{1}$$
needs neither derivation nor justification; instead, it acts as an Ansatz for the electric field and as a definition for the pair of functions
$$A(\mathbf{r},t)e^{-i\phi(\mathbf{r},t)}:=E(\mathbf{r},t)e^{-i\langle\omega\rangle t}.\tag{1'}$$
Now, one of the weirder quirks of mathematics (and the mathematics of physics) is that you're free to define whatever you want, no matter how weird it might seem a priori. The only requirement is that you then go on to use those definitions to do something useful.
In this specific case, $E(\mathbf r,t)$ is a complex-valued function of both position and time, and since complex numbers have an amplitude and a phase, proposing an Ansatz of the form $E(\mathbf{r},t)=C(\mathbf{r},t)e^{-i\varphi(\mathbf{r},t)}$ would carry pretty much zero new information.
Your Ansatz in $(1)$, however, is different, because you're saying something nontrivial about the structure of the phase, namely that it's of the form $\varphi(\mathbf{r},t) = \phi(\mathbf{r},t)-⟨\omega⟩t$, where the variation in $\phi$ is much smaller than the central frequency. Here's the first core point: this is not guaranteed, i.e. there's plenty of imaginable waveforms for which there is no frequency $⟨\omega⟩$ such that that holds. (For examples, try superpositions of quasi-monochromatic waves at different central frequencies, or short pulses with a broad bandwidth and a strong chirp.) Similarly, there's no guarantee that you're going to be able to bound the time variation of the amplitude with reference to the central frequency. (Again, for examples, look to ultrashort pulses.)
Now, none of this is a problem, because we're not here to build a formalism that will handle every imaginable waveform. Instead, building on the auxiliary definitions in $(1)$, the bit that really does the work is the condition that
$$\left[\frac{1}{A} \frac{\partial A}{ \partial t},\frac{\partial \phi}{ \partial t} \right] \ll \langle\omega\rangle, \tag 3$$
and it is this that acts as the definition of quasi-monochromatic waves. Again, all you've done thus far is define things (in this case, the term quasi-monochromatic), so again, you don't actually need to justify anything*. Instead, if the author does their job correctly, the justification will come from showing that quasi-monochromatic waves, defined in this fashion, have useful properties (which they do).
*(extended) footnote:
OK, so maybe I lied a little at that point. You don't really need to justify the stuff you define, but if you're re-using terms that have previous connotations (or which partially overlap with such terms) then you do need to show that you're not radically changing those terms. For the case of quasi-monochromatic waves, you do need to show that your definition agrees with the intuitive understanding of the term.
There's two components of this, and they're both mathematical.
- One is a link between the relevant time derivatives, $\frac{1}{A} \frac{\partial A}{ \partial t}$ and $\frac{\partial \phi}{ \partial t}$, and the width of the power spectrum of the function $A(\mathbf{r},t)e^{-i\phi(\mathbf{r},t)}$.
- The other is the fact that multiplying a function by $e^{-i\langle\omega\rangle t}$ on the time domain is equivalent to shifting its frequency-domain representation by $⟨\omega⟩$, which is a corollary of the convolution theorem.
Both can be shown and they are relatively reasonable theorems, but I don't think the technicalities are that important here. Ultimately, those mathematical facts allow you to link your definition, $(3)$ with the physical fact that the bandwidth of your waveform is much smaller than its central frequency, which is about as close to the intuitive concept of quasi-monochromatic as you can reasonably get.
So, in a way, this last bit is the justification for the definition.
A general solution can be presented as
\begin{equation}
\int d^3 \vec k \Big[ g_1(\vec k)e^{i|v||\vec k|t-i\vec k\cdot \vec x}+g_2(\vec k)e^{-i|v||\vec k|t+i\vec k\cdot \vec x}\Big] \qquad (1)
\end{equation}
with arbitrary $g_1(k)$ and $g_2(k)$.
It follows from your formulas, after you express $\omega$ in terms of $k$. Two terms correspond to two ways of taking the square root. If you additionally require that the solution is real, then you need to impose
\begin{equation}
\Big(g_1(\vec k) \Big)^*=g_2(\vec k).
\end{equation}
In (1) the integral goes over spatial $k$ only. In principle, you can make it over all $k$ in space-time (I mean, including $\omega$). Then, as it is clear form what you wrote, you need to solve
\begin{equation}
g(\vec k,\omega)(\omega^2 - |v|^2 |\vec k|^2)=0.
\end{equation}
Its general solution is
\begin{equation}
g(\vec k,\omega) = f(\vec k,\omega) \delta(\omega^2 - |v|^2 |\vec k|^2), \qquad (2)
\end{equation}
where $\delta$ is the Dirac delta. Note that in (2) only values of $f(\vec k, \omega)$ at $\omega^2 - |v|^2 |\vec k|^2=0$ are relevant. Indeed, if we replace
\begin{equation}
f(\vec k,\omega)\to f(\vec k,\omega)+(\omega^2 - |v|^2 |\vec k|^2)\alpha (\vec k,\omega)
\end{equation}
with any $\alpha$, then (2) remains unchanged due to $\delta(x)x=0$. The above replacement allows to change values of $f$ arbitrarily everywhere away from $\omega^2 - |v|^2 |\vec k|^2=0$. Thus, only values of $f$ on $\omega^2 - |v|^2 |\vec k|^2=0$ actually contribute to (2).
The general solution than acquires the form
\begin{equation}
\int d^3\vec k d\omega f(\vec k,\omega) \delta(\omega^2 - |v|^2 |\vec k|^2) e^{-i\omega t + i\vec k\vec x}. \qquad (3)
\end{equation}
If you need a real solution, you need to impose
\begin{equation}
f(-\vec k,-\omega) = \Big( f(\vec k,\omega)\Big)^*.
\end{equation}
The solution in this form can be easily connected to (1) if one eliminates the $\omega$ integral using the delta function.
Best Answer
The entire phase factor is a function of $\mathbf{k}$. The center of the wave packet is where the individual waves interfere the most constructively. In order to maximize this constructive interference, the phase factors should be as close to each other as possible, and that occurs precisely at a stationary point, where the derivative of the phase factor with respect to $\mathbf{k}$ is zero. As we know from elementary calculus, there is no first-order change in a function at its stationary point.