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)


Is there any mathematical justification as to why (higher degree) polynomial functions overfit the data?

Sure. As others mentioned, particularly stachyra and fblundun, it's about the complexity of the hypothesis class relative to the amount of data you have. A highly complex model will always find a way to explain a small amount of data, regardless of whether that explanation generalizes correctly. A simple model won't be able to fit the training data unless it actually explains the underlying relationship.
Imagine that my data distribution looks like this: for each $x$, I pick $y$ uniformly at random from $\{0,1\}$. No model in the world will be able to predict the next data point. It's totally random.
But suppose you plan to draw $10$ training data points and fit a degree-$10$ polynomial. I can already tell you in advance what will happen: you will be able to fit the data perfectly. In other words, your approach can't tell whether you're going to generalize well or not. It always has zero training error even when the next prediction will be very bad.
Whereas with a degree-$3$ polynomial, you will immediately notice a poor fit and conclude that you will not generalize. A simple model can detect whether or not it is fitting the data. A complex one will always fit the data regardless of how it's generated.

Intuition like this is formalized by VC-dimension, a measure of complexity of hypothesis classes for binary classification (but there are versions for regression as well, psuedo-dimension). The theory promises that for a simple model class, if we draw relatively few data points, then the fit to the training data is representative of the fit to the actual generative model. Whereas for a more complex model class, it may fit the training data significantly better than the test data.
