Substituting for Taylor series So my question is simple:
Why is substitution valid? I mean it seems counter-intuitive to me mainly because of the chain rule.
For example: 
The Taylor series of $e^{x^2}$  is  simply done by substituting $x^2$ wherever $x$ goes in the original sum of the Taylor series for $e^x$.
But if I do the Taylor series manually for the function $e^{x^2}$  I have to apply the chain rule and so I get other terms and it's not immediately obvious to me why it is that the same series comes up.
For what reason is this substitution valid? Thanks.
 A: Intuitively, the taylor series of $f$, denoted by $T_f$, evaluated at $x$ is equal to the function $f$ evaluated at $x$.
$$T_f(x)=f(x)$$
Thus evaluating it at some point $x=g(t)$ yields
$$T_f(g(t))=f(g(t))\tag{1}\label{1}$$
which is what we wanted. We will treat this intuitive conclusion as our hypothesis and prove it in a way that addresses your concerns regrading chain rule.
As a preliminary, note that due to the operator
$$e^{h\partial}=\sum\limits_{n=0}^{\infty}\frac{(h\partial)^n}{n!}\tag{2}\label{eq:1}$$
where $\partial=\frac{d}{dx}$, taylor series expansion is an operator on functions. We will use it to prove the hypothesis
$$e^{h\partial}(f\circ g)=e^{h\partial}(f)\circ g\tag{3}\label{hypo}$$
meaning that the taylor series of $f\circ g$ is equivalent to composing the taylor series of $f$ with the function $g$, namely $T_{f\circ g}(x)=T_f(g(x))$.
Taylor series of compositions is expressed as
$$e^{h\partial}(f\circ g)=\exp(\partial_fe^{h\partial_g})(g,f)\tag{4}\label{eq:2}$$
where $\exp(\partial_fe^{h\partial_g})$ is an operator on pairs of maps $(g,f)$, where $\partial_g$ is applied to $g,$ and $\partial_f$ to $f$.
The proof of $\eqref{eq:2}$ can be found here [Theorem 5.6.] and it contains the answer to your question. The lines of the proof pertaining to your inquiry regrading chain rule are 
$$\exp(\partial_fe^{h\partial_g})=\exp\left(\partial_f\sum\limits_{i=0}^{\infty}\frac{(h\partial_g)^i}{i!}\right)=\prod_{i=1}^{\infty}e^{\partial_f\frac{(h\partial_g)^i}{i!}}\Big(e^{\partial_f}\Big)$$
 $$\implies$$
 $$\exp(\partial_fe^{h\partial_g})(g,f)=\sum\limits_{\forall_n}\frac{h^n}{n!}\sum\limits_{\lambda(n)}n!\prod\limits_{k\cdot l\in\lambda}\Big(\frac{\partial_f\partial_g^l(g)}{l!}\Big)^k\frac{1}{k!}\Big(\Big(e^{\partial_f}\Big)f\Big)$$
 where $\lambda(n)$ stands for the partitions of $n$. Taking into consideration the fact that $e^{\partial_f}(f)$ evaluated at a point $x$ is the same as evaluating $f$ at $x$, we see
$$\exp(\partial_fe^{h\partial_g})(g,f)=\sum\limits_{\forall_n}\frac{h^n}{n!}\sum\limits_{\lambda(n)}n!\prod\limits_{k\cdot l\in\lambda}\Big(\frac{\partial_f\partial_g^l(g)}{l!}\Big)^k\frac{1}{k!}\Big(f(g)\Big)$$
which is percisely the taylor series expansion of a composition (noting Faà di Bruno's formula in the innermost sum, which resolves your issue with higher order chain rule). Thus
$$e^{h\partial}(f\circ g)=e^{h\partial}(f)\circ e^{h\partial}(g)$$
meaning that taylor series of $f\circ g$ is equivalent to composing the taylor series of $f$ with the taylor series of $g$. Again, noting that $e^{\partial_g}(g)$ evaluated at a point $x$ is the same as evaluating $g$ at $x$, we see 
$$e^{h\partial}(f)\circ e^{h\partial}(g)=e^{h\partial}(f)\circ g$$ 
$$\implies$$
$$e^{h\partial}(f\circ g)=e^{h\partial}(f)\circ g$$
meaning that the taylor series of $f\circ g$ is equivalent to composing the taylor series of $f$ with the function $g$, proving \eqref{hypo}, as inquired in the question.
A: Def: A function $f$ is $o(x^n)$ for $x \to 0$ if $$\lim_{x \to 0} \frac{f(x)}{x^n} = 0$$
If $f$ is $o(x^n)$ for $x \to 0$, and if $m < n$, then $f$ is $o(x^m)$ too.
There is a theorem, due to Peano, that if $f(x) = P(x) + o(x^n)$, for $P(x)$ a polynomial of degree $n$, then $P(x)$ is the Taylor polynomial of order $n$ of $f$, around zero. Clearly the definition I gave right at the beginning can be generalized for another center, $x_0$. 
I'll find the Taylor expansion around zero of $\ln (1 + x^2)$, of order $6$. It is known that: $$\ln (1 + x) = x - \frac{x^2}{2} + \frac{x^3}{3} + o(x^3)$$
Using $x^2$ instead of $x$, we get: $$\ln (1 + x^2) = x^2 - \frac{x^4}{2} + \frac{x^6}{3} + o(x^6)$$
Since the remainder is $o(x^6)$, that is the polynomial desired, for Peano's theorem.
Bigger examples might require that some terms are "swallowed" by the remainder $o(x^n)$, if the degree of said terms gets more than $n$. See that $o(x^n)$ is a notation, so things like $o(x^n) + o(x^n) = o(x^n)$ shouldn't bother you. At the end of the day, the $o(x^n)$ you have is not the same you began with. For example, $x^4 + x^5 + o(x^3) = o(x^3)$. 
So, it is not simple as just making substitutions, you always have to pay attention to the remainder, or else you might get a polynomial, that is not the best approximation. 
A: We know that for any real number $x,$
$$e^x=\sum_{n=0}^{\infty}\frac{x^n}{n!}.$$
Now if $x$ is a real number, then so is $x^2!$ It follows that
$$\tag 1 e^{x^2}=\sum_{n=0}^{\infty}\frac{(x^2)^n}{n!} = \sum_{n=0}^{\infty}\frac{x^{2n}}{n!}.$$
We didn't use the word "substitution" or say "let $u$ equal …" anywhere. That's all that is going on here.
How do we know $(1)$ is the Taylor series of $e^{x^2}$ at $0?$ That follows from this theorem: Suppose $f(x)=\sum_{n=0}^{\infty} a_nx^n,$ where the power series converges in $(-r,r)$ for some $r>0.$ Then $f$ is infinitely differentiable in $(-r,r),$ and $a_n = D^nf(0)/n!$ for all $n.$
Now to your comment on this easy approach vs. the long awful slog of computing the Taylor coefficients of $e^{x^2}$ by hand: I'm not sure why you're bothered by this. It is just one of many examples of a problem that has two solutions, one simple and direct, the other long, difficult, and possibly tedious. Happens all the time.
A: If you have $\displaystyle f(x)=\sum_{n=0}^\infty a_n x^n$ for any $x\in (-r,r)$ then you have the same equality for $f(x^2), x\in(-\sqrt{r},\sqrt{r}),$ for $f(\sqrt{x}) ,x\in(-r^2,r^2),$ $\cdots$ It is just a change of the name of the variable. The only important thing is that the new variable belongs to the  domain where equality holds.
Of course, you can get the Taylor series of $f(x^2)$ or $f(\sqrt{x})$ (if it exists, in the second case), but if you know the Taylor series of $f(x)$ is just to make a substitution.
A: We consider a more general scenario of the substitution rule.

Let's say we want to expand a function $f$ around $x=0$ as Taylor series. We assume $f$ is the composition of a function $g$ with the power $(x-x_0)^m, m\in \mathbb{N}$.
  \begin{align*}
  f(x)=g((x-x_0)^m)
  \end{align*}
If the function $g$ can  be expanded in a Taylor series around $x=0$ with radius of convergence $R$,
  \begin{align*}
  g(x)=\sum_{n=0}^\infty\frac{g^{(n)}(0)}{n!}x^n  \qquad\qquad |x|<R\tag{1}
  \end{align*}
  we obtain the Taylor series of $f$ around $x=x_0$ by substituting $x$ with $(x-x_0)^m$,
\begin{align*}
  f(x)=\sum_{n=0}^\infty\frac{g^{(n)}(0)}{n!}(x-x_0)^{nm}\qquad\qquad |x-x_0|< \sqrt[m]{R}\tag{2}
  \end{align*}
Note: The validity of the representation as Taylor series of $g$ in (1) is given for $|x|<R$. The validity of the substitution $x\rightarrow(x-x_0)^m$ is therefore assured for $|x-x_0|^m<R$, resp. by the equivalent condition stated in (2).

If the function $e^{x^{2}}$ is to be expanded around $x=0$, we consider the special case
\begin{align*}
g(x)=e^x=\sum_{n=0}^\infty\frac{x^n}{n!}\qquad\qquad\text{and}\qquad\qquad (x-0)^2
\end{align*}
and obtain
\begin{align*}
  e^{x^{2}}=\sum_{n=0}^\infty\frac{x^{2n}}{n!}\qquad\qquad x\in\mathbb{R}
  \end{align*}
