Z. Keszthelyi, S. A. Brands, A. de Koter, N. Langer, J. Puls

The cluster R136 in the LMC contains a population of stars in excess of 100 M$_\odot$, including R136a1, the most massive star known. Very Massive Stars (VMSs) play an influential role in feedback processes and may potentially produce exotic supernova types and black holes of tens of solar masses. The evolutionary history and final fate of the three most luminous stars, R136a1, R136a2, and R136a3, has been a puzzling issue. We aim to resolve this by rotating single-star MESA models. We produce interpolated model grids and apply a Markov-Chain Monte Carlo analysis to compare our models with observations. The nature of supernova progenitors strongly depends on mass loss and the AM coupling schemes. We predict no pair-instability and no GRB progenitors from our fiducial model grid at LMC metallicity. The onset of Wolf-Rayet-type mass-loss rates on the main sequence leads to a rapid decrease in stellar mass and luminosity. The mass turnover implies that the evolutionary history can only be inferred if additional constraints are available. We utilise the surface helium abundance, which poses a conundrum: R136a1, the most luminous star, is less enriched in helium than R136a2 and R136a3. We propose that this can be explained if both R136a2 and R136a3 were initially more massive than R136a1. From a rigorous confrontation of our models to spectroscopically-derived observables, we estimate an initial mass of 346$\pm41$ M$_\odot$ for R136a1, and $\gtrsim$500 M$_\odot$ for R136a2 and R136a3. Even though VMSs are only present in the youngest clusters below 2 Myr of age, our study strengthens their role in local and galaxy evolution. At LMC metallicity, they will be observable as helium-enriched massive stars after their drastic mass loss, produced via single-star evolution. If the core collapse leads to a supernova, it will be of Type Ib/c. [abridged]