We present the first Hubble diagram of superluminous supernovae (SLSNe) out to a redshift of two, together with constraints on the matter density, $\Omega_{\rm M}$, and the dark energy equation-of-state parameter, $w(\equiv p/\rho)$. We build a sample of 20 cosmologically useful SLSNe~I based on light curve and spectroscopy quality cuts. We confirm the robustness of the peak decline SLSN~I standardization relation with a larger dataset and improved fitting techniques than previous works. We then solve the SLSN model based on the above standardisation via minimisation of the $\chi^2$ computed from a covariance matrix which includes statistical and systematic uncertainties. For a spatially flat $\Lambda$CDM cosmological model, we find $\Omega_{\rm M}=0.38^{+0.24}_{-0.19}$, with a rms of 0.27 mag for the residuals of the distance moduli. For a $w_0w_a$CDM cosmological model, the addition of SLSNe~I to a `baseline' measurement consisting of Planck temperature together with type Ia supernovae, results in a small improvement in the constraints of $w_0$ and $w_a$ of 4\%. We present simulations of future surveys with 868 and 492 SLSNe I (depending on the configuration used) and show that such a sample can deliver cosmological constraints in a flat $\Lambda$CDM model with the same precision (considering only statistical uncertainties) as current surveys that use type Ia supernovae, while providing a factor 2-3 improvement in the precision of the constraints on the time variation of dark energy, $w_0$ and $w_a$. This paper represents the proof-of-concept for superluminous supernova cosmology, and demonstrates they can provide an independent test of cosmology in the high-redshift ($z>1$) universe.
Comment: 16 pages, 9 figures. MNRAS accepted