Double covariant derivative in coordinates: Why does this work? Lets take a vector field $X$ on some Riemannian manifold $(\mathcal{M},g)$ with Levi-Civita connection $\nabla$. Then, the components of its covariant derivative in coordinates are
$$\nabla_{\alpha}X^{\beta}=\partial_{\alpha}X^{\beta}+\Gamma_{\alpha\gamma}^{\beta}X^{\gamma}$$
Now, since $\nabla X$ is again tensor field, we can apply a second covariant derivative, which in coordinates yields
$$\nabla_{\gamma}\nabla_{\alpha}X^{\beta}:=\nabla_{\gamma}(\nabla_{\alpha}X^{\beta})=\partial_{\gamma}(\nabla_{\alpha}X^{\beta})-\Gamma_{\gamma\alpha}^{\delta}\nabla_{\delta}X^{\beta}+\Gamma_{\gamma\delta}^{\beta}\nabla_{\alpha}X^{\delta}=...$$
Now, instead of writing it like this, let us formally change the order in which the covariant detivative act, i.e. let us write
$$\nabla_{\gamma}\nabla_{\alpha}X^{\beta}=\nabla_{\gamma}(\nabla_{\alpha}X^{\beta})=\nabla_{\gamma}(\partial_{\alpha}X^{\beta}+\Gamma_{\alpha\gamma}^{\beta}X^{\gamma})=\nabla_{\gamma}(\partial_{\alpha}X^{\beta})+\nabla_{\gamma}(\Gamma_{\alpha\delta}^{\beta}X^{\delta})=...$$
Now, mathematically speaking, the two terms on the right-hand side are ill-defined and do not make sense, since $\nabla$ is an operation acting on tensor fields and neither $\partial_{\alpha}X^{\beta}$ nor $\Gamma_{\alpha\delta}^{\beta}X^{\delta}$ are the components of a tensor field. However, if we treat these two terms as if they were rank (1,1) tensor fields, i.e. elements of $\Gamma^{\infty}(T\mathcal{M}\otimes T^{\ast}\mathcal{M})$, and use the standard formula for the connections, we will find the same and correct result for the components of $\nabla^{2}X$.
Now, this seem to work in general, i.e. for arbitrary rank tensors and an arbitrary amount of covariant derivatives (at least I never have seen a counter example).
My question, or lets say, my curiosity, is:

Why does the second "approach", in which we produce ill-defined terms in the steps in-between, work?

Is there any mathematical reason? For example, maybe one can extend the covariant derivative to more general "objects with indices" in a unique way, such that the steps in between become well-defined. For example, one can extend the covariant derivative to a map acting on tensor densities (i.e. sections of tensor product of a tensor bundle and a density bundle). Maybe there is a similar and more general notion which also includes objects like the partial derivative and the Christoffel symbols.
 A: The notion that Christoffel symbols are not tensor fields is not entirely true, or at least somewhat deceptive. If we choose a set of local coordinates $x^\alpha$, then the corresponding Chrisotoffel symbols $\Gamma^\alpha{}_{\beta\gamma}$ represent a perfectly well-defined local tensor field. Issues only arise when we try to interpret the objects involved as global, coordinate-independent objects, since a different set of coordinates will not generally result in the same tensor field, even on their common domain.
In this particular computation, we are fixing a coordinate system at the start, and finding an expression for $\nabla\nabla X$ in terms of the Christoffel symbols of those coordinates, so we can interpret expressions like $\nabla_\alpha(\Gamma^{\beta}{}_{\gamma\delta})$ as a covariant derivative of a tensor field without issue.
A: As an alternative point of view, let me point out the following: In your second approach, if you keep the term $(\partial_{\alpha}X^{\beta}+\Gamma_{\alpha\gamma}^{\beta}X^{\gamma})$ in brackets and use linearity only once there is no more covariant derivative, everything stays perfectly well-defined:
$$\nabla_{\gamma}\nabla_{\alpha}X^{\beta}=\nabla_{\gamma}(\partial_{\alpha}X^{\beta}+\Gamma_{\alpha\gamma}^{\beta}X^{\gamma})=\\=\partial_{\gamma}(\partial_{\alpha}X^{\beta}+\Gamma_{\alpha\gamma}^{\beta}X^{\gamma})+\Gamma_{\gamma\delta}^{\beta}(\partial_{\alpha}X^{\delta}+\Gamma_{\alpha\gamma}^{\delta}X^{\gamma})-\Gamma_{\gamma\alpha}^{\delta}(\partial_{\delta}X^{\beta}+\Gamma_{\delta\gamma}^{\beta}X^{\gamma})$$
Now, as soon as you arrive here, you can use linearity, since $\partial_{\gamma}$ acts on every $C^{\infty}$-function and the terms with Christoffel symbols can also be clearly separated.
