Cubic hypergeometric integrals of motion in affine Gaudin models

We construct cubic Hamiltonians for quantum Gaudin models of affine types $\hat{\mathfrak{sl}}_M$. They are given by hypergeometric integrals of a form we recently conjectured in arXiv:1804.01480. We prove that they commute amongst themselves and with the quadratic Hamiltonians. We prove that their vacuum eigenvalues, and their eigenvalues for one Bethe root, are given by certain hypergeometric functions on a space of affine opers.