How can we explain the "bad reputation" of higher-order polynomials? We all must have heard it by now - when we start learning about statistical models overfitting data, the first example we are often given is about "polynomial functions" (e.g., see the picture here):

We are warned that although higher-degree polynomials can fit training data quite well, they surely will overfit and generalize poorly to the test data.
Why does this happen? Is there a mathematical justification as to why (higher-degree) polynomial functions overfit the data? The closest explanation I could find online was something called "Runge's phenomenon", which suggests that higher-order polynomials tend to "oscillate" a lot - does this explain why polynomial functions are known to overfit data?
I understand that there is a whole field of "regularization" that tries to fix these overfitting problems (e.g., penalization can prevent a statistical model from "hugging" the data too closely) - but just using mathematical intuition, why are polynomials known to overfit the data?
In general, "functions" (e.g., the response variable you are trying to predict using machine learning algorithms) can be approximated using older methods like Fourier series, Taylor series and newer methods like neural networks. I believe that there are theorems that guarantee that Taylor series, polynomials and neural networks can "arbitrarily approximate" any function. Perhaps neural networks can promise smaller errors for simpler complexity?
But are there mathematical reasons behind higher-order polynomials (e.g., polynomial regression) being said to have a bad habit of overfitting, to the extent that they have become very unpopular? Is this solely explainable by Runge's phenomenon?
Reference:
Gelman, A. and Imbens, G. (2019) Why high order polynomials should not be used in regression discontinuity designs. Journal of Business and Economic Statistics 37(3), pp. 447-456. (An NBER working paper version is available here)

Before I attempt to answer this, I'd like to just point out that what you are observing here (overfitting in a linear regression) is just a specific example of a more general phenomenon, bias-variance trade-off, which is also observed in more "modern" machine learning contexts as well as the "classical" setting of regression fitting.  Your question could probably be expanded and rephrased more generally as something like, "why do high complexity models (e.g., those with larger numbers of adjustable parameters, such as higher degree polynomials) usually exhibit higher variance?"  I've never seen a really good general mathematical explanation or "proof" of this phenomenon; most discussion that I've encountered seems to treat it as more of an empirical observation rather than theoretical result.  This may be related to the fact that the field has a huge number of alternative model selection criteria currently in widespread use, indicating that there isn't really a good theoretical consensus yet on how to properly quantify a model's complexity (nor the closely related bias-variance trade-off) in the first place.
With those preliminary remarks out of the way, what do we mean by the term "overfitting" in the context of a statistical regression against a high degree polynomial?  I submit that it may be broken down into two separate phenomena:

*

*As the degree of the polynomial in the regression increases, the resulting curve fits each data point ever more closely

*Simultaneously, the oscillations between the data points become more numerous and larger in amplitude--meaning that if we later acquire new data points by making additional observations, those new data points will generally tend to fit the previously fitted curves more and more poorly, as the degree of the fit model increases

How might we explain the above two points, mathematically?
For a polynomial curve fit model specifically, one way of understanding bullet one mathematically is to observe that a lower degree polynomial model is just a higher degree polynomial in which most of the fitted coefficients have been artificially pre-constrained to be equal to zero.  For example, if we model a data set using a polynomial of degree 2:
$$
f(x) = C_{0} + C_{1} x + C_{2} x^{2}
$$
adjusting the parameters $C_{0}, C_{1}, C_{2}$ to result in the best possible fit, it is in some sense equivalent to saying that we have modeled it using a degree 100 polynomial:
$$
f(x) = C_{0} + C_{1} x + C_{2} x^{2} + C_{3} x^{3} + ... + C_{100} x^{100}
$$
in which the upper 98 coefficents have all been constrained such that:
$$
C_{3} = C_{4} = ... C_{99} = C_{100} = 0
$$
Now, looking at it from that point of view, consider what happens when we add a degree to our polynomial model function, increasing from degree 2 to degree 3: it is in some sense like removing a constraint: where $C_{3}$ was previously forced to be $0$, now the fit algorithm is free to adjust it.  Essentially, adding another adjustable parameter gives the fitter an additional "knob" that can be adjusted to allow the resulting fitted curve to more closely track the underlying data.  Another crucial point to bear in mind: by removing the constraint that $C_{3} = 0$, you also allow the fitter a slightly wider latitude to adjust the original three parameters, $C_{0}, C_{1}, C_{2}$, across a wider range of potential values, because there are more opportunities for additional new adjustments to compensate or partially cancel each other out.  So bullet number one (the phenomenon that higher degree polynomials more closely reproduce the observation data) may be understood as a simple result of the fact that fewer constraints and more "adjustable dials" means more opportunities to adjust all the parameter values "just so" in order to obtain a nearly perfect fit to the existing observed data.
But what about all of the wildly gyrating oscillations--how may we understand those, mathematically?  Well, consider that the locations of the local minima / maxima for a polynomial of degree $n$ are obtained by solving the following equation:
$$
\frac{d f(x)}{dx} = 0
$$
or equivalently,
$$
C_{1} + 2 C_{2} x + 3 C_{3} x^{2} + ... + n C_{n} x^{n-1} = 0
$$
In general, a polynomial of degree $n$ may have up to $n-1$ unique local minima / maxima, because the algebraic equation above may have up to $n-1$ unique, real-valued (i.e., non-complex) roots.  Because each local minimum / maximum in this case corresponds to one oscillation, it means that in the absence of a regularization constraint, increasing the degree of the fit model will naturally tend to increase the number of oscillations.
Additional bonus comment: if you program in R, there is an example script published with this question which facilitates experimenting with the practical effects of fitting real data with very high degree polynomials.
