What is the commutator of a Poisson bracket and the covariant derivative? Consider a classical vector field $V^\mu$ on a curved background. We make a 3+1 split of coordinates into $t,x^i$, where $x^i$ are coordinates on spatial hypersurfaces $\Sigma$ and $t$ the parameter labeling them.
Now consider a canonically conjugate 
$$\tilde{\pi}_\mu = \frac{\partial \tilde{\mathcal{L}}}{\partial (\partial_0 \phi^\mu)},$$
where $\tilde{\mathcal{L}}$ is the Lagrangian scalar (the Lagrangian density would be $\mathcal{L} = \tilde{\mathcal{L}} \sqrt{-g}$). Then the following Poisson bracket holds
$$\{V^\alpha(x^i,t),\tilde{\pi}_\beta (y^i,t)\} = \frac{1}{\sqrt{d}}\delta^{\alpha}_{\;\beta} \delta^{(3)} (x^i - y^i) $$
where $\sqrt{d}$ is the square root of the determinant of the metric $d_{ij}$ induced on $\Sigma$. 
Now consider the total momentum
$$P_\mu(t) \equiv \int_\Sigma (\tilde{\pi}_\alpha(x^i,t) V^{\alpha}_{\;\; ;\mu}(x^i,t) - \delta^t_\mu \tilde{\mathcal{L}}(x^i,t)) \sqrt{d} \,d^3\! x$$
where $V^\alpha_{\;\;;\mu} = \nabla_\mu V^\alpha$ is the covariant gradient and $d\Sigma = \sqrt{d} \, d^3 \! x$ is a covariant spatial volume element on $\Sigma$. 
I would now like to evaluate brackets such as $\{V^\alpha(y^i,t), P_\mu (t)\},\{P_\mu,P_\nu\}$ or $\{V^{\alpha}_{\;\; ;\mu},P_\nu\}$ to explore this algebra further. The problem is, however, that the Poisson bracket and $\nabla_\mu$ obviously do not commute because if I swap their order in different ways, I seem to get different results in every case. So, how does one get the covariant derivative outside the Poisson bracket?
In other words, I am looking for $A^\alpha_{\; \beta \mu}$ and $B^\alpha_{\; \beta \nu}$ such that
$$\{V^\alpha_{\;\;;\mu}(x^i,t), \tilde \pi_\beta(y^i,t)\} = \nabla^{(x)}_\mu\{V^\alpha (x^i,t), \tilde \pi_\beta(y^i, t)\} + A^\alpha_{\; \beta \mu}$$
$$\{V^\alpha(x^i,t), \tilde  \pi_{\beta;\nu}(y^i,t)\} = \nabla^{(y)}_\nu\{V^\alpha (x^i,t), \tilde  \pi_\beta(y^i, t)\} + B^\alpha_{\; \beta \nu}$$
What are they and how can I find them?
 A: There was an error in the foundational assumption of the question. (I have edited the original question so that the notation agrees and you can find definitions of quantities there if not defined here.)
Much of the difficulties with the computation come from the fact that the "proper" canonically conjugate momentum is actually defined as
$$\pi_\mu = \frac{\partial (\tilde{\mathcal{L}} \sqrt{-g})}{\partial (\partial_0 V^\mu)} = \frac{\partial \mathcal{L}}{\partial (\partial_0 V^\mu)}\,.$$
This momentum $\pi_\mu = \tilde{\pi}_\mu \sqrt{-g}$ is then itself a vector density on $\Sigma$ and the Poisson bracket reads
$$\{V^\alpha(x^i,t), \pi_\beta (y^i,t)\} = \delta^\alpha_\beta \delta^{(3)}(x^{i} - y^i)\,.$$
Note that in the original question the $\{V^\alpha(x^i,t), \tilde \pi_\beta (y^i,t)\}$ bracket is stated incorrectly, because we see that we have 
$$\{V^\alpha(x^i,t), \tilde \pi_\beta (y^i,t)\} =  \frac{1}{\sqrt{-g}}\delta^\alpha_\beta \delta^{(3)}(x^{i} - y^i)\,.$$
Using the properties of the delta function it is then easy to show that the following Poisson bracket can be defined for the functionals of the fields: 
$$
A(t)[\pi,\phi] = \int \mathcal{A}(\pi_\alpha,\pi_{\alpha,i}, V^\alpha, V^\alpha_{\;,i},x^i,t) d^3 x $$ $$
B(t)[\pi,\phi] = \int \mathcal{B}(\pi_\alpha,\pi_{\alpha,i}, V^\alpha, V^\alpha_{\;,i},x^i,t) d^3 x $$ $$
\{A(t),B(t)\} = \int \frac{\delta \mathcal{A}}{\delta V^\alpha} \frac{\delta \mathcal{B}}{\delta \pi_\alpha} - \frac{\delta \mathcal{B}}{\delta V^\alpha} \frac{\delta \mathcal{A}}{\delta \pi_\alpha} d^3 x  $$
where $\mathcal{A},\mathcal{B}$ are densities on $\Sigma$ and $\delta \mathcal{F}/\delta f$ is the variational derivative
\begin{equation}
\frac{\delta \mathcal{F}}{\delta f} = \frac{\partial \mathcal{F}}{\partial f} - \frac{\partial \;}{\partial x^i} \frac{\partial \mathcal{F}}{\partial (f_{,i})}\,,
\end{equation}
where we have assumed that $\mathcal{F}$ is dependent only on $f$ and its first-order gradients (for higher order gradients we get a series of analogous terms of varying sign). This addresses only spatial gradients on $\Sigma$ and not the temporal ones. 
The temporal gradients can be eliminated by $field_{,0} = \{field,\mathcal{H}\}$, some first order gradients will be possible to eliminate by substituting
$$\pi_\mu V^{\mu}_{\;,0} - \mathcal{L} = \mathcal{H}\,.$$
As for the commutators of gradients of fields, with these better variables it is easy to compute
$$\{V^\alpha_{\;;\mu}(x^i,t), \pi_\beta (y^i,t)\} = (\delta^\alpha_\beta \partial_{\mu(x)} + \Gamma^\alpha_{\mu\beta})\delta^{(3)}(x^{i} - y^i)$$
$$\{V^\alpha(x^i,t), \pi_{\beta|\mu} (y^i,t)\} = (\delta^\alpha_\beta \partial_{\mu(y)} - \Gamma^\alpha_{\mu\beta})\delta^{(3)}(x^{i} - y^i)\,,$$
where the symbol $\pi_{\beta|\mu}$ stands for a pseudo-covariant derivative $\pi_{\beta|\mu} = \pi_{\beta,\mu} - \Gamma^\gamma_{\mu\beta} \pi_\gamma$ (remember that $\pi_\mu$ is not a covariant quantity) and of course $\partial_0 \delta^{(3)} = 0$. This derivative proves to be very useful in the computation of many brackets.
