  We present a comprehensive analytic model of relativistic jet propagation in
expanding homologous media (ejecta). This model covers the entire jet evolution
as well as a range of configurations that are relevant to binary neutron star
mergers. These include low and high luminosity jets, unmagnetized and mildly
magnetized jets, time-dependent luminosity jets, and Newtonian and relativistic
head velocities. We also extend the existing solution of jets in a static
medium to power-law density media with index $\alpha<5$. Our model provides
simple analytic formulae (calibrated by 3D simulations) for the jet head
propagation and breakout times. We find that the system evolution has two main
regimes: strong and weak jets. Strong jets start their propagation immediately
within the ejecta. Weak jets are unable to penetrate the ejecta at first, and
breach it only after the ejecta expands significantly, thus their evolution is
independent of the delay between the onset of the ejecta and the jet launching.
After enough time, both strong and weak jets approach a common asymptotic
phase. We find that a necessary, but insufficient, criterion for the breakout
of unmagnetized [weakly magnetized] jets is $E_{j,{\rm iso,tot}} \gtrsim
3[0.4]E_{ej,{\rm tot}}\left({\theta_j}/{0.1{\rm~rad}}\right)^2$, where
$E_{j,{\rm iso,tot}}$ is the jet total isotropic equivalent energy, $\theta_j$
is its opening angle, and $E_{ej,{\rm tot}}$ is the ejecta energy. Applying our
model to short GRBs, we find that there is most likely a large diversity of
ejecta mass, where mass $ \lesssim 10^{-3}~{\rm M}_{\odot} $ (at least along
the poles) is common.
