  We calculate the three bulk viscosity coefficients as arising from the
collisions among phonons in superfluid neutron stars. We use effective field
theory techniques to extract the allowed phonon collisional processes, written
as a function of the equation of state of the system. The solution of the
dynamical evolution of the phonon number density allows us to calculate the
bulk viscosity coefficients as function of the phonon collisional rate and the
phonon dispersion law, which depends on the neutron pairing gap. Our method of
computation is rather general, and could be used for different superfluid
systems, provided they share the same underlying symmetries. We find that the
behavior with temperature of the bulk viscosity coefficients is dominated by
the contributions coming from the collinear regime of the $2\leftrightarrow3$
phonon processes. For typical star radial pulsation frequencies of $\omega \sim
10^{4} s^{-1}$, we obtain that the bulk viscosity coefficients at densities $n
\gtrsim 4n_0$ are within 10% from its static value for $T \lesssim 10^9 K$ and
for the case of strong neutron superfluidity in the core with a maximum value
of the $^3P_2$ gap above 1 MeV, while, otherwise, the static solution is not a
valid approximation to the bulk viscosity coefficients. Compared to previous
results from Urca and modified Urca reactions, we conclude that at $T\sim
10^9$K phonon collisions give the leading contribution to the bulk viscosities
in the core of the neutron stars, except for $n \sim 2n_0$ when the opening of
the Urca processes takes place.
