Multi-messenger gravitational wave (GW) astronomy has commenced with the detection of the binary neutron star merger GW170817 and its associated electromagnetic counterparts. The almost coincident observation of both signals places an exquisite bound on the GW speed $|c_g/c-1|leq5cdot10^{-16}$. We use this result to probe the nature of dark energy (DE), showing that a large class of scalar-tensor theories and DE models are highly disfavored. As an example we consider the covariant Galileon, a cosmologically viable, well motivated gravity theory which predicts a variable GW speed at low redshift. Our results eliminate any late-universe application of these models, as well as their Horndeski and most of their beyond Horndeski generalizations. Three alternatives (and their combinations) emerge as the only possible scalar-tensor DE models: 1) restricting Horndeskis action to its simplest terms, 2) applying a conformal transformation which preserves the causal structure and 3) compensating the different terms that modify the GW speed (to be robust, the compensation has to be independent on the background on which GWs propagate). Our conclusions extend to any other gravity theory predicting varying $c_g$ such as Einstein-Aether, Hov{r}ava gravity, Generalized Proca, TeVeS and other MOND-like gravities.
Download