docs: def _spectrum(self, photon_energy):
          """
          Compute differential spectrum from pp interactions using Eq.71 and
          Eq.58 of Kelner, S.R., Aharonian, F.A., and Bugayov, V.V., 2006
          PhysRevD 74, 034018 (`arXiv:astro-ph/0606058
          <http://www.arxiv.org/abs/astro-ph/0606058>`_).
  
          Parameters
          ----------
          photon_energy : :class:`~astropy.units.Quantity` instance
              Photon energy array.
          """
  
          outspecene = _validate_ene(photon_energy)
  
          with warnings.catch_warnings():
              warnings.simplefilter("ignore")
              self.nhat = 1.0  # initial value, works for index~2.1
              if np.any(outspecene < self.Etrans) and np.any(
                  outspecene >= self.Etrans
              ):
                  # compute value of nhat so that delta functional matches
                  # accurate calculation at 0.1TeV
                  full = self._calc_specpp_hiE(self.Etrans)
                  delta = self._calc_specpp_loE(self.Etrans)
                  self.nhat *= (full / delta).decompose().value
  
              self.specpp = np.zeros(len(outspecene)) * u.Unit("1/(s TeV)")
  
              for i, Egamma in enumerate(outspecene):
                  if Egamma >= self.Etrans:
                      self.specpp[i] = self._calc_specpp_hiE(Egamma)
                  else:
                      self.specpp[i] = self._calc_specpp_loE(Egamma)
  
          density_factor = (self.nh / (1 * u.Unit("1/cm3"))).decompose().value
  
          return density_factor * self.specpp.to("1/(s eV)")