  While the GeV $\gamma$-ray emission of starburst galaxies (SBG) is commonly
thought to arise from hadronic interactions between accelerated cosmic rays and
interstellar gas, the origin of the TeV $\gamma$-ray emission is more
uncertain. One possibility is that a population of pulsar wind nebulae (PWNe)
in these galaxies could be responsible for the TeV $\gamma$-ray emission. In
this work, we first synthesize a PWNe population in the Milky Way, and assessed
their contribution to the $\gamma$-ray emission of the Galaxy, using a
time-dependent model to calculate the evolution of the PWN population. Such
synthetic PWN population can reproduce the flux distribution of identified PWNe
in the Milky Way given a distribution of the initial state of the pulsar
population. We then apply it to starburst galaxies and quantitatively calculate
the spectral energy distribution of all PWNe in the SBG NGC 253 and M82. We
propose that TeV $\gamma$-ray emission in starburst galaxies can be dominated
by PWNe for a wide range of parameter space. The energetic argument requires
that $\eta_e \times v_{\rm SN} > 0.01 {\rm yr}^{-1}$, where $\eta_e$ is the
fraction the spin-down energy going to electrons and $v_{\rm SN}$ is the
supernova rate. By requiring the synchrotron emission flux of all PWNe in the
galaxy not exceeding the hard X-ray measurement of NGC 253, we constrain the
initial magnetic field strength of PWNe to be $< 400 \mu$G. Future observations
at higher energies with LHAASO or next-generation neutrino observatory
IceCube-Gen2 will help us to understand better the origin of the TeV
$\gamma$-ray emission in SBGs.
