You can go from definition of linearity.
If $x$ and $y$ are solutions of differential equation defined by linear operator $H$, then $x+y$ is also solution.
In other words, if $Hx=0$ and $Hy=0$, then it should follow that $H(x+y)=0$
Another property of linear operators is scaling: $H(\alpha x)=\alpha H(x)$
In your case it is easy to see that operator defining equation is non-linear. Let's say that $x$ and $y$ are solutions to the equation. Let's assume that $(x+y)$ is also solution:
${\frac {\mathrm {d} ^{2}(x+y)}{\mathrm {d} t^{2}}}+2\zeta \omega _{0}{\frac {\mathrm {d} (x+y)}{\mathrm {d} t}}+\omega _{0}^{\,2}(x+y)=(-2\zeta \omega _{0}{\frac {\mathrm {d} (x+y)}{\mathrm {d} t}})^2$
If you expand it:
${\frac {\mathrm {d} ^{2}x}{\mathrm {d} t^{2}}}+{\frac {\mathrm {d} ^{2}y}{\mathrm {d} t^{2}}} +
2\zeta \omega _{0}{\frac {\mathrm {d} x}{\mathrm {d} t}}+
2\zeta \omega _{0}{\frac {\mathrm {d} y}{\mathrm {d} t}}+
\omega _{0}^{\,2}x+
\omega _{0}^{\,2}y
=4(\zeta \omega _{0}{\frac {\mathrm {d} x}{\mathrm {d} t}})^2+
4(\zeta \omega _{0}{\frac {\mathrm {d} y}{\mathrm {d} t}})^2+
8(\zeta \omega_{0})^2 {\frac {\mathrm {d} x}{\mathrm {d} t}}{\frac {\mathrm {d} y}{\mathrm {d} t}}$
You will see requirement for $(x+y)$ to be solution:
$8(\zeta \omega_{0})^2 {\frac {\mathrm {d} x}{\mathrm {d} t}}{\frac {\mathrm {d} y}{\mathrm {d} t}}=0$
Which is not automatically follows from $x$ and $y$ being solutions.
Hence, original operator (or original equation) is non-linear.