Abstract Computing the return times of extreme events and assessing the impact of climate change on such return times is fundamental to extreme event attribution studies. However, the rarity of such events in the observational record makes this task a challenging one, even more so for ‘record-shattering’ events that have not been previously observed at all. While climate models could be used to simulate such extremely rare events, such an approach entails a huge computational cost: gathering robust statistics for events with return time of centuries would require a few thousand years of simulation. In this study, we use an innovative tool, rare event algorithm, that allows us to sample numerous extremely rare events at a much lower cost than direct simulations. We employ the algorithm to sample extreme heatwave seasons, corresponding to large anomalies of the seasonal average temperature, in a heatwave hotspot of South Asia using the global climate model Plasim. We show that the algorithm estimates the return levels of extremely rare events with much greater precision than traditional statistical fits. It also enables the computation of various composite statistics, whose accuracy is demonstrated through comparison with a very long control run. In particular, our results reveal that extreme heatwave seasons are associated with an anticyclonic anomaly embedded within a large-scale hemispheric quasi-stationary wave-pattern. Additionally, the algorithm accurately represents the intensity-duration-frequency statistics of sub-seasonal heatwaves, offering insights into both seasonal and sub-seasonal aspects of extreme heatwave seasons. This innovative approach could be used in extreme event attribution studies to better constrain the changes in an event’s probability and intensity with global warming, particularly for events with return times spanning centuries or millennia.