11institutetext: Faculty of Aerospace Engineering, Delft University of
Technology, Delft, The Netherlands
11email<EMAIL_ADDRESS>22institutetext: Centro de Astrobiología (CSIC-
INTA), Ctra. de Ajalvir, km 4, Torrejón de Ardoz, 28850 Madrid, Spain
33institutetext: Leiden Observatory, Leiden University, P.O. Box 9513, NL 2300
RA Leiden, The Netherlands 44institutetext: Department of Physics, National
Central University, Jhongli City, Taoyuan County 32054, Taiwan
55institutetext: Max Planck Institute for Extraterrestrial Physics, Postfach
1312, 85741 Garching, Germany

# Characterization of carbon dioxide on Ganymede and Europa supported by
experiments: Effects of temperature, porosity, and mixing with water

L. Schiltz 11 B. Escribano 22 G. M. Muñoz Caro 22 S. Cazaux 1133 C. del Burgo
Olivares 22 H. Carrascosa 22
I. Boszhuizen 11 C. González Díaz 22 Y.-J. Chen 44 B. M. Giuliano 55 P.
Caselli 55

(Received XXXXXX; accepted XXXX)

###### Abstract

Context. The surfaces of icy moons are primarily composed of water ice that
can be mixed with other compounds, such as carbon dioxide. The carbon dioxide
(CO${}_{2})$ stretching fundamental band observed on Europa and Ganymede
appears to be a combination of several bands that are shifting location from
one moon to another.

Aims. We investigate the cause of the observed shift in the CO2 stretching
absorption band experimentally. We also explore the spectral behaviour of CO2
ice by varying the temperature and concentration.

Methods. We analyzed pure CO2 ice and ice mixtures deposited at 10 K under
ultra-high vacuum conditions using Fourier-transform infrared (FTIR)
spectroscopy and temperature programmed desorption (TPD) experiments.
Laboratory ice spectra were compared to JWST observation of Europa’s and
Ganymede’s leading hemispheres. The simulated IR spectra were calculated using
density functional theory (DFT) methods, exploring the effect of porosity in
CO2 ice.

Results. Pure CO2 and CO2-water ice show distinct spectral changes and
desorption behaviours at different temperatures, revealing intricate CO2 and
H2O interactions. The number of discernible peaks increases from two in pure
CO2 to three in CO2-water mixtures.

Conclusions. The different CO2 bands were assigned to $\tilde{\nu}_{3,1}$
(2351 cm-1, 4.25 $\mu$m) caused by CO2 dangling bonds (CO2 found in pores or
cracks) and $\tilde{\nu}_{3,2}$ (2345 cm-1, 4.26 $\mu$m) due to CO2 segregated
in water ice, whereas $\tilde{\nu}_{3,3}$ (2341 cm-1, 4.27 $\mu$m) is due to
CO2 molecules embedded in water ice. The JWST NIRSpec CO2 spectra for Ganymede
and for Europa can be fitted with two Gaussians attributed to
$\tilde{\nu}_{3,1}$ and $\tilde{\nu}_{3,3}$. For Europa, $\tilde{\nu}_{3,1}$
is located at lower wavelengths due to a lower temperature. The Ganymede data
reveal latitudinal variations in CO2 bands, with $\tilde{\nu}_{3,3}$
dominating in the pole and $\tilde{\nu}_{3,1}$ prevalent in other regions.
This shows that CO2 is embedded in water ice at the poles and it is present in
pores or cracks in other regions. Ganymede longitudinal spectra reveal an
increase of the CO2 $\tilde{\nu}_{3,1}$ band throughout the day, possibly due
to ice cracks or pores caused by large temperature fluctuations.

###### Key Words.:

Carbon dioxide ices – methods: laboratory – JWST observations – infrared
spectroscopy – icy moons

##  1 Introduction

Mixed ices of CO2 and water are known to be present in the frozen nuclei of
comets (Crovisier, 2006b, a) in various satellites of our Solar System (e.g.
Grundy et al., 2003; Buratti et al., 2005) and as components of dust
interstellar particles (e.g Dartois et al., 2005; Draine, 2003). On Ganymede,
the CO2 stretching fundamental band, which is ordinarily at 2341 cm-1 (4.27
$\mu$m) (Falk, 1987), has been observed with the Near Infrared Mapping
Spectrometer (NIMS) aboard the Galileo spacecraft (McCord et al., 1998). The
location of this absorption band is displaced slightly to a shorter wavelength
(2345 cm-1, 4.26 $\mu$m) indicating that the CO2 is bound, or trapped, in a
host. Recent James Webb Space Telescope (JWST) observations (Bockelée-Morvan
et al., 2024) have shown variations in the latitude and longitude of the CO2
bands on Ganymede’s surface. In the boreal region of the leading hemisphere,
the CO2 band is dominated by the 2341 cm-1 (4.27 $\mu$m) band, consistent with
CO2 trapped in amorphous water ice, while at equatorial latitudes (and
especially on dark terrains) the observed band is broader and located around
2345 cm-1 (4.26 $\mu$m), suggesting CO2 adsorbed on non-icy materials, such as
minerals or salts. On Europa, CO2 was detected with the Galileo/NIMS
instrument at 2351 cm-1 (4.25 $\mu$m) and mostly located on the anti-Jovian
and trailing sides (Hansen & McCord, 2008); however, observations in the near-
infrared (NIR) were not able to confirm the presence of CO2 (Mishra et al.,
2021). Recent JWST observations by Villanueva et al. (2023) have highlighted
the presence of CO2 on Europa’s surface, with bands located at 2351 cm-1 (4.25
$\mu$m) and 2341 cm-1 (4.27 $\mu$m). These observations suggest that CO2 is
mixed with other compounds and that carbon is sourced from within Europa,
probably from the liquid sub-surface ocean.

Prompted by these findings, solid H2O/CO2 ice mixtures, as laboratory
analogues of astrophysical objects, have been studied for many years, mainly
by infrared spectroscopy and mass spectrometry, in controlled warm-up
experiments (i.e. Sandford & Allamandola, 1990; Ehrenfreund et al., 1997;
Bernstein et al., 2005; Kumi et al., 2006; Malyk et al., 2007). From these
studies, the interaction between CO2 and H2O on a molecular level has been
shown to cause significant changes in the position and profile of CO2 peaks in
the IR. However, a concise answer to why the CO2 stretching fundamental band
on icy moons, like Europa and Ganymede, shows frequency shifts remains
elusive.

In the present work, we attempt to find an explanation for the different CO2
bands to characterise the carbon dioxide on Europa and Ganymede. This is
achieved by investigating the influence of temperature on the spectral
properties of laboratory CO2 ice, and in particular the behaviour of the CO2
stretching mode, which (in its pure form) peaks around 2344 cm-1 (4.27
$\mu$m). We note that CO2 is first studied in its pure form and afterwards co-
deposited with water so that the H2O:CO2 deposition ratio is varied. In
particular, we considered low (¡5$\%$) and high (¿25$\%$) carbon dioxide
concentrations in the ice. Experiments were carried out at temperatures
ranging between 10 K and 160 K and ultra high vacuum conditions. This research
uses a combination of the following two techniques: temperature programmed
desorption (TPD) and Fourier-transform infrared spectroscopy (FTIR) in
transmittance, across a variety of experimental conditions. Additionally,
Gaussian deconvolution was applied to the spectra obtained experimentally.
This technique was first applied to pure CO2 ice and afterwards to the ice
mixture H2O:CO2. It was found that deconvolution of the asymmetric stretching
band $\nu_{3}$ requires two Gaussians for pure CO2 and three when mixed with
water. The evolution of the integrated absorbance areas, band shape and
position as a function of temperature can be compared to the TPD curve in
order to assign the bands to different molecular interactions, thus
characterising the type of CO2 ice. The synergy between laboratory studies and
data gathered by JWST enriches our comprehensive approach to understanding CO2
behaviour on icy worlds, enabling a thorough exploration of the fundamental
properties of CO2 ice in various conditions.

The structure of this paper is organized as follows. In section 2, we outline
the experimental setup and the determination of the gas mixture ratio. In
section 3, we focus on the spectroscopic behaviour of pure CO2 ice,
elucidating the shape and location of the $\nu_{3}$ stretching band, alongside
the TPD curve. In section 4 and section 5, we describe the acquired laboratory
spectra for low and high CO2 concentrations in water ice as a function of
temperature, accompanied by explanations of integrated absorbance areas and
TPD curves. In section 6 and section 7, we discuss the results and assign the
IR vibrational modes from experiments and simulations. Finally, section 8
delves into the application of findings to icy moons, specifically Europa and
Ganymede, involving a comparison of new experimental results with recent
observations.

##  2 Experimental method

###  2.1 ISAC set-up

The experiments reported in this work have been performed with the
Interstellar Astrochemistry Chamber (ISAC), described in more detail in Muñoz
Caro et al. (2010). ISAC is an ultra-high vacuum (UHV) chamber with a base
pressure of 4 $\times$ 10-11 mbar, designed to simulate the conditions present
in the interstellar medium (ISM), regarding temperature, pressure and
ultraviolet (UV) radiation field. A closed-cycle He cryostat allows for the
tip of the cold finger to cool down to 8 K, where a sample holder with a
potassium bromide (KBr) window acting as the substrate for ice deposition is
located. A schematic representation is shown in Figure 1. A Lakeshore
temperature controller 331 with 0.1 K accuracy was used. The gas line, which
permits the introduction of gas species with regulated compositions, is linked
to the main chamber through a leak valve. The identification of the species
present in the chamber is done via a quadrupole mass spectrometer (QMS,
Pfeiffer Vacuum, Prisma QMS 200). The valve is opened during deposition, and
the gas is delivered to the cold substrate through a deposition tube. This
tube’s end is roughly placed at a distance of 3 cm from the substrate. FTIR in
transmittance mode using a Bruker Vertex 70 at a working spectral resolution
of 2 cm-1 is used to record infrared spectra during deposition, and later
during the TPD using a heating ramp in K/min. ISAC uses laser interferometry
at 632.8 nm to assess changes in ice thickness (González Díaz et al. (2022)).
A 5.0 mW He-Ne red laser and a Silicon Photodiode Power Sensor to measure the
optical power of the laser light (model S120C) are mounted at an angle of 6∘,
as seen in Figure 1.

In general, ices were grown at 10 K by opening the leak valve of gas line 2 in
Figure 1, feeding a mixture of H2O and CO2 into the main chamber, at a
deposition pressure of 2$\times$10-7 mbar. To investigate the influence of
deposition temperature and pressure on the spectral properties of CO2 ice (in
particular, the behaviour of the CO2 stretching mode), some experiments were
also carried out at 70 K and higher deposition pressure. The experiments
presented in Table 1 were repeated for reproducibility or for suspected
contamination. During the deposition of the gas mixture to form the ice,
infrared spectra at 45∘ incidence of the beam relative to the substrate were
collected every 300 s. Warming up the samples after deposition was done at 0.2
K min-1 until a temperature of 210 K was attained, ensuring that both the CO2
and H2O were entirely desorbed. IR measurements were taken every 300 s again,
yielding a spectrum every 1 K temperature increment.

Figure 1: ISAC cross-section at the ice sample level. The measuring devices
and sensors, including the laser interferometry, are displayed. The FTIR
source is positioned on one side and the FTIR detector on the opposite one.
The UV spectrometer is placed directly across from the vacuum-UV lamp, see
González Díaz et al. (2022) for more details.

The key parameters measured during the experiments using various devices
include: the laser intensity to measure the ice thickness (González Díaz et
al., 2022), pressure inside the chamber, the ion currents of the relevant
molecules (via QMS), and temperature. The temperature of the ice sample was
measured using a silicon diode sensor connected to the sample holder and
placed directly above the ice sample. A Bayard-Alpert gauge positioned about
23 cm below the plane where deposition occurs was used to measure the pressure
in the main chamber of ISAC. Laser interference was monitored continuously,
resulting in an interference pattern during ice deposition and a new
interference pattern during TPD. The ice column density was calculated using
IR spectroscopy from the areas calculated by integration of the absorption
bands. This is covered in more detail in the following section.

###  2.2 Estimation of the ice mixture ratio

The integration of the infrared absorption band yields the column density $N$
of the ice layer accreted on the cold substrate in molecules cm-2 with the
following formula:

$N=\int_{band}\frac{\tau_{\nu}d\nu}{A},$ (1)

with $A$ the band strength in cm molecule-1, $\tau_{\nu}$ the optical depth of
the band and d$\nu$ the wavenumber differential in cm-1. The integrated
absorbance is equal to 0.43 $\times$ $\tau$, where $\tau$ is the integrated
optical depth of the band. The ice mixture ratio is calculated by dividing the
column density of the carbon dioxide anti-symmetric stretching band to the
water stretching band. The IR band strength of the CO2 str. feature in the
H2O:CO2 = 24:1 ice mixture is about 94% with respect to the pure CO2 ice (from
Table 3 of Gerakines et al. (1995)), and for higher CO2 concentrations this
difference is expected to be smaller. We therefore used the band strengths of
pure ices at 10 K: A(H2O)=2.0$\times$10-16 cm molecule-1 (Hagen et al. (1981))
and A(CO2)=7.6$\times$10-17 cm molecule-1 (Bouilloud et al. (2015)). To
compute the integrated intensities of the H2O stretching peaks, the measured
peak positions and the integration limits in cm-1($\mu$m) are taken between
3000 (3.33) and 3600 (2.78). The same method is applied to the CO2 anti-
symmetric stretching mode in the 2330 (4.29) to 2360 (4.24) spectral range. A
table summarizing the experimental parameters can be found in Table 1, where
$p_{dep}$ is the total pressure in the ISAC chamber during deposition of the
ice, $T_{dep}$ is the temperature at which the ice is deposited, $d_{ice}$ is
the ice thickness, and $dT/dt$ is the heating rate in Kelvin per minute.

Table 1: Experimental parameters Composition |  | $p_{dep}$
---
(mbar)
| $T_{dep}$
---
(K)
| $d_{ice}$
---
($\mu$m)
| $\frac{dT}{dt}$
---
(K/min)
Pure CO2 | 2$\times$10-7 | 10 | 0.62 | 0.5
H2O:CO2=1:0.04 | 2$\times$10-7 | 10 | 0.54 | 0.2
H2O:CO2=1:0.25 | 2$\times$10-7 | 10 | 0.67 | 0.2

##  3 Pure CO2 ice

Pure CO2 ice was deposited at 10 K and warmed up until 100 K. The study of the
CO2 ice stretching band, commonly known as $\nu_{3}$, as well as the TPD are
reported in Sections 3.1 and 3.2, respectively.

###  3.1 $\nu_{3}$ stretching band shape and position

Figure 2 shows the CO2 ($\nu_{3}$) asymmetric stretching fundamental for 12CO2
when deposited at 10 K, at a pressure of 2$\times$10-7 mbar and warmed-up with
a rate of 0.5 K/min up to 90 K. All the FTIR spectra have a resolution of 2
cm-1 and are shifted vertically for clarity. The 12CO2 asymmetric stretching
fundamental ($\nu_{3}$) is located at $\sim$2345 cm-1 (4.26 µm), and is
redshifted from the gas-phase value located at 2348 cm-1 (4.26 $\mu$m) due to
interactions with the surrounding matrix environment (e.g. Isokoski et al.
(2013)). The profile is asymmetric with a prominent blue shoulder around 2350
cm-1 (4.25 $\mu$m). When the CO2 ice is warmed up, the main peak at 2345 cm-1
(4.26 $\mu$m) shifts to lower wavenumbers. In addition, the peak intensity
increases with increasing temperature as a result of the decrease in
bandwidth. Figure 3 illustrates the effect of temperature on the band’s shape
and position.

In general, each spectrum can be deconvolved into two Gaussian distributions,
which we refer to as $\tilde{\nu}_{3,1}$ at 2351 cm-1 (4.25 $\mu$m) and
$\tilde{\nu}_{3,2}$ at 2345 cm-1 (4.26 $\mu$m), both for 10 K (see Table 2).
The figure clearly shows the redshift of both Gaussians during warm-up as well
as the decrease in FWHM for $\tilde{\nu}_{3,2}$ (2345 cm-1, 4.26 $\mu$m),
going from 5.7 cm-1 at 10 K to 3.8 cm-1 at 80 K.

Figure 2: Spectra over the 2380-2320 cm-1 (4.20-4.31 $\mu$m) range of a pure CO2 ice deposited at 10 K and 2$\times$10-7 mbar, later warmed up to 90 K. All spectra are measured at 2 cm-1 resolution, and at temperatures indicated in each graph. Spectra are shifted vertically for clarity. Figure 3: Evolution of both Gaussian distributions during warm-up from 10 K to 80 K for a pure CO2 ice. The arrows indicate the direction of evolution. Thicker lines indicate the initial and final states. Table 2: Band positions ($\tilde{\nu}$) for all three fitted Gaussians as a function of temperature. This is done for a pure CO2 ice, a low CO2 concentration and a high CO2 concentration. | Pure CO2 | Low CO2 | High CO2
---|---|---|---
T (K) |  | $\tilde{\nu}_{3,1}$
---
(cm-1)
| $\tilde{\nu}_{3,2}$
---
(cm-1)
| $\tilde{\nu}_{3,3}$
---
(cm-1)
| $\tilde{\nu}_{3,1}$
---
(cm-1)
| $\tilde{\nu}_{3,2}$
---
(cm-1)
| $\tilde{\nu}_{3,3}$
---
(cm-1)
| $\tilde{\nu}_{3,1}$
---
(cm-1)
| $\tilde{\nu}_{3,2}$
---
(cm-1)
| $\tilde{\nu}_{3,3}$
---
(cm-1)
10 | 2351.3 | 2345.0 | - | 2350.5 | 2345.0 | 2340.8 | 2353.5 | 2344.3 | 2337.3
30 | 2351.0 | 2344.8 | - | 2350.2 | 2344.8 | 2340.8 | 2353.7 | 2344.0 | 2337.0
50 | 2350.0 | 2344.3 | - | 2349.8 | 2344.5 | 2340.7 | 2353.8 | 2344.0 | 2337.1
70 | 2349.4 | 2344.2 | - | 2349.8 | 2344.2 | 2340.6 | 2353.7 | 2344.1 | 2337.2
80 | 2349.0 | 2344.0 | - | 2347.5 | - | 2340.7 | 2351.5 | 2344.1 | 2339.0
90 | - | - | - | 2347.0 | - | 2340.7 | 2348.0 | - | 2339.7
100 | - | - | - | 2347.0 | - | 2340.8 | 2348.3 | - | 2340.0
110 | - | - | - | 2347.0 | - | 2340.7 | 2347.8 | - | 2339.9
130 | - | - | - | 2346.9 | - | 2340.8 | 2347.8 | - | 2339.8
140 | - | - | - | 2347.0 | - | 2340.8 | 2347.8 | - | 2339.9
150 | - | - | - | 2346.8 | - | 2340.7 | 2349.0 | - | 2340.3
160 | - | - | - | 2348.0 | - | 2340.6 | 2349.3 | - | 2340.5
170 | - | - | - | 2348.0 | - | - | 2349.1 | - | -

###  3.2 Thermal desorption of CO2

A TPD curve of a pure CO2 ice layer is shown in Figure 4. The QMS data
corresponding to the molecular mass of CO2 is plotted as the ion current in
Ampere versus increasing temperature. This figure shows one desorption peak at
85 K. It is at this temperature that CO2 thermally desorbs. The desorption
rate, expressed in molecules cm-2 s-1, of CO2 ice is described by the Polanyi-
Wigner equation:

$\frac{dN_{g}(CO_{2})}{dt}=\nu_{i}[N_{s}(CO_{2})]^{i}\text{exp}\left(-\frac{E_{d}(CO_{2})}{T}\right),$
(2)

where $N_{g}(CO_{2})$ is the column density of CO2 molecules that desorb from
the ice surface (cm-2), $\nu_{i}$ is the frequency factor (molecules1-i
cm${}^{-2(1-i)}s^{-1})$ for desorption order, $i$, $N_{s}(CO_{2})$ is the
column density of CO2 molecules on the surface at time $t$, $E_{d}(CO_{2})$ is
the binding energy expressed in K, and $T$ the surface temperature in K. The
TPD data can be fitted using Equation 2 and the relation:

$\frac{dN_{g}(CO_{2})}{dt}=\frac{dT}{dt}\frac{dN}{dT},$ (3)

where $\frac{dT}{dt}$ is the heating rate, 0.5 K min-1 in the reported
experiments. Figure 5 shows the QMS data, column densities during warm-up,
Polanyi-Wigner equation, and calculated zero-order coverage. The coverage is
calculated using 1-$\frac{dN_{g}(CO_{2})}{dt}$. The TPD data were fitted using
the parameter values $\nu_{0}$ = 0.4 $\times$ 1031 molecules cm-2 s-1 and
$E_{d}$ = 2555.78 K for the desorption of CO2 ice. We find that the coverage
values determined by transmittance FTIR (green triangles) and the coverage
curve computed by fitting the TPD curve recorded by QMS (light blue line) are
in good agreement. There are only two remaining triangle datapoints that are
poorly fitted. The issue is that a zero-order fit by definition fails when
there are fewer than 4 $\times$ 1016 molecules per square centimeter of
coverage, as demonstrated by Figure 5.

Figure 4: TPD curve of pure CO2 ice layer deposited at 10 K and heated at 0.5
K/min. The ion current (A) is plotted on a logarithmic scale for a better
appreciation of the curve profile and roughly corresponds to partial pressure
in mbar. Figure 5: Column density of CO2 ice (triangles) calculated from IR
measurements matches the calculated coverage (light blue line). TPD curve of
CO2 from QMS measurements (dark blue) agrees well with the fit corresponding
to the Polanyi-Wigner equation (red line).

##  4 Low CO2 concentration in water ice

A low CO2 content in the ice mixtures, up to 5$\%$, is considered in this
section. The important parameters of the presented experiments can be found in
Table 1. Figure 6 depicts the spectra for a H2O:CO2=1:0.04 ratio deposited at
10 K. The experiment is performed with a deposition pressure of 2$\times$10-7
mbar as reported in Table 1. At 10 K, the 12CO2 asymmetric stretching
fundamental ($\nu_{3}$) can be deconvolved into three Gaussians in these
experiments, $\tilde{\nu}_{3,1}$ around 2351 cm-1 (4.25 $\mu$m),
$\tilde{\nu}_{3,2}$ at 2345 cm-1 (4.26 $\mu$m), and $\tilde{\nu}_{3,3}$ at
2341 cm-1 (4.27 $\mu$m). The integrated absorbance area as a function of
temperature is shown in the bottom panel of Figure 7. Blue dots represent
$\tilde{\nu}_{3,1}$ (2351 cm-1, 4.25 $\mu$m) starting at 10 K, which undergoes
a red shift at elevated temperatures. Yellow triangles represent
$\tilde{\nu}_{3,2}$ (2345 cm-1, 4.26 $\mu$m), which is red shifted by 1 cm-1
at 80 K right before thermal desorption. Green crosses representing
$\tilde{\nu}_{3,3}$ (2341 cm-1, 4.27 $\mu$m) indicate no observable changes
with increasing temperature.

Figure 6: Spectra over the 2380-2320 cm-1 (4.20-4.31 $\mu$m) range of the
H2O:CO2=1:0.04 ice mixture deposited at 10 K and 2$\times$10-7 mbar. All
spectra are measured with 2 cm-1 resolution at the temperatures indicated in
each graph. Spectra are shifted vertically for clarity.

As mentioned before, pure CO2 ice thermally desorbs around 85 K. For the ice
mixture H2O:CO2=1:0.04, the top panel of Figure 7 shows the TPD curves from 10
K to 180 K, with three desorption peaks for carbon dioxide: one around 80 K,
another one around 146 K, and a final one at 162 K, while only one peak
appears for water, also at 162 K. There is a smaller CO2 desorption peak at
150 K which is not reproducible in repeated experiments and could be caused by
a minor rearrangement of the water structure. The following explanation is
given for the three main desorption peaks of carbon dioxide. Firstly, pure
carbon dioxide ice thermally desorbs around 80 K, meaning that only CO2
molecules bonded with other CO2 molecules are desorbing. In second place, the
peak at 140 K is caused by the transformation from cubic to hexagonal water
ice (Martín-Doménech et al., 2014). During this process, carbon dioxide is
pushed out of the water matrix. This peak is known as volcano desorption. The
third peak between 160 and 170 K is the co-desorption of all remaining carbon
dioxide along with water.

It is worth noting that carbon dioxide constantly desorbs between 80 and 146
K, as evidenced by the QMS signal being greater after the initial peak than
before (see blue line in top panel of Figure 7). This continuous desorption is
caused by amorphous water ice rearrangement, which produces a change in the
size and form of the porous structure (pores coalesce) when heated (Cazaux et
al., 2015). Carbon dioxide molecules may diffuse through the pores and
sublimate continuously between 80 and 146 K due to the increasing ice
temperature and favored by this rearrangement. The argument of diffusing
carbon dioxide molecules can be verified by looking at the integrated
absorbance areas of the three Gaussian peaks in the bottom panel of Figure 7.
The data points (triangles) represent $\tilde{\nu}_{3,2}$ (2345 cm-1, 4.26
$\mu$m) and confirm that pure CO2 desorbs thermally at 80 K. During this
desorption of CO2 molecules starting at 75 K, the integrated area of
$\tilde{\nu}_{3,3}$ (2341 cm-1, 4.27 $\mu$m) is increasing significantly,
pointing towards the diffusion of the carbon dioxide molecules through the
water ice structure. Once temperatures of 140 K are reached, the integrated
areas of both $\tilde{\nu}_{3,1}$ (2351 cm-1, 4.25 $\mu$m) and
$\tilde{\nu}_{3,3}$ (2341 cm-1, 4.27 $\mu$m) bands are continuously decreasing
during annealing, confirming the amorphous water ice rearrangements. The cubic
to hexagonal water-ice phase transition that occurs at 146 K further
accelerates the desorption of CO2 molecules. Finally, any remaining CO2 will
co-desorb with water around 162 K.

Figure 7: TPD curves of CO2 and H2O for a H2O:CO2=1:0.04 ice heated at 0.2
K/min (top). Integrated absorbance areas of the three Gaussian distributions
as a function of temperature for the same experiment (bottom).

##  5 High CO2 concentration in water ice

In this section, a high CO2 content in the ice mixtures, above 20$\%$, is
considered. The parameters of the conducted experiments can be found in Table
1. Figure 8 shows the spectra for H2O:CO2=1:0.25 ice deposited at 10 K and
2$\times$10-7 mbar. The broadening of the band is a consequence of the
increased CO2 concentration, which causes particle aggregation and thus
creates additional trapping sites between CO2 and H2O (Ehrenfreund et al.
(1997)). As with the low CO2 concentration, the spectra can be deconvolved
into three Gaussian distributions, the band positions at 10 K are
$\tilde{\nu}_{3,1}$ = 2354 cm-1 (4.25 $\mu$m), $\tilde{\nu}_{3,2}$ = 2344 cm-1
(4.27 $\mu$m), and $\tilde{\nu}_{3,3}$ = 2336 cm-1 (4.28 $\mu$m), at 10 K. Due
to the broadening of the CO2 band, the shoulder on the left is blueshifted by
2 cm-1, compared to the low CO2 concentration spectra at 10 K. In addition,
the shoulder on the right is redshifted by 4 cm-1.

Figure 8: Spectra over the 2380-2320 cm-1 (4.20-4.31 $\mu$m) range of the
H2O:CO2=1:0.25 ice mixture deposited at 10 K and 2$\times$10-7 mbar. All
spectra are measured with 2 cm-1 resolution at the temperatures indicated in
each graph. Spectra are shifted vertically for clarity.

The TPD curves and integrated absorbance area are shown in Figure 9. Again,
three desorption peaks for carbon dioxide can be identified, namely at 82 K,
146 K, and 161 K. The same explanation as for a low CO2 concentration is valid
for this experiment. The small bump in the TPD curve around 30 K is not
reproducible in repeated experiments and might be due to a small amount of air
in the gas line. The TPD curve of carbon dioxide shows the main desorption at
82 K, corresponding to pure CO2 desorption in these experiments, as is
confirmed by the drastic drop of the integrated absorbance area of
$\tilde{\nu}_{3,2}$ at that temperature. The bottom panel of Figure 9 clearly
shows, as in the case of low CO2 concentration, the diffusion through the
water ice structure after the first desorption peak, indicated by the increase
in integrated area for $\tilde{\nu}_{3,3}$ The desorption of CO2 is enhanced
by the phase transition from cubic to hexagonal water ice at 146 K and any
remaining CO2 co-desorbs with water at 161 K.

Figure 9: TPD curves of CO2 and H2O for a H2O:CO2=1:0.25 ice layer heated at
0.2 K/min (top). Integrated absorbance areas of the three Gaussian
distributions as a function of temperature for a H2O:CO2=1:0.25 ice (bottom).

##  6 Infrared band assignments

Table 2 provides an overview of the positions of the fitted Gaussians of the
asymmetric stretching bands, denoted as $\tilde{\nu}_{3,1}$,
$\tilde{\nu}_{3,2}$, and $\tilde{\nu}_{3,3}$. This was done for a pure CO2
ice, low or high CO2 concentration in water ice. These measurements span from
10 K (the deposition temperature) to 170 K, which corresponds to the
desorption temperature of water in these experiments.
For pure CO2 ice, the positions of $\tilde{\nu}_{3,1}$ and
$\tilde{\nu}_{3,2}$, respectively 2351 cm-1 (4.25 $\mu$m) and 2345 cm-1 (4.26
$\mu$m) at 10 K, are affected by temperature, shifting towards the red as the
temperature increases (as also demonstrated in Figure 3). A similar behaviour
is observed both when low and high CO2 concentrations are mixed with water. We
note that the $\tilde{\nu}_{3,2}$ (at 10 K it falls at 2345 cm-1 or 4.26
$\mu$m) shift for high CO2 concentration is much smaller than for the low
concentration. At 10 K, a result of $\tilde{\nu}_{3,1}$ (2351 cm-1, 4.25
$\mu$m) for a low CO2 concentration closely aligns with the position in pure
CO2 ice, whereas a high concentration results in a blueshift of over 2 cm-1.
In addition, $\tilde{\nu}_{3,2}$ at high concentration (2344 cm-1, 4.27
$\mu$m) is redshifted by 0.7 cm-1 compared to the low concentration and pure
CO2 ice. The $\tilde{\nu}_{3,3}$ appears only when mixed with water. The peak
position of $\tilde{\nu}_{3,3}$ (2337 cm-1, 4.28 $\mu$m) in high CO2
concentration experiences a redshift of nearly 4 cm-1, compared to low
concentration. These shifts are due to the band broadening (Ehrenfreund et al.
(1997)), and cause a shift of the Gaussian fits that is observed when the CO2
concentration is increased in water ice.

When mixed with water, the thermal evolution of the CO2 $\tilde{\nu}_{3,1}$,
$\tilde{\nu}_{3,2}$, and $\tilde{\nu}_{3,3}$ are reported in Figure 10. This
figure depicts the three Gaussians from 10 K to 170 K. The top panel shows a
low CO2 concentration whereas the bottom panel corresponds to a high
concentration. The Gaussians for 10 K are highlighted for clarity representing
the beginning of the experiment. The arrows indicate the direction of
evolution with increasing temperature, clearly illustrating how
$\tilde{\nu}_{3,3}$ (2337 cm-1, 4.28 $\mu$m) shifts to the blue only for a
high CO2 concentration. For a low concentration, the peak intensity of
$\tilde{\nu}_{3,3}$ (2341 cm-1 at 10 K, or 4.27 $\mu$m) first increases after
80 K, which becomes CO2 embedded in the water ice matrix, as illustrated in
the top panel of Figure 10. At 80 K, near the desorption temperature of pure
CO2, $\tilde{\nu}_{3,1}$ shifts to lower wavenumbers for both low and high CO2
concentrations, their respective positions at this temperature are about 2348
cm-1 (4.26 $\mu$m) and 2352 cm-1 (4.25 $\mu$m). At 140 K, the integrated
absorbance areas decrease due to the transition from cubic to hexagonal ice
structure (Martín-Doménech et al. (2014)). It is worth noting that in Table 2,
for temperatures above 160 K, $\tilde{\nu}_{3,1}$ shifts back towards the
blue, ultimately resulting in a single peak at around 2348 cm-1 (4.259 $\mu$m)
at 170 K for low and 2349.1 cm-1 (4.257 $\mu$m) for high CO2 concentration.

Figure 10: Three Gaussian distributions plotted from 10 K to 170 K, for a low
and high CO2 concentration in the ice. Blue represents $\tilde{\nu}_{3,1}$,
orange $\tilde{\nu}_{3,2}$, and green $\tilde{\nu}_{3,3}$. The bold graphs
indicate the beginning of the experiment at 10 K. The arrows indicate the
direction of evolution.

Furthermore, we explored the influence of deposition pressure on the band
positions. To investigate this aspect, an experiment was conducted with a
deposition pressure of 8$\times$10-7 mbar and a 4$\%$ CO2 concentration. Upon
analyzing the band positions at 10 K, it was observed that the positions of
$\tilde{\nu}_{3,2}$ (2345 cm-1, 4.26 $\mu$m) and $\tilde{\nu}_{3,3}$ (2341
cm-1, 4.27 $\mu$m) were consistent with the experiment conducted at a lower
deposition pressure. The higher deposition pressure is known to induce changes
in the ice morphology, rendering it more porous (Maté et al., 2008; Bossa et
al., 2015). The DFT simulations detailed in section 7 show the effect porosity
has on amorphous CO2 ice, allowing the CO2 dangling bond to occur. As the ice
produced at higher pressure is more porous, it undergoes greater
reorganization during the warm-up process. At 10 K, $\tilde{\nu}_{3,1}$ has a
small $\sim$1 cm-1 blue shift compared to the low deposition pressure. When
the temperature reached 50 K, CO2 became crystalline and the reorganization
was complete. At this point, $\tilde{\nu}_{3,1}$, $\tilde{\nu}_{3,2}$, and
$\tilde{\nu}_{3,3}$ were all peaking at approximately 2349.8 cm-1 (4.256
$\mu$m), 2344.5 cm-1 (4.265 $\mu$m), and 2340.6 cm-1 (4.272 $\mu$m),
respectively. These band positions are similar to the low deposition pressure
experiment. The assignment of the different band components is discussed in
the remainder of this section. As previously mentioned, when CO2 is mixed with
water, three distinct peaks emerge.

The $\tilde{\nu}_{3,2}$ band found at approximately 2345 cm-1 (4.26 $\mu$m)
can be attributed to bulk CO2 ice. This is corroborated by the desorption
temperature of pure CO2, which falls in the range of 80-85 K for both low and
high CO2 concentrations. This desorption is illustrated in Figure 4, with band
positions shown in Table 2. This band is also referred to as the CO2-ext where
the interaction with water is quite weak as expected if CO2 is superficially
adsorbed, as proposed by Gálvez et al. (2007).

The $\tilde{\nu}_{3,3}$ band arises when CO2 is mixed with water and peaks
around 2340 cm-1 (4.27 $\mu$m). This band is the consequence of CO2 molecules
interacting with the water ice matrix. This peak is also referred to as
CO2-int in previous work (Gálvez et al., 2007; Maté et al., 2008), where
individual CO2 molecules are trapped in the amorphous water ice and this
interaction results in a slight weakening of the C-O bond, which produces a
redshift on the IR spectrum (Sandford & Allamandola, 1990). When increasing
the temperature the integrated absorbance area of $\tilde{\nu}_{3,3}$ near
2340 cm-1 (4.27 $\mu$m) diminishes as the water ice evolves to a hexagonal
structure and the CO2 molecules desorb from the water ice matrix.

Here, $\tilde{\nu}_{3,1}$ was found at approximately 2351 cm-1 (4.254 $\mu$m)
for low CO2 concentration and around 2353 cm-1 (4.250 $\mu$m) for high
concentration. This band is present in pure CO2 ice, meaning that it is not
caused by the interaction with water but is instead a consequence of the
morphology of the ice. The band has been previously observed in ice growth
experiments with high deposition rates (Gálvez et al., 2008), although it was
not observable at lower deposition rates (Falk, 1987). We propose that this
band is caused by pores in the ice where the binding is particularly weak, as
it occurs on the top surface of the ice, and the CO2 molecules are dangling
from the pore surface. This preliminary assignment is supported by a
vibrational mode that is connected to porosity, namely, a degree of freedom to
vibrate that is not allowed in a compact CO2 ice structure; furthermore, the
2351 cm-1 position relative to that of the fundamental C=O stretch in CO2 ice
is analogous to the relative position of the O-H dangling with respect to the
O-H stretch in water ice (Matsuda et al., 2018).This phenomenon is further
supported by the discussion in section 7. Therefore, we assigned
$\tilde{\nu}_{3,1}$ to the CO2 dangling bonds. When CO2 is mixed with water,
this peak is also present; however, in this case, it remains above the 80 K
temperature, which implies that following the desorption of bulk CO2 ice,
carbon dioxide molecules can remain trapped in the pores of the amorphous
solid water.

A similar explanation for this band location, involving molecules bonded to
-OH dangling bonds in the water ice, was suggested by Matsuda et al. (2018)
for CO molecules. In that case, two bands are assigned, namely CO molecules
interacting with the OH dangling groups (CO-dangling OH) and CO molecules
interacting with the oxygen atoms of the surface water molecules (CO-bonded
OH). It was found that the area of the CO-dangling band increases when cycling
to higher temperatures, meaning that CO molecules partially diffuse to CO-
dangling OH groups on the surface of the amorphous water ice pores.

The position of this band, blueshifted from the pure CO2 stretching band,
could be easily confused with a similar band observed in clathrate hydrates
(Dartois & Schmitt, 2009). In our case we can discard these cage-like
structures as they are normally formed at much slower deposition rates and
higher temperatures.

Other proposed explanations for this band location is that carbon dioxide is
in a “complexed” form with water, as suggested by Chaban et al. (2007). A blue
shift of 5 cm-1 compared to the pure CO2 ice is observed when bonded with one
water molecule, and 10 cm-1 with two. Furthermore, the symmetric O-H
stretching band undergoes a redshift when complexed with two water molecules
and the intensity increases significantly. In this work, for a 4$\%$ CO2
concentration, there was no increase in intensity nor a redshift at 10 K for
the symmetric O-H stretching band when compared to that of a pure H2O ice,
ruling out the possibility of a complexed CO2.

##  7 Assignment of $\nu_{3.1}$ with DFT calculations

In order to better understand the influence of porosity in amorphous CO2,
simulations were carried out using density functional theory (DFT) (Hohenberg
& Kohn, 1964; Kohn & Sham, 1965) using the CASTEP program (Clark et al., 2005)
for geometry optimization and vibrational spectrum prediction.

Starting with a CO2 crystal, the amorphous ice models are created by applying
temperature using molecular dynamics (MD) in the NPT ensemble until the
structure melts. This was done with the Andersen barostat (Andersen, 1980) at
constant pressure and the Nose thermostat (Martyna et al., 1992) at rising
temperatures. The porous model is then formed by removing CO2 molecules from
the amorphous model’s core. DFT was then used to geometrically optimize the
systems using the generalized gradient approximation (GGA) and functionals by
Perdew-Burke-Ernzerhof (PBE). Infrared spectra were simulated using density
functional perturbation theory (Refson et al., 2006) based on these improved
models.

The calculated spectra are not intended to perfectly duplicate the measured
wavenumbers for vibrational bands, but should be near within a margin of
error. Most significantly, they are repeatable and can aid in understanding
the influence of experimentally uncontrollable factors.

Figure 11: DFT calculated spectra for bulk amorphous CO2 ice (blue) and a
simulated pore in the same system (orange) as a function of wavenumber.

The simulated IR spectra of an amorphous CO2 model (blue line) and a porous
CO2 model (orange line) are shown in Figure 11. The morphology of the bulk ice
model (right) and the porous model (left) are included for comparison. In the
case of the bulk amorphous ice, we can see that the spectrum is very similar
to the experimental case, with some expected offset for the wavenumber. In the
porous model there is a clear blueshift of appproximately 4 cm-1 for the main
peak. Most importantly, a new secondary peak emerges at higher wavenumbers.
This pattern is very similar to the spectra seen in laboratory ice, as well as
the recent observations of this same CO2 band in Ganymede. It must be stressed
that the goal of this simulation is to simulate the effect of porosity on the
peak location and shape. Therefore, finding the mechanism behind the shift
rather than matching perfectly the band location with the same wavenumbers.

##  8 Applications to icy moons

In this section, we use JWST ERS data from Ganymede (Bockelée-Morvan et al.,
2024) and Europa (Villanueva et al., 2023) to study the position and evolution
of the CO2 bands with varying latitude and longitude. The data were
reprocessed from the Mikulski Archive for Space Telescopes (MAST) as described
in Bockelée-Morvan et al. (2024). Our goal is to assign the CO2 bands observed
to the physical state of CO2 ice. On Ganymede, the amount of CO2 relative to
water in mass is of the order of 1$\%$ (Bockelée-Morvan et al., 2024), which
implies that our low concentration experiments are more representative of
these icy moons’ surfaces.

###  8.1 Ganymede

This section offers a discussion of the latitude and longitudinal spectra of
Ganymede’s leading hemisphere, as similar trends are observed on the trailing
hemisphere.

Figure 12 shows the CO2 band observed with JWST NIRspec on Ganymede as a
function of latitude. The solid lines are JWST spectra and the colours
indicate different latitudinal ranges, in good agreement with Bockelée-Morvan
et al. (2024). The latitudinal variations of the CO2 band and its asymmetric
band shape support the view that different CO2 ice physical states are
involved. The 2342 cm-1 (4.27 $\mu$m) peak is stronger across the northern
latitudes, which are colder and more enriched in water ice (Bockelée-Morvan et
al., 2024; Trumbo et al., 2023). Across the equatorial latitudes, the 2353
cm-1 (4.25 $\mu$m) band becomes stronger. Ganymede has an average surface
temperature around 100 K, with a maximum of 160 K at the equator (Bockelée-
Morvan et al., 2024). Each latitude spectrum can be deconvolved into two
Gaussian distribution that are represented with dotted lines and plotted for
each range of latitudes. When analyzing the results, we can see that one of
the Gaussians is located at 2341.7 cm-1 (4.270 $\mu$m) at the poles and shifts
to lower wavelengths for other regions (supposedly warmer). This is in
agreement with the $\tilde{\nu}_{3,3}$ band observed experimentally,
attributed to CO2 interacting with H2O molecules in the ice, which shifts
slightly to the blue as the temperature increases. Moving closer to the
equator, this peak shifts even more to the blue, reaching its maximum at
2342.5 cm-1 (4.269 $\mu$m). The second Gaussian is located around 2350.3 cm-1
(4.255 $\mu$m) at the North pole and around 2350.9 cm-1 (4.254 $\mu$m) at
other latitudes. This is also in agreement with the $\tilde{\nu}_{3,1}$ band
observed experimentally, which is blue-shifting from temperature above 150 K,
which coincides with the surface temperature at Ganymede’s equator. In line
with the experimental results, the presence of the $\tilde{\nu}_{3,1}$ band is
identified as CO2 dangling bonds, which are an indication of ice porosity.

Figure 12: Latitude spectra of Ganymede’s leading hemisphere acquired with
JWST. Colours are used to indicate different latitude ranges. The dotted lines
represent the Gaussian fit of each individual spectrum which is deconvolved in
two Gaussians.

The longitudinal spectra for Ganymede’s leading hemisphere are shown in Figure
13. The solid lines are JWST spectra and the colours indicate different
longitudinal ranges. The variations among the CO2 band position and intensity
are shown in the figure for latitudes between 60∘ north and 60∘ south. The
morning limb is represented by the 120-150 W spectrum (blue), whereas the
evening limb is represented by the 30-60 W spectrum (red). Each spectrum is
deconvolved into two Gaussian distributions, which are shown by dotted lines.
The band at 2351 cm-1 (4.25 $\mu$m) has been assigned to the
$\tilde{\nu}_{3,1}$ band being an indication of the CO2 dangling bonds and of
porosity. When comparing morning (blue) to evening (red), the
$\tilde{\nu}_{3,1}$ band increases, suggesting that porosity increases during
a day on Ganymede. This increase in porosity could be attributed to cracks in
the ice that are forming due to temperature fluctuations. The effect of UV
irradiation could also play an important role in this change in porosity that
should also be considered.

Figure 13: Longitudinal spectra of Ganymede’s leading hemisphere acquired with
JWST. The longitudes are taken between 60∘ North and 60∘ South in latitude.
The spectrum for 120-150 W represents the morning limb of the leading edge and
the evening limb by the 30-60 W spectrum. Colors are used to indicate
different longitudinal ranges. The dotted lines represent the Gaussian fit of
each individual spectrum which is deconvolved into two Gaussians.

###  8.2 Europa

Figure 14 shows the CO2 band observed on Europa as a function of the latitude
and exhibits similar variation as that reported by Villanueva et al. (2023).
The band ca. 2342 cm-1 (4.27 $\mu$m) is stronger across the northern
latitudes, which are colder and more enriched in water ice (Trumbo & Brown
2023). Across the equatorial latitudes, the 2353 cm-1 (4.25 $\mu$m) peak
becomes stronger. The temperature on Europa ranges from 60 K at the poles to
around 110 K at the equator, which is colder than on Ganymede. The Gaussian
fits of the data are represented with dotted lines, and plotted for each range
of latitudes. Each latitude spectrum can be deconvolved into two Gaussian
distributions (as seen in Figure 14). The band around 2342 cm-1 (4.27 $\mu$m)
feature shifts when moving closer to the equator, from 2340.9 cm-1 (4.272
$\mu$m) for 30∘ to 60∘ North to 2341.2 cm-1 (4.271 $\mu$m) at the equator.
This is in agreement with the $\tilde{\nu}_{3,3}$ band position observed
experimentally, which is located around 2340.8 cm-1 (4.272 $\mu$m) and shifts
to the blue with increasing temperatures. The second Gaussian
$\tilde{\nu}_{3,1}$, located at 2353.5 cm-1 (4.249 $\mu$m), is shifted to the
blue by around 3 cm-1 when compared to Ganymede. This is consistent with our
experimental results showing the $\tilde{\nu}_{3,1}$ band shifts to higher
wavelength with decreasing temperature.

Figure 14: Latitude spectra of Europa’s leading hemisphere acquired with JWST.
Colours are used to indicate different latitude ranges. The dotted lines
represent the Gaussian fit of each individual spectrum, each consisting of two
Gaussians.

The main conclusion of this section is that CO2 observations of both Ganymede
and Europa can be explained by the presence of two components, the
$\tilde{\nu}_{3,1}$ band and $\tilde{\nu}_{3,3}$ band. The $\tilde{\nu}_{3,3}$
band located at 2342 cm-1 (4.27 $\mu$m) represents CO2 trapped in water ice
and can account for the CO2 band observed in Ganymede’s and Europa’s spectra
of the leading hemisphere, and is stronger in the boreal regions especially
for Ganymede. The $\tilde{\nu}_{3,1}$ band at 2353 cm-1 (4.25 $\mu$m) for
Europa and at 2350.9 cm-1 (4.254 $\mu$m) for Ganymede is stronger in other
regions for both Ganymede and Europa. This band is associated with CO2
dangling bonds which depend on the pores present in the ice, as CO2 molecules
on the pore walls can exhibit such dangling bonds. The band assignment of
$\tilde{\nu}_{3,1}$ to CO2 in pores has been suggested, based on our
experimental data combined with DFT simulations. It should be noted that the
2353 cm-1 (4.25 $\mu$m) band depends on the concentration of CO2 in the water
ice, which will make the band shift.

##  9 Conclusions

In this study, we combined FTIR measurements and TPD experiments to explain
the complexity of the CO2 stretching band. The examination of CO2 ice
behaviour on icy moons has offered intriguing revelations about its
interaction within these distinct environments. Moreover, aligning laboratory
findings and DFT calculations with JWST data from Ganymede and Europa, bridges
the gap between theory, experiments, and celestial icy surfaces. These
findings can help enrich our understanding of these moons’ environments.

  * •

The C–O stretching band observed for pure CO2 ice at 10 K consists of two
distinct vibrations: $\tilde{\nu}_{3,1}$ which is caused by porosity and peaks
at 2351.3 cm-1 (4.253 $\mu$m), and $\tilde{\nu}_{3,2}$ which corresponds to
bulk CO2 ice and peaks at 2345 cm-1 (4.264 $\mu$m). The band’s shape evolves
during warm-up from 10 K to 80 K, exhibiting a redshift in both Gaussians. The
TPD data indicates a desorption peak at 85 K, suggesting a thermal desorption
of CO2 from the ice layer.

  * •

CO2 mixed with water ice (at low and high concentration) exhibits multiple
peaks, of which two are the same as for pure CO2. The presence of three peaks,
$\tilde{\nu}_{3,1}$ at 2351 cm-1 (4.25 $\mu$m), $\tilde{\nu}_{3,2}$ at 2345
cm-1 (4.26 $\mu$m), and $\tilde{\nu}_{3,3}$ at 2341 cm-1 (4.27 $\mu$m) at 10 K
demonstrates the complexity of CO2 interaction within water ice. The different
CO2 bands are assigned as follows: $\tilde{\nu}_{3,1}$ (2351 cm-1, 4.25
$\mu$m) is associated to CO2 dangling bonds in which CO2 molecules are found
in pores or cracks, $\tilde{\nu}_{3,2}$ (2345 cm-1, 4.26 $\mu$m) is due to CO2
trapped in the water ice (segregated) forming pockets of pure CO2, and
$\tilde{\nu}_{3,3}$ (2341 cm-1, 4.27 $\mu$m) is related to CO2 molecules
embedded in the ice interacting with water molecules.

  * •

Thermal desorption analysis allowed for distinct peaks in the desorption
curves to be assigned, indicating different desorption behaviours for CO2
within the ice. The low and high concentrations display similar desorption
patterns. Initial CO2 desorption occurs around 80 K, followed by peaks at 146
K and 162 K, indicating various CO2-water interactions during ice phase
transitions. Continuous desorption between 80 and 146 K is observed, driven by
the amorphous water ice rearrangement and thereby facilitating CO2 diffusion
through pores in the ice matrix.

  * •

The JWST NIRSPEC spectra of Ganymede’s leading hemisphere reveals significant
variations in the CO2 band profile across different latitudes. Two main CO2
bands have been observed, that we assigned to $\tilde{\nu}_{3,1}$ (2353 cm-1,
4.25 $\mu$m) and $\tilde{\nu}_{3,3}$ (2342 cm-1, 4.27 $\mu$m). The dominance
of $\tilde{\nu}_{3,1}$ (2353 cm-1, 4.25 $\mu$m) in northern latitudes,
associated with colder regions enriched in water ice, is in contrast with the
prevalence of $\tilde{\nu}_{3,3}$ (2342 cm-1, 4.27 $\mu$m) in equatorial
latitudes. Gaussian fits of the spectra suggest two distinct physical states
of CO2 that shift with temperature, which confirm the assignments from our
laboratory findings. The CO2 band shifts observed on Ganymede with latitude
could be attributed to an increase in temperature, showing that in the poles
CO2 is well embedded in water ice.

  * •

The longitudinal spectra for Ganymede’s leading hemisphere indicates
variations in CO2 spectra from morning to evening. The band is assigned to the
$\tilde{\nu}_{3,1}$ band (2351 cm-1, 4.25 $\mu$m) that links to the porosity
of the ice. This increase in porosity could be attributed to cracks created in
the ice due to variation of temperatures.

Similarly, the JWST NIRSPEC spectra of Europa’s leading hemisphere can be
deconvolved into two Gaussian fits: $\tilde{\nu}_{3,3}$ (2342 cm-1, 4.27
$\mu$m) assigned to CO2 molecules trapped in the water ice and
$\tilde{\nu}_{3,1}$ (2353 cm-1, 4.25 $\mu$m) assigned to dangling CO2
molecules in pores or cracks. Temperatures in Europa are colder than Ganymede
and as expected show a blueshift for $\tilde{\nu}_{3,1}$ and a redshift for
$\tilde{\nu}_{3,3}$ when compared to the warmer moon. A similar effect is
observed when comparing spectra from the warmer equator to the colder poles,
indicating higher porosity near the equator.

###### Acknowledgements.

The authors would like to thank I. de Pater, T. Fouchet, and D. Bockelée-
Morvan for the use of the JWST ERS Ganymede data and assistance with
processing, as well as G. Villanueva for the use of the JWST Europa data. This
research has been funded by project PID2020-118974GB-C21 by the Spanish
Ministry of Science and Innovation. B.E. acknowledges support by grant
PTA2020-018247-I by the Spanish Ministry of Science and Innovation/State
Agency of Research MCIN/AEI.

## References

  * Andersen (1980) Andersen, H. C. 1980, J. Chem. Phys., 71, 2384
  * Bernstein et al. (2005) Bernstein, M. P., Cruikshank, D. P., & Sandford, S. A. 2005, Icarus, 179, 527
  * Bockelée-Morvan et al. (2024) Bockelée-Morvan, D., Lellouch, E., Poch, O., et al. 2024, Astronomy & Astrophysics, 681, A27
  * Bossa et al. (2015) Bossa, J.-B., Maté, B., Fransen, C., et al. 2015, The Astrophysical Journal, 814, 47
  * Bouilloud et al. (2015) Bouilloud, M., Fray, N., Bénilan, Y., et al. 2015, Monthly Notices of the Royal Astronomical Society, 451, 2145
  * Buratti et al. (2005) Buratti, B., Cruikshank, D. P., Brown, R. H., et al. 2005, The Astrophysical Journal, 622, L149
  * Cazaux et al. (2015) Cazaux, S., Bossa, J.-B., Linnartz, H., & Tielens, A. 2015, Astronomy & Astrophysics, 573, A16
  * Chaban et al. (2007) Chaban, G. M., Bernstein, M., & Cruikshank, D. P. 2007, Icarus, 187, 592
  * Clark et al. (2005) Clark, S. J., Segall, M. D., Pickard, C. J., et al. 2005, Z. Kristall., 220, 567
  * Crovisier (2006a) Crovisier, J. 2006a, Faraday discussions, 133, 375
  * Crovisier (2006b) Crovisier, J. 2006b, Molecular Physics, 104, 2737
  * Dartois et al. (2005) Dartois, E., Pontoppidan, K., Thi, W.-F., & Caro, G. M. 2005, Astronomy & Astrophysics, 444, L57
  * Dartois & Schmitt (2009) Dartois, E. & Schmitt, B. 2009, Astronomy & Astrophysics, 504, 869
  * Draine (2003) Draine, B. T. 2003, Annual Review of Astronomy and Astrophysics, 41, 241
  * Ehrenfreund et al. (1997) Ehrenfreund, P., Boogert, A., Gerakines, P., Tielens, A., & Van Dishoeck, E. 1997, Astronomy and Astrophysics, 328, 649
  * Falk (1987) Falk, M. 1987, The Journal of Chemical Physics, 86, 560
  * Gálvez et al. (2008) Gálvez, O., Maté, B., Herrero, V. J., & Escribano, R. 2008, Icarus, 197, 599
  * Gálvez et al. (2007) Gálvez, O., Ortega, I. K., Maté, B., et al. 2007, Astronomy & Astrophysics, 472, 691
  * Gerakines et al. (1995) Gerakines, P. A., Schutte, W. A., Greenberg, J. M., & van Dishoeck, E. F. 1995, A&A, 296, 810
  * González Díaz et al. (2022) González Díaz, C., Carrascosa, H., Muñoz Caro, G. M., Satorre, M., & Chen, Y.-J. 2022, Monthly Notices of the Royal Astronomical Society, 517, 5744
  * Grundy et al. (2003) Grundy, W., Young, L., & Young, E. 2003, Icarus, 162, 222
  * Hagen et al. (1981) Hagen, W., Tielens, A., & Greenberg, J. 1981, Chemical Physics, 56, 367
  * Hansen & McCord (2008) Hansen, G. B. & McCord, T. B. 2008, Geophysical Research Letters, 35
  * Hohenberg & Kohn (1964) Hohenberg, P. & Kohn, W. 1964, Phys. Rev., 136, B864
  * Isokoski et al. (2013) Isokoski, K., Poteet, C., & Linnartz, H. 2013, Astronomy & Astrophysics, 555, A85
  * Kohn & Sham (1965) Kohn, W. & Sham, L. J. 1965, Phys. Rev., 140, A1133
  * Kumi et al. (2006) Kumi, G., Malyk, S., Hawkins, S., Reisler, H., & Wittig, C. 2006, The Journal of Physical Chemistry A, 110, 2097
  * Malyk et al. (2007) Malyk, S., Kumi, G., Reisler, H., & Wittig, C. 2007, The Journal of Physical Chemistry A, 111, 13365
  * Martín-Doménech et al. (2014) Martín-Doménech, R., Caro, G. M., Bueno, J., & Goesmann, F. 2014, Astronomy & Astrophysics, 564, A8
  * Martyna et al. (1992) Martyna, G. L., Klein, M. L., & Tuckerman, M. 1992, J. Chem. Phys., 97, 2635
  * Matsuda et al. (2018) Matsuda, S., Yamazaki, M., Harata, A., & Yabushita, A. 2018, Chemistry Letters, 47, 468
  * Maté et al. (2008) Maté, B., Gálvez, O., Martín-Llorente, B., et al. 2008, The Journal of Physical Chemistry A, 112, 457
  * McCord et al. (1998) McCord, T. B., Hansen, G. B., Clark, R. N., et al. 1998, J. Geophys. Res., 103, 8603
  * Mishra et al. (2021) Mishra, I., Lewis, N., Lunine, J., et al. 2021, The Planetary Science Journal, 2, 183
  * Muñoz Caro et al. (2010) Muñoz Caro, G., Jiménez-Escobar, A., Martín-Gago, J., et al. 2010, Astronomy & Astrophysics, 522, A108
  * Refson et al. (2006) Refson, K., Tulip, P. R., & Clark, S. J. 2006, Phys. Rev. B, 73, 155114
  * Sandford & Allamandola (1990) Sandford, S. & Allamandola, L. 1990, Astrophysical Journal, 355, 357
  * Trumbo & Brown (2023) Trumbo, S. K. & Brown, M. E. 2023, Science, 381, 1308
  * Trumbo et al. (2023) Trumbo, S. K., Brown, M. E., Bockelée-Morvan, D., et al. 2023, Science Advances, 9, eadg3724
  * Villanueva et al. (2023) Villanueva, G., Hammel, H., Milam, S., et al. 2023, Science, 381, 1305


:
LEAVE UNSET 2021 LEAVE UNSET LEAVE UNSET Machine Learning for Health (ML4H)
2021

#  Optimizing Binary Symptom Checkers via Approximate Message Passing

Mohamed<EMAIL_ADDRESS>
AIP Labs  Hungary  Faouzi Bellili, Amine Mezghani{Faouzi.Bellili,
<EMAIL_ADDRESS>
University of Manitoba  Canada  Hayet Amdouni<EMAIL_ADDRESS>
Self-employed Dermatologist  Tunisia

###### Abstract

Symptom checkers have been widely adopted as an intelligent e-healthcare
application during the ongoing pandemic crisis. Their performance have been
limited by the fine-grained quality of the collected medical knowledge between
symptom and diseases. While the binarization of the relationships between
symptoms and diseases simplifies the data collection process, it also leads to
non-convex optimization problems during the inference step. In this paper, we
formulate the symptom checking problem as an underdertermined non-convex
optimization problem, thereby justifying the use of the compressive sensing
framework to solve it. We show that the generalized vector approximate message
passing (G-VAMP) algorithm provides the best performance for binary symptom
checkers.

###### keywords:

Symptom checkers, approximate message passing.

##  1 Introduction

With the growing shortage of doctors around the world, the geographical reach
of the online symptom assessment tools, particularly symptom checkers, have
provided alternative channels to better support patients at the very beginning
of their diagnostic journey. The challenges faced by symptom checkers can be
categorized into three classes:

Medical knowledge data collection: It is usually a tedious task to define the
medical knowledge represented by the matrix $\mathbf{A}$ in Fig. 1 since it
requires a manual labor-intensive process to define the probabilistic
relationships $p(s_{i}|d_{j})$ between each symptom $s_{i}\in\mathcal{S}$ and
disease $d_{j}\in\mathcal{D}$. While few studies (Rotmensch et al. (2017))
explored the possibility to automate this process by learning high-quality
knowledge bases associating diseases to symptoms directly from electronic
medical records, the data-driven medical knowledge suffer from underestimating
rare disease medical knowledge and requires an expert assessment.

Sparse underdertermined inference: The performance of most symptom checkers is
limited by the number of observed symptoms which is usually a binary input,
i.e., either the patient has the symptom or not. Additionally, the accuracy of
the patient input is usually biased due to the medical vocabulary gap,
incomplete information, and the inability to differentiate between highly
correlated medical concepts. Moreover, the common design of symptom checkers
usually require the number of conditions to be higher than the number of
symptoms, i.e., $M<N$. This makes the problem in (1) an underdertermined
inference problem, which is precisely the typical assumption for recovering
the sparse disease vector $\bm{d}$.

Unbalanced datasets: Symptom checkers based on machine learning algorithms
such as supervised learning (Choi et al. (2016)) and reinforcement learning
(Akrout et al. (2019)) are biased toward unbalanced datasets whose real-world
relative frequencies of common vs. rare diseases are highly disproportional
(Fraser et al. (2017)).

While early diagnosis is the most effective way to detect early diseases and
reduce mortality (Coleman (2017)), symptom checkers111See WebMD and
Symptomate’s symptom checkers as examples. still face three main challenges to
accurately infer the correct disease $d_{i}^{*}$ among a set of $N$ supported
diseases $\mathcal{D}=\\{d_{1},\,\dots,\,d_{N}\\}$ given $K$ observed symptoms
subset $\mathcal{S}_{\textrm{obs}}$ among a set of $M$ supported symptoms
$\\{s_{1},\,\dots,\,s_{M}\\}$, i.e.,
$\mathcal{S}_{\textrm{obs}}\subset\mathcal{S}$. That is to say:

$\small\widehat{\mathbf{d}}=\operatorname*{argmin}_{\mathbf{d}\in\mathbb{R}^{N}}\frac{1}{2}\|\mathbf{s}-\phi(\mathbf{A}\,\mathbf{d})\|_{2}^{2}+f(\mathbf{d}).$
(1)

where $\phi(\cdot)$ is a non-linear function applied component-wise while
$f(\mathbf{d})$ is a regularization function chosen to promote a desired
structure in the unknown disease vector $\mathbf{d}$. For example, letting
$f(\mathbf{d})=\lambda\,\|\mathbf{d}\|_{1}$ with $\lambda>0$ promotes sparsity
in $\widehat{\mathbf{d}}$.

Figure 1: The symptom checking model: a patient is represented by a disease
vector $\mathbf{d}$ sampled from a predefined disease density $p(\mathbf{d})$,
whose product with the matrix $\bm{A}$ containing the probabilistic
relationships between symptoms and diseases leads to the dense symptom vector
$\mathbf{z}$. The observed symptom vector $\mathbf{s}$ is a sparse version of
$\mathbf{z}$ i.e. $\mathbf{s}=\phi(\mathbf{z+\mathbf{w}})$ where
$\mathbf{w}\sim\mathcal{N}(\mathbf{w};\bm{0},\gamma_{w}^{-1}\,\mathbf{I})$.

In this work, we propose a new methodology to jointly tackle the two
challenges of symptom checkers by advocating the use of message passing
algorithms for compressive sensing (Donoho et al. (2009)) to draw the
following potential benefits:

Binary medical knowledge: While collecting the probabilistic medical knowledge
for each symptom-disease pair of the matrix $\mathbf{A}$ is a tedious task,
collecting their binary version is more practical in terms of data collection
speed and expert consensus convergence. In fact, most of the expert consensus
issues lie in the magnitude of the relationship between any symptom-disease
pair, which is automatically removed via the binarization.

Efficient algorithms for sparse inference: The fact that the vector
$\mathbf{d}$ is sparse in practice urges the need to use low-complexity
compressive sensing algorithms such as GAMP (Rangan (2011)) and VAMP (Rangan
et al. (2019)) which are well-suited for the under-complete problems in (1)
under non-convex sparsity-inducing penalty functions. The performance of
symptom checkers based on sparse Bayesian inference are model-based algorithms
which only rely on a predefined prior on the unknown disease vector
$\mathbf{d}$.

##  2 Dermatology case study

To evaluate the performance of the proposed symptom checker based on message
passing, we describe in this section our the data collection process for our
dermatological study case. While many possible medical specialities were
possible to explore, we decided to focus on dermatology due to the
accessibility of the visual skin information, which will allow us to integrate
the proposed symptom checker with visual classifiers in future works.

###  2.1 Medical knowledge collection

We collaborated with three practicing dermatologists to elicit the binary
medical knowledge matrix $\mathbf{A}$ required for the symptom checker model
in (1). The participating doctors were encouraged to use the available medical
literature and their experience in the process of creating the medical
knowledge matrix. We explicitly asked them to provide the binary conditional
probabilities $p(s_{i}|d_{j})$, between each symptom $s_{i}\in\mathcal{S}$ and
disease $d_{j}\in\mathcal{D}$ covering the following $M=27$ symptoms and
$N=31$ diseases:

Set of diseases $\mathcal{D}=\\{$psoriasis, seborrheic dermatitis, lichen,
acne, atopic dermatitis, keloids, cheilitis, condylomata, candidiasis,
dermatophyties, stasis dermatitis, dysidrosis, erysipelas, bedsores,
folliculitis, hidradenitis suppurativa, cutaneous leishmaniasis, lupus
erythematosus, melanoma, noevus, rosacea, toxidermia, ulcer venous, urticaria,
varicella, herpes, zoster, sarcoidosis$\\}$,

Set of symptoms $\mathcal{S}=\\{$redness, dander, vesicles, bubbles,
pigmentation, swelling, pustule, macule, plate, nodule, papule, crusts,
hypochromia, atrophy, fever, pain, pruritus, oozing, hyperkeratosis, cracks,
ulceration, ulcer, edema, induration, necrosis, infiltration,
telangiectasias$\\}$.

$p_{\mathsf{\mathbf{d}}}(\mathbf{d})$
$\mathbf{\mathsf{\mathbf{{z}}}=\mathbf{A}\,\mathsf{\mathbf{d}}}$
$\mathsf{\mathbf{{{s}}}}=\phi(\mathsf{\mathbf{{{z}}}}+\mathsf{\mathbf{{{w}}}})$
extextextext$\mathbf{\widehat{d}}^{+}_{\mathsf{p}}$$\gamma_{\bm{\mathsf{d}}_{\mathsf{p}}}^{+}$$\mathbf{\widehat{d}}^{+}_{\mathsf{e}}$$\gamma_{\bm{\mathsf{d}}_{\mathsf{e}}}^{+}$$\mathbf{\widehat{z}}^{+}_{\mathsf{p}}$$\gamma_{\bm{\mathsf{z}}_{\mathsf{p}}}^{+}$$\mathbf{\widehat{z}}^{+}_{\mathsf{e}}$$\gamma_{\bm{\mathsf{z}}_{\mathsf{e}}}^{+}$$\mathbf{\widehat{z}}^{-}_{\mathsf{p}}$$\gamma_{\bm{\mathsf{z}}_{\mathsf{p}}}^{-}$$\mathbf{\widehat{z}}^{-}_{\mathsf{e}}$$\gamma_{\bm{\mathsf{z}}_{\mathsf{e}}}^{-}$$\mathbf{\widehat{d}}^{-}_{\mathsf{p}}$$\gamma_{\bm{\mathsf{d}}_{\mathsf{p}}}^{-}$$\mathbf{\widehat{d}}^{-}_{\mathsf{e}}$$\gamma_{\bm{\mathsf{d}}_{\mathsf{e}}}^{-}$
Figure 2: Block diagram of the generalized VAMP algorithm with its three
modules: the prior modules $p_{\mathsf{\mathbf{{{d}}}}}(\bm{d})$, the linear
minimum mean square estimator (LMMSE) module
$\mathbf{\mathsf{\mathbf{{{z}}}}=\mathbf{A}\,\mathsf{\textbf{{{d}}}}}$ and the
non linear block
$\mathsf{\mathbf{{s}}}=\phi(\mathsf{\mathbf{{{z}}}}+\mathsf{\mathbf{{{w}}}})$.
Each module uses the turbo principle (widely known in iterative decoding
literature) to pass the extrinsic messages to its adjacent module.

###  2.2 Evaluation of vignettes collection

Assessing the performance of the proposed symptom checker requires the
collection of a real-world test set containing the medical consultation
information, a.k.a., a “vignette”, of patients represented by:

  * •

The list of symptoms, $\mathcal{S}_{\textrm{obs}}$, of the patient described
by the medical consultation, which is different from the list of symptoms
reported by the dermatologist. This choice guarantees that the symptom vector
$\mathbf{s}$ takes into account the partial observability of the set of
manifested symptoms of the patient due to his/her limited medical knowledge.

  * •

The ground truth disease, $d_{i}^{*}$, of the patient confirmed by the
dermatologist after one or multiple medical consultations. Therefore, the
target disease vector $\mathbf{d}$ contains one non-zero element equal to 1 in
the $i$th position.

At the end of the process, under the patients’ consent, we collected 200
vignettes covering the 31 diseases, i.e., $\sim$ 7 vignettes per disease.

##  3 Approximate message passing optimization

Approximate message passing (AMP)-based computational techniques have gained a
lot of attention since their introduction within the compressed sensing
framework (Donoho et al. (2009)) to solve optimization problems involving a
linear mixing of the vector to be optimized, as in (1). In particular, vector
AMP (VAMP) (Rangan et al. (2019)) and generalized VAMP (G-VAMP) (Schniter et
al. (2016)) algorithms exhibit faster convergence for a broad (non i.i.d.)
class of sensing matrices that are right orthogonally-invariant, unlike the
generalized AMP (GAMP) algorithm (Rangan (2011)).

$\bm{g}_{\bm{\mathsf{d}}}\left(\mathbf{\widehat{d}},\gamma_{\bm{\mathsf{d}}}\right)=\frac{{\int}\bm{d}\,\,p_{\bm{\mathsf{d}}}\left(\bm{d}\right)\,\mathcal{N}\big{(}\bm{d};\mathbf{\widehat{d}},\gamma_{\bm{\mathsf{d}}}^{-1}\bm{I}\big{)}\,\textrm{d}\bm{d}}{{\int}p_{\bm{\mathsf{d}}}\left(\bm{d}\right)\,\mathcal{N}\big{(}\bm{d};\mathbf{\widehat{d}},\gamma_{\bm{\mathsf{d}}}^{-1}\bm{I}\big{)}\,\textrm{d}\bm{d}}$
(2)

$\bm{g}_{\bm{\mathsf{z}}}\left(\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}},\mathbf{s}\right)=\frac{{\int}\bm{z}\,\,p_{\bm{\mathsf{s}}|\bm{\mathsf{z}}}\left(\bm{s}|\bm{z}\right)\,\mathcal{N}\big{(}\bm{z};\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}}^{-1}\bm{I}\big{)}\,\textrm{d}\bm{z}}{{\int}p_{\bm{\mathsf{s}}|\bm{\mathsf{z}}}\left(\bm{s}|\bm{z}\right)\,\mathcal{N}\big{(}\bm{z};\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}}^{-1}\bm{I}\big{)}\,\textrm{d}\bm{z}}$
(3)

$\bm{f}_{\bm{\mathsf{r}}}\left(\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}},\mathbf{\widehat{d}},\gamma_{\bm{\mathsf{d}}}\right)=\frac{{\int}{\int}\bm{r}\,\,\mathcal{N}\big{(}\bm{d};\mathbf{\widehat{d}},\gamma_{\bm{\mathsf{d}}}^{-1}\bm{I}\big{)}\,\mathcal{N}\big{(}\bm{z};\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}}^{-1}\bm{I}\big{)}\,\mathcal{N}\big{(}\bm{z}|\bm{d};\mathbf{\bm{A}\,\bm{d}},\gamma_{\bm{\mathsf{z}}}^{-1}\bm{I}\big{)}\,\textrm{d}\bm{z}\,\textrm{d}\bm{d}}{{\int}{\int}\mathcal{N}\big{(}\bm{d};\mathbf{\widehat{d}},\gamma_{\bm{\mathsf{d}}}^{-1}\bm{I}\big{)}\,\mathcal{N}\big{(}\bm{z};\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}}^{-1}\bm{I}\big{)}\,\mathcal{N}\big{(}\bm{z}|\bm{d};\mathbf{\bm{A}\,\bm{d}},\gamma_{\bm{\mathsf{z}}}^{-1}\bm{I}\big{)}\,\textrm{d}\bm{z}\,\textrm{d}\bm{d}},~{}\text{with}~{}\bm{\mathsf{r}}=\\{\bm{\mathsf{d}},\bm{\mathsf{z}}\\}.\vspace{-0.7cm}$
(4)

One could also consider the use of alternating direction method of multipliers
(ADMM) (Boyd et al. (2011)), which is a penalty-based iterative algorithm for
solving constrained optimization problems. However, it is known that VAMP
converges faster than ADMM because of its automatic tuning of the augmented
Lagrangian parameter $\rho$ in ADMM (Manoel et al. (2018)). In fact, one of
the variable splitting strategies employed by ADMM, known as Douglas-Rachford
splitting (Douglas and Rachford (1956)), corresponds to VAMP by letting
$\rho=\gamma_{\bm{\mathsf{d}}_{\mathsf{p}}}^{+}/2$ where
$\gamma_{\bm{\mathsf{d}}_{\mathsf{p}}}^{+}$ is the precision estimate of
$\mathbf{d}$.

VAMP algorithm for symptom checkers: For better illustration, the block
diagram of the generalized VAMP (G-VAMP) algorithm and its use toward solving
the symptom checkers optimization problem in (1) is depicted in Fig. 2. There,
we show its different constituent blocks, namely the different denoisers as
they interact with the so-called LMMSE module. All the exchanged extrinsic
information is performed using expectation propagation (a.k.a., turbo)
principle which approximates the posterior messages by Gaussian distributions.
The algorithmic steps of the G-VAMP are summarized in Algorithm 1 where each
line number is colored similarly to the corresponding module in Fig. 2. The
subscript $t$ stands for the iteration index and subscripts $\mathsf{p}$ and
$\mathsf{e}$ are used to distinguish “posterior” and “extrinsic” variables,
respectively.

Algorithm 1 Generalized VAMP for symptom checkers

1:$\mathbf{Input:}$ observed symptom vector $\bm{s}\in\mathbb{R}^{M}$; medical
knowledge matrix $\mathbf{A}$; precision tolerance $\xi=10^{-6}$; number of
iterations $T_{\textrm{max}}$; noise precision $\gamma_{w}$.

2:$\mathbf{Initialize}$

3:$t\leftarrow 1$

4:$\triangleright$ COMMENTinitialize posterior and extrinsic means and
precisions

5:$\bm{\widehat{d}}^{-}_{\mathsf{e},0}$,
$\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},0}^{-}$,
$\bm{\widehat{z}}^{-}_{\mathsf{e},0}$,
$\gamma_{\bm{\mathsf{z}}_{\mathsf{e}},0}^{-}$

6:Repeat

7:$\triangleright$ COMMENTdenoising $\bm{d}$

8:$\widehat{\bm{d}}_{\mathsf{p},t}^{+}=\bm{g}_{\bm{\mathsf{d}}}\left(\mathbf{\widehat{d}}^{-}_{\mathsf{e},t-1},\gamma_{\bm{\mathsf{d}}_{\mathsf{e},t-1}}^{-}\right)$

9:$\alpha_{\bm{\mathsf{d}}_{\mathsf{p}},t}^{+}=\left\langle\bm{g}_{\bm{\mathsf{d}}}^{\prime}\left(\mathbf{\widehat{d}}^{-}_{\mathsf{e},t-1},\gamma_{\bm{\mathsf{d}}_{\mathsf{e},t-1}}^{-}\right)\right\rangle$

10:$\gamma_{\bm{\mathsf{d}}_{\mathsf{p}},t}^{+}=\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t-1}^{-}/\alpha_{\bm{\mathsf{d}}_{\mathsf{p}},t}^{+}$,
$\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{+}=\gamma_{\bm{\mathsf{d}}_{\mathsf{p}}}^{+}-\gamma_{\bm{\mathsf{d}}_{\mathsf{e},t-1}}^{-}$

11:$\mathbf{\widehat{d}}^{+}_{\mathsf{e},t}=\left(\gamma_{\bm{\mathsf{d}}_{\mathsf{p},t}}^{+}\,\widehat{\bm{d}}_{\mathsf{p},t}^{+}-\gamma_{\bm{\mathsf{d}}_{\mathsf{e},t-1}}^{-}\,\mathbf{\widehat{d}}^{-}_{\mathsf{e},t-1}\right)/\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{+}$

12:$\widehat{\bm{d}}_{\mathsf{p},t}^{-}=\bm{f}_{\bm{\mathsf{d}}}\left(\mathbf{\widehat{z}}^{-}_{\mathsf{e},t-1},\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t-1}}^{-},\mathbf{\widehat{d}}^{+}_{\mathsf{e},t},\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{+}\right)$

13:$\alpha_{\bm{\mathsf{d}}_{\mathsf{p}},t}^{-}=\left\langle\bm{f}_{\bm{\mathsf{d}}}\left(\mathbf{\widehat{z}}^{-}_{\mathsf{e},t-1},\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t-1}}^{-},\mathbf{\widehat{d}}^{+}_{\mathsf{e},t},\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{+}\right)\right\rangle$

14:$\gamma_{\bm{\mathsf{d}}_{\mathsf{p}},t}^{-}=\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{+}/\alpha_{\bm{\mathsf{d}}_{\mathsf{p}},t}^{-}$,
$\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{-}=\gamma_{\bm{\mathsf{d}}_{\mathsf{p}}}^{-}-\gamma_{\bm{\mathsf{d}}_{\mathsf{e},t}}^{+}$

15:$\mathbf{\widehat{d}}^{-}_{\mathsf{e},t}=\left(\gamma_{\bm{\mathsf{d}}_{\mathsf{p},t}}^{-}\,\widehat{\bm{d}}_{\mathsf{p},t}^{-}-\gamma_{\bm{\mathsf{d}}_{\mathsf{e},t}}^{+}\,\mathbf{\widehat{d}}^{+}_{\mathsf{e},t}\right)/\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{-}$

16:$\triangleright$ COMMENTdenoising $\bm{z}$

17:$\widehat{\bm{z}}_{\mathsf{p},t}^{+}=\bm{f}_{\bm{\mathsf{z}}}\left(\mathbf{\widehat{z}}^{-}_{\mathsf{e},t-1},\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t-1}}^{-},\mathbf{\widehat{d}}^{+}_{\mathsf{e},t},\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{+}\right)$

18:$\alpha_{\bm{\mathsf{z}}_{\mathsf{p}},t}^{+}=\left\langle\bm{f}_{\bm{\mathsf{z}}}^{\prime}\left(\mathbf{\widehat{z}}^{-}_{\mathsf{e},t-1},\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t-1}}^{-},\mathbf{\widehat{d}}^{+}_{\mathsf{e},t},\gamma_{\bm{\mathsf{d}}_{\mathsf{e}},t}^{+}\right)\right\rangle$

19:$\gamma_{\bm{\mathsf{z}}_{\mathsf{p}},t}^{+}=\gamma_{\bm{\mathsf{z}}_{\mathsf{e}},t-1}^{-}/\alpha_{\bm{\mathsf{z}}_{\mathsf{p}},t}^{+}$,
$\gamma_{\bm{\mathsf{z}}_{\mathsf{e}},t}^{+}=\gamma_{\bm{\mathsf{z}}_{\mathsf{p}}}^{+}-\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t-1}}^{-}$

20:$\mathbf{\widehat{z}}^{+}_{\mathsf{e},t}=\left(\gamma_{\bm{\mathsf{z}}_{\mathsf{p},t}}^{+}\,\widehat{\bm{z}}_{\mathsf{p},t}^{+}-\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t-1}}^{-}\,\mathbf{\widehat{z}}^{-}_{\mathsf{e},t-1}\right)/\gamma_{\bm{\mathsf{z}}_{\mathsf{e}},t}^{+}$

21:$\triangleright$ COMMENTEstimation of $\bm{z}$ from $\bm{s}$

22:$\widehat{\bm{z}}_{\mathsf{p},t}^{-}=\bm{g}_{\bm{\mathsf{z}}}\left(\mathbf{\widehat{z}}^{+}_{\mathsf{e},t},\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t}}^{+},\mathbf{s}\right)$

23:$\alpha_{\bm{\mathsf{z}}_{\mathsf{p}},t}^{-}=\left\langle\bm{g}_{\bm{\mathsf{z}}}^{\prime}\left(\mathbf{\widehat{z}}^{+}_{\mathsf{e},t},\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t}}^{+},\mathbf{s}\right)\right\rangle$

24:$\gamma_{\bm{\mathsf{z}}_{\mathsf{p}},t}^{-}=\gamma_{\bm{\mathsf{z}}_{\mathsf{e}},t}^{+}/\alpha_{\bm{\mathsf{z}}_{\mathsf{p}},t}^{-}$,
$\gamma_{\bm{\mathsf{z}}_{\mathsf{e}},t}^{-}=\gamma_{\bm{\mathsf{z}}_{\mathsf{p}}}^{-}-\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t-1}}^{+}$

25:$\mathbf{\widehat{z}}^{-}_{\mathsf{e},t}=\left(\gamma_{\bm{\mathsf{z}}_{\mathsf{p},t}}^{-}\,\widehat{\bm{z}}_{\mathsf{p},t}^{-}-\gamma_{\bm{\mathsf{z}}_{\mathsf{e},t}}^{+}\,\mathbf{\widehat{z}}^{+}_{\mathsf{e},t}\right)/\gamma_{\bm{\mathsf{z}}_{\mathsf{e}},t}^{-}$

26:$t\leftarrow t+1$

27:Until
$\Big{(}\big{|}\\!\big{|}\bm{\widehat{d}}_{\mathsf{p},t+1}^{+}-\bm{\widehat{d}}_{\mathsf{p},t}^{+}\big{|}\\!\big{|}^{2}_{\textrm{F}}\Big{)}$
$\leq\xi\Big{(}\big{|}\\!\big{|}\bm{\widehat{d}}_{\mathsf{p},t}^{+}\big{|}\\!\big{|}^{2}_{\textrm{F}}$)
or ($t>T_{\textrm{max}}\Big{)}$

28:return $\widehat{\bm{d}}_{\mathsf{p},T_{\textrm{max}}+1}^{+}$

Denoising $\bm{d}$ and $\bm{z}$: The denoising functions of $\mathbf{d}$ and
$\mathbf{z}$ in lines 8 and 17 are given in (2) and (4), respectively. The
function $\bm{g}_{\bm{\mathsf{d}}}\left(\cdot,\cdot\right)$ in (2) acts as a
“denoiser” of the additive white Gaussian noise-corrupted pseudo-measurement
$\mathbf{\widehat{d}}^{-}_{\mathsf{e}}=\bm{d}^{*}+\mathbf{w}$ where
$\mathbf{w}\sim\mathcal{N}(\bm{w};\bm{0},\bm{I}/\gamma_{\bm{\mathsf{d}}_{\mathsf{e}}}^{-})$,
using the prior knowledge on the disease ground truth $\bm{d}^{*}$. The
denoising bloc in purple in Fig. 2 denoises both $\mathbf{d}$ and $\mathbf{z}$
under the pseudo-prior
$\mathcal{N}(\bm{d};\mathbf{\widehat{d}}^{+}_{\mathsf{e}},\bm{I}/\gamma_{\bm{\mathsf{d}}_{\mathsf{e}}}^{+})$
and
$\mathcal{N}(\bm{z};\mathbf{\widehat{z}}^{-}_{\mathsf{e}},\bm{I}/\gamma_{\bm{\mathsf{z}}_{\mathsf{e}}}^{-})$,
respectively, using the function
$\bm{f}_{\bm{\mathsf{r}}}\left(\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}},\mathbf{\widehat{d}},\gamma_{\bm{\mathsf{d}}}\right)$
given in (4) where $\bm{\mathsf{r}}=\\{\bm{\mathsf{d}},\bm{\mathsf{z}}\\}$.

Estimating $\bm{z}$ from symptoms $\bm{s}$: The fact that the generative model
of the observed symptom vector $s$ has the form
$\mathsf{\mathbf{{s}}}=\phi(\mathsf{\mathbf{{{z}}}}+\mathsf{\mathbf{{{w}}}})$
allows the calculation of the posterior mean
$\mathbf{\widehat{z}}^{-}_{\mathsf{p}}$ using the function
$\bm{g}_{\bm{\mathsf{z}}}(\cdot)$ given in (3). There, the density
$p_{\bm{\mathsf{s}}|\bm{\mathsf{z}}}\left(\bm{s}|\bm{z}\right)$ represents the
non-linear element-wise channel through which the dense symptom vector
$\bm{z}$ is transformed in a sparse symptom vector $\bm{s}$. For binary
symptom checkers, we let $\phi(x)=\textrm{sign}(x)$ which corresponds to the
following denoising function:

$\bm{g}_{\bm{\mathsf{z}}}\left(\mathbf{\widehat{z}},\gamma_{\bm{\mathsf{z}}},\mathbf{s}\right)=\begin{cases}\frac{v}{u}+\frac{u}{\gamma_{\bm{\mathsf{z}}}}\cdot\frac{\phi(v)}{\Phi(v)}&,\textrm{if}\quad
s_{i}=1\\\
\frac{v}{u}-\frac{u}{\gamma_{\bm{\mathsf{z}}}}\cdot\frac{\phi(v)}{1-\Phi(v)}&,\textrm{if}\quad
s_{i}=-1\end{cases}$ (5) where $\displaystyle
u=\frac{1}{\sqrt{\gamma_{\bm{\mathsf{z}}}^{-2}+\gamma_{w}^{-2}}},\quad
v=\mathbf{\widehat{z}}\cdot u$
$\displaystyle\phi(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}x^{2}},\quad\Phi(x)=\int_{-\infty}^{x}\phi(t)\,\textrm{d}t$

##  4 Evaluation & Results

We evaluate the NRMSE performance of G-VAMP (Rangan et al. (2019)) on the
$P=200$ collected vignettes. We benchmark it with other optimization
techniques, namely: $i)$ sparse sensing (Rozenberg et al. (2018)) which works
by scanning all possible support bases for the disease vector $\bm{d}$ when it
has either 1 or 2 non-zero elements, $ii)$ ADMM (Boyd et al. (2011)) for a
range values of augmented Lagrangian penalty $\rho\in[0,5]$ with a step size
of 0.1, and $iii)$ the standrad underdetermined least square (ULS) closed-form
solution (Datta (2010)). We also evaluate the symptom checker using the top-K
accuracy metric for $K\in\\{1,2,3\\}$, which is well known in computer vision
classification scenarios to successfully assess the differential diagnosis
cases. All G-VAMP experiments were run with a signal-to-noise ratio (SNR) of
25 dB.

Results are reported in Table 1. We observe that G-VAMP outperforms the three
baseline algorithms significantly. While ADMM exhibits a performance close to
G-VAMP, this has been achieved by reporting its best performance among all 51
experiments by varying its parameter $\rho\in[0,5]$ with a step size of 0.1.
Clearly, the auto-tuning feature of G-VAMP not only provides a better
performance but also reduces the time complexity of symptom checkers which
must provide a real-time result for the end user.

Table 1: NRMSE and top-K metrics of the different algorithms on the binary symptom checker problem (1). | ULS | Sparse sensing | ADMM | G-VAMP
---|---|---|---|---
NRMSE | 0.70 | 0.56 | 0.53 | 0.44
TOP 1 | 0.39 | - | 0.78 | 0.81
TOP 2 | 0.41 | - | 0.81 | 0.83
TOP 3 | 0.43 | - | 0.84 | 0.87

##  5 Future work & Conclusion

To simplify the medical knowledge data collection, we have formulated binary
symptom checkers as an optimization problem and advocated AMP algorithms to
solve it by showing a significant improvement in terms of NRMSE and top-K
accuracy performance. Many avenues for further investigation remain open. When
the prior knowledge of the disease density on $\bm{d}$ is not available, it is
possible to relax this requirement using data-driven methods to infer the
prior from the available datasets. It is also possible to substitute the
denoising function $\mathbf{g}_{\bm{\mathsf{d}}}(\cdot)$ with a separate
optimization procedure via Plug-and-Play priors (Venkatakrishnan et al.
(2013)). One can also investigate the bilinear problem version of (1) to
jointly recover both the disease vector $\mathbf{d}$ and the medical knowledge
matrix $\mathbf{A}$ by employing bilinear message passing algorithms such as
BiG-VAMP (Akrout et al. (2020)) and BiG-AMP (Parker et al. (2014)). This will
enable the comparison of the reconstructed matrix $\mathbf{A}$ with the one
designed by medical experts.

## References

  * Akrout et al. (2019) Mohamed Akrout, Amir-massoud Farahmand, Tory Jarmain, and Latif Abid.  Improving skin condition classification with a visual symptom checker trained using reinforcement learning.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 549–557. Springer, 2019.
  * Akrout et al. (2020) Mohamed Akrout, Anis Housseini, Faouzi Bellili, and Amine Mezghani.  Bilinear generalized vector approximate message passing.  _arXiv preprint arXiv:2009.06854_ , 2020.
  * Boyd et al. (2011) Stephen Boyd, Neal Parikh, and Eric Chu.  _Distributed optimization and statistical learning via the alternating direction method of multipliers_.  Now Publishers Inc, 2011.
  * Choi et al. (2016) Edward Choi, Mohammad Taha Bahadori, Andy Schuetz, Walter F. Stewart, and Jimeng Sun.  Doctor ai: Predicting clinical events via recurrent neural networks.  In _Proceedings of the 1st Machine Learning for Healthcare Conference_ , volume 56 of _Proceedings of Machine Learning Research_ , pages 301–318, Northeastern University, Boston, MA, USA, 18–19 Aug 2016. PMLR.
  * Coleman (2017) Cathy Coleman.  Early detection and screening for breast cancer.  In _Seminars in oncology nursing_ , volume 33, pages 141–155. Elsevier, 2017.
  * Datta (2010) Biswa Nath Datta.  _Numerical linear algebra and applications_ , volume 116.  Siam, 2010.
  * Donoho et al. (2009) David L Donoho, Arian Maleki, and Andrea Montanari.  Message-passing algorithms for compressed sensing.  _Proceedings of the National Academy of Sciences_ , 106(45):18914–18919, 2009.
  * Douglas and Rachford (1956) Jim Douglas and Henry H Rachford.  On the numerical solution of heat conduction problems in two and three space variables.  _Transactions of the American mathematical Society_ , 82(2):421–439, 1956.
  * Fraser et al. (2017) Hamish SF Fraser, Susan Clamp, and Carol J Wilson.  Limitations of study on symptom checkers.  _JAMA internal medicine_ , 177(5):740–741, 2017\.
  * Manoel et al. (2018) Andre Manoel, Florent Krzakala, Gaël Varoquaux, Bertrand Thirion, and Lenka Zdeborová.  Approximate message-passing for convex optimization with non-separable penalties.  _arXiv preprint arXiv:1809.06304_ , 2018.
  * Parker et al. (2014) Jason T Parker, Philip Schniter, and Volkan Cevher.  Bilinear generalized approximate message passing—part i: Derivation.  _IEEE Transactions on Signal Processing_ , 62(22):5839–5853, 2014.
  * Rangan (2011) Sundeep Rangan.  Generalized approximate message passing for estimation with random linear mixing.  In _2011 IEEE International Symposium on Information Theory Proceedings_ , pages 2168–2172. IEEE, 2011.
  * Rangan et al. (2019) Sundeep Rangan, Philip Schniter, and Alyson K Fletcher.  Vector approximate message passing.  _IEEE Transactions on Information Theory_ , 65(10):6664–6684, 2019.
  * Rotmensch et al. (2017) Maya Rotmensch, Yoni Halpern, Abdulhakim Tlimat, Steven Horng, and David Sontag.  Learning a health knowledge graph from electronic medical records.  _Scientific reports_ , 7(1):1–11, 2017.
  * Rozenberg et al. (2018) Igal Rozenberg, Yuval Beck, Yonina C Eldar, and Yoash Levron.  Sparse estimation of faults by compressed sensing with structural constraints.  _IEEE Transactions on Power Systems_ , 33(6):5935–5944, 2018.
  * Schniter et al. (2016) Philip Schniter, Sundeep Rangan, and Alyson K Fletcher.  Vector approximate message passing for the generalized linear model.  In _2016 50th Asilomar Conference on Signals, Systems and Computers_ , pages 1525–1529. IEEE, 2016.
  * Venkatakrishnan et al. (2013) Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg.  Plug-and-play priors for model based reconstruction.  In _2013 IEEE Global Conference on Signal and Information Processing_ , pages 945–948. IEEE, 2013.