<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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