docs: def _spectrum(self, photon_energy):
          """Compute intrinsic synchrotron differential spectrum for energies in
          ``photon_energy``
  
          Compute synchrotron for random magnetic field according to
          approximation of Aharonian, Kelner, and Prosekin 2010, PhysRev D 82,
          3002 (`arXiv:1006.1045 <http://arxiv.org/abs/1006.1045>`_).
  
          Parameters
          ----------
          photon_energy : :class:`~astropy.units.Quantity` instance
              Photon energy array.
          """
  
          outspecene = _validate_ene(photon_energy)
  
          from scipy.special import cbrt
  
          def Gtilde(x):
              """
              AKP10 Eq. D7
  
              Factor ~2 performance gain in using cbrt(x)**n vs x**(n/3.)
              Invoking crbt only once reduced time by ~40%
              """
              cb = cbrt(x)
              gt1 = 1.808 * cb / np.sqrt(1 + 3.4 * cb ** 2.0)
              gt2 = 1 + 2.210 * cb ** 2.0 + 0.347 * cb ** 4.0
              gt3 = 1 + 1.353 * cb ** 2.0 + 0.217 * cb ** 4.0
              return gt1 * (gt2 / gt3) * np.exp(-x)
  
          log.debug("calc_sy: Starting synchrotron computation with AKB2010...")
  
          # strip units, ensuring correct conversion
          # astropy units do not convert correctly for gyroradius calculation
          # when using cgs (SI is fine, see
          # https://github.com/astropy/astropy/issues/1687)
          CS1_0 = np.sqrt(3) * e.value ** 3 * self.B.to("G").value
          CS1_1 = (
              2
              * np.pi
              * m_e.cgs.value
              * c.cgs.value ** 2
              * hbar.cgs.value
              * outspecene.to("erg").value
          )
          CS1 = CS1_0 / CS1_1
  
          # Critical energy, erg
          Ec = (
              3
              * e.value
              * hbar.cgs.value
              * self.B.to("G").value
              * self._gam ** 2
          )
          Ec /= 2 * (m_e * c).cgs.value
  
          EgEc = outspecene.to("erg").value / np.vstack(Ec)
          dNdE = CS1 * Gtilde(EgEc)
          # return units
          spec = (
              trapz_loglog(np.vstack(self._nelec) * dNdE, self._gam, axis=0)
              / u.s
              / u.erg
          )
          spec = spec.to("1/(s eV)")
  
          return spec