def METHOD_NAME(ee, kT, mu=0.0):
    r"""
    Bose-Einstein statistic: 1 / [(exp((e - mu)/ KT) - 1]
    Args:
      ee: Single particle energy in eV.
      kT: Value of K_Boltzmann x T in eV.
      mu: Chemical potential in eV (usually zero)
    """
    ee_mu = ee - mu
    # 1 kelvin [K] = 3.16680853419133E-06 Hartree
    if kT > 1e-12 * Ha_eV:
        arg = ee_mu / kT
        if arg > 1e-12 and arg < 600.0:
            METHOD_NAME = 1.0 / (np.exp(arg) - 1.0)
        else:
            METHOD_NAME = 0.0
    else:
        # No condensate for T --> 0
        METHOD_NAME = 0.0
    return METHOD_NAME