The evolution of the solar wind

How has the solar wind evolved to reach what it is today? In this review, I discuss the long-term evolution of the solar wind, including the evolution of observed properties that are intimately linked to the solar wind: rotation, magnetism and activity. Given that we cannot access data from the solar wind 4 billion years ago, this review relies on stellar data, in an effort to better place the Sun and the solar wind in a stellar context. I overview some clever detection methods of winds of solar-like stars, and derive from these an observed evolutionary sequence of solar wind mass-loss rates. I then link these observational properties (including, rotation, magnetism and activity) with stellar wind models. I conclude this review then by discussing implications of the evolution of the solar wind on the evolving Earth and other solar system planets. I argue that studying exoplanetary systems could open up new avenues for progress to be made in our understanding of the evolution of the solar wind.

The big picture: evolution of winds of cool dwarf stars. As the star ages, its rotation and magnetism decrease, causing also a decrease in angular momentum removal. In the images, I highlight some of the areas to which wind-rotation-magnetism interplay is relevant. Image reproduced with permission from Vidotto (2016a) Fig. 2 An overview of mass-loss rates (colour coded) in the cool-star HR diagram. Winds of cool stars evolve from hot ( $ 10 6 K) and tenuous ('hot corona') to cold ( $ 10 4 K) and denser ('no corona'). Evolved low-mass stars that show weak or sporadic signatures of a hot corona are denoted 'warm/hybrid'. The zero-age main sequence is shown by the grey line. Image reproduced with permission from Cranmer and Winebarger (2019), copyright by Annual Reviews have hotter winds (on the order of 10 6 K) and low mass-loss rates, while on the top right corner of the HR diagram cool stars have colder winds that can reach temperatures of 10 4 K and maybe even lower, and high mass-loss rates. There is a relatively smooth transition between these two groups, with stars that show signs of weak/warm coronae belonging to an intermediate 'hybrid' group, and they perhaps have a combination of wind driving mechanisms (thermal and mechanical, Ó Fionnagáin et al. 2021).
To derive an accurate picture of the solar wind evolution, a combined observational-theoretical approach is more suitable. Ideally, we would like to inform models by measuring winds of Sun-like stars at different evolutionary stages, thus sampling the physical properties of winds from young to old suns. In practice, however, this is not a simple task, as observing winds of cool dwarfs is currently very challenging. I start this review with an overview of observations of winds of solar-like stars (Sect. 2). Using these results, I attempt to derive an evolutionary sequence of the solar wind mass-loss rate in Sect. 3. One cannot talk about the solar wind evolution without talking about evolution of the key physical ingredients of the solar wind: rotation, magnetism and activity. Thus, in Sect. 4, I introduce some key observations that have allowed us to infer the evolution of these three parameters. As we will see, data sampling activity, rotation and even magnetism are more abundant in the literature than observational data sampling wind evolution. In Sect. 5, I present the results of some models investigating the solar wind evolution. In Sect. 6, I discuss some implications of solar wind evolution on the past of the solar system, as well as on other extra-solar systems. I finish this review with a summary and a discussion of open questions in the field (Sect. 7).
From now on, I will refer to winds of isolated, main-sequence, Sun-like stars simply as winds or stellar winds. Otherwise, whenever I refer to winds of other types of stars, or at different evolutionary phases, I will specify (e.g., winds of red giants, or winds of pre-main sequence stars).
2 How to detect winds of solar-like stars decelerates from supersonic to subsonic velocities. Further out, the astropause separates the stellar wind and the ISM flow. Surrounding the astrosphere, if the relative motion between the stellar wind and the ISM is supersonic, a bow shock forms. The enhancement in hydrogen density, known as the hydrogen wall, occurs between the astropause and the bow shock. Therefore, after Ly-a photons are emitted from the star, they have a long journey until they reach us. They first cross the stellar wind and the hydrogen wall around the astrosphere. Then comes the long journey through the ISM itself. Finally, they cross the hydrogen wall around the heliosphere and traverse the interplanetary space (solar wind) until they reach us and we can observe them (in space, above the Earth's atmosphere). During each stage of this journey, neutral hydrogen present in these different regions absorb the Ly-a spectral line at different velocities. As a consequence, the original Ly-a line is severely altered. The hydrogen wall, in particular, is key for quantifying the stellar wind. In this region, ISM and stellar wind materials mix together, which can lead to charge-exchange between the ionised stellar wind material and the neutral component of the ISM. As a result, stellar wind particles are neutralised, but they still retain their high velocity and temperature. This hot neutral hydrogen, which causes substantial absorption in the Ly-a line, provides the signature required to indirectly quantify stellar winds.
By modelling these different absorption components, one can then infer the column density of neutral hydrogen in the wall surrounding an astrosphere (Wood and Linsky 1998). With assistance from hydrodynamical models of the astrosphere-ISM interaction, it is then possible to estimate the wind ram pressure at the site of the interaction (Wood et al. 2001). Given that the wind ram pressure is linked to the  Wood (2004) 2.2 Radio free-free emission Denser winds and/or radio flares  and Güdel (1992) 2.3 Exoplanets as probes Evaporating planet Vidotto and Bourrier (2017) 2.4 Prominences in H-a Fast rotation Jardine and Collier Cameron (2019) 2.5 Detection of CMEdominated winds Fast CME associated with a flare Crosley and Osten (2018) 2.6 Propagation of radio emission Point source (planet?) within wind Vidotto and Donati (2017) 2.7 X-ray emission from the stellar wind Hot coronal winds  and Llama et al. (2013) 2.8 Charge-exchange induced X-rays Partially neutral ISM Wargelin and Drake (2001) and Wargelin and Drake (2002) 2.9 Accretion onto white dwarfs Close binary with cool dwarf secondary Debes (2006) wind mass-loss rate as P ram / _ Mu 1 , it is then possible to derive _ M if u 1 is known. In the case of astrospheric measurements, a unique u 1 ¼ 400 km/s has been assumed for all observed stars, allowing the mass-loss rate of the observed star to be inferred. For more details, I refer the reader to the method review in Wood (2004).
The ''astrosphere method'' has been the most successful method for measuring mass-loss rates of winds of solar-like stars, with nearly 20 measurements to date (Wood 2018), from stars ranging from spectral types M to F on the main sequence, and a few more evolved cool stars as well. The reason this method has not provided measurements for a larger sample of stars is that it has a sweet spot for optimal performance. Firstly, the ISM needs to be neutral, or at least partially neutral, for charge-exchange to take place, and thus the astrospheric absorption signature to occur in the Ly-a line. This means that once the ISM becomes fully ionised, beyond 10-15 pc, the method becomes inapplicable (Redfield and Linsky 2008). Secondly, if the stars are beyond $ 10 pc from us, the ISM column density can be too large and the ISM absorption obscures the astrosphere absorption. For these reasons, Ly-a observations of nearly all stars beyond 10 pc have turned out to be non-detections (Wood et al. 2005b). Figure 4 summarises the detections of winds of solar-like stars (filled red circles) using the astrosphere method. This figure shows that the mass-loss rate per unit surface area ( _ M=R 2 H ) varies as a function of X-ray flux F X as _ M=R 2 H / F 1:34 X (shaded line). Because X-ray flux is a measure of stellar activity, F X can be used as a rough Fig. 3 Sketch of the ISM (left) interacting with a stellar wind (right) and giving rise to an astrosphere, which is surrounded by a hydrogen wall and possibly a bow shock. The hydrogen wall, which is an enhancement of hydrogen density, is located between the bow shock and the astropause. Figure adapted from Ó Fionnagáin (2020) proxy for age-stars with relatively large X-ray fluxes tend to be younger than stars with lower X-ray fluxes. For stars with F X J10 6 erg s À1 cm À2 , Wood (2004) proposed the existence of a wind dividing line, beyond which the power-law fit (shaded line) would cease to be valid. For solar-like stars, this X-ray flux of 10 6 erg s À1 cm À2 roughly corresponds to an age of 600 Myr. Would there be something happening at around this age that makes the wind mass-loss rate drop more than 2 orders of magnitude? It has been suggested that the magnetic topology of stars with F X J10 6 erg s À1 cm À2 suffers an abrupt change that could affect their winds (Wood et al. 2005a), an idea that is backed-up by solar wind observations (see Sect. 3 for the solar wind analogy). However, magnetic field reconstructions of young stars do not show abrupt changes in their field topology, but rather a smooth transition from a magnetic field with an important toroidal component at young ages (fast rotation) to a topology that is dominated by a poloidal field at older ages Folsom et al. 2018;See et al. 2015;Vidotto et al. 2016).
p 1 UMa, a young solar-analogue star, is at the centre of the discussion on the existence of the wind dividing line. As can be seen in Fig. 4, the derived astrospheric mass-loss rate is quite low, and yet, the star has a high F X (Wood et al. 2014). p 1 UMa would, therefore, support the idea of the existence of a wind dividing line. However, an alternative explanation is that p 1 UMa might be in a Fig. 4 Summary of mass-loss rates for low-mass stars derived from the astrosphere method. For GK dwarfs (red filled circles), the mass-loss rate per unit surface area varies as a function of X-ray flux as / F 1:34 X (shaded line). It has been suggested that active (and overall younger) stars with F X J10 6 erg s À1 cm À2 would have reduced mass-loss rates, thus giving rise to a 'wind dividing line'. Some new mass-loss rates measurements indicate that mass-loss rates of young Suns could actually remain large (cf. Sect. 3). Image reproduced with permission from Wood (2018), copyright by the author strongly ionised ISM (the star is at a distance of 14 pc) and the measured Ly-a absorption could be taking place in a cloud between the star and the observer instead of in the hydrogen wall (Wood et al. 2014). Indeed, spin down models (see Sect. 5.3) predict that solar analogues with p 1 UMa's rotation period and the age should have a mass loss rate that is about one order of magnitude higher than the present-day solar wind mass-loss rate (Johnstone et al. 2015a). One way to clarify this contradiction is to use multiple wind-detection methods applied to the same star. In fact, radio observations to detect the wind of p 1 UMa have been conducted (Fichtinger et al. 2017). In the next section, I will discuss this method in more details.

Detecting free-free radio emission from winds
What is unfortunate about the astrosphere method is that the non-detections do not allow us to derive upper limits for mass-loss rates. Although a non-detection could indeed be due to very low wind mass-loss rates, it could also be due to a fully ionised ISM around the astrosphere or an ISM that causes too much absorption, rendering the Ly-a astrospheric signature undetectable. On the other hand, the method that I will discuss now, based on thermal radio emission of winds, is a method that can extract meaningful information even in cases of non-detection of stellar winds. The main current disadvantage of this method is that no wind of lowmass star has ever been detected and, thus, all it has provided so far are upper limits on mass-loss rates. With the advance of radio technologies and construction of more sensitive radio telescopes, this could change in the very near future.
The stellar wind is a thermal, ionised plasma and, as such, it emits continuum bremsstrahlung (free-free) radiation across the electromagnetic spectrum. In particular, the densest region in a stellar wind (its innermost region) can emit at radio wavelengths, thus providing a way to directly detect the wind in radio (Güdel 2002). For large enough densities, the innermost regions can become optically thick to radio wavelengths, creating a radio photosphere Kavanagh et al. 2019), that, if detected, can allow us to quantify the mass-loss rate of the wind. In this case, the underlying non-thermal radio emission from the star cannot be seen. On the other hand, a low-density wind would be radio transparent, and radio (non-thermal) emission from the stellar surface can pass through the wind unattenuated. Non-thermal radio emission, such as radio flares, has been detected in M dwarfs and solar-like stars (e.g., Fichtinger et al. 2017). These radio-transparent winds can nonetheless provide important upper limits to wind mass loss rates. Regardless of whether the wind is optically thin or thick, the wind itself produces radio emission and that may provide a direct detection method.
The idea that winds of early-type (hot) stars could emit in radio was proposed in the seminal works of Panagia and Felli (1975), Wright and Barlow (1975) and Olnon (1975), which initially evolved from the theory of continuum emission in HII regions. The main modification to the previous theory was that, contrary to a constant density approximation adopted in HII regions, stellar winds were assumed to have a spherically-symmetric density profile of the form q / r À2 , where r is the radial coordinate (this implicitly assumes a constant wind speed, which, as I will show next, is incorrect in the innermost region of stellar winds). This modification led to a change in the shape of the spectrum. In the asymptotic limit of an optically thick plasma (s m ) 1 at low frequencies m), the flux density of an HII region behaves as S m / m 2 , while for a stellar wind, it was then found that S m / m 0:6 . In the optically thin asymptotic limit (s m ( 1 at high m), the radio flux density of a stellar wind has the same frequency-dependence as that of an HII region: S m / m À0:1 . Lateron, Reynolds (1986) dropped the assumption of spherically symmetric winds with constant speeds and demonstrated that spectral indices could range between À 0:1 and 2. Indeed, in the low-frequency, optically thick regime, wind models predict a range of slopes-Ó Fionnagáin et al. (2019) derive indices from 1.2 to 1.6 in 3D wind simulations.
Although currently there has been no detection of free-free emission originating in winds of solar-type stars, radio observations have provided upper limits of massloss rates for a number of objects (Güdel 1992;Gaidos et al. 2000;Villadsen et al. 2014;Fichtinger et al. 2017;Vidotto and Donati 2017). Even for the pre-main sequence phase, when stellar winds are expected to have higher mass-loss rates, only upper limits have been derived (e.g., 3 Â 10 À9 M yr À1 for the weak-lined TTauri star V830 Tau, Vidotto and Donati 2017), indicating that this is indeed a tricky detection with current instrumentation. For stars identified as good solar analogues, like p 1 UMa, b Com, j Ceti, EK Dra and n 1 Ori, upper limits of mass-loss rates are as low as 3 Â 10 À12 M yr À1 and as high as 1:3 Â 10 À10 M yr À1 (Gaidos et al. 2000;Fichtinger et al. 2017 complication is that, if the wind has a small mass-loss rate, its weak emission competes against the expected thermal free-free radio emission from the chromosphere and gyromagnetic emission from active regions. For such weak detections, the distinction between a wind signature and the thermal component from close to the stellar surface requires extremely accurate radio spectra (Drake et al. 1993;Villadsen et al. 2014). In a number of simulations of winds of solar-like stars, Ó Fionnagáin et al. (2019) predicted that the radio flux density should increase for higher stellar rotation rates X H (Fig. 5). In particular, at m ¼ 1 GHz, the radio flux density S 1 GHz ' 0:7 lJy X H X ! 0:7 10 pc d ! 2 ' 0:7 lJy 4:6 Gyr t where d is the distance, t the age, and X is the present-day solar rotation rate. In the equation above, I used the fact that X H / t À0:56 , which applies to ages J 600 Myr (Delorme et al. 2011, see Sect. 4.2), i.e., younger stars rotate faster. This implies that an analogue of our present-day Sun, when placed a 10 pc, would emit a flux density of $ 0:7 lJy at 1 GHz. 2 This is below detection limits of current radio telescopes. In fact, the distance-square decay in Eq.
(1) plays a more important role than the increase in rotation rate (or decrease in age) of the stars. Young suns, such as j Ceti and p 1 UMa, are also the closest to us in the simulated sample of Ó Fionnagáin et al. (2019)-for these objects, the predicted 1-GHz flux density is still only a few lJy. At the time of writing (Spring 2020), there exist plans to upgrade the existing VLA system, which would increase instrument sensitivity by a factor of 10 ). The expected sub-lJy sensitivity level of the future Square Kilometre Array (SKA) means that winds of close-by young suns (such as j Ceti and p 1 UMa) are potentially detectable below 1 GHz (see also discussion in the previous paragraph). One consequence of the free-free emission of stellar winds is the creation of a radio photosphere, which is paraboloidal in shape (Kavanagh et al. 2019). Radio emission produced from sources inside the photosphere might be absorbed by the stellar wind. An expected source of low-frequency radio emission are magnetised close-in exoplanets (Farrell et al. 1999). If a planet is embedded in the radio photosphere of their host star, then it is possible that most of the planetary radio emission gets absorbed and does not escape Kavanagh and Vidotto 2020). Although this might be problematic for planet detection in radio, close-in exoplanets can be a useful tool to probe the inner regions of stellar winds, as I explain next.

Using exoplanets to probe stellar winds
Close-in, giant planets form the majority of the exoplanet population known nowadays. With orbital distances significantly smaller than Mercury's orbit (and as small as . 0:02 au), these exoplanets are embedded in a stellar wind regime that is unprecedented for solar system planets. For comparison, the innermost planet in the solar system is Mercury, with a semi-major axis of 0.4 au. Although the solar wind speed at Mercury is similar to the wind speed at Earth's orbit (roughly 400 km/s), the local solar wind density at Mercury is about 1=0:4 2 times larger. This is a consequence of the r 2 -decay of density for winds at asymptotic terminal speeds. However, going even closer to the Sun, one reaches the acceleration zone of the solar wind, in which the wind speed is still increasing with distance. Due to mass conservation, densities increase nearly exponentially towards small heliocentric distances.
Not only the local stellar wind densities around close-in planets are orders of magnitude larger than the $ 5 cm À3 at Earth's orbit, the magnetic field embedded in the stellar wind is also expected to be several orders of magnitude larger . Additionally, the close orbital distances imply that close-in planets have high Keplerian velocities (v kep / r À1=2 ). This means that, in the planet's reference frame, the stellar wind particles could still arrive at large velocities, even though the local stellar wind itself might still not have reached terminal speed (Vidotto et al. 2010a(Vidotto et al. , 2011. Altogether, this shows that close-in planets face extreme wind environments . Although an extreme stellar wind environment could be detrimental for life formation Khodachenko et al. 2007;Scalo et al. 2007;Vidotto et al. 2013), it is precisely the extreme wind conditions around close-in planets that amplify signatures of the wind-planet interaction, potentially leading to detection of such signatures. Signatures of wind-planet interaction have been detected in the hot-Jupiters HD209458b, HD189733b and the warm-Neptune GJ436b, that orbit main-sequence stars of spectral types F8, K2 and M3, respectively. These planets show strong atmospheric escape, which can only be interpreted when considering the interaction with the winds of their host stars (e.g., Holmström et al. 2008;Bourrier et al. , 2016Villarreal D'Angelo et al. 2014;Kislyakova et al. 2014). Given that the escaping atmospheres of these gas giants are hydrogen dominated, their outflows are detected in the Ly-a line during planetary transits, in a technique known as transmission spectroscopy or spectroscopic transit observations (for a review on the topic see, e.g., Kreidberg 2018).
By modelling the stellar Ly-a line profile that is transmitted through the planetary atmosphere, one can derive the conditions of the stellar wind surrounding the exoplanet. The stellar wind has two effects that leave their fingerprints in the Ly-a line. Firstly, the wind shapes the outflow of the planet, in a similar way as the ISM shapes the astrosphere. This causes an asymmetry in the planetary outflow, as the interaction occurs preferentially on one side of the planet, causing lightcurve asymmetries (Villarreal D'Angelo et al. 2018a;McCann et al. 2019). Secondly, the ionised wind exchanges charge with the neutral hydrogen escaping the planetary atmosphere, creating a population of high-velocity neutrals (previously, these were high velocity ions from the stellar wind that became neutralised during the chargeexchange process). This leads to a high-velocity, blue-shifted component of the stellar Ly-a line (e.g., Holmström et al. 2008). Therefore, modelling these fingerprints left in the Ly-a line, one can then obtain the local densities and speeds of the stellar wind (Bourrier et al. 2016), and, consequently, the wind mass-loss rates as well (Kislyakova et al. 2014, Vidotto andBourrier 2017).
In the case of GJ436b, a warm Neptune orbiting an M dwarf, models of the spectroscopic transit in Ly-a predict local wind speeds of 85 km/s and proton densities of 2 Â 10 3 cm À3 (Bourrier et al. 2016), which translates to wind mass-loss rates of 1:2 Â 10 À15 M yr À1 (Vidotto and Bourrier 2017). For the solar-like star HD209458, Kislyakova et al. (2014) were able to model the Ly-a line with a wind that resembles that of the present-day Sun: with the same mass-loss rate, but with a scaled-up density at the orbit of HD209458b. This results in local wind speeds of 400 km/s and wind densities of 5 Â 10 3 cm À3 . For HD189733, the derived local stellar wind speed at the orbit of the planet of 190 km/s and density of 3 Â 10 3 cm À3 (Bourrier and Lecavelier des Etangs 2013) yielded a mass-loss rate of approximately 4 Â 10 À15 M yr À1 . One caveat to keep in mind, though, is that the values derived above are likely not unique, and some degeneracy might exist (Mesquita and Vidotto 2020).
Although the wind properties can be derived as a by-product of a method to detect escape in exoplanets, there are disadvantages to conduct transmission spectroscopy in the Ly-a line. The biggest of which is that the line falls in the ultraviolet part of the electromagnetic spectrum, whose observations require space instrumentation, which are very expensive. Currently, only the Hubble Space Telescope is able to observe in the ultraviolet. For this reason, the recent detection of escaping atmospheres in lines that can be detected with ground-based instrumentation, such as H-a or the HeI triplet at 10,830 Å Nortmann et al. 2018;Spake et al. 2018;Allart et al. 2019), can present new opportunities for probing escaping atmospheres, and likely the interaction region with stellar winds (Oklopčić and Hirata 2018;Villarreal D'Angelo et al. 2021).
The typical wind conditions around exoplanets are soon to become better understood, with experiments that are taking place, right now, in our own solar system. Until very recently, no in-situ measurements of the solar wind plasma at close distances to the Sun existed. Measurements of the first encounter of the NASA's Parker Solar Probe (PSP) were recently released, probing the solar wind at distances as close as 36 R ' 0:17 au (Bale et al. 2019). PSP is set to remain in the heliospheric equator and, after a series of Venus flybys, it will get closer and closer to the Sun, reaching a highly elliptical orbit with a perihelion of 0.046 au-typical of hot Jupiters. A second complementary spacecraft, ESA's Solar Orbiter, was launched in February 2020 and will also study the solar wind at close distances, down to 0.29 au, albeit at high heliospheric latitudes (polar regions). Together, these two spacecrafts will allow in situ measurements of the solar wind at unprecedented close heliospheric distances. They will provide more information about the physical mechanism that accelerates the solar wind and will provide information about the environment at the orbits of close-in exoplanets.

Using prominences to probe stellar winds
The last method for detecting winds of solar-like stars that I would like to discuss involves using observations of slingshot prominences to derive wind mass-loss rates. Slingshot prominences occur on fast rotating stars, and are very extended (Jardine and van Ballegooijen 2005), in contrast to solar prominences. They are detected as absorption features that travel in the H-a stellar line profile as the star rotates. Their observed velocities indicate that slingshot prominences occur at or beyond the co-rotation radius, which are several radii above the stellar surface (Collier Cameron and Robinson 1989a, b). Jardine and van Ballegooijen (2005) suggested that prominences are formed at the top of long magnetic loops, which are filled with mass from the stellar wind. The material in the prominence cools to a few 10 4 K-thus, if the material has enough optical depth, slingshot prominences are seen in H-a Doppler maps in absorption, when they pass in front of the stellar disc.
The dynamical support and lifetime of slingshot prominences depend on the relative position between where they are formed (i.e., the co-rotation radius) and the Alfvén and sonic surfaces of a stellar wind (Fig. 6). The Alfvén and sonic surfaces are defined as the surface where the wind speed reaches the Alfvén and sound speeds, respectively (more about this will be discussed in the modelling Sect. 5). Two main conditions should be kept in mind with regards to the existence of slingshot prominences. Firstly, given that these prominences occur at loop tops, they can only exist when the co-rotation radius lies below the Alfvén surface, as beyond the Alfvén surface, magnetic field lines are open (Jardine and Collier Cameron 2019). Villarreal D'Angelo et al. (2018b) demonstrated that for older solar-like stars, the co-rotation radius is above the Alfvén radius, and these stars cannot support these types of prominences. Young, fast rotating stars, on the other hand, belong to the other group, in which prominences occur within the Alfvén surface.
The second important condition to consider is the relative location between a prominence (co-rotation point) and the sonic point. Information about what is happening in the wind cannot be passed back to the star for any event that takes place above the sonic point (Del Zanna et al. 1998;Vidotto and Cleary 2020). Thus, if the prominence, formed at the co-rotation radius, is formed above the sonic point, the star keeps loading the prominence with stellar wind material and the loop top becomes denser and denser, until it eventually erupts, and the cycle starts again. This represents the 'limit-cycle regime' proposed in Jardine and Collier Cameron (2019). If the site of prominence formation at the co-rotation radius, on the other hand, occurs below the sonic point, the mass-loading from the surface gets readjusted-although prominences could still be formed in the 'hydrostatic regime', they would erupt on an occasional basis only.
Altogether, these conditions imply that the formation of slingshot prominence occur in the 'limit-cycle regime' (Fig. 6). Such a condition is more easily met in faster rotators, as these stars have smaller co-rotation radii and thus more easily to be formed in the sub-Alfvénic regime. Their winds need to be relatively hotter, so that the sonic point occurs at lower heights. As I will discuss later on, we expect that active stars have hotter winds. These conditions are more easily met in fast and ultrafast rotators, in agreement with extended prominences seen in the ultrafast rotators such as Speedy Mic (Jeffries 1993;Dunstone et al. 2006), HQ Lup (Donati et al. 2000) and, more recently, in V530 Per (Cang et al. 2020), all of which with rotation periods \0:4 d. Slingshot prominences have been predicted to occur in stars that are rotating at less extreme rates as well, such as in the young Sun HIP 12545 (Villarreal D'Angelo et al. 2018b), which shows a rotation period of nearly 5 days.
From H-a observations, one can measure the prominence lifetime and the amount of mass contained in the prominence, giving the rate of mass that upflows into the prominence. To convert this mass-loading rate, which occurs in a localised region of the stellar surface, to a wind mass-loss rate, the prominence surface coverage is needed. Jardine and Collier Cameron (2019) estimated that about 1% of the surface of a fast rotating star would be covered in prominences. These authors then derived mass-loss rates of 350, 4500, and 130 times the solar-wind mass loss rate for the young solar-type stars AB Dor, LQ Lup and Speedy Mic, respectively, all of which rotate with a period shorter than half a day. 3 Due to their youth, these three solar-type stars are X-ray luminous. When placing them together in the diagram shown in Fig. 4, we notice that the trend of mass-loss rates with X-ray flux extends for a larger dynamical range in X-ray fluxes, compared Fig. 6 The formation of slingshot prominences occurs when in the 'limit-cycle regime'. In this case, the site of prominence formation, i.e., on loop-tops at the corotation radius, occurs above the sonic point (the star is ''unaware'' that the prominence has formed and thus keeps loading it with stellar wind material) and below the Alfvén radius (beyond the Alfvén radius all magnetic field lines will be open). By observing slingshot prominences, one can estimate the rate at which mass is loaded into the loop tops and thus derive mass-loss rates of stellar winds. Image reproduced with permission from Jardine and Collier Cameron (2019), copyright by the authors to the work of Wood (2018). In Sect. 3, I will discuss how the outcomes of these detection methods can be used to derive an observed evolutionary sequence for mass-loss rates.

Detecting coronal mass ejections (CMEs) through type II radio bursts
Although not all solar CMEs have a flare counterpart (and vice-versa), solar CMEs are often associated to flares. Aarnio et al. (2011) showed that stronger flares lead to more massive CMEs and that the CME mass scales with the X-ray flare flux as M CME / F 0:7 flare or, in terms of the flare energy, as M CME / E 0:68 flare (Aarnio et al. 2012). If one were to extend these solar empirical relations to active stars, which show higher flare energies, one would expect that the mass contained in stellar CMEs could become substantially large. As a result, active stars, due to their higher flare rates, could have CME-dominate winds with mass-loss rates that could be substantially larger than solar.
In the Sun, CMEs contribute, on average, to a mass-loss rate of ' 4 Â 10 À16 M yr À1 (Vourlidas et al. 2010), which is about only a few percent of the present-day solar wind mass-loss rate. In contrast are the predictions for young stars-using flare-CME empirical scalings, Aarnio et al. (2012), Drake et al. (2013) and Osten and Wolk (2015) estimated CME mass-loss rates for T Tauri stars in the range $ 10 À12 -10 À9 M yr À1 , which are several orders of magnitude larger than the present-day solar wind. However, for CME-dominated winds with rates J 10 À10 M yr À1 , Drake et al. (2013) argued that they would have kinetic energies that could amount to 10% of the stellar bolometric luminosity, which is a rather substantial energy budget associated to CMEs. These authors questioned that instead the solar flare-CME relation could not be extrapolated indefinitely to higher flare energies. For example, the relationship might flatten out towards higher flare energies or that observed stellar flares might not always be accompanied by a CME. This could happen if, for example, CMEs on active stars were more strongly confined, and so fewer of them would be produced for any given number of X-ray flares.
This stronger confinement could be caused by the noticeable differences between the solar and stellar magnetic field characteristics. Indeed, more active stars seem to have more toroidal large-scale magnetic field topologies See et al. 2015;Vidotto et al. 2016), which could indeed result in more confined CMEs. Using numerical simulations, Alvarado-Gómez et al. (2018) showed that a stronger overlying large-scale dipolar (poloidal) magnetic field of 75 G could prevent a typical solar CME from erupting and that only CMEs with magnetic energies J 30 times larger than those in typical solar CMEs were able to escape. Those that escaped, however, were not able to accelerate efficiently due to the strong overlying magnetic field.
While Aarnio et al. (2011) and Drake et al. (2013) concentrated on flare energies in the X-ray band, Osten and Wolk (2015) investigated empirical solar flare-CME relations using the bolometric energy from flares. These authors provided energy partition calculations that can be used to relate the amount of radiated flare energy in one bandpass (e.g., in white light or X-rays) to the total bolometric energy of flares. White-light flares, for example, are frequently seen in high-cadence observations of missions such as Kepler, K2 and TESS. The long-term monitoring provided by these missions allows us to build more complete statistical studies of flares from stars of different spectral types and ages (Maehara et al. 2012(Maehara et al. , 2015Davenport et al. 2014), which can then be used to study stellar CMEs.
An alternative approach to estimate mass-loss rates of CME-dominated winds was proposed by Cranmer (2017). He used relations between solar magnetic energy flux and the kinetic energy flux of the solar wind/CME outflows to predict the evolution of CME mass-loss rates of solar-like stars. His models predict that both the wind and CME mass-loss rates are larger at younger ages, with CMEs dominating the mass loss process in the first 0.3 Gyr of the lifetime of a solar-like star. At these younger ages, Cranmer (2017) estimated that CME mass-loss rates are a factor of 10-100 higher than the quiescent (overlying) stellar wind. At ages J 1 Gyr, the mass loss in the quiescent wind dominates that in CMEs.
Although these models qualitatively agree that CMEs can contribute significantly to the total mass-loss rates at younger ages, observationally confirming the presence and amount of stellar CMEs is challenging and their detection has remained elusive (see, e.g., Leitzinger et al. 2020, for an updated census). Given that type II radio bursts are related to CMEs in the Sun (not all CMEs produce type II bursts though), one possible way to detect stellar CMEs is through observations of radio bursts (see, e.g., Crosley et al. 2016, and references therein). Solar type II radio bursts are believed to originate in the shocked material created as the CME propagates outwards at super-Alfvénic velocities (Gopalswamy et al. 2019). The plasma emission from this shocked material is seen in the dynamic spectrum as a burst drifting in time (related to the speed of the shock/CME) and frequency (related to the density of the corona). Given that the coronal density decreases with height, type II radio bursts drifts towards lower frequencies as the CME propagates outwardsstarting (maximum) frequencies of $ 100 MHz corresponds to electron densities of $ 10 8 cm À3 (see Eq. 5 that will be discussed below).
Not all CMEs should produce type II bursts though. For example, if stellar CMEs are not sufficiently accelerated to become super-Alfvénic (see, e.g., slowly accelerated CMEs as seen in the models from Alvarado-Gómez et al. 2018), a shock wave would not be formed and thus a type II burst would not be seen. In spite of this, given that type II bursts are indicative of solar CMEs, by analogy, detecting stellar type II radio bursts could provide a lower limit on stellar CMEs (Crosley et al. 2016), which can be used to determine CME occurrence rates and characterise their properties. Crosley et al. (2016) relates the shock speed v CME to the frequency drift rate _ f and coronal scale height H as where f is the frequency. Although H has to be modelled, both _ f and f are quantities that are directly obtained from type II radio burst observations. With this, one could derive the CME speed. One further step is required to convert the CME speed into CME mass. For that, it is assumed that energy equipartition between the bolometric radiated energy of the flare (E flare; bol ) and the kinetic energy of its associated CME holds. Thus where Eq.
(2) was used in the last equality of Eq. (3). E flare; bol is not a directly observed quantity, but energy partition calculations can be used to convert from observed flare energies at a given bandpass to E flare; bol (Osten and Wolk 2015). Thus, while no type II burst has yet been observed associated with a stellar flare, if one were detected, the rate of frequency change would provide a very useful measure of CME speed (Eq. 2) and mass (Eq. 3, cf. Crosley and Osten 2018). With an observationally derived flare frequency distribution, this could then be used to estimate the total mass loss in CME-dominated winds. More recently, the prospects of using coronal dimming have been suggested as a means to identify stellar CMEs (Harra et al. 2016;Jin et al. 2020). In the Sun, coronal dimming is seen in certain EUV coronal lines, tracing the post-CME evacuation of the corona. From the depth of the light curve, one can infer how much mass has been blown out, and from the slope of the light curve, one can determine the CME speed (Mason et al. 2016). The suggestion of detecting the stellar equivalent of solar coronal dimming is still in its first steps and will be explored with future instrumentations (France et al. 2019). Some points still need to be elucidated, such as whether dimming could be detected in the case where CMEs are happening all the time, as could potentially be the case of young stars, or whether dimming could change the ''basal'' level of the stellar EUV emission to such an extent that one would not be able to disentangle particular CMEs. To answer these and other open questions, more detailed modelling studies are needed (see, e.g., Jin et al. 2020).
2.6 Propagation/suppression of radio emission from a point source (planet) embedded in a stellar wind As discussed in Sect. 2.2, a stellar wind, if sufficiently dense, can emit free-free emission at radio frequencies. We can define a boundary around the star wherein the bulk of the wind emission comes from. Here, we define this boundary, also known as the radio photosphere, as the isosurface where the frequency-dependent optical depth is s m . The value of s m ¼ 0:399, for example, delineates the region within which 50% of the radiation is absorbed by the stellar wind and 50% of the emission escapes (Panagia and Felli 1975). The left panel in Fig. 7 illustrates the radio photosphere at m ¼ 30 MHz (dashed line) and the density profile of the wind of a sun-like star with a mass-loss rate of 2 Â 10 À12 M yr À1 (Kavanagh and Vidotto 2020). In this figure, we are seeing a 2D cut of the radio photosphere, where the observer sees the system from the negative x axis. In three dimensions, the radio photosphere resembles a 'tea cup' (see Fig. 5 in Kavanagh et al. 2019), which means that in the plane of the sky, this isocontour would have an approximately circular shape (the shape is only perfectly circular for a spherically symmetric wind though).
We can now imagine a situation in which a point-source radio emitter is embedded in the wind, which is optically thick at radio frequencies. The point source could be, for example, an exoplanet. Exoplanets, if they are magnetised, are believed to emit at radio frequencies, analogously to the magnetised planets in the solar system (Farrell et al. 1999). Their emission is cyclotronic and thus takes place at cyclotron frequencies where B p is the planet's magnetic field, c the speed of light, and e and m e are the electron charge and mass. We explore here two different scenarios. In the first scenario, we investigate whether the emission of the point source could propagate though the wind of the host star, while, in the second scenario, we assume the emission can propagate, but it is attenuated. The planetary radio emission can only propagate in the stellar wind plasma if the cyclotron frequency of emission f c is larger than the stellar wind plasma frequency f p everywhere along the propagation path, where Left: Number density of the wind of a solar-like star with a mass-loss rate of 2 Â 10 À12 M yr À1 . A planet is considered to orbit at 0.02 au, with the observer looking towards the system from the negative xdirection. The dashed line shows the radio photosphere where 50% of a 30-MHz wind emission is produced. The hypothetical 30-MHz radio emission of this planet is increasingly more attenuated after the planet ingresses the radio photosphere and the attenuation peaks at orbital phase / ¼ 0:5. Its emission is least attenuated at / ¼ 0. Given that the position of the radio photosphere is linked to the stellar wind properties, monitoring of planetary radio emission could allow one to derive stellar wind properties. Right: The situation on the left panel only occurs if the plasma emission f p is below the cyclotron frequency f c of the planetary emission (white area). If the wind of the host star has a high mass-loss rate and the planet has a weak magnetic field, such that f p [ f c , then the planetary radio emission cannot propagate through the wind of the host star (grey area). Detections of planetary radio emission can thus place an upper limit on the mass-loss rate of the star (Eq. 8). Note that the wind parameters used to produce this figure is based on models of the weak-lined T Tauri star V830 Tau. Images reproduced with permission from [left] Kavanagh and Vidotto (2020), copyright by the authors; and from [right] Vidotto and Donati (2017), copyright by ESO f p ¼ n e e 2 pm e 1=2 ¼ 9 n e 1 cm À3 1=2 kHz: The local electron density n e of a fully ionised hydrogen stellar wind is n e ¼ n=2, where n is the total particle density. The condition f c [ f p is met when The density is related to the mass-loss rate of a stellar wind by the mass continuity equation. Therefore, the condition expressed in Eq. (6) can be translated into where we assumed a mass-loss rate of a steady, spherically symmetric, fully ionised hydrogen wind: _ M ¼ 4pr 2 orb m p n e ðr orb Þuðr orb Þ, with uðr orb Þ being the wind speed at the orbital distance r orb . The right panel in Fig. 7 shows the region of parameter space where radio emission of the exoplanet V830 Tau b could propagate through the wind of the host star . As we can see, the propagation depends on the combination of planetary and stellar wind characteristics. As such, if one can detect radio emission from the planet, the planetary magnetic field can be derived from the frequency of the emission (Eq. 4), and an upper limit of the wind mass-loss rates can be estimated (Eq. 8).
Assuming the planetary radio emission can propagate through the wind of the host star, there is still the question of by how much the planetary emission is attenuated by the wind. To answer this question, we turn our attention back to the left panel in Fig. 7 and more specifically to the location of the planet through its orbit with respect to the radio photosphere (dashed curve). The planetary radio emission is least attenuated (i.e., least absorbed by the wind of the host star) when the planet is between the observer and the star at orbital phase / ¼ 0. As the planet ingresses in the region where the wind is optically thick, the emission from the planet will get increasingly more attenuated, due to an increase in optical depth. In the case shown in the left panel of Fig. 7, the reduction in radio emission from the planet starts at orbital phase / ¼ 0:36, it is maximum when the planet is the furthest from the observer (at / ¼ 0:5) and decreases back again until the planet emerges from the radio photosphere (at / ¼ 0:64). Given that the size of the radio photosphere depends on the physical properties of the stellar wind, if radio emission is detected from the planet for a certain fraction of the orbit, one can then use this information to further constrain the stellar wind properties (Kavanagh and Vidotto 2020).

Constraining stellar winds (upper limits) though X-ray emission
As I will discuss further on, in the Sun, the bulk of the X-ray emission originates from active regions, whose magnetic loops can confine hot plasma. Nevertheless, a hot stellar wind, optically thin in X-rays, can also contribute to a fraction of the stellar X-ray emission. Assuming that the total stellar luminosity is L X ¼ L AR X þ L wind X , with the superscript 'AR' denoting the contribution from active regions, we have that L wind X \L X . Therefore, we can use the observed stellar luminosity L X to infer the upper limit of the wind contribution L wind X . As I will show next, this can then place a constraint (upper limit) in the wind mass-loss rate.
The X-ray emission of the wind can be estimated as where m is the X-ray emissivity integrated over the X-ray frequency range ½m 1 ; m 2 and V is the volume of the emitting wind. Considering the wind X-ray emission is caused by free-free radiation, the emissivity of an optically thin fully ionised hydrogen plasma is given by (Rybicki and Lightman 1986) where g ff is the Gaunt factor of the order of unit (Karzas and Latter 1961), h is the Planck constant and k B is the Boltzmann constant. Note that, at temperatures from below 1 MK up to several MK, lines dominate the X-ray spectrum. Therefore, when assuming that the emissivity is entirely due to free-free radiation, we are providing a lower limit on the emission, given the probable temperature range of coronal winds. Let us first estimate the wind emissivity in the X-ray range where we assumed g ff ' 1. Solving the integral analytically, one gets With this, we get that the wind contribution to the X-ray luminosity is where we assumed the wind is isothermal and the emission measure is where we used the fact that for a fully ionised hydrogen wind, we have n i ¼ n e ¼ n=2.
Equation (13) shows that the X-ray emission coming from the wind depends on the density profile n of the wind. In general, the relationship between density and wind mass-loss rate is not straightforward, as _ M depends also on the velocity structure u of the wind ( _ M / ur 2 n). However, for non-magnetised, thermally-driven winds at a given temperature, it can be shown that the mass-loss rate scales linearly with n, as the velocity structure only depends on the temperature of the wind (see Sect. 5.2.1 and in particular Eq. 42). Therefore, Eq. (13) indicates that isothermal winds with higher densities and thus higher mass-loss rates would have higher L wind X . To compute L wind X , we can use stellar wind models to predict how the density of the wind varies with r and plug this in Eq. (14). Here, however, we proceed with a simplified approximation. For an isothermal wind, we can approximate the density structure below the sonic point by a hydrostatic density structure where n 0 is the wind base density, H H ¼ k B T=ðg H m p =2Þ is the scale height of a fully ionised hydrogen wind and g H ¼ GM H =R 2 H is the stellar gravity at the surface. Thus where x ¼ r=R H . The limits of the integral above should be from the stellar surface x ¼ 1 to the observer x max ! 1. However, given the validity of the hydrostatic density structure, the value of x max should not be larger than the sonic point. Numerically, we can see that the integrand decays quite fast with x, having its maximum of 1 at x ¼ 1. For the sake of simplicity, we will take x max ' 2R H . For isothermal winds of solar-like stars, this is actually not a bad assumption, as the bulk of the emission for winds with temperatures ranging from 1.0 to 2.5 MK occurs within 2R H . For the present-day solar wind, with a temperature of ' 1:5 Â 10 6 K and base density of n 0 ¼ 10 8 cm À3 , solving the integral in Eq. (16) numerically gives ' 0:1 and thus EM ' 10 48 cm À3 . Substituting this in Eq. (13), for an X-ray range from 0.2 to 10 keV, our estimated X-ray emission of the solar wind is quite low, on the order of 10 À10 L . The total X-ray luminosity of the Sun varies during the solar cycle from about 7:0 Â 10 À8 to 1:2 Â 10 À6 L (Peres et al. 2000), which is 700-12,000 times larger than our estimated X-ray wind emission.
As we can see, the present day solar wind indeed contributes to a small fraction of the total solar X-ray luminosity. Although we know the base density of the solar wind, this parameter cannot be easily derived for winds of solar-like stars. We can then ask ourselves, if the Sun were instead a distance star, what would be the base density and, thus mass-loss rates, we would derive from its X-ray luminosity? Given that L wind X / n 2 0 , we would obtain a density value that is larger by a factor of ffiffiffiffiffiffiffi ffi 700 p to ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 12;000 p the n 0 value we adopted before. For the same wind temperature of ' 1:5 Â 10 6 K, this implies that our estimated wind mass-loss rates would be larger by a factor of $ 26 to 110 and the derived maximum mass-loss rate of the Sun-as-astar would have been . 5 Â 10 À13 to 2 Â 10 À12 M yr À1 . The Sun-as-a-star experiment illustrates that stellar X-ray emission can only provide upper limits for wind mass-loss rates.
There has not been many examples in the literature where the technique presented here has been used to measure mass-loss rates in cool dwarf stars.  used soft X-ray observations of the M dwarf star YZ CMi to derive an upper limit of its mass-loss rate. Llama et al. (2013) used the X-ray observations of the K dwarf HD 189733 to constrain the base density n 0 , which is a free parameter in their models.

Detecting charge-exchange induced X-ray emission
The interaction between an ionised stellar wind with a neutral ISM can lead to charge-exchange, in which the ISM neutral atom transfers charge to a solar wind ion during a collision. In this process, a highly-charged ion, particularly oxygen, is excited to a high excitation state, which is then followed by a single or a cascade of radiative decays. These decays lead to emission in the X-ray range. Contrary to the charge-exchange process discussed in the context of Ly-a observations (Sect. 2.1) that only occurs at the site of the interaction between wind and ISM, the chargeexchange induced X-ray emission should take place throughout the stellar wind, as neutral ISM atoms penetrate in the astrosphere. This gives rise to an X-ray 'halo'.
This method was first idealised in Wargelin and Drake (2001), using the solar wind as an example. Considering that low-mass stars have winds that are embedded in a partially neutral ISM, Wargelin and Drake (2001) further suggested that stars with mass-loss rates not much greater than solar could also produce chargeexchange induced X-ray emission. This emission has a distinct profile with a steep rise closer to the star and then a slow decay (see Figure 1 of Wargelin and Drake 2001 for a Sun-as-a-star wind emission). In a follow-up study, the same authors observed Proxima Centauri (Wargelin and Drake 2002), an M dwarf star, but unfortunately the signature of charge-exchange induced X-ray emission was not detected. With their non-detection, they were able to place an upper limit for the mass-loss rate of Proxima Centauri of about 14 times that of the present-day Sun.

Accretion onto white dwarfs as probe of the wind of secondary companion
The last method I would like to present here consists of observing signatures of mass accretion in close binary systems, in which the primary is a white dwarf. The assumption is that material from the wind of the secondary is accreted onto the white dwarf and thus, if one can measure the accretion rate, the wind mass-loss rate of the secondary can be inferred. To the best of my knowledge, this method has not yet been used in systems with solar-like stars, and estimates of mass-loss rates so far exist for some M dwarf stars that are members of eclipsing binary systems (Debes 2006;Parsons et al. 2012).
For the interacting stars where only mass-accretion rates onto the white dwarf companions have been reported (Parsons et al. 2012), these accretion rate values can be used as lower limits for the mass-loss rates of the M dwarf stars. These are considered lower limits of wind mass-loss rates, because not necessarily all the mass lost in the stellar wind of the secondary will be accreted into the primary white dwarf. Debes (2006) finds that M dwarf mass-loss rates are about 15-100 times larger than white dwarf mass-accretion rates.
As I mentioned before, this method has been used to model mass-loss rates in M dwarfs only (see Table 1 in Vidotto and Bourrier 2017 for a compilation of such values). At the moment, it is unclear if the coronae and winds of M dwarf stars in these binary systems are similar to those of isolated/non-interacting M dwarf stars.
From the nine methods presented in Sect. 2, the first four are some of the key methods of wind detection, with the other proposed methods mostly used as case studies and applied on individual stars.
3 An observed evolutionary sequence for mass-loss rates? Figure 8 summarises the stellar wind measurements derived using the methods discussed here. Colour indicates the method used in the derivation: blue for exoplanets, orange for astrospheres, green for prominences, and the Sun is indicated in black (values at minimum and maximum of cycle are provided). Note that for the radio observations, only upper limits have been extracted and these are represented by the arrows pointing down (a few upper limits obtained from other methods are also included). The larger symbols are the stars that are relevant for this reviewthey are main-sequence, solar-like stars, while the smaller symbols are either evolved stars or M dwarfs. The solid line is a fit through the larger circles. Figure 8 shows that, overall, mass-loss rates increase with X-ray flux F X and, given that young stars have higher F X , this implies that young stars have overall higher massloss rates. The recent stellar wind measurements from Jardine and Collier Cameron (2019, greencircles), when contextualised with other measurements, suggest that mass-loss rate continues to increase with F X and that a wind dividing line (cf. Sect. 2.1) is no longer required.
Quantifying the wind of the same star through the use of different methods is an ideal approach to verify and validate measurements. Multiple measurements already exists for some stars in Fig. 8. One example is p 1 UMa-in Fig. 8, the astrospheric measurement of p 1 UMa can be identified at F X ¼ 1:7 Â 10 6 erg cm À2 s À1 and _ M=R 2  (Fichtinger et al. 2017). The spin down model of Johnstone et al. (2015a) predicts a ten times larger than solar mass-loss rate for p 1 UMa, i.e., at _ M=R 2 H $ 11 _ M =R 2 . Although the upper limit provided by the radio observation does not contradict the astrospheric measurement nor the spin down modelling, the astrospheric measurement for this star is in contradiction with spindown tracks that are statistically observed for this type of star. As I discussed in Sect. 2.1, it is possible that p 1 UMa is in a high-ionisation region, which might have affected the result obtained in the astrospheric method (see Wood et al. 2014). Another star in Fig. 8 with multiple wind measurements is the M dwarf Proxima Centauri (F X ¼ 1:4 Â 10 6 erg cm À2 s À1 ): the astrospheric method (Wood et al. 2001), the charge-exchange induced X-ray emission (Wargelin and Drake 2002) and the free-free radio emission  all produced upper limits of . 4 Â 10 À15 , . 2:8 Â 10 À13 and . 7 Â 10 À12 M yr À1 , respectively. Unfortunately, in this case, the three measurements are less stringent, as they are all upper limits. Figure 8 also shows that there is a significant spread in the observations. Vidotto et al. (2016) noted that the more active stars can show a significant variation (cyclic or not) of stellar magnetism in timescales of the order of years. Such variations could increase the spread in stellar wind properties, that is worsened when the X-ray flux and mass-loss rates are not contemporaneously derived. Note for example in Fig. 8 Summary of derived mass-loss rates for low-mass stars combining results from the different methods discussed in Sect. 2. The y-axis is given in solar values, i.e., _ M =R 2 , with _ M ¼ 2 Â 10 À14 M yr À1 . Colour indicates the method used in the derivation: blue for exoplanets, orange for astrospheres, green for prominences, black for the Sun at minimum/maximum of its sunspot cycle. Grey arrows indicate upper limits, which are mostly derived from radio observations. The solid line is a power-law fit through the larger circles. The smaller symbols are either evolved stars or M dwarfs, which were not included in the fit and neither were the stars for which only upper limits exist (arrows). The values used in this plot were compiled from the following works: Drake et al. (1993); ; Gaidos et al. (2000); Wood et al. (2001Wood et al. ( , 2002Wood et al. ( , 2005aWood et al. ( , 2014; Wood and Linsky (2010) The mass-loss rate of the Sun, however, does not seem to vary significantly during the cycle (Cohen 2011), with Finley et al. (2019 estimating that the solar wind mass-loss rate varies in the range $ ½1:6; 3:2 Â 10 À14 M yr À1 as shown by the black line in Fig. 9 (see also Wang 1998). However, other wind properties, like density and velocity and their latitudinal distribution vary more significantly during the cycle (Fig. 10). When the Sun is at minimum activity, its large-scale magnetic field is dominated by a dipolar field, whose axis is roughly aligned with the rotation axis (Sanderson et al. 2003;DeRosa et al. 2012;Vidotto et al. 2018a). In this case, the solar wind velocity distribution is organised into a faster stream ( $ 700-800 km/s) emerging from high latitude coronal holes and a slower stream ( $ 400 km/s) around the equatorial plane (McComas et al. 1998). Conversely, at solar maximum, the large-scale magnetic field geometry is more complex, with a dipolar component vanishing and giving rise to higher-order fields, in particular, even-mode components become more important, such as the quadrupole component (DeRosa et al. 2012;Vidotto et al. 2018a). The solar wind speed in this case has a more complex distribution, as shown in the right panel of Fig. 10. The solar case illustrates the dependence of wind properties on stellar magnetism and later on in this review I will come back to the evolution of stellar magnetism.
I would like now to come back to the fits in Figs. 4 and 8, which represent the fits to the observations of mass-loss rates for solar-like stars. Using measurements derived from the astrosphere method (Fig. 4), Wood et al. (2005a) found that _ M=R 2 Ã / F 1:34AE0:18 X , with the fit being valid for F X .10 6 erg cm À2 s À1 . On the other hand, the solid line in Fig. 8 fits through all the solar-like star measurements (larger symbols) derived from several of the methods discussed in Sect. 2. Note that this fit shows no break within the range of F X shown in the x-axis. For those stars, I find a less steep dependence with F X , with a significant (20%) 1-r uncertainty in the slope of the fit 4 I can now derive an approximate evolutionary sequence for the solar wind from this. I am only focusing on stars older than $ 600 Myr, after which their rotational evolution has converged (see Sect. 4.2). These are also stars that are no longer in the saturated regime, i.e., those for which X-ray increases with rotation (see Sect. 4.3).
There are several relations in the literature that link X-ray flux or luminosity with age or rotation rate (Ribas et al. 2005;Güdel 2007;Reiners et al. 2014). Here, I adopt the relation from Güdel (2007), who found that L X / t À1:5AE0:3 , for solar-like stars in the non-saturated regime. From this, Eq. (17) becomes Wood et al. (2005a) did this exercise, albeit using a different X-ray-age relationship. Using their fit, the line presented in Fig. 4, Wood et al. (2005a) derived that _ M / t À2:33 , which is steeper than our derivation. These two power-laws are shown in Fig. 11. Both of these power-laws indicate that the solar wind had higher mass loss rate in the past. The dispersion is however quite significant, due to all uncertainties in the involved power-laws. For example, the uncertainty in the slope led Wood et al. (2005a) to estimate a solar-wind mass-loss rate in the range ½5; 30 Â 10 À12 M yr À1 for when the Sun was 600 Myr old.
Regardless of which of the two fits we use ( _ M / t À2:33 or _ M / t À0:99 ) for mainsequence stars, there is still a great unknown on how the solar wind has evolved from early stages in the main sequence (or final stages of the pre-main sequence) until the age of $ 600 Myr. Can we extrapolate the t À0:99 dependence to early times  (cycle 22), the solar magnetic field resembles an aligned dipole and the solar wind shows a bimodal velocity distribution, with faster streams emerging from high-latitude coronal holes and slower streams remaining in the equatorial plane. Right: At cycle maximum (cycle 23), the solar magnetic field geometry is more complicated, which is reflected in the solar wind velocity distribution. The background image shows a zoom-in of the solar corona extending out to a few solar radii, while the polar plot shows the solar wind speed measured by Ulysses at several au from the Sun. Colour indicates the magnetic field polarity (red for outward, blue for inward). Images reproduced with permission from McComas et al. (2003), copyright by AGU in the main sequence? It has also been suggested that mass-loss rates increase with rotation (towards lower ages) but then saturate. This suggestion was discussed in Johnstone et al. (2015a) and is necessary to avoid the most rapidly rotating stars to spin down too fast. I will come back to this in Sect. 5.3. In Fig. 11, I represent the mass-loss rate saturation only schematically. Figure 11 also shows the mass-loss rate derived for the weak-lined TTauri star V830 Tau, which is considered to be a 2 Myr-old 'baby Sun'. With a mass of about 1M , this star still has an inflated radius of 2R . The upper limit for its mass-loss rate is estimated to be \3 Â 10 À9 M yr À1 , with likely values ranging between ½1 Â 10 À12 ; 1 Â 10 À10 M yr À1 .
Summarising the discussion presented in this Section, it seems natural that the solar wind has evolved from a high mass-loss rate at early ages until today. However, observations so far do not give us very stringent limits on how precisely this evolution took place. radiation in X-rays and ultraviolet, among others. All these magnetic proxies evolve with stellar rotation, the clock that tracks the passage of time in solar-like stars.
This clock is ultimately regulated by stellar winds, which remove angular momentum from the star and thus force single solar-like stars to spin down with time. In this section, I discuss the observational point of view of the evolution of stellar rotation, activity and magnetism as they play important roles in the theory of winds of solar-like stars, which will be discussed in Sect. 5. For a recent review on the theory of stellar magnetic field generation and links with stellar rotation, I point the reader to Brun and Browning (2017).

Evolution of magnetism
Stellar magnetism can be observationally probed with different techniques. A technique that has been particularly successful in imaging stellar magnetic fields is the Zeeman Doppler Imaging (ZDI, Donati and Brown 1997; for a review see Donati and Landstreet 2009). Through a series of circularly polarised spectra (Stokes V) distributed over one or more rotational cycles, this technique has allowed the reconstruction of the large-scale surface fields (intensity and orientation) of more than one hundred cool stars to date (e.g., Donati et al. 2006;Morin et al. 2008;Petit et al. 2008;Marsden et al. 2011;Mengel et al. 2016;Folsom et al. 2016). The top panel in Fig. 12 illustrates an output of this technique (Fares et al. 2009), in which the three components of the stellar magnetic field are reconstructed: radial, azimuthal (West-East) and meridional (North-South). For comparison, I show in the middle panel a synoptic map from the Sun, produced using SOLIS data (Gosain et al. 2013). The difference between both sets of maps is striking. Firstly, the solar magnetic field strength is an order of magnitude larger than the stellar map. Secondly, we see all this salt-and-pepper structure in solar maps that are not seen in the stellar maps. The reason for these differences lies ultimately in the resolution of the data. Due to its proximity, synoptic magnetic maps of the Sun can be reconstructed down to much smaller scales than stellar maps.
If we decompose these maps using spherical harmonics, a stellar map would typically reach up to a maximum harmonic order 5 ' max $ 10. For the solar map, on the other hand, the high resolution allows harmonics of order ' max ¼ 192 to be achieved. One can estimate how the maximum order is linked to angular resolution as Dh ' 180 =' max , which means that the surface of the Sun can be mapped down to $ 1 , while for a star similar to the one shown in the top panel in Fig. 12, the resolution is on the order of $ 23 . Because of the lower resolution of ZDI maps, magnetic fields of opposite polarities that fall within an element of resolution cancel out. Given that the small-scale field (e.g., concentrated in spots and active regions) have high intensities and that these fields cancel out, the ZDI maps of older solarlike stars do not reach the large strength fields observed in the Sun and allow only the large-scale field to be reliably reconstructed Arzoumanian et al. 2011;Lang et al. 2014).
In the bottom panel of Fig. 12, I show how the Sun would look like if it were observed as a star. This is done by filtering out the small-scale field of the solar synoptic map. In practice, I decompose the solar map shown in the middle panel using spherical harmonics as presented in Vidotto (2016b). This allows me to calculate the spherical harmonics coefficients for each given harmonic order '. Note that the smaller scale structure is described by increasingly larger values of '. Thus, to filter out the small-scale component, I only use the derived coefficients up to a maximum harmonic order of ' max ¼ 5 for reconstructing the large-scale field. This method is commonly used in the literature to separate the large-scale field, i.e., the low harmonic orders (e.g., DeRosa et al. 2012;Petrie 2013;Vidotto et al. 2018a;Lehmann et al. 2018Lehmann et al. , 2019. With this, we can now more easily compare the largescale field of the Sun (bottom panel of Fig. 12) to the ZDI map of a solar-like star (top panel).
A thorough investigation of the limitations of ZDI was recently presented in Lehmann et al. (2019), a study that incorporated simulated data in the ''ZDI machinery''. The authors then compare what the ZDI method derives against the controlled input data. This is done for simulated stars with two different inclination Fig. 12 Top: The reconstructed stellar magnetic field using the ZDI technique for the F-type star s Boo (Fares et al. 2009). Middle: Synoptic map of the Sun plotted with data from Gosain et al. (2013). The solar map, due to its increased resolution, shows a lot more structure in the surface magnetic field. Bottom: To compare the solar map with the stellar map, I filtered out the small-scale magnetic field structure of the Sun (Vidotto 2016b). At the top of each panel, I show the maximum harmonic order ' max for each map and its typical spatial resolution angles of the spin axis with respect to the plane-of-sky (20 and 60 ) and three rotational periods (17, 19 and 27 days). Lehmann et al. (2019) showed that for these low-activity stars, ZDI recovers relatively well the large-scale field morphologies, but magnetic energies can be underestimated by up to one order of magnitude (roughly a factor of ffiffiffiffiffi 10 p $ 3 in field strengths). They also showed that the reconstruction of the field geometry is less reliable if the star is viewed more poleon (lower inclination angles). For this reason, over the years, ZDI studies have favoured mapping the magnetic field of stars that have intermediate inclinations.
One of the major advantages of ZDI maps is that they allow the derivation of the three components of the stellar magnetic field topology. Even though this topology is only restricted to the large-scale field, it can still be useful in stellar wind modelling. As I will show in Sect. 5, surface maps have indeed been increasingly adopted in 3D simulations of stellar winds (e.g., Vidotto et al. 2012;Alvarado-Gómez et al. 2016b;Boro Saikia et al. 2020) and the limited resolution of ZDI magnetograms has been demonstrated not to affect stellar wind models Boro Saikia et al. 2020). The reason for this is the following. If you consider a multipolar field, the dipolar field (' ¼ 1) decays with r À3 ; a quadrupolar field (' ¼ 2) with r À4 ; an octupolar field (' ¼ 3) with r À5 ; and so on. Generalising, a field with a harmonic order ' will decay with r Àð'þ2Þ . This means that small-scale fields (i.e., with large ' values) have very short reach. Thus, it is the large-scale field that is embedded in winds of stars. Stellar winds flow through open magnetic field lines, which ultimately are formed by the stellar large-scale fields that were blown open by the wind outflow. The small-scale field can also affect stellar winds, through, for example, determining the ''micro-physics'' of wind acceleration, which is key for understanding the heating of the solar lower atmosphere (Suzuki and Inutsuka 2005;Cranmer et al. 2007;Shoda et al. 2019).
While the ZDI is blind to the small-scale magnetic field, the Zeeman broadening technique is not. This technique uses Zeeman-induced line broadening of unpolarized light (Stokes I) to derive the average unsigned surface magnetic field hjB I ji (Johns-Krull 2007; Shulyak et al. 2019). When we compare the unsigned average field strength hjB V ji from ZDI maps with hjB I ji, we see that hjB V ji is only a fraction ( $ 10%) of hjB I ji (Reiners and Basri 2009;Morin et al. 2010). Note that here I use the subscripts I and V to differentiate between Stokes I and Stokes V derived fields. For a review on this technique, see Reiners (2012).
hjB I ji is the product of the intensity-weighted surface filling factor of active regions f and the mean unsigned field strength in these regions, which is roughly assumed to be the same as the equipartition field B eq : hjB I ji ¼ fB eq (Cranmer and Saar 2011;See et al. 2019b). The equipartition field is found by balancing the thermal and magnetic pressure at the photosphere of the star. Recent results (Cranmer and Saar 2011;See et al. 2019b;Kochukhov et al. 2020) suggest that B eq itself does not change significantly from star to star, but the filling factor f increases for fast rotators (younger stars). Eventually, the surface of the star becomes covered in active regions and the filling factor saturates for the very fast rotating stars at f $ 1. This leads to a saturation in the magnetic flux of fast rotating stars (Fig. 13a,  Reiners 2012). Figure 13 shows a compilation of a few trends derived empirically considering measurements from these two aforementioned techniques: Zeeman broadening (left panels) and ZDI (right panels). The top panels show how magnetism varies with Rossby number. 6 The Rossby number Ro is defined as the ratio between rotation period P rot and convective turnover time s c Three panels of this figure show magnetic fluxes, defined as where the unsigned average field strength from ZDI observations is given by Here, h and u are the surface colatitude and longitude, and jB ZDI ðh; uÞj is the unsigned surface field strength of the ZDI map at ðh; uÞ. Note that in Fig. 13b and d, only the radial component was used to calculate / V . The trend we see in Fig. 13a is also seen in other tracers of stellar activity, such as in X-ray luminosity versus rotation. I will come back to this later on this Section. This trend consists of two parts: a flat, or saturated, part that is independent of stellar rotation, and a power-law that shows magnetism/activity decreasing with the increase of Ro (i.e., towards slow rotators). Saturation occurs in fast rotating stars and the break occurs roughly at Ro $ 0:1. Here, we are not focusing on the red squares in panel a nor on the open circles in panel b, as those represent mid to late-M dwarfs. Panels a and b in Fig. 13 show that both the large-scale field hjB V ji and the total field hjB I ji behave similarly and are summarised as follows: -I computed a rough fit for RoJ0:1 for panel Fig. 13a and found a power-law slope of À 1:41 AE 0:22 (cf. with À 1:2 from Saar 2001), which is consistent with the slope found in ZDI studies of À 1:19 AE 0:14 (panel b). -The saturation in panel a occurs at about hjB I ji sat $ 2400 G, or, assuming an average stellar radius of 1 R for the solar-like stars, a magnetic flux of / I;sat $ 1:4 Â 10 26 Mx. Saturation for ZDI measurements is currently not clearthe indicative grey dashed line in panel b is at / V;sat $ 1:8 Â 10 24 Mx (see discussion in Vidotto et al. 2014a) Another very interesting relation is between stellar X-ray luminosities L X and magnetic fluxes, as shown in the bottom panels of Fig. 13. Combining solar data and Zeeman broadening measurements, Pevtsov et al. (2003) reported a large dynamical range spanning about 12 orders of magnitude in L X and / I . If we focus only on dwarf stars (crosses in Fig. 13c), the dynamical range is more modest and spans about 2 orders of magnitude, which is also similar to the dynamical range in / V (Fig. 13d). Comparing panels c and d, we see that in both cases, the large-scale field and the total field increase with X-ray luminosity. This is an indication that magnetic fields power coronal activity, in the form of X-rays. Considering only the dwarf stars, Pevtsov et al. (2003) found a power-law slope of 0:98 AE 0:19 for Zeeman broadening measurements, while ZDI measurements show a slope of 1:80 AE 0:20 (Vidotto et al. 2014a). Technically, these slopes are similar within 3r, but the differences in slopes could also be real. If this is the case, this could indicate different 'efficiencies' in producing large-and small-scale fields Morin et al. 2008). One needs to keep in mind that the X-ray values and magnetic field values shown in each panel are not contemporaneous. For example, similar to the Sun, stars have cycles, which affect both field strengths and X-ray luminosities. Thus, non-contemporaneous X-ray and magnetic measurements would lead to increased scatter in the relations shown in Fig. 13. The scatter is seen in all four  et al. (2003), copyright by AAS panels for stellar measurements (it is less evident in panel c, due to the large dynamical range of the plot). Figure 13 shows that overall, hjB I ji and hjB V ji follow activity trends that are similar to those reported in chromospheric and coronal activity proxies (e.g., Skumanich 1972;Wright et al. 2011). An interesting question still remains: how does the magnetic field evolve with age itself? Figure 14 shows that the average magnetic field intensity decays with age À0:65AE0:04 (Vidotto et al. 2014a). This trend is valid over three orders of magnitude in field intensity (from G to kG fields) and four orders of magnitude in ages (from Myr to several Gyr). Although the trend is clear, it also presents a large scatter. Part of this scatter is due to short-term evolution of magnetism, including magnetic cycles. Additionally, in the case of stars younger than $ 600 Myr, part of the scatter is also caused by stars of similar ages and masses, but different rotation rates, showing different levels of magnetisation (see Figure 7a in Folsom et al. 2016, where we see a factor of 10 in hjB V ji for solar-mass stars at 120 Myr). It is fascinating, nevertheless, to see that the relation in Fig. 14 is in accordance to the observed decrease of X H with approximately the square-root of age. I will discuss age-rotation relations in Sect. 4.2.

Evolution of stellar rotation
Stellar clusters are laboratories for stellar evolution models. In particular, observations of rotation rates of stars in clusters at different ages allow us to draw an evolutionary history for stellar rotation. Figure 15 shows age-rotation diagrams from Gallet and Bouvier (2015). The observations of rotation rates of stars with masses $ 1M are shown as pluses, the Sun is depicted as an open circle and rotational velocity dispersion of old field stars in the Galactic disc is represented by the rectangle labeled 'OD' (old disc). The right panel shows, in addition, model tracks as solid (describing the rotation of the envelope) and dotted (same, but for the core) lines. The observations shown in these diagrams are from stars in open clusters, which have good age estimates. The variety of stellar rotation rates found in each cluster results in the vertical alignment of the pluses. For each cluster, the blue, green and red diamonds represent the 25th, 50th, and 90th rotational percentiles, respectively, and the associated models are shown in curves with the same colours.
Two clear behaviours are seen in Fig. 15. Firstly, it is noted that stars, even those with different rotation histories at the early stages of their lives, converge to a unique rotation-age relation approximately at the age of the Hyades open cluster (' 625 Myr, Delorme et al. 2011). From then on, observations show that their rotation rates evolve as where b is known as the magnetic braking index. The equation above forms the basis of the gyrochronology dating method (Barnes 2003), whereby a measurement of rotation rate allows one to infer the age of the star. Skumanich (1972)  Myr, as rotation rates are relatively simple to measure and allows us to derive the age of stars. The problem of deriving ages using this method really starts for stars younger than $ 600 Myr. The second behaviour seen in Fig. 15 is that younger clusters show a wider dispersion of rotation rates, as opposed to older ones. Given that the dispersion of rotational velocities is already seen at very young, Myr-old open clusters, the dispersion is attributed to different 'initial' conditions acquired at star formation. It is believed that a star born as, say, a fast rotator will remain as a fast rotator during its evolution, i.e., it will evolve along the upper envelope shown by the blue solid line in Fig. 15b. Broadly speaking, the rotational evolution models represent three phases. Starting from the young ages, we find the first phase, when the profile shows a constant X H . This is the disc-locking phase, in which the accretion disc regulates stellar rotation. Meanwhile, the star is contracting onto the main sequence. Once the disc is dispersed, disc locking ceases to exist and only contraction continues to take place and, mostly due to angular momentum conservation, the star then spins up. This is the only phase in which rotation increases and we see the peak in rotation happening at around 20-30 Myr. From there on, rotation decreases with age, which means that single stars on the main sequence will spin down with time. The main driver for this spin down, as we already discussed in the introduction, is the magnetised stellar wind, which carries away angular momentum. The particular model I present in Fig. 15 (Gallet and Bouvier 2015) includes more physical aspects, such as the evolution of the moment of inertia (mass, radius) and coupling between the rotations of the core and of the envelope. Other models exist in the literature-they differ, for example, on which stellar evolution code is used, different treatments for the core-envelope decoupling, or different prescriptions for the 'stellar wind braking law', among others (e.g., Kawaler 1988;MacGregor and Brenner 1991;Keppens et al. 1995;Spada et al. 2011;Reiners and Mohanty 2012;Epstein and Pinsonneault 2014;Johnstone et al. 2015a;Matt et al. 2015;Amard et al. 2019). The braking law describes how angular momentum loss rate _ J depends on the stellar properties, such as mass, radius, stellar magnetism, rotation, surface metallicity, and stellar wind mass-loss rate. The constructions of braking laws are based on theoretical models, empirical scalings, or a combination of both. I refer the reader to Amard et al. (2016) for a study on how different braking laws can affect rotational evolution models and thus the goodness of fits to observed rotation rates. Reviews on the topic include the extensive lecture notes from Palacios (2013) and Bouvier (2013). I will come back to angular momentum losses in Sect. 5.3.

Evolution of coronal activity
The third wind ''ingredient'' I would like to discuss is stellar activity. In particular, I would like to focus on coronal activity. The corona starts above the chromosphere (and above the transition region) and it is characterised by a high temperature (' 10 6 K) and low density plasma. The emission from the corona falls in the X-ray and UV part of the electromagnetic spectrum. The solar wind is, in a sense, a continuation of the solar corona, expanding into the interplanetary space. However, in the cool-star community, it is common to separate the corona as the region of complex, closed magnetic field lines and the wind as belonging to the region of open field lines. For this reason, the corona is sometimes referred as the 'closed corona' Fig. 15 Rotational evolution from solar-mass stars. Pluses are observations of rotation rates of open cluster stars, circle is the present-day Sun and the rectangle shows the velocity dispersion of old disc field stars. The grey shaded bars in a highlight the velocity dispersion at each age. The red, green, and blue diamonds represent the 25th, 50th, and 90th rotational percentiles, respectively. In b, models are overplotted to these observations. Image adapted and reproduced with permission from Gallet and Bouvier (2015), copyright by ESO and it is the source of X-ray and UV emission, with the wind-bearing regions being X-ray dark (see Sect. 2.7). There is still no consensus at how precisely the corona and the wind are heated to such high temperatures and it is possible that a combination of different mechanisms heat the plasma in different regions and at different times, as discussed in the recent review on solar coronal and wind heating by Cranmer and Winebarger (2019).
Coronal activity is observed in other stars through X-ray and UV observations. Figure 16 shows how the normalised X-ray luminosity R X ¼ L X =L bol evolves with Rossby number. Similar to what was seen in Fig. 13a and b, the emission saturates at R X ' 10 À3 for fast rotating stars with Ro.0:1 and a decrease in X-ray emission is seen with Ro À2 , or $ X 2 H , in the unsaturated part (Pizzolato et al. 2003;Reiners et al. 2014). As discussed in Sect. 4.2, rotation can be used as a proxy for age and thus this decay in X-ray can be attributed to age evolution. Three factors contribute to the large spread seen in Fig. 16: inaccuracies in X-ray measurements, variability caused by stellar cycles and intrinsic differences between stars (some being more active than others for same mass, age and rotation, Johnstone et al. 2021). Figure 17 shows the evolution of high-energy radiation of solar-type stars. The left panel divides the high-energy flux into 6 bands: the [1, 20] Å band corresponds to hard X-rays, [20, 100] Å to the ROSAT band, [100, 360] Å to the EUVE band and [920, 1180] Å corresponds to FUSE (Far Ultraviolet Spectroscopic Explorer) bands. There is a gap between 360 and 920 Å for solar-like stars, so the fit for this Fig. 16 The X-ray-rotation activity relation shows how the normalised X-ray luminosity R X ¼ L X =L bol varies as a function of Rossby number (fast rotators on the left, slow rotators on the right of the x axis). Image adapted from Reiners et al. (2014), copyright by AAS band shown in the figure relies on interpolations from neighbour wavelength bands. The slopes in the power-law fits are shown in Table 5 of Ribas et al. (2005)-they range from the steepest value of À 1:92 for the hard X-rays to À 0:85 for the less energetic band. The fluxes are normalised to the present-day Sun values. We see that the contribution of the X-ray part of the spectrum is quite strong at young ages, compared to present-day solar values.
The EUV part of the spectrum (100-920 Å ) is relevant in studies of atmospheric escape in planets and exoplanets. Sanz-Forcada et al. (2011) estimated that the EUV luminosity L EUV is related to the X-ray luminosity (5-100 Å ) as with luminosities given in erg s À1 . Using the X-ray-rotation activity relation (similar to shown in Fig. 16), the equation above and rotational evolution models (similar to shown in Fig. 15), Tu et al. (2015) derived an evolutionary model to describe the evolution of L X and L EUV that is shown in the right panel of Fig. 17. The spread in rotation rates observed for young stars implies a spread in XUV (Xray plus UV) luminosities, such that slow rotators at a given age will have lower XUV luminosities than a fast-rotating star at the same age (Tu et al. 2015).  Fig. 17 Evolution of high-energy irradiation of solar type stars. Left: The observed evolution of highenergy flux at a given orbital distance at different bands. Right: The X-ray (left axis) and extreme ultraviolet (right axis) luminosities derived from rotational evolution tracks (e.g., Fig. 15) as a function of stellar age. Images reproduced with permission from (left) Ribas et al. (2005), copyright by AAS; and (right) Tu et al. (2015), copyright by ESO than the Sun. To compare to the 'total' EUV band (100-920 Å ) used in Eq. (23), one should then sum F ½100À360 EUV and F ½360À920 EUV , obtaining F EUV ¼ 10 2:04 F 0:681 X þ 0:31F 0:626 Comparing Eqs. (23) and (26), we notice a smaller exponent is the latter. This is because one relation uses flux and the other luminosity. Johnstone et al. (2021) demonstrate that using surface fluxes reduce spread in empirical relations between X-ray, Ly-a and overall EUV emissions.
In the Sun, X/EUV flux is very well correlated with sunspot number and the solar surface magnetic flux (Figs. 1 and 8 in Hazra et al. 2020). Figure 18 shows how the solar XUV flux varies with surface magnetic flux. The XUV flux is calculated over the range 5-912 Å , using data from the Solar Extreme Ultraviolet Experiment (SEE) instrument on the NASA Thermosphere Ionosphere Mesosphere Energetics Dynamics (TIMED) spacecraft (Woods et al. 2005). The solar synoptic maps used to compute the magnetic fluxes are from the Helioseismic Magnetic Imager (HMI) on board the Solar Dynamics Observatory (Scherrer et al. 2012). The empty symbols consider the HMI data in its full resolution, while the filled symbols only considers magnetic field components with spherical harmonics orders below ' max ¼ 10. This represents the Sun-as-a-star case. Figure 18 shows the relatively tight correlation between surface XUV flux and magnetic fluxes considering both the full resolution maps as well as the low resolution maps (large-scale fields). This tight correlation led Hazra et al. (2020) to  Hazra et al. (2020) suggest that the large-scale fields observed in other stars (e.g., Fig. 12) could be used to infer the stellar X/EUV fluxes of stars through the relation with / V defined in Eq. (20). Note, however, that if one assumes a lower harmonics order (e.g., ' max $ 5) for the 'Sun as a star' case, the scatter in the relation increases substantially. The reason why XUV emission is better related with higher harmonics order is that the solar XUV emission mostly comes from compact active regions, which are characterised by smaller-scale fields, i.e., with higher harmonics orders (see discussion at the beginning of Sect. 4.3).
5 Where are we in terms of modelling?
The three observed ingredients discussed in Sect. 4 are often used to inform stellar wind models, which developed from models of the present-day solar wind. Here, I do not go into details of the present-day solar wind models. I refer the reader to some recent reviews in the area, such as Gombosi et al. (2018)  This section is focused on models of winds of solar-type stars and, in particular, on creating an evolutionary sequence that considers and reproduces the stellar observations available to us, namely rotation, magnetism and activity (Sect. 4). In Sect. 5.1, I provide a brief overview of the main types of models used to describe winds of solar-type stars. One important physical process that is not frequently dealt with in solar wind models is angular momentum loss. This is mostly because solar wind studies have focused on the properties of the Sun today: its open magnetic flux, mass flux (densities and velocities) impinging on the solar system objects, and its mass-loss rate. The solar wind is usually treated as a snapshot in evolutionary timescales, and rotation, although considered in solar wind models, has its feedback neglected on the Sun's long-term evolution. For stellar astronomers, rotation is the key for the evolution of the Sun and its wind properties. Thus, in Sect. 5.3, I focus on modelling angular momentum losses, after presenting a brief overview of magnetohydrodynamics (MHD) wind models in Sect. 5.2.

Overview on the different treatments used in stellar wind models
Studies of the solar wind have provided guidance for models of winds of low-mass stars. The detailed physical mechanism(s) that provides the heating and acceleration of solar-like winds is not known, similarly to the solar wind itself (e.g., Cranmer 2009). These winds are believed to be magnetically-driven, in which coupling between stellar magnetism and convection transports free magnetic energy, which in turn is converted into thermal energy in the upper atmosphere of stars (Matsumoto and Suzuki 2014), giving rise to a hot corona and a hot wind (J 10 6 K).
The X-ray emitting stellar corona, set by energy input, varies with the properties of the star (Jardine 2004;Güdel 2004), as do the stellar wind properties .
To the best of my knowledge, the currently-available stellar wind models are based on fluid descriptions (magnetohydrodynamics, MHD from now on). Although used to study the solar wind itself, I am not aware of models that have adopted kinetic theory to study the evolution (or even a snapshot) of winds of solar-type stars. In terms of the fluid description, two modelling approaches are used in the study of the hot coronal winds of low-mass stars. The first one explicitly considers atmospheric heating that is deposited at the photospheric (or sometimes the chromospheric) level, while for the second one, a million Kelvin-temperature wind is assumed (i.e., heating is assumed a priori and thus is ''implicit''). I discuss these different treatments next.
(i) Explicitly including heating from the photosphere: A possible scenario to convert magnetic into thermal energy in the atmosphere of stars involves the dissipation of waves and turbulence, based on scenarios developed for the solar wind (e.g., Holzer et al. 1983;Cranmer 2008;Cranmer and Saar 2011;Suzuki et al. 2013;Matsumoto and Suzuki 2014). In addition to depositing energy, the gradient of wave pressure provides a volumetric force that accelerates the wind, similar to thermal pressure gradients. The modelling approach that explicitly considers heating from the photosphere involves a more rigorous computation of the wave energy and momentum transfer, i.e., the computations are done from ''first principles'' (e.g., Hollweg 1973;Holzer et al. 1983;Vidotto and Jatenco-Pereira 2006;Cranmer 2008;Suzuki et al. 2013;Shoda et al. 2019). In these models, the increase in temperature from the colder photosphere to the hotter corona arises naturally in the solution of the equations as does the wind acceleration. Most of the models that treat the stellar wind acceleration starting from the photosphere have focused on the wind dynamics along a single open magnetic flux tube, as, depending on the level of details of the physics involved in the wind acceleration/heating mechanisms, models can become computationally intensive. In particular, a challenging numerical aspect is the large contrast between the dense and cold photosphere and the hot and rarefied corona (e.g., Matsumoto and Suzuki 2012).
Models validated for the solar wind have been applied to stellar wind studies. Suzuki et al. (2013) presented a parametric study of Alfvén-wave driven winds applied to solar-like stars. 7 The authors simulated a wide range of input wave fluxes, starting from the photosphere, by adopting different values of magnetic field strengths and turbulent velocities, representative of stellar values. The left panel of Fig. 19 shows a summary of their simulation results, where the output kinetic luminosity of the stellar wind ( _ Mu 2 r =2, y-axis) is plotted against the the input wave luminosity (x-axis). Different values of average magnetic field strengths are shown by the different symbols. Here, B r;0 represents the magnetic field intensity at the bottom of a magnetic flux tube and f 0 represents the fraction of the star that is covered in open flux tubes. Note the vertical scatter, indicating that similar input wave energies give rise to different wind kinetic energies. The right panel of Fig. 19 shows a similar plot, but instead of the input wave luminosity, the x-axis now shows the wave luminosity as measured at the transition region (top of the chromosphere). Note that these values are smaller than the input wave luminosity at the photosphere, because a fraction of the waves is lost as they reflect back downward, being mostly radiated away. Suzuki et al. (2013) found that a fraction between 1 and 30% of the original input wave energy is able to get into the chromosphere. Of this survival wave energy, part of the energy is lost in the form of radiation and part is used to do work against the gravitational force-thus, only a small amount of the input wave energy at the photosphere is actually used to drive the stellar winds. This study shows that waves cannot easily penetrate in the corona/wind. Overall, for a given magnetic field, an increase in the wave input generates a larger output kinetic energy, until saturation occurs. The consequence of this is that, for a given value of B r;0 f 0 , the mass-loss rate of the stellar wind increases with input wave energy up to a saturation value _ M sat . Beyond this point, the increase in input wave energy seems to no longer increase mass-loss rate. The saturation is also seen to occur at higher values for larger average fields: _ M sat ¼ 7:86 Â 10 À12 ðB r;0 f 0 =½1GÞ 1:62 M yr À1 . This result indicates that stars more magnetically active would have higher 'maximum' mass-loss rates. It is interesting to note that this model links not only the magnetic field of an open flux tube, but the product between the magnetic field of a flux tube and the amount of open flux tubes in the surface.
In Sect. 4.1, I discussed surface filling factor of active regions, f, mentioning that recent results suggest that the magnetic field in active regions does not change significantly from star to star, but the filling factor f increases for fast rotators/ younger stars (Cranmer and Saar 2011;See et al. 2019b;Kochukhov et al. 2020). The filling factor in active regions is likely linked to the filling factor of open flux tubes-if one assumes that the field lines in active regions will remain closed, at the very least we could imagine a scenario where f ' 1 À f 0 , i.e., the star has only f 0 Unfortunately, to determine the amount of open magnetic flux, one needs models, such as potential field extrapolations (e.g., Jardine et al. 2002) or stellar wind models (e.g. Vidotto et al. 2014b;Réville et al. 2015). Nevertheless, even without models, there is an observationally-derived proxy that could tell us a bit more about how open magnetic flux, i.e., the one embedded in the wind, varies with stellar activity. A large fraction ( $ 90%, Reiners and Basri 2009;Morin et al. 2010) of the unsigned field derived in Zeeman broadening measurements can be missed in ZDI measurements due to flux of different polarities cancelling out within an element of resolution (Lang et al. 2014). The bulk of the ''missed'' field is thus in small-scale fields (in regions that are likely smaller than active-region sizes), which are mostly connected to each other forming closed-field structures. The smaller the loop size, and thus the larger its harmonic order ', the faster is its decay with radial distance, given that the magnetic field of order ' decays with r Àð'þ2Þ . These small loops, thus, decay very fast with distance and they do not thread the wind, which represents a more global-scale quantity. The magnetic field that will thread the wind is, on the contrary, a large-scale field, with See et al. (2018) showing that the open flux is largely dominated by the dipolar field alone. ZDI measures large-scale fields, a fraction of which will open up and contribute to the open magnetic flux. Of course, one cannot tell which fraction of the ZDI field will open up and contribute to the wind without a model, but in light of the two available observational quantities, namely the total field (small?large scales) derived in Zeeman broadening measurements and the large-scale field derived in ZDI measurements, the latter represents a better proxy for the wind-bearing magnetic field lines. Thorough discussions on the fraction of open field in the present-day Sun and on other cool dwarfs can be found in Cranmer (2017) and See et al. (2019b). Given that the magnetic flux in ZDI fields increases with X-ray luminosities (see Fig. 13d), this probably means that the magnetic flux in open field lines increase for more active stars. Thus, even if the closed-field region increases in coverage for fast rotators, the remaining area of open field lines, even occupying a smaller portion of the star, harbours a larger amount of magnetic flux, which ultimately pertains to the stellar wind.
Theory and numerical simulations predict that wind quantities such as mass-and angular momentum-loss rates depend on the amount of open magnetic flux (e.g., Washimi and Shibata 1993;Vidotto et al. 2014b;Réville et al. 2015;See et al. 2019a), thus estimating the fraction f 0 of open field lines is important for stellar wind models. This is also seen in the work of Suzuki et al. (2013)- Figure 19 (left), for example, shows that the output kinetic luminosity of the wind depends on the fraction f 0 of the photosphere covered in open flux tubes, along which Alfvén waves with luminosity L A;0 can propagate. (The wave luminosity itself depends on the $ kG-level magnetic field intensity at the bottom of a magnetic flux tube: L A;0 / v A;0 / B r;0 .) (ii) Implicitly assuming heating in wind models: The models that explicitly solve for the energy input in the upper atmosphere of Sun-like stars, such as from Suzuki et al. (2013), reproduce the rapid increase in temperature from the photosphere to the corona. One of their limitations is that, because they are computationally expensive, they are computed along a flux tube, and usually do not consider rotation (however, see Shoda et al. (2020) for first efforts in including rotation in these types of model). These two limitations can be overcome in the second approach for modelling winds of solar-like stars that I discuss now (albeit other limitations appear). This approach adopts a simplified energy equation, usually assuming that the thermal variables are connected by a polytropic relation. Typically, a power-law between the thermal pressure P and the density q or temperature T are used: P / q C , or T / q CÀ1 where C is known as the polytropic index. A value of C ¼ 1, for example, indicates that the temperature is independent of the density, i.e., mimicking an isothermal relation. In the solar wind, the measured effective polytropic index in the corona is around 1.1 (Van Doorsselaere et al. 2011), increasing to 1.46 near Earth (Totten et al. 1995) and decreasing with distance in the outer heliosphere (Elliott et al. 2019). Usually, polytropic wind models assume constant values of C ranging from 1 to 1.15, although a few models treating a distance-dependent C exist for solar wind models (Roussev et al. 2003;Cohen et al. 2007) and stellar wind models (Vidotto 2009;Johnstone et al. 2015b, a).
In polytropic wind models (which can be magnetised or not), the computation often starts at the point where the temperature has already reached coronal values $ 10 6 K (e.g., Pneuman and Kopp 1971;Washimi and Shibata 1993;Keppens and Goedbloed 2000;Matt et al. 2012;Vidotto et al. 2009b;Réville et al. 2015). This approach ignores the physical reason of what led temperatures to increase from photospheric to coronal values. Technically, it is as if the wind base occurred at the corona-if we take the height of the solar corona to be at about 2000 km (e.g., Gary 2001;Yang et al. 2009), this means that the base of the corona, and of these wind models, occurs at $ 1:003 R . In practice, models usually start at 1 R and this small difference in neglected in simulations.
By assuming a polytropic wind, the energy equation is simplified and no additional equation (e.g., describing wave propagation) is computed. This allows us to not only perform three-dimensional numerical simulations, but also to extend the simulation box to much larger distances from the star. For this reason, I call these ''global'' wind models, or large-scale wind models. These models usually consider rotation, allowing us to calculate angular momentum losses (cf. Sect. 5.3). The obvious drawback of such models is that they do not consider the physics of what happens within 0:003 R above the photosphere, which is where the bulk of the heating takes place. Therefore, the temperature, and also the density, at the wind base are free parameters of polytropic wind models. These free parameters are usually determined from complementary observations, such as X-ray observations, radio observations and observationally-derived mass-loss rates (e.g., Holzwarth and Jardine 2007;Johnstone et al. 2015b, a;Réville et al. 2016;Ó Fionnagáin and Vidotto 2018;Vidotto et al. 2018b;Ó Fionnagáin et al. 2021).
To reduce the uncertainty in these free parameters, some wind models have used X-ray properties to constrain, for example, the density and temperature at the wind base (e.g., Holzwarth and Jardine 2007;Réville et al. 2016;Ó Fionnagáin et al. 2019). In the Sun, coronal holes are X-ray dark, and the closed field line region is X-ray bright, with the more active/denser region showing higher temperatures (Peres et al. 2004). Given that the wind flows through coronal holes, relating wind base properties with X-ray properties is at the end an assumption that the thermal properties of coronal holes are proportional to the thermal properties of the closed corona. Such an assumption implies that, given that the corona of a fast rotating star is hotter (Güdel et al. 1997;Telleschi et al. 2005), their wind should be proportionally hot (corona and wind temperatures are not necessarily the same though). Likewise, stars with denser corona would have denser winds in these models. These wind models thus predict that mass-loss rates are higher for young and fast rotating stars.
An alternative approach that has been used to constrain the wind base temperature in other stars is by assuming that the thermal speed (and thus the temperature) at the base of the wind is a fixed fraction (about 22%) of the surface escape velocity (Matt and Pudritz 2008;Matt et al. 2012). Because the activity level is not considered in this approach, two stars with same mass and radius, but with different magnetic activity, would have winds with the same temperature. One uncertainty is that, because the magnetic field flux at the base of open-field regions increases with magnetic activity, the heating efficiency may not be the same during the spin evolution. Given that the surface escape velocity does not change considerably during the main sequence, these wind temperatures change by a smaller amount during stellar spin evolution than the approach using X-ray observations.
There are two advantages of global wind models that I would like to highlight here. Firstly, because the numerical grid can extend out to large distances, it is possible to characterise the stellar wind conditions around exoplanets (Vidotto et al. 2009a(Vidotto et al. , 2010b(Vidotto et al. , 2012Cohen et al. 2011a, b;Llama et al. 2013;Nicholson et al. 2016;Vidotto and Donati 2017;Strugarek et al. 2019). The characterisation of the local conditions of the stellar wind is important to quantify the wind (magnetic and particles) effects on exoplanets, as I will discuss later on in this review. Secondly, the global wind models can also incorporate more complex magnetic field topologies, including synoptic maps derived from ZDI studies. The magnetic map is imposed as boundary condition at the stellar wind base. Usually, the process involves extrapolating the surface field into the computational domain initially assuming the field is in its lowest energy state (e.g., a potential field). After the interaction with the stellar wind flow, the magnetic field becomes stressed. The selfconsistent interaction between magnetic field lines and stellar wind are allowed to evolve, until a relaxed solution is found (for more details, see e.g., Vidotto et al. 2014b). Figure 20 illustrates the initial (left panel) and final (right panel) conditions of a 3D magnetohydrodynamics (MHD) stellar wind simulation, that adopts a polytropic heating and incorporates observed surface magnetic fields at its inner boundary.
Being able to include the diversity of observed magnetic field geometries is important. In the case of the Sun, observations of the solar wind with the Ulysses spacecraft revealed that the geometry of the solar magnetic field affects the velocity distribution of the solar wind (McComas et al. 2008). When the Sun is in its activity minimum and has a simpler magnetic field topology, close to a dipole, the wind Fig. 20 a Magnetic field line extrapolations, assuming a potential field source surface model. The surface magnetic field is derived from ZDI observations . The potential field model assumes the stellar magnetic field is in its minimum energy state. b The field lines are stressed after interaction with stellar wind flow. Based on the simulations presented in Vidotto et al. (2014b) Fig. 21 Left: Velocity of the simulated stellar wind (in km/s) of the planet-hosting star HD 189733, at the position of the orbit of the exoplanet HD 189733b. The blue circles denote an inward component of the magnetic field and the green diamonds denote an outward component. The computed free-free X-ray emission of the coronal wind is shown in the background. The structure that the stellar magnetic field imposes on the X-ray corona is correlated with the structure of the stellar wind. Right: Histogram of wind velocities for 3D MHD simulations of the wind of the young solar-like star HII 296. Three velocity components related to magnetic field geometry can be identified: a slow one at 250 km/s emanating from helmet streamers, a fast component at 500 km/s (dashed-magenta line) from expanded flux tubes and an intermediate velocity of 400 km/s. Images reproduced with permission from [left] Llama et al. (2013), copyright by the authors; and from [right] Réville et al. (2016), copyright by AAS structure is bi-modal, with larger wind velocities in the (high-latitude) coronal holes than in the low-latitude region. On the other hand, when the Sun is at its activity maximum, its magnetic field geometry becomes more complex, which is then reflected in the wind structure. Numerical simulations indeed show that the geometry of the stellar magnetic field affects the wind velocity (e.g., Vidotto et al. 2009b;Llama et al. 2013;Réville et al. 2016;Ó Fionnagáin et al. 2019). Figure 21 shows results of two 3D MHD simulations of winds of solar-like stars that include the observed surface magnetic fields. The left panel of Fig. 21 shows the solution of the stellar wind model of the planet-hosting star HD 189733 from Llama et al. (2013) that was computed using the observationally-derived ZDI magnetic map from Fares et al. (2010). This panel has a similar format as Fig. 10. The background shows the computed X-ray emission of the hot, quiescent corona of the star due to thermal free-free radiation. Note that coronal X-ray emission should come mainly from smaller scale active regions. As the small-scale magnetic structure is not resolved in ZDI observations, the X-ray emission computed in Llama et al. (2013) captures only the emission of the background coronal wind and, as such, provide a lower limit for the emission. Overlaid to the computed X-ray image is the velocity of the stellar wind in km/s at the position of the orbit of the exoplanet HD 189733b. Regions of fast wind correspond to X-ray dark regions in the corona. Slow wind regions tend to lie over the largest helmet streamers in the coronal field. This is a result of the nature of the magnetic force, which generates meridional flows that bring wind from open-line regions (coronal holes) to the top of the streamers (Vidotto et al. 2009b). Therefore the structure that the stellar magnetic field imposes on the X-ray corona is correlated with the structure of the stellar wind. The right panel of Fig. 21 shows a histogram of wind velocities derived in another 3D MHD simulation of the wind of the young solar-like star HII 296, where three velocity components can be identified (Réville et al. 2016). A slow component around 250 km/s originates from wind emanating from helmet-streamer structures of the magnetic field, a fast component around 500 km/s originates due to magnetocentrifugal effects and flux tube expansion (dashed magenta line), and an intermediate component at 400 km/s, which coincides with the non-magnetised polytropic wind solution (dashed black line). Réville et al. (2016) showed that the result of the interaction of the flow with strong, non-axisymmetric field generates a complex wind-speed distribution, which is wider in the case of stars with higher magnetisation and faster rotation (see also (iii) Hybrid approach: Recently, there have been efforts in developing a hybrid model that combines the two approaches described above to study winds of lowmass stars. In these hybrid models, a phenomenological approach of the (solarbased) wave heating mechanism is implemented in three-dimensional simulations of solar/stellar winds, starting from the upper chromosphere (van der Holst et al. 2010(van der Holst et al. , 2014Sokolov et al. 2013;Garraffo et al. 2016;Alvarado-Gómez et al. 2016a, b;Cohen 2017;Boro Saikia et al. 2020;Ó Fionnagáin et al. 2021). These models present a step forward in the 3D modelling of winds of low-mass stars, as they, for example, do not need to impose a polytropic index to mimic energy deposition in the wind. However, the full computation of the heating process from the photosphere is still not available.
In the hybrid models, there are also free parameters that need to be considered, such as the energy flux of waves at the inner boundary and its dissipation length scale (Sokolov et al. 2013;Ó Fionnagáin et al. 2021). In the case of the solar wind, these free parameters can be constrained from solar observations (Sokolov et al. 2013;van der Holst et al. 2014). In the case of winds of other stars, it is less clear how to constrain these free parameters. Of particular relevance is the wave energy flux, or Poynting flux, at the inner boundary of these models (chromosphere). Alvarado-Gómez et al. (2016b) and Cohen (2017), for example, followed the approach of Sokolov et al. (2013) who used empirical relation between X-ray and magnetic fluxes to constrain the wave input Poynting flux. While in a polytropic formalism the heating is not directly related to the surface magnetic field, in their hybrid approach, the input wave flux at the chromosphere is taken to be proportional to the local surface magnetic field strength (Sokolov et al. 2013). By assuming that the Alfvén wave turbulent dissipation is the main source of heating of the solar wind/corona, Sokolov et al. 's approach implies that the total heating power scales with the magnetic flux (Eq. 1 in Sokolov et al. 2013), with a fraction of this total power generating the observed X-ray luminosity seen in scalings from Pevtsov et al. (2003). Sokolov et al. (2013) derived the constant of proportionality between the Poynting flux and the local surface magnetic field as F A =B ' 44:7ð4pR 2 H Þ 0:1488 ' 1:1 Â 10 5 ðR H =R Þ 0:2976 erg cm À2 s À1 G À1 . Note the stellar size dependence, which should be considered when this implementation is applied to other stars (Cohen 2017). This is not the only approach to determine the input wave energy flux. An alternate one was discussed in Boro Saikia et al. (2020), who constrained the wave flux using far-ultraviolet spectral lines that are formed in the upper chromosphere or transition region of stars. For completeness, I also highlight alternative approaches for determining wave energy fluxes used in 1D models (which start from the photosphere, contrary to the hybrid studies discussed in this section, that start in the chromosphere). Saar (2011) andCranmer (2017) constrained the Poynting flux used in their models from the magneto-convection models of Musielak and Ulmschneider (2002), who studied the dependence on the emerging wave energy flux with stellar gravity under the mixing-length theory. The gravity dependence can be particularly relevant for stars that evolve off the main sequence (Ó Fionnagáin et al. 2021). Dependence with mass and metallicity of the emerging wave energy flux has also been considered in the work of Suzuki (2018), when investigating winds of population III stars. In this case, the amplitude of the wave velocity fluctuations dv (note that F A / qhdv 2 i) at the photosphere is determined from the surface convective flux, which is proportional the stellar luminosity: qdv 3 / L H (see their Eq. 15).
By considering a Poynting flux dependent on the local magnetic field strength, the Poynting flux varies across the stellar surface. As a result, hybrid 3D models show stronger contrast between the wind properties and the large-scale magnetic field geometry than those based on a purely polytropic formalism. This technique is recent in comparison to those discussed in (i) and (ii) and I expect it will be further developed in the years to come.

Brief overview of MHD stellar wind theory
Each of the different treatments adopted in wind models have their advantages and their drawbacks. When it comes to modelling rotational evolution, models that explicitly treat heat deposition in the wind are less suitable, as most of these models neglect rotation. To study angular momentum losses, it is more common to adopt polytropic models (i.e., with implicit energy deposition). The first polytropic model for the solar wind was developed by Parker (1958). This model is, in fact, isothermal, which is a special case of a polytropic wind with C ¼ 1. Here, I do not present the derivation of the MHD equations that are often used in stellar wind theory-those can be found in textbooks (e.g., Mestel 1999; Goedbloed and Poedts 2004). Here, I only present some key MHD equations that I will use to discuss stellar wind models and angular momentum losses.

Thermally-driven winds
The first important equation in wind models is that of mass continuity where q is the mass density, and ũ is the wind velocity. In the absence of sinks or sources, mass is conserved in a stellar wind. For the momentum equation, it is useful to use the Lagrangian formalism. The Lagrangian time derivative D/Dt, sometimes also known as the total derivative, is given by The D/Dt operator can be applied to scalar Q or vector Q quantities. The first term on the right hand side is the rate of change of Q computed at a fixed location (i.e., the Eulerian time derivative) and the second term arises because the fluid element with velocity ũ has moved to a new location where the quantity Q has a different value. This term is also known as the 'convective derivative'. In the Lagrangian description, Newton's second law can be simply written as The term of the right hand side represents the sum of all volumetric forces acting on the wind. Once derived, however, the momentum equation is more convenient to use from the Eulerian viewpoint. For a polytropic (or an isothermal) stellar wind with only gravitational forces and pressure gradients acting on the wind, the momentum equation, using the two previous equations, is wherer is the unit vector along the radial direction. This equation can be further simplified by assuming steady state (o=ot ¼ 0) and spherical symmetry (o=oh ¼ 0, Under these conditions, the continuity Eq. (28) can be written as Taking into account that the stellar wind emerges from the 4p steradians of stellar surface, we then define the mass-loss rate as which is constant with radial distance in a 1D stellar wind.
Equations (32) and (34) form the basis of the thermally-driven wind model derived by Parker (1958). The thermal pressure can be written using the ideal gas law where k B is the Boltzmann constant, m p the proton mass, T the wind temperature and l is the mean ''molecular'' weight. In the case of hot stellar winds, no molecules exist, and the mean molecular weight is simply the average mass of the particles. For a purely hydrogen plasma, fully ionised (50% protons, 50% electrons), the mean mass of the particles are lm p ¼ 0:5m p þ 0:5m e ' 0:5m p , where m e is the electron mass, so l ¼ 0:5. In the equation above, I wrote the isothermal sound speed as a ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k B T=ðlm p Þ p . Thus, differentiating the ideal gas law along the radial coordinate, we get In Parker's model, the wind is considered isothermal, thus the terms with dT/dr are null in the previous equation. I am going to assume here that the wind is polytropic, which implies that where I used the ideal gas law to convert between the two equations. From this, the temperature gradient is and Eq. (36) is written as where I used the mass conservation Eq. (33) to write dq=dr and wrote the sound speed for a polytropic gas as Substituting the previous equation in the momentum Eq. (32) and rearranging, we arrive at the momentum equation of a thermally-driven wind: The equation above reduces to Parker's solar wind equation if one assumes the wind is isothermal and thus C ¼ 1 and hence c s ¼ a.
There are two important points to highlight here. Firstly, note that the density cancels out. Thus, in a thermally-driven wind, the acceleration of the wind does not depend on the density of the wind, which leaves the density as a scaling parameter that is usually adopted to match the required mass-loss rate (Eq. 34). Secondly, this equation has a singularity at u r ¼ c s , i.e., when the stellar wind radial velocity reaches sound speed, the denominator goes to zero. The momentum equation for a polytropic wind is, regardless of that, continuous, as long as the numerator and the denominator go to zero simultaneously. For this to happen, the wind speed must reach the sonic speed at r ¼ r c : u r ðr c Þ ¼ c s . At this critical point, known as the sonic point, the numerator is written as This is actually a very important point in the theory of stellar winds. Equation (42) has an infinite number of mathematical solutions, but the only physical solution for a stellar wind is the transonic solution, i.e., the one that starts subsonic and becomes supersonic. Figure 22 illustrates solutions for the wind radial velocity, assuming isothermal winds (C ¼ 1) at four different temperatures: 0.75, 1.5, 3 and 4 MK (from bottom to top curves). I assumed here a star with 1 M and 1 R . The curves show the typical wind radial velocity profile, which is characterised by a relatively quick increase in the velocity at inner regions and then an asymptotic solution, which tends to the terminal wind velocity, u 1 . Note how sensitive u 1 is to the temperature, going from 400 km/s for a 0.75-MK wind to 1200 km/s for a 4-MK wind. The zoom-in panel in the right shows the position of the sonic point. Because r c / c À2 s / 1=T, the sonic point happens closer to the star for higher temperature winds. Likewise, because c s / T 1=2 , the speed at which the sonic point is reached is increasingly larger for higher temperatures. As discussed in page 43, the base temperature, along with the base density, are free parameters in polytropic wind models.
Using the mass conservation equation and the solutions for u r in Fig. 22, the density profiles can be found: q ¼ _ M=ð4pr 2 u r Þ. Note here that the density depends on the choice of mass-loss rate _ M, which is related to the discussion above on the density in thermally-driven wind being a scaling factor. For larger distances, when the wind has reached asymptotic speeds, q / r À2 . In this case, the ram pressure of the wind P ram ¼ qu 2 r also decreases with r À2 . This is relevant for calculating the size of the heliosphere, and how it evolved. I will discuss this in Sect. 6.3. Fig. 23 Left: initial setup of the split monopole. Right: the trailing spiral structure that is created once stellar rotation is set

The magneto-rotator wind
The thermally-driven wind discussed before does not consider two physical parameters that are observed in cool stars: rotation and magnetic fields. I am going to show now an important effect when considering magnetised, rotating winds: the removal of angular momentum from the star. Consider the initial setup shown in Fig. 23a: a split monopole magnetic field in the equator of the star. The magnetic field lines are anchored on the stellar surface, and thus, once stellar rotation is set, they must co-rotate with the star (i.e., at same rotation rate X H ). I will consider the stellar rotation axis pointing out of the page (rotation is counter-clockwise).
Magnetic field geometry: I will first derive the geometry of the magnetic field lines, using the magnetic induction equation for a perfectly conducting plasma (null resistivity) I will keep assuming steady state (o=ot ¼ 0) and axi-symmetry (o=ou ¼ 0). Given that the equation is solved in the equatorial plane, o=oh ¼ 0. Thus the magnetic induction equation becomes The wind velocity vector in spherical coordinates is ũ ¼ ½u r ; 0; u u and similarly the magnetic field is B ¼ ½B r ; 0; B u . Notice that even though I presented in Fig. 23a a split monopole (B u ¼ 0), the rotating star will generate an azimuthal component for the magnetic field in the wind (Fig. 23b). Using the curl in spherical coordinates, and the fact that ũ Â B ¼ ½B r u u À u r B u ĥ, the last equation is written as 1 r o or ðr½B r u u À u r B u ûÞ ¼ 0: This implies that the term inside the parenthesis is constant for different radial distances: rðB r u u À u r B u Þ ¼ constant. To find this constant, I analyse the wind physical properties at the wind base, where r ¼ R H . There, the wind radial velocity is very small u r ðR H Þ ' 0 and the azimuthal velocity is the rotational velocity of the star, u u ðR H Þ ¼ X H R H . The magnetic field anchored on the surface of the star is still mostly radial, thus B u ðR H Þ ' 0. Hence, the constant computed at R H is: This can be further generalised by using magnetic flux conservation: B r r 2 ¼ B r ðR H ÞR 2 H . The magnetic induction equation then becomes Rewriting this equation, I arrive at the equation describing the magnetic field geometry of the wind of a rotating magnetised star Close to the star (small r), the rotational velocity is u u ðR H Þ ¼ X H R H and thus B u ' 0. Far from the star (large r), u u ( X H R H and thus B u =B r ' ÀX H R H =u r . This implies that not only an azimuthal component of the magnetic field is generated at large distances, but that this component has opposite sign to the radial field: B u =B r \0. This is illustrated in the sketch on the right panel of Fig. 23. As a result, the rotating wind generates a trailing spiral: if B r [ 0 (pointing outward), then B u \0 (clockwise). On the contrary, if B r \0 (pointing inward), then B u [ 0 (counter-clockwise).
The momentum equation of a magneto rotator wind: The momentum equation of a magnetic rotator is similar to the equation used before, except that now the Lorentz force must also be considered where I assumed steady state and F B is the Lorentz force per unit volume where J is the current density and I used the Ampère's law (J ¼ c=ð4pÞðr Â BÞ) in the second equality. Using the curl in spherical coordinates, it is straightforward to demonstrate that r Â B ¼ Àr À1 oðrB u Þ=orĥ and thus Before substituting Eq. (51) in the momentum Eq. (49), it is necessary to write the convective derivative of ũ in spherical coordinates: Substituting Eqs. (52) and (51) into (49), I finally arrive at the momentum equation of a magneto-rotator wind. It is convenient to split it into the radial and azimuthal components, which are respectively and q u r ou u or þ u r u u r ¼ B r 4pr These are the equations derived by Weber and Davis (1967).
The radial component of the momentum equation of the magneto-rotator wind (53), compared to that of the thermally-driven wind (Eq. 32) has two additional terms: the slinging effect of the rotating magnetic field (ÀB u =½4pr o½rB u =or) and the centrifugal force (qu 2 u =r). These two terms originate from the rotation of the star. Similar to the thermally-driven wind, the magneto-rotator wind also has critical points. One important critical point for angular momentum losses is the Alfvén radius, 8 r A , which is defined as the position where the wind radial velocity reaches Alfvén velocity, v A : u r ðr A Þ ¼ v A . The Alfvén velocity is given by . This means that the wind of a magnetised, rotating star is trans-Alfvénic. In the sub-Alfvénic part of the stellar wind (closer to the star), the kinetic energy is smaller than the magnetic energy: qu 2 r =2 ( B 2 =ð8pÞ. In this case, the wind is too weak to modify the structure of magnetic field and the wind is forced to flow along the magnetic field lines. In the super-Alfvénic part of the wind, the opposite is true and qu 2 r =2 ) B 2 =ð8pÞ. In this case, the magnetic field lines are dragged with the flow. Angular momentum losses: The angular momentum carried away by the stellar wind can be found from the azimuthal component of the momentum equation of the magneto-rotator (54), which can be written as Rearranging this equation, one has where I use the facts that the mass-loss rate ( _ M ¼ 4pr 2 qu r , Eq. 34) and the magnetic flux (B r r 2 ) do not depend on r. Thus, which implies that the term in parenthesis is constant with radial distance where L is the specific angular momentum (angular momentum per unit mass). Figure 24 shows the two terms as computed for the present-day solar wind. The first term on the left hand side is the specific angular momentum carried by the gas (the one that we are most familiar with) and the second term is the angular momentum carried by magnetic stresses. We notice from Fig. 24 that even a star that has a modest magnetic field like the Sun and a relatively weak wind, the magnetic term dominates the angular momentum losses. This can be understood in light of the size of the Alfvén surface. Rewriting Eq. (58) with the help of Eq. (48), one finds Given that L is a constant, we can calculate the previous equation at r ¼ r A , in which case u r ¼ v A and the equation simplifies to This means that, if we know where the Alfvén radius is, then we can calculate the specific angular momentum of the wind by using the previous equation. I will show later on that this is however not an easy task. It is curious to see that Eq. (60) resembles so much the angular momentum of a solid body with a radius r A rotating at the angular speed of the star. Because of this, sometimes this equation is incorrectly interpreted as a wind rotating as a solid body out to r A . In fact, the wind is forced to co-rotate with the magnetic field (which is anchored on the star) out to a distance that is much smaller than r A . Finally, I can now calculate the angular momentum loss rate. The flux of angular momentum is the product of the specific angular momentum L by the mass flux qũ. Integrating the flux of angular momentum over a surface area S, one finds the angular momentum loss rate where the subscript 'A' indicates values computed at the Alfvén radius. Note that all this calculation was done for the equatorial plane, but in reality, the wind is Fig. 24 The specific angular momentum of the solar wind calculated by Weber and Davis (1967). The sum of these two terms is constant with distance and simplifies to L ¼ X H r 2 A (see Eqs. 58 and 60), where r A is the distance to the Alfvén point. Image reproduced with permission, copyright by AAS The previous equation has to be integrated to derive X H ðtÞ, i.e., the rotational evolution of the star. I already assumed that M H and R H are independent of age, and thus, independent of rotation. Two further assumptions need to be made before integrating the equation above. The first one is related on how B H depends on X H . I will assume they are related to each other in a linear way (i.e., a linear-type dynamo). This is not too far from what is observed, as I showed in Fig. 13b. The second assumption, and the most tricky one, is on how v A depends on X H . The reason why this is tricky is that v A ¼ B= ffiffiffiffiffiffiffiffi 4pq p depends on the magnetic field and density at the Alfvén point, which in turn depends on the thermal and magnetic properties of the wind. Thus, there is no way to derive the location of the Alfvén point, and its velocity, before computing a wind model. I will assume, for the sake of simplicity, that v A is independent of X H in the main sequence for tJ600 Myr, but I will come back to this in the next section. With these assumptions, the previous equation can be written as Integrating the previous expression results in t / X À2 H , arriving at the following relation which surprisingly agrees with observations! In summary, two main conclusions can be extracted from this simple evolutionary model. First, it demonstrates that rotating, magnetised stellar winds carry away angular momentum from the star. As angular momentum leaves the star, the star spins down, explaining the observed trend that stars spin down as the square-root of age. Second, angular momentum extraction is enhanced by the magnetic field: even for a star with a $ weak magnetic field as the Sun, the angular momentum carried away by magnetic stresses is the dominant contributor to the loss rate.
In a more sophisticated treatment, one should not use the moment of inertia of a spherical solid body nor the assumption of constant mass and radius through the main sequence. These quantities can be extracted from stellar evolution models and thus coupling the wind angular momentum loss prescription to stellar evolution models would provide a more sophisticated modelling for X H ðtÞ. The other naive assumption made above refers to the velocity at the Alfvén point (or surface). I assumed that the Alfvén velocity is independent of stellar rotation. However, this is a very rough oversimplification of the problem, as v A depends on the structure of the wind, which depends on the magnetic, thermal and rotational properties of the wind. Thus, there is no way to derive the location of the Alfvén surface, and its velocity, without computing a wind model. I will explore in the next section how this can be achieved.

Magnetic braking laws: prescriptions for angular momentum-loss rate
One of the biggest challenge in stellar rotational evolution studies is on determining the size of the Alfvén surface, so that the angular momentum loss rate _ J can be calculated. Equation (62) looks deceptively simple! The Alfvén radius, in a 1D geometry, or surface, in 3D, can be interpreted as the lever arm of the wind torquethe ''position'' where the torque is applied in order to change the angular rotation of the star. In fact, there is not a unique radius where the torque is applied in a stellar Fig. 25 In 3D, the Alfvén radius takes the form of an Alfvén surface (grey surface). Because of the complex stellar surface magnetic field geometry (shown by the colour bar), the Alfvén surface can be highly asymmetric. Image reproduced with permission from Vidotto et al. (2014b), copyright by the authors thus constants that were previously assumed (e.g., a, n, etc) can be derived from numerical studies. Matt et al. (2012) performed a range of 2D simulations for solartype stars with different rotation rates X H and dipolar magnetic field strengths B dip;H . The authors found a generalised fit for their simulations that can be written in terms of constants k 1 , k 2 and m as with k 1 ¼ 1:30, k 2 ¼ 0:0506, and m ¼ 0:2177. The surface escape velocity is v esc; With all values given in cgs units, then the angular momentum loss rate has units of g cm 2 s À2 , which is an energy unit and thus _ J is frequently expressed in erg.
The angular-momentum loss prescriptions shown, for example, in Eqs. (70) and (71) can then be implemented in stellar evolution codes. Using Eq. (71) is an advantage over the use of Eq. (70), as discussed in Johnstone et al. (2015a). However, stellar evolution models still need to adjust the value of k 1 so that they can be properly calibrated for the Sun (Gallet and Bouvier 2013;Johnstone et al. 2015a;See et al. 2019a). One challenge in such prescriptions for _ J is that they require a value of _ M, which, as discussed in Sect. 2, is difficult to observationally derive. In particular, the mass-loss rates are dependent of the thermal properties of the wind. In wind models that implicitly assume heating, the thermal properties at the wind base are free parameters (see page 43). Models that assume, for example, that wind temperatures increase with stellar activity (e.g., Ó Fionnagáin and Vidotto 2018) predict that mass-loss rates are higher for young and fast rotating stars.
Would then mass-loss rates increase indefinitely with rotation? As I presented in Sect. 3, although it is believed that younger stars should have higher _ M, there is still no clear picture from observations (see Fig. 11). Recent observations, using prominences to probe stellar winds, suggest that mass-loss rates seem to keep increasing towards young ages (Figs. 8, 11)-AB Dor, for instance, is a $ 120 Myr-old star, with 350 times the present-day solar mass-loss rate, or the even younger LQ Lup ( $ 25 Myr), with _ M $ 4500 times the solar value (Jardine and Collier Cameron 2019). These findings challenge the results from the astrosphere method, which indicated that mass-loss rate might actually reduce at ages below $ 600 Myr, with p 1 UMa ( $ 500 Myr) showing half the present-day solar massloss rate (Wood et al. 2014). I believe that the answer to this should come from a combined approach between modelling and observations. 6 Implications for the evolution of the solar system planets and exoplanets In this section, I briefly discuss possible implications of the solar wind evolution on planets in the solar system and on exoplanets orbiting Sun-like stars.

The evolution of Earth's magnetosphere
One frequent argument in favour of the larger mass-loss rates of the young Sun is the lack of a substantial atmosphere in Mars. Given that Mars lacks a large-scale protective magnetic field similar to that of Earth, the interaction between the planet and the solar wind is not shielded by the planet's magnetic field. As a result, the young Martian atmosphere is believed to have been eroded by a stronger solar wind (e.g., Kulikov et al. 2007). This reasoning has been recently challenged in light of ion escape measurements in terrestrial planets (Brain et al. 2013). Escape rates of O þ in Earth, Mars and Venus are reasonably similar, yet these terrestrial planets have very different properties, such as gravity, orbital distances and the presence of a large-scale magnetic field. The fact that there are not two planets in the solar system that are identical, except for their magnetic fields, makes it very difficult to isolate how, or if, a large-scale magnetic field can be an effective atmospheric shield (Brain et al. 2013).
There are two arguments relating to whether a planetary magnetic field might, or not, protect planetary atmospheres.
In favour: The magnetosphere of a planet acts to deflect the stellar wind particles, preventing its direct interaction with planetary atmospheres and, therefore, potential atmospheric erosion. If the magnetosphere is too small, part of the atmosphere of the planet becomes exposed to the interaction with stellar wind particles. Lammer et al. (2007) suggested that, for a magnetosphere to protect the atmosphere of planets, planets are required to have magnetospheric sizes J 2 times the planet's radius. This can have important effects on the habitability of exoplanets, including terrestrial type planets orbiting inside the habitable zones of their host stars (e.g., Lammer et al. 2007;Vidotto et al. 2013;See et al. 2014;Ribas et al. 2016). Against: However, it is not only the size of the magnetosphere that matters. Depending on the orientation of the stellar wind (i.e., interplanetary) magnetic field relative to the planetary one, magnetic reconnection can occur and the interplanetary magnetic field lines can connect to auroral zones around the planet's magnetic pole. Because these field lines are open, stellar wind plasma can be channeled towards the polar caps, causing local heating that could ultimately increase atmospheric losses through polar/auroral flows (Moore and Horwitz 2007. This is potentially more important for heavier species like O þ than for lighter ones like H þ , Brain et al. 2013.) The amount of stellar wind plasma that is channeled towards the polar caps depends on the collecting area of the magnetosphere (i.e., its cross section $ pr 2 M ), which is larger for higher magnetic fields (equations will be derived below).
Thus, while a magnetosphere larger than the extent of the upper atmosphere can protect atmospheric losses, larger magnetospheres also have greater 'collecting area' for stellar wind plasma, which could lead to enhanced atmospheric losses. It is likely that these two processes coexist and that their contributions differ throughout the life of the planet. Blackman and Tarduno (2018) suggested that the competition between the strength of the stellar wind and the large collecting areas of magnetospheres indicates whether a planetary magnetic field has a protective effect on the planetary atmosphere or not. In the case of Earth, they suggested that our planet's magnetic field has played a crucial role in protecting our atmosphere. However, different solar wind characteristic in the future could reverse this protective effect, where a reduced inflow speed (related to the speed of magnetic reconnection) and larger planetary magnetospheres would conspire to generate more atmospheric losses than what an unmagnetised Earth would have presented.
In a first order approximation, we can estimate the magnetospheric stand-off distance by setting pressure balance between the incoming stellar wind and the planetary magnetic pressures. At the planet-stellar wind interaction zone, pressure balance between the stellar wind (left-hand side) and planetary magnetosphere (right-hand side) can be written as where q w and u w are the density and speed of the stellar wind at the planetary orbit, and B p;r M is the planetary magnetic field intensity at a distance r M from the planet centre. Assuming that the planet's magnetic field is dipolar, with B p;eq its surface magnetic field at the equator, the previous equation can be rewritten so that magnetospheric size of the planet is given approximately by Equation (73) shows that a large stellar wind pressure acts to reduce the size of planetary magnetospheres for a given planetary magnetisation. This is the Chapman-Ferraro equation firstly derived for Earth's magnetosphere. Note that this equation needs to be modified for exoplanets in closer orbits, as in those cases, the magnetic pressure and thermal pressure of the stellar wind, as well as ram pressure of the planetary atmosphere, might have to be incorporated in Equations (72) and (73) Vidotto and Cleary 2020).
Recently, Carolan et al. (2019) studied the evolution of Earth's magnetosphere during the solar main sequence. For that, they first performed simulations of the evolution of the solar wind using a magneto-rotator wind model (1.5D MHD). With that, they were able to model the conditions of the evolving solar wind at the orbit of the Earth. These conditions were then implemented as external boundary conditions into 3D MHD simulations of the Earth's magnetosphere. Carolan et al. (2019) showed that, as the Sun spun down and the solar wind density and velocity values at Earth's orbit evolved, the Earth's magnetosphere gradually increased following a power-law r M / X À0:27 H up until the Sun reached a rotation of X H $ 1:4X . After that (i.e., for lower rotations X H .1:4X ), their model predicted a more rapid increase of the Earth's magnetosphere with r M / X À2 H . The piece-wise evolution they predict for the Earth's magnetosphere lies on assumptions for the temperature of their stellar wind models. Carolan et al. (2019) used a simplified stellar wind model, which allowed them to predict the evolution of the solar wind for a large interval of time. More sophisticated wind models have opted for a more focused approach. For example, Sterenborg et al. (2011) used modified solar magnetograms as input for their 3D MHD simulations for how the solar wind was 3.5 Gyr ago. Their wind results were then implemented in simulations of Earth's ''paleomagnetosphere'', in which the strength and orientation of Earth's magnetic field were varied. Depending on the specific conditions of the Earth's magnetic field, Sterenborg et al. (2011) showed that the paleomagnetosphere stand-off distance could have been as small as 4.25 Earth radii.
As the solar wind interacts with Earth's magnetosphere, not only the magnetospheric stand-off distance changes, but also does the area of Earth's polar cap containing open magnetic field lines. Sterenborg et al. (2011) showed that at 3.5Ga the Earth's open-field polar cap area was larger for several conditions they explored. The only exception was for when Earth's dipolar field had a substantial tilt, or when the solar wind magnetic field was parallel to the Earth's magnetic field at the interaction zone. In this case, the Earth's magnetosphere was closed, without any open-field polar cap area. Appendix B in Carolan et al. (2019) discusses the differences between open and closed planetary magnetospheres, which can more clearly be seen in their Fig. B1.
Another focused work on the young Sun-Earth study was presented in do Nascimento et al. (2016), who developed a 3D study of the magnetohydrodynamic environment surrounding j 1 Cet, which has been recognised as a good proxy for the young Sun at the age when life arose on Earth. According to the authors, the massloss rate of j 1 Cet is about 50 times larger than the present-day Sun (see also Airapetian and Usmanov 2016). Due to this larger mass-loss rate, the ram pressure of the young solar wind impacting on the magnetosphere of the young Earth was larger than the present-day values. do  showed that the magnetospheric sizes of the young Earth should have been reduced to approximately half to a third of the value it is today ( $ 11 Earth radii). The relatively large size of the early Earth's magnetosphere may have been the reason that prevented the volatile losses from Earth and created conditions to support life (Tarduno et al. 2010;. This is opposed to the fate of a young, nonmagnetised Mars, which could have lost its atmosphere due to a stronger younger solar wind ).

Implication for exoplanets and new lessons for the solar system
There are several ways in which the solar wind evolution, and the linked activity evolution, can affect planets. One way, already mentioned above, refers to the direct interaction between the evolving solar wind and the planetary magnetosphere, in the case of magnetised planets. In the case of unmagnetised or weakly magnetised planets, the magnetic pressure on the right-hand side of Eq. (72) should be replaced by the thermal and ram pressure of the planetary atmosphere. In particular, planets that receive a large amount of high-energy irradiation can have evaporating atmospheres. This happens in either planets that orbit old stars at close distances, or planets with larger semi-major axis orbiting more active (young) stars. These two conditions lead to enhanced XUV energy flux deposition on the top layers of planet atmospheres, which leads to an increased heating, causing atmospheres to inflate and to more easily evaporate (e.g., Lammer et al. 2003;Baraffe et al. 2004;Sanz-Forcada et al. 2011). Since the XUV emission of solar-type stars is observed to evolve (Sect. 4.3), exoplanetary evaporation is also expected to evolve in time.
This implies that the history of the host star plays an essential role in the evolution of atmospheric evaporation of their planets. For example, a young host star rotating more slowly than other stars at the same age would generate less EUV radiation and thus their planets would evaporate less and end up with a higher atmospheric mass fraction. On the contrary, a young host star that is a fast rotator would emit more EUV radiation, which could potentially lead to a planet with a reduced atmosphere. Tu et al. (2015) showed that a 0.5-Earth-mass planet orbiting at 1 au around a solar-mass star will present an atmosphere with a very different hydrogen content at 4.5 Gyr, depending on whether the host star was a slow or fast rotator during its youth. Starting with an initial hydrogen atmosphere of 0.005 Earth mass, they showed that the atmosphere of the terrestrial planet is entirely lost before an age of 100 Myr if the host star begins its main-sequence evolution as a rapid rotator, while if the host star was a slow rotator, at 4.5 Gyr, the planet would still have about 45% of its initial atmosphere.
Because the 'final' (i.e., observed) planetary atmospheric fraction depends on the EUV history of the star, Kubyshkina et al. (2019a) suggested that detecting the atmospheres of planets today can provide a way to probe the evolutionary path of the host star. Kubyshkina et al. (2019b) applied this model to the Kepler-11 system, formed by a star that is very similar to the Sun (in mass and age) and that has six known transiting planets with orbital distances between 0.09 and 0.5 au. Given that these six close-in planets seem to have retained hydrogen-dominated atmospheres, these authors concluded that this Sun-like star was a slow rotator at its youth. In a similar reasoning, models using instead the Earth and Venus atmospheric composition were conducted by Lammer et al. (2020), who showed that one possible explanation for the noble gas composition in Earth's and Venus' atmospheres could be due to the Sun being a slow to moderate rotator at young ages.
Works like these demonstrate that modelling star-exoplanet interactions could open up new avenues for us to progress in our understanding of the solar system, and the Sun and solar wind evolutions.
Another important consideration for the young solar system is that the wind of the young Sun might have had an increased rate of CMEs and even might have been CME-dominated at an early point in time (see Sect. 2.5). Although CME rates decrease with age, due to geometrical effects, close-in planets interact much more frequently with stellar CMEs than planets at larger orbits (Kay et al. 2016). The CME itself might also be highly affected by the presence of a close-in planet, given that the CME does not have enough time to expand before it encounters the close-in planet (Cohen et al. 2011b). CMEs cause perturbations in stellar wind conditions, which can then lead to changes in the magnetospheres and extended atmospheres of exoplanets Cherenkov et al. 2017). CME-planet interaction might be more common in planets orbiting M dwarfs, as these stars remain active for a long part of their lives (West et al. 2008;Irwin et al. 2011) and have high flare and CME rates (Vida et al. 2016(Vida et al. , 2017. M dwarf stars are the prime targets for detecting terrestrial planets in their habitable zone. However, the consequence of their high rate of energetic events might affect habitability of M dwarf planets . Although stellar winds constantly interact with exoplanets, exoplanets do not interact with every stellar CME, because the planet may not lie along certain CME trajectories. In the Sun, CMEs are observed within AE 30 latitude, corresponding to the active region belt (Gopalswamy et al. 2019). However, in young, fast rotating stars, active regions are often located at high, polar-region latitudes (Marsden et al. 2006;Jeffers et al. 2007). If CMEs originate in active regions, it is therefore more likely that CMEs in young stars could emerge at high latitudes and be directed above/below the equatorial plane. Considering planet formation in the equatorial plane, this could mean that young planets might at the end not be strongly affected by CMEs, even though their host stars are expected to have higher CME rates. Future research in this area might shed more light into the effects of CMEs in young planetary systems.
Regardless of whether exoplanets lie along CME trajectories, they can nevertheless be affected by the almost-instantaneous increase in irradiation caused by an eventual CME-associated flare, since this emission covers a wider solid angle. An increase in irradiation input caused by a flare can lead to increase in evaporation in close-in gas giants (Cherenkov et al. 2017;Bisikalo et al. 2018;Hazra et al. 2020). Transit observations of the hot Jupiter HD189733b showed an enhancement of atmospheric evaporation that took place 8h after the host star flared. Lecavelier des Etangs et al. (2012) suggested that the violent event and the enhanced evaporation are potentially related to each other. The remaining open question is whether it is the increase in XUV emission caused by the flare or the interaction between an unseen associated CME and the exoplanet (or even both) that caused the increase in the observed evaporation rate.

Evolution of the heliospheric size and implications for the penetration of Galactic cosmic rays
The solar wind expands into the interstellar medium (ISM), giving rise to the heliosphere. Analogously to the Chapman-Ferraro equation for calculating the Earth's magnetospheric stand-off distance, one can use pressure balance between the solar wind and ISM to derive the size of the heliosphere. The ram pressure of the solar wind can be rewritten as P ram ¼ qu 2 ¼ _ Mu 1 =ð4pr 2 Þ, where I assumed that for larger r, the wind has reached terminal speed u 1 . Thus, pressure balance leads to _ Mu 1 4pr 2 ¼ q ISM u 2 where the subscript 'ISM' refer to ISM 'wind' quantities. If one assumes that the ISM density and velocity have not changed significantly over the past 4.5 Gyr, the larger mass-loss rates of the young Sun implies that the heliosphere was larger in the past. Rodgers-Lee et al. (2020) estimated that heliospheric size of the solar wind extended out to 950 au when the Sun had an age of 1 Gyr and could have ranged between 1300 and 1700 au at t ¼ 600 Myr, i.e., more than a factor of ten larger than the present-day value of 122 au. The size of the heliosphere, as well as the solar wind conditions, have consequences for the propagation of cosmic rays. Galactic cosmic rays are suppressed as they travel through the solar wind all the way to the Earth-this ''suppression'' is usually referred to as cosmic ray ''modulation''. I refer the reader to Potgieter (2013) for a comprehensive review on the modulation of cosmic rays in the heliosphere. In the framework of a diffusive transport equation for the Galactic cosmic rays, there are three main competing processes in place: the diffusion of cosmic rays into the solar system, and their advection, in momentum and in space, which acts to suppress the flux of Galactic cosmic rays reaching Earth. These processes depend on the level of turbulence present in the solar wind magnetic field and on the wind velocity, both of which evolve in time.
Using the solar wind evolution model of Carolan et al. (2019), Rodgers-Lee et al. (2020) showed that, at the Earth's orbit, the advective processes (pushing cosmic rays out) are relatively more important than diffusive processes (allowing them in) in the Sun's past, for GeV and MeV cosmic rays. This means that cosmic rays of these energies are more suppressed and thus fewer of those have reached Earth in the past, compared to present-day values. This is seen in Fig. 26, which shows the modelled evolution of Galactic cosmic ray spectrum at Earth for a range of solar rotation rates. The black solid line is the cosmic ray spectrum at the local interstellar medium (Vos and Potgieter 2015). The black dashed-line is the present-day value (which was calibrated to present-day measurements) at Earth's orbit, at approximately solar minimum. The red lines and shaded area represent a range of rotation values that the Sun could have had at 600 Myr, i.e., the predicted flux of Galactic cosmic rays can change by more than one order of magnitude depending on whether the Sun was a fast/slow rotator. Rodgers-Lee et al. (2020) considered a 1D propagation model for cosmic rays. An interesting future study would be to use a 3D cosmic ray propagation model, in 3D wind simulations that include ZDI observations, similar to the approach adopted in Cohen et al. (2012).
With the knowledge of the spectrum of cosmic rays arriving at the top of Earth's atmosphere, further modelling can then be conducted to investigate variations of high-energy particle fluxes into the young Earth's atmosphere (e.g., Grenfell et al. 2007;Stadelmann et al. 2010;Rimmer and Helling 2013;Grießmeier et al. 2016). Interestingly, it has been suggested that the variation of cosmic ray flux over the past 3 billion years coincide with periods of glaciations at Earth (Svensmark 2006, see also Shaviv 2002). Another interesting application, not discussed in this review, is to investigate cosmic ray propagation through the winds of M dwarf stars (Grießmeier et al. 2005(Grießmeier et al. , 2015. These stars are currently among the main targets in search for habitable planets.

Final remarks and open questions
Over the years, we have collected a large amount of information on the solar wind, in an incredible level of detail. All this information, however, can only tell us about how the solar wind looks like now. To understand the past, and future, evolution of the solar wind, we rely on observations from other suns in the Universe, i.e., stars that resemble the Sun and are at different evolutionary stages. In theory, connecting the observations of winds of different solar-like stars under one unique evolutionary sequence would tell us how the wind of our own star has evolved. In practice, however, this is not straightforward.
One reason for this is that detecting winds of solar-like stars is a very challenging task. Direct detection of winds of other Suns has not been possible and the Sun still remains the only star for which the wind has been probed in situ. However, over the years, a small number of clever techniques have been developed to indirectly derive stellar wind properties, such as mass-loss rates and speeds. Overall, these techniques Fig. 26 Evolution of the Galactic cosmic ray spectrum at Earth's orbit. The figure above shows Galactic cosmic ray differential intensity as a function of their kinetic energies. The different curves refer to different solar rotation rates normalised to present-day value (here denoted by X 0 ). Each curve thus represents a different age, which ranges from about 6-1 Gyr from the top to the bottom dashed lines. The black dashed line is thus for the present-day value (1X 0 ) at Earth's orbit. The black solid line is the cosmic ray spectrum at the local interstellar medium. Cosmic rays are suppressed (i.e., modulated) as they travel through the solar wind and, given the solar wind evolution, the spectrum of Galactic cosmic rays at Earth changes with time. The red lines and red shaded area represent a range of rotation values that the Sun could have had at 600 Myr. Image reproduced with permission from Rodgers-Lee et al. (2020), copyright by the authors have shown that the solar wind had a higher mass-loss rate for the past $ 3.5-4 billion years. Before that, from the time the Sun started its evolution in the main sequence until about 600-800 Myr, clues point in different directions.
It has been suggested that the young solar wind had a lower mass loss rate, probed in the lower-than-expected mass-loss rate of p 1 UMa, a 500-Myr-old solar analogue, which has only half the present-day solar mass-loss rate _ M (Wood et al. 2014). In particular, p 1 UMa does not fit the trend of _ M-F X relation built from the astrosphere wind detection method, that older solar-like stars seem to fall on. However, a recent method for detecting winds of rapidly rotating stars, using prominences to probe the wind mass-loss rate, showed that young Suns could have much higher mass-loss rates, with the solar analogue AB Dor ( $ 120 Myr) presenting a mass-loss rate of 350 _ M and the even younger LQ Lup ( $ 25 Myr), with 4500 _ M (Jardine and Collier Cameron 2019). One way to shed light on the solar wind evolution is to look at physical characteristics known to be intimately related to winds of solar-like stars: rotation, magnetism and activity. These three 'ingredients' are related to the evolution of the solar wind in a feedback loop that I summarised in Fig. 1. The stellar wind carries away angular momentum, which leads stars to spin down. A decrease in rotation rate implies that the dynamo weakens, generating weaker magnetic fields. This affects the stellar wind and the amount of angular momentum it can carry away. Then, the loop restarts. As a result, rotation, magnetism and its proxy in the form of stellar activity decrease with time.
These ingredients, which are better observationally constrained than wind properties, can be implemented in theoretical models. Models of winds of solar-like stars have been mainly developed in the magnetohydrodynamics (MHD) framework, which treats the wind as a magnetised fluid. Within this framework, models can be broadly divided into ones that calculate explicitly heating deposition (due to dissipation of MHD waves mostly) and those that assume the plasma has already been heated to MK temperatures, and thus energy deposition is treated implicitly. Each model set has its advantages and drawbacks. The first group has a more sophisticated treatment of the physical wind driving mechanism, but is usually computed along a magnetic flux tube, neglects rotation, and extend only out to a few solar radii. The second group has a less sophisticated treatment of the physics, which implies that they can be computed over the entire stellar surface, extend out to larger stellar distances and include rotation.
For studying the evolution of angular momentum, the second group of models are usually adopted. In that framework, it is fascinating to see that the early model of Weber and Davis (1967), which considers a simple split-monopole magnetic field embed in a thermally-driven wind of a rotating Sun already does a good job in reproducing the observed relation X H / t À1=2 for older stars (Sect. 5.3). In this sense, stellar evolution models, calibrated to the observed rotation distribution of stars in open clusters, have become an important aid for constraining the solar wind evolution.
In this review, I focused on the solar wind evolution from the early main sequence until today. That does not mean that what lies ahead is less interesting! The Sun is about half-way through its main-sequence journey, and will have thus another $ 4 Gyr worth of evolution in a relatively calm way. When the hydrogen fuel in its core is finally exhausted, the Sun will evolve off the main sequence and become a red giant. The timescales of this phase and subsequent ones become shorter and thus the Sun will pass through these phases very quickly in comparison to the long evolution during its main sequence. Eventually, the Sun will end its life as a white dwarf and cool down indefinitely. Certainly, the solar wind will change radically. In the red giant phase alone, the solar wind will likely be cooler (. 10 4 K), with lower velocities, but higher densities, leading to mass-loss rate that will be several magnitudes higher than today (Wood et al. 2016). Solar rotation will change considerably due to its radius expansion. Given that rotation also seems to be related to magnetism in evolved stars (Aurière et al. 2015), the Sun's magnetic field will change as well.
The evolution of the solar wind has applications that goes beyond understanding stellar rotational evolution-it affects planets and exoplanets, which at the end of the day are embedded in the outward-streaming plasma from their host stars. A strong solar wind in the past is often cited to explain the lack of a substantial atmosphere in Mars today-planetary atmospheric erosion can be enhanced under strong stellar wind conditions. Such strong stellar wind conditions are also faced by close-in exoplanets, even when orbiting old and less active stars. Additionally, the Earth's magnetosphere is believed to have been smaller than it is today, as a result of a strong stellar wind. Similarly, the strong solar wind in the past was able to push further the boundary with the ISM-the heliosphere was likely larger in the past. A different solar wind condition and much more expanded heliosphere would then have likely reduced the amount of Galactic cosmic rays reaching Earth. Both of these (Earth's magnetosphere and cosmic ray flux) are believed to be important factors for life formation and sustainability in this planet.
I hope to have shown in this review that significant progress in understanding the evolution of the solar wind is being made. Nevertheless this is currently a ''living'' area of research and, as such, still many open questions remain. Understanding the physical processes heating and driving the present-day solar wind could provide a step further in modelling the solar wind at different ages. Missions like NASA's Parker Solar Probe and ESA's Solar Orbiter will obtain in situ measurements of the solar wind at unprecedented close heliospheric distances. These regions are where the bulk of the acceleration of the solar wind takes place, and thus these missions are expected to soon shed light on the physical mechanisms that accelerate the solar wind. As an added bonus, these close heliospheric distances coincide with the regions where many close-in exoplanets are located-these missions will thus provide information about the environment at the orbits of close-in exoplanets, which will aid in further understanding how extreme wind conditions can affect exoplanets, and potentially, link this to how the stronger solar wind in the past interacted with the young solar system planets.