A Headwind-like Helium Outflow from TOI-3884b

Corresponding author: Joe P. Ninan

Introduction

Atmospheric escape drives the evolution and composition of planetary atmospheres, from primordial atmospheres to the final atmospheric composition of planets that we observe. Direct measurement of this phenomenon is crucial for understanding how planetary atmospheres evolve, shaping the observed properties of planets and their demographics. To measure the extended but tenuous exospheres of planets via transmission spectroscopy, one needs to select absorption lines with the largest cross sections. Since neutral hydrogen is abundant in giant planetary atmospheres, Ly$\alpha$ provides the strongest absorption signal. Vidal-Madjar et al. [53] detected exospheric outflow from HD 209458b via excess absorption in the Ly$\alpha$ wings. A few other planetary exospheres have also been detected using Ly$\alpha$ ([23, 22, 16, 11, 15, 25]). While Ly$\alpha$ is the strongest absorber, it suffers from the drawback that the interstellar medium is also full of neutral hydrogen. This, combined with geocoronal contamination, completely absorbs the zero-velocity component of Ly$\alpha$, leaving only the high-velocity wings detectable. Seager & Sasselov [42]; Turner et al. [52]; Oklopčić & Hirata [35] identified the He I 10830 Å absorption line of metastable helium as a potential alternative. If the conditions allow the formation of metastable helium atoms, the absorption line formed by the resonant scattering of 10830 Å photon is ideal, because it is not affected by the interstellar medium and has minimal contamination from the Earth’s atmosphere. This implies that ground-based, high-resolution spectrographs on large telescopes can detect this signature. Initial measurements of the He I 10830 Å signature were made around F-, G-, and K-type star planets ([46, 40, 28, 2, 34, 5] (Allart et al. 2019)). For these spectral types, the stellar XUV flux is sufficient to ionize helium and to form metastable helium atoms during the recombination cascade through triplet states. While the XUV flux from a late-type M-dwarf photosphere is not sufficient to ionize helium, the active chromospheres of some M-dwarfs produce enough XUV flux. The first detection of a planetary helium exosphere around an M-dwarf (GJ 3470b) was reported by Ninan et al. [33]; Palle et al. [37]. These detections enabled measurements of planetary outflow rates around stars beyond the traditional FGK-type hosts.

The Habitable-zone Planet Finder (HPF) has a long-running campaign to detect these exospheres ([33, 56, 18, 36]). TOI-3884b is a super-Neptune in a polar orbit around an active M4 dwarf with persistent spots ([4, 24]). Since it has a low bulk density of 0.67 g cm$^{-3}$, it is a good candidate for atmospheric mass loss. The transits of TOI-3884 show persistent starspot crossings, and several studies have confirmed that this spot complex lies near the stellar pole ([24, 50, 30, 12]), and is partially occulted by the planet in the first half of its transit. Being fully convective and a low-mass star, its photosphere is not hot enough to produce sufficient XUV radiation to ionize helium. However, XUV radiation from active regions and the chromosphere could potentially produce enough metastable helium atoms to enable detection of the outflow. Irrespective of whether the outflow is core-powered or photoevaporation-driven, measurement of the outflow morphology can help constrain the mass-loss mechanisms in these late M-dwarf systems. In § 2 we outline our in- and out-of-transit observations of TOI-3884b and the data reduction. § 3 discusses the detection and modelling of the observed outflow rate. Finally, in § 4 we summarise the main results and provide insights into what is driving the outflow in these late M-dwarf systems.

Observations and Data Reduction

We observed TOI-3884 with the Habitable Zone Planet Finder (HPF) spectrograph on the 10-m Hobby-Eberly Telescope (HET) at McDonald Observatory, USA. HPF is an ultra-stabilized extreme precision radial velocity spectrograph, operating in near-infrared (8079-12786 Å) with a resolving power of $R \sim55{,}000$ ([27, 29, 47]). HET has a fixed-altitude primary mirror and a movable prime-focus unit that tracks the star. This limits the observable altitude range of HET to 55 ± 8 degrees ([38, 44]). For TOI-3884 at declination +12.5 degrees, this implies a continuous maximum track length of 4080 seconds. We observed TOI-3884 across multiple nights spread over 4 years. We obtained 49 in-transit spectra across 5 transits, and 64 out-of-transit spectra randomly spread across the out-of-transit phase. Table 1 lists our in-transit and out-of-transit observations of TOI-3884b.

Observation JD	Exposure Time (s)	SNR	In/Out Tranist
2459550.990452	969.1	64	Out
2459551.001865	969.1	73	Out
2459564.962886	969.1	60	Out
2459564.973614	969.1	65	Out
2459569.946142	969.1	54	Out
2459569.957349	969.1	60	Out
2459570.940155	969.1	52	Out
2459570.951242	969.1	37	Out
2459571.942175	969.1	63	Out
2459571.952683	969.1	61	Out
2459575.932589	969.1	37	Out
2459575.943844	969.1	65	Out
2459596.868282	969.1	49	Out
2459596.88009	969.1	69	Out
2459597.877765	969.1	72	Out
2459597.889182	969.1	72	Out
2459601.868763	969.1	60	In
2459601.88035	969.1	51	In
...	...	...	...

We optimized our scheduling to preferentially observe when the observer barycentric velocity shifts the telluric OH emission line and the water absorption line near 10833Å away from the stellar and planetary rest frame position of He I 10830Å.

The observed data were reduced using the standard PSU HPF pipeline. Briefly, the steps include conversion from raw up-the-ramp data to slope images (after bias-noise, and nonlinearity correction) using HxRGProc ([32]). The slope image contains the photoelectrons collected in each pixel per second, as well as an estimated uncertainty for each pixel value. This error includes both detector and photon noise, and it is propagated through all the subsequent steps of the pipeline. The spectra are extracted using the QTH lamp continuum spectrum in the flat-relative optimal extraction algorithm ([55]). The extracted spectra are wavelength calibrated using Laser Frequency Comb (LFC) spectra obtained throughout the night. This wavelength calibration procedure is detailed in Stefansson et al. [48]. The telluric absorption in the spectra were corrected using a grid of telluric models generated with LBLRTM ([13]) interpolated to each frame’s best fit precipitable water vapor column and convolved with an accurate representation of HPF’s instrumental line spread function (Krolikowski et al. in prep.). The simultaneous sky spectrum from the sky fiber is interpolated onto a common wavelength grid and subtracted from the stellar science spectra. This subtraction is not perfect and leaves some residuals in the data. The OH sky-emission line strengths are correlated across the sky spectrum. Hence, we could remove the residuals from sky-subtraction, without overfitting the OH lines near 10830 Å by modelling all the residuals of the sky-emission subtraction across the echelle order containing the 10830 Å. We first created a stellar template spectrum by averaging all the spectra in barycentric velocity frame. This stellar spectrum was shifted to the telluric frame and subtracted from each epoch to obtain the residuals from the sky-subtraction. We performed a Principal Component Analysis (PCA) to obtain eigenvectors that represent the correlated residual structures in the OH sky-emission lines across the echelle order. Based on the PCA explained variance ratio, we selected six principal components. These vectors were linearly fit to each spectrum to remove residuals from the OH sky-emission line subtraction across the echelle order. In addition, in the observations taken just before sunrise (after 3:30 am CST, 9:30 am UTC), we found an excess He I 10830Å emission signature in the telluric rest frame. This emission originates from the Earth’s exosphere ([43]). We inflated the error bars in this narrow emission-line region in Earth’s rest frame for these epochs to render any residual signal from this feature statistically insignificant.

To generate an out-of-transit template spectrum, we used all spectra taken outside a 7-hour window around the center of transit. This conservative out-of-transit window was chosen to avoid any potential helium escape beyond the white light transit duration of TOI-3884b (1.58 hours). All out-of-transit spectra were averaged together, weighted by the variance, to obtain a high-S/N out-of-transit template spectrum. This template spectrum was then used to divide each spectrum, yielding the ratio of in-transit to out-of-transit, i.e., the transmission spectra. A second-order polynomial across the echelle order containing 10830Å was used to normalize epoch-to-epoch variations in the continuum that could arise from instrumental effects and changing sky conditions.

Since our in-transit spectra were observed across multiple transits under different sky conditions, barycentric velocity, and S/N ratios, we used the following procedure to measure the excess transit depth for each epoch. First, all in-transit spectra were averaged together in planet’s reference frame, weighted by the variance, to obtain the in-transit/out-of-transit He I 10830Å transmission signature. This is shown in Figure 1. This average profile was then fit, using a weighted least-squares method (with pixel-by-pixel uncertainties) on individual epoch spectrum, to obtain the peak transit depth at each epoch. Since the shape of the absorption profile is fixed, the area inside the curve (i.e., equivalent width) also scales proportional to the transit depth. This procedure provides a robust and efficient estimate of the transit depth and corresponding equivalent widths for heteroscedastic transmission spectra. One major drawback of this technique is that it assumes the absorption profile shape is not variable and is the same across all transit epochs. Within the uncertainties of our single-epoch-spectra we do not see any statistically significant variations. The TOI-3884b orbital and planetary parameters used in this work are summarized in Table 2. The transit depth and corresponding scaled equivalent width for all epochs are shown in Figure 2.

Parameter	Value	Description	Reference
$\gamma$	8.3444 $\pm$ 0.0015	Absolute stellar RV (km s$^{-1}$)	Almenara et al. [4]
$M_{*}$	0.298 $\pm$0.018	Stellar mass (M$_{\odot}$)	Libby-Roberts et al. [24]
$T_{0}$	2459556.51661753 $\pm$0.00025	Transit midpoint (BJD$_{\mathrm{TDB}}$)	Libby-Roberts et al. [24]
$T_{14}$	0.06609 $\pm$0.002	Transit duration (days)	Libby-Roberts et al. [24]
$P$	4.54458889 $\pm$0.0000098	Period (days)	Libby-Roberts et al. [24]
$e$	0.06 $\pm$0.05	Eccentricity	Libby-Roberts et al. [24]
$\omega$	$-1.95^{+4.28}_{-0.04}$	Argument of periastron (radians)	Libby-Roberts et al. [24]
$K$	$28.03^{+6.06}_{-6.23}$	RV semi amplitude (m s$^{-1}$)	Libby-Roberts et al. [24]
$M_{\mathrm{pl}}$	$32.59^{+7.31}_{-7.38}$	Mass of the planet (M$_{\oplus}$)	Libby-Roberts et al. [24]
$R_{\mathrm{pl}}/R_{*}$	0.197$\pm$0.002	Scaled radius of the planet (R$_{*}$)	Libby-Roberts et al. [24]
$b$	$0.089^{+0.082}_{-0.061}$	Impact Parameter	Libby-Roberts et al. [24]
a	0.0361 $\pm$0.0008	Semimajor axis (au)	Libby-Roberts et al. [24]

Magnetic field measurement

To model the photospheric magnetic field, we modeled a set of eight Ti I lines in the multi-epoch-averaged HPF spectra of TOI-3884 (with a combined signal to noise ratio of ~300). The details of the modeling are **Note—**This table is available in its entirety in machine-readable form online.

described in Getman et al. [17]; Amateis et al. [6]. Briefly, the spectral synthesis of these magnetically sensitive lines was performed by including the effects of Zeeman broadening and polarized radiative transfer. The observed spectra were modeled as a linear combination of synthetic models with field strengths of 0 kG, 2 kG, 4 kG, and 6 kG, with variable filling factors. Figure 3 shows the inferred filling factors for different magnetic field strengths. The average field strength obtained for TOI-3884 is $\langle B \rangle= 1.40^{+0.17}_{-0.18}$ kG. The small increase in the filling factor for 6 kG relative to 4 kG (see Figure 3) might be due to the large active spot region.

Results and Discussion

In-transit Signal

We detect a strong He I 10830Å signature, with an excess transit depth of 6.8$\pm$0.27% and an equivalent width of 0.101$\pm$0.004Å on weighted average across the five transits we observed. The width of the profile (see Figure 1), combined with the relatively small change in the planet’s velocity during transit, makes it difficult to estimate any velocity variation during the transit. As detailed in § 3.6, our confidence in the in-transit signal from the planet is strengthened by the comparable levels of pre-transit signal seen in this system. We detect significant variations in the strength of the He I 10830Å signal across the five transits. While the signal was present in all the transits, its amplitude varied from transit to transit (see the scatter in the inset of Figure 2). Further follow-up observations are needed to determine whether this variability is due to the transit light source effect (TLSE) or to fluctuations in the exoplanet’s exospheric outflow.

1D model

To constrain the mass outflow rate, we modeled the exospheric outflow using an isothermal 1D Parker wind. While this model cannot explain the pre-transit signature, its simplicity allows us to constrain the mass outflow rate using the equivalent-width measurement during transit. The modeling of metastable helium atoms and the line radiative transfer were carried out following the procedure developed in Ninan et al. [33]; Oklopčić & Hirata [35]. TOI-3884, the host star, is an M4 dwarf star ($P_{\rm rot} \sim11$ days; Mori et al. [30]). The XUV spectrum irradiating the planet is one of the most critical inputs for modeling the exosphere. For instance, our model using the XUV spectrum of another M4 dwarf star, GJ876, from the X-exoplanets archive (Sanz-Forcada et al. [41]) produced no He I 10830Å signature. The XUV spectrum of the moderate to low activity GJ 876 ($P_{\rm rot} \sim97$ days; Rivera et al. [39]) was insufficient to ionize helium rapidly enough to drive an outflow with any detectable He I 10830Å signature. Since TOI-3884 is an active star with a persistent spot (Libby-Roberts et al. [24]), we instead adopted the XUV spectrum of an active M4 dwarf star, YZ CMi ($P_{\rm rot} \sim2.77$ days; Donati et al. [14]), to irradiate the planetary atmosphere in our model (estimated from FUV spectroscopy with HST; GO 17428, PI - K. France). YZ CMi has a $\log(R'_{HK}) = -4.078$ (Astudillo-Defru et al. [10]), while TOI-3884 has a $\log(R'_{HK}) = -4.79 \pm0.13$ (Almenara et al. [4]). This model produced a significant outflow and a He I 10830Å signature consistent with the observations. Since the Parker-wind-based 1D model has two parameters, $\dot{M}$ for the mass-loss rate and $T$ for the temperature of the upper exosphere, we compared models across a grid using the equivalent width measurement. Figure 4(a) shows the chi-square valley in the $\dot{M}$-$T$ space that matches the observed equivalent width. The right panel of Figure 4(b) shows how the modeled profile compares with the observed data. Within the caveats of this 1D model, TOI-3884b has an evaporating atmosphere with a mass-loss rate of $\sim1-4 \times10^{10}$ g s$^{-1}$ for a typical isothermal model with $T = 4000-6000$ K.

Pre-transit

In two epochs when we observed before the optical transit started, we see a clear pre-transit He I 10830Å absorption signature in the data (see Figure 2). However, we do not see any indication of a post-transit cometary tail signature. As far as we know, the only other system that has shown this behavior is HAT-P-67b (Gully-Santiago et al. [18]). The host star of HAT-P-67b is an F5 subgiant, and the 3D modeling of HAT-P-67b (Nail et al. [31]) suggests that this kind of dayside-only outflow can occur if the system is sufficiently cool and the exosphere easily overflows the Roche lobe from the dayside. It is interesting that we see a very similar outflow structure in this planet, which orbits an M4 dwarf star. Another system with a longer pre-transit signal than post-transit tail is HAT-P-32b (Zhang et al. [56]). In all three cases, the absorption signal is strong: 10% in HAT-P-67b, 8.2% in HAT-P-32b, and 6.8% in TOI-3884b. However, the mass-loss rate in TOI-3884b is three orders of magnitude smaller. Since the inactive photosphere and chromosphere of an M4 dwarf does not produce enough XUV flux to drive an outflow, we hypothesize that the persistent active regions on this host as well as an active chromosphere could be driving this dayside-only outflow. Since a core-powered mass outflow is symmetric, this strong asymmetry detected between the pre- and post-transit signals further supports irradiation driven mass loss over core-powered mass outflow. The pre-transit signature starts $\sim0.25$ days before the optical transit. This implies that the headwind stream extends out to about 23 times the planetary diameter (7.8 times the Roche lobe radius).

Bubble versus Stream morphology

MacLeod et al. [26] proposed a classification of exoplanetary outflow morphology into bubble and stream systems. When the planet has a hot (i.e., faster) outflow, the wind escapes symmetrically and forms a bubble morphology, resulting in a nearly isotropic absorption signature during both pre- and post-transit observations. When the planet has a colder (i.e., slower) outflow, the wind is affected by the two-body gravitational potential beyond the Hill radius, and it primarily escapes as a stream via the L1 and L2 Lagrange points. Stream morphology can lead to the detection of extended pre- or post-transit tails. This morphological distinction can be quantified by the Rossby number. The Roche lobe radius of TOI-3884b is $r_{\mathrm{h}} = 5.8R_{\mathrm{p}}$. For our best-fit 1D Parker wind model, the outflow velocity at the Roche lobe is $v_{\mathrm{h}} = 7.5\ \mathrm{km\ s}^{-1}$ and the sound speed is $c_{\mathrm{s}} = v_{\mathrm{s}} = 7.9\ \mathrm{km\ s}^{-1}$. This implies a hydrodynamic escape parameter of $\lambda_{\mathrm{H}} = 0.9$ and a Rossby number of $Ro_{\mathrm{H}} = 1.8$. In the idealized MacLeod et al. [26] model, a Rossby number greater than 1 implies a bubble morphology rather than a stream morphology. However, similar to the case of HAT-P-67b, despite the slightly higher Rossby number, the pre-transit absorption indicates a stream-like structure.

Magnetically channeled outflow

When a planet orbits inside a star’s Alfvén radius, that is, where the Alfvén speed exceeds the stellar wind speed, the planet will not form a magnetic bow shock; instead, it threads through the stellar magnetic field and forms magnetic-field line channels to the star (Strugarek & Shkolnik [49]). The exosphere outflow can then be channeled along open field lines toward the star through the polar regions (Adams [1]).

Since the dipole dominates at large distances, we assume a $1/r^{3}$ radial falloff of the field strength starting with $\langle B \rangle= 1.40^{+0.17}_{-0.18}\ \mathrm{kG}$ at the photosphere estimated in § 2.1. This implies a field strength of $B = 8 \times10^{-5}\ \mathrm{kG}$ at the orbital distance of 25.9 $R_{\odot}$ for TOI-3884b. The density of the stellar wind depends on the mass outflow of the host star. Mass outflow rates

from M-dwarf stars are poorly constrained. From mass conservation arguments, the density scales as the inverse square of radius $1/r^2$, and inversely with velocity $1/v$. In an isothermal Parker wind or Alfvén wave driven wind, the velocity increase is nominal outside the critical point. We approximate the $1/r^2$ regime in the outer regions of the stellar wind density profile presented in Figure 2 of Shoda et al. [45] for Sun-like stars, with $\sim1 \times10^{-18}/r^2\ g\ cm^{-3}$, where $r$ is the orbital radius in units of stellar radii. This implies a mass density of $\rho\sim1.49 \times10^{-21}\ g\ cm^{-3}$ in the stellar wind at TOI-3884b’s orbital distance. The stellar wind velocity at this distance is $\sim400\ km\ s^{-1}$ [45]. The Alfvén speed is given by the formula $\frac{B(r)}{\sqrt{\mu_0\rho(r)}}$, where $\mu_0$ is the permeability of the vacuum. Hence, at the planet’s location, the Alfvén speed is $\sim5887\ km\ s^{-1}$. This Alfvén speed is significantly greater than the stellar wind speed at TOI-3884b, and therefore we expect sub-Alfvénic interaction between the planet and the star. Since TOI-3884b is in a polar orbit (see Figure 9 in Libby-Roberts et al. [24]), the sub-Alfvénic interaction opens up an alternate hypothesis for the pre-transit signature. It could be a magnetically channeled outflow. The planet’s magnetic fields could connect to the polar region of the stellar hemisphere visible to us. Before transit begins, this funnel connecting the planet to the polar region could be channeling the outflow whose signature we detect in He I 10830 Å. After the transit, this channel will be weaker as the planet moves away from this pole, and toward the far hemisphere pole behind the star, resulting in the lack of a post-transit tail. A complete MHD simulation of this interaction is beyond the scope of this manuscript, but it remains a viable hypothesis to explore.

Stellar activity contamination

Since TOI-3884b transits consistently over the active spot region [24], the observations could be susceptible to the transit light source effect (TLSE), i.e., if there are bright He I 10830Å-emitting regions around these stellar spots that are selectively blocked by the planet during transit, it can result in an extra absorption dip that masquerades as the planetary signal. We could not rule out the TLSE using velocity variations during transit, because of the small expected velocity shifts and the insufficient S/N. If we assume that the active regions correspond to the three spots estimated in Libby-Roberts et al. [24], then for the planet to produce a 6.8% extra transit depth, the excess He I emission per unit area in the active regions must be at least 4 times that of the rest of the stellar disc. We can confidently rule out the hypothesis that the detected signal is entirely due to TLSE, since that scenario cannot explain the strong pre-transit absorption signal. However, we cannot rule out some level of stellar contamination in the amplitude of the in-transit signal or its variations reported in § 3.1. One of the out-of-transit outliers in Figure 2, epoch marked by E, can be explained by a strong stellar activity event.

At that epoch spectrum, we could clearly see a change in stellar rest frame He 10830 Å line strength.

Conclusion

We report the detection of a strong exospheric outflow from TOI-3884b, despite its host star being an M4 dwarf. In our isothermal 4000 - 6000 K model, the planet is losing mass at a rate of $\sim1 - 4 \times10^{10}$ g s$^{-1}$. Since XUV radiation from an inactive M4 star such as GJ 876 cannot drive the atmospheric outflow or ionize Helium efficiently enough to produce metastable helium atoms, we hypothesize that the detected outflow is driven by XUV radiation from chromospheric active regions on the TOI-3884 host star, analogous to those in YZ Cmi. The combination of a detected pre-transit headwind and the lack of a post-transit tail is unusual. Its behavior is similar to that of HAT-P-67b, a hot Saturn orbiting an F5 dwarf star. TOI-3884b is orbiting a much cooler M4 dwarf. Such an extended headwind, which is longer than $\sim23$ times the planetary diameter, will need to be explained by advanced exospheric models of star-planet interactions. These discoveries also caution against using spectra taken close to the transit as out-of-transit baselines, since the exosphere can be many times larger than the exoplanet.

J.P.N. & S.M. thank Harish Vedantham and Sagnik Saha for discussion on the sub-Alfvénic interaction at ICTS-TIFR Strange New Worlds School in 2026 April.

The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. Computations for this research were performed on the Pennsylvania State University’s Institute for CyberScience Advanced CyberInfrastructure (ICS-ACI).

These results are based on observations obtained with the Habitable-zone Planet Finder Spectrograph on the Hobby-Eberly Telescope. We thank the Resident astronomers and Telescope Operators at the HET for the skillful execution of our observations of our observations with HPF. The Hobby-Eberly Telescope is a joint project of the University of Texas at Austin, the Pennsylvania State University, Ludwig-Maximilians-Universität München, and Georg-August Universität Gottingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly. The HET collaboration acknowledges the support and resources from the Texas Advanced Computing Center.

S.P.H acknowledges that this research was carried out, in part, at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

This research has made use of NASA’s Astrophysics Data System Bibliographic Services. Perplexity AI was used for checking grammar and getting critical feedback on language quality. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the publication.

J.P.N. acknowledge the support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI 4002. J.P.N acknowledges support in part by a generous donation (from the Murty Trust) aimed at enabling advances in astrophysics through the use of machine learning. The Murty Trust, an initiative of the Murty Foundation, is a not-for-profit organization dedicated to the preservation and celebration of culture, science, and knowledge of systems born out of India. The Murty Trust is headed by Mrs. Sudha Murty and Mr. Rohan Murty. O.K. acknowledges support by the Swedish Research Council (grant agreement no. 2023-03667) and by the Swedish National Space Agency. S.M., M.D., P.R. acknowledge support from NASA through the Exoplanet (XRP) research program grant 80NSSC24K0155. S.M., P.R., G.S. and J.L-R. are part of NASA’s CHAMPs team, supported by NASA under grant No. 80NSSC23K1399 issued through the Interdisciplinary Consortia for Astrobiology Research (ICAR) program.

Facilities: HPF (HET)

Software: HxRGProc ([32]), barycorrpy ([21]), astropy ([7, 8, 9]), numpy ([19]), scipy ([54]), matplotlib ([20]), GNU parallel ([51])