11institutetext: Anton Pannekoek Institute for Astronomy, University of
Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

# Expanding the inventory of spectral lines used to trace atmospheric escape
in exoplanets

D.C. Linssen E-mail<EMAIL_ADDRESS>A. Oklopčić

(Received 3 April 2023 / Accepted 12 June 2023)

Escaping exoplanet atmospheres have been observed as deep transit signatures
in a few specific spectral lines. Detections have been made in the hydrogen
Ly-$\alpha$ line, the metastable helium line at 10830 Å and some UV lines of
metallic species. Observational challenges, unexpected non-detections and
model degeneracies have generally made it difficult to draw definitive
conclusions about the escape process for individual planets. Expanding on the
suite of spectral tracers used may help to mitigate these challenges. We
present a new framework for modeling the transmission spectrum of
hydrodynamically escaping atmospheres. We predict FUV to NIR spectra for
systems with different planet and stellar types and identify new lines that
can potentially be used to study their upper atmospheres. Measuring the radius
in the atmosphere at which the strongest lines form puts them into context
within the upper atmospheric structure. Targeting a set of complementary
spectral lines for the same planet will help us to better constrain the
outflow properties.

##  1 Introduction

Planets that orbit close to their host stars can experience extreme levels of
high-energy stellar irradiation, which can drive hydrodynamic escape of their
atmospheres (see review by Owen, 2019). At the same time, heat released by the
cooling core of a planet can also supply energy to drive the escape process
(Ginzburg et al., 2018). Hot Jupiters undergoing escape have been shown to be
stable over their lifetimes (e.g. Murray-Clay et al., 2009; Vissapragada et
al., 2022). Less massive planets such as hot Neptunes might however lose a
significant fraction of their atmosphere in this way (Ionov et al., 2018),
which can affect their composition and stability. This process has been
invoked to explain the bimodal radius distribution of sub-Neptune planets
around $\sim$1.6 $R_{\oplus}$ (e.g. Owen & Wu, 2017; Ginzburg et al., 2018),
as well as the lack of sub-Jovian planets at short orbital distances (e.g.
Owen & Wu, 2013; Lundkvist et al., 2016).

Atmospheric escape has been observed for a growing sample of exoplanets in a
few specific spectral lines that trace the low-density gas in the outflow.
Absorption in the wings of the Ly-$\alpha$ line has provided direct evidence
of gas escaping beyond the planet’s Hill radius (e.g. Vidal-Madjar et al.,
2003). This line allows us to study the regions where the planetary wind
interacts with the stellar environment, but the wind-launching regions are
obscured (Owen et al., 2022). The metastable helium triplet can be observed at
high resolution in the near-infrared (e.g. Nortmann et al., 2018; Allart et
al., 2018). This line probes deeper into the planetary outflow than
Ly-$\alpha$ and therefore does allow us to study the wind-launching mechanism
(Oklopčić & Hirata, 2018). For a small sample of exoplanets, strong absorption
from metal species in the UV has also indicated the presence of a planetary
outflow. This includes the Fe II UV1 and UV2 multiplets (Sing et al., 2019;
Cubillos et al., 2020), the Mg II doublet (Fossati et al., 2010; Sing et al.,
2019), the C II $\lambda$ 1335 lines (Vidal-Madjar et al., 2004; García Muñoz
et al., 2021; Ben-Jaffel et al., 2022), the O I $\lambda$ 1305 lines (Vidal-
Madjar et al., 2004; Ben-Jaffel & Ballester, 2013) and the Si III $\lambda$
1206 line (Linsky et al., 2010).

Each of these spectral lines have their unique observational challenges. For
example, the ISM and geocorona render the core of the Ly-$\alpha$ line
inaccessible (e.g. Lecavelier des Etangs et al., 2010), and the helium line
preferentially forms in planets around K-type stars (Oklopčić, 2019). Still,
there are many exoplanets for which non-detections of atmospheric-escape
tracers have been reported, even though atmospheric mass loss is expected
based on e.g. the received irradiation (see dos Santos 2022 for an overview).
Similarly for some planets, ongoing atmospheric mass loss has been observed in
one or more spectral tracers, but has eluded detection in others. For example,
in HD 189733 b, Ly-$\alpha$ (Lecavelier des Etangs et al., 2010), metastable
helium (Salz et al., 2018) and O I (Ben-Jaffel & Ballester, 2013) have been
detected, but Mg II absorption has not (Cubillos et al., 2023). The
interpretation of (non-)detected signatures has also not been straightforward.
For example, different types of models have been used to interpret helium
observations (e.g. a particle code as employed by Allart et al. (2018), a 3D
hydrodynamics code as employed by Wang & Dai (2021) and a 1D isothermal Parker
wind model as employed by Mansfield et al. (2018)), with varying predictions
and inferences about some of the core properties of the outflow, such as the
mass-loss rate.

Finding new spectral lines that probe escaping atmospheres may mitigate these
challenges and help to enhance our understanding of atmospheric mass loss. A
step in this direction was taken by Turner et al. (2016), who predicted
spectral signatures for a uniform slab of material at conditions similar to
that of escaping planetary gas. In this work, we present a new framework for
producing synthetic transmission spectra of evaporating exoplanet atmospheres
from FUV to NIR wavelengths and we use it to identify new spectral lines that
can potentially be used to study atmospheric escape. We measure the
atmospheric radii probed by the most prominent spectral lines in order to
create a more complete picture of the upper atmosphere.

##  2 Methods

###  2.1 Modeling the escaping atmosphere with Cloudy

We model the escaping planet atmosphere as a Parker wind (Parker, 1958; Lamers
& Cassinelli, 1999) made with the p_winds code (dos Santos et al., 2022). This
is a parametrized outflow model with an isothermal temperature and mass-loss
rate as parameters, as well as a chemical composition that we keep fixed at
90% hydrogen and 10% helium by number. We employ our iterative algorithm
presented in Linssen et al. (2022) to simulate the Parker wind density profile
with Cloudy v17.02111https://gitlab.nublado.org/cloudy/cloudy/-/wikis/home
(Ferland et al., 1998, 2017). In short, our algorithm solves for a (non-
isothermal) temperature structure based on the radiative heating and cooling
rates as calculated by Cloudy, while also including the hydrodynamic thermal
effects of expansion cooling and heat advection. Since the mass-loss rate and
constant temperature are free parameters of the Parker wind model, many
different models can be constructed for a given planet. The chosen temperature
may however be unrealistic based on the stellar flux and density and velocity
of the outflow. We leverage this to identify ‘self-consistent’ Parker wind
profiles through a comparison between the chosen isothermal temperature and
Cloudy’s converged temperature structure. There is freedom in how to
practically do this comparison and the most useful way depends on the aim of
the modeling effort. In Linssen et al. (2022) we interpreted metastable helium
detections and hence constructed a characteristic temperature in the helium
line-forming region from the Cloudy simulations, which we compared to the
original isothermal temperature. In this work we aim to simulate a whole suite
of spectral lines probing different parts of the upper atmosphere. It is thus
impossible to find one characteristic temperature for the line-forming regions
of all these lines, and hence we choose to simply calculate the radial average
of Cloudy’s temperature profile between 1 and 2 planetary radii and compare it
to the isothermal value. At those radii the planetary wind is launched and
respecting self-consistency in this region generally results in a selection of
Parker wind models whose density and velocity profiles are most similar to
self-consistent outflow models (see also Sec. 4.4). In this way we focus on
obtaining a sensible overall outflow structure instead of focusing on one part
of the outflow where a particular spectral line forms. Still, we reiterate
that the choice is ultimately arbitrary. After checking the self-consistency
of the models, for one value of the mass-loss rate we tend to find one Parker
wind model with a self-consistent temperature and consider just this model. We
thus effectively reduce the free parameters of the Parker wind model to only
the mass-loss rate (and the composition). In Sec. 3 we explain which specific
Parker wind model parameters we use in this work.

Once we have found a self-consistent Parker wind model, we calculate its
transmission spectrum. We do not use Cloudy’s own predicted spectrum for this,
but take its thermal and chemical atmospheric profiles and post-process these
to calculate the spectral lines found in the NIST
database222https://www.nist.gov/pml/atomic-spectra-database. Appendix B
describes the procedure in more detail.

###  2.2 Selecting suitable spectral lines

We calculate mid-transit transmission spectra over the wavelength range
911-20,000 Å (see Fig. 2 for an example spectrum). Below 911 Å, the continuum
opacity of hydrogen (and helium) ionization renders other spectral lines
inaccessible. Above 20,000 Å, there are almost no strong atomic spectral
lines. In each transmission spectrum, we then search for spectral lines that
are strong enough to be potentially observable. For every atomic and ionic
species, we find the strongest spectral line in the UV (¡4000 Å) and the
strongest line in the VIS/NIR (¿4000 Å). If the strongest UV line has more
than 1% excess absorption at line center, we include it in our sample of
suitable tracers, and then do the same for the strongest VIS/NIR spectral line
if it is stronger than 0.1%. We chose to split the wavelength domain in two
and use different thresholds because the UV and VIS/NIR are typically observed
by different instruments with different precision. For some species such as
e.g. Fe+, there are many strong spectral lines that could be used as upper
atmospheric tracers (and it may be possible to use cross-correlation
techniques), but we include here only the strongest spectral line to keep the
analysis clear and straightforward. We emphasize that we are modeling
hydrodynamic upper atmospheres that are completely atomic/ionic, and thus any
molecular absorption bands originating from the ‘lower’ atmosphere are absent
from our spectra.

###  2.3 Calculating the radius where spectral lines form

From looking at the transmission spectrum, one might naively assume that
deeper spectral lines originate at higher altitudes in the atmosphere. While
this is true in many cases, it need not always be, as a spectral line of a
given depth might originate from a confined region close to the planet with
relatively high optical depth or from a more diffuse, extended region with
relatively low optical depth, thereby probing a different part of the outflow.
In this work, we want to distinguish between these cases and measure the
radius of the atmospheric region in which each spectral line forms.

Figure 1: Illustration of how the line-forming radius is defined for the Mg II
line at 2796 Å. We consecutively keep adding deeper atmospheric layers and
track the increase in equivalent width (EW) of the spectral line. Top panel:
2D projection of the atmosphere separated into different layers (indicated by
different colors). We show 18 layers for clarity, but our real calculations
consist of 60 layers. The planet is indicated by the black circle. The star is
to the left and the observer to the right. Middle panel: the ‘observed’
spectral line when only including the part of the atmosphere with the
corresponding color as well as lighter shaded layers (i.e. more yellow
layers). Bottom panel: contributions to the EW of the line from each
atmospheric layer. The black point indicates the line-forming radius and
standard deviations that are calculated based on the histogram and defined in
Eqs. 1 and 2 (the plotted $\Delta$EW value of the black point is meaningless).

To this end, we devise a measure for the ‘line-formation radius’ of a given
spectral line, based on the contribution to the equivalent width (EW) of the
line of each atmospheric layer (illustrated in Fig. 1). We divide our
atmosphere into a dayside ($x<0$) and nightside ($x>0)$ and further into
layers logarithmically spaced in radius. We then start by calculating the
transit spectrum and its EW while only including the material of the uppermost
layer on the dayside (the yellow layer between 8$R_{p}$ and 6.3$R_{p}$ in the
figure). We then add a deeper layer to the atmosphere (the orange layer
between 6.3$R_{p}$ and 5$R_{p}$ in the figure) and re-calculate the transit
spectrum and its EW. The increase in EW is the contribution of the second
layer, and these contributions per layer are shown in the bottom panel of Fig.
1. We repeat this process until all dayside layers are included, after which
we add the lowest nightside layer (the pink layer between $1R_{p}$ and
$\sim$1.3$R_{p}$ in the figure) and work our way outwards again. We finally
co-add the EW contributions of the dayside and nightside layers of equal
radius, which results in a histogram like the bottom panel of Fig. 1.

The reason we treat the dayside and nightside separately is that this is also
the order in which a light ray originating from the star would probe the
planetary atmosphere. This is important when the dayside layer at, for
example, $8R_{p}$ absorbs the first light and thus contributes to the EW, but
then the atmosphere becomes optically thick at lower altitudes after which the
nightside layer at $8R_{p}$ does not contribute significantly to the EW
anymore. If we instead had started by calculating the EW using the uppermost
layer of both day- and nightside, the nightside would have contributed to the
EW as the gas would appear to be optically thin still.

From the histogram, we finally define the line-forming radius $r_{line}$ as
the average of the radial bins $r$ weighted by their EW contribution
($\Delta$EW):

$r_{line}=\frac{\sum_{r}r\cdot\Delta\mathrm{EW}(r)}{\sum_{r}\Delta\mathrm{EW}(r)}$
(1)

and define a corresponding standard deviation that gives an indication of the
spread of radii probed

$\sigma_{r}=\sqrt{\frac{\sum_{r^{\prime}}(r^{\prime}-r_{line})^{2}\cdot\Delta\mathrm{EW}(r^{\prime})}{\sum_{r^{\prime}}\Delta\mathrm{EW}(r^{\prime})}}$
(2)

which we do separately for radii above and below $r_{line}$ (i.e.
$r^{\prime}>r_{line}$ and $r^{\prime}<r_{line}$, respectively). Performing
this exercise for each identified strong spectral line allows us to make a
plot that resembles a spectrum but presents the radii at which the lines form
instead of the transit depth.

##  3 Results

###  3.1 Model grid

Our goal is to find the spectral lines that trace upper exoplanet atmospheres
and assess where these spectral lines form. The results generally depend on
the parameters of the system, hence we investigate a grid of systems with
different stellar spectral types and planet gravities. In terms of the
spectral type of the host star, we explore the A0 and A7 models described in
Sec. 2.2 of Oklopčić (2019), which are based on the work of Fossati et al.
(2018), a solar (G2) spectrum from combined TIMED/SEE (Woods et al., 2005) and
SORCE (Rottman, 2005) data, and HD 97658 (K1), HD 85512 (K6) and GJ-832 (M2)
spectra from the MUSCLES survey (France et al., 2016; Youngblood et al., 2016;
P. Loyd et al., 2016). The spectral energy distributions (SEDs) are shown in
Fig. 8. In terms of the planet gravity, we explore a hot Jupiter planet with
the mass and radius akin to HD 209458 b ($R=1.38R_{J}$, $M=0.73M_{J}$,
$\phi=GM/R=10^{12.98}$ erg g-1) and a hot Neptune planet with the mass and
radius akin to HAT-P-11 b ($R=0.389R_{J}$, $M=0.0736M_{J}$, $\phi=10^{12.53}$
erg g-1). We use Parker wind models with mass-loss rates of $\dot{M}=10^{10}$
g s-1 as in Linssen et al. (2022) we retrieved similar values when fitting a
sample of metastable helium observations, but we do explore lower and higher
mass-loss rates in Sec. 4.1. The Parker wind temperatures are found following
the methodology described in Sec. 2.1 and are listed in the captions of the
figures. We fix the semi-major axis of all models to $a=0.05$ AU (changing to
$a=0.025$ AU only had a minor impact on the results). We keep the composition
of the gas fixed at solar. Fig. 2 shows the transmission spectrum of the hot
Jupiter planet around a G2 star. The other transmission spectra are not shown
in this work but are available for download333The transmission spectra and a
reproduction package for all figures of this manuscript are available for
download at https://doi.org/10.5281/zenodo.8013349. Figs. 3 and 4 show the
resulting line-forming radii. Table 1 lists all unique spectral lines
identified. Our model also predicts the hydrogen Ly-$\alpha$ line, but we
caution that our model does not include interactions with the stellar wind
which will shape the observable wings of the line. This makes the calculated
Ly-$\alpha$ transit depth and formation radius unreliable.

Figure 2: Transmission spectrum of the hot Jupiter planet around a G2 star,
for a Parker wind profile with $\dot{M}=10^{10}$ g s-1 and $T=8000$ K. The
strongest line of each atomic/ionic species is marked if it exceeds our
thresholds of 1% at $\lambda<4000$ Å and 0.1% at $\lambda>4000$ Å and its
line-forming radius is shown in the top right panel of Fig. 3. Figure 3: Line-
forming radii for the hot Jupiter planet around stars of different spectral
types. The outflow models used are Parker winds with $\dot{M}=10^{10}$ g/s and
(self-consistent) temperatures of 7000, 7500, 8000, 8000, 8000, 8000 K for the
A0, A7, G2, K1, K6, M2 spectral types, respectively. The marker size indicates
the excess absorption at line center. The figures show all spectral lines
deeper than 1% at $\lambda<4000$ Å and deeper than 0.1% at $\lambda>4000$ Å
(Table 1 lists the wavelength of each line). In every panel, the gray shaded
area demarcates the 4000 Å boundary, the dotted line indicates the Hill radius
and the solid line shows the stellar radius. Gas at larger radii than the
stellar radius can contribute to the transmission spectrum if it is located
‘along the direction of the stellar light rays’, instead of at high impact
parameter. Figure 4: Similar to Fig. 3 but for the hot Neptune planet. The
outflow models used are Parker winds with $\dot{M}=10^{10}$ g/s and (self-
consistent) temperatures of 5000, 5000, 5500, 5000, 5000, 5500 K for the A0,
A7, G2, K1, K6, M2 spectral types, respectively.

###  3.2 General trends with stellar spectral type

We generally find more spectral lines that reach our excess absorption
threshold around later type stars. This is in large part due to the transit
geometry, since the 1% and 0.1% thresholds are reached more easily around late
type stars as their radii are smaller. Fainter lines may however be easier to
observe around earlier type stars as they have a higher luminosity and thus
better signal-to-noise ratio. In terms of the different SEDs, the XUV flux
differs by up to two orders of magnitude (even more for the A0 star) and the
maximum is reached in the spectrum of the K1 star (Fig. 8). At higher XUV
flux, the degree of ionization increases, and we tend to see stronger lines
originating from (highly) ionized species and weaker lines originating from
atomic species. For example, the C IV $\lambda$ 1548 and Si IV $\lambda$ 1393
lines are strongest around a K1 host star. The different XUV flux levels also
result in slightly different self-consistent temperatures and mean molecular
weight structures found for the Parker wind models. The A0 and M2 stars have
the lowest XUV flux, resulting in lower hydrogen and helium ionization and
hence higher $\mu$. At fixed mass-loss rate, a Parker wind model with lower
$T/\mu$ has a denser outflow profile, and thus generally stronger lines.

###  3.3 Dependence of specific lines on stellar spectral type

In the VIS/NIR wavelength range, the strongest tracers are the metastable
helium triplet and the calcium infrared triplet. The metastable helium line
preferably forms around K-type stars as found in Oklopčić (2019) and
corroborated by the current study. Fig. 3 shows that the transit depth of the
calcium infrared triplet remains roughly constant with stellar spectral type,
which means that the total amount of absorption must increase around earlier
type stars (as their radii are larger). This can be explained by the fact that
it originates from excited Ca+ ions for which a main population channel is
absorption of a $\sim$730 nm photon by the ground-state. Earlier type stars
will have a higher flux of such NIR photons. This may make the Ca II infrared
triplet a highly complementary tracer to the metastable helium line by
favoring ‘opposite’ stellar types.

The H$\alpha$ line originates from hydrogen atoms that are excited (among
other processes) through absorption of a stellar Ly-$\alpha$ photon. The A0
stellar SED template has the strongest Ly-$\alpha$ flux and this is why (for
the hot Jupiter planet) we see stronger H$\alpha$ absorption there. Around the
M2 star we see H$\alpha$ due to the favorable radius contrast as discussed
above.

After Ly-$\alpha$, the strongest UV line is C III $\lambda$ 977 around all
stars except A0. Our models predict transit depths roughly half that of
Ly-$\alpha$, yet this line has never been observed. Other lines that are
strong irrespective of stellar type are Si II $\lambda$ 1264, the Mg II
doublet, Si III $\lambda$ 1206 and C II $\lambda$ 1335\. Both silicon lines
have not been reliably detected (see Table 1 and references therein).

###  3.4 Relation between transit depth and line-forming radius

Our results show that while there is certainly a trend in stronger spectral
lines forming at larger radii, this is definitely not true for all lines. A
good example is the N III $\lambda$ 992 line, which forms at somewhat similar
radii as the C III $\lambda$ 977 line, yet is much weaker. For most lines, the
radii probed are a lot larger than one would infer by doing a simple
conversion of transit depth to apparent planet size. Such a conversion
effectively assumes that the atmosphere is completely opaque up until said
planet size but becomes completely transparent above it. Instead of this
scenario, we find that for many lines there is important contribution from
optically thin high altitude layers. Even though the density is low at these
layers, noticeable absorption still occurs due to the total absorbing volume
increasing with $r^{3}$ in a spherically symmetric atmosphere (we discuss this
in more detail in Sec. 4.3). Additionally, the degree of ionization increases
with radius (see Fig. 9) so that the relative (or even absolute) density of
ion species rises with radius. This pushes their line-forming regions higher,
but these lines are not necessarily strong (for example the C IV line at 1548
Å).

##  4 Discussion

###  4.1 Dependence on the mass-loss rate

Some spectral lines that are generally thought to trace the upper atmosphere,
such as e.g. H$\alpha$ and K I $\lambda$ 7667 (e.g. Keles et al., 2020; dos
Santos, 2022), are missing in many panels of Figs. 3 and 4. The Parker wind
profiles do not reach the column density needed for these lines to become
stronger than 0.1%. For the hot Jupiter planet around a G2 star, we also
simulated a Parker wind model with $\dot{M}=10^{11}$ g s-1 and $T=8000$ K. As
Fig. 5 shows, we do see many more lines appearing when using a higher mass-
loss rate, including the aforementioned lines in the optical wavelength range.
Their line-forming radii indicate that most of these lines do in fact trace
the upper atmosphere, but their absence from the sample at lower mass-loss
rates indicates that they may not always be strong enough to be observed. The
figure also shows the results for a Parker wind model with $\dot{M}=10^{9}$ g
s-1 and $T=7000$ K for the same system. Decreasing the mass-loss rate by one
order of magnitude in this way resulted in a spectrum with weaker spectral
lines and thus fewer tracers above our threshold. At the same time, the lines
probe smaller radii compared to the case with $\dot{M}=10^{10}$ g s-1. This
behaviour is remarkably similar for different spectral lines: there is an
almost monotonously increasing relation between the line-forming radii and the
mass-loss rate. The same is true for the transit depth. This means that we do
not identify any spectral lines that are particularly (in)sensitive to the
mass-loss rate.

Figure 5: Line-forming radii for the hot Jupiter planet around a G2 star for
different mass-loss rates. In the left panel we used $\dot{M}=10^{9}$ g s-1
and $T=7000$ K. In the middle panel we used $\dot{M}=10^{10}$ g s-1 and
$T=8000$ K (identical to the upper right panel of Fig. 3). In the right panel
we used $\dot{M}=10^{11}$ g s-1 and $T=8000$ K.

###  4.2 Using an SED of a young star

Thus far we have used SEDs of mature stars. Young stars have much higher XUV
flux levels, thereby driving stronger planetary outflows through
photoevaporation (e.g. Lammer et al., 2003). To investigate the observable
signatures of such a system, we simulated the hot Jupiter planet around a
young K1 type star. We took the MUSCLES spectrum of HD 97658 and multiplied
the X-ray and EUV parts by factors of 2014 and 32, respectively, to match
observations of the XUV flux of the young V1298 Tau system by Maggio et al.
(2023) (the SED is shown in Fig. 8). We assumed a mass-loss rate of
$\dot{M}=10^{12}$ g s-1 and found a self-consistent temperature of $T=12000$
K. The line-forming radii are shown in Fig. 6.

The high stellar XUV flux results in an outflow with a much higher degree of
ionization compared to the mature K1 case (bottom left panel of Fig. 3). This
is reflected in the tracers we find, as the C IV, O VI, N V and Si IV are
among the strongest lines in the spectrum. At the same time, it means that
lines originating from atomic or singly ionized species are not much stronger
than in the mature K1 case, even though the total density is higher due to the
higher mass-loss rate. The high $T/\mu$ also means that the outflow velocity
is higher by a factor of a few, and the spectral tracers should thus be highly
broadened. At fixed mass-loss rate, such an increase in velocity would mean a
decrease in density and thus generally weaker spectral lines. Therefore, the
fact that we still see strong spectral lines in Fig. 6 is because we used a
higher mass-loss rate.

Figure 6: Line-forming radii for the hot Jupiter planet around a young K1
star, using a Parker wind model with $\dot{M}=10^{12}$ g s-1 and $T=12000$ K.
The stellar spectrum was constructed by scaling the X-ray and EUV portions of
the MUSCLES spectrum of HD 97658 to match the observed high-energy flux of
V1298 Tau.

###  4.3 Dependence on the computational domain and outflow geometry

It is important to consider the impact of the upper boundary of the
computational domain for lines that form at high altitudes. Our Cloudy
simulations extend to 8 planetary radii, even if the Hill radius of the planet
is smaller. As the results show, many lines have significant contributions
from material beyond the Hill sphere, where the assumption of a 1D spherically
symmetric outflow is likely not valid. We chose to still include the material
beyond the Hill sphere to allow for a qualitative comparison between the line-
forming radii of different spectral lines, as well as between different
stellar SED shapes. In Fig. 7 we show line-forming radii for the hot Jupiter
planet around a G2 star when we stop our calculations at 4 planetary radii (a
bit below the Hill radius of 4.6$R_{p}$). For the most part, the middle two
panels of Fig. 7 appear as ‘vertically squashed’ versions of the corresponding
panels in Fig. 5, i.e. the lines that form significantly below the Hill radius
are hardly affected, while the high-altitude lines are less strong and are now
predicted to form at smaller radii. This means that for these lines, similar
to the Ly-$\alpha$ line, the 3D geometry of the outflow will presumably have a
strong influence on the line shape, strength, and where it forms. The 3D
outflow geometry and resulting observational signatures can depend on many
factors, including the stellar wind (McCann et al., 2019), radiation pressure
(Bourrier & Lecavelier des Etangs, 2013), flares and coronal mass ejections
(Hazra et al., 2021), as well as the planet day- and nightside differences
(Carroll-Nellenback et al., 2017) and magnetic field (Khodachenko et al.,
2021). These effects may shape the outflow into highly non-spherical
configurations such as an up- and downstream two-arm structure or a cometary
tail. We note that the upper boundary of the atmosphere is a free parameter in
our simulations, but observations of high-altitude spectral lines should be
able to constrain the extent of the absorbing material, for example by
measuring the time of ingress. To provide a rough estimate, the hot Jupiter
planet at a distance of 0.05 AU from a G2 star has an orbital velocity of
around 1$R_{p}$ per 0.2 hours, so that lines forming around 5$R_{p}$ may show
significant absorption already an hour before transit.

Figure 7: Line-forming radii for the hot Jupiter planet around a G2 star for
different outflow profiles and with calculations truncated at 4$R_{p}$. The
top panel shows the density and velocity profiles (dashed). The second panel
shows the result for a Parker wind profile with $\dot{M}=10^{10}$ g s-1 and
$T=8000$ (same system as in the upper right panel of Fig. 3, but now stopped
below the Hill radius). The third panel shows the result for a Parker wind
profile with $\dot{M}=10^{11}$ g s-1 and $T=8000$. The bottom panel shows the
results for ATES, which predicted a mass-loss rate of $\dot{M}=10^{11}$ g s-1.
We find additional lines at small radii ($\lesssim 1.1R_{p}$) when using ATES
because its density profile increases steeply in the lower atmosphere.

For our simulations of the hot Neptune planet (Fig. 4) most lines seem to
cluster around a line-forming radius of $\sim$4.5$R_{p}$ with large standard
deviations. For this planet, the Parker wind density profiles are less ‘steep’
than for the hot Jupiter, and thus the relative contribution of material at
larger radii increases. In fact, the bulk density decreases by roughly three
orders of magnitude from $1R_{p}$ to $8R_{p}$, but the volume of a spherical
shell increases by a factor $8^{3}\approx 500$, so that the total absorbing
gas mass hardly changes as a function of radius. This makes the spectral line
formation particularly sensitive to the upper boundary of the simulation,
especially if it originates from a species whose relative abundance does not
vary strongly with radius. The histogram in the bottom panel of Fig. 1 will
look flat for such a spectral line and it will thus tend to form around the
middle of the radial domain ($4.5R_{p}$ for our simulations) with little
dependence on the line strength. Similarly to the high-altitude lines
discussed in the paragraph above, the outflow geometry will ultimately shape
the lines for this type of planet.

###  4.4 Using an ATES outflow profile

In Sec. 4.1, we investigated how the spectral tracers depend on the mass-loss
rate of the Parker wind profile. Here, we study the effect of using an outflow
model generated with the photoionization hydrodynamics code ATES (Caldiroli et
al., 2021). The main difference with a Parker wind model is that ATES
additionally solves the equation of energy conservation such that it self-
consistently calculates the temperature structure instead of assuming an
isothermal profile (see also Sec. 2 in Linssen et al. 2022). In doing so, ATES
predicts the mass-loss rate instead of assuming it, but its atmospheric
structure is generally similar to a Parker wind model with the same mass-loss
rate and a comparable temperature (see top panel of Fig. 7). ATES does however
have the planet equilibrium temperature and a manually chosen density as lower
boundary conditions at 1$R_{p}$. The equilibrium temperature is much lower
than the thermospheric temperature and the density at planet radius is
typically chosen as $10^{14}$ cm-3 (Caldiroli et al., 2021, 2022). Together,
this results in a much steeper and denser lower atmospheric ($\lesssim
1.1R_{p}$) density structure than that of Parker wind models. Additionally,
the velocity profile of ATES tends to be steeper at smaller radii, as it
transitions into a hydrostatic (i.e. velocity$\sim$0) lower atmosphere, while
a Parker wind model does not necessarily do so. Precise measurements of
spectral line broadening might in this way be able to distinguish between
these two outflow profiles.

We ran ATES for our hot Jupiter planet around a G2 star and post-processed the
resulting density and temperature structures with Cloudy in order to make
synthetic spectra. ATES predicts a mass-loss rate of $10^{11}$ g s-1 for this
system. We compare the spectral tracers in Fig. 7. As expected, the results
are very similar. The Parker wind model predicts a few lines to be stronger at
radii of roughly 1.5$R_{p}$ such as e.g. Cu II $\lambda$ 1359 and Cr I
$\lambda$ 3580\. The ATES model predicts a few additional lines at small radii
at optical wavelengths, such as e.g. V I $\lambda$ 4578 and Li I $\lambda$
6710\. Both these differences can be explained by a higher density of the
considered model at these radii. The sub-1.1$R_{p}$ radii of the extra lines
found by ATES indicate that they do not really probe the escaping atmosphere,
but rather the transition region from hydrostatic ‘lower’ atmosphere to
hydrodynamic ‘upper’ atmosphere. We caution that this modeled region is
somewhat uncertain; on our end because the lines form in only the few deepest
computational atmospheric layers and hence the calculated radii are not very
precise, but also because one has to assume and fix a density at the lower
boundary of the computational domain with ATES, which allows arbitrarily
scaling the strength of these spectral lines if a well-motivated or
observationally constrained value is lacking. Salz et al. (2016) showed in
their Fig. 4 that this ‘reference’ density does not influence the structure of
the outflow above $\gtrsim 1.1R_{p}$ (for TPCI simulations, which produce
similar results to ATES).

As explained in Sec. 2.1, we chose to assess the self-consistency of the
Parker wind models in the 1-2$R_{p}$ region. As the results throughout this
work show, we find many spectral lines that originate at radii significantly
larger than this, where we in principle may deviate from the self-consistent
temperature and hence density and velocity structure, as these are coupled.
However, as the top panel of Fig. 7 shows, even at these larger radii, the
structure of Parker wind models are relatively similar to ATES models of
comparable mass-loss rate. We thus do not consider the fact that we have to
make a choice about where to assess the self-consistency of the Parker wind
models to be a major limitation of our approach.

###  4.5 Dependence on the composition

In this work, we assumed a solar composition atmosphere. A change in the
metallicity may have important consequences for the outflow structure and
observed spectrum. The structure will change because the metallicity affects
the mean molecular weight, but also because metals can be efficient coolants
in the upper atmosphere (Fossati et al., 2021) and hence influence the
temperature structure of the outflow. The spectrum will change as a result of
the different outflow profile, but also because the strength of spectral lines
originating from metal species simply scales with the metal abundance,
although not necessarily linearly. Especially smaller planets can be expected
to have (highly) supersolar metallicities (e.g. Fortney et al., 2013;
Thorngren et al., 2016; Wakeford et al., 2017) and the observed signatures may
depend heavily on it. At the same time it may prove important to consider a
gas composition that is not constant with radius. Heavier elements are more
difficult to drag along in the planetary wind and they may preferentially
reside at lower altitudes. Specific elements that are easily condensed into
clouds in the lower atmosphere may also be less abundant in the upper
atmosphere. We will study non-solar composition outflows in a future work.

###  4.6 Targeting complementary spectral lines

In general, observing multiple lines for the same planet will put stronger
constraints on models and model parameters, as a given model then has to be
able to fit these lines simultaneously. Huang et al. (2023) recently presented
such a multi-line analysis for combined UV and optical observations of
different escaping species. They highlighted the fact that combining several
features helped to break the degeneracy between their model parameters. Some
spectral lines complement each other well and are particularly useful in
constraining theoretical models. The line-forming radii that we calculate in
this work allow us to explore such specific synergies between different
spectral tracers. For example, lines that form at similar radii should also
have similar temperatures and outflow velocities, which will both contribute
to spectral line broadening. Since the thermal broadening also depends on the
mass of the absorbing species, observing two spectral lines from species with
different masses may help to reveal the dominant broadening mechanism of the
lines (in absence of other line broadening effects). Namely, such lines
forming at similar radii should have similar broadening if the velocity
dominates, but have different broadening due to the different mass if the
temperature dominates. An example of this is the (heavy) Fe II $\lambda$ 2382
multiplet combined with the (lighter) He I $\lambda$ 10833 triplet or O I
$\lambda$ 1302 lines.

It will also be useful to observe lines that form at different radii, which
may help to constrain the temperature and/or velocity profile of the outflow.
An example of this is combining the C I $\lambda$ 1657, C II $\lambda$ 1335, C
III $\lambda$ 977 and C IV $\lambda$ 1548 lines. These lines span a large
range of radii and little assumptions would have to be made about the
composition of the upper atmosphere as they all originate from the same
element. This series of lines would additionally allow one to study the drag
of carbon to high altitudes in the outflow (which is currently included in our
model through a constant composition with radius). The same could be done for
a heavier species by targeting the Si I $\lambda$ 1850, Si II $\lambda$ 1264,
Si III $\lambda$ 1206 and Si IV $\lambda$ 1393 series.

##  5 Summary

Escape of exoplanet atmospheres has thus far been observed in a handful of
spectral lines, with the Ly-$\alpha$ line and metastable helium triplet being
the most prolific. Expanding on this suite of spectral lines would help to
better constrain the many different models used to interpret observations and
get a more complete view of upper atmospheres. In this work, we aim to find
additional spectral lines tracing upper atmospheres that could be targeted
with observations. We do this by making synthetic transmission spectra of
hydrodynamically escaping atmospheres and identifying the strongest spectral
lines.

We describe the upper atmosphere as a solar composition 1D Parker wind and use
the NLTE photoionization code Cloudy to obtain the radial density profiles of
each absorbing species. We post-process these density profiles with the
available line coefficients in the NIST database to compute the transmission
spectrum from the FUV to the NIR. We make transit spectra for generic hot
Jupiter and hot Neptune planets orbiting stars of spectral types ranging from
A0 to M2. For every transit spectrum, we construct a sample of suitable escape
tracers that includes the strongest spectral line of each atomic/ionic species
if it has a transit depth stronger than 1% or 0.1% in the UV or VIS/NIR,
respectively. In this way, we identify many spectral lines that have not yet
been observed in an exoplanet atmosphere, most of which are at ultraviolet
wavelengths.

We want to know which parts of the outflow are probed by a given spectral line
and thus define a line-forming radius based on the contribution of each radial
atmospheric layer to the ‘observed’ spectral line. This illustrates that the
different tracers we find probe different line-forming radii, from the ‘lower’
regions where the outflow is launched to the ‘upper’ regions where
interactions with the stellar environment are expected to shape the lines. The
information of where the lines form helps to find potential synergies between
them, and we give a few examples of how one might leverage observations of
multiple lines to better constrain the structure of the escaping atmosphere.

###### Acknowledgements.

We are grateful to A. Caldiroli and the ATES team, as well as G.J. Ferland and
the Cloudy team, for making their codes public and providing support. We
appreciate the helpful conversations with M. MacLeod and F. Nail. We thank
SURFsara (www.surfsara.nl) for the support in using the Lisa Compute Cluster.
AO gratefully acknowledges support from the Dutch Research Council NWO Veni
grant.

## References

  * Allart et al. (2018) Allart, R., Bourrier, V., Lovis, C., et al. 2018, Sci, 362, 1384
  * Ballester & Ben-Jaffel (2015) Ballester, G. E. & Ben-Jaffel, L. 2015, ApJ, 804, 116
  * Ben-Jaffel & Ballester (2013) Ben-Jaffel, L. & Ballester, G. E. 2013, A&A, 553, A52
  * Ben-Jaffel et al. (2022) Ben-Jaffel, L., Ballester, G. E., Muñoz, A. G., et al. 2022, Nat Astron, 6, 141
  * Bourrier & Lecavelier des Etangs (2013) Bourrier, V. & Lecavelier des Etangs, A. 2013, A&A, 557, A124
  * Caldiroli et al. (2021) Caldiroli, A., Haardt, F., Gallo, E., et al. 2021, A&A, 655, A30
  * Caldiroli et al. (2022) Caldiroli, A., Haardt, F., Gallo, E., et al. 2022, A&A, 663, A122
  * Carroll-Nellenback et al. (2017) Carroll-Nellenback, J., Frank, A., Liu, B., et al. 2017, MNRAS, 466, 2458
  * Casasayas-Barris et al. (2019) Casasayas-Barris, N., Pallé, E., Yan, F., et al. 2019, A&A, 628, A9
  * Charbonneau et al. (2002) Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377
  * Cubillos et al. (2023) Cubillos, P. E., Fossati, L., Koskinen, T., et al. 2023, A&A, 671, A170
  * Cubillos et al. (2020) Cubillos, P. E., Fossati, L., Koskinen, T., et al. 2020, AJ, 159, 111
  * dos Santos (2022) dos Santos, L. A. 2022, arXiv:2211.16243
  * dos Santos et al. (2022) dos Santos, L. A., Vidotto, A. A., Vissapragada, S., et al. 2022, A&A, 659, A62
  * Ferland et al. (1998) Ferland, G., Korista, K., Verner, D., et al. 1998, PASP, 110, 761
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzman, F., et al. 2017, Rev. Mex. Astron. Astrofis., 53, 385
  * Fortney et al. (2013) Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80
  * Fossati et al. (2010) Fossati, L., Haswell, C. A., Froning, C. S., et al. 2010, ApJ, 714, L222
  * Fossati et al. (2018) Fossati, L., Koskinen, T., Lothringer, J. D., et al. 2018, ApJ, 868, L30
  * Fossati et al. (2021) Fossati, L., Young, M. E., Shulyak, D., et al. 2021, A&A, 653, A52
  * France et al. (2016) France, K., Loyd, R. O. P., Youngblood, A., et al. 2016, ApJ, 820, 89
  * García Muñoz et al. (2021) García Muñoz, A., Fossati, L., Youngblood, A., et al. 2021, ApJ, 907, L36
  * Ginzburg et al. (2018) Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS, 476, 759
  * Hazra et al. (2021) Hazra, G., Vidotto, A. A., Carolan, S., Villarreal D’Angelo, C., & Manchester, W. 2021, MNRAS, 509, 5858
  * Huang et al. (2023) Huang, C., Koskinen, T., Lavvas, P., & Fossati, L. 2023, arXiv:2304.07352
  * Ionov et al. (2018) Ionov, D. E., Pavlyuchenkov, Y. N., & Shematovich, V. I. 2018, MNRAS, 476, 5639
  * Jensen et al. (2012) Jensen, A. G., Redfield, S., Endl, M., et al. 2012, ApJ, 751, 86
  * Keles et al. (2020) Keles, E., Kitzmann, D., Mallonn, M., et al. 2020, MNRAS, 498, 1023
  * Khodachenko et al. (2021) Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., et al. 2021, MNRAS, 507, 3626
  * Lamers & Cassinelli (1999) Lamers, H. J. G. L. M. & Cassinelli, J. P. 1999, Introduction to Stellar Winds (Cambridge University Press)
  * Lammer et al. (2003) Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121
  * Lecavelier des Etangs et al. (2010) Lecavelier des Etangs, A., Ehrenreich, D., Vidal-Madjar, A., et al. 2010, A&A, 514, A72
  * Linsky et al. (2010) Linsky, J. L., Yang, H., France, K., et al. 2010, ApJ, 717, 1291
  * Linssen et al. (2022) Linssen, D. C., Oklopčić, A., & MacLeod, M. 2022, A&A, 667, A54
  * Lundkvist et al. (2016) Lundkvist, M. S., Kjeldsen, H., Albrecht, S., et al. 2016, Nat Commun, 7, 11201
  * Maggio et al. (2023) Maggio, A., Pillitteri, I., Argiroffi, C., et al. 2023, arXiv:2305.06931
  * Mansfield et al. (2018) Mansfield, M., Bean, J. L., Oklopčić, A., et al. 2018, ApJ, 868, L34
  * McCann et al. (2019) McCann, J., Murray-Clay, R. A., Kratter, K., & Krumholz, M. R. 2019, ApJ, 873, 89
  * Murray-Clay et al. (2009) Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23
  * Nortmann et al. (2018) Nortmann, L., Pallé, E., Salz, M., et al. 2018, Sci, 362, 1388
  * Oklopčić (2019) Oklopčić, A. 2019, ApJ, 881, 133
  * Oklopčić & Hirata (2018) Oklopčić, A. & Hirata, C. M. 2018, ApJ, 855, L11
  * Owen (2019) Owen, J. E. 2019, Annu. Rev. Earth Planet. Sci., 47, 67
  * Owen et al. (2022) Owen, J. E., Murray-Clay, R. A., Schreyer, E., et al. 2022, MNRAS, 518, 4357
  * Owen & Wu (2013) Owen, J. E. & Wu, Y. 2013, ApJ, 775, 105
  * Owen & Wu (2017) Owen, J. E. & Wu, Y. 2017, ApJ, 847, 29
  * P. Loyd et al. (2016) P. Loyd, R. O., France, K., Youngblood, A., et al. 2016, ApJ, 824, 102
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664
  * Rottman (2005) Rottman, G. 2005, Solar Physics, 230, 7
  * Salz et al. (2018) Salz, M., Czesla, S., Schneider, P. C., et al. 2018, A&A, 620, A97
  * Salz et al. (2016) Salz, M., Czesla, S., Schneider, P. C., & Schmitt, J. H. M. M. 2016, A&A, 586, A75
  * Schlawin et al. (2010) Schlawin, E., Agol, E., Walkowicz, L. M., Covey, K., & Lloyd, J. P. 2010, ApJ, 722, L75
  * Sing et al. (2011) Sing, D. K., Désert, J.-M., Fortney, J. J., et al. 2011, A&A, 527, A73
  * Sing et al. (2019) Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, AJ, 158, 91
  * Spake et al. (2018) Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nat, 557, 68
  * Thorngren et al. (2016) Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64
  * Turner et al. (2016) Turner, J. D., Christie, D., Arras, P., Johnson, R. E., & Schmidt, C. 2016, MNRAS, 458, 3880
  * Vidal-Madjar et al. (2003) Vidal-Madjar, A., des Etangs, A. L., Désert, J.-M., et al. 2003, Nat, 422, 143
  * Vidal-Madjar et al. (2004) Vidal-Madjar, A., Dsert, J.-M., Etangs, A. L. d., et al. 2004, ApJ, 604, L69
  * Vidal-Madjar et al. (2013) Vidal-Madjar, A., Huitson, C. M., Bourrier, V., et al. 2013, A&A, 560, A54
  * Vissapragada et al. (2022) Vissapragada, S., Knutson, H. A., Greklek-McKeon, M., et al. 2022, AJ, 164, 234
  * Wakeford et al. (2017) Wakeford, H. R., Sing, D. K., Kataria, T., et al. 2017, Sci, 356, 628
  * Wang & Dai (2021) Wang, L. & Dai, F. 2021, ApJ, 914, 98
  * Woods et al. (2005) Woods, T. N., Eparvier, F. G., Bailey, S. M., et al. 2005, JGR, 110
  * Yan et al. (2019) Yan, F., Casasayas-Barris, N., Molaverdikhani, K., et al. 2019, A&A, 632, A69
  * Youngblood et al. (2016) Youngblood, A., France, K., Loyd, R. O. P., et al. 2016, ApJ, 824, 101
  * Zhang et al. (2022) Zhang, Y., Snellen, I. A. G., Wyttenbach, A., et al. 2022, A&A, 666, A47

##  Appendix A Additional figures

Figure 8: The stellar SEDs used throughout this work. References for these
spectra are in Sec. 3.1 and in Sec. 4.2 for the young K1 spectrum.

##  Appendix B Calculating synthetic transmission spectra

In this work, we do not use Cloudy’s own predicted spectrum, as it has a
spectral resolution that is too low for our purposes and furthermore does not
include the effects of Doppler broadening due to the outflow velocity.
Instead, we use Cloudy to obtain the densities of the different energy levels
of all atoms and ions up to six times ionized (higher degrees of ionization
are typically not reached in our simulations) as a function of radius. For
most metal species, Cloudy resolves the lowest $\sim$15 energy levels, which
is satisfactory since we are interested in strong lines that usually originate
from the ground state or the lower excited states. Fig. 9 shows a subset of
the density structures of different carbon atomic/ionic energy levels. For
each of these atomic and ionic species, we download a table of spectral lines
from the NIST database. For every spectral line, these tables specify the
wavelength, the energy level of the atom/ion it originates from, the
oscillator strength (used to calculate the cross-section) and the Einstein A
coefficient (used to calculate the natural broadening).

Figure 9: Density structures of different carbon ions and their energy levels
as returned by Cloudy for the hot Jupiter planet around a G2 star. The blue,
pink, green and yellow colors correspond to atomic, single, double and triple
ionized carbon, respectively. For each of these, the density structures of the
four lowest energy levels are shown in different line styles. We post-process
such density structures with the coefficients from the NIST database to
produce synthetic transmission spectra. For example, the C I $\lambda$ 1657
multiplet is a blend of lines originating from the ${}^{3}P_{0}$,
${}^{3}P_{1}$ and ${}^{3}P_{2}$ levels, while the C III $\lambda$ 977 line
originates solely from the ground-state (i.e. ${}^{1}S_{0}$).

Ultimately, the Parker wind velocity structure, Cloudy’s temperature structure
and density structures of the energy levels of each atomic/ionic species,
together with the spectral line coefficients from the NIST database, allow us
to perform radiative transfer calculations using the equations in Section 3.4
of Oklopčić & Hirata (2018). In summary, we project the 1D density, velocity
and temperature profiles onto a 2D semi-circle and calculate the optical depth
$\tau$ at frequency $\nu$ along rays with different impact parameters $y$ from
the planet as

$\tau_{\nu}(y)=\int n(x)\sigma_{0}\Phi(\Delta\nu)\mathrm{d}x,$ (3)

where $n$ is the number density of the absorbing species, $\sigma_{0}$ the
line cross-section, $\Phi$ the Voigt line profile and $\Delta\nu$ the offset
from line center which depends on $\nu$ and the line-of-sight velocity. The
extinction of each ray is $1-\mathrm{exp}(-\tau_{\nu}(y))$ and we find the
total transit spectrum $F_{in}/F_{out}$ by adding the extinction of all rays
multiplied with their projected surface on the stellar disk (which would point
‘out of the paper’ in Fig. 1 so that we effectively have a 3D atmosphere). We
assume a zero transit impact parameter and no stellar limb darkening.

##  Appendix C Identified lines

Table 1 lists all unique spectral lines that appear in Figs. 3 and 4 and
whether they have been previously detected. Our simulations with Cloudy
included all elements up until zinc. We post-processed these simulations to
make spectra of all atoms and ions up to 6+. However, for some of these ions,
the NIST database does not report any spectral lines with the coefficients
needed to produce spectra. Thus, we were not able to produce spectra for the
following atomic and ionic species: F, Ti4+, Ti6+, V3+, V4+, V5+, V6+, Cr2+,
Cr3+, Cr4+, Cr5+, Cr6+, Mn2+, Mn3+, Mn4+, Mn5+, Mn6+, Fe3+, Fe5+, Co3+, Co4+,
Co5+, Co6+, Ni3+, Ni5+, Ni6+, Cu2+, Cu3+, Cu4+, Cu5+, Cu6+, Zn+, Zn2+, Zn3+,
Zn4+, Zn5+, Zn6+.

Table 1: All unique identified spectral lines tracing upper atmospheres of
Figs. 3, 4, 5 and 6 (thus excluding the additional lines found in Fig. 7 when
using an ATES model). We give only one example reference for lines that have
been previously detected.

Wavelength | Species | Previously detected
---|---|---
(vac) (Å) |  |
912.7 | S II | No
916.7 | N II | No
933.4 | S VI | No
977.0 | C III | No
991.6 | N III | No
1021.3 | S III | No
1031.9 | O VI | No
1048.2 | Ar I | No
1073.0 | S IV | No
1122.5 | Fe III | No
1199.5 | N I | No
1206.5 | Si III | e.g. Linsky et al. (2010)11$1$Contested by Ballester & Ben-Jaffel (2015)
1215.7 | H I | e.g. Vidal-Madjar et al. (2003)
1238.8 | N V | No
1264.7 | Si II | No
1302.2 | O I | e.g. Vidal-Madjar et al. (2004)
1310.7 | P II | No
1335.7 | C II | e.g. Vidal-Madjar et al. (2004)
1358.8 | Cu II | No
1393.8 | Si IV | e.g. Schlawin et al. (2010)22$2$Tentative detection
1425.0 | S I | No
1548.2 | C IV | No
1657.0 | C I | No
1670.8 | Al II | No
1850.7 | Si I | No
1854.7 | Al III | No
2056.3 | Cr II | No
2139.2 | Zn I | No
2217.2 | Ni II | No
2276.2 | Ca I | No
2286.9 | Co II | No
2382.8 | Fe II | e.g. Sing et al. (2019)
2408.0 | Co I | No
2576.9 | Mn II | No
2796.4 | Mg II | e.g. Fossati et al. (2010)
2853.0 | Mg I | e.g. Vidal-Madjar et al. (2013)33$3$Contested by Cubillos et al. (2020)
3248.5 | Cu I | No
3384.7 | Ti II | No
3579.7 | Cr I | No
3889.7 | He I | No
3934.8 | Ca II | e.g. Yan et al. (2019)
4227.9 | Ca I | No
4255.5 | Cr I | No
5170.5 | Fe II | e.g. Casasayas-Barris et al. (2019)
5185.0 | Mg I | e.g. Zhang et al. (2022)
5891.6 | Na I | e.g. Charbonneau et al. (2002)
6564.5 | H I | e.g. Jensen et al. (2012)
7667.0 | K I | e.g. Sing et al. (2011)
8544.4 | Ca II | e.g. Casasayas-Barris et al. (2019)
10833.3 | He I | e.g. Spake et al. (2018)

444


* (87) P. Artoisenet, R. Frederix, O. Mattelaer and R. Rietkerk, _Automatic spin-entangled decays of heavy resonances in Monte Carlo simulations_ , _JHEP_ 03 (2013) 015 [1212.3460].
  * (88) S. Catani, B. R. Webber and G. Marchesini, _QCD coherent branching and semiinclusive processes at large x_ , _Nucl. Phys._ B349 (1991) 635.
  * (89) A. Banfi, G. P. Salam and G. Zanderighi, _Principles of general final-state resummation and automated implementation_ , _JHEP_ 03 (2005) 073 [hep-ph/0407286].
  * (90) R. Bonciani, S. Catani, M. L. Mangano and P. Nason, _Sudakov resummation of multiparton QCD cross-sections_ , _Phys. Lett._ B575 (2003) 268 [hep-ph/0307035].
  * (91) S. Frixione, P. Nason and C. Oleari, _Matching NLO QCD computations with Parton Shower simulations: the POWHEG method_ , _JHEP_ 11 (2007) 070 [0709.2092].
  * (92) G. Rodrigo and F. Krauss, _Resummed jet rates for heavy quark production in e+ e- annihilation_ , _Eur. Phys. J._ C33 (2004) S457 [hep-ph/0309325].
  * (93) G. Luisoni, P. Nason, C. Oleari and F. Tramontano, _$HW^{\pm}$ /HZ + 0 and 1 jet at NLO with the POWHEG BOX interfaced to GoSam and their merging within MiNLO_, _JHEP_ 10 (2013) 083 [1306.2542].
  * (94) K. Hamilton, T. Melia, P. F. Monni, E. Re and G. Zanderighi, _Merging WW and WW+jet with MINLO_ , _JHEP_ 09 (2016) 057 [1606.07062].
  * (95) S. Frixione, Z. Kunszt and A. Signer, _Three jet cross-sections to next-to-leading order_ , _Nucl. Phys._ B467 (1996) 399 [hep-ph/9512328].
  * (96) N. Hansen and A. Ostermeier, _Completely derandomized self-adaptation in evolution strategies_ , _Evolutionary Computation_ (2001) 159.
  * (97) A. Banfi, P. F. Monni, G. P. Salam and G. Zanderighi, _Higgs and Z-boson production with a jet veto_ , _Phys. Rev. Lett._ 109 (2012) 202001 [1206.4998].
  * (98) NNPDF collaboration, R. D. Ball et al., _Parton distributions for the LHC Run II_ , _JHEP_ 04 (2015) 040 [1410.8849].
  * (99) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J._ C75 (2015) 132 [1412.7420].
  * (100) Particle Data Group collaboration, C. Patrignani et al., _Review of Particle Physics_ , _Chin. Phys._ C40 (2016) 100001.
  * (101) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An Introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * (102) B. Cabouat and T. Sjöstrand, _Some Dipole Shower Studies_ , _Eur. Phys. J._ C78 (2018) 226 [1710.00391].
  * (103) A. Ballestrero et al., _Precise predictions for same-sign W-boson scattering at the LHC_ , 1803.07943.
  * (104) F. Krauss and G. Rodrigo, _Resummed jet rates for e+ e- annihilation into massive quarks_ , _Phys. Lett._ B576 (2003) 135 [hep-ph/0303038].
  * (105) S. Catani, D. de Florian and M. Grazzini, _Universality of nonleading logarithmic contributions in transverse momentum distributions_ , _Nucl. Phys._ B596 (2001) 299 [hep-ph/0008184].
  * (106) D. de Florian and M. Grazzini, _The Structure of large logarithmic corrections at small transverse momentum in hadronic collisions_ , _Nucl. Phys._ B616 (2001) 247 [hep-ph/0108273].


# MTRGL: Effective Temporal Correlation Discerning through Multi-modal
Temporal Relational Graph Learning

###### Abstract

In this study, we explore the synergy of deep learning and financial market
applications, focusing on pair trading. This market-neutral strategy is
integral to quantitative finance and is apt for advanced deep-learning
techniques. A pivotal challenge in pair trading is discerning temporal
correlations among entities, necessitating the integration of diverse data
modalities. Addressing this, we introduce a novel framework, Multi-modal
Temporal Relation Graph Learning (MTRGL). MTRGL combines time series data and
discrete features into a temporal graph and employs a memory-based temporal
graph neural network. This approach reframes temporal correlation
identification as a temporal graph link prediction task, which has shown
empirical success. Our experiments on real-world datasets confirm the superior
performance of MTRGL, emphasizing its promise in refining automated pair
trading strategies. ††footnotetext: *Corresponding author: Shan Wu,
<EMAIL_ADDRESS>To appear as a conference paper at ICASSP 2024

Index Terms—  Graph, finance, multimodal

##  1 Introduction

Pair trading, a pivotal investment strategy, capitalizes on price disparities
between related assets or markets [1]. Traders identify temporary price
variances, executing simultaneous long and short positions on correlated
assets. In a long position, a trader anticipates a price rise, while in a
short position, they expect a decline, aiming to repurchase the asset at a
reduced price later. This strategy, illustrated in Figure 1, leverages market
inefficiencies, enabling traders to profit from transient price deviations. By
doing so, pair trading not only enhances market efficiency but also fosters
improved returns and risk management in portfolios.

The essence of pair trading lies in identifying temporal correlations among
assets or markets. Recognizing pairs with synchronized price movements, which
tend to converge or diverge, is complex due to the market’s dynamic nature and
the vast number of potential pairs. This complexity demands robust
quantitative analysis and advanced statistical techniques, further complicated
by the evolving nature of asset relationships influenced by market shifts,
regulations, and macroeconomic events.

Fig. 1: Illustration of the Pair Trading Strategy. Entities A and B are two
correlated entities whose price tends to move together. The circled section
indicates a temporal divergence of the price of the two entities due to market
fluctuation or inefficiency. During the divergent period, the trader can
simultaneously take a long position on entity B and a short position on entity
A, and benefit from the divergence period when the movement of A and B
converge again.

Historically, experts manually identified these correlations. While
traditional methods focused on statistical inferences, they often captured
only basic data trends. The rise of machine learning has ushered in a new era,
demonstrating its prowess in handling vast datasets, recognizing non-linear
relationships, and outperforming traditional methods. However, the potential
of machine learning in pair trading remains underexplored, necessitating
further research.

Applying machine learning to pair trading introduces challenges. Firstly,
while abundant data exists in finance, its simplicity can limit machine
learning optimization. Advanced models thrive on complex, multi-dimensional
inputs, surpassing traditional methods. Thus, integrating diverse data
modalities is crucial. Secondly, the ever-changing financial markets,
influenced by various factors, can drastically alter asset relationships. This
demands models that can adapt to these temporal correlations, requiring
advanced algorithms.

In this study, we examine the role of machine learning in the stock market by
focusing on pair trading, introducing a novel framework, the multi-modal
temporal relation graph learning (MTRGL). The contributions of this paper are
highlighted as follows:

$\bullet$ We probe the challenges associated with the application of machine
learning in finance. We highlight the pivotal importance of integrating
information from diverse modalities.

$\bullet$ The introduction of MTRGL, a framework that seamlessly integrates
various modalities. It features two key components: a mechanism for
constructing a dynamic graph encapsulating both time series data (e.g., price
trends) and discrete feature information (e.g., sector classifications), and a
neural model powered by a memory-based dynamic graph neural network, an
efficacious tool in temporal graph learning [2].

$\bullet$ Empirical analysis on real-world data, showcasing MTRGL’s
superiority in identifying temporal correlations. An ablation study further
validates its effectiveness, offering insights for future machine learning
applications in finance.

##  2 RELATED WORK

The rise in computational capabilities has enabled neural-network-based
machine learning models to derive insights from complex datasets, including
images [3, 4, 5], languages [6, 7, 8], and networks [9, 10, 11, 12]. Financial
markets, with their intricate dynamics, have become a focal point for machine
learning applications. The allure of financial gains has driven increased
research interest in this domain, with a comprehensive review available in
[13]. Yet, much of the existing work has been confined to stock price data, a
relatively straightforward time series, restricting the full exploitation of
machine learning’s capabilities.

Recent efforts [14, 15] have explored the capabilities of large language
models like ChatGPT [16, 17] to extract features from financial news. However,
the combination of these features with stock price data often remains
rudimentary. In this study, we address this integration challenge, emphasizing
temporal correlations for pair trading with MTRGL. Our approach presents a
refined method to integrate multi-modal information, creating structure within
financial markets.

(a) Temporal Graph Construction Process of MTRGL

(b) Training Pipeline of MTRGL

Fig. 2:  Figure 2(a) displays how the time series data spanning from $[0,T]$
is partitioned into a sequence of intervals $\\{\Delta_{1},...,\Delta_{n}\\}$.
The latter portion of this figure presents how the information from these time
series segments is amalgamated with feature data to create the sequence of
temporal investment graphs $\\{G_{1},...,G_{n}\\}$. On the other hand, Figure
2(b) provides a visualization of MTRGL’s training process. The constructed
graphs are processed as event batches and employed for model updates. The loss
of event batches $B_{i}$, denoted as $\mathcal{L}(B_{i},\theta_{k}^{(j)})$, is
calculated using the parameter $\theta_{k}^{(j)}$ obtained from the $k$-th
epoch and $j$-th iteration. The parameters derived from the final iteration of
each epoch serve as initial values for the subsequent epoch, denoted as
$\theta_{k}^{(K-1)}=\theta_{k+1}^{(0)}$.

##  3 Methodology

Consider $\mathcal{A}=\\{A_{1},A_{2},...,A_{n}\\}$ as a collection of $n$
entities under examination, where entities may represent companies, assets, or
markets. Each entity $A_{i}$ is linked with a dynamic feature vector
$X_{i}(t)=(x^{(1)}(t),x^{(2)}(t),...,x^{(m)}(t))$, with each entry
$x^{(j)}(t)$ signifying the $j$-th feature of entity $A_{i}$ at a given time
$t$ (for instance, sector name, trading volume, market share, and so on). We
use $\mathcal{X}(t)$ to denote the collection of feature vectors associated
with $\mathcal{A}$ at time $t$. Moreover, $P_{i}(t):t\mapsto\mathbb{R}$ is the
trading price time series of entity $A_{i}$, and we denote $P_{i}([0,t])$ as
the time series information of entity $A_{i}$ in the time interval $[0,t]$ for
a certain $t>0$. Let $\mathcal{P}=\\{P_{i}(t)\\}$ denote the set of time
series linked to the entity set $\mathcal{A}$. Similarly,
$\mathcal{P}([0,t])=\\{P_{i}([0,t])|P_{i}(.)\in\mathcal{P}\\}$ represents the
set of all entities’ time series associated with the time interval $[0,t]$.

As previously discussed, a crucial aspect of pair trading is to discern the
temporal correlation among entities based on their historical and feature
information. More formally, if we denote $S(P_{i}(t),P_{j}(t))$, e.g., the
normalized historical difference (NHD) [1], as a predetermined correlation
measure of two time series and $\gamma$ as a defined threshold, the aim is,
given a set of entities $\mathcal{A}$, their feature vectors $\mathcal{X}$,
and their time series information $\mathcal{P}([0,T])$, to identify pairs of
entities $A_{i},A_{j}$ such that
$S(p_{i}([T,T+\delta]),p_{j}([T,T+\delta]))\geq\gamma$ for some $\delta>0$.

###  3.1 Temporal Graph Construction

Given data from the time interval $[0,T]$, we initially divide this interval
into $K$ smaller intervals, $\Omega=\\{\Delta_{i}\\}_{i=1,...,K}$. We assume
that each time interval $\Delta_{i}$ has the same size of
$|\Delta_{i}|=\delta=\frac{T}{K}$, i.e.,
$\Omega=\\{[0,\delta),[\delta,2\delta),...,[(K-1)\delta,K\delta)\\}$ for
simplification.

Assuming $S(P_{i}(t),P_{j}(t))\mapsto[0,1]$ is a given measure function (for
simplicity, we use the NHD) for the correlation of two time series, and
$\gamma$ is a given threshold, we construct the $k-th$ temporal graph $G_{k}$
associated with the time interval $\Delta_{k}$ as follows.

  * •

We add a vertex $v_{i}$ to the graph for each entity $A_{i}$.

  * •

An edge $e_{ij}(t)$ is added between $v_{i}$ and $v_{j}$ only if
$S(P_{i}(\Delta_{k}),P_{j}(\Delta_{k}))\geq\gamma$.

  * •

We assign $t=(k-\frac{1}{2})\delta$ to the edge above.

The overall procedure of constructing the temporal graph is summarized in
Algorithm 2. We denote $\mathcal{G}=\\{G_{i}\\}_{i=1,...,K}$ as the temporal
graph sequence given by the above construction process. The set of vertices
(entities) in $\mathcal{G}$ is denoted as $\mathcal{V}=\\{1,...,n\\}$.
Additionally, following the practice in dynamic graph learning, we refer to an
edge in $\mathcal{G}$ as an event.

###  3.2 Temporal Graph Neural Networks

Temporal graph neural networks (TGNNs) have shown substantial potency as
neural models, particularly for forecasting temporal graphs [2, 18, 19].
Memory-based TGNNs (MTGNNs), a subset of TGNNs, offer superior outcomes over
memory-less alternatives [2]. One salient feature of MTGNNs is the integration
of a memory module, which serves as a filter, iteratively processing and
refining data from new and historical graph events. As a result, MTGNNs
capture extensive dependencies efficiently and deliver leading performance on
a wide variety of dynamic graph-related tasks [20].

In our study, we employ and adapt MTGNNs as our neural model. Adhering to the
settings in [20][21], we structure our MTGNN with an encoder-decoder
configuration. The encoder of MTGNN comprises three primary modules: msg
(message), mem (memory), and emb (embedding) modules. The output from the
embedding module is fed into a decoder (in this case, a simple two-layer MLP)
to execute the inference task. The details of each module in the MTGNN are
presented in the following sections.

####  3.2.1 Encoder

The encoder of MTGNN houses three modules (message, memory, and embedding),
introduced separately for clarity.

Message. For each event (edge) that involves node $i$, a message is computed
to refresh $i$’s memory. Given an interaction event $e_{ij}(t)$ between source
node $i$ and target node $j$ at time $t$, two messages can be formulated:

$\begin{split}m_{i}(t)=\mathrm{msg}(s_{i}(t^{-}),s_{j}(t^{-}),e_{ij}(t),\psi(t-t_{i}^{\prime})),\\\
m_{j}(t)=\mathrm{msg}(s_{j}(t^{-}),s_{i}(t^{-}),e_{ij}(t),\psi(t-t_{j}^{\prime})).\end{split}$
(1)

Here, $s_{i}(t^{-})$ and $s_{j}(t^{-})$ are the memory states of node $i$ and
$j$ just before the event at time $t$, $\mathrm{msg}(.)$ is the message
function, and $t^{\prime}_{j}$ is the timestamp of the last event involving
node $j$. The time encoding function $\psi(.)$ [22] maps the time interval
into a $d$-dimension vector. For simplification, we use the popular identity
message function that returns the concatenation of the input vectors [20][21].

Memory. The memory of a node is updated for each event involving that node:

$\begin{split}s_{i}(t)=\mathrm{mem}(m_{i}(t),s_{i}(t^{-})).\end{split}$ (2)

In the case of interaction events that involve two nodes $i$ and $j$, the
memories of both nodes are updated post the event. Here, $\mathrm{mem}$ is a
learnable memory update function. In practice, we use the gated recurrent unit
(GRU) [23]. The memory module functions as a filtering mechanism, learning and
distilling information from both new and historical temporal graph data
iteratively. Thus, it effectively captures long-range dependencies.

Algorithm 1 Training Procedure of MTRGL

Input: $\mathcal{G}=\\{G_{1},...,G_{k}\\}$ // {The temporal graph sequence}

Input: $T$ // {Number of epoches}

Initialization: $S_{0}\leftarrow 0$ {Initialize memory vectors to be zero}

Initialization: $\theta_{0}^{(0)}\leftarrow 0$ {Initialize model parameter}

for t=1 to T do

for $G_{i}\in G_{2},...,G_{K}$ do

$B_{i}^{+}\leftarrow\mathcal{E}_{i}$, $B_{i}^{-}\leftarrow\text{Sample
negative events}$

$\bar{B}_{i}=B_{i}^{-}\bigcup B_{i}^{+}$ {Use the temporal batch from last
iteration to update the memory and embedding}

$\bar{B}_{i-1}\leftarrow\text{Temporal batch from last iteration}$

$M_{i}=\mathrm{msg}(S_{i-1},\bar{B}_{i-1})$ {Compute messages for events}

$S_{i}=\mathrm{mem}(S_{i-1},M_{i})$ {Update memory the message}

$H_{i}=\mathrm{emb}(S_{i},A_{i}),$ {Compute the embedding}

Compute the loss $\mathcal{L}(H_{i},B_{i})$(given in Eq. (5)) and update the
model parameter with training algorithm (e.g., backpropagation and Adam)

end for

end for

Embedding. The aim of the temporal embedding module is to generate pre-jump
representations $z(t^{-})$ at any time $t$ before the arrival of the next
event. We specifically employ an $L$-layer temporal graph attention network to
gather neighbourhood information.

We first enhance the memory vector of entity $A_{i}$ with the node feature
$z_{i}^{(0)}(t)=s_{i}(t)+X_{i}(t)$. This allows the model to utilize both the
current memory $s_{i}(t)$ and the temporal node features $X_{i}(t)$. Then for
each layer $1\leq l\leq L$, we employ the multi-head attention [22] to
aggregate neighbourhood information as follows:

$\begin{split}z_{i}^{(l)}&=\mathrm{mlp}^{(l)}(z^{(l-1)}||\tilde{z}_{k}^{(l)}),\\\
\tilde{z}_{k}^{(l)}&=\mathrm{mha}^{(l)}(q^{(l)}_{i}(t),K^{(l)}_{i}(t),V^{(l)}_{i}(t)),\\\
q^{(l)}(t)&=z_{i}^{(l-1)}||\psi(0),\\\
K_{i}^{(l)}(t)=V^{(l)}_{i}(t)&=\begin{bmatrix}z_{i}^{(l-1)}||e_{\pi_{i}(1)}(t_{\pi_{i}(1)})||\psi(t-t{\pi_{i}(1)})\\\
...\\\
z_{i}^{(l-1)}||e_{\pi_{i}(N)}(t_{\pi_{i}(N)})||\psi(t-t_{\pi_{i}(N)})\\\
\end{bmatrix}.\end{split}$ (3)

Here, $||$ denotes vector concatenation, $\mathrm{mlp}(.)$ are one-layer MLP
with hidden/output dimensions of $d$, and $\mathrm{mha}(.)$ are multi-head
attention functions with queries $q(.)$, keys $K(.)$ and values $V(.)$. We
limit the receptive fields to the most recent $N$ events for each node $i$,
represented by
$\pi_{i}=\\{e_{\pi_{i}(1)}(t_{\pi_{i}(1)}),e_{\pi_{i}(1)}(t_{\pi_{i}(2)}),...,e_{\pi_{i}(N)}(t_{\pi_{i}(N)})\\}$.
In this context, $\pi$ is a permutation, and $\pi_{i}(.)$ are temporal
neighbours of node $i$. The subscript $i$ in $\pi_{i}(.)$ is omitted,
referring to an event from $i$ to $\pi_{i}(.)$ or vice versa.

####  3.2.2 Decoder

In our study, we simplify the problem of discerning temporal correlation into
the temporal link prediction task. Given $z_{i}(t^{-})$ and $z_{j}(t^{-})$ as
produced by the encoder, the decoder computes the probability of the event
$e_{ij}(t)$ utilizing a two-layer MLP followed by the sigmoid function
$\sigma(.)$,

$\hat{p}_{ij}(t)=\sigma(\mathrm{MLP}(z_{i}(t^{-})||z_{j}(t^{-}))).$ (4)

###  3.3 Model Training and Inference

The training process for the model leverages binary cross-entropy as the loss
function and contrastive learning:

$\mathcal{L}=-\sum_{e_{ij}(t)\in\mathcal{E}}[\log\hat{p}_{ij}(t)+\log(1-\hat{p}_{ik}(t))],$
(5)

where a negative destination node is randomly sampled as $k$, and
$\hat{p}_{ik}(t)$ is calculated in a similar manner using $z_{i}(t^{-})$ and
$z_{k}(t^{-})$. This facilitates the neural model in learning contrastive
signals derived from the generated graphs.

Moreover, we implement a temporal batch training method for the sake of
efficiency. Successive events are segmented into a temporal batch, within
which events are processed concurrently. For nodes associated with multiple
events within a batch, we resort to the simplest recent message aggregator
that retains only the most recent message for each node in the batch,
following the practice from [20]. In order to avoid information leakage
(preventing the use of a batch’s information to predict its own events), we
employ a lag-one scheme in which temporal batch $B_{i-1}$ is utilized to
refresh the memory state and create embeddings for predicting $B_{i}$. The
comprehensive training procedure is provided in Algorithm 2 and Figure 2(b)
gives a graphical illustration.

Algorithm 2 Multi-modal Temporal Relation Graph Construction

Input: $\mathcal{A}=\\{A_{1},...,A_{n}\\}$ // { set of entities}

Input: $K$ // {Number of the time interval to be constructed}

Input: $\gamma$ // {The threshold for correlation }

Input: $S(.,.)$ // {a given covariance measure }

Initialization: $\mathcal{G}\leftarrow[.]$ {Initialize an empty list for
storing graph}

Divide $\mathcal{A}$ into $K$ intervals $\\{\Delta_{1},..,\Delta_{K}\\}$.

for $\Delta_{i}\in\\{\Delta_{1},..,\Delta_{K}$ } do

Initialize $G_{i}$.

Create a vertex $v_{j}$ for each $A_{j}\in\mathcal{A}$ with feature $X_{j}$.

Create an edge $e_{ij}(t)$ between $v_{i}$ and $v_{j}$ if and only if
$S(P_{i}(\Delta_{k}),P_{j}(\Delta_{k}))\geq\gamma.$

Assign $t=(k-\frac{1}{2})\delta$ and append $G_{i}$ to $\mathcal{G}$

end for

return $\mathcal{G}$.

##  4 Experiments

This section presents an empirical evaluation of the proposed method, with a
focus on two key research questions. Q1: How does the performance of our
approach compare with existing baselines? Q2: Does incorporating multi-modal
information enhance the solution?

####  4.0.1 Data and Preprocessing

The evaluation is carried out using publicly available financial data from
Yahoo Finance and Naver Finance. The study primarily targeted three indices:
the Korea Composite Stock Price Index (KOSPI), the Standard & Poor’s 500 Index
(S&P 500), and Heng Seng Index (HSI), thus offering a comparative study across
different markets. We intentionally focus on the pre-pandemic period,
specifically the daily closing price from 2015 to 2019, to circumvent the
”black swan” effect induced by the COVID-19 pandemic. Information such as
market capitalization, and sector information are used as elements in the
feature vector of each entity.

We divide data into training, validation and testing in a chronological
manner. The data is split in a $60/20/20$ (%) ratio. The first three years
serve as the training data, the following year is used for validation, and the
final year is used for testing. This arrangement ensures an evaluation setup
that realistically mimics the chronological nature of financial data and
allows for a rigorous evaluation

####  4.0.2 Benchmark Methods

Our approach is compared with four different methods which act as our
experimental baselines. The first one utilizes a simple 2-layer MLP, depending
solely on static entity features. The second method leverages Long Short-Term
Memory (LSTM) [24], which uses historical price data to identify correlations.
The third method projects the time series into a time-frequency domain,
utilizing a Convolutional Neural Network (CNN) for learning [25]. We also
include a traditional statistical method based on cointegration (COINT) [24].
All machine learning methods (including ours) are trained until convergence or
a maximum of $50$ epochs, using the Adam optimizer with a learning rate of
0.01 and a weight decay of $0.0001$.

####  4.0.3 Evaluation Measures

Our performance evaluation relies on two key measures: average precision (AP)
and mean absolute percentage error (MAPE). AP is a standard measure for models
predicting categorical outcomes, indicating the proportion of correctly
identified correlated pairs in our case. Conversely, MAPE measures the
precision in estimating the actual correlation value, giving us insight into
how well the model predicts the correlation degree between identified pairs.

###  4.1 Results

Effectiveness. The effectiveness of our proposed method, MTRGL, is evaluated
first. We selected a time interval of one month for both training and
inference, enabling the capture of short-term volatility and long-term trends.
The results of this experiment are presented in Table 1 and the values are
averaged from five independent trials. Bolded entries represent the best
performance. $\uparrow$ indicates that higher values are better, and
$\downarrow$ denotes that lower values are better. As Table 1 shows, MTRGL
significantly outperforms the existing baselines in predicting future
correlations among entities.

Ablation Study. A ablation study is conducted to validate the design decisions
of MTRGL, examining the use of feature information and the structural
information from the constructed graph. The impact of feature information is
summarized in Table 2. In MTRGL-one-hot, feature vectors are replaced with a
unique one-hot identifier for each entity. The results reveal a significant
performance drop when feature information is removed. The impact of structural
information is summarized in Table 3. In MTRGL-edgeless, the memory module is
removed, and each entity’s memory vector is used as input for the decoder.
This eliminates the concept of temporal neighbourhood and hence the structural
information. As Table 3 illustrates, removing structural information also
leads to a performance reduction in MTRGL.

Table 1: Performance comparison of MTRGL with baselines.

| KOSPI | S&P 500 | HSI
---|---|---|---
|  AP(%) $\uparrow$ |  MAPE $\downarrow$ |  AP(%) $\uparrow$ |  MAPE $\downarrow$ |  AP(%) $\uparrow$ |  MAPE $\downarrow$
COINT | 55.8 ± 0.5 | 43.6 ± 1.2 | 56.6 ± 0.9 | 40.2 ± 1.3 | 51.4 ± 0.7 | 32.8 ± 1.2
MLP | 48.2 ± 0.4 | 45.2 ± 1.6 | 46.8 ± 0.3 | 43.2 ± 1.4 | 44.2 ± 0.4 | 34.4 ± 1.3
CNN | 62.7 ± 0.4 | 32.8 ± 1.3 | 64.8 ± 0.3 | 37.2 ± 1.5 | 64.2 ± 0.8 | 25.8 ± 1.4
LSTM | 61.4 ± 0.3 | 30.6 ± 1.5 | 65.6 ± 0.5 | 34.8 ± 1.3 | 61.5 ± 0.5 | 23.3 ± 1.2
MTRGL (ours) | 72.8 ± 0.4 | 24.2 ± 1.0 | 74.2 ± 0.4 | 27.8 ± 1.3 | 69.8 ± 0.7 | 16.8 ± 1.4

Table 2: Performance of MTRGL w./w.o feature information

| KOSPI | S&P 500 | HSI
---|---|---|---
|  AP(%) $\uparrow$ |  MAPE $\downarrow$ |  AP(%) $\uparrow$ |  MAPE $\downarrow$ |  AP(%) $\uparrow$ |  MAPE $\downarrow$
MTRGL-one-hot | 63.2 ± 0.9 | 32.2 ± 1.2 | 64.3 ± 0.2 | 34.6 ± 1.5 | 60.8 ± 0.5 | 21.6 ± 1.9
MTRGL | 72.8 ± 0.4 | 24.2 ± 1.0 | 74.2 ± 0.4 | 27.8 ± 1.3 | 69.8 ± 0.7 | 16.8 ± 1.4

Table 3: Performance of MTRGL w./w.o structural information

| KOSPI | S&P 500 | HSI
---|---|---|---
|  AP(%) $\uparrow$ |  MAPE $\downarrow$ |  AP(%) $\uparrow$ |  MAPE $\downarrow$ |  AP(%) $\uparrow$ |  MAPE $\downarrow$
MTRGL-edgeless | 57.2 ± 1.1 | 36.2 ± 2.4 | 58.0 ± 0.9 | 38.2 ± 2.1 | 53.6 ± 0.8 | 26.2 ± 1.4
MTRGL | 72.8 ± 0.4 | 24.2 ± 1.0 | 74.2 ± 0.4 | 27.8 ± 1.3 | 69.8 ± 0.7 | 16.8 ± 1.4

##  5 CONCLUSION

In this paper, we have embarked on an exploration of the application of deep
learning in the realm of pair trading, a well-regarded quantitative investment
strategy. This journey has led to the creation of a unique approach, MTRGL,
explicitly designed to amalgamate descriptive and time series data, thereby
optimizing the process of discerning temporal correlations in pair trading.
Our empirical evidence shows that MTRGL is highly effective in automatically
identifying correlated pairs, surpassing the performance of traditional
baselines that rely exclusively on either descriptive or time series data.

Moreover, the novel integration of multi-modal information in our approach
extends beyond the scope of this study and pair trading. Its potential impact
is significant for other quantitative finance-related problems, and we look
forward to seeing how this innovation could transform these areas in the
future.

Acknowledgements. This study was funded by the Yangtze River Delta Joint
Research and Innovation Community (Grant No.2022CSJGG1200), National Natural
Science Foundation of China (Grant No. 42302326) and the Central University
Basic Research Business Special Fund Assistance (No. JZ2023HGTA0178).

## References

  * [1] Evan Gatev, William N Goetzmann, and K Geert Rouwenhorst,  “Pairs trading: Performance of a relative-value arbitrage rule,”  The Review of Financial Studies, vol. 19, no. 3, pp. 797–827, 2006\.
  * [2] Farimah Poursafaei, Shenyang Huang, Kellin Pelrine, , and Reihaneh Rabbany,  “Towards better evaluation for dynamic link prediction,”  in Neural Information Processing Systems (NeurIPS) Datasets and Benchmarks, 2022.
  * [3] Kwok Tai Chui, Brij B Gupta, Hao Ran Chi, Varsha Arya, Wadee Alhalabi, Miguel Torres Ruiz, and Chien-Wen Shen,  “Transfer learning-based multi-scale denoising convolutional neural network for prostate cancer detection,”  Cancers, vol. 14, no. 15, pp. 3687, 2022.
  * [4] Quang Thanh Bui, Van Long Vu, Xuan Linh Nguyen, and Van Manh Phan,  “Combination of uav images and dsm for land cover classification using convolutional neural network,”  VNU Journal of Science: Earth and Environmental Sciences, vol. 38, no. 1, 2022.
  * [5] Jana Wäldchen and Patrick Mäder,  “Machine learning for image based species identification,”  Methods in Ecology and Evolution, vol. 9, no. 11, pp. 2216–2225, 2018.
  * [6] Ka Wing Wan, Chun Hoi Wong, Ho Fung Ip, Dejian Fan, Pak Leung Yuen, Hoi Ying Fong, and Michael Ying,  “Evaluation of the performance of traditional machine learning algorithms, convolutional neural network and automl vision in ultrasound breast lesions classification: A comparative study,”  Quantitative imaging in medicine and surgery, vol. 11, no. 4, pp. 1381, 2021.
  * [7] Wahab Khan, Ali Daud, Jamal A Nasir, and Tehmina Amjad,  “A survey on the state-of-the-art machine learning models in the context of nlp,”  Kuwait journal of Science, vol. 43, no. 4, 2016.
  * [8] Dan Ofer, Nadav Brandes, and Michal Linial,  “The language of proteins: Nlp, machine learning & protein sequences,”  Computational and Structural Biotechnology Journal, vol. 19, pp. 1750–1758, 2021.
  * [9] Joon-Woo Lee, HyungChul Kang, Yongwoo Lee, Woosuk Choi, Jieun Eom, Maxim Deryabin, Eunsang Lee, Junghyun Lee, Donghoon Yoo, Young-Sik Kim, et al.,  “Privacy-preserving machine learning with fully homomorphic encryption for deep neural network,”  IEEE Access, vol. 10, pp. 30039–30054, 2022.
  * [10] Feng Xia, Ke Sun, Shuo Yu, Abdul Aziz, Liangtian Wan, Shirui Pan, and Huan Liu,  “Graph learning: A survey,”  IEEE Transactions on Artificial Intelligence, vol. 2, no. 2, pp. 109–127, 2021.
  * [11] Yanjie Gao, Xianyu Gu, Hongyu Zhang, Haoxiang Lin, and Mao Yang,  “Runtime performance prediction for deep learning models with graph neural network,”  in 2023 IEEE/ACM 45th International Conference on Software Engineering: Software Engineering in Practice (ICSE-SEIP). IEEE, 2023, pp. 368–380.
  * [12] Junwei Su, Difan Zou, Zijun Zhang, and Chuan Wu,  “Towards robust graph incremental learning on evolving graphs,”  in International Conference on Machine Learning. PMLR, 2023, pp. 32728–32748.
  * [13] Ahmet Murat Ozbayoglu, Mehmet Ugur Gudelek, and Omer Berat Sezer,  “Deep learning for financial applications: A survey,”  Applied Soft Computing, vol. 93, pp. 106384, 2020.
  * [14] Alejandro Lopez-Lira and Yuehua Tang,  “Can chatgpt forecast stock price movements? return predictability and large language models,”  arXiv preprint arXiv:2304.07619, 2023.
  * [15] Riaz Qureshi, Daniel Shaughnessy, Kayden AR Gill, Karen A Robinson, Tianjing Li, and Eitan Agai,  “Are chatgpt and large language models “the answer” to bringing us closer to systematic review automation?,”  Systematic Reviews, vol. 12, no. 1, pp. 72, 2023.
  * [16] Zihan Chen, Lei Nico Zheng, Cheng Lu, Jialu Yuan, and Di Zhu,  “Chatgpt informed graph neural network for stock movement prediction,”  arXiv preprint arXiv:2306.03763, 2023.
  * [17] Brady D Lund, Ting Wang, Nishith Reddy Mannuru, Bing Nie, Somipam Shimray, and Ziang Wang,  “Chatgpt and a new academic reality: Artificial intelligence-written research papers and the ethics of the large language models in scholarly publishing,”  Journal of the Association for Information Science and Technology, vol. 74, no. 5, pp. 570–581, 2023.
  * [18] Emanuele Rossi, Ben Chamberlain, Fabrizio Frasca, Davide Eynard, Federico Monti, and Michael Bronstein,  “Temporal Graph Networks for Deep Learning on Dynamic Graphs,”  in Proceedings of International Conference on Learning Representations, 2021.
  * [19] Srijan Kumar, Xikun Zhang, and Jure Leskovec,  “Predicting dynamic embedding trajectory in temporal interaction networks,”  in Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, 2019, pp. 1269–1278.
  * [20] Yao Zhang, Yun Xiong, Yongxiang Liao, Yiheng Sun, Yucheng Jin, Xuehao Zheng, and Yangyong Zhu,  “Tiger: Temporal interaction graph embedding with restarts,” 2023.
  * [21] Emanuele Rossi, Ben Chamberlain, Fabrizio Frasca, Davide Eynard, Federico Monti, and Michael Bronstein,  “Temporal graph networks for deep learning on dynamic graphs,”  arXiv preprint arXiv:2006.10637, 2020.
  * [22] Da Xu, Chuanwei Ruan, Evren Korpeoglu, Sushant Kumar, and Kannan Achan,  “Inductive representation learning on temporal graphs,”  arXiv preprint arXiv:2002.07962, 2020.
  * [23] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio,  “Learning phrase representations using rnn encoder-decoder for statistical machine translation,”  arXiv preprint arXiv:1406.1078, 2014.
  * [24] Simão Moraes Sarmento and Nuno Horta,  “Enhancing a pairs trading strategy with the application of machine learning,”  Expert Systems with Applications, vol. 158, pp. 113490, 2020.
  * [25] Bendong Zhao, Huanzhang Lu, Shangfeng Chen, Junliang Liu, and Dongya Wu,  “Convolutional neural networks for time series classification,”  Journal of Systems Engineering and Electronics, vol. 28, no. 1, pp. 162–169, 2017.