# Space Weather-driven Variations in Ly$\alpha$ Absorption Signatures of
Exoplanet Atmospheric Escape: MHD Simulations and the Case of AU Mic b

Ofer Cohen Lowell Center for Space Science and Technology, University of Massachusetts Lowell, 600 Suffolk Street, Lowell, MA 01854, USA  Julián D. Alvarado-Gómez Leibniz Institute for Astrophysics Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany  Jeremy J. Drake Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Laura M. Harbach Astrophysics Group, Department of Physics, Imperial College London, Prince Consort Rd, London, SW7 2AZ, UK  Cecilia Garraffo Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Federico Fraschetti Dept. of Planetary Sciences-Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ, 85721, USA  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

###### Abstract

We simulate the space environment around AU Microscopii b and the interaction
between the magnetized stellar wind with a planetary atmospheric outflow for
ambient stellar wind conditions and Coronal Mass Ejection (CME) conditions. We
also calculate synthetic Ly$\alpha$ absorption due to neutral hydrogen in the
ambient and the escaping planetary atmosphere affected by this interaction. We
find that the Ly$\alpha$ absorption is highly variable due to the highly-
varying stellar wind conditions. A strong Doppler blue-shift component is
observed in the Ly$\alpha$ profile, in contradiction to the actual escape
velocity observed in the simulations themselves. This result suggest that the
strong Doppler blue-shift is likely attributed to the stellar wind, not the
escaping neutral atmosphere, either through its advection of neutral planetary
gas, or through the creation of a fast neutral flow via charge exchange
between the stellar wind ions and the planetary neutrals. Indeed, our CME
simulations indicate a strong stripping of magnetospheric material from the
planet, including some of the neutral escaping atmosphere. Our simulations
show that the pressure around close-in exoplanets is not much lower, and may
be even higher, than the pressure at the top of the planetary atmosphere.
Thus, the neutral atmosphere is hydrodynamically escaping with a very small
velocity ($<15~{}km~{}s^{-1}$). Moreover, our simulations show that an MHD
treatment is essential in order to properly capture the coupled magnetized
stellar wind and the escaping atmosphere, despite of the atmosphere being
neutral. This coupling should be considered when interpreting
Ly$\alpha$observations in the context of exoplanets atmospheric escape.

stellar winds, outflows— planets and satellites: atmospheres —
magnetohydrodynamics (MHD) — planets and satellites: magnetic fields —
planet–star interactions

##  1 Introduction

The orbital period vs. planetary radius distribution of confirmed Neptune-size
exoplanets reveals a gap at periods less than about 4 days where no “hot
Neptunes” are found. This “Neptunian Desert” has been explained by a very high
mass loss rate that Neptune-size planets experience at these orbits, leading
to the complete evaporation of their gaseous envelope and leaving behind a
bare terrestrial-size core. Planets with a much larger size at these orbital
distances (hot Saturns and hot Jupiters) seem to have sufficient mass to
sustain their gaseous envelope despite their high mass loss rates (e.g., Szabó
& Kiss, 2011; Beaugé & Nesvorný, 2013; Batygin et al., 2016; Lundkvist et al.,
2016; Matsakos & Königl, 2016; Mazeh et al., 2016; Owen & Lai, 2018). The
evaporation of the gaseous envelopes in hot Neptunes occurs so quickly that it
is very hard to find observational evidence to support the proposed mechanism.
The recently discovered exoplanet AU Microscopii b (AU Mic b hereafter,
Plavchan et al., 2020; Martioli et al., 2021) offers a unique opportunity to
investigate a short-orbit, Neptune-size planet as it undergoes strong
evaporation on its way to losing most of its envelope (AU Mic b also co-exists
with a debris disk).

Evaporation of gaseous envelopes (the primary, mostly hydrogen atmosphere) of
short-orbit planets has been investigated using H I Ly$\alpha$ observations of
several hot Jupiters in absorption during transits. These observations
revealed conspicuous blue Doppler shifts in the line, indicating potentially
strong outflows of neutral hydrogen from the planet (e.g., Vidal-Madjar et
al., 2003; Lecavelier des Etangs et al., 2004; Linsky et al., 2010). The
evaporative outflows have been attributed to strong stellar ionizing radiation
that leads to hydrodynamic escape (“Parker Wind”, Parker, 1958) of the neutral
hydrogen (e.g., Erkaev et al., 2007; Murray-Clay et al., 2009; Owen & Jackson,
2012; Shaikhislamov et al., 2014). More recently, similar signatures have been
observed in the Helium 10830Å line (e.g., Spake et al., 2018; Oklopčić &
Hirata, 2018a, b). The importance of the latter is that this line can be
observed from the ground and it does not get absorbed by the ISM, in contrast
to the Ly$\alpha$ line which can only be observed from space due to geocoronal
absorption by neutral H in the Earth’s atmosphere.

The Ly$\alpha$ and He observations indicate very high mass loss rates of
$10^{8}-10^{10}~{}g~{}s^{-1}$. In some cases the blue Doppler shift in
Ly$\alpha$ line shows speeds of $100~{}km~{}s^{-1}$ or more. Even if the
hydrodynamic escape is very strong due to copious EUV irradiation, it still
cannot fully explain the high velocity of the escaping neutral hydrogen, and
it is even harder to explain in the case of the heavier helium atoms (e.g.,
Oklopčić & Hirata, 2018a).

The high Doppler shift velocity has recently been investigated by McCann et
al. (2019). Using a hydrodynamic simulation of the interaction between the
stellar wind and the escaping atmosphere, they showed that the stellar wind
can sweep up the escaping material and accelerate it to the observed
velocities. They also discussed the possibility that the high velocity can be
attributed to charge-exchange between the ionized, fast stellar wind ions and
the neutral, slow escaping hydrogen atoms, but they did not include this
process in their simulation. The interaction between the escaping atmosphere
and a generic stellar wind has also recently been studied by Shaikhislamov et
al. (2016); Carolan et al. (2020); Kubyshkina et al. (2022) using a
hydrodynamic model. They have found that a stronger stellar wind reduces the
mass loss rate by a factor of 2. Hazra et al. (2021) used a similar
hydrodynamic approach to study the impact on the escaping atmosphere of HD
189733 by stellar CME and flares. They found that the high blue-shifted
velocity may be attribute to CMEs.

Both the simulations of McCann et al. (2019) and Carolan et al. (2020)
described the stellar wind - planetary outflow interaction using hydrodynamic
models, where the magnetic field and electromagnetic forces have been
neglected. However, it is known from the interaction of the solar wind with
planetary atmospheres in the solar system that this interaction is controlled
by the magnetic field - plasma interaction, even in the case of non-magnetized
planets, where a magnetosphere can be induced by the stellar wind pushing
against the planet’s upper atmosphere (Kivelson & Russell, 1995; Schrijver &
Siscoe, 2009; Russell, 2009). Therefore, the ability of hydrodynamic
simulations to properly capture the wind - atmosphere interaction is limited.
Ekenbäck et al. (2010); Kislyakova et al. (2014, 2019) used direct Monte-carlo
simulations that included the effect of the magnetic field. Kislyakova et al.
(2019) used an MHD treatment to set the lower boundary conditions. These
detailed simulations, which included photoionization and charge exchange,
showed that the observed high velocity Doppler shift could be attributed to
Energetic Neutral Atoms (ENA) created via stellar wind - escaping neutrals
charge exchange (see also Holmström et al., 2008). The focus in these
simulations was to fit a particular set of parameters to the observed signal
at the time of the observation.

More recently, Harbach et al. (2021) performed MHD simulations of the stellar
wind interacting with an atmospheric outflow of a TRAPPIST 1e-like planet.
They showed how the Ly$\alpha$ line profile can change dramatically on
timescales as short as an hour as a function of the stellar wind conditions at
different orbital phases. The stellar wind conditions themselves can vary by a
few orders of magnitude along the orbits of short-period exoplanets (e.g.,
Garraffo et al., 2016, 2017).

In this paper, we extend the work by Harbach et al. (2021) to predict how the
Ly$\alpha$ observation varies as a function of the stellar wind in AU Mic b.
We simulate the stellar wind conditions using an MHD model for the stellar
corona and stellar wind, and use the results to drive an MHD model that
simulates the stellar wind - magnetosphere interaction. The results of the
latter are used to investigate the Ly$\alpha$ line profile as a function of
the planetary orbital phase. In addition to the effect of the ambient Stellar
Wind (SW), we also investigate the effect of the Coronal Mass Ejection (CME)
conditions, obtained from detailed CME simulations by Alvarado-Gómez et al.
(2022), on the planetary outflow and its associated Ly$\alpha$ signature.

We discuss the model setup and the input parameters in Section 2, and describe
the results in Section 3. In Section 4, we discuss the consequences of our
modeling work for future observations of atmospheric escape in AU Mic b, and
we conclude our findings in Section 5.

##  2 Model Description

###  2.1 Background Stellar Wind Model

We use the Alfvén Wave Solar Model (AWSoM, van der Holst et al., 2014) to
obtain the SW conditions near AU Mic b. The model solves the non-ideal
magnetohydrodynamic (MHD) equations, which include Alfvén waves coronal
heating and SW acceleration, as well as thermodynamic terms, on a spherical
grid between the base of the stellar corona and $200R_{\star}$. The assumed
stellar parameters in the simulations are: $R_{\star}=0.75R_{\odot}$,
$M_{\star}=0.5M_{\odot}$, and $P_{rot}=4.85~{}days$ (Kiraga, 2012; Plavchan et
al., 2020), and the smallest grid size near the inner boundary is
$0.025R_{\star}$ in the radial direction and $1.4~{}degree$ in the
latitudinal/azimuthal directions.

The inner boundary conditions (which essentially define the final solution)
are constrained by surface magnetic field data in the form of a “magnetogram”.
Magnetogram data for AU Mic has been made available from Kochukhov & Reiners
(2020) and from Klein et al. (2021), where the former also proposed the
presence of an additional dipole component of about $2~{}kG$. We investigated
test simulations for both Kochukhov & Reiners (2020) and Klein et al. (2021)
magnetograms but found the results to be similar, especially in terms of the
variation of conditions along the orbit of AU Mic b (see Alvarado-Gómez et
al., 2022, for a more complete description of the magnetograms and resulting
simulated wind conditions). Since the main focus of this work is the variation
of the absorption profile due to SW variations, we focus here on the SW
parameters obtained from the simulations driven by the magnetogram derived by
Klein et al. (2021).

Once the inner boundary conditions are specified, the initial condition for
the three-dimensional magnetic field is calculated using the potential field
approximation (Altschuler & Newkirk, 1969). The model then solves the non-
ideal MHD equations, where it accounts for the Alfvén waves coronal heating
and wind acceleration, as well as coronal thermodynamics, electron heat
conduction and radiative cooling, until a steady-state is obtained. The final
steady-state represents the state of the stellar corona and SW during the time
that the magnetogram data had been obtained. This approach to obtain the
stellar corona and SW solution has been extensively used by e.g., Cohen et al.
(2011, 2014); Alvarado-Gómez et al. (2016); do Nascimento et al. (2016);
Garraffo et al. (2016, 2017); Dong et al. (2017); Vidotto & Donati (2017);
Dong et al. (2018); Kavanagh et al. (2021). We refer the reader to those
references and to van der Holst et al. (2014) for the complete model
description. From the three-dimensional steady-state coronal solution, the
parameters along the orbit of the planet are extracted and are used to drive
the global magnetosphere model (see Section 2.3). The parameters are extracted
with a resolution of 3 degrees along the circular orbit (about 3.5 hours).

###  2.2 CME Model

In addition to the steady-state SW conditions, we calculate Ly$\alpha$
absorption profile variations due to a CME that hits the planet AU Mic b. The
CME is obtained in the coronal model by superimposing an unstable flux-rope
onto the steady state solution. Due to its instability, the flux-rope erupts
and its evolution, propagation, and interaction with the ambient SW is
simulated in a time dependent manner. The CME conditions near the planet are
recorded as a function of time, and they are also used to drive the global
magnetosphere model (similar to the orbital SW conditions). Here, we simulate
90 minutes of the CME as it passes through the location of AU Mic b.

The prescription for the unstable flux-rope is given by Titov & Démoulin
(1999) and it has been used for simulations of CMEs at Earth (e.g., Jin et
al., 2013) and for simulations of stellar CMEs (Alvarado-Gómez et al., 2019).
The bipolar flux-rope parameters include its location on the photosphere, tilt
angle, separation length, internal radius, total plasma mass, and stored
energy. The latter is controlled by a current term in the flux-rope formalism.
The location of the flux-rope was determined using data from Wisniewski et al.
(2019), while the CME mass, $M_{CME}=10^{20}~{}g$, and energy,
$E_{CME}=10^{36}~{}ergs$, were fit to the best current CME candidate event on
AU Mic (Moschou et al., 2019, see also Katsova et al. 1999; Cully et al.
1994). The complete description of the CME simulation can be found in
Alvarado-Gómez et al. (2022).

###  2.3 Global Magnetosphere Model

In order to study the impact of the stellar environment conditions on the
Ly$\alpha$ profile, we use the BATSRUS Global Magnetosphere (GM) MHD code
(Powell et al., 1999; Tóth et al., 2005, 2012). The model is constructed with
a spherical inner boundary at the planetary surface, and a cartesian grid that
extends between $+100R_{p}$ to $-200R_{p}$ along the star-planet line
($X_{GSE}$, day to night side), and between $+130R_{p}$ to $-130R_{p}$ in the
other directions ($Y_{GSE}$ and $Z_{GSE}$). Here, GSE is the Geocentric solar
ecliptic coordinate system (Kivelson & Russell, 1995). We adopt a planetary
radius $R_{p}=25,930$ km ($0.35R_{J}$, $4R_{E}$, Martioli et al., 2021). The
grid size near the inner boundary is $\Delta x=0.1R_{p}$.

The space environment conditions for the orbital SW or the time-dependent CME
propagation near the planet are imposed on the outer face of the GM domain
which faces the star. The conditions are given as time-dependent upstream
conditions, and the time-dependent GM simulation provides the response of the
escaping atmosphere and magnetosphere to variations in these upstream
conditions.

We are interested in capturing the variations of the escaping atmosphere and
its observable signature (i.e., the Ly$\alpha$ profile) due to the changing
upstream conditions. Thus, we impose a prescribed inner boundary condition for
the plasma density and temperature which results in a strong outflow due to
the pressure gradient. We stress that we do not attempt to simulate an actual
photoevaporative outflow and its acceleration at lower parts of the
atmosphere. We do not even specify an outflow velocity (set to zero) at the
boundary, and only specify temperature and density. Our approach enables us to
obtain an outflow that is strong enough to interact with the incoming SW/CME,
as well as to produce a Ly$\alpha$ absorption signature. The prescribed
boundary conditions are number density of $n=3\cdot 10^{8}~{}cm^{-3}$ and
temperature of $T=10^{4}~{}K$, resulting in an outflow mass-loss rates shown
in Table 1. The mass-loss rates are calculated over a sphere just slightly
above the spherical inner boundary (located at $1.2R_{p}$), and are of the
order of $10^{10}~{}g~{}s^{-1}$. This falls within the low range of mass-loss
rates of $10^{10}-10^{12}~{}g~{}s^{-1}$ predicted in hot Jupiters (Owen,
2019). Of course, the mass loss rate is controlled mostly by the base density
value used in the simulations.

We stress that our goal here is to drive a strong outflow to investigate the
impact of the stellar wind on the outflow and its Ly$\alpha$ signature. We do
not aim to predict an actual observed profile for AU Mic b. Thus, we apply a
prescribed pressure at the inner boundary which leads to an outflow from it.
This is different from one-dimensional models of hydrodynamic escape that use
the energy-limited stellar radiation to heat the base of the atmosphere and to
create the high base pressure along with a very low pressure at the top
boundary. Kubyshkina et al. (2018) performed a grid of hydrodynamic models for
different exoplanets. AU Mic b is probably most similar to the $Pd2$ case from
Kubyshkina et al. (2018), except that here we set the base temperature to
$10^{4}~{}K$, with a mass-loss rate of $10^{10}~{}g~{}s^{-1}$, while Case
$Pd2$ has a temperature of $4370~{}K$ and a mass-loss rate of $2\cdot
10^{9}~{}g~{}s^{-1}$. As a sanity check, we re-run our Case 1 using the
temperature of $4370~{}K$ and we obtained a mass-loss rate of $4\cdot
10^{9}~{}g~{}s^{-1}$, which is more or less consistent with Kubyshkina et al.
(2018). Since we aim to demonstrate that the Ly$\alpha$ signature is sensitive
to stellar wind variations, overestimating the outflow here does not matter,
since if the mass-loss rate would have been lower, it would be even more
sensitive to the stellar wind variations. Figure 1 compares the Ly$\alpha$
absorption images and line profiles for the two prescribed mass-loss rates. It
can be seen that a lower mass-loss rate results in a less visible absorption
(see the following section for a detailed description of how these images are
produced).

A planetary magnetic field of $B=-0.3~{}G$ (Earth-like) is assumed in all
cases. We leave the study of the impact of different planetary magnetic field
strengths on the planetary outflow for a future investigation (see Section
4.4).

###  2.4 Synthetic Ly$\alpha$ Profile

We follow the formalism by Tasitsiomi (2006), Khodachenko et al. (2017), and
Harbach et al. (2021) in order to obtain the Ly$\alpha$ synthetic observables.
First, we define a line-of-sight (LOS) which points in the positive $X_{GSE}$
direction, which extends from the back of the GM domain (behind the planet)
towards the star. We calculate the intensity of the Ly$\alpha$ emission,
$I_{Ly}(v)$, which comes out from the star and passes through the planetary
atmosphere, for a given velocity, $v$:

$I_{Ly}(v)=\int\int I_{Ly0}(y,z,v)\exp{[-\tau(y,z,v)]}dydz,$ (1)

with $I_{Ly0}(y,z,v)$ being the background stellar intensity, and
$\tau(y,z,v)$ being the optical depth. A strong absorption leads to an
increase in the optical depth, and a reduction in the Ly$\alpha$ intensity at
the observation point. The optical depth is obtained by integrating the column
density and the Ly$\alpha$ cross section, $\sigma_{Ly(v)}$:

$\tau(y,z,v)=\int_{LOS}0.5\cdot N_{H}(x)\sigma_{Ly}(v)dx.$ (2)

Here $N_{H}(x)$ is the local density at the cell along the LOS, and we assume
a 50% ionization of the gas, hence the factor of 0.5 (Owen & Mohanty, 2016).

The cross section is approximated as:

$\displaystyle\sigma_{Ly}(v)=5.8\times
10^{-14}\sqrt{10^{4}K/T(K)}\cdot\exp{(-b^{2})}$ (3) $\displaystyle+2.6\times
10^{-19}\left[\frac{100km~{}s^{-1}}{v-v_{x}}\right]^{2}q(b^{2}),$ (4)

where

$b=\frac{v-v_{x}}{\sqrt{2kT/m_{p}}}$ (5)

with $k$ is the Boltzmann constant, and $m_{p}$ is the proton mass. The
coefficient $q$ is given by:

$q(b^{2})=\begin{cases}\frac{21+b^{2}}{1+b^{1}}z(b^{2})(0.1117+z(b^{2}))&z(b^{2})>0\\\
\times[4.421+z(b^{2})(5.674z(b^{2})-9.207)]&\\\ 0&z(b^{2})<0,\end{cases}$ (6)

with

$z(b^{2})=\frac{b^{2}-0.855}{b^{2}+3.42}.$ (7)

The stellar background Ly$\alpha$ flux, $I_{Ly0}$, is assumed to be a Gaussian
with a full width at half maximum (FWHM) of $148.6~{}km~{}s^{-1}$, centered at
$1215.67\AA$, and peaking at $10^{-11}~{}ergs~{}cm^{-2}~{}s^{-1}$ (Schneider
et al., 2019).

For each of the GM MHD solutions, we perform the LOS integration for each line
on the grid in the $X_{GSE}$ direction for a given velocity. The velocity grid
runs from $-300~{}km~{}s^{-1}$ to $+300~{}km~{}s^{-1}$ with a resolution of
$10~{}km~{}s^{-1}$. Each LOS integration provides a pixel in a
$Y_{GSE}$-$Z_{GSE}$ image. At the end of the procedure, we obtain an intensity
image for every velocity, $I_{Ly}(y,z,v)=exp{(-\tau)}$. The image for each
velocity is multiplied by a background stellar disk that covers a circle of
$18R_{p}$ (AU Mic’s radius in units of AU Mic b’s radius), and the result is
integrated over the image to give the total flux at a given velocity (or
Doppler shift), i.e., Eq. 1. Multiplying this flux with the background stellar
flux Gaussian gives the absorption at each velocity, and the Ly$\alpha$ line
Doppler profile.

##  3 Results

###  3.1 Background Stellar Wind Conditions

Figure 2 shows the SW conditions along the orbit of AU Mic b. It shows the
orbital variations of the number density, $n$, the SW speed, $u$, the SW
magnetic field strength, $B$, the SW Alfvénic Mach number, $M_{A}=u/u_{A}$,
the SW dynamic pressure, $p$ (normalized to typical solar wind conditions at 1
AU with $n=5~{}cm^{-3}$ and $u=500~{}km~{}s^{-1}$), and plasma temperature,
$T$. The local Alfvén speed is $u_{A}=B/\sqrt{4\pi\rho}$, with $\rho=nm_{p}$
being the mass density, and $m_{p}$ the proton mass. The SW velocity and
magnetic field vectors, as well as the number density and temperature are used
to drive the GM model. The stellar wind conditions vary strongly with
azimuthal angle. In particular, AU Mic b will reside in a sub-Alfvénic regime
for most of the orbit, but transitions into a super-Alfvénic regime twice
during the crossing of the astrospheric current sheet.

The SW speeds are not so different from those observed near the Earth.
However, the number density and the magnetic field strength are orders of
magnitude higher due to the close proximity of the orbit to the star (typical
values for the Earth are $n=1-100~{}cm^{-3}$ and $B=1-500~{}nT$, for ambient
solar wind and CME conditions111https://cdaweb.gsfc.nasa.gov/istp_public/).
The orbital period of AU Mic b is about $8.5~{}d$ (Plavchan et al., 2020;
Martioli et al., 2021). Thus, the crossing of the astrospheric current sheet
(low to high density and high to low magnetic field/wind speed) takes about
two days, and the transition between the sub- to super-Alfvénic SW takes less
than a day. Cohen et al. (2014) have shown that such quick transitions may
lead to a rapid global change in the magnetospheric structure from an Alfvén
wings (Neubauer, 1980, 1998), Io-like topology to a stretched, Earth-like
topology, which leads to strong variations in the SW-planetary outflow
interaction.

Here, we focus on three azimuthal angles in the simulations that we use as
representative points along the orbit of AU Mic b. We investigate a case of
the sub-Alfvénic region at the beginning of the orbit (Case 1), a case of the
super-Alfvénic point (Case 2), and a case of the second sub-Alfvénic region
(Case 3). We assume that all other orbital locations behave similarly to these
locations.

Figure 3 shows the results for Cases 1–3. The top panel shows the mid-transit
absorption ($e^{-\tau}$) images integrated over all velocities. It can be seen
that the Ly$\alpha$ absorption corresponds to the cooler (under 50,000K),
denser material that occupies the inner regions of the planetary
magnetosphere. The middle panel of Figure 3 shows the Ly$\alpha$ Doppler
profile for the three cases (left), and the ratio of their absorbed to the
non-absorbed (stellar background) profiles (right). These plots show
significant absorption in the blue wing of the profile, ranging between
10-20%, and peaking at velocities between $-40$ to $-100~{}km~{}s^{-1}$
($1215.25-1215.5$Å).

In the bottom panel of Figure 3, we show the corresponding three-dimensional
density and magnetic field structure of the planetary magnetosphere for the
three orbital epochs. It shows that initially (Case 1), the planet resides in
a sub-Alfvénic SW with two magnetospheric lobes extended towards and away from
the star, with an angle that depends on the particular SW magnetic field
vector at that time. A relatively high density is seen close to the planet due
to the planetary outflow from the inner boundary. When the planet moves to the
super-Alfvénic SW (Case 2, current sheet crossing), the SW density increases,
stretching the magnetospheric lobes behind the planet, forming a bow shock at
the day side of the magnetosphere, and pushing against the planetary outflow.
When the planet moves back to the lower density, sub-Alfvénic SW region (Case
3), the initial structure is recovered (even though Cases 1 and 3 are not
identical).

Figure 4 shows the predicted Ly$\alpha$ lightcurve. It shows the Ly$\alpha$
flux integrated over all velocities as the planet passes in front of the
stellar disk. Our predicted transit lasts for about 4-5h, which is consistent
with e.g., Plavchan et al. (2020). Figure 4 and the middle panel of Figure 3
clearly shows that the Ly$\alpha$ transit changes significantly along the
orbit due to the different SW conditions that the planet experiences. The blue
wing and the overall absorption is less than 10% in Case 1, then it increases
to almost 20% in Case 2, then it decreases again to about 13-14% in Case 3.
These absorption variations occur on a time scale of less than a day, and they
are clearly due to the SW variations as our planetary outflow driver at the
inner boundary is uniform and steady. It is also clear from the top panel of
Figure 3 that the SW conditions change the overall (density) structure of the
inner magnetosphere, which is the region that is also responsible for the
absorption of the stellar Ly$\alpha$ radiation.

###  3.2 CME Conditions

Figure 5 shows the SW conditions for the CME event hitting AU Mic b. These
conditions were extracted with a 1 minute cadence over the course of 90
minutes. The CME front is clearly seen in the Alfvénic Mach number plot, where
it arrives at the planet after only about 20 minutes. This is not surprising
due to the very fast CME speed of over $8000~{}km~{}s^{-1}$. The CME front is
also indicated by the sudden increase in density, magnetic field, and speed.
The post-front CME conditions occupy the vicinity of the planet for the rest
of the simulation.

Figure 6 is similar to Figure 3 but it shows the results for the CME event at
$t=0$, $24$, and $90$ minutes. The choice of $t=24$ min is due to the fact
that the CME front arrives at the edge of the GM simulation domain after 20
minutes, but it takes 4 more minutes to see the impact of the CME on the
absorption profile when the CME reaches close proximity to the planet. It can
be seen that prior to the arrival of the CME at the initial state, a strong
absorption occurs in a similar manner to Case 1, as the SW conditions at that
time are sub-Alfvénic. When the CME arrives after 24 minutes, and throughout
the CME event up to $t=90~{}min$, the signature completely disappears from
both the absorption images and the Ly$\alpha$ Doppler profile. The three-
dimensional plots at the bottom panel of Figure 6 clearly show that the CME
plasma takes over most of the simulation domain, replacing the magnetospheric
material that occupied the domain prior to the CME arrival.

##  4 Discussion

The results presented above highlight the conclusions of Harbach et al.
(2021): for close-in planets the stellar wind is instrumental in shaping
atmospheric outflows and can dictate the form and variability of absorption
signatures. Such variability can be large. Different stellar wind and magnetic
field conditions with azimuthal angle imply that strong variations in outflow
signatures will occur over an orbit, potentially leading to observable
variations in absorption signatures during a single transit on timescales as
short as a few hours.

Regardless of any secular evolution of the stellar surface magnetic field that
could lead to changes in the wind conditions and atmospheric outflow
absorption signatures, the occurrence of transit-to-transit variations can be
expected.

###  4.1 Stellar-wind Driven Variations of the Ly$\alpha$ Observables

One of the aims of this study is to shed light on the relation between the
three-dimensional physical system and the one-dimensional, Ly$\alpha$ signal
that is often the only observed diagnostic. In interpreting observed
Ly$\alpha$ profiles as an indicator of a strong atmospheric outflow, many
assumptions are made regarding the space environment and stellar wind near the
planet, the particular modeling approach and assumptions, and the radiative
transfer processes. It is important to stress that these interpretations are
generally model dependent (or “assumptions dependent” on the practical level).
Our simulations provide a number of important features that could constrain
these assumptions, in the context of the three-dimensional system.

As first proposed by Vidal-Madjar et al. (2003), an excess of blue-shifted
Ly$\alpha$ absorption during transit as compared with the out of transit flux
indicates a flow of neutral hydrogen in the direction of the observer.
Assuming that the neutral hydrogen is concentrated in the planetary
atmosphere, the observed blue-shift of about $100~{}km~{}s^{-1}$ indicates a
strong outflow of neutral hydrogen from the planet with that speed. The strong
absorption also requires a quite high density, with values that are likely to
be possible only at the top of the planetary atmosphere (below an altitude of
a few thousand km). In contrast, models for neutral hydrogen hydrodynamic
escape from hot Jupiters (e.g., Murray-Clay et al., 2009; Owen & Jackson,
2012; Tripathi et al., 2015; McCann et al., 2019) predicted escape velocities
in the range of $5-20~{}km~{}s^{-1}$.

Despite the high mass-loss rate we obtain in our simulations, the outflow
speed near the planet is much less than $20~{}km~{}s^{-1}$, which is
consistent with the hydrodynamic escape models referenced above. This is far
lower than the negative $100~{}km~{}s^{-1}$ speed associated with the
Ly$\alpha$ observations. We draw a particular attention to an important aspect
of the nature of the hydrodynamic escape process. In its essence, hydrodynamic
escape is driven by a pressure gradient between the planet and space, and the
greater this pressure gradient is the greater the outflow. Indeed, in his
original paper on the solar wind, magnetic field and interstellar medium,
Parker (1958) showed that the low pressure estimated for the interstellar
medium was consistent with hydrodynamic outflow of a $3\times 10^{6}$ K solar
corona. Such a setting is assumed in many models for hydrodynamic escape
(mentioned in Section 1). In our simulations here, we impose a high pressure
at the inner boundary via high density and temperature that were obtained from
the literature. If our simulated exoplanet would be located far from the star,
then a high escape rate would be self-consistently generated in the simulation
via the pressure gradient. However, in the case of close-in exoplanets, space
is not that empty (i.e., the pressure is not necessarily low). The ambient SW
density near AU Mic b is 3-4 orders of magnitude higher than that near the
Earth, and the temperature is approaching the coronal temperature of more than
1MK. Thus, the pressure at “infinity” is not much smaller, and in some
situations may actually be higher than that of the planet inner boundary.

Figure 7 shows the radial velocity close to the planet (magnitude and
streamlines displayed on a $x-z$ slice), together with the $U_{x}$ velocity
component (the velocity component towards the observer) and the thermal
pressure profile as a function of distance from the planet. The latter were
extracted from the simulation at the back of the planet (night side) in the
negative $\hat{x}$ direction. It can be seen that the velocity close to the
planet is much less than even $15~{}km~{}s^{-1}$, and it gets to higher
negative values far from the planet, indicating a SW plasma, not planetary
material. It also shows that the thermal pressure near the planet does not
drop much: by an order of magnitude or so in the sub-Alfvénic cases (1 and 3)
within 3-4 $R_{p}$, and by an even smaller factor in the super-Alfvénic case
(2). The weak pressure gradient may be even weaker if we consider the total
pressure, which is the sum of the thermal, dynamic, and magnetic pressures.

Figures 8 and 9 show the different pressures extracted along the planet-star
line (sub-stellar line) for Cases 1-3 and the CME epochs, respectively. The
magnetic pressure, which is very small in the SW near the Earth, is much more
significant in the case of AU Mic (and close-in planets in general). It is
actually the dominant pressure term in some parts. This emphasises the
importance of including magnetic effects in studies of planetary atmospheric
escape. When considering the total pressure, the pressure gradient becomes
even more moderate than that of the thermal pressure alone.

The reduction of the pressure far from the planet is taken to the extreme in
our CME simulations (Figures 5 and 9). Here, it is clearly seen how the
increased SW pressure completely shuts down and suppresses the escape of
planetary material via reduction of the pressure gradient. On the one hand,
this is similar to the escape suppressing by a strong SW, as proposed by
Shaikhislamov et al. (2016); Vidotto & Cleary (2020). However, the escape
suppression clearly occurs due to the overall reduction in the pressure
gradient surrounding the planet. In general, our simulations show that for
close-in exoplanet simulations of hydrodynamic escape, the upper boundary
conditions for the pressure (or the pressure at ”infinity”) need to account
for the rather high SW pressure, including the magnetic pressure.

###  4.2 Mass-loss Due to the CME Event

One important aspect of CME events on AU Mic is their possible impact on
atmospheric escape. Our GM simulation cannot capture any actual impact on the
escape in terms of atmospheric heating or an acceleration of escaping
particles (see e.g., Cohen et al., 2014; Garcia-Sage et al., 2017). However,
we can still look at the effect on the outflow of the CME as it passes the
planet.

Our planetary outflow is driven by setting a constant, high pressure at the
inner boundary, which is higher than that of the surrounding space, creating a
mass-loss rate of about $10^{10}~{}g~{}s^{-1}$. However, as the CME arrives,
the total pressure in the vicinity of the planet becomes larger than the
pressure at the inner boundary, so the outflow is shut down.

Table 1 shows the mass-loss rates during the CME event. Initially, the outflow
is $2\times 10^{10}~{}g~{}s^{-1}$. However, when the CME arrives at
$t=24~{}min$, we see an inflow due to the incoming CME material pushing
against the outflowing planetary material. Finally, after $90~{}min$, the CME
penetrates through the sphere and the inflowing mass flux through it is
completely associated with the CME itself. Since the planetary outflow fills
the region near the planet and perhaps parts of its magnetosphere (if it
exists), it is expected that a significant amount of mass will be stripped by
the CME. Estimating this amount of mass-loss requires taking into account the
time scales involved in filling the magnetosphere and stripping it, as well as
some more consistent escape mechanisms than the one we use here. Thus, we
leave this investigation for future studies.

###  4.3 The Escaping Neutral Atmosphere

Although the escaping neutral hydrogen should not be affected by the magnetic
field carried by the stellar wind, one can ask whether the neutral hydrogen
atoms travel freely without interacting with the ionized stellar wind and
magnetospheric plasma.

The stellar wind particle density from our simulations in the vicinity of the
planet (Figure 8) is in the range $n\sim 10^{4}$-$10^{5}$ cm-3. The cross-
section for collisions between protons and neutral hydrogen atoms is dominated
by the elastic collision term and is of the order $\sim 10^{-14}$ cm2 for the
temperatures of interest here (Hunter & Kuriyan, 1977). The mean free path for
the low kinetic energy neutrals near the planet between collisions with
stellar wind protons is then of the order of 10,000-100,000 km, or similar to
the planetary radius of $R_{p}=25,930$ km (Martioli et al., 2021). Thus, it is
reasonable to believe that the neutral material will follow the magnetic field
structure to some extent via the collisions with the ions. The inclusion of
the magnetized environment (i.e., MHD approach) is then crucial for rigorous
study of the interaction between the escaping atmosphere and the space
environment. McCann et al. (2019) estimated that the timescale for ionization
of the escaping neutrals is about 4 hours. With a maximum escape velocity of
$20~{}km~{}s^{-1}$, these neutrals could reach a distance of about 10
planetary radii before they get ionized. However, it is more likely that their
escape velocity is much less than the above upper limit, and that the neutrals
will collide with other neutrals/ions much closer to the planet. All of the
above suggest that the scenario of a neutral escaping atmosphere that freely
expands to large distances from the planet is very unlikely, especially in a
close-in orbits, where the photoevaporation is expected to be very strong. It
is important to note here that instead of a pure hydrodynamic evaporation,
other escape mechanisms may lead to the significant loss of the planetary
gaseous envelop (see recent review by Gronoff et al., 2020).

The Ly$\alpha$ absorption we obtain using our assumption of a 50% ionization
fraction of our single-fluid plasma provides an absorption profile with 10-20%
flux reduction of the blue wing. In order to further investigate the size and
shape of the escaping neutral atmosphere, we repeat our simulation for Cases
1-3 and include a neutral hydrogen fluid. The neutral fluid is imposed at the
inner boundary in the same manner as the single-fluid plasma, and it expands
from the inner boundary. It is important to note that in this numerical
experiment, the plasma and neutral fluids are completely decoupled, and no
interaction between them is included (including the collisional analysis
described above). Of course, in this formalism the neutral fluid is not
affected by the magnetic field and it is purely hydrodynamic.

The results of our simulations for the neutral hydrogen fluid are shown in
Figure 10. The display is the same as Figure 3. One can see that the neutrals
expand significantly far from the planet, beyond 10 planetary radii as
predicted by the conventional view of an extended neutral atmosphere in close-
in gas planets. However, when using this extended neutral atmosphere to
calculate the Ly$\alpha$ absorption profile, we find that the blue wing is
almost completely absorbed (while the red wing is not). Even though the
neutral gas should not be affected by the magnetic field, it is relatively
coupled to the stellar wind. Thus, in the positive x direction (red-shift,
away from the observer), the neutrals’ motion is stopped by the incoming
stellar wind, while in the negative x direction (blue-shift, towards the
observer), the neutrals are carried by the ion motion, as well as their own,
slow motion. Our neutral experiment suggests that a very extended neutral
atmosphere will be manifested in a greater blue-wing absorption. Of course,
our simulation here is over-simplified. A better treatment would include a
detailed ion-neutral interaction, self-consistent ionization, and ion escape
in a multi-species manner. We plan to implement these processes in future
extension of this work.

The low escape velocities we obtain are consistent with the concept that the
Ly$\alpha$ absorption profile is blue-shifted mostly due to the stellar wind -
outflow interaction and not due to an extreme planetary outflow by itself
(Murray-Clay et al., 2009; McCann et al., 2019; Carolan et al., 2020).
Moreover, our MHD simulations have even stronger coupling between the SW and
the planetary magnetosphere than these past hydrodynamic simulations. In our
simulations, the fast SW (and even faster CME) compress near the planet,
leading to a high enough density to contribute to the synthetic absorption
line profile. In reality, the ionized SW may need to undergo charge-exchange
in order to contribute to the Ly$\alpha$ absorption as proposed by Holmström
et al. (2008); Murray-Clay et al. (2009); McCann et al. (2019); Debrecht et
al. (2022). Such charge-exchange is seen in many bodies in our own solar
system, and it is expected to be stronger in close-in exoplanets, where the SW
density is orders of magnitude higher than the SW density around solar system
bodies.

The possibility that the observed Ly$\alpha$ absorption is due to the SW and
not the planetary escaping atmosphere is potentially problematic for observing
the planetary atmosphere to extend our knowledge about atmospheric loss and,
by extension, exoplanet evolution. However, this opens the door to use these
observations to constrain the SW itself. Our work also raises the possibility
of a future study about how the planetary magnetic field affects the
Ly$\alpha$ profile (here we only use a single planetary field value of
$0.3~{}G$). If Ly$\alpha$ observations could provide information about the
planetary magnetic field, then its contribution to exoplanet evolution may
still be significant.

###  4.4 the Possible Role of the Planetary Magnetic Field Strength

The possible role of the planetary magnetic field on the escaping neutrals and
their pattern is not clear. It is reasonable to assume that a stronger
planetary field would probably reduce the sensitivity of the outflow to the
stellar wind conditions since a stronger field would keep a more steady,
dipolar structure close to the planet, where the outflow operates (see, e.g.,
Owen & Adams, 2014; Owen, 2019; Khodachenko et al., 2021).

The impact of the magnetic field on the escape itself is a completely open
question. The intuitive assumption is that a stronger planetary field would
protect the atmosphere from stellar wind stripping. However, at the upper
atmosphere, ion chemistry dominates, and the impact of the planetary field is
stronger. It has been shown that in some cases the strong field can actually
enhance the escape via wave-particle interaction heating (Strangeway et al.,
2010a, b, 2019). This is a very broad subject that is beyond the scope of our
paper, especially due to the fact that our model here cannot capture any of
these processes. This should be studied by a model for the upper atmosphere,
which includes ion chemistry, electrodynamics, and the planetary field.

A stronger field could potentially impact the amount of neutral Hydrogen,
which absorbs the Ly alpha flux. If the Ly alpha signature is produced by
charge-exchanged stellar wind, then a stronger field would push the incoming
stellar wind further from the planet, where the planetary neutral density is
much lower. We would expect that this would lead to a reduction in the Ly
alpha absorption.

##  5 Summary & Conclusions

In this paper, we simulate the response of the atmospheric outflow from AU Mic
to the time-varying conditions of the stellar wind along the orbit, and to a
stellar CME. We also calculated the predicted Ly$\alpha$ transit absorption
signature of this interaction. We find that:

  1. 1.

The Ly$\alpha$ emission is highly variable due to the highly-varying stellar
wind conditions (ranging between 10-20% absorption in our simulations) through
an orbit. Thus, Ly$\alpha$ observations represent a particular epoch, which
may change even over a fraction of the exoplanet orbit. Thus, one epoch alone
should not be used to deduce information about the global system, e.g., total
escape rate and the planetary magnetic field strength. Instead, multi-epoch
Ly$\alpha$ profiles could be useful for characterizing these parameters.

  2. 2.

The pressure gradient, which drives planetary atmosphere hydrodynamic escape,
may be weak in close-in exoplanets due to the relatively high pressure of the
nearby space environment. During CME events, the pressure gradient may be
completely removed, thereby suppressing atmospheric escape.

  3. 3.

Our results show a strong impact of the magnetic field and magnetic pressure
on the escaping atmospheric material and the Ly$\alpha$ profile. Thus, our
results show that an MHD, and not hydrodynamic treatment is crucial for the
study of escaping atmospheres and their interaction with the nearby space
environment. Ideally, a multi-fluid approach that couples ions and neutrals
would be the most complete.

  4. 4.

Our results support the claim that the high Doppler blue-shifted velocity is
likely attributed to the stellar wind sweeping the escaping neutral material,
or it might be due to charge-exchanged, neutralized fast stellar wind. It is
unlikely that the atmospheric neutral material reaches such velocities.

  5. 5.

Our CME simulations indicate a potential strong stripping of magnetospheric
material from the planet, including some of the neutral escaping atmosphere.
However, without a consistent model for the atmospheric escape, our
simulations lack the ability to properly quantify the rate of mass-loss per
CME event.

## acknowledgments

We than an unkown referee for her/his useful comments. This work is supported
by NASA XRP grant 80NSSC20K0840. FF was supported, in part, by NASA through
Chandra Theory Award Number $TM0-21001X$ issued by the Chandra X-ray
Observatory Center, which is operated by the Smithsonian Astrophysical
Observatory for and on behalf of NASA under contract NAS8-03060. LMH has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (grant agreement No.
853022, PEVAP). JJD was supported by NASA contract NAS8-39073 to the Chandra
X-ray Center and thanks the Director, P. Slane, for continuing advice and
support. Simulation results were obtained using the (open source) Space
Weather Modeling Framework, developed by the Center for Space Environment
Modeling, at the University of Michigan with funding support from NASA ESS,
NASA ESTO-CT, NSF KDI, and DoD MURI. The simulations were performed on NASA’s
Pleiades cluster (under SMD-20-52848317) and on the Massachusetts Green High
Performance Computing Center (MGHPCC) cluster.

## References

  * Altschuler & Newkirk (1969) Altschuler, M. D., & Newkirk, G. 1969, Sol. Phys., 9, 131, doi: 10.1007/BF00145734
  * Alvarado-Gómez et al. (2019) Alvarado-Gómez, J. D., Drake, J. J., Moschou, S. P., et al. 2019, ApJ, 884, L13, doi: 10.3847/2041-8213/ab44d0
  * Alvarado-Gómez et al. (2016) Alvarado-Gómez, J. D., Hussain, G. A. J., Cohen, O., et al. 2016, A&A, 588, A28, doi: 10.1051/0004-6361/201527832
  * Alvarado-Gómez et al. (2022) Alvarado-Gómez, J. D., Cohen, O., Drake, J. J., et al. 2022, arXiv e-prints, arXiv:2202.07949.  https://arxiv.org/abs/2202.07949
  * Batygin et al. (2016) Batygin, K., Bodenheimer, P. H., & Laughlin, G. P. 2016, ApJ, 829, 114, doi: 10.3847/0004-637X/829/2/114
  * Beaugé & Nesvorný (2013) Beaugé, C., & Nesvorný, D. 2013, ApJ, 763, 12, doi: 10.1088/0004-637X/763/1/12
  * Carolan et al. (2020) Carolan, S., Vidotto, A. A., Plavchan, P., Villarreal D’Angelo, C., & Hazra, G. 2020, MNRAS, 498, L53, doi: 10.1093/mnrasl/slaa127
  * Cohen et al. (2014) Cohen, O., Drake, J. J., Glocer, A., et al. 2014, ApJ, 790, 57, doi: 10.1088/0004-637X/790/1/57
  * Cohen et al. (2011) Cohen, O., Kashyap, V. L., Drake, J. J., et al. 2011, ApJ, 733, 67, doi: 10.1088/0004-637X/733/1/67
  * Cully et al. (1994) Cully, S. L., Fisher, G. H., Abbott, M. J., & Siegmund, O. H. W. 1994, ApJ, 435, 449, doi: 10.1086/174827
  * Debrecht et al. (2022) Debrecht, A., Carroll-Nellenback, J., Frank, A., et al. 2022, MNRAS, doi: 10.1093/mnras/stac112
  * do Nascimento et al. (2016) do Nascimento, J. D., J., Vidotto, A. A., Petit, P., et al. 2016, ApJ, 820, L15, doi: 10.3847/2041-8205/820/1/L15
  * Dong et al. (2018) Dong, C., Jin, M., Lingam, M., et al. 2018, Proceedings of the National Academy of Science, 115, 260, doi: 10.1073/pnas.1708010115
  * Dong et al. (2017) Dong, C., Lingam, M., Ma, Y., & Cohen, O. 2017, ApJ, 837, L26, doi: 10.3847/2041-8213/aa6438
  * Ekenbäck et al. (2010) Ekenbäck, A., Holmström, M., Wurz, P., et al. 2010, ApJ, 709, 670, doi: 10.1088/0004-637X/709/2/670
  * Erkaev et al. (2007) Erkaev, N. V., Kulikov, Y. N., Lammer, H., et al. 2007, A&A, 472, 329, doi: 10.1051/0004-6361:20066929
  * Garcia-Sage et al. (2017) Garcia-Sage, K., Glocer, A., Drake, J. J., Gronoff, G., & Cohen, O. 2017, ApJ, 844, L13, doi: 10.3847/2041-8213/aa7eca
  * Garraffo et al. (2016) Garraffo, C., Drake, J. J., & Cohen, O. 2016, ApJ, 833, L4, doi: 10.3847/2041-8205/833/1/L4
  * Garraffo et al. (2017) Garraffo, C., Drake, J. J., Cohen, O., Alvarado-Gómez, J. D., & Moschou, S. P. 2017, ApJ, 843, L33, doi: 10.3847/2041-8213/aa79ed
  * Gronoff et al. (2020) Gronoff, G., Arras, P., Baraka, S., et al. 2020, Journal of Geophysical Research (Space Physics), 125, e27639, doi: 10.1029/2019JA027639
  * Harbach et al. (2021) Harbach, L. M., Moschou, S. P., Garraffo, C., et al. 2021, ApJ, 913, 130, doi: 10.3847/1538-4357/abf63a
  * Hazra et al. (2021) Hazra, G., Vidotto, A. A., Carolan, S., D’Angelo, C. V., & Manchester, W. 2021, MNRAS, doi: 10.1093/mnras/stab3271
  * Holmström et al. (2008) Holmström, M., Ekenbäck, A., Selsis, F., et al. 2008, Nature, 451, 970, doi: 10.1038/nature06600
  * Hunter & Kuriyan (1977) Hunter, G., & Kuriyan, M. 1977, Proceedings of the Royal Society of London Series A, 353, 575, doi: 10.1098/rspa.1977.0051
  * Jin et al. (2013) Jin, M., Manchester, W. B., van der Holst, B., et al. 2013, ApJ, 773, 50, doi: 10.1088/0004-637X/773/1/50
  * Katsova et al. (1999) Katsova, M. M., Drake, J. J., & Livshits, M. A. 1999, ApJ, 510, 986, doi: 10.1086/306587
  * Kavanagh et al. (2021) Kavanagh, R. D., Vidotto, A. A., Klein, B., et al. 2021, MNRAS, 504, 1511, doi: 10.1093/mnras/stab929
  * Khodachenko et al. (2021) Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., et al. 2021, MNRAS, 507, 3626, doi: 10.1093/mnras/stab2366
  * Khodachenko et al. (2017) —. 2017, ApJ, 847, 126, doi: 10.3847/1538-4357/aa88ad
  * Kiraga (2012) Kiraga, M. 2012, Acta Astron., 62, 67.  https://arxiv.org/abs/1204.3825
  * Kislyakova et al. (2014) Kislyakova, K. G., Johnstone, C. P., Odert, P., et al. 2014, A&A, 562, A116, doi: 10.1051/0004-6361/201322933
  * Kislyakova et al. (2019) Kislyakova, K. G., Holmström, M., Odert, P., et al. 2019, A&A, 623, A131, doi: 10.1051/0004-6361/201833941
  * Kivelson & Russell (1995) Kivelson, M. G., & Russell, C. T. 1995, Introduction to Space Physics, 586
  * Klein et al. (2021) Klein, B., Donati, J.-F., Moutou, C., et al. 2021, MNRAS, 502, 188, doi: 10.1093/mnras/staa3702
  * Kochukhov & Reiners (2020) Kochukhov, O., & Reiners, A. 2020, ApJ, 902, 43, doi: 10.3847/1538-4357/abb2a2
  * Kubyshkina et al. (2022) Kubyshkina, D., Vidotto, A. A., Villarreal D’Angelo, C., et al. 2022, MNRAS, 510, 2111, doi: 10.1093/mnras/stab3594
  * Kubyshkina et al. (2018) Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018, A&A, 619, A151, doi: 10.1051/0004-6361/201833737
  * Lecavelier des Etangs et al. (2004) Lecavelier des Etangs, A., Vidal-Madjar, A., McConnell, J. C., & Hébrard, G. 2004, A&A, 418, L1, doi: 10.1051/0004-6361:20040106
  * Linsky et al. (2010) Linsky, J. L., Yang, H., France, K., et al. 2010, ApJ, 717, 1291, doi: 10.1088/0004-637X/717/2/1291
  * Lundkvist et al. (2016) Lundkvist, M. S., Kjeldsen, H., Albrecht, S., et al. 2016, Nature Communications, 7, 11201, doi: 10.1038/ncomms11201
  * Martioli et al. (2021) Martioli, E., Hébrard, G., Correia, A. C. M., Laskar, J., & Lecavelier des Etangs, A. 2021, A&A, 649, A177, doi: 10.1051/0004-6361/202040235
  * Matsakos & Königl (2016) Matsakos, T., & Königl, A. 2016, ApJ, 820, L8, doi: 10.3847/2041-8205/820/1/L8
  * Mazeh et al. (2016) Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75, doi: 10.1051/0004-6361/201528065
  * McCann et al. (2019) McCann, J., Murray-Clay, R. A., Kratter, K., & Krumholz, M. R. 2019, ApJ, 873, 89, doi: 10.3847/1538-4357/ab05b8
  * Moschou et al. (2019) Moschou, S.-P., Drake, J. J., Cohen, O., et al. 2019, The Astrophysical Journal, 877, 105, doi: 10.3847/1538-4357/ab1b37
  * Murray-Clay et al. (2009) Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23, doi: 10.1088/0004-637X/693/1/23
  * Neubauer (1980) Neubauer, F. M. 1980, J. Geophys. Res., 85, 1171, doi: 10.1029/JA085iA03p01171
  * Neubauer (1998) —. 1998, J. Geophys. Res., 103, 19843, doi: 10.1029/97JE03370
  * Oklopčić & Hirata (2018a) Oklopčić, A., & Hirata, C. M. 2018a, ApJ, 855, L11, doi: 10.3847/2041-8213/aaada9
  * Oklopčić & Hirata (2018b) —. 2018b, ApJ, 855, L11, doi: 10.3847/2041-8213/aaada9
  * Owen (2019) Owen, J. E. 2019, Annual Review of Earth and Planetary Sciences, 47, 67, doi: 10.1146/annurev-earth-053018-060246
  * Owen & Adams (2014) Owen, J. E., & Adams, F. C. 2014, MNRAS, 444, 3761, doi: 10.1093/mnras/stu1684
  * Owen & Jackson (2012) Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931, doi: 10.1111/j.1365-2966.2012.21481.x
  * Owen & Lai (2018) Owen, J. E., & Lai, D. 2018, MNRAS, 479, 5012, doi: 10.1093/mnras/sty1760
  * Owen & Mohanty (2016) Owen, J. E., & Mohanty, S. 2016, MNRAS, 459, 4088, doi: 10.1093/mnras/stw959
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664, doi: 10.1086/146579
  * Plavchan et al. (2020) Plavchan, P., Barclay, T., Gagné, J., et al. 2020, Nature, 582, 497, doi: 10.1038/s41586-020-2400-z
  * Powell et al. (1999) Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., & Zeeuw, D. L. D. 1999, Journal of Computational Physics, 154, 284 , doi: http://dx.doi.org/10.1006/jcph.1999.6299
  * Russell (2009) Russell, C. T. 2009, Space Physics: An Introduction
  * Schneider et al. (2019) Schneider, A. C., Shkolnik, E. L., Barman, T. S., & Loyd, R. P. 2019, ApJ, 886, 19, doi: 10.3847/1538-4357/ab48de
  * Schrijver & Siscoe (2009) Schrijver, C. J., & Siscoe, G. L. 2009, Heliophysics: Plasma Physics of the Local Cosmos
  * Shaikhislamov et al. (2014) Shaikhislamov, I. F., Khodachenko, M. L., Sasunov, Y. L., et al. 2014, ApJ, 795, 132, doi: 10.1088/0004-637X/795/2/132
  * Shaikhislamov et al. (2016) Shaikhislamov, I. F., Khodachenko, M. L., Lammer, H., et al. 2016, ApJ, 832, 173, doi: 10.3847/0004-637X/832/2/173
  * Spake et al. (2018) Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nature, 557, 68, doi: 10.1038/s41586-018-0067-5
  * Strangeway et al. (2010a) Strangeway, R. J., Russell, C. T., & Luhmann, J. G. 2010a, in European Planetary Science Congress 2010, 334
  * Strangeway et al. (2019) Strangeway, R. J., Russell, C. T., Luhmann, J. G., Moore, T. E., & Barabash, S. V. 2019, in AGU Fall Meeting Abstracts, Vol. 2019, SM52B–05
  * Strangeway et al. (2010b) Strangeway, R. J., Russell, C. T., Luhmann, J. G., et al. 2010b, in AGU Fall Meeting Abstracts, Vol. 2010, SM33B–1893
  * Szabó & Kiss (2011) Szabó, G. M., & Kiss, L. L. 2011, ApJ, 727, L44, doi: 10.1088/2041-8205/727/2/L44
  * Tasitsiomi (2006) Tasitsiomi, A. 2006, ApJ, 645, 792, doi: 10.1086/504460
  * Titov & Démoulin (1999) Titov, V. S., & Démoulin, P. 1999, A&A, 351, 707
  * Tóth et al. (2005) Tóth, G., Sokolov, I. V., Gombosi, T. I., et al. 2005, Journal of Geophysical Research (Space Physics), 110, A12226, doi: 10.1029/2005JA011126
  * Tóth et al. (2012) Tóth, G., van der Holst, B., Sokolov, I. V., et al. 2012, Journal of Computational Physics, 231, 870, doi: 10.1016/j.jcp.2011.02.006
  * Tripathi et al. (2015) Tripathi, A., Kratter, K. M., Murray-Clay, R. A., & Krumholz, M. R. 2015, ApJ, 808, 173, doi: 10.1088/0004-637X/808/2/173
  * van der Holst et al. (2014) van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81, doi: 10.1088/0004-637X/782/2/81
  * Vidal-Madjar et al. (2003) Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J.-M., et al. 2003, Nature, 422, 143, doi: 10.1038/nature01448
  * Vidotto & Cleary (2020) Vidotto, A. A., & Cleary, A. 2020, MNRAS, 494, 2417, doi: 10.1093/mnras/staa852
  * Vidotto & Donati (2017) Vidotto, A. A., & Donati, J. F. 2017, A&A, 602, A39, doi: 10.1051/0004-6361/201629700
  * Wisniewski et al. (2019) Wisniewski, J. P., Kowalski, A. F., Davenport, J. R. A., et al. 2019, ApJ, 883, L8, doi: 10.3847/2041-8213/ab40bf

Case |  Mass-loss Rate [$g~{}s^{-1}$]
---|---
1 |  $-2\cdot 10^{10}$ (outflow)
2 |  $-2\cdot 10^{10}$ (outflow)
3 |  $-2\cdot 10^{10}$ (outflow)
CME Time [m] |  Mass-loss Rate [$g~{}s^{-1}$]
0 |  $-1.8\cdot 10^{10}$ (outflow)
24 |  $7\cdot 10^{8}$ (inflow)
90 |  $3\cdot 10^{9}$ (inflow)
Table 1: Ion and neutral mass-loss rates during the CME event. Figure 1:
Ly$\alpha$ line absorption (left) and absorption images (middle and right) for
planetary mas-loss rate of $4\cdot 10^{9}$ and $10^{10}~{}g~{}s^{-1}$. Figure
2: SW parameters as extracted from the AWSOM solution along the orbit of AU
Mic b and used to drive the GM model. Plots are for the SW number density
(top-left), SW speed (top-middle), SW magnetic field strength (top-right), and
SW Alfvénic Mach number (bottom-left), SW dynamic pressure (bottom-middle),
and SW temperature (bottom-right). SW conditions for cases 1-3 are marked by
the vertical black lines (left to right, respectively). Figure 3: Top:
Absorption ($e^{-\tau}$) images for the three cases integrated over all
velocities. Middle: Ly$\alpha$ Doppler-shift profile for the three cases
(left), and the ratio of these profiles to the un-absorbed profile (right).)
Bottom: the three-dimensional solution near the planet for Cases 1-3 (left,
middle, right, respectively). Color contours are for the number density where
selected magnetic field lines are also shown. The solid white line marks the
Alfvén surface (seen only in the super-Alfvénic Case 2). Figure 4: Synthetic
transit profile of the Ly$\alpha$ flux for Cases 1-3. Figure 5: SW parameters
as extracted from the AWSOM solution near the planet during the CME event
(display is similar to Figure 2). In this case, the exoplanet (and the GM
domain) are assumed to be located at an orbital distance from the star along
the star-planet line, facing the center of the CME (similar to an L1 point
satellite near the Earth measuring an incoming CME). The CME arrival time to
the edge of the GM solution occurs around $t=20~{}min$.  Figure 6: Similar
display as in Figure 3 but for the times of the CME event, showing the CME has
pressure-stripped the escaping envelope, with the Ly$\alpha$ absorption
greatly reduced in the $T=24$m and $T=90$m panels. Figure 7: Top: color
contours of the radial velocity displayed on the $x-z$ plane for Cases 1-3.
Velocity three-dimensional streamlines are represented by the black arrow
lines. Bottom: Line plots show the LOS ($u_{x}$) velocity component (solid)
and the thermal pressure (dashed) profiles along the negative $\hat{x}$ (night
side) line starting from the planet up to $20~{}R_{p}$. Figure 8: Profiles of
the different pressure components as a function of radial distance extracted
along the sub-stellar line (positive $\hat{x}$, day side) from the planet up
to $20~{}R_{p}$ (right column), and zoomed up to $3~{}R_{p}$ (left column).
Plots are for Cases 1-3. Solid lines represent the different pressure
components while the dashed line represents the total pressure. Figure 9:
Similar to Figure 8 but showing the profile of the various pressures vs.
radial distance from the planet for the different times during the CME event
close to the planet (left column) and out to $20~{}R_{p}$ (right column). The
pressure increase due to the CME can be seen at $T=24$m at about $8R_{p}$.
Figure 10: Similar display as in Figure 3 but for the neutral hydrogen fluid.


# Partial sensitivity analysis in differential privacy

Tamara T. Müller †‡, Alexander Ziller †‡, Dmitrii Usynin †‡§, Moritz Knolle
†‡§,
Friederike Jungmann †‡, Daniel Rueckert ‡§, Georgios Kaissis ‡§¶
‡Technical University of Munich, §Imperial College London,
†equal contribution<EMAIL_ADDRESS>
###### Abstract

Differential privacy (DP) allows the quantification of privacy loss when the
data of individuals is subjected to algorithmic processing such as machine
learning, as well as the provision of objective privacy guarantees. However,
while techniques such as individual Rényi DP (RDP) allow for granular, per-
person privacy accounting, few works have investigated the impact of each
input feature on the individual’s privacy loss. Here we extend the view of
individual RDP by introducing a new concept we call partial sensitivity, which
leverages symbolic automatic differentiation to determine the influence of
each input feature on the gradient norm of a function. We experimentally
evaluate our approach on queries over private databases, where we obtain a
feature-level contribution of private attributes to the DP guarantee of
individuals. Furthermore, we explore our findings in the context of neural
network training on synthetic data by investigating the partial sensitivity of
input pixels on an image classification task.

##  1 Introduction

The handling and processing of sensitive data inherently entails the risk of
compromising individual privacy and exposing personal information. The
differential privacy (DP) framework dwork2014algorithmic  allows to reason
about privacy loss in the setting of data processing and provide privacy
guarantees to individuals while allowing to draw conclusions from the dataset
as a whole. However, certain attributes of the individual which serve as input
variables, may contribute more to their identifiability in a dataset than
others. Moreover, the individuals themselves may perceive some of their
attributes as being more sensitive and thus more worthy of protection from
exposure. Transparent and trustworthy data processing systems must therefore
be capable of not only accounting for individual privacy loss (e.g. shown in
work on individual Rényi DP (RDP) feldman2020individual ), but to also
selectively apply privacy-preserving mechanisms to specific private attributes
alaggan2015heterogeneous . So far, few works have investigated the
contribution of individual attributes to the query function’s sensitivity, and
therefore to overall privacy loss. This is, in part, because in existent
machine learning frameworks, it is impossible to disentangle the interaction
between input features and model components such as weights. To address this
issue, recent work has introduced symbolic automatic differentiation
(AD)-based sensitivity analysis ziller2021sensitivity . The main contribution
of our work is the introduction of partial sensitivity, representing the
fractional contribution of each input feature to the gradient norm, which
determines privacy loss in DP, and can be derived using symbolic AD. We
demonstrate this technique in the contexts of queries over private databases
and neural network training on synthetic data. Moreover, we showcase the
markedly different influence of non-private stochastic gradient descent (SGD)
on partial sensitivity compared to DP-SGD.

##  2 Related Work

Typically, DP guarantees are allotted homogeneously across individuals in the
dataset, only considering the worst case scenario (i.e. where privacy
accounting is carried out using global sensitivity). As, however, the realised
gradient norm of a query may be considerably lower, the outlook of this type
of privacy accounting may be unnecessarily pessimistic. Individual privacy
accounting, and specifically the work by Feldman et al. feldman2020individual
, proposes an alternative formulation which separately accounts for each
individual’s influence on the outcome of a computation via the actual gradient
norm. The authors express privacy guarantees using Rényi DP mironov2017renyi .
We note that the utilisation of the gradient norm for individual privacy
accounting, which we will be using in this work, is performed in conjunction
with a privacy filter or odometer presented in feldman2020individual .
Heterogeneous DP alaggan2015heterogeneous  captures the variation of privacy
expectations across information held by one individual and uses a modified
Laplace mechanism. Phan et al. Phan2019 ; Phan2017  address the distribution
of noise between different input variables depending on their importance to
the network output using the Laplace or the Heterogeneous Gaussian Mechanism.
Our work is complementary to these techniques as partial sensitivity captures
the influence of individual input attributes on the gradient norm and can
therefore be employed to guide noise addition to selectively perturb them
using one of these heterogeneous noise mechanisms.

##  3 Background

We briefly introduce key terminology, and assume familiarity with the concepts
of differential privacy, the Gaussian mechanism dwork2014algorithmic  and
(DP-)SGD abadi2016deep . We will use the following notations: $D$ and
$D^{\prime}$ for two adjacent datasets, whereby $\sim$ denotes add/remove one
adjacency. We will use $f$ for a function/query and $\Delta_{2}(f)$ for its
$L_{2}$-sensitivity.

######  Definition 1 ($L_{2}$-sensitivity).

The $L_{2}$-sensitivity of a function $f$ is defined as the maximum change in
$L_{2}$-norm in the image of $f$ over all pairs of adjacent inputs:

$\Delta_{2}(f)=\max_{D,D^{\prime}:D\sim D^{\prime}}\|f(D)-f(D^{\prime})\|_{2}$
(1)

The $L_{2}$-sensitivity is used to calibrate the noise addition in the
Gaussian Mechanism. If $f$ is Lipschitz continuous, its Lipschitz constant is
equivalent to its sensitivity raskhodnikova2016lipschitz :

######  Definition 2 (Lipschitz continuity).

Let $f$ be a function $X\rightarrow Y$ with associated metrics $d_{X}$ and
$d_{Y}$. $f$ is Lipschitz continuous if and only if there exists a real
constant $K\geq 0$ such that, $\forall x_{1},x_{2}\in X$:

$d_{Y}(f(x_{1}),f(x_{2}))\leq Kd_{X}(x_{1},x_{2})$ (2)

The Lipschitz constant of a real-valued function is equal to the supremum of
the norm of its gradient: $K(f)=\sup_{x}\|\nabla f(x)\|_{2}$.

##  4 Theoretical Results

######  Definition 3 (Partial sensitivity of a function).

Let $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{1},m>1$ be a real-valued function.
The partial sensitivity of $f$ at an input $X$ is defined as the gradient of
the gradient norm of $f$ with respect to the input variables:

$\Delta_{p}(f(X)):=\nabla_{X}\|\nabla f(X)\|_{2}$ (3)

The definition of partial sensitivity is derived from the following property:
Let $\mathbf{x}=\left(x_{1},x_{2},\dots,x_{d}\right)^{T}$ be the vector of
inputs. Since $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{1}$, $f(\mathbf{x})$ is
a scalar. As such, $\nabla f(\mathbf{x})=\left(\frac{\partial
f(\mathbf{x})}{\partial x_{1}},\dots,\frac{\partial f(\mathbf{x})}{\partial
x_{d}}\right)$. Furthermore, let $\left(\frac{\partial f(\mathbf{x})}{\partial
x_{i}}\right)=a_{i},i\leq d$ and let $\varphi$ be the $L_{2}$-norm of
$\nabla{f(\mathbf{x})}$, given by
$\varphi=\|\nabla{f(\mathbf{x})}\|_{2}=\sqrt{\sum_{i=1}^{d}a_{i}^{2}}$ which
is also a scalar. As we are interested in the rate of change of the
sensitivity with respect to the individual input attributes,
$(x_{1},\dots,x_{d})$, we obtain the gradient of $\varphi$:

$\displaystyle\nabla\varphi=\left(\frac{\partial}{\partial
x_{1}}\sqrt{\sum_{i=1}^{d}a_{i}^{2}},\dots,\frac{\partial}{\partial
x_{d}}\sqrt{\sum_{i=1}^{d}a_{i}^{2}}\right)$ (4)

The $j_{th}$ component of $\nabla\varphi$, $\nabla^{j}\varphi$, where $j\leq
d$ is given by:

$\displaystyle\frac{\partial}{\partial
x_{j}}\sqrt{\sum_{i=1}^{d}a_{i}^{2}}=\frac{1}{2}\frac{2a_{j}}{\sqrt{\sum_{i=1}^{d}a_{i}^{2}}}=\frac{a_{j}}{\varphi}=\frac{a_{j}}{\|\nabla{f(\mathbf{x})}\|_{2}}=\frac{\frac{\partial
f(\mathbf{x})}{\partial x_{j}}}{\|\nabla{f(\mathbf{x})}\|_{2}}$ (5)

Thus:

$\displaystyle\Delta^{(x_{1},\dots,x_{d})}_{p}(f(\mathbf{x}))=\left(\frac{\frac{\partial
f(\mathbf{x})}{\partial
x_{1}}}{\|\nabla{f(\mathbf{x})}\|_{2}},\dots,\frac{\frac{\partial
f(\mathbf{x})}{\partial x_{d}}}{\|\nabla{f(\mathbf{x})}\|_{2}}\right)$ (6)

Partial sensitivity therefore represents the fractional contribution of the
individual input attributes to the gradient norm of the function. Its symbolic
representation, which can be easily obtained using a symbolic AD system
ziller2021sensitivity , is independent of the actual input data and can
therefore be used to interpret the impact of individual input attributes on
the gradient norm, used for (individual) privacy accounting, as shown in the
next section.

##  5 Experiments

###  5.1 Partial sensitivity analysis of multivariate database queries

We assume a scenario where an analyst wants to construct a statistical query
$f(a,b)$ over a database containing private data from a population. They
specify reasonable ranges on the private attributes
$a\in[1.0,2.0],b\in[0.5,3]$ (e.g. based on prior knowledge but not on the
actual data values). Let $f(a,b)=a^{2}+e^{2b-a}$ and
$\phi:X\rightarrow\mathbb{R}^{1}$ be a linear aggregation function with
sensitivity $1$. To obtain the Rényi DP guarantee for one individual $i$ as
specified in Section 2, the Lipschitz constant $(L_{\phi})_{i}$ of $\phi$ and
the norm of the gradient of $f$ of each individual $i$, $\|\nabla_{i}f\|_{2}$,
are required. Following Definition 2.7 from feldman2020individual  and
Definition 2, the linearity of $\phi$ means that the same Lipschitz constant
$L_{\phi}=\sup\|\nabla\phi\|_{2}=1$ can be used for all individuals. The
privacy loss of individual $X_{i}$ can then be calculated as follows:
$(\alpha,\frac{\alpha L_{\phi}^{2}\|\nabla
f_{i}\|_{2}^{2}}{2\sigma^{2}})=(\alpha,\frac{\alpha\|\nabla
f_{i}\|_{2}^{2}}{2\sigma^{2}})$, where $\alpha$ is the Rényi divergence order
and $\sigma^{2}$ the variance of the Gaussian noise added by the mechanism.
Furthermore, compiling the symbolic expression for the gradient norm of $f$
and maximising its value using a suitable technique such as simplicial
homology global optimisation endres2018simplicial , the global sensitivity
$\Delta_{2}(f)$ can be calculated given the pre-defined ranges of the input
variables $a$ and $b$. This has to be done data-independently, and
$\Delta_{2}(f)$ is a constant (in this example $\Delta_{2}(f)=66.19$).
Finally, the partial sensitivity of $f$, $\Delta^{(a,b)}_{p}(f)$, can be
analysed and plotted to permit visual reasoning about the effects of $a$ and
$b$ on $\|\nabla_{i}f\|_{2}$, and thus on Rényi DP guarantees. Figure 1 shows
how the two variables contribute to the sensitivity of $f$. This knowledge
allows to assign a different weight to each attribute and calibrate noise
addition to selectively protect specific attributes.

Figure 1: Partial sensitivity plot for section 5.1. The horizontal axes show
the partial sensitivity of $f$ with respect to inputs $a$ and $b$. The
vertical axis shows $\|\nabla f_{i}\|_{2}$. $\Delta_{2}(f)$ given the input
ranges is marked in red.

###  5.2 Partial sensitivity modelling in (DP-)SGD

In this section, we investigate the impact of individual input attributes on
the privacy loss of a specific neural network training example. To illustrate
this effect, we performed binary classification on synthetic $5\times 5$-pixel
images of vertical and horizontal bars using SGD and DP-SGD using a two-layer
neural network. The symbolic representation for the partial sensitivity of the
loss function, $\Delta_{p}^{(x_{1,1},\dots,x_{5,5})}(\mathcal{L})$ can again
be derived through symbolic AD. By substitution of the weights, class label
and input pixel values, it is possible to obtain the partial sensitivity of
the input pixels at any moment in training. At the end of training, where we
consider the weights to be fixed, this can be used to reason over the relevant
input features leading the network to assign an image to one of the two
classes. For the vertical bar class, the maximum partial sensitivity values
were observed across the horizontal of the image, and vice versa for the
horizontal bar class (Figure 2A). Evidently, the presence of features
characteristic for the opposite class renders the individual an outlier in the
specific distribution, resulting in a large change in the function’s gradient
norm and thus, a high individual privacy loss. Interestingly, training with
DP-SGD led to a substantial reduction in the partial sensitivities of the
corresponding pixels, while increasing the values of unrelated pixels in the
input space. To further illustrate this effect, we generated histograms of the
partial sensitivities for each pixel over one thousand samples. We found the
distributions of the partial sensitivities to be highly concentrated around
specific values in the case of SGD, while being substantially more dispersed
and centered around zero for DP-SGD (Figure 2B). From this we conjecture that
the effect of DP-SGD training is the homogenisation of the partial sensitivity
across the space of inputs, which may suppress the memorisation of strongly
identifying private attributes. 111Source code at:
https://github.com/tamaramueller/Deuterium_Partial_Sensitivity

Figure 2: Pixel-wise privacy loss analysis of neural networks. A shows input
images of the vertical and horizontal bar class, the maximum partial
sensitivities for each pixel for SGD and DP-SGD. B shows histograms of the
partial sensitivity of each pixel using SGD and DP-SGD.

##  6 Discussion and Conclusions

Our work extends sensitivity analysis to the level of individual input
attributes by introducing partial sensitivity. Individuals may wish to weight
their specific attributes in a personalised manner as they consider them to
have a different privacy impact, or even allocate fractions of their privacy
budget to different analysts based on such feature-level granularity. Analysts
on the other hand, may be interested to attenuate the specific features
leading to high privacy expenditure, as demonstrated in our case study on DP-
SGD. Our results may help to better understand the privacy-utility trade-offs
of this algorithm dp_impact_accuracy ; dp_bad_2 ; avent2019automatic , its
impact on fairness fioretto2021decision , and the link between memorisation of
private features, model overfitting, and generalisation feldman2020does .
Moreover, partial sensitivity exhibits strong links to gradient-based model
introspection techniques gradcam . We intend to explore these topics, and the
utilisation of partial sensitivity to guide targeted noise addition via
heterogeneous mechanisms, in future work.

## References

  * [1] Cynthia Dwork, Aaron Roth, et al.  The algorithmic foundations of differential privacy.  Foundations and Trends in Theoretical Computer Science, 9(3-4):211–407, 2014.
  * [2] Vitaly Feldman and Tijana Zrnic.  Individual privacy accounting via a renyi filter.  arXiv preprint arXiv:2008.11193, 2020.
  * [3] Mohammad Alaggan, Sébastien Gambs, and Anne-Marie Kermarrec.  Heterogeneous differential privacy.  arXiv preprint arXiv:1504.06998, 2015.
  * [4] Alexander Ziller, Dmitrii Usynin, Moritz Knolle, Kritika Prakash, Andrew Trask, Rickmer Braren, Marcus Makowski, Daniel Rueckert, and Georgios Kaissis.  Sensitivity analysis in differentially private machine learning using hybrid automatic differentiation.  arXiv preprint arXiv:2107.04265, 2021.
  * [5] Ilya Mironov.  Rényi differential privacy.  In 2017 IEEE 30th Computer Security Foundations Symposium (CSF), pages 263–275. IEEE, 2017.
  * [6] NhatHai Phan, Minh Vu, Yang Liu, Ruoming Jin, Dejing Dou, Xintao Wu, and My T Thai.  Heterogeneous gaussian mechanism: Preserving differential privacy in deep learning with provable robustness.  arXiv preprint arXiv:1906.01444, 2019.
  * [7] NhatHai Phan, Xintao Wu, Han Hu, and Dejing Dou.  Adaptive laplace mechanism: Differential privacy preservation in deep learning.  In 2017 IEEE International Conference on Data Mining (ICDM), pages 385–394. IEEE, 2017.
  * [8] Martin Abadi, Andy Chu, Ian Goodfellow, H Brendan McMahan, Ilya Mironov, Kunal Talwar, and Li Zhang.  Deep learning with differential privacy.  In Proceedings of the 2016 ACM SIGSAC conference on computer and communications security, pages 308–318, 2016.
  * [9] Sofya Raskhodnikova and Adam Smith.  Lipschitz extensions for node-private graph statistics and the generalized exponential mechanism.  In 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), pages 495–504. IEEE, 2016.
  * [10] Stefan C Endres, Carl Sandrock, and Walter W Focke.  A simplicial homology algorithm for Lipschitz optimisation.  Journal of Global Optimization, 72(2):181–217, 2018.
  * [11] Eugene Bagdasaryan and Vitaly Shmatikov.  Differential Privacy Has Disparate Impact on Model Accuracy.  arXiv preprint arXiv:1905.12101, 2019.
  * [12] Peter Kairouz, Keith Bonawitz, and Daniel Ramage.  Discrete distribution estimation under local privacy.  arXiv preprint arXiv:1602.07387, 2016.
  * [13] Brendan Avent, Javier Gonzalez, Tom Diethe, Andrei Paleyes, and Borja Balle.  Automatic discovery of privacy-utility pareto fronts.  arXiv preprint arXiv:1905.10862, 2019.
  * [14] Ferdinando Fioretto, Cuong Tran, and Pascal Van Hentenryck.  Decision making with differential privacy under a fairness lens.  arXiv preprint arXiv:2105.07513, 2021.
  * [15] Vitaly Feldman.  Does learning require memorization? a short tale about a long tail.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing, pages 954–959, 2020.
  * [16] Ramprasaath R. Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  International Journal of Computer Vision, 128(2):336–359, Oct 2019\.

##  Appendix A Appendix

#### Description of computational resources used

All experimentation was conducted on a single workstation computer with 18
physical CPU cores and 512 GB of physical memory, as well as 2 NVidia Quadro
RTX 8000 GPUs running Ubuntu Linux v. 20.04 LTS.

#### Model architectures

For the partial sensitivity modelling in DP-SGD, we designed a feedforward
neural network consisting of an input layer of dimensionality $25\times 8$, an
intermediate layer of dimensionality $8\times 8$ with a bias vector of length
$8$ and an output layer of dimensionality $8\times 1$. All layers were
followed by logistic sigmoid activation functions. Training was performed
using the binary cross-entropy loss function.

#### Synthetic data generation

Synthetic data generation was conducted as follows:
Two base images with white background were constructed, where one consisted of
a vertical bar of black pixels and the other of a horizontal line. From these,
$1000$ images per class were created by the addition of random Gaussian noise
with mean zero and a standard deviation of $0.2$. Evaluation was performed on
a test set of $100$ images per class. A fixed random seed was used to
deterministically generate the images in the experiments of the main
manuscript.

#### Hyperparameters

We trained all networks to convergence using SGD with a batch size of $2000$
(batch gradient descent) and a learning rate $\eta=0.1$. We repeated this
procedure using DP-SGD with the same learning rate, a noise multiplier of
$5.0$ and an $L_{2}$-bound of $0.1$.

#### Details on the attached code

Symbolic automatic differentiation was performed using the Deuterium framework
[4], whose source code is available alongside the experiment code. Compilation
performed by Deuterium relies on a suitable compiler for the C or Fortran
programming language or for the LLVM tool-chain.


derived from the track fitting algorithm. In this paper, the track fitting was
done with the custom software implementation of the General Broken Lines
algorithm [20].

Minimization of the $\chi^{2}$ of Eq.(47) leads the partial derivative with
respect to each ($l$-th) local track parameter of each ($i$-th) track,
$\Delta{q}_{il}$, to equal zero:

$\begin{split}\frac{\partial{\chi}^{2}}{\partial
q_{il}}=&2\sum_{j=1}^{n_{meas}}\Bigl{(}{\frac{\partial\boldsymbol{\varepsilon}_{ij}}{\partial
q_{il}}}\Bigr{)}^{\mathsf{T}}\mathbf{V}_{ij}^{-1}\boldsymbol{\varepsilon}_{ij}+2\sum_{j=2}^{n_{scat}-1}\Bigl{(}{\frac{\partial\boldsymbol{\beta}_{ij}}{\partial
q_{il}}}\Bigr{)}^{\mathsf{T}}\mathbf{W}_{ij}^{-1}\boldsymbol{\beta}_{ij}=0\\\
{\simeq}&2\sum_{j=1}^{n_{meas}}\Bigl{(}{\frac{\partial\boldsymbol{\varepsilon}_{ij}}{\partial
q_{il}}}\Bigr{)}^{\mathsf{T}}\mathbf{V}_{ij}^{-1}\Bigl{[}\boldsymbol{\varepsilon}_{ij}(\boldsymbol{q}_{i}^{0},\boldsymbol{p}^{0})+\sum_{l^{\prime}}{\frac{\partial\boldsymbol{\varepsilon}_{ij}}{\partial
q_{il^{\prime}}}}{\Delta}q_{il^{\prime}}+\sum_{g^{\prime}}{\frac{\partial\boldsymbol{\varepsilon}_{ij}}{\partial
p_{g^{\prime}}}}{\Delta}p_{g^{\prime}}\Bigr{]}\\\
&+2\sum_{j=2}^{n_{scat}-1}\Bigl{(}{\frac{\partial\boldsymbol{\beta}_{ij}}{\partial
q_{il}}}\Bigr{)}^{\mathsf{T}}\mathbf{W}_{ij}^{-1}\sum_{l^{\prime}}{\frac{\partial\boldsymbol{\beta}_{ij}}{\partial
q_{il^{\prime}}}}{\Delta}q_{il^{\prime}}\end{split}$ (B.6)

Eq.(B.6) can be simplified in matrix form as:

$\boldsymbol{b}^{i}=(\mathbf{G}^{i})^{\mathsf{T}}{\Delta}\boldsymbol{p}+\mathsf{\Gamma}^{i}{\Delta}\boldsymbol{q}_{i}$
(B.7)

where $\boldsymbol{b}^{i}$ is a vector whose $l$-th element is given by:

$b_{l}^{i}=-\sum_{j=1}^{n_{meas}}\Bigl{(}{\frac{\partial\boldsymbol{\varepsilon}_{ij}}{\partial
q_{il}}}\Bigr{)}^{\mathsf{T}}\mathbf{V}_{ij}^{-1}\boldsymbol{\varepsilon}_{ij}(\boldsymbol{q}_{i}^{0},\boldsymbol{p}^{0})$
(B.8)

$(\mathbf{G}^{i})^{\mathsf{T}}$ is the transpose of the matrix
$\mathbf{G}^{i}$ which is defined in Eq.(B.5), and $\mathsf{\Gamma}^{i}$ is a
matrix whose $(l,l^{\prime})$ entry is given by:

$\Gamma_{ll^{\prime}}^{i}=\sum_{j=1}^{n_{meas}}\Bigl{(}{\frac{\partial\boldsymbol{\varepsilon}_{ij}}{\partial
q_{il}}}\Bigr{)}^{\mathsf{T}}\mathbf{V}_{ij}^{-1}\frac{\partial\boldsymbol{\varepsilon}_{ij}}{\partial
q_{il^{\prime}}}+\sum_{j=2}^{n_{scat}-1}\Bigl{(}{\frac{\partial\boldsymbol{\beta}_{ij}}{\partial
q_{il}}}\Bigr{)}^{\mathsf{T}}\mathbf{W}_{ij}^{-1}{\frac{\partial\boldsymbol{\beta}_{ij}}{\partial
q_{il^{\prime}}}}$ (B.9)

The partial derivatives of the scattering angle with respect to the local
track parameters,
${\partial\boldsymbol{\beta}_{ij}}/{\partial\boldsymbol{q}_{i}}$, are derived
from the track fitting algorithm.

Combining Eq.(B.2) and Eq.(B.7), all the global alignment parameters,
$\Delta\boldsymbol{p}$, and all the local track parameters,
$\Delta\boldsymbol{q}$, can be solved simultaneously from following matrix
equation:

$\begin{pmatrix}\sum_{i}\mathbf{C}^{i}&\mathbf{G}^{1}&\ldots&\mathbf{G}^{j}&\ldots&\mathbf{G}^{N}\\\
(\mathbf{G}^{1})^{\mathsf{T}}&\mathsf{\Gamma}^{1}&\ldots&\mathbf{0}&\ldots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
(\mathbf{G}^{j})^{\mathsf{T}}&\mathbf{0}&\ldots&\mathsf{\Gamma}^{j}&\ldots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
(\mathbf{G}^{N})^{\mathsf{T}}&\mathbf{0}&\ldots&\mathbf{0}&\ldots&\mathsf{\Gamma}^{N}\\\
\end{pmatrix}\begin{pmatrix}\Delta\boldsymbol{p}\\\
\Delta\boldsymbol{q}_{1}\\\ \vdots\\\ \Delta\boldsymbol{q}_{j}\\\ \vdots\\\
\Delta\boldsymbol{q}_{N}\\\ \end{pmatrix}=\\\
\begin{pmatrix}\sum_{i}\boldsymbol{d}^{i}\\\ \boldsymbol{b}^{1}\\\ \vdots\\\
\boldsymbol{b}^{j}\\\ \vdots\\\ \boldsymbol{b}^{N}\\\ \end{pmatrix}$

## References

  * [1] V. Bindi, et al., Nuclear Instruments and Methods in Physics Research Section A 743, 22 (2014).  https://doi.org/10.1016/j.nima.2014.01.002
  * [2] T. Kirn, et al., Nuclear Instruments and Methods in Physics Research Section A 706, 43 (2013).  https://doi.org/10.1016/j.nima.2012.05.010
  * [3] M. Aguilar, et al., Nuclear Instruments and Methods in Physics Research Section A 614, 237 (2010).  https://doi.org/10.1016/j.nima.2009.12.027
  * [4] C. Adloff, et al., Nuclear Instruments and Methods in Physics Research Section A 714, 147 (2013).  https://doi.org/10.1016/j.nima.2013.02.020
  * [5] P. von Doetinchem, et al., Nuclear Physics B-Proceedings Supplements 197, 15 (2009).  https://doi.org/10.1016/j.nuclphysbps.2009.10.025
  * [6] K. Luebelsmeyer, et al., Nuclear Instruments and Methods in Physics Research Section A 654(1), 639 (2011).  https://doi.org/10.1016/j.nima.2011.06.051
  * [7] M. Aguilar, et al., Physical Review Letters 114(17), 171103 (2015).  https://doi.org/10.1103/PhysRevLett.114.171103
  * [8] W. Burger, et al., Nuclear Instruments and Methods in Physics Research Section A 512(3), 517 (2003).  https://doi.org/10.1016/S0168-9002(03)01866-7
  * [9] G. Ambrosi., V. Choutko, C. Delgado, A. Oliva, Q. Yan, and Y. Li, Nuclear Instruments and Methods in Physics Research Section A 869, 29 (2017).  https://doi.org/10.1016/j.nima.2017.07.014
  * [10] T. Siedenburg, et al., Nuclear Physics B-Proceedings Supplements 113(1-3), 154 (2002).  https://doi.org/10.1016/S0920-5632(02)01835-2
  * [11] S. Borsini, et al., SAE Technical Paper pp. 01–3023 (2007).  https://doi.org/10.4271/2007-01-3023
  * [12] Q. Yan, V. Choutko, A. Oliva, and M. Paniccia, Nuclear Physics A 996, 121712 (2020).  https://doi.org/10.1016/j.nuclphysa.2020.121712
  * [13] M. Stoye, Calibration and Alignment of the CMS Silicon Tracking Detector.  Ph.D. thesis, U. Hamburg, Germany (2007).  URL https://cds.cern.ch/record/1047047
  * [14] C. Thiel, J. Voss, R.J. Martin, and P. Neff, International Journal of Non-Linear Mechanics 112, 57 (2019).  https://doi.org/10.1016/j.ijnonlinmec.2018.10.002
  * [15] V. Blobel, Nuclear Instruments and Methods in Physics Research Section A 566(1), 5 (2006).  https://doi.org/10.1016/j.nima.2006.05.157
  * [16] G. Flucke, P. Schleper, G. Steinbrück, and M. Stoye, Journal of Instrumentation 3(09), P09002 (2008).  https://doi.org/10.1088/1748-0221/3/09/P09002
  * [17] S. Viret, C. Parkes, and M. Gersabeck, Nuclear Instruments and Methods in Physics Research Section A 596(2), 157 (2008).  https://doi.org/10.1016/j.nima.2008.07.153
  * [18] ALICE collaboration, Journal of Instrumentation 5(03), P03003 (2010).  https://doi.org/10.1088/1748-0221/5/03/P03003
  * [19] ATLAS Collaboration, European Physical Journal C 80, 1194 (2020).  https://doi.org/10.1140/epjc/s10052-020-08700-6
  * [20] V. Blobel, C. Kleinwort, and F. Meier, Computer Physics Communications 182(9), 1760 (2011).  https://doi.org/10.1016/j.cpc.2011.03.017
  * [21] G.R. Lynch and O.I. Dahl, Nuclear Instruments and Methods in Physics Research Section B 58(1), 6 (1991).  https://doi.org/10.1016/0168-583X(91)95671-Y
  * [22] C. Kleinwort, Nuclear Instruments and Methods in Physics Research Section A 673, 107 (2012).  https://doi.org/10.1016/j.nima.2012.01.024
  * [23] V. Blobel and C. Kleinwort, arXiv preprint hep-ex/0208021 (2002).  https://doi.org/10.48550/arXiv.hep-ex/0208021
  * [24] V. Blobel, Institut für Experimentalphysik, Universität Hamburg (2007).  URL https://www.desy.de/~sschmitt/blobel/Mptwo.pdf
  * [25] DESY gitlab.  Millepede-II.  URL https://gitlab.desy.de/claus.kleinwort/millepede-ii.git
  * [26] D. Leonardo and R. Menon, Computing in Science and Engineering 5(1), 46 (1998).  https://doi.org/10.1109/99.660313
  * [27] R. Brun and F. Rademaker, Nuclear Instruments and Methods in Physics Research Section A 389(1-2), 81 (1997).  https://doi.org/10.1016/S0168-9002(97)00048-X
  * [28] A.J. Peters, E.A. Sindrilaru, and G. Adde, in _Journal of Physics: Conference Series_ , vol. 664 (2015), vol. 664, p. 042042.  https://doi.org/10.1088/1742-6596/664/4/042042
  * [29] S. Agostinelli, et al., Nuclear Instruments and Methods in Physics Research Section A 506(3), 250 (2003).  https://doi.org/10.1016/S0168-9002(03)01368-8
  * [30] M. Aguilar, et al., Physical Review Letters 115(21), 211101 (2015).  https://doi.org/10.1103/PhysRevLett.115.211101
  * [31] C.H. Reinsch, Numerische Mathematik 10(3), 177 (1967).  https://doi.org/10.1007/BF02162161
  * [32] D.F. Smart and M.A. Shea, Advances in Space Research 36(10), 2012 (2005).  https://doi.org/10.1016/j.asr.2004.09.015
  * [33] E. Thébault, et al., Earth, Planets and Space 67(1), 1 (2015).  https://doi.org/10.1186/s40623-015-0228-9
  * [34] J. Berdugo, V. Choutko, C. Delgado, and Q. Yan, Nuclear Instruments and Methods in Physics Research Section A 869, 10 (2017).  https://doi.org/10.1016/j.nima.2017.07.012
  * [35] M. Aguilar, et al., Physics Reports 894, 1 (2021).  https://doi.org/10.1016/j.physrep.2020.09.003
  * [36] CMS collaboration, Journal of Instrumentation 9(06), P06009 (2014).  https://doi.org/10.1088/1748-0221/9/06/P06009