How to justify this small angle approximation $\dot{\theta}^2=0$? Suppose the equation of motion for some oscillating system takes the following form:
$$\ddot{\theta}+\dot{\theta}^2\sin\theta+k^2\theta\cos\theta=0$$
Applying small angle approximation to $\theta$ gives $\sin\theta\approx\theta$ and $\cos\theta\approx1$,
$$\ddot{\theta}+\dot{\theta}^2\theta+k^2\theta=0$$
But I am wondering if it is alright to set $\dot{\theta}^2=0$ so that the equation simplifies to
$$\ddot{\theta}+k^2\theta=0$$
If so, how can this approximation be justified?
 A: By drawing a phase plot you can more visually see where the approximations are valid. Define
\begin{align}
x&\equiv \theta\\
y&\equiv\dot\theta
\end{align}
This transforms your initial equation into the following first order differential equations:
\begin{align}
\pmatrix{\dot x\\\dot y}=\pmatrix{y\\-y^2\sin x-k^2x\cos x}
\end{align}
We can plot the vector on the right as a velocity field to see the dynamics. You are probably familiar with this. I use the following Mathematica code to do this:
range = 3;
expressions = {
   {{y, -k^2 x}, {y, -k^2 Sin[x]}},
   {{y, -y^2 x - k^2 x}, {y, -y^2 Sin[x] - k^2 x Cos[x]}}
   };
Grid@Table[
  StreamPlot[
   expressions[[i, j]] /. k -> 1, {x, -range, range}, {y, -range, 
    range}, PlotLabel -> MatrixForm@{expressions[[i, j]]}],
  {i, 1, 2}, {j, 1, 2}
  ]


I added the small angle approximation of the harmonic oscillator as well for comparison.
My interpretation: as long as $x$ and $y$ are both small the approximation is valid. The region of validity is fairly small though and outside the unit circle the approximation becomes pretty bad. The valid region will probably change for different $k$.
A: To think about the problem physically, we can imagine differential equations with the same dimensionality as linear force problems.  Your idealized approximation $\ddot\theta≈-k^2\theta$ is the frequently useful simple harmonic oscillator, given by a Hooke’s Law-type of linear restoring force. In linear coordinates, we would write this force like $m\vec a\propto - \vec x$.
What are some interactions that correspond to the term you’d like to approximate as zero?  Consider that
$$
m\vec a\propto -\vec v
$$
is usually considered a “drag” force, acting to bring the velocity towards zero.  In one sense we can already stop there: if you approximated $\dot\theta\to0$, you are eliminating some kind of a drag term from your equation, and turning a damped oscillator into an undamped oscillator.
Yours is not the usual damped oscillator, because the directionality comes from the position, not from the squared velocity:
$$
m\vec a ≈ -v^2\vec x
$$
A related remark: when you’re approximating by throwing away terms, you don’t generally want to say “this term should go to zero.” It’s better hygiene to say “this term is much smaller than that term,” which is a reason to keep only “that term.” You’ve already done this with your small angle approximation, which is
$$
\cos x = 1-\frac{x^2}{2!}+\frac{x^4}{4!}-\cdots≈1
$$
and which is justified in the limit $x^2\ll 1$, whether $x\to0$ or not.  The distinction between “much smaller than” versus “near zero” is especially important if you do some algebra and find that you actually have $1-\cos x$, in which case you need to be prepared to go back and re-do your approximation with more terms.
In this case, you have (in the linear-force language),
$$
m\vec a = -(v^2+k^2)\vec x
$$
and the relevant comparison is whether $v^2\ll k^2$ (or vice-versa).  You might remember that a simple harmonic oscillator has its minimum $v$ when $x$ is large, and its maximum $v$ as $x$ passes through the origin.  So the question of whether you can ignore $\dot\theta^2$ and treat your oscillator as undamped is a question which depends both on your spring constant $k$ and on your amplitude.  The damping becomes negligible if the oscillations are small and the spring is stiff.
Specifically, since the maximum speed for an undamped oscillator is $\dot\theta_\text{max}=k\theta_\text{max}$, you have
$$
\frac{\dot\theta_\text{max}^2}{k^2} ≈ \theta_\text{max}^2 
$$
which you have already assumed was small, in order to use the small-angle approximation.
An exercise for the reader: while you can ignore the damping over a single oscillation, you can’t ignore the damping over many oscillations. How many? My ballpark guess is $1/\theta_\text{max}$ or $1/\theta_\text{max}^2$, but I would expect to learn something from constructing the argument. The square would make a big difference in estimating the dissipation time of small oscillations.
A: There is no mathematical justification for setting $\dot \theta^2$ to $0$. If this equation represents some specific physical system and if it is known for physical reasons that $\dot \theta^2 \approx 0$ then it could be justified. But that will require a specific physical argument for the specific system in question, not a generic mathematical justification.
