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)

The ringing is an artifact of using uniformly spaced points, because the lagrange polynomials for this spacing are not tightly concentrated around the points they are trying to fit. E.g. for 11 evenly-spaced points on the interval [0,1], here is the degree 10 polynomial that is used to fit the value at x=0.4 (is zero at all other sample points):

Clearly nudging the value at x=0.4 will wildly change the fitted function at unrelated locations. As the Chebfun package (and Trefethen's other work) shows, this problem disappears when using well chosen sampling points over the interval, based on the roots of Chebyshev polynomials. The lagrange polynomials (i.e. basis functions) in this case naturally look closer to smooth kernels. With this choice of fitting points/polynomial basis the polynomial attempting to fit x=0.4 is in fact maximized at this location.

The tradeoff is that you have to fix your desired fitting interval in advance. As we can see in the example the 'concentration' property drops off very sharply outside of [0,1].
So the wild overfitting isn't really about polynomials: the issue is that on evenly spaced points the mapping between curve-distance over the interval and sample-distance on the sampled points is near-degenerate. While we intuitively think of high-degree polynomials as smooth interpolants, this is not automatically true of high-order polynomials and you have forgotten to inform your fitting procedure about what a normal measure of function-distance looks like. This is completely fixable by switching basis, and once this is done you can do very reliable numerical interpolation with 100+ degree polynomials. Of course, this solution is better suited to working numerically with in-principle continuous quantities, rather than for discrete data where we don't get to choose our sample points.
