def poles_normalization(driver):
    """
    Return P(k,ell=0) monopole normalization, using the current
    values of `f` and `b1`, and the properly normalized
    no-wiggle linear power spectrum

    Returns
    -------
    callable :
        the P(k,ell=0) function callable
    """
    f = driver.theory.fit_params['f'].value
    b1 = driver.theory.fit_params['b1'].value
    beta = f / b1
    m = driver.theory.model
    power_norm = (m.sigma8_z / m.cosmo.sigma8()) ** 2
    return lambda k: (1.0 + 2.0 / 3 * beta + 1.0 / 5 * beta ** 2
        ) * b1 ** 2 * power_norm * m.power_lin_nw(k)
