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?