  We construct a model for cosmic ray acceleration from protostellar accretion
shocks and calculate the resulting cosmic ray ionization rate within
star-forming molecular clouds. We couple a protostar cluster model with an
analytic accretion shock model to calculate the cosmic ray acceleration from
protostellar surfaces. We present the cosmic ray flux spectrum from keV to GeV
energies for a typical low-mass protostar. We find that at the shock surface
the spectrum follows a power-law trend across 6 orders of magnitude in energy.
After attenuation, the spectrum at high energies steepens, while at low
energies it is relatively flat. We calculate the cosmic ray pressure and cosmic
ray ionization rate from relativistic protons at the protostellar surface and
at the edge of the core. We present the cosmic ray ionization rate for
individual protostars as a function of their instantaneous mass and final mass.
The protostellar cosmic ray ionization rate is $\zeta \approx 0.01 - 1$
s$^{-1}$ at the accretion shock surface. However, at the edge of the core, the
cosmic ray ionization rate drops substantially to between $\zeta \approx
10^{-20}$ to $10^{-17}$ s$^{-1}$. There is a large spatial gradient in the
cosmic ray ionization rate, such that inner regions may experience cosmic ray
ionization rates larger than the often assumed fiducial rate, $\zeta =
3\times10^{-17}$ s$^{-1}$. Finally, we calculate the cosmic ray ionization rate
for protostellar clusters over 5 orders of magnitude of cluster size. We find
that clusters with more than approximately 200 protostars produce a higher
cosmic ray ionization rate within their natal cloud than the fiducial galactic
value.
