def symmetrize_monomial(monomial, symmetries):
    """Symmetrize monomial by averaging over all symmetry images under symmetries.

    Parameters
    ----------
    monomial : Model or BlochModel object
        Hamiltonian term to be symmetrized
    symmetries : iterable of PointGroupElement objects
        Symmetries to use for symmetrization. `symmetries` must form a closed group.

    Returns
    -------
    Model or BlochModel object
        Symmetrized term.
    """
    return sum([sym.apply(monomial) for sym in symmetries]) * (1/len(symmetries))