def pl_psd(f, Si, fL):
    """
    Return the power spectral density P(f) = (f / yr) ^ {-Si} * Theta(f - fL)
    @param f:   Frequency
    @param Si:  Spectral index
    @param fL:  Low-frequency cut-off
    """
    rv = f**(-Si)
    rv[f < fL] = 0.0
    return rv