While the GeV $\gamma$-rays 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$-rays emission in SBGs.