I never encountered negative mass anywhere in physics, and would say it must be due to a mistake (or numerical instability) in the calculations. If it was numerical instability, repeat the calculation with slightly different inputs (at the order of 10*roundoff) and check whether this significantly affects the results. If yes, the diagnosis is correct, and you'd need to look at your algorithms to find out where the culprit is. If no, rounding errors are irrelevant and there must be another mistake.
On the other hand, negative mass squared (imaginary mass, equivalently tachyonic states) mean that the spectrum was computed in an unstable (and unphysical) pseudo-vacuum, which would undergo a rapid phase transition.
This means that the fields must be shifted (e.g., broken symmetry in Higgs) and/or linearly transformed by a Bogoliubov transformation (e.g., superconductivity) before they are Fourier expanded in terms of creation and annihilation operators used to find the spectrum.