Metastable Helium Absorptions with 3D Hydrodynamics and Self-Consistent Photochemistry I: WASP-69b, Dimensionality, XUV Flux Level, Spectral Types, and Flares

The metastable Helium (He*) lines near $10830~{\rm A}$ are ideal probes of atmospheric erosion--a common phenomenon of close-in exoplanet evolution. A handful of exoplanet observations yielded well-resolved He* absorption features in transits, yet they were mostly analyzed with 1D isothermal models prescribing mass-loss rates. This work devises 3D hydrodynamics co-evolved with ray-tracing radiative transfer and non-equilibrium thermochemsitry. Starting from the observed stellar/planetary properties with reasonable assumptions about the host's high energy irradiation, we predict from first principle the mass loss rate, the temperature and ionization profiles, and 3D outflow kinematics. Our simulations well reproduce the observed He* line profiles and light curves of WASP-69b. We further investigate the dependence of He* observables on simulation conditions and host radiation. The key findings are: (1) Simulations reveal a photoevaporative outflow ($\sim 0.5~M_{\oplus}~{\rm Gyr}^{-1}$) for WASP-69b without a prominent comet-like tail, consistent with the symmetric transit shape (Vissapragada et al. 2020). (2) 3D simulations are mandatory for hydrodynamic features, including Coriolis force, advection, and kinematic line broadening. (3) EUV ($>13.6~{\rm eV}$) photons dominate photoevaporative outflows and populate He* via recombination; FUV is also detrimental by destroying He*; X-ray plays a secondary role. (4) K stars hit the sweet spot of EUV/FUV balance for He* line observation, while G and M stars are also worthy targets. (5) Stellar flars create characteristic responses in the He* line profiles.


INTRODUCTION
One of the most exciting discoveries in exoplanetary sciences in recent years is that the radii of sub-Neptune planets have a bimodal distribution (Fulton et al. 2017). The prevailing explanation is the atmospheric erosion by either photoevaporation (e.g. Owen & Wu 2013 or core-powered mass loss (Ginzburg et al. 2018;Gupta & Schlichting 2019. In any case, the prominence of the radius gap implies that atmospheric erosion is probably a stage of evolution that close-in exoplanets very commonly go through. Since the discovery of the Lyman α (Lyα hereafter) transit of hot Jupiter HD209458 b (Vidal-Madjar et al. 2003), Lyα has been a workhorse for studying atmospheric erosion (e.g. Lecavelier Des Etangs et al. 2010;Kulow et al. 2014). However, Lyα has some unavoidable limitations, namely, it is heavily contaminated by geocoronal emission, and interstellar absorption saturating the very center of the line (e.g. Ehrenreich et al. 2015). Moreover, one has to go space to observe this UV transition. These effects significantly limit the number of systems for which we can study Lyα transits.
Besides Lyα , helium lines are emerging as promising outflow indicators. The 2 3 S state of helium is often called the "metastable state" (He * hereafter), because the 2 3 S→ 1 1 S transition is a magnetic dipole process with a slow spontaneous decay rate of A 1.3×10 −3 s −1 (Drake 1971(Drake , 2006. Meanwhile, the transition between the lower 2 3 S and the upper 2 3 P J states of helium consists of three lines with A > 10 7 s −1 , whose wavelengths in vacuum are 10832.08Å (for the J = 0 upper state), 10833.24Å (J = 1), and 10833.33Å (J = 2) respectively. These lines are often referred to as "He i 10830Å lines" or the "metastable helium lines", as they are radiatively decoupled from the ground state. The abundance of helium, the absence of geocoronal or intersteallar containmination, the longevity of metastable state, and observability from the ground together enabled the He * lines as an excellent probe of ionized flows in various scenarios of astrophysics, including quasars (e.g. Leighly et al. 2011), stellar atmospheres and outflows (see Edwards et al. 2003, and references to the article), and T Tauri stars (Kwan et al. 2007).
Over the years, researchers have proposed the He * lines as a tracer of mass loss of close-in exoplanets (Seager & Sasselov 2000;Turner et al. 2016;Oklopčić & Hirata 2018). It was the secure detection of Spake et al. (2018) that revitalized interest in this unique transition. At the time of writing this paper, several close-in exoplanets have transmission He * line profiles resolved by ground-based spectrographs (e.g. Allart et al. 2018;Nortmann et al. 2018;Salz et al. 2018;Kirk et al. 2020;Ninan et al. 2020). More recently, Vissapragada et al. (2020) custom made a narrow band filter specifically for the He * transitions on the diffuser-based photometric system on Palomar/WIRC. The resultant precise light curves of the He * is complementary to the line profiles from the near-infradred spectrographs. A lot of information about atmospheric outflow is hiding in these line profiles and light curves waiting to be interpreted. The models that are commonly used in the literature to interpret these He * observations are 1D spherical symmetric models that are isothermal (Oklopčić & Hirata 2018;Oklopčić 2019;Palle et al. 2020), or have prescribed heating efficiency (Lampón et al. 2020). The model is widely recognized for its simplicity and effectiveness, however it has to prescribe, rather than predict, a mass loss rate or a temperature profile.
In this work, we build upon our previous model that conducts hydrodynamics, self-consistent thermochemistry, and ray-tracing radiative transfer to study the photoevaporation of sub-Neptune planets (Wang & Dai 2018, WD18 hereafter). We have streamlined the code so that it is sufficiently fast to run in 3D to fully capture outflow dynamics, and added various processes that are relevant to the (de)population of He * . We will show in this paper that using the observed stellar/planetary properties and making reasonable estimate of the high energy spectral energy distribution (SED) about the host star, our model can predict mass loss rate, the temperature profile, the ionization states, and synthesize the observed He * line observables from first principles.
In this first work of a series, we focus on WASP-69b, which is one of the first detections of He * line absorption with well-resolved line profile (Nortmann et al. 2018). Acknowledging the many limitations of a 1D isothermal model, Nortmann et al. (2018) did not tie their He * line observation with a theoretical model. Instead they only reported what the data showed directly. Notably, Nortmann et al. (2018) reported an asymmetric transit profile characterized by a longer-than-expected egress that could be caused by a comet-like tail associated with the mass loss. However, Vissapragada et al. (2020) suggests a symmetric shape of transit using better-sampled light curves with higher precision and signal-to-noise ratio (SNR). Another interesting point about WASP-69b is the apparent temporal variability of the He * transit depth seen in Nortmann et al. (2018). We will try to understand these observations of WASP-69b in the framework of our 3D hydrodynamic simulations. Afterwards, we will use WASP-69 as a fiducial case to investigate the impact of dimensionality, XUV flux level, and host spectral types on the observables of the He * lines.
This paper is structured as follows. §2 describes our methods of numerical simulations and synthetic observations. In §3, we present a fiducial model of WASP-69b that well reproduces all current observations. Based on this model, §4 studies how various system parameters impact the rate of photoevaporation and He * observables. §5 explores the possibility that stellar flare may cause the observed temporal variability of WASP-69b. §6 summarizes the findings of this paper.

Basic Setup
We conceptually divided a planet into four regions: (1) a dense core, (2) a convective inner atmosphere, (3) a quasi-isothermal outer atmosphere with equilibrium temperature T eq and (4) an outflowing region irradiated by high energy photons (e.g. Rafikov 2006;Ginzburg et al. 2016;Owen & Wu 2016). The equilibrium temperature T eq satisfies, where L * is the bolometric luminosity of the star, and a is the semi-major axis of the planetary orbit. Our simulations will focus on regions (3) and (4), whereas the structure of regions (1) and (2) provide the correct boundary conditions crucial for correctly reproducing the measured mass and radius of the planet. We refer the reader to Appendix A for the details of how we set up the internal structures of our planet and resultant boundary conditions for our simulations. We characterize the high energy radiation spectral energy distribution (SED) of the host star with 5 different characteristic energy bins: (1) hν = 2 eV for infrared, optical and near ultraviolet (NUV) photons, (2) hν = 7 eV for "soft" far ultraviolet (FUV) photons that can photoionize He * , (3) hν = 12 eV for the Lyman-Werner band FUV photons that can photodissociate molecular hydrogen, (4) hν = 20 eV for "soft" extreme ultraviolet (soft EUV) photons that can ionize hydrogen but not helium, (5) hν = 40 eV for hard EUV photons that ionize hydrogen and helium, and (6) hν = 3 keV for the X-ray.
Our simulation combines ray-tracing radiative transfer, real-time non-equilibrium thermochemistry, and full hydrodynamics calculations (based on a higher order Godunov hydrodynamic solver Athena++; Stone et al. 2020). The simulation is mostly based on our WD18 work with a few modifications and improvements added for the higher dimensionality and the inclusion of He * .

Geometry and Boundary Conditions
The density distribution, temperature profile and the dynamics of outflowing atmosphere all play a part in the He * observables. To capture the outflow dynamics accurately, simulations should include the gravity of the star and the planet and the effects of orbital motion: the centrifugal and Coriolis forces. Therefore, 3D simulations are required. Given its short orbital period and observed radial velocities (Anderson et al. 2014), we assume that WASP-69b is tidally locked and circularized. Our simulation is run in a co-rotating frame centered on the planet. We adopt a spherical polar coordinate (r, θ, φ) with θ = 0 pointing towards the host star and φ = 0 pointing in the direction of orbital motion.
The mesh covers the domain (r, θ, φ) ∈ [r in , r out ] ⊗ [0, π] ⊗ [0, π]. Planet-specific radial boundaries r in and r out usually extend from the the base of the quasiisothermal layer to a relatively large radii (150 R ⊕ in this case) such that the density/opacity drops low enough. The radial grids are placed logarithmically to strategically capture the steep change of density, while latitudinal and azimuthal grids are spaced evenly. Reflecting boundary conditions are enforced at the r = r in , φ = 0 and φ = π boundaries, while the r = r out boundary is an outflowing boundary. The θ = 0, θ = π boundaries are polar wedges to avoid coordinate singularity. The whole mesh, with its polar axis always pointing towards the host star, co-rotates with the orbital motion and the rotation of the tidally-locked planet.
In a 3D spherical polar mesh, the grids near the polar axis are narrow in the azimuthal direction (δx φ r cc sin θ cc δφ; subscripts "cc" stand for "cell center"). This can result in highly non-unitary aspect ratios (δx θ r cc δθ), and a stringent Courant-Friedrichs-Lewy (CFL) condition. We thus introduce an adaptive "mesh coarsening" technique for the azimuthal grids near the poles. Without any violations of conservation laws, the effective aspect ratio of the high-latitude zones become close to one and the timestep constraints imposed by CFL condition is not as severe (see also Nakamura et al. 2019;Müller et al. 2019). This helps to greatly speed up our model.

He * in non-LTE Thermochemistry
Our simulation includes a non-local-thermodynamicequilibrium reaction network that coevolves with the hydrodynamics (see WD18 for detail). With the addition of the metastable state of neutral helium and all relevant reactions that populate and de-populate this state (see Oklopčić & Hirata 2018 and references therein). Our reaction network now has 26 thermochemical "species" (24 chemical species in WD18, He * , and internal energy density) and 135 reactions such as ionization, recombination, collisional (de-)excitation, photodissociation, and cooling and heating processes. The ordinary differential equations (ODEs) of the thermochemical network were solved efficiently using the semi-implicit method specially optimized for the graphics processing units (GPUs). The resultant efficiency allows us to coevolve the hydrodynamics with the thermodynamics, rather than including thermodynamics as a post-processing step that is often done in the literature. Again, we refer interested readers to WD18 for more details.

Synthetic Observations
We synthesize both the transmission line profiles (Nortmann et al. 2018) and the narrow-band light curves (Vissapragada et al. 2020) of He * transitions using our simulations. At each wavelength λ and a particular orbital phase, the optical depth along a line of sight (LoS hereafter) is given by, (2) where we have transformed from our planet-centered coordinate systems in the simulations to a star-centered coordinate system for the synthetic observations. Thus x and v are the position and velocity vector measured from the host star. The integration goes along the designated LoS, and the summation index i runs over He * 's three lines with different upper state quantum number J. The cross-section σ i is assumed to be a Voigt profile which convolves the intrinsic Lorentz profile (γ = A/4π, A = 1.0216 × 10 7 s −1 ; see Drake 2006) with a Gaussian profile from thermal broadening at temperature T (x). This Voigt profile is shifted by the local bulk velocity including orbital motion and outflow kinematics and the projection onto the LoSn LoS · v(x).
This integration is performed for all relevant LoS that originate from the surface of the host star: where (λ) is the relative extinction at wavelength λ, S(λ) is the normalized surface brightness ( dΣS(λ) = 1) of the star after accounting for limb darkening and stellar rotation. The integral runs through the entire projected stellar surface. (λ; Φ) is the absorption line profile as a function of wavelength and orbital phase (time). We mimic what observers often do in He * observations i.e. time averaging (λ; Φ) over the entire transit event from nominal ingress to egress (t ii through t iii ) 1 . The outcome ∆ (λ) is a line profile of excess absorption to be compared with observations directly.
We report a number of summary statistics including the equivalent widths W λ ≡ ∆ (λ) dλ, the radial velocity shift of the absorption peak ∆v peak and the fullwidth-half-maximum (FWHM) of the absorption line profile. These summary statistics help us to compare between models and observations efficiently and are less prone to measurement uncertainty, bad pixels and other instrumental effects.
Finally, we integrated (λ; Φ) multiplied by a filter bandpass function over λ. The result is a transit light curve near the He * transitions. In this work, we use the bandpass function provided by Vissapragada et al. (2020) for a direct comparison with their results.

FIDUCIAL MODEL OF WASP-69b
In this section, we will show how we arrived at a fiducial model that gives rather remarkable agreement with the observed He * line profiles (Nortmann et al. 2018) and light curves (Vissapragada et al. 2020). We note that our 3D hydrodynamic model is not fast enough 2 for a full exploration of the parameter space with techniques such Markov Chain Monte Carlo or even simple gradient descent. Instead we had to rely on the reported system parameters and making reasonable assumption as well as hand tuning the high energy SED of the host star. We will see shortly, without much tuning, we can arrive at a fiducial model that fits various observations of WASP-69b very well.  We set up our simulations to match the reported system properties of WASP-69b (Anderson et al. 2014). The host star is K star with M * = 0.826 M , R * = 0.813 R and T eff = 4715 K. WASP-69b has a circular orbit with semi-major axis a = 0.04525 AU. The equilibrium temperature is estimated using Eqn. 1 T eq = 965 K. The planet has an optical transit radius of R p 1.057 R Jup and a mass of M p 0.26 M Jup from radial velocity follow-ups. Details of the fiducial model are presented in Table 1.
The interior of our WASP-69b model is set up as described in Appendix A. The core size, the equation of state and other details of the interior of a giant planet is still subject to a lot of uncertainties even in the case of Jupiter (see e.g. Wahl et al. 2017). However, the details of the interior should not affect the outflowing region of the envelope which is what we are interested in this work. We set the inner boundary of our simulation at 11.37 R ⊕ so that we capture several scale heights be-low the optical transit radius r eff R p at 1.057 R Jup . The outer boundary is located at 150 R ⊕ . For simplicity, we assumed an atmospheric metallicity often seen in protoplanetary disk (WD18) which is slightly below solar value (Table 1). We will explore any metallicity dependence in a future work.
Optical and infrared fluxes, represented by the hν = 2 eV photon energy bin, are calculated using the host star radius and effective temperature. For the high energy photons, more uncertainties are involved depending on the age/activity of the host star; while direct measurements are also lacking. We note that WASP-69 is moderately active indicated by the Ca ii H and K lines log R H,K = −4.54 (Anderson et al. 2014). As we will see later in §4.2, the He * absorption line profile depends critically on the shares taken by various high-energy radiation bins. After gaining intuition on how each energy bin affects the He * line profiles (again §4.2), we varied the high energy SED of WASP-69 until we achieved a reasonable agreement with both the line profile and light curve measurements. The resultant high energy SED is quite typical of a K5 star when compared to observational constraint of FUV and EUV flux of the MUSCLES survey Youngblood et al. 2016;Loyd et al. 2016;Youngblood et al. 2017), the X-ray flux according to Gudel (1992), and the more comprehensive compilation of Oklopčić (2019).

A Photoevaporative Outflow on WASP-69b
Before analyzing our simulations, we ensured that quasi-steady states have been achieved. This usually involved running the simulations for many dynamical timescales, specifically we set t sim 30 τ dyn . The dynamical timescale τ dyn is estimated by the sound crossing time of the Bondi radius: Here, c s is the sound speed. For a typical T = 10 4 K outflow we see for WASP-69b, τ dyn ∼ 10 4 s. Moreover, we also check explicitly if the simulation has settled down to a quasi-steady state by comparing key hydrodynamic/thermodynamic properties in neighboring dump files.
Our fiducial model for WASP-69b shows clear signs of a photoevaporative outflow. In Figure 1, we show 2D slices of the density, temperature and LOS velocity distributions of our 3D simulations looking down the North pole of the planet. We see a T 10 4 K hot ionized supersonic outflow originating at a wind base of r 13 R ⊕ which eventually accelerates to ∼ 23 km s −1 when leaving the domain of our simulation. This outflow disperses the planet atmosphere at a mass-loss rate oḟ M 5.5×10 −10 M ⊕ yr −1 . Since we assumed a constant high energy radiation output from the host star, the mass-loss rates and the hydrodynamic/thermodynamic profiles remain nearly constant after reaching the quasisteady state. We also note that we did not put in any stellar wind from the host star, as the current windless model fits the data reasonably well and is preferred by Occam's Razor. However in a companion paper on WASP-107 we will show that stellar winds may generate Kelvin-Helmholtz instability that leads to fluctuations of a photoevporative outflow.
Which mechanisms control the population of the He * state? The bottom row in Figure 2 compares the rates of different (de-)population processes along the two particular streamlines (thickened curves in Figure 1). We compute the rate of ionization, recombination, spontaneous decay, collisional excitation and de-excitation as well as an advection attenuation term |v · ∇n(He * )|. Along the representative streamline presented by the left column, the abundance of He * is determined by the relatively stiff balance between the recombination (He + + e − → He * ) and the collisional de-excitation at small radii 30 R ⊕ . As expected these two processes are efficient at higher densities consistent with the law of mass action. Photoionization of He * by soft FUV starts to take over the destruction channel of He * where the density of free electrons declines at higher altitudes. The other channels have negligible importance: e.g. collisional excitation from 1 1 S to the metastable state is more than five orders of magnitude slower than recombination. On the right column of Figure 2, we show an interesting streamline that crosses in the shadow of the planet. The number density of He * soars in the shadow because the photoionization of He * by soft FUV vanishes here.
Beneath the base of the photoevaporative outflow (r 13 R ⊕ ), the temperature gradient between the day-side and the night-side generates a slow "zonal" circulation (∼ 0.1 km s −1 ). However, this region has little observational effect on the overall He * * observables which are mostly controlled by the much more extended low density regions of the outflowing atmosphere. We will return to this point shortly. Moving to higher altitude, this day-night advection continues, amounting to a 2 − 3km s −1 blueshift at about 20 − 40 R ⊕ . Going further away from the planet, the Coriolis effect starts to shape the outflow streamlines into spiral curves resulting in blueshifts on the leading edge and redshifts on the trailing edge. Considering that the outflow is still primarily radial, increments in the latitudinal velocity |∆v θ | after traveling through a radial distance ∆r can v K Ω K Figure 1. Profiles of the simulation for WASP-69b (fiducial model 69-0) in its quasi-steady-state. Stellar radiation comes from the left of the plot, and the orbital angular velocity ΩK points out of the paper plane; the Keplerian motion of the plaent is upwards vK. Colormaps show the mass density ρ (upper left panel), temperature T (upper middle), line-of-sight velocity vLoS (upper right; the value is measured at mid-transit), He * number density n(He * ) (lower left), inverse timescale of recombination He * formation (defined as formation rate normalized by n(He * ); lower middle), and free electron number density n(e − ) (lower right). White streamlines (projected to the orbital plane) are overlaid on each panel; two neighbor streamlines are separated in such a way that they are ∆θ = π/16 apart when they reach the outer radial boundary (r = 150 R⊕). The heavy streamline are the reference lines on which the profiles are plotted in Figure 2). Black solid lines indicate the sonic surface.
be estimated by, This estimation is confirmed by the velocity profile in the top panels of Figure 2: if we compare the values at r 40 R ⊕ and r 100 R ⊕ , the difference in v LoS (approximately equals to v θ for this streamline) is ∆v LoS ∼ 11 km s −1 , and eq. (5) yields ∼ 10.8 km s −1 . This effect re-distributes He * atoms in the velocity space and broadens the observed He * line profiles as we will see shortly.  Figure 1). The top abscissa is the radial coordinate r corresponding to the curve length along the streamline on the bottom abscissa. The top panels contain the profiles of scaled mass density ρr 2 , temperature T , velocity magnitude |v| and the line-of-sight velocity vLoS. Dashed part of the vLoS curve indicate negative values. The middle panels present the abundances of free electrons (relative to total hydrogen nucleus density nH) and helium in different forms (relative to the total helium nucleus density nHe; note that n(He * ) is multiplied by 10 5 for clarity). Inverse timescales of He * formation (in solid curves; note that the collisional excitation rate is multiplied by 10 5 ) and destruction (in dashed curves) processes are shown in the bottom panels. In all panels, the vertical dotted line indicates the sonic critical point, and the vertical dash-dotted line shows the location of Roche radius. Figure 3 shows our synthetic observations of both the line profiles and light curves of WASP-69b (Nortmann et al. 2018;Vissapragada et al. 2020). We have binned the light curve data from Vissapragada et al. (2020) for better clarity and the uncertainty represents the standard deviation within each phase bin. Our fiducial model of WASP-69 seems to fit both the spectroscopic and photometric observations well simultaneously. In particular, the synthetic line profile reproduced the subtle blueshift of the peak absorption, the overall line depth, and the relative ratio between the lines of this triplet. Numerically, Nortmann et al. (2018) reported a net blueshift of −3.58 ± 0.23 km s −1 . This blueshift was based on fitting Gaussians to the observed line profiles; however, as we argued in the previous section, kinematic shift of the outflow introduces significant distortion of the spectral shape. Instead of fitting Gaussian to our line profile, we compared our simulations directly to the line profile itself which shows great agreement. We also use a different way to measure the blueshift: we report the blueshift of the peak absorption relative to a lineratio-weighted average of the rest-frame line center for the two longer wavelength transistions that are usually blended together. Nortmann et al. (2018) hinted at the possibility of a comet-like tail trailing behind WASP-69b. The basis of their suggestion is that additional He * absorption can  still be seen ∼ 20 min after the nominal egress of the planet. The higher precision, better temporally sampled photometric data from Vissapragada et al. (2020) nontheless favors a symmetric transit. The symmetric transit shape (lower panel of Figure 3) does not support an extended comet-like tail. Our simulation seems to be more consistent with Vissapragada et al. (2020), the photoevaporative outflow of WASP-69b in our fiducial model is largely symmetric between the leading and trailing edge, hence it produces a more symmetric transit shape. We note that an extended comet-like tail will also introduce significant distortion to He * line profile (see our Companion paper on WASP-107b for strong comet-like tail generated by strong stellar wind in that system). Here in the case of WASP-69b, our fiducial model produces a good fit the resolved line profile (Nortmann et al. 2018) while it does not need to invoke a prominent comet-like tail. Another point we would like to emphasize is that a significant part of the He * absorption for WASP-69 seems to be produced by an extended, optically thin (τ < 1) outer layer of the photoevaporative outflow. In Figure 4, we show the mid-transit extinction (1 − exp −τ ) at three different wavelengths near the He * transitions. The outer regions (10s of R ⊕ ) contribute significantly to the overall extinction thanks to their extended area and the slow decrease of He * number density in the outflow. Because of the unsaturated optical depth, the line ratios between the He * triplet are close to 1 : 3 : 5 i.e. their quantum degeneracies. More accurately, the line ratios are close to 1 : 8 as the longer two lines are blended by kinematics and thermal broadening. This suggests that the line ratios between the He * triplet can be a diagnostic of the number density in the outflow. If most of the He * absorption is due to higher-density region where one of the line may saturate first, the line ratio will deviate from the quantum degeneracy ratio; this would tell us about the density of the outflow region in a model-independent way (see also discussions in Salz et al. 2018). This does not seem to be the case for WASP-69b, as most He * absorption happens in lower density regions.

Comparison with Observations
It is worth noting that, due to heavy computational cost of our 3D simulations, we could not afford to numerically fit the data with multiple simulation runs. Instead, our fiducial model serves as a validation our self-consistent 3D hydrodynamic simulations: with reasonable assumptions of the planetary/stellar properties and high energy SED, we can at least qualitatively reproduce the various He * observables. That said, the degree to which our simulation agrees with observations is quite encouraging if not remarkable. After this validation of model, we stand at a position to perturb our fiducial model and investigate how the He * observables are impacted by various factors that control the underlying photoevaporative outflow and He * population in the following section.

PARAMETRIC STUDY
How do the photoevaporative outflow and the resultant He * observables depend on key parameters in our simulations? In this section, we explore the impact of simulation dimensionality, XUV flux levels, host star spectral type, and planet surface gravity. This is done by perturbing validated fiducial model in these parameters. Before that, we did a further convergence test. We reran the fiducial model with a much finer simulation grid of N log r × N θ × N φ = 192 × 192 × 128 (versus the fiducial model, N log r × N θ × N φ = 144 × 128 × 64). The resultant He * observables in the convergence test are almost identical (Figure 3) to the much faster fiducial model. This gives us confidence that the adopted spatial grid is fine enough to resolve the photoevaporative outflow on WASP-69b.

Dimensionality
To test how our model depends on the spatial dimensions of the simulations, we ran a 2D model with axisymmetry (about the φ axis i.e. N φ = 1) while keep all system parameters the same as the fiducial model. The Coriolis forces is not captured in this 2D simulation while stellar gravity and orbital centrifugal forces are still involved. A reference 1D spherical symmetric model is also implemented using the θ = π/2, φ = π/2 radial line and removing the θ and φ components of the velocity. Figure 3 compares the synthesized line profiles and light curves with all three dimensionality models. The 1D spherically symmetric model suffers from the loss of all non-radial kinematic information. It is clearly inconsistent with the observed line profile with no blueshift and limited kinematic broadening. The 2D axisymmetric test is able to capture the day-night advection. It shows good agreement with the observed line profile. The 3D fiducial model further modifies the line profile by including the Coriolis force. In this case of WASP-69b such a modification is quite subtle, which again tes-tifies that the photoevaporative outflow on WASP-69b is largely symmetric between the leading and trailing edge. Again see our companion paper on WASP-107b for how this symmetry is broken by the inclusion of stellar winds.
The three models from 1D through 3D have almost identical equivalent width ( W λ ∼ 3.1Å) and light curves. This stresses the importance of spectrally resolving the He * line profiles which are seen to vary the most between dimensions. The mass-loss rate are again quite similar between 3D and 2D models at abouṫ M 5.5×10 −10 M ⊕ yr −1 . The mass-loss rate in our 1D model is off (Ṁ 6.9 × 10 −10 M ⊕ yr −1 ) because it assumes perfect spherical symmetry. However, the streamlines in Figure 1 are clearly non-radial. We also compare our results with that from a 1D isothermal model (Oklopčić & Hirata 2018;Vissapragada et al. 2020) of 9.5 × 10 −10 M ⊕ yr −1 ( 3 × 10 −3 M Jup Gyr −1 ) at an assumed temperature of 12000K. The results are in rough agreement with differences arising from more careful treatment of the hydrodynamics, thermodynamics and radiative transfer.

XUV Flux Intensity
Photoevaporative outflows are driven by high energy radiation from the host star. Moreover, the population of He * states are also controlled by the critical balance high energy photons of different energy bins. We examine the impact of high energy radiation in each energy bin by perturbing the fidual model. The amount of high energy radiation a star outputs is variable depending on the evolution stage, activity and spectral types of the host star. Direct measurements are also lacking as the XUV measurements have to be performed in space. Note-For simplicity, F N (hν) ≡ F (hν)/(10 N cm −2 s −1 ), calibrated for the value at the planet orbit without any extinction. Vertical dotted lines mark the boundary between different energy bands: "Hard EUV" for hν > 24.6 eV photons that can ionize helium; "Soft EUV" for 13.6 < (hν/eV) < 24.6 photons that can ionize hydrogen ; "LW" (short for Lyman-Werner) for 11.2 < (hν/eV) < 13.6 photons that can photodissociate H2; "Soft FUV" for hν < 7 eV photons.
Therefore, in our Models 69-1 to 69-5 (soft FUV for 69-1, LW for 69-2, soft EUV for 69-3, hard EUV for 69-4, X-ray for 69-5), we bump up the flux level in each high energy bin by a whole order of magnitude to reflect the intrinsic variation in high energy flux level. We summarize key He * observables in Table 3; we also show the synthetic line profiles and light curves in Figure 6 and the relative abundance of He * as a function of radius in Figure 8. The He * line profiles are controlled mainly by the FUV (adverse effect) and the EUV bands (positive effect); the LW and X-ray bands only play minor roles under  Table 3. The fiducial model 69-0 is also included for reference.
the "typical" host star conditions. More specifically, the relative abundance of He * is suppressed by soft FUV because it photoionizes the He * states, thus significantly reducing its population and the He * absorption depth. In Model 69-1 (×10 soft FUV flux), the over-all massloss rate is enhanced by a few percent thanks to extra energy deposited into the atmosphere. However the stronger soft FUV flux significantly lowered the number density of He * at a larger radial extent (r 20 R ⊕ ). The He * line profile depth and the light curve depth are both reduced by about a factor of two (3.16Å to 1.75Å). The line width also decreased (FWHM from 17.3Å down to 14.0Å) because the high-altitude region with a higher velocity dispersion contribute less to the He * extinction.
In Model 69-2 (×10 LW flux) the He * observables are mostly unaltered from the fiducial model. This is because the LW band is intrinsically narrow thus only amount to a very small fraction of the overall high energy radiation flux. Moreover, most molecular H 2 are already dissociated at the ∼ 10 4 K in our simulations.   Table 3. Note fiducial model 69-0 is a K star. Figure 8. Relative abundance profiles of He * (relateive to total helium nucleus density nHe), measured along the radial line with θ = π/4, φ = 0 in each simulation domain, for models described in Table 3. The upper panel shows models 69-0 and 69-1 through 69-6, while the lower panel specifically compares the results of different types of host stars (Model 69-0 for the K star WASP-69, and Models 69-F, 69-G and 69-M for F, G, M stars respectively).
At higher EUV fluxes (Model 69-3 and 69-4), the much faster outflows not only bring more He * into the exosphere but also spread them out in velocity space effectively broadening the line profiles. This confirms our earlier picture that EUV flux are most effective in driving photoevaporative outflows (WD18).
Finally, X-ray seems to play a secondary role in photoevaporation and He * observables. Model 69-5 (×10 X-ray flux) has very similar observables as the fiducial model. Although X-rays photons are very energetic, they also have much smaller cross sections than EUV photons. As a result, X-ray penetrate deeper into the atmosphere where collisional processes and dust particles quickly convert the X-ray energies to infrared radiation that escapes easily. This limits the heating potential of X-ray. We reached a similar conclusion in WD18.

Host Spectral Type
Among the handful of reported He * detections, 4 out of 6 are planets around K-type hosts (Spake et al. 2018;Venzmer & Bothmer 2018;Allart et al. 2018;Salz et al. 2018;Ninan et al. 2020;Alonso-Floriano et al. 2019). Oklopčić (2019) confirmed that K-stars, at least in 1D isothermal models, may be at the sweet spot of FUV/EUV flux balance that best promote the He * population and thus observablity. This section re-evaluates such a claim with our 3D hydrodynamic simulation coupled with self-consistent thermodynamics and radiative transfer.
We set up three additional models 69-F, 69-G, and 69-M, whose luminosities in different high energy bins emulate typical F-type, G-type, and M-type main-sequence stars based on the compilation of Oklopčić 2019. We remind the reader the fiducial model 69-0 has a K-star SED. The flux in each high energy bin is summarized in Table 2). Broadly speaking, F-type and G-type stars output similar levels of soft and hard EUV fluxes as K-type stars, however their FUV luminosities are significantly higher by a factor of ∼ 3000 and ∼ 130 respectively. For a typical M star, fluxes in all high energy bands are lower by about one order of magnitude.
He * observables of these models are summarized in Table 3 for their mass-loss rates and a few key diagnostics etc. We show the line profiles and light curves in Figure 6 and the radial distribution of He * in Figure 8. With a soft FUV flux ∼ 3000 times higher than our fiducial model, F-type stars significantly suppress the population and observability of He * with an equivalent width reduced by almost two orders of magnitude (3.16 to 0.04 A). However, the mass loss rate of the photoevaporative outflow is similar to that of the fiducial model. Again, photoevaporation is driven mostly by EUV which have similar flux levels between F and K stars. On the other hand, for a typical M star host, whose higher energy flux levels are weaker in all bands, the mass loss rate is reduced by one order of magnitude. Nonetheless, the equivalent width of He * in the transmission spectrum only decline by a factor of 2 (3.16 to 1.4Å). Again this is because its much weaker soft FUV flux allows proportionally more He * to exist in the outflow (Figure 8). G-star represents some middle ground, its ∼ 130 times higher soft FUV flux suppress the equivalent width He * by a factor of 6.
In summary, our findings suggest that K-star planet hosts are indeed favorable targets for He * observations consistent with the suggestion of Oklopčić (2019). The high-energy SED, nonetheless, is expected to change significantly as a function of host star age and activity level. The suppression factor of He * around G and M type stars are often only a factor of a few. We encourage observers to keep them in their target list particularly the young and active ones. We also predict that there will be more reports of He * detection around G and M type hosts soon. Another important point we would like to stress is that the depth of He * line profile can not be translated to the underlying mass loss rate without knowing the high energy SED of the host star. In other words, measuring the XUV SED of the host star directly is crucial for correctly interpreting the He * observations.

Surface gravity
The mass-loss rate of photoevaporation depends quite strong on the depth of gravitational potential well of the planet. A shallower potential allows faster outflow with the same high energy irradiation. In Model 69-6 we adjusted the planet interior such that the planet mass is reduced by half while keeping the transit radius the same. This effectively lowers the surface gravity of the planet by a factor of two. The mass loss rate in Model 69-6 increases toṀ 9.3×10 −10 M ⊕ yr −1 which is ∼ 70 % larger than the fiducial model. This larger mass loss rate can be decomposed into an increase in the terminal outflow velocity by ∼ 10 % and an increase of the outflow density by ∼ 60 %.
The He * line profile depth responds to this increase of mass loss rate sub-linearly. In Figure 6, the line profile has ∼ 40 % larger depth than that in the fiducial model, while the equivalent width increases by about ∼ 30 %. However, the He * line profile maintains a similar morphology with the peak velocity shift and the FWHM unchanged from the fiducial model (Table 3). In short, puffier, low surface gravity planets are more likely to undergo strong photoevaporative mass loss and should prove great target for He * observations.  Relative Soft FUV Flux Figure 10. Slimilar to the lower panel of Figure 3, showing the light curves following soft-FUV-only flares that start at mid-transit (∆t = 0) and terminating after 0.5 hr. Different flare intensities are indicated by colors. The light curves in solid lines should read the left ordinate, while the dashed lines are the flare shapes and should be read with the right ordinate. Soft-FUV-only destroys He * while does not significantly perturb the photoevaporative outflow, thus only produces a decrease in He * line depth that is seen in WASP-69b. Nontheless, we do not believe this is the explanation for the observed variability (see §5 for detail).
drop in magnitude that lasted for about 20 minutes. This variability could well be instrumental in origin, but here we explore an alternative explanation that it is generated by a stellar flare on the host star. Solar/stellar flares are associated with the surface magnetic activity of the host star. Their amplitudes can range below a percent to even orders of magnitude in extreme cases ("superflares" Günther et al. 2020, and references therein). The sudden rise of luminosity is often followed by exponential decays to nominal flux level on minutes or hours timescale. To investigate the consequences of flaring events on He * observables, we first inject a simple flare model in which fluxes across all energy bands energy bands increase by a factor 10 and 100 times which then quickly decay to the quiescent state exponentially with a timescale of 500 s.
The temporal response of our fiducial model to these flares are shown in Figure 9. Independent of the energy injected, the common mode of response is as follows. Before the dynamics of the outflowing atmosphere can fully respond to the flares, the photoionization of He * by soft FUV photons reduces the number density of He * and cause a significant decrease in the equivalent width. It is only after the dynamical timescale t − t flare ∼ τ dyn ∼ hours that the flare-generated surge of photoevaporative mass-loss reach the higher altitudes where most of He * absorption occurs. As a result, the equivalent width of He * increases after the first hour or so and remains high for several hours. Looking at the mass-loss rate across the outer boundary of the simulation (lower panel of Figure 9), it experiences several oscillations on dynamical timescales as the systems response to the increased flux from the flare. The equivalent width (and other observables) however are the spatially integrated quantities, thus it effectively smears out most of these oscillations and has a much smoother variation (upper panel of Figure 9). Comparing with variability seen in WASP-69b (Nortmann et al. 2018), a flare that simultaneously raises all high energy radiation does not appear to be a good explanation. This is because it should be observed as a decrease followed by an increase of He * absorption rather than the decrease only in the observations (Nortmann et al. 2018).
We therefore explore a different flare model in which only the soft FUV band. We do not have observational support that flares of this kind exist. We explore this rather contrived scenario just out of curiosity. Remember from §4.2 that the soft FUV primarily suppresses the He * abundance by photoionization without significantly changing the overall kinematics. A soft-FUV-only flare may produce the observed decrease of He * line depths. We setup three extra simulation runs, again based on the fiducial model of WASP-69b. We put in soft FUV flares (10, 100 and 1000 times the nominal value) that start at the middle of the transits and last for 30 min. In Figure 10, we can see that the light curves respond to these soft-FUV-only flares quickly. In order to reproduce the ∼30% variation, a soft FUV flare between 10 − 100× the nominal level is required. This should be readily observable in the Ca ii H, K lines of the CARMENES spectra in Nortmann et al. (2018). Nevertheless, enhanced activity was not observed in the spectra during or preceding the observed light curve variation (private communications, Nortmann). Afterall, stellar flares do not seem to be a viable explanation of the temporal variability seen the He * line profiles of Nortmann et al. (2018); instrumental effect is perhaps a better solution.

SUMMARY
In this work, we simulate the ionized mass loss of close-in exoplanets and the metastable helium absorption during the planetary transit. We produce synthetic spectrally resolved line profiles and the light curves in a narrow filter band around the He * transitions. Dynamics of such synthesis requires 3D hydrodynamic simulations of photoevaporating planetary atmospheres; nonequilibrium thermochemistry and ray-tracing radiative transfer are co-evolved with the hydrodynamics. The processes that populate and depopulate the metastable state of neutral helium are included in a thermochemical network and solved efficiently on GPUs.
With reasonable assumptions about the system parameters and high energy SED of WASP-69, we find a plausible model that launches a photoevaporative outflow with a mass-loss rate ofṀ 5.5 × 10 −10 M ⊕ yr −1 . The model yields a spectrum and a light curve that are in remarkable agreement with the observations in terms equivalent width, line-ratios, blueshift and line broadening (Nortmann et al. 2018;Vissapragada et al. 2020). Inside this outflow, metastable helium is formed almost solely by recombination. Its destruction is mainly due to collisional de-excitation at small radii where the density is high and photoionization by FUV photons at outer lower-density regions.
With this fiducial model of WASP-69b, we investigated how the photoevaporative outflow and He * observables depend on various input parameters. 3D simulations are crucial for capturing the full hydrodynamics including Coriolis force and advection. These effects are needed for producing correct line profile including the line ratios, kinematic broadening and the overall blueshift. We found that EUV photons are most efficient in driving the photoevaporation dynamics and in producing He + as the progenitors of recombination excitation of He * . The soft FUV photons that can ionize He * often play an adverse effect on the He * observability. X-ray photons, having much lower interaction cross section, are of secondary importance. Surface gravity also determines the effectiveness of photoevaporative outflows with puffier planets experiencing significantly stronger outflows, but the response of He* equivalent width is sub-linear.
K-stars are at a sweet spot of FUV/EUV balance that maximize the detectability of He * lines. F or earlier type stars have excessive FUV fluxes that suppresses the He * lines by orders of magnitude. G and M dwarfs represent a middle ground: He * lines should still be observable particularly for the younger and more active ones. In any case, the depth He * line profiles cannot be translated to a mass-loss rate without knowing the host star high energy SED.
We also investigated whether stellar flares could explain some of the temporal variability of WASP-69b (Nortmann et al. 2018). We found that a flare which enhances all high energy radiations initially suppresses He * lines due to FUV fluxes ionizing the He * before the whole system can adjust to higher mass-loss state after some dynamical timescales (usually hour-timescale). This characteristic shape is not consistent with the observed temporal variability of WASP-69b which only shows a decline of He * line depth before returing to nominal levels. Stellar flares are unlikely to be the explanation for this type of variability.