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 \leftrightarrow 3 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 4 n _ { 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 ^ { 3 } P _ { 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 2 n _ { 0 } when the opening of the Urca processes takes place .