The following instance of your identity is well known in physics (and is sometimes called an "operator disentangling" identity)
$:\exp\left[\left(e^W-1\right)_{ij}a_i^\dagger a_j\right]:\;=\exp\left(W_{ij}a^\dagger_i a_j\right)$
where $::$ denotes normal ordering, $a_i$ and $a^\dagger_j$ are canonical Bose annhilation and creation operators satisfying $[a_i,a^\dagger_j]=\delta_{ij}$, and W is an arbitrary matrix (summation implied).
For general (rather than just quadratic) $\phi$ the formula is completely new (indeed remarkable) to me.
Frustratingly, it's hard to track down the origins of the above formula. Here's a recent discussion that includes the above version for a single boson mode (Eq. 30):
Combinatorics and Boson normal ordering: A gentle introduction American Journal of Physics, 75 (7), pp. 639 (2007)
The authors' comments after Eq. 30 seem to imply that the formula doesn't generalize simply.
EDIT: I realized that Sylvester initially states the quadratic form above, and then limits his generalization to functions $\phi$ "linear quantic in $\delta_x$, $\delta_y$, $\delta_z$,...". Still, this generalization appears to contradict Eq. 31 of the above article.
