  A remarkable discovery of NASA's Kepler mission is the wide diversity in the
average densities of planets of similar mass. After gas disk dissipation, fully
formed planets could interact with nearby planetesimals from a remnant
planetesimal disk. These interactions would often lead to planetesimal
accretion due to the relatively high ratio between the planet size and the hill
radius for typical planets. We present calculations using the open-source
stellar evolution toolkit MESA (Modules for Experiments in Stellar
Astrophysics) modified to include the deposition of planetesimals into the H/He
envelopes of sub-Neptunes (~1-20 MEarth). We show that planetesimal accretion
can alter the mass-radius isochrones for these planets. The same initial planet
as a result of the same total accreted planetesimal mass can have up to ~5%
difference in mean densities several Gyr after the last accretion due to
inherent stochasticity of the accretion process. During the phase of rapid
accretion these differences are more dramatic. The additional energy deposition
from the accreted planetesimals increase the ratio between the planet's radius
to that of the core during rapid accretion, which in turn leads to enhanced
loss of atmospheric mass. As a result, the same initial planet can end up with
very different envelope mass fractions. These differences manifest as
differences in mean densities long after accretion stops. These effects are
particularly important for planets initially less massive than ~10 MEarth and
with envelope mass fraction less than ~10%, thought to be the most common type
of planets discovered by Kepler.
