  We present a search for the flavor-changing neutral current decay B0_s -> phi
mu+ mu- using about 0.45 fb^-1 of data collected in p \bar p collisions at
sqrt{s} =1.96 TeV with the D{\O}detector at the Fermilab Tevatron Collider. We
find an upper limit on the branching ratio of this decay normalized to B0_s ->
J/psi phi of B(B0_s -> phi mu+ mu-)/B(B0_s -> J/psi phi) < 4.4\times 10^{-3} at
the 95% C.L. Using the central value of the world average branching fraction of
B0_s -> J/psi phi, the limit corresponds to B(B0_s -> phi mu+ mu-) < 4.1 \times
10^{-6} at the 95% C.L., the most stringent upper bound to date.
