Is there a sense in which one can "Taylor Series'' a distribution? For example, is it true that $$ \Theta(x - \epsilon) = \Theta(x) - \epsilon \delta(x) + \frac{\epsilon^2}{2} \delta'(x) + \ldots $$ I have the feeling the above is not well-defined. However, the above expression seems to make sense under an integral sign.
To see why, let $f(x)$ be some smooth test function defined on the real line, and let $F$ be a primitive of $f$ such that $F'(x) = f(x)$. Suppose we integrate this distribution with $f$ for $a < 0 < \epsilon < b$, then we have $$ \int_a^b dx\ f(x) \Theta(x - \epsilon) = \int_\epsilon^b f(x) = F(b) - F(\epsilon) $$
Using the above "Taylor Series" we have at the same time $$ \int_a^b dx\ f(x) \bigg[ \Theta(x) - \epsilon \delta(x) + \frac{\epsilon^2}{2} \delta'(x) \bigg] = [ F(b) - F(0) ] - \epsilon f(0) - \frac{\epsilon^2}{2} f'(0) - \ldots $$ But $- F(0) - \epsilon f(0) - \frac{\epsilon^2}{2} f'(0) - \ldots = - F(\epsilon)$ and so we seem to have equality here.
Is there a way in which the above makes sense?
I believe so. Let $\Phi\in\mathcal{D}'(\Bbb{R})$ be any distribution on $\Bbb{R}$, and $f\in \mathcal{D}(\Bbb{R})$ any test function. Note that by definition, to translate the distribution by $\epsilon$, we define its action as \begin{align} \langle \tau_{-\epsilon}\Phi,f\rangle&:=\langle \Phi,\tau_{\epsilon}f\rangle \end{align} where $(\tau_{\epsilon}f)(x):=f(x-\epsilon)$. Now, we can use the usual Taylor's theorem here to write for any $k\in\Bbb{N}$, a $k^{th}$ order Taylor expansion \begin{align} f(x-\epsilon)&=\sum_{j=0}^j(-1)^j\frac{f^{(j)}(x)}{j!}\epsilon^j + R_{k,f}(x;\epsilon) \end{align} So, plugging this above, we see that \begin{align} \langle \tau_{-\epsilon}\Phi,f\rangle&=\left\langle \Phi, \sum_{j=0}^j(-1)^j\frac{f^{(j)}(\cdot)}{j!}\epsilon^j + R_{k,f}(\cdot,\epsilon)\right\rangle\\ &=\left\langle \sum_{j=0}^k\frac{\Phi^{(j)}}{j!}\epsilon^j,f\right\rangle + \langle \Phi, R_{k,f}(\cdot,\epsilon)\rangle \end{align} In other words, we can "Taylor expand" the distribution $\tau_{-\epsilon}\Phi$ (i.e $\Phi(x+\epsilon)$ in more convenient but sloppy notation) just as for functions: \begin{align} \tau_{-\epsilon}\Phi &=\sum_{j=0}^k\frac{\Phi^{(j)}}{j!}\epsilon^j + \rho_{k,\Phi}(\epsilon), \end{align} where $\rho_{k,\Phi}(\epsilon)$ is the "remainder term" of the "Taylor Polynomial" of $\Phi$, whose action on a test function $f$ is $\langle \rho_{k,\Phi}(\epsilon),f\rangle= \langle \Phi,R_{k,f}(\cdot, \epsilon)\rangle$.
To make this really mimic the standard Taylor's theorem, it would be nice to know that the remainder term is "small" to the $k^{th}$ order in $\epsilon$, i.e \begin{align} \lim\limits_{\epsilon\to 0}\frac{\rho_{k,\Phi}(\epsilon)}{\epsilon^k}&=0, \end{align} where the limit is to be understood in the (weak* topology) of the space of distributions $\mathcal{D}'(\Bbb{R})$. This is indeed the case, because of the Taylor's theorem (say with Lagrange remainder) for single-variable functions (I quickly wrote down a proof before writing my answer, and I believe it's right, but I'd highly suggest you try to do it yourself, and maybe afterwards, if need be I'll write up my proof, which is really just definition chasing and applying the usual Taylor's remainder formulae).
Of course, we can now apply this in particular for $\Phi = \Theta$, and since its derivatives are given by $\delta$ and its derivatives, the answer to your question is yes, provided you interpret everything appropriately.