We develop and analyze a surface integral equation (SIE) whose solution pertains to numerical simulations of propagating time-harmonic electromagnetic waves in three-dimensional dielectric media. The formulae to evaluate the far-field pattern and propagation of the electric and magnetic fields in the interior and exterior of a dielectric body, through surface integrals, require the solution of a 2 × 2 system of weakly-singular SIEs for the two unknown electric and magnetic fields at the interface surface of the dielectric body. The SIE is governed by an operator that is of the classical identity plus compact form. The tangential surface currents and normal surface charges of the dielectric model can be easily computed from the surface electric and magnetic fields. A novel feature of that second-kind SIE is that when augmented with two stabilization scalar equations, the corresponding reformulation does not suffer from spurious resonances or low-frequency breakdown. We analyze this novel feature of the SIE in both the Sobolev and continuous function space settings (without any tangential constraints) by proving the uniqueness of the stabilized solution of the SIE, for all positive incident wave frequencies and physically appropriate real/complex dielectric media constants. Using numerical examples, we also demonstrate that without our stabilization constraints, the 2 × 2 second-kind SIE suffers from some spurious resonances and incorporating the two additional scalar equations removes these spurious resonances. At the zero frequency limit, the two stabilization scalar equations are automatically satisfied and the 2 × 2 formulation nicely decouples to a well-posed system at the zero frequency limit. Our analysis to establish the equivalence of the classical dielectric model (governed by the Maxwell equations) and the all-frequency SIE based representation provides a natural decomposition of the surface unknowns into their tangential currents and normal charges. The salient features of the representation (and well-posedness analysis) provide a tool to develop a new class of efficient and high-order algorithms (with numerical analysis) for simulation of the three-dimensional dielectric media model from low- to high-frequencies.