i'm trying to understand the derivation of the Einstein Field Equations through the variational principle.
I have found on most sources online about the Einstein-Hilbert action, that the variation of the Riemann tensor is
$ \delta R^{\rho}_{\sigma\mu\nu} = \nabla_{\mu}(\delta\Gamma^{\rho}_{\nu\sigma}) - \nabla_{\nu}(\delta\Gamma^{\rho}_{\mu\sigma}) $
and we justify taking the covariant derivative of $\delta\Gamma^{\rho}_{\nu\sigma}$ by saying that the difference of connections is a tensor.
How exactly is $\delta \Gamma^{\rho}_{\nu\sigma}$ a difference of connections?
My understanding of the variation is
$ \delta\Gamma^{\rho}_{\mu\nu} := \lim_{\epsilon \to 0} \frac{\Gamma^{\rho}_{\mu\nu}(g^{\alpha\beta}_\epsilon)-\Gamma^{\rho}_{\mu\nu}(g^{\alpha\beta})}{\epsilon}$
where $g^{\alpha\beta}_\epsilon$ is a curve $ \mathbb{R} \rightarrow \Gamma(S^0_2(M))$ in the space of the metrics on $M$, where the connection is defined.
So isn't $\delta\Gamma^{\rho}_{\mu\nu}$ a $\textit{limit} $ of a difference of connections? Does that still count?
Let $\epsilon\mapsto g_{\epsilon}$ be a smoothly-varying $1$-parameter family of metrics. Let $B_{\epsilon}$ denote the $(1,2)$ tensor field corresponding to the difference of connections with metrics $g_{\epsilon}$ and $g_0$, i.e $(B_{\epsilon})^{a}_{bc}=\Gamma^{a}_{bc}(g_{\epsilon})-\Gamma^a_{bc}(g_0)$. You can easily prove that $\epsilon\mapsto B_{\epsilon}$ is a smoothly varying 1-parameter family of $(1,2)$ tensor fields, and that $B_0=0$ is the zero tensor field. Therefore, the derivative $\frac{dB_{\epsilon}}{d\epsilon}\bigg|_{\epsilon=0}$ exists (by the smoothness) and is again a $(1,2)$ tensor field; and this derivative is exactly the limit you're talking about. This tensor field is what we call $\delta\Gamma$, the variation of the Christoffel symbols/connection. Thus, the (smooth) variation of Christoffel symbols really is a smooth $(1,2)$ tensor field.
In general, you're right to be worried that a limit of a certain type of object may fail to be that type of object (e.g a limit of rational numbers could be an irrational number), but the thing is that in our case, we have sufficient smoothness to rule out such a possibility.
One remark: it is rather unwieldy to define the appropriate topologies on the spaces of metrics/ tensor fields so talking about smoothness with respect to $\epsilon$ may seem weird. However, for our purposes, let us just say that $\epsilon\mapsto \zeta_{\epsilon}$, defined for $\epsilon$ in an interval $I\subset \Bbb{R}$, is a smoothly varying 1-parameter family of $(r,s)$ tensor fields if when written with respect to some (hence every) chart $(U,x)$ the component functions $(\zeta_{\epsilon})^{i_1\dots i_r}_{j_1\dots j_s}$ are smooth functions on $I\times U$.