docs: def _sigma_2(gam, eps):
          """
          gam and eps in units of m_e c^2
          Eq. A3 of Baring et al. (1999)
          Return in units of cm2 / mec2
          """
          s0 = r0 ** 2 * alpha / (3 * eps) / mec2_unit
  
          s1_1 = 16 * (1 - eps + eps ** 2) * np.log(gam / eps)
          s1_2 = -1 / eps ** 2 + 3 / eps - 4 - 4 * eps - 8 * eps ** 2
          s1_3 = -2 * (1 - 2 * eps) * np.log(1 - 2 * eps)
          s1_4 = 1 / (4 * eps ** 3) - 1 / (2 * eps ** 2) + 3 / eps - 2 + 4 * eps
          s1 = s1_1 + s1_2 + s1_3 * s1_4
  
          s2_1 = 2 / eps
          s2_2 = (4 - 1 / eps + 1 / (4 * eps ** 2)) * np.log(2 * gam)
          s2_3 = -2 + 2 / eps - 5 / (8 * eps ** 2)
          s2 = s2_1 * (s2_2 + s2_3)
  
          return s0 * np.where(eps <= 0.5, s1, s2) * heaviside(gam - eps)