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)

It's much worse than just overfitting.
The problems with polynomials don't become clear in examples with only 10 or 20 parameters, so I'll examine a function that we want to fit with 200,000 parameters, where 200,000 parameters really is the correct number - we aren't overfitting.
Our function is a 10 second audio clip, with pressure sampled 20,000 times per second. If we use linear interpolation, then we can interpret our function as a linear combination of 200,000 little bumps. All these bumps behave basically the same, and if we need more parameters we can just make more. We are dealing with lots of parameters, but adding 10 more parameters isn't that much harder than adding the first 10. A 5 cent computer embedded in a happy meal toy can handle this to play back a cow mooing.
Another classic approach would be to represent our function as a linear combination of sine waves. Finding 200,000 reasonable sine waves is easy enough, and using a fourier transform we can do this computation on any phone processor, without much fuss. sin(2x) isn't really that different from sin(200,000x.)
Now, if we want to model our audio clip as a polynomial, then we have to model it as a linear combination of 200,000 monomials. While finding 200,000 well behaved bumps was easy, and finding 200,000 chill sine waves was a breeze, finding 200,000 well behaved monomials is much harder! $1$ and $x$ were easy enough to work with. We start to run into trouble with $x^{20}$, but with some careful work we could wrangle it into representing audio. By the time we are dealing with $x^{1000}$ we are pulling our hair out- when we put in 0, 1, 2, it puts out 0, 1, 1.071509e+301, which barely even fits in a 64 bit floating point number. $x^{20,000}$ is not a reasonable function to work with, and yet when we resort to recruiting it, we have 180,000 functions still to go. $x^{200,000}$ is right out.
tldr: When representing a complicated thing as a linear combination of building blocks, you need a large collection of reasonable building blocks. There aren't very many reasonable monomials.
