  We develop path-integral Monte Carlo simulations for a parabolic
two-dimensional (2D) quantum dot containing $N$ interacting electrons in the
presence of Dresselhaus and/or Rashba spin-orbit couplings. Our method solves
in a natural way the spin contamination problem and allows for numerically
exact finite-temperature results at weak spin-orbit coupling. For $N<10$
electrons, we present data for the addition energy, the particle density, and
the total spin $S$ in the Wigner molecule regime of strong Coulomb
interactions. We identify magic numbers at N=3 and N=7 via a peak in the
addition energy. These magic numbers differ both from weak-interaction and
classical predictions, and are stable with respect to (weak) spin-orbit
couplings.
