  We present a calculation of the rare decay modes D^0 -> e^+ e^- gamma and D^0
-> mu^+ mu^- gamma in the framework of Standard Model. For the short distance
part, we have derived QCD corrections to the Wilson coefficients involved,
including C_9. The latter is found to be strongly suppressed by the
corrections, leading to diminished values for the c -> u l^+l^- branching
ratios in the 10^(-10) range. Within SM the exclusive decays are dominated by
long distance effects. Nonresonant contributions are estimated using heavy
quark and chiral symmetries to be at the level of 10%, compared to the
contributions arising from D^0-> V gamma->l^+ l^- gamma, with V=rho, omega,
phi. The total SM branching ratio is predicted to be in the range (1-3) x
10^(-9). We also consider contributions coming from MSSM with and without R
parity conservation. Effects from MSSM are significant only for the R-parity
violating case. Such contributions enhance the branching ratio D^0-> mu^+ mu^-
gamma to 0.5 x 10^(-7), based on appropriately allowed values for C_9 and
C_{10}. This selects D^0->mu^+ mu^- gamma as a possible probe of new physics.
