I have a continuous transformation given by $$\vec r_i(t)\to\tilde{\vec r_i}(t)=\frac{1}{s}\vec r_i(s^2t), \ \ \ i=1,2$$ Taking $s=1+\epsilon$ we aim to find the corresponding infinitesimal transformation given by $$\vec r_i(t)\to\tilde{\vec r_i}(t)=\vec r_i(t)+\epsilon\zeta_i(r_i,\dot{\vec r_i},t), \ \ \ i=1,2$$ Solution
Beginning with $$\frac{1}{1+\epsilon}\vec r_i((1+\epsilon)^2t)$$ and approximating, utilising the fact that $\epsilon^2\approx0$ we have $$(1-\epsilon)\vec r_i(t+2\epsilon t)=(1-\epsilon)(\vec r_i(t)+\vec r_i(2\epsilon t))+\mathcal{O}(\epsilon^2)$$ Which gives $$(1-\epsilon)(\vec r_i(t)+2\epsilon t\dot{\vec r_i}(t))+\mathcal{O}(\epsilon^2)$$ Utilising $\epsilon^2\approx0$ once more we obtain $$\vec r_i(t)+\epsilon(2 t\dot{\vec r_i}(t)-\vec r_i(t))+\mathcal{O}(\epsilon^2)$$ finally we have $$\zeta_i=2 t\dot{\vec r_i}(t)-\vec r_i(t)$$ I am unsure of the validity of $\vec r_i(2\epsilon t)=2\epsilon t\dot{\vec r_i}(t)$, would anyone be able to advise me on this?
The infinitesimal transformation just means we're considering $f:\Bbb{R}\to\Bbb{R}^2$ defined as \begin{align} f(\epsilon)&=\frac{1}{1+\epsilon}\mathbf{r}((1+\epsilon)^2t), \end{align} and we want to calculate the derivative $f'(0)$ at the origin; i.e this is what you're calling $\zeta$. From the product/chain/quotient rules, we have \begin{align} f'(\epsilon)&=\left(-\frac{1}{(1+\epsilon)^2}\right)\mathbf{r}((1+\epsilon)^2t)+ \frac{1}{1+\epsilon}\dot{\mathbf{r}}((1+\epsilon)^2t)\cdot (2(1+\epsilon)t) \end{align} Hence, evaluating at $\epsilon=0$ yields \begin{align} f'(0)&=-\mathbf{r}(t)+2t\,\dot{\mathbf{r}}(t) \end{align} In other words, $\zeta(\mathbf{r},\mathbf{v},t)=-\mathbf{r}+2t\mathbf{v}$ is the "infinitesimal transformation".
Also, just to emphasize: we're keeping $t$ (which we like to imagine as the "curve parameter") fixed. We're letting $\epsilon$ (the so called "variation-parameter") vary. Differentiating with respect to $\epsilon$ at $0$ thus describes how (to first order) things are changing when the modify the curve itself. Hence, it is called the infinitesimal (which means first order) transformation.
Note that if you want to arrive at this via a "Taylor expansion", then you're essentially going to be reproving the basic rules of differential calculus. Strictly speaking, we can only use the big-oh $\mathcal{O}(\epsilon^2)$ if we assume the functions are $C^1$ (which is actually an unnecessary assumption); first derivatives deal with little-oh $o(\epsilon)$. Anyway, if you don't care about such technicalities, and would like to see how to arrive at the infinitesimal transformation the "physics way", here it is: \begin{align} f(\epsilon)&=\frac{1}{1+\epsilon}\mathbf{r}((1+\epsilon)^2t)\\ &=(1-\epsilon + \mathcal{O}(\epsilon^2))\cdot \mathbf{r}\left(t+ 2\epsilon t+ \mathcal{O}(\epsilon^2)\right)\tag{i}\\ &=(1-\epsilon + \mathcal{O}(\epsilon^2))\cdot \left[\mathbf{r}(t)+ \left(2\epsilon t+\mathcal{O}(\epsilon^2)\right)\,\dot{\mathbf{r}}(t) + \mathcal{O}\left(\left(2\epsilon t+\mathcal{O}(\epsilon^2)\right)^2\right)\right]\tag{ii}\\ &=(1-\epsilon + \mathcal{O}(\epsilon^2))\cdot \left[\mathbf{r}(t)+2\epsilon t\, \dot{\mathbf{r}}(t) + \mathcal{O}(\epsilon^2)\right]\\ &=\mathbf{r}(t)+ \epsilon\left[-\mathbf{r}(t)+2t\dot{\mathbf{r}}(t)\right]+\mathcal{O}(\epsilon^2) \end{align} hence the thing in square brackets is once again what we found above for $f'(0)$. Note that in going from (i) to (ii), what I did was to define the quantity $\eta=2\epsilon t+ \mathcal{O}(\epsilon^2)$, so that \begin{align} \mathbf{r}(t+2\epsilon t + \mathcal{O}(\epsilon^2))=\mathbf{r}(t+\eta)=\mathbf{r}(t)+\eta \,\dot{\mathbf{r}}(t)+\mathcal{O}(\eta^2). \end{align} Then, I plugged in the definition of $\eta$, simplified the resulting algebra, and used some basic properties of big-oh, such as a function being $\mathcal{O}(\eta^2)$ implying that it is also $\mathcal{O}(\epsilon^2)$ and so on.
You mention
Well, you're right to be suspicious because this is completely false. You can't bring scalars out of a function and replace it with a derivative. Remember that the if a function $g$ is differentiable at a point $a$, then $g(a+\epsilon)=g(a)+\epsilon \cdot g'(a)+o(\epsilon)$ (or if you assume $C^1$ then you can put $\mathcal{O}(\epsilon^2)$ for the remainder). But I hope the above explanation is sufficient to understand where the $2\epsilon t$ is actually coming from.