# Optimized experimental optical tomography of quantum states of room-
temperature alkali-metal vapor

Marek Kopciuch 1,2  Magdalena Smolis2  Adam Miranowicz3  Szymon Pustelny2
1Doctoral School of Exact and Natural Sciences, Jagiellonian University,
Faculty of Physics, Astronomy and Applied Computer Sciences, Łojasiewicza 11,
30-348 Kraków, Poland
2Institute of Physics, Jagiellonian University in Kraków, Łojasiewicza 11,
30-348 Kraków, Poland
3Institute of Spintronics and Quantum Information, Faculty of Physics, Adam
Mickiewicz University, 61-614 Poznań, Poland
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We demonstrate a novel experimental technique for quantum-state tomography of
the collective density matrix. It is based on measurements of the polarization
of light, traversing the atomic vapor. To assess the technique’s robustness
against errors, experimental investigations are supported with numerical
simulations. This not only allows to determine the fidelity of the
reconstruction, but also to analyze the quality of the reconstruction for
specific experimental parameters (light tuning and number of measurements). By
utilizing the so-called conditional number, we demonstrate that the
reconstruction can be optimized for a specific tuning of the system
parameters, and further improvement is possible by selective repetition of the
measurements. Our results underscore the potential high-fidelity quantum-state
reconstruction while optimizing measurement resources.

††journal: opticajournal††articletype: Research Article

##  1 Introduction

Quantum technology is built on the precise manipulation and reconstruction of
quantum states. When dealing with single microscopic quantum objects, the
reconstruction of states becomes challenging. This stems from a (often)
destructive nature of the reconstruction and small amplitudes of recorded
signals. To address these difficulties, some researchers have turned their
focus towards studying ensembles of quantum objects, which display a
collective quantum behavior [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15].

Atomic vapors serve as a prime example of a medium utilized for the
engineering of collective quantum states. In their ultracold form, they allow
for precise quantum control through light and other external fields, albeit
implementation of the control requires complex experimental setups. On the
other hand, room-temperature vapors can be studied using simpler apparatuses,
but they simultaneously present challenges in terms of theoretical
understanding [3]. Despite these problems, however, the room-temperature
atomic vapors were used to demonstrate various quantum-mechanical effects
including coherent population trapping [4], spin squeezing [5, 6], macroscopic
entanglement [7, 8], spin waves [8, 9], squeezed light generation [10, 11, 5]
and entanglement of light modes[12]. Rubidium vapor was also used to construct
an on-demand quantum memory [13, 14]. These experiments revived the interest
in such media, while also necessitated the development of reliable quantum-
state tomography (QST) methods.

In this work, we demonstrate the first experimental implementation of recently
proposed QST method [15]. The method enables the reconstruction of a
collective density matrix of a room-temperature atomic vapor and is based on
the illumination of the vapor with an off-resonant probing light and
monitoring properties of the light after traversing a medium subjected to an
external magnetic field. This enables to reconstruct a collective quantum
state of 87Rb atoms residing in the $F=1$ ground state (qutrit).

To evaluate the efficiency of the tomographic technique, we used the so-called
conditional number [16, 17, 18]. Previously, the parameter was used for a
comprehensive comparison of tomographic methods of two polarization qubits
[17], NMR tomography of two ${}^{1}H$ spins-1/2 (two qubits) [19], and a
single nuclear spin-3/2 (a quartit) in a semiconductor quantum well [20]. We
demonstrate that by an appropriate tuning of the probing light, the
conditional number can be minimaized (corresponding to an optimized
reconstruction) and as small as 2.25 can be achieved. We also discuss means of
further improvement of the reconstruction efficiency by repeating specific
measurements.

##  2 Principles of the optical tomography

We begin with a brief overview of the QST technique developed in Ref. [15].
This method relies on measuring the polarization rotation of linearly
polarized probe light traversing a medium (e.g., room-temperature alkali metal
atoms) subjected to a longitudinal magnetic field. We assume that the
amplitude of the light is low, which allows us to describe its interaction
with atoms using perturbation theory at the lowest order (linear interaction).
At the same time, unlike previous approaches (see, e.g. Ref. [21, 22]), we do
not assume a significant detuning of the light from the optical transition.
This enables us to consider not only the vector contributions to a
polarization rotation [23, 21, 24], but also tensor one [25], and hence
reconstruct the collective density matrix of the atoms. It is noteworthy that
this reconstruction is achieved without full control over the system, as
successive magnetic sublevels are equally splitted due to a weak magnetic
field (under the conditions of the linear Zeeman effect) [26].

In Ref. [15], the relation between time-dependent polarization rotation
$\delta\alpha(t)$ and operators $\hat{\alpha}_{R,I}$ and $\hat{\beta}$ was
introduced. The operators are associated with coherences and population
difference of specific magnetic sublevels hence provide access to specific
density-matrix elements. In this work, we employ a slightly modified version
of that relationship, i.e.,

$\delta\alpha(t;\Delta)=\eta(\Delta)\left(e^{-\gamma_{1}t}\left[\expectationvalue{\hat{\alpha}_{R}}\sin(2\Omega_{L}t)+\expectationvalue{\hat{\alpha}_{I}}\cos(2\Omega_{L}t)\right]-\zeta(\Delta)e^{-\gamma_{2}t}\expectationvalue{\hat{\beta}}\right),$
(1)

where $\eta(\Delta)=\chi V_{R}(\Delta)$ and
$\zeta(\Delta)=V_{I}(\Delta)/V_{R}(\Delta)$ are the so-called global and local
scaling factors associated with real $V_{R}$ and imaginary $V_{I}$ parts of
the Voigt profile, and $\chi$ is related to experimental parameters such as
atomic density and transition frequency (for more details see the Supplemental
Information – SI). As shown in Eq. (1), the time dependence of the
polarization rotation is determined by the Larmor frequency $\Omega_{L}$ and
the relaxation rates $\gamma_{1}$ and $\gamma_{2}$.

Since a single measurement described by Eq. (1) allows us to extract only
limited information about the system (specifically the population difference
and coherence between magnetic sublevels with $\Delta m_{F}=2$) it is
necessary to expand the set of measured signals to obtain a more comprehensive
information. To achieve this, we introduce a series of unitary operations
known as control pulses, which systematically manipulate a given state in the
Hilbert space. This provides the access to other density-matrix elements and
hence offer a complete characterization of the system [15]. In turn, the
reconstruction problem can be presented as

$\mathbb{O}\rho_{V}=\mathbf{b},$ (2)

where $\mathbb{O}$ represents the coefficient matrix determined by the set of
observables, and
$\rho_{V}=\left[\rho_{\bar{1}\bar{1}}^{R},\rho_{\bar{1}0}^{R},\rho_{\bar{1}0}^{I},\ldots\right]^{T}$
(where $\rho_{mn}^{R}={\rm Re}\\{\rho_{mn}\\}$, $\rho_{mn}^{I}={\rm
Im}\\{\rho_{mn}\\}$ and $\bar{1}=-1$) is the vectorized form of a standard
form density matrix $\rho$ with entrances $\rho_{ij}$ (see SI for more
information), $\mathbf{b}$ is the observation vector and contains the values
of the measured values of the observables. In a typical experimental scenario,
the set of measurements given in Eq. (2) is often overdetermined, and it is
advantageous to rescale it to a more suitable form

$\mathbb{C}\rho_{V}=\tilde{\mathbf{b}},$ (3)

where $\mathbb{C}=\mathbb{O}^{\dagger}\mathbb{O}$ and
$\tilde{\mathbf{b}}=\mathbb{O}^{\dagger}\mathbf{b}$. This rescaling enables
the calculation of the density operator by simply inverting the aforementioned
linear problem.

##  3 Experimental details

###  3.1 Experimental setup

|
---|---
(a) | (b)
Figure 1: (a) Simplified scheme of the experimental setup used for the
quantum-state generation and tomography. SAS – saturation absorption
spectroscopy, DAVLL – dichroic atomic vapor laser lock, AOM – acusto-optic
modulator, PD – photodiode, Pol – Glan-Thompson polarizer, $\lambda/4$ –
quater-wave plate, 87Rb – parafin-coated vapor cell filled with 87Rb, Wol –
Wollaston prism, BPD - balanced photodetector. (b) Experimental sequence used
in our method. The initial state of atoms is prepared with the pump light
turned on for about 200 ms (red trace). After the preparation period, we apply
a sequence of magnetic pulses to modify the state of the atoms (green trace).
Here, the CYCLOPS pulses (see Sec. 3.4) are first used and then the control
pulses are implemented for a total time of about 1.5 ms. Finally, the probe
light is turned on, alongside with the longitudinal magnetic field, for about
1 s (blue) and the polarization rotation signal is recorded.

The heart of our experimental system is a 3 cm diameter paraffin-coated
spherical cell, containing an istopically enriched sample of 87Rb atoms. The
cell is heated up to $50^{\circ}$C and is placed inside a cylindrical magnetic
shield made of three layers of mumetal and a qubic innermost ferrite layer.
Apart from the cell, the shield additionally contains a set of magnetic-field
coils, which enables residual-field compensation and generation of magnetic-
field pulses in the $\mathbf{x}$, $\mathbf{y}$, and $\mathbf{z}$ directions.
Light used for the illumination of the rubidium atoms is provided by three
diode lasers, where the pump and probe lasers are distributed-feedback lasers
(DFBs), and the repump laser is the Fabry-Perot laser (ECDL). All lasers are
independently tuned, and the repump laser wavelength is frequency-stabilized
using a Dichroic Atomic Vapor Laser Lock (DAVLL) [27]. The wavelengths of the
other two lasers are passively maintained due to their inherent temporal
stability. Performance of all lasers is monitored using a wavemeter, while the
pump and probe lasers are additionally monitored through saturated absorption
spectroscopy (SAS). The intensities of the laser beams are dynamically
controlled by three acusto-optical modulators (AOMs). To generate a specific
quantum state in the vapor, the pump-light polarization is set by polarizers
(POLs) and quarter-wave plates ($\lambda/4$), while the repump light is
linearly polarized orthogonal to the pump-light propagation direction (i.e.,
along $\mathbf{y}$-axis). To determine the local scaling factor (see
discussion below), the intensity of the probe light is monitored and its
$\mathbf{y}$ linear polarization prior to the shield is provided by a Glan-
Thomson polarizer. Finally, the polarization rotation of the probe light is
measured after the cell using a balanced polarimater consisting of a
Wollastron prism (WOL) and a balanced photodetector (BPD). The schematic of
the setup is shown in Fig. 1(a).

###  3.2 Experimental sequence

The experimental sequence utilized in our measurements is shown in Fig. 1(b).
The sequence begins with a pumping period during which a specific quantum
state is engineered. This stage typically consists of a 200 ms light pulse
(optical pumping), which is applied simultaneously with the repumping that
prevent the atoms from escaping into the dark ($F=2$) state, followed by a few
short ($\approx 100$ $\mu$s) magnetic-field pulses, enabling generation of a
desired complex state. Subsequently, a series of magnetic-field pulses is used
to mitigate technical problems (see Sec. 3.4), which is followed by a set of
control pulses. Once the pulses are completed, a constant magnetic field along
the $\mathbf{z}$-direction, ranging from 10-100 nT, is established. At the
same time, a probe light beam, propagating along $\mathbf{z}$ with intensity
of 1-10 $\mu$W/cm2, is turned on. In order to improve the signal-to-noise
ratio, the intensity of the probe light is modulated at a frequency of 200 kHz
and the polarimeter signal is detected using a lock–in amplifier.

###  3.3 Global and local scaling factor

An important element of the reconstruction of the density matrix is the
determination of the global scaling factor $\eta(\Delta)$ [see Eq. (1)]. This
can be done by measuring the light absorption in an unpolarized vapor. Using
the absorption relationship derived in the SI, the factor can be identified by
comparing the absorption of the probe light, tuned to the same wavelength as
that during the tomography measurements (i.e., blue-detuned from
$f=1\rightarrow F=2$ by 50–400 MHz), with the absorption of far-detuned light
(>15 GHz).

$\eta(\Delta)=\dfrac{27}{16}\left(\sqrt{\frac{U_{2}(\Delta)/U_{1}(\Delta)}{U_{2}(\infty)/U_{1}(\infty)}}-1\right),$
(4)

where $U_{1}$ is the voltage measured at the transimpedance photodetector
placed in front of and $U_{2}$ after the medium (see Fig. 1(a) and the SI for
more details) with $\Delta$ indicating the probe light tuned for QST and
$\infty$ far-detuned light.

Experimental determination of the local scaling factor $\zeta(\Delta)$ [see
Eq. (1)] presents a greater challenge. It requires preparation of an
anisotropic, yet well-defined quantum state. In this work, we select
“stretched” states that are generated along the $\mathbf{x}$\- and
$\mathbf{z}$-axes. The first state can be created by illuminating the atoms
with a circularly polarized pump light propagating along the
$\mathbf{x}$-axis. The preparation of the second state is more involved and
requires the application of an additional magnetic-field pulse after the
pumping, which rotates the atomic $\mathbf{x}$-polarization to the
$\mathbf{z}$-direction (we have experimentally verified that this process did
not introduce dephasing, as evidenced by the unchanged signal amplitude for a
many-$\pi$ pulse). Employing this procedure allows us to mitigate potential
systematic errors arising from varying polarization levels achieved with the
pump light propagating along different directions, while simultaneously
simplifying the experimental setup. The formulas for the light polarization
rotation corresponding to these two states are (see the SI for more details)

$\displaystyle\delta\alpha^{(z)}(t;\Delta)$ $\displaystyle=$
$\displaystyle-\dfrac{5(1-\epsilon)}{24}\eta(\Delta)\zeta(\Delta)e^{-\gamma_{2}t},$
(5a) $\displaystyle\delta\alpha^{(x)}(t;\Delta)$ $\displaystyle=$
$\displaystyle-\dfrac{(1-\epsilon)}{48}\eta(\Delta)e^{-\gamma_{1}t}\cos(2\Omega_{L}t),$
(5b)

where $\epsilon$ is the remaining isotropic part of state. This allows one to
calculate the local scaling factor

$\zeta(\Delta)=\dfrac{1}{10}\dfrac{\delta\alpha^{(z)}(0;\Delta)}{\delta\alpha^{(x)}(0;\Delta)}.$
(6)

###  3.4 CYCLOPS-like measurement

Equation (1) shows that our reconstruction method is sensitive to the initial
phase of the measured signal. As uncontrollable phase delays are present in
every experiment, the identification of the quadrature components of the
signal becomes difficult. To address this issue, we adapt the CYCLically
Ordered Phase Sequence (CYCLOPS) method, commonly utilized in nuclear magnetic
resonance experiments [28, 29]. In our approach, we leverage the fact that the
$\pi$-rotation of the state around the $\mathbf{y}$-axis leads to a sign
reversal of $\expectationvalue{\hat{\alpha}_{I}}$ and
$\expectationvalue{\hat{\beta}}$ (for more information, see the SI). At the
same time, by applying the pulse rotating the state by $\pi/2$ around the
$\mathbf{z}$-axis and next the pulse rotating the state around the
$\mathbf{y}$-axis by $\pi$ (see the SI) the signs of the
$\expectationvalue{\hat{\alpha}_{R}}$ and $\expectationvalue{\hat{\beta}}$ are
reversed. By subtracting these two transformed states from the initial signal,
we obtain

$\displaystyle\left(\delta\alpha-\delta\alpha^{(Y)}\right)(t;\Delta)=2\eta(\Delta)\left[-\zeta(\Delta)e^{\gamma_{2}t}\expectationvalue{\hat{\beta}}+e^{-\gamma_{1}t}\expectationvalue{\hat{\alpha}_{I}}\cos(2\Omega_{L}t+\varphi)\right],$
(7a)
$\displaystyle\left(\delta\alpha-\delta\alpha^{(ZY)}\right)(t;\Delta)=2\eta(\Delta)\left[-\zeta(\Delta)e^{\gamma_{2}t}\expectationvalue{\hat{\beta}}+e^{-\gamma_{1}t}\expectationvalue{\hat{\alpha}_{R}}\sin(2\Omega_{L}t+\varphi)\right],$
(7b)

where $\varphi$ is an unknown phase shift originating from the experimental
apparatus. In our CYCLOPS-like measurements, the problem of unknown phase is
alleviated, as the final signals [see Eqs. (7)] depend only on one quadrature
(via either sine or cosine time dependence) and, thus, $\varphi$ becomes
insignificant. The procedure also allows us to remove systematic shifts of the
signals associated with the imbalance of the polarimeter (for more details,
see the SI).

##  4 Reconstruction of states

To perform QST, we conducted the above-described nine measurements, consisting
of three sets of CYCLOPS-like pulses for each of three control pulses. To
ensure the self-consistency of our reconstruction procedure, we simultaneously
fit all of the polarization-rotation signals with shared parameters such as
the global phase, relaxation rates, and oscillation frequency. The fitting
values are then used to determine the observables and reconstruct the qutrit
density-matrix elements using the linear inversion method given in Eq. (3).
However, as this method does not guarantee the reconstructed matrices to be
positive semidefinite, we utilize the maximum likelihood method with the
Euclidean norm [20, 30] to find the closest physical realization of the
reconstructed matrix.

To validate our tomography technique, we compare the reconstructed density
matrices with numerical simulations of the state obtained during the pumping
stage. For the simulations, we assume the interaction of an appropriately
polarized light with a Doppler-broadened medium consisting of atoms of the
energy-level structure similar to that of the $D_{1}$ line in 87Rb. As in the
real experiment, we assume that there are two distinct regions between which
the atoms can freely move. In the first region, the atoms evolve in a
homogeneous magnetic field and relax to thermal equilibrium due to the
collisions with vapor-cell walls and between one another. This corresponds to
the atoms residing outside of the light beams. In the second region, the atoms
still interact with the magnetic field but also with the pump and repump
light. Moreover, we neglect the wall relaxation in this region. The latter
region corresponds to the atoms inside the light beams. All parameters used in
the simulations match the parameters of our experimental setup.

|  |
---|---|---
(a) |  | (b)
Figure 2: Comparison of two experimentally reconstructed density matrix
elements (blue bars), with simulated ones (red bars), including their
amplitude (upper plots) and phase (lower plots). Here we chose simple pumping
schemes with (a) circularly polarized and (b) linearly $\mathbf{z}$–polarized
pumps propagating along $\mathbf{x}$–axis. The fidelity achieved between the
experimental results and simulations in both cases exceeds 0.99.

As representative examples for our reconstruction, we consider two states that
can be easily generated experimentally and simulated theoretically. The first
state can be pumped with a strong, circularly polarized pumping light,
propagating along the $\mathbf{x}$-axis [Fig. 2(a)]. The state has a
nonuniform population distribution and its all coherences are nonzero. This
allows us to demonstrate that our method can reconstruct not only different
coherences but also determine their amplitudes and phases with a high
accuracy. The results of the experimental reconstruction and simulations are
presented in Fig. 2(a). As seen, the results are in a very good agreement
revealing a reconstruction fidelity of 0.995. As the second example, we
considered a state pumped with the $\pi$-polarized light, propagating along
the $\mathbf{x}$-axis. In the ideal case (without experimental artefacts),
this scheme leads to the total depletion of the $m_{F}=\pm 1$ states and no
coherences between any sublevels. As shown in Fig. 2(b), our measurements
demonstrate a good agreement with numerical simulations, demonstrating the
fidelity of 0.998. Nonetheless one can notice that a very small amplitude of
the coherences can lead to the deterioration of the phase reconstruction. The
very high quality of the reconstruction of these two representative states
demonstrates the usefulness of our QST technique.

##  5 Conditioning and optimization of quantum state tomography

###  5.1 Condition number in linear inversion

As mentioned above, the condition number $\kappa$ is a useful parameter to
evaluate the reliability of a QST method [see Eq. (3)]. Specifically, to
quantify the ability to tolerate errors or sensitivity to the errors, we use
the condition number of a (nonsingular) matrix $\mathbb{C}$, which, assuming
the spectral norm $\|\dots\|_{2}$, can be defined as [31, 32, 33]

$\kappa(\mathbb{C})=\|\mathbb{C}\|_{2}\;\|\mathbb{C}^{-1}\|_{2}=\max[{\rm
svd}(\mathbb{C})]\max[{\rm svd}(\mathbb{C}^{-1})]=\dfrac{\max[{\rm
svd}(\mathbb{C})]}{\min[{\rm svd}(\mathbb{C})]}\geq 1,$ (8)

where ${\rm svd}(\mathbb{C})$ denotes the singular values of $\mathbb{C}$. The
significance of this error-robustness parameter explains well the Gastinel-
Kahan theorem [32], which states that a relative distance of a nonsingular
square matrix $\mathbb{C}$ from the set of singular matrices corresponds to
the inverse of a condition number. Utilizing the error
$\delta\tilde{\mathbf{b}}$ in the observation vector $\tilde{\mathbf{b}}$ and
the condition number $\kappa(\mathbb{C)}$, one can estimate the error
$\delta\rho_{V}$ in the reconstructed density matrix $\rho_{V}$ from the so-
called Atkinson inequalities [31]

$\dfrac{1}{\kappa(\mathbb{C})}\dfrac{\|\delta\tilde{\mathbf{b}}\|}{\|\tilde{\mathbf{b}}\|}\leq\dfrac{\|\delta\rho_{V}\|}{\|\rho_{V}\|}\leq\kappa(\mathbb{C})\dfrac{\|\delta\tilde{\mathbf{b}}\|}{\|\tilde{\mathbf{b}}\|}.$
(9)

When the condition number approaches 1, it becomes apparent that small
relative variations in the observation vector $\tilde{\mathbf{b}}$ result in
correspondingly small relative changes in the reconstructed state $\rho_{V}$.
In order to account for errors $\delta\mathbb{C}$ present in the coefficient
matrix $\mathbb{C}$, these inequalities can be expanded according to the
formulation derived in Ref. [31], giving rise to the expression

$\dfrac{\|\delta\rho_{V}\|}{\|\rho_{V}\|}\leq\dfrac{\kappa(\mathbb{C})}{1-\kappa(\mathbb{C})\|\delta\mathbb{C}\|/\|\mathbb{C}\|}\left[\dfrac{\|\delta\tilde{\mathbf{b}}\|}{\|\tilde{\mathbf{b}}\|}+\dfrac{\|\delta\mathbb{C}\|}{\|\mathbb{C}\|}\right].$
(10)

By referring to the inequalities in Eqs. (9) and (10), we can infer that the
quality of a QST method, in terms of its error sensitivity or robustness, can
be assessed through its condition number $\kappa(\mathbb{C})$, which
characterizes the degree to which small (large) changes in the observation
vector $\tilde{\mathbf{b}}$ lead to relatively small (large) changes in the
reconstructed state $\rho_{V}$. Thus, if $\kappa(\mathbb{C})$ is small
(large), the QST method is well-conditioned (ill-conditioned), indicating the
robustness (sensitivity) of the method to errors in the observation vector
$\tilde{\mathbf{b}}$. In the case of ill-conditioned QST, even slight errors
in $\tilde{\mathbf{b}}$ can cause significant errors in the reconstructed
$\rho_{V}$. In short, the smaller the condition number the stronger robustness
of a given linear-inversion-based QST method against errors. Thus, one can
refer to an optimal method in this respect if $\kappa(\mathbb{C})=1$.
Numerical examples of ill-conditioned QST problems can be found in Refs. [31,
17].

###  5.2 Optimization via probe light tuning

In order to optimize a QST process, it is desired to make the coefficient
matrix $\mathbb{C}$ more isotropic, which means that each measurement brings
an equal amount of information about the system. A simple example of such an
optimized problem is when each measurement brings information about only a
specific density-matrix element, with all measurements having the same weight
[17, 18]. In this case, the coefficient matrix $\mathbb{C}$ is proportional to
identity. Even though such optimization is intuitive, it is often unpractical,
as experimental transformations required to achieve a desired scheme are very
complex. Instead, here we propose a scheme, where a single experimental
parameter is adjusted. In our case, this parameter is the probing light
detuning, which, incorporated in Eq. (1) through $\zeta(\Delta)$, makes one of
the observable detuning dependent.

It is important to note that our method does not guarantee an optimal
tomography process, $\kappa(\mathbb{C})=1$. Therefore, to explore the limit of
the method, we calculate the eigenvalues of coefficient matrix with
$\zeta(\Delta)$ as a free parameter. In our case, the eigenvalues of
$\mathbb{C}$ can be analytically calculated, taking the values:
$\left\\{\dfrac{1}{100},\dfrac{1}{150},\dfrac{1}{225},\dfrac{1}{225},\dfrac{1}{225},\dfrac{\zeta^{2}}{18},\dfrac{\zeta^{2}}{9},\dfrac{\zeta^{2}}{9}\right\\}$.
From this, we obtain the dependence of $\kappa(\mathbb{C})$ on $\zeta(\Delta)$
[Fig. 3(a)] and a minimal possible conditional number of 2.25 is determined.

|
---|---
(a) | (b)
Figure 3: (a) Condition number of a state measured versus the probing-light
detuning from the center of the Doppler-broadened $f=1\rightarrow F=2$
transition. The red points indicate the values calculated based on
experimental measurements (horizontal uncertainty comes from the uncertainty
of detuning, the evaluated vertical errors are small hence not visible), while
the blue line shows the theoretical dependence calculated from the absorption
measurements. Green dashed line indicates the smallest
$\kappa(\mathbb{C})=2.25$ achievable using this approach. (b) Relative
uncertainty of the linear inversion [see Eq. (9)] as a function of the
condition number $\kappa(\mathbb{C})$. Red points correspond to the
reconstruction of the state pumped with circularly polarized light propagated
along the $\mathbf{x}$-axis [see Fig. 2(a)] and blue points correspond to the
reconstruction of the state generated with linearly $\mathbf{z}$-polarized
pump light, propagating along the $\mathbf{x}$-axis [see Fig. 2(b)]. The
uncertainty of the condition number is significantly smaller than the size of
the data points and solid lines are added for clarity.

To further illustrate the effect of the probing-light detuning on the
reconstruction uncerainty and, hence, demonstrate the potential of this
approach, we perform a series of reconstructions of a state generated under
the same conditions but reconstructed using different probing-light detunings.
In our experiment, the detuning is changed from 50 to 270 MHz. The results of
these investigations are shown in Fig. 3(b). They demonstrate that the
reliability of the reconstruction deteriorates with the detuning and it
achieves the minimum in the vicinity of the center of the Doppler-broadened
$f=1\rightarrow F=2$ transition. This agrees with our theoretical prediction
of the condition-number detuning dependence, which we calculate assuming that
$\zeta(\Delta)=V_{I}(\Delta)/V_{R}(\Delta)$.

####  5.2.1 Conditional number versus the number of measurements

The repetition of specific measurements offers a straightforward and versatile
method for optimizing the relative weights of the observables used in the
state-reconstruction procedure. This approach allows for achieving an
arbitrarily small condition number, making it particularly valuable when the
previous method is infeasible or when the condition number is desired to be
smaller than the detuning-optimized bound (e.g., 2.25). However, it should be
noted that this technique is associated with a potential drawback; the number
of repetitions required to attain $\kappa(\mathbb{C})=1$ is typically
substantial, especially when dealing with initially high condition numbers, as
illustrated in Fig. 4.

Figure 4: Minimum achievable condition number versus a number of repetitions
of specific measurements. Three sets of points correspond to the condition
number at different detunings: 300 MHz (green triangles), 200 MHz (red
squares), and 60 MHz (blue diamonds). Solid line used for eye guidance.

##  6 Conclusions

In this study, we presented the first experimental implementation of a
quantum-state tomography technique, which was originally proposed in Ref.
[15]. The technique enabled us the successful reconstruction of collective
quantum states of a qutrit in room-temperature rubidium vapor at the $f=1$
ground state with a fidelity of 0.99. To overcome experimental challenges of
the reconstruction, we adapted the CYCLOPS technique, which allowed us to
achieve reliable reconstruction by mitigating a problem of unknown phase
delays present in measured signals. Additionally, we presented a comprehensive
analysis of the technique by introducing the conditional number, which
quantifies the reliability of the reconstruction. This parameter was
investigated versus different experimental factors, including tuning of the
probing light used for the reconstruction. We demonstrated that by appropriate
tuning of the light, the conditional numbers as low as 2.25 can be achieved
(where the conditional number of 1 refers to ideal reconstruction). We also
demonstrated that further improvement of the reconstruction (lowering the
conditional number) can be achieved by the repetition of the specific
measurements.

The successful implementation of the presented QST technique opens up avenues
for measuring a range of fundamental properties of qutrits. In future, we plan
to focus on exploring different measures of nonclassicality and establishing
their ordering for various classes of quantum states. We also plan on a
further development of the technique to demonstrate quantum-process
tomography, expanding the method capabilities in the characterization of
quantum operations and transformations. Finally, the ability to accurately
reconstruct the quantum states of atomic ensembles allows for experimental
optimization of generation of metrologically appealing quantum states. This is
the research direction that we currently pursuit in our work.

##  7 Acknowledgements

The authors would like to thank Arash D. Fard for his help in experimental
measurements. The work was supported by the National Science Centre, Poland
within the SONATA BIS programme (Grant No. 2019/34/E/ST2/00440). MK would like
to acknowledge support from the Excellence Initiative – Research University of
the Jagiellonian University in Kraków. A.M. is supported by the Polish
National Science Centre (NCN) under the Maestro Grant No.
DEC-2019/34/A/ST2/00081.

## References

  * [1] K. Mouloudakis, G. Vasilakis, V. G. Lucivero, J. Kong, I. K. Kominis, and M. W. Mitchell, “Effects of spin-exchange collisions on the fluctuation spectra of hot alkali-metal vapors,” Phys. Rev. A 106, 023112 (2022).
  * [2] R. Shaham, O. Katz, and O. Firstenberg, “Strong coupling of alkali-metal spins to noble-gas spins with an hour-long coherence time,” Nature Physics 18, 506–510 (2022).
  * [3] W. Happer, “Optical pumping,” Rev. Mod. Phys. 44, 169–249 (1972).
  * [4] O. Schmidt, R. Wynands, Z. Hussein, and D. Meschede, “Steep dispersion and group velocity below c/3000 in coherent population trapping,” Physical Review A 53, R27 (1996).
  * [5] K. Hammerer, K. Mølmer, E. S. Polzik, and J. I. Cirac, “Light-matter quantum interface,” Phys. Rev. A 70, 044304 (2004).
  * [6] H. Krauter, C. A. Muschik, K. Jensen, W. Wasilewski, J. M. Petersen, J. I. Cirac, and E. S. Polzik, “Entanglement generated by dissipation and steady state entanglement of two macroscopic objects,” Phys. Rev. Lett. 107, 080503 (2011).
  * [7] O. Katz, R. Shaham, E. S. Polzik, and O. Firstenberg, “Long-lived entanglement generation of nuclear spins using coherent light,” Phys. Rev. Lett. 124, 043602 (2020).
  * [8] K. Jensen, W. Wasilewski, H. Krauter, T. Fernholz, B. M. Nielsen, M. Owari, M. B. Plenio, A. Serafini, M. M. Wolf, and E. S. Polzik, “Quantum memory for entangled continuous-variable states,” Nature Physics 7, 13–16 (2011).
  * [9] M. Dąbrowski, M. Parniak, and W. Wasilewski, “Einstein–Podolsky–Rosen paradox in a hybrid bipartite system,” Optica 4, 272–275 (2017).
  * [10] E. E. Mikhailov and I. Novikova, “Low-frequency vacuum squeezing via polarization self-rotation in rb vapor,” Opt. Lett. 33, 1213–1215 (2008).
  * [11] I. H. Agha, G. Messin, and P. Grangier, “Generation of pulsed and continuous-wave squeezed light with 87Rb vapor,” Opt. Express 18, 4198–4205 (2010).
  * [12] W. Wasilewski, T. Fernholz, K. Jensen, L. S. Madsen, H. Krauter, C. Muschik, and E. S. Polzik, “Generation of two-mode squeezed and entangled light in a single temporal and spatial mode,” Opt. Express 17, 14444–14457 (2009).
  * [13] A. Heifetz, A. Agarwal, G. C. Cardoso, V. Gopal, P. Kumar, and M. Shahriar, “Super efficient absorption filter for quantum memory using atomic ensembles in a vapor,” Optics Communications 232, 289–293 (2004).
  * [14] M. Dąbrowski, R. Chrapkiewicz, and W. Wasilewski, “Magnetically tuned, robust and efficient filtering system for spatially multimode quantum memory in warm atomic vapors,” Journal of Modern Optics 63, 2029–2038 (2016).
  * [15] M. Kopciuch and S. Pustelny, “Optical reconstruction of the collective density matrix of a qutrit,” Phys. Rev. A 106, 022406 (2022).
  * [16] Y. I. Bogdanov, G. Brida, M. Genovese, S. P. Kulik, E. V. Moreva, and A. P. Shurupov, “Statistical estimation of the efficiency of quantum state tomography protocols,” Phys. Rev. Lett. 105, 010404 (2010).
  * [17] A. Miranowicz, K. Bartkiewicz, J. Peřina Jr., M. Koashi, N. Imoto, and F. Nori, “Optimal two-qubit tomography based on local and global measurements: Maximal robustness against errors as described by condition numbers,” Phys. Rev. A 90 (2014).
  * [18] K. Bartkiewicz, A. Černoch, K. Lemr, and A. Miranowicz, “Priority choice experimental two-qubit tomography: Measuring one by one all elements of density matrices,” Sci. Rep. 6, 19610 (2016).
  * [19] S. S. Roy and T. Mahesh, “Density matrix tomography of singlet states,” Journal of Magnetic Resonance 206, 127–133 (2010).
  * [20] A. Miranowicz, S. K. Özdemir, J. Bajer, G. Yusa, N. Imoto, Y. Hirayama, and F. Nori, “Quantum state tomography of large nuclear spins in a semiconductor quantum well: Optimal robustness against errors as quantified by condition numbers,” Phys. Rev. B 92, 075312 (2015).
  * [21] I. H. Deutsch and P. S. Jessen, “Quantum control and measurement of atomic spins in polarization spectroscopy,” Opt. Commun. 283, 681–694 (2010).
  * [22] C. Lovecchio, S. Cherukattil, B. Cilenti, I. Herrera, F. S. Cataliotti, S. Montangero, T. Calarco, and F. Caruso, “Quantum state reconstruction on atom-chips,” New Journal of Physics 17, 093024 (2015).
  * [23] K. Hammerer, A. S. Sørensen, and E. S. Polzik, “Quantum interface between light and atomic ensembles,” Rev. Mod. Phys. 82, 1041–1093 (2010).
  * [24] Y. Takahashi, K. Honda, N. Tanaka, K. Toyoda, K. Ishikawa, and T. Yabuzaki, “Quantum nondemolition measurement of spin via the paramagnetic faraday rotation,” Phys. Rev. A 60, 4974–4979 (1999).
  * [25] G. Colangelo, R. J. Sewell, N. Behbood, F. M. Ciurana, G. Triginer, and M. W. Mitchell, “Quantum atom–light interfaces in the gaussian description for spin-1 systems,” New Journal of Physics 15, 103007 (2013).
  * [26] D. d’Alessandro, _Introduction to quantum control and dynamics_ (CRC Press, 2021).
  * [27] S. Pustelny, V. Schultze, T. Scholtes, and D. Budker, “Dichroic atomic vapor laser lock with multi-gigahertz stabilization range,” Review of Scientific Instruments 87 (2016). 063107\.
  * [28] R. Freeman, _A Handbook of Nuclear Magnetic Resonance_ (Longman Scientific & Technical, 1987).
  * [29] F. A. Bonk, E. R. deAzevedo, R. S. Sarthour, J. D. Bulnes, J. C. C. Freitas, A. P. Guimarães, O. I. S., and T. J. Bonagamba, “Quantum logical operations for spin 3/2 quadrupolar nuclei monitored by quantum state tomography,” Journal of Magnetic Resonance 175, 226–234 (2005).
  * [30] M. Paris and J. Řehaček, _Quantum state estimation_ , vol. 649 (Springer Science & Business Media, 2004).
  * [31] K. Atkinson, _An Introduction to Numerical Analysis_ (Wiley, 1989).
  * [32] N. Higham, _Accuracy and Stability of Numerical Algorithms_ (Society for Industrial and Applied Mathematics, 1996).
  * [33] G. Golub and C. Van Loan, _Matrix Computations_ (Johns Hopkins University Press, 1996).

See pages - of supplementary_information.pdf


# Reactivity of HCO with CH3 and NH2 on Water Ice Surfaces. A Comprehensive
Accurate Quantum Chemistry Study

Joan Enrique-Romero  [  Albert Rimola  [<EMAIL_ADDRESS>Cecilia
Ceccarelli  [<EMAIL_ADDRESS>Piero Ugliengo  [
Nadia Balucani  [  Dimitrios Skouteris  [

###### Abstract

Interstellar complex organic molecules (iCOMs) can be loosely defined as
chemical compounds with at least six atoms in which at least one is carbon.
The observations of iCOMs in star-forming regions have shown that they contain
an important fraction of carbon in a molecular form, which can be used to
synthesize more complex, even biotic molecules. Hence, iCOMs are major actors
in the increasing molecular complexity in space and they might have played a
role in the origin of terrestrial life. Understanding how iCOMs are formed is
relevant for predicting the ultimate organic chemistry reached in the
interstellar medium. One possibility is that they are synthesized on the
interstellar grain icy surfaces, via recombination of previously formed
radicals. The present work focuses on the reactivity of HCO with CH3/NH2 on
the grain icy surfaces, investigated by means of quantum chemical simulations.
The goal is to carry out a systematic study using different computational
approaches and models for the icy surfaces. Specifically, DFT computations
have been bench-marked with CASPT2 and CCSD(T) methods, and the ice mantles
have been mimicked with cluster models of 1, 2, 18 and 33 H2O molecules, in
which different reaction sites have been considered. Our results indicate that
the HCO + CH3/NH2 reactions, if they actually occur, have two major
competitive channels: the formation of iCOMs CH3CHO/NH2CHO, or the formation
of CO + CH4/NH3. These two channels are either barrierless or present
relatively low ($\leq 10$ kJ/mol equal to about 1200 K) energy barriers.
Finally, we briefly discuss the astrophysical implications of these findings.

###### keywords:

Interstellar medium, Astrochemistry, DFT, iCOMs, grains

Université Grenoble Alpes] Univ. Grenoble Alpes, Institut de Planétologie et
Astrophysique de Grenoble (IPAG), 38000 Grenoble, France [Universitat Autònoma
de Barcelona] Departament de Química, Universitat Autònoma de Barcelona,
08193, Catalunya, Spain Universitat Autònoma de Barcelona] Departament de
Química, Universitat Autònoma de Barcelona, 08193, Catalunya, Spain Université
Grenoble Alpes] Univ. Grenoble Alpes, Institut de Planétologie et
Astrophysique de Grenoble (IPAG), 38000 Grenoble, France Università degli
Studi di Torino]Dipartamento di Chimica and Nanostructured Interfaces and
Surfaces (NIS), Università degli Studi di Torino, Via P. Giuria 7, 10125
Torino, Italy Università degli Studi di Perugia]Dipartamento di Chimica,
Biologia e Biotecnologie, Università degli Studi di Perugia, Via Elce di Sotto
8, I-06123 Perugia, Italy [INAF]Osservatorio Astrofisico di Arcetri, Largo E.
Fermi 5, 50125 Firenze, Italy [UniversitèGrenoble Alpes]Univ. Grenoble Alpes,
Institut de Planètologie et Astrophysique de Grenoble (IPAG), 38000 Grenoble,
France Master-Up, Perugia, Italy]Master-Up, Strada Vicinale Sperandio 15,
I-06123 Perugia Italia Perugia, Italy ISM,ASW

##  1 1\. Introduction

It has been long demonstrated that star forming regions are places with a rich
organic chemistry (e.g. 1, 2, 3, 4, 5, 6). Although there are no proofs that
organic molecules formed in the interstellar medium (ISM) did play a role in
the emergence of terrestrial life, there is mounting evidence that they were
inherited by the small objects of the Solar System: for example, carbonaceous
chondrites and comets contain a wide variety of organic molecules, some of
them probably being a direct heritage of the ISM based on their relative
abundances and ratio of deuterated versus hydrogenated species (e.g. 7, 4, 8,
9).

Knowing how detected interstellar complex organic molecules (iCOMs: C-bearing
molecules with at least six atoms 10) are formed and destroyed is not only
important per se but also because it is the only way to understand the
ultimate organic complexity present in the ISM. Indeed, there is an intrinsic
limit to the detection of large iCOMs (excluding linear chains), which is
caused by the fact that the larger the molecule the larger the number of
rotational transitions (because of the larger number of functional groups of
the iCOM) and, consequently, the weaker the intensity of the lines. As a
result, the numerous and weak lines of large iCOMs produce a “grass” of lines
in the spectra, which makes the identification of a large molecule eventually
impossible. Therefore, there is a limit to the largest detectable iCOM, and
this has a direct consequence: we need to rely on our knowledge of the
processes to predict which large molecules are synthesized in the ISM.

How iCOMs are formed is a question that has baffled astrochemists for decades.
In the 90s it was thought that gas-phase reactions were the dominant formation
processes (e.g. 11, 12). However, subsequent astronomical observations (e.g.
13, 2), laboratory experiments (e.g. 14), and theoretical calculations (e.g.
15) challenged this synthetic route model. As a result, a new paradigm was
proposed, which assumes that most (if not all) iCOMs are formed on the
surfaces of icy interstellar grains.16, 17 According to this paradigm,
radicals are created by the UV photons impinging on the ice grain mantles,
which are formed during the cold cloud/prestellar phase and contain simple
hydrogenated species (mostly water but also other species, such as methanol
and ammonia, in smaller quantities). With the warming up of the environment
caused by the birth of a protostar, the radicals trapped in the ice can
diffuse on the grain surfaces and react between them to form iCOMs. The
reaction process is assumed to be a direct combination of radicals, which are
considered as “lego-like” blocks.

This paradigm is usually assumed in current gas-grain astrochemical models
(e.g. 18, 19, 20, 21, 22). However, there are still many uncertainties that
makes the paradigm not fully validated. First, it is not clear that the lego-
like radical combination is actually an effective process (refs: e.g. 23, 24).
Second, it appears that the role of gas-phase reactions has been
underestimated in the past. 25, 26, 27, 28, 29, 9, 30, 31, 32 Finally, models
based on the “exclusive grain-surface”paradigm are unable to reproduce the
observed abundance of several iCOMs (e.g. 33, 34, 35).

The numerous laboratory experiments reported in the literature are extremely
useful to study the possible processes occurring on ices illuminated by UV
photons and/or irradiated by energetic particles (e.g. 36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46, 47). However, investigating surface-induced iCOMs
formation by means of experimental techniques is extremely difficult, if not
impossible, since the actual interstellar conditions cannot exactly be
reproduced in terrestrial laboratories. For example, the UV irradiation used
in laboratory experiments is, for practical reasons, more than million times
larger than the one impinging on the interstellar grain mantles: this causes
an instantaneous injection of energy of several orders of magnitude larger
than that absorbed by the interstellar grain mantles and, consequently, likely
introduces different reaction routes (e.g. energy barriers are probably
overcome in laboratory experiments while they are not in the cold ISM).
Similarly, the H-flux in experiments is also extremely larger compared with
the actual ISM conditions, causing a dramatic amplification of the results
towards a full saturation of the compounds. A suitable alternative, and
certainly a complementary method for such investigations, is the use of
theoretical simulations based on quantum chemical calculations. Indeed, these
calculations provide a description of the surface reactions from an atomistic
perspective, providing unique, relevant energetic information, and,
accordingly, they can be useful to assess the validity of the “exclusive
grain-surface” paradigm.

In this work, we focus on two reactions occurring on amorphous solid water
(ASW) surfaces: HCO + CH3 and HCO + NH2. In the “exclusive grain-surface”
paradigm mentioned above, it is expected that the radical coupling produces
two iCOMs: CH3CHO (acetaldehyde) and NH2CHO (formamide). However, previous
works by some of us 24, 23 showed that other reactive channels can compete
with the iCOM formation. Specifically, the two reactions can lead to the
formation of CO + CH4 and CO + NH3, respectively, in which the H atom of HCO
is transferred to the radical partner. Similar processes were identified
computationally when HCO reacts with CH3 on surfaces of CO-pure ices.48

The goal of the present work is to carry out a systematic study of two
reactions, i.e., HCO + CH3 and HCO + NH2, considered here as prototype
reactions for the formation of iCOMs, using different approximations for the
calculations and models of ASW with the aim to: (1) understand how the
different methods and models affect the results, (2) individuate the most
convenient methods and models to use in future calculations of other similar
radical-radical systems, and (3) identify are the products of the reactions
for different conditions. To this end, the present work focuses on the
following three points:

  1. 1.

Methodology benchmark: (i) The energy barriers for reactions between the two
couples of radicals in the presence of 1 and 2 water molecules are computed
with two DFT methods (B3LYP and BHLYP) and compared to the values calculated
with the multi-reference CASPT2 method; (ii) The interaction of the three
radicals (CH3, HCO and NH2) with 1 and 2 water molecules is studied and bench-
marked taking as reference the binding energies computed at the CCSD(T) level.

  2. 2.

Radical-surface binding enthalpies: We study the binding of the three radicals
(CH3, HCO and NH2) to an ASW cluster model of 18 water molecules and to a
larger cluster of 33 water molecules sporting two different morphological
sites, a cavity and its side.

  3. 3.

Radical-radical reactivity: The reactivity of the two radical couples is
studied (i) on the 18 water molecules cluster, and (ii) on the two different
morphological sites of the 33 water molecules cluster.

It is worth mentioning that the systems to deal with here sport an additional
complexity from an electronic structure point of view. That is, two radicals
interacting with water ice, in which the unpaired electrons have opposite
signs, constitute a singlet biradical system. Describing this electronic
situation with quantum chemistry calculations is delicate. Ab initio multi-
reference methods are a good choice to describe biradical systems but they are
extremely expensive and, accordingly, unpractical for large systems.
Alternatively, a good compromise between accuracy and computational cost is
the DFT broken-symmetry approach.49. By using this method, however, our own
experience24, 23 indicates that one has to be sure that the initial guess wave
function corresponds to the actual singlet biradical state (i.e., the unpaired
electrons being localized on the corresponding radicals and with opposite spin
signs), as it may well happen that the initial wave function represents
wrongly a closed-shell-like situation, in which the unpaired electrons are 50%
distributed among the two radicals. Results derived from one or the other
situation are dramatically different.

This article is organized as follows. First the adopted methods are presented
(§ 2), then the results are provided following the 3 objectives described
above (§ 3) and finally a discussion, including the astrophysical
implications, and the conclusions are presented in § 4 and 5, respectively.

##  2 2\. Methods

All DFT and CCSD(T) calculations were performed using the GAUSSIAN0950
software package, while the multi-reference calculations were carried out with
the OpenMolcas 18.0951, 52, 53, 54 program.

Stationary points of the potential energy surfaces (PESs) were fully optimized
using two hybrid density functional theory (DFT) methods: B3LYP and BHLYP.
These methods have the same Lee, Yang, Parr correlation functional (LYP)55 but
differ on the exchange functionals: the Becke’s three parameter (B3), which
includes a 20% of exact exchange in its definition,56 and the Becke’s half-
and-half (BH), which mixes the pure DFT and the exact exchange energy in a 1:1
ratio.57. For B3LYP calculations, the Grimme’s D3 dispersion term58 was
accounted for during the geometry optimizations. In contrast, for BHLYP, both
the D259 and D3 dispersion terms were included in a posteriori way onto the
pure BHLYP optimized geometries.

A calibration study was first carried out for (i) the NH2 \+ HCO and CH3 \+
HCO reactivities in the presence of 1 and 2 water molecules (W1 and W2,
respectively) and (ii) the interaction of each radical with W1 and W2. For
this calibration study, the DFT methods were combined with the Pople’s
6-311++G(2df,2pd) basis set. As reference for the reactivity results we used
single point energy calculations at the multi-reference CASPT2 level combined
with the Dunnning’s cc-pVTZ basis using as initial guess the orbitals
generated at CASSCF(2,2) level. In the same way, single point energy
calculations at the CCSD(T)/aug-cc-pVTZ level were also carried out in order
to compare them to CASPT2 values. Regarding the interaction energies, the same
DFT methods were compared to CCSD(T)/aug-cc-pVTZ level. All the single point
energy calculations for this benchmark study were carried out on the B3LYP-D3
optimized geometries.

Radical-radical reactivity was also studied more realistically on two
amorphous solid water (ASWs) ices modelled by molecular clusters consisting of
18 (W18) and 33 (W33) water molecules, which were also used in previous
works.60, 24, 23 The optimized structures are shown in Figure 1.
Interestingly, W33 exhibits a hemispherical cavity and, accordingly, we
studied the surface processes considering both this cavity and an extended
side of the ice surface (see Figure 1(b)), as they exhibit different surface
properties. In order to make the calculations computationally affordable, for
these cases the DFT methods were combined with the Pople’s 6-31+G(d,p) basis
set.

Figure 1: Structures of the 18 and 33 water molecules clusters, (a) and (b)
respectively, optimised at the BHLYP/6-31+G(d,p) level.

All stationary points were characterized by the analytical calculation of the
harmonic frequencies as minima (reactants, products and intermediates) and
saddle points (transition states). Intrinsic reaction coordinate (IRC)
calculations at the same level of theory were carried out when needed to
ensure that the transition states connect with the corresponding minima.
Thermochemical corrections to the potential energy values were carried out
using the standard rigid rotor/harmonic oscillator formulas61 to compute the
zero point energy (ZPE) corrections.

Adsorption energies of HCO, CH3 and NH2 on W18 and W33 were refined with
single point energy calculations combining the DFT methods with the extended
6-311++G(2df,2pd) basis set and corrected for the basis set superposition
error (BSSE). Considering A as the adsorbate and B as the surface cluster
model, the BSSE-non-corrected adsorption energy was computed as $\Delta
E_{ads}=E^{AB}_{AB}(AB)-E^{A}_{A}(A)-E^{B}_{B}(B)$, where superscripts denote
the basis set used and the subscripts the geometry at which the calculation
was done. BSSE-corrected energies were calculated as:

$\Delta E_{ads}^{CP}(AB)=\Delta
E_{ads}+{BSSE}(A)+{BSSE}(B)+\delta^{A}(A)+\delta^{B}(B)$ (1)

where the BSSE values were calculated following the Boys and Bernardi
counterpoise correction method ($BSSE(A)=E^{AB}_{AB}(A)-E^{A}_{AB}(A)$, 62),
and where the deformation of each monomer was also accounted for
($\delta^{A}(A)=E^{A}_{AB}(A)-E^{A}_{A}(A)$).

Inclusion of ZPE corrections allowed us to obtain adsorption enthalpies at 0
K:

$\Delta H_{ads}(AB)=\Delta E_{ads}^{CP}(AB)+\Delta ZPE$ (2)

In the sign convention followed in this work, the adsorption energy is the
negative of the binding energy, i.e. $\Delta E_{ads}=-\Delta E_{bind}$.

The systems containing two radical species were first optimized in the triplet
electronic state, which was then optimized in the singlet state to describe
the biradical system. Structures with doublet and triplet electronic states
were simulated with open-shell calculations based on the unrestricted
formalism. Singlet biradicals systems were calculated adopting the
unrestricted broken symmetry (BS) approach, in which the most stable initial
wave function was found using the stable=opt keyword in Gaussian09.

Following the International System Units, all energy units are given in kJ
mol-1, whose conversion factor to K is 1 kJ mol${}^{-1}=$ 120.274 K.

##  3 3\. Results

In this section, results of the calibration study devoted to the reactivity
and interaction of HCO with CH3 and NH2 in the presence of 1 and 2 water
molecules are first presented. Then, results on adsorption properties and the
radical-radical reactions on W18 and W33 are reported.

###  3.1 3.1. Methodology benchmark

This section aims to be a calibration study to check the reliability of the
B3LYP and BHLYP methods for (i) the radical-water interactions and (ii) the
activation energy of the reactions of CH3/HCO and NH2/HCO, in both cases with
one and two water molecules (W1 and W2 hereafter) as minimal models
representing an ice surface. The references used for the study of the binding
and activation energies are CCSD(T) and CASPT2, respectively. In all cases,
single point energies were computed onto the B3LYP-D3 optimised geometries. In
the supporting information (SI) section the reader can find these geometries
as well as the relative errors of the data presented in the following.

Table 1 contains the binding energies for the three radicals involved in this
study (namely CH3, NH2 and HCO) interacting with W1 and W2. The systems where
a radical interacts with W1 have been based on those presented by Wakelam et
al. 63, where the initial structures were built in a chemical-wise manner
following the ability of each component of the radical to establish a hydrogen
bond to a single water molecule, e.g. in the cases where NH2 and HCO interact
with W1 two possibilities were considered: the radicals acting as either
H-bond donors (through one of their H atoms), or acting as H-bond acceptors
(through the N and O atoms, respectively). The initial geometries of the
systems with W2 were built similarly, with the radicals having the maximum
number of hydrogen bonds to the two water molecules.

Table 1: ZPE- and BSSE-non-corrected binding energies computed with different methods for the interaction of CH3, HCO and NH2 with 1 and 2 water molecules (W1 and W2) at the B3LYP-D3 and BHLYP levels. None, D2 (accounting for 2-body interactions) and D3 (accounting for 2- and 3-body interactions) dispersion corrections were considered for the latter. BHLYP-based and CCSD(T) values were calculated as single points on the B3LYP-D3 optimised geometries. Energy units are in kJ/mol. | B3LYP-D3 | BHLYP | BHLYP-D2 | BHLYP-D3 | CCSD(T)
---|---|---|---|---|---
CH3/W1 | 9.5 | 5.2 | 11.3 | 8.1 | 6.9
HCO/W1 (H) | 12.9 | 10.9 | 13.6 | 13.1 | 12.0
HCO/W1 (O) | 11.7 | 9.7 | 12.5 | 12.1 | 11.5
NH2/W1 (H) | 14.1 | 12.0 | 14.9 | 14.4 | 12.8
NH2/W1 (N) | 23.4 | 21.5 | 24.3 | 23.6 | 21.6
CH3/W2 | 8.6 | 3.1 | 11.3 | 8.0 | 7.1
HCO/W2 | 14.8 | 9.2 | 16.1 | 15.1 | 13.2
NH2/W2 | 41.6 | 36.6 | 44.5 | 41.9 | 38.3

It can be seen that CH3 is the species having the weakest interaction with the
water molecules (6.9-7.1 kJ/mol at the CCSD(T) level). HCO and NH2, instead,
can form intermediate and strong H-bonds respectively, and, accordingly, they
show higher binding energies (11.5-13.2 kJ/mol for the former and 12.8-38.3
kJ/mol for the latter at the CCSD(T) level). This trend is in agreement with
that found by Wakelam et al. 63, in which the binding energies of these
species with W1 were 13, 19-23 and 23-38 kJ/mol for CH3, HCO and NH2
respectively, computed at the M062X/aug-cc-PVTZ level. Note, however, that the
quoted values are not the final values used in the model by Wakelam et al. 63,
as they also use binding energies from other sources in the literature in
order to fit experimental data.
Regarding the performance of the methods, the best ones are B3LYP-D3 and
BHLYP-D3, with a relative error of 2-39% and 5-18% respectively (see SI
section). The pure BHLYP method systematically underestimates the binding
energies, providing strong deviations for CH3 and HCO, specially on W2.
BHLYP-D2 dramatically overestimates the binding energies of CH3-containing
systems (relative errors of $\sim$ 60%), indicating that the D2 term probably
accounts for dispersion in excess for this kind of weakly bound complexes.

Table 2 shows the calculated energy barriers of CH3/HCO and NH2/HCO on W1 and
W2. Three different possible reactions have been identified: i) radical-
radical coupling, leading to formation of the iCOMs (i.e. CH3CHO and NH2CHO);
ii) direct hydrogen abstraction, in which the H atom of HCO is transferred to
the other radicals, forming CO + CH4 and CO + NH3, respectively; and iii)
water assisted hydrogen transfer, the same as ii) but the H transfer is
assisted by the water molecules adopting a H relay mechanism. These reactions
will be referred to as Rc, dHa and wHt, respectively, along the work.

Table 2: ZPE-non-corrected energy barriers ($\Delta E^{\ddagger}$), computed with different methods, for radical-radical coupling (Rc), direct hydrogen abstraction (dHa) and water assisted hydrogen transfer (wHt) reactions of NH2 \+ HCO and CH3 \+ HCO in the presence of 1 and 2 water molecules computed at the B3LYP-D3 and BHLYP-based levels. For the latter two dispersion corrections have been considered: D2 (which considers all 2-body interactions) and D3 (which considers 2- and 3-body interactions). BHLYP-based, CCSD(T) and CASPT2 values were calculated as single point energy calculations on the B3LYP-D3 optimised geometries. Energy units are in kJ mol-1. NB stands for “No Barrier” and means that the process is found to be barrierless. System  | Process |  B3LYP -D3 | BHLYP |  BHLYP -D2 |  BHLYP -D3 | CCSD(T) | CASPT2
---|---|---|---|---|---|---|---
NH2/HCO $\cdots$ W1 | Rc | NB | NB | NB | NB | NB | NB
dHa | 3.5 | 6.1 | 4.1 | 5.1 | 5.3 | 4.2
wHt | 10.9 | 48.6 | 46.2 | 48.8 | 43.8 | 48.3
CH3/HCO $\cdots$ W1 | Rc | NB | NB | NB | NB | NB | NB
dHa | 3.1 | 6.9 | 3.3 | 5.1 | 5.0 | 1.5
wHt | 23.5 | 59.4 | 55.5 | 58.6 | 60.0 | 52.2
NH2/HCO $\cdots$ W2 | Rc | 6.5 | 8.5 | 6.1 | 6.8 | 6.3 | 6.1
dHa | NB | NB | NB | NB | NB | NB
wHt | 15.5 | 65.8 | 62.5 | 65.1 | 52.3 | 52.9
CH3/HCO $\cdots$ W2 | Rc | NB | NB | NB | NB | NB | NB
dHa | 1.6 | 4.0 | 1.5 | 1.6 | 1.6 | 1.0
wHt | 30.7 | 65.8 | 62.5 | 65.1 | 64.0 | 52.0

According to these values, two general trends are observed: (i) Rc and dHa are
either barrierless (meaning that the initial biradical structures were not
stable) or have relatively low energy barriers (1.0-6.1 kJ/mol for CASPT2),
and (ii) wHt are the processes presenting the highest energy barriers (around
50 kJ/mol for CASPT2).

For those Rc and dHa reactions having an energy barrier, one can see that the
worst performance is given by the pure BHLYP method compared to CASPT2. The
rest of the methods perform similarly, BHLYP-D2 and B3LYP-D3 being the best
ones, which are followed by BHLYP-D3. Regarding the wHt reactions, B3LYP-D3
dramatically underestimates the energy barriers (errors of $\sim$ 40-80%),
while BHLYP-based methods perform reasonably well with errors below $\sim$
25%. CCSD(T) performs similarly to the BHLYP-based methods, presenting
moderate energy barrier deviations compared to the CASPT2 values. All errors
can be found in the SI section.

In summary, BHLYP-D3 provides the most reliable binding energies, B3LYP-D3,
BHLYP-D2 and BHLYP-D3 show good performances for energy barriers related to Rc
and dHa processes, while BHLYP-D2 and BHLYP-D3 perform good for wHt energy
barriers. Accordingly, and with the aim to be consistent, in the following
sections (devoted to the binding and potential energy profiles on the W18 and
W33 cluster models) we provide the results at the BHLYP-D3 level of theory,
while results based on BHLYP-D2 (on both W18 and W33) and on B3LYP-D3 (on W18)
are reported in the SI.

###  3.2 3.2. Radical-surface binding enthalpies

Here we present the ZPE- and BSSE-corrected binding energies of the radicals
on the W18 and W33 cluster models and their comparison with recent literature.
These values are important because binding energies are essential parameters
in astrochemical modelling studies.

As it was stated in the previous section, the binding of the radicals on the
water ice surfaces are mainly dictated by hydrogen bonds (H-bonds) and
dispersion interactions. The clusters exhibit several potential binding sites
(Figure 1) and, accordingly, different radical/surface complexes can exist.
Here, for the sake of simplicity, we choose those complexes in which the
inter-molecular interactions between the two partners are maximized. The
underlying assumption is that radicals on the water ice surfaces have had
enough time to thermalize and establish the largest number of inter-molecular
interactions with the surface, as very probably is the case in the
interstellar conditions. Figure 2 shows the BHLYP-D3 fully optimised complexes
on W18 and on the two differentiated W33 sites: its side (W33-side) and its
cavity (W33-cav). The corresponding BSSE-corrected binding enthalpies (at 0 K)
are also shown. The same information at BHLYP-D2 and B3LYP-D3 is available in
SI.

Figure 2: Fully optimized geometries for the binding of CH3, HCO and NH2 on
(a) W18, (b) W33-side and (c) W33-cav at BHLYP-D3/6-31+G(d,p). Bond distances
are in Å. Binding enthalpies (at 0 K) values (in kJ/mol) are corrected for
BSSE and are shown below their respective structures.

The binding enthalpy trend is NH2 $>$ HCO $>$ CH3, irrespective of the cluster
model and surface morphology, in agreement with the results of the previous
section. Specifically, the binding enthalpy ranges are: 22.4-34.4, 16.9–22.2
and 5.0–9.5 kJ/mol for NH2, HCO and CH3, respectively. These values compare
well with those previously computed by Enrique-Romero et al. 24 for HCO and
CH3 on the W18 cluster (19 and 6 kJ/mol) and by Rimola et al. 23 for NH2 and
HCO on the W33-cav cluster (33.5 and 17.5 kJ/mol). In the model presented in
Wakelam et al. 63, the authors reported binding energies of 27, 20 and 13
kJ/mol for NH2, HCO and CH3, respectively, in rough agreement with our
results. Sameera et al. 64 have recently computed the binding energies of HCO
and CH3 adsorbed on hexagonal ice slabs111Authors used either full DFT with a
double-$\zeta$ basis set or the ONIOM approach combining DFT with force
fields. They reported several values owing to the different adsorption modes
of these two species adopting the two computation approaches. This allowed
them to report the wide range of binding energies reported here.. They found
HCO binding energies ranging between 12–40 kJ/mol, and 11–25 kJ/mol for HCO
and CH3, respectively. These values are similar or moderately larger than the
values found in our W33-cav cluster.

Finally, it is worth mentioning that these complexes do not present hemi-
covalent bonds, as it was the case for CN in Rimola et al. 23. Attempts to
identify this type of binding in the current systems have been made but the
initial structures collapsed to the complexes presented here upon geometry
optimization. Therefore, for these systems, the interaction of the radicals
with the ice surfaces is essentially dictated by H-bonds and dispersion
forces.

###  3.3 3.3. Radical-radical reactivity

In this section, the reactivity of HCO with CH3 and NH2 on W18, W33-side and
W33-cav are investigated. Given that the chemical environment between W18 and
W33-side is similar (note that the surface morphology of W33-side is very
similar to that of W18 because W33 was built up by joining two W18 units23),
the comparison of the results between these two models will allow us to assess
the effects introduced by the ASW model size. On the other hand, comparison
between the reactivity on W33-side and W33-cav will allow us to assess the
effects due to different surface environments, namely the presence of a higher
number of radical/surface interactions. Please note that the cavity sites are
probably more representative of the interstellar conditions than the side
sites, as the vast majority of radicals are trapped in the bulk of the ice.
Figures 3, 4 and 5 show the PESs of the reactivity of CH3/HCO and NH2/HCO
calculated at BHLYP-D3 on top of W18, W33-side and W33-cav, respectively.
Based on the results for the reactivity in the presence of W1 and W2 (see
§3.1), the same three reaction paths, i.e., Rc, dHa and wHt, have been
investigated. However, as the later processes are those exhibiting the highest
activation enthalpies (as high as 100.6, 78.8, 92.5 and 79.1 kJ/mol for
CH3/HCO$\cdots$W18, CH3/HCO$\cdots$W33-cav, NH2/HCO$\cdots$W18 and
NH2/HCO$\cdots$W33-cav, respectively), for the sake of clarity and with the
aim of focusing only on the reactions that might play a role in interstellar
chemistry, all results related to wHt can be found in SI, this section only
showing the Rc and dHa reactions calculated at BHLYP-D3. In the same way, the
reader can also find in SI the results for all systems computed at BHLYP-D2
and, for W18, also at B3LYP-D3.

The initial structures of these systems were built according to the
interaction patterns present in the single adsorbed radical complexes (see
reactant structures of Figures 3-5). The Rc and dHa reactions, forming
acetaldehyde or formamide and CO+CH4 or CO+NH3, respectively, both take place
through a single step. That is, the bond formation between the two radicals
for Rc, and the H transfer for dHa. The only exception is the dHa reaction of
CH3/HCO$\cdots$W18, which displays first a submerged activation energy step
where HCO breaks its H-bonds with the surface to facilitate the H transfer
(see dHa-TS1 of Figure 3(a)). By comparing the Rc and dHa reactions on these
three cluster models, the general trend is that Rc activation enthalpies are
lower than the dHa ones, between 0.5 - 4 kJ/mol lower. This was already
observed in the presence of W1 and W2, where most of the Rc reactions are
barrierless while the dHa present physical (although low) energy barriers. It
should be noticed that the HCO/CH3 reactions on W33-side (Figure 4a) are
barrierless due to the ZPE correction, as opposed to the the W1 and W2 cases
where the lack of barrier is because of the unstable singlet evolving directly
towards products. Indeed, all the stationary points shown in Figure 4(a) have
been identified as stable structures in the pure PESs (i.e., without ZPE
corrections) but after introduction of ZPE corrections, the transition state
becomes lower in energy than the reactants, and hence the barrierless
character.

Another interesting point is that results for W18 and W33-side present some
differences, particularly for the HCO + CH3 reactions: on W33-side, both Rc
and dHa are barrierless while on W18 they present energy barriers. This is
indicative of the fact that the size of the cluster for this radical-radical
reactivity is important. Finally, no clear trends related to the effect of the
ASW ice morphology can be obtained by comparing the W33-side results with the
W33-cav. That is, for HCO+CH3 reactions, activation enthalpies are higher on
W33, while the opposite occurs for HCO+NH2 ones.
In the following section, a comprehensive discussion of all these results is
provided.

Figure 3: ZPE-corrected Rc (solid black line) and dHa (dashed green lines)
PESs for (a) HCO/CH3 $\cdots$W18 and (b) HCO/NH2 $\cdots$W18 calculated at the
BHLYP-D3/6-31+G(d,p). Energy units are in kJ/mol and distances in Å. Notice
that dHa-TS1 for (a) HCO/CH3 $\cdots$W18 lies below the intermediate dHa-I,
due to the ZPE correction, setting effectively a single energy barrier: dHa-
Ts2.

Figure 4: ZPE-corrected Rc (solid black lines) and dHA (green dashed lines)
PES for the HCO/CH${}_{3}\cdots$W33-side (a) and HCO/NH${}_{2}\cdots$W33-side
(b) systems optimized at the BHLYP-D3 theory level. Energy units are in kJ/mol
and distances in Å. Notice that both Rc and dHa TS for CH3/HCO lie below the
energy of reactants due to the ZPE correction.

Figure 5: ZPE-corrected Rc (solid black lines) and dHA (green dashed lines)
PES for the HCO/CH${}_{3}\cdots$W33-cav (a) and HCO/NH${}_{2}\cdots$W33-cav
(b) systems optimized at the BHLYP-D3 theory level. Energy units are in kJ/mol
and distances in Å. BHLYP-D2 values can be found in the SI.

##  4 4\. Discussion

###  4.1 4.1 Reaction channels and competition to iCOMs formation

Our computations show that if HCO + CH3 and HCO \+ NH2 react on top of an ASW
ice, they have two sets of possible reaction products: i) formation of iCOMs
(Rc process), where the two radical species meet and couple, and ii) the
formation of the hydrogenated CH4/NH3 species, where the H atom of HCO
migrates to CH3/NH2, which can happen either directly (dHa process) or through
the ice water molecules adopting a H transfer relay mechanism (wHt process).
The energetics associated with each process at the BHLYP-D3 level are
summarised in Table 3.

Table 3: Highest activation enthalpies (at 0 K) for Rc, dHa and wHt reactions on W18, W33-side and W33-cav at the BHLYP-D3 level. Note that NB stands for “No Barrier”. Units are in kJ/mol. | W18 | W33-side | W33-cav
---|---|---|---
Sys | Rc | dHa | wHt | Rc | dHa | wHt | Rc | dHa | wHt
HCO + CH3 | 2.6 | 5.2 | 100.6 | NB | NB | - | 6.8 | 10.2 | 78.8
HCO + NH2 | 3.1 | 5.0 | 92.5 | 5.1 | 5.4 | - | 2.1 | 1.6 | 79.1

The Rc and dHa reactions show, in all the studied systems, similar energetic
features, i.e., they are either barrierless or exhibit relatively low energy
barriers. The highest pair of energy barriers concerns the HCO + CH3 reactions
on W33-cav, i.e., 6.8 and 10.2 kJ/mol for Rc and dHa, respectively. As a
general trend, dHa reactions have slightly higher activation energies than Rc
(by as much as 3.4 kJ/mol for HCO/CH3 and 1.9 kJ/mol for HCO/NH2). In some
cases, like HCO/CH${}_{3}\cdots$W33-side, HCO/NH${}_{2}\cdots$W33-side and
HCO/NH${}_{2}\cdots$W33-cav both Rc and dHa can be considered as competitive
reactions given the small activation energy differences. In contrast, the
lowest energy barriers for wHt HCO + CH3 and HCO + NH2 reactions are $\sim$ 80
kJ/mol, respectively. These values are larger than any Rc and dHa energy
barrier and, accordingly, wHt reactions cannot be considered by any means as
competitive channels. The explanation of these energetic differences is
provided by the reaction mechanisms. Rc and dHa reactions take place, in most
of the cases, in a concerted way, in which the radicals, in essence, have to
partly break the interactions with the surface to proceed with the reaction.
In contrast, most of the wHt reactions adopt a multi-step mechanism since the
H transfer, which is assisted by different ice water molecules, involves
different breaking/formation bonds. In these cases, high energy intermediates
consisting of the coexistence of HCO and an OH radical are involved (see SI).

In previous works by some of us, (e.g. 24) the wHt reactions between CH3 \+
HCO on the W18 cluster model were observed to spontaneously occur during
geometry optimisation, i.e. they were found to be barrierless. The difference
with the computations presented in this work resides on the fact that, in the
previous work24, the initial wave function did not describe a singlet
biradical system but a metastable singlet closed-shell-like one, and hence the
spontaneous evolution to form CO+CH4. In this work, as well as in Rimola et
al. (2018) 23, the initial wave function is actually describing a singlet
biradical situation, which leads to a significant stabilization of the
reactants and hence the presence of high energy barriers.

Finally, some words related to the chemical role played by the ice on these
reactions are here provided. The Rc and dHa processes in the gas-phase
(namely, in the absence of the icy grain) are, in both cases, barrierless. In
contrast as explained above, in the presence of the surface, they exhibit,
although low, energy barriers. Accordingly, from a rigorous chemical kinetics
standpoint, the grains slow down the reactions. This leads us to think that a
major role played by the grains is as that of third bodies, by quickly
absorbing the nascent energy associated with the reactions, hence stabilising
the products. This aspect is particularly appealing in the iCOMs formation
processes via radical recombination since in the gas phase iCOMs can
redissociate back to reactants if they are not stabilised through three-body
reactions. As we will discuss below, the water morphology plays a role in the
reaction energetics, although it does not change the essence of our
conclusions.

###  4.2 4.2 Influence of the water ice surface model

Clear differences arise when comparing W33-side and W33-cav, as radicals on
the latter exhibit more inter-molecular interactions with the surface. This
can be seen for example in the binding energies (higher on W33-cav than on
W33-side, see Figure 2). But also on the Rc and dHa energy barriers, for which
a different behaviour is observed. Indeed, for HCO + CH3 these reactions are
barrierless on W33-side, while on W33-cav they present energy barriers of 7
(Rc) and 10 (dHa) kJ/mol. On the other hand, for HCO + NH2 the opposite
behaviour is observed: the energy barriers are higher on W33-side than on
W33-cav (see Table 3). This might be indicating that the different polarity of
the radicals, i.e., CH3/apolar and NH2/polar, is important when several polar-
based inter-molecular interactions surround the reaction sites, like it is the
case of W33-cav.

Finally, the size difference between the W18 and W33 models does not seem to
provide a consistent trend, neither for binding energies nor for the PESs.
This is probably due to the modest energetics of reactions of interest, which
are the result of many intermingled effects, i.e. H-bond and dispersion
interaction strength, small charge transfer and polarization. All these
components are affected by the nature and size of the water adopted clusters
without a definite and predictable structure-properties relationship.

###  4.3 4.3 Astrophysical implications

A major goal of this study is to understand whether iCOMs can be formed on the
icy grain surfaces by the direct combination of radicals, a process assumed to
be efficient in the majority of current astrochemical models (see
Introduction). The present computations show that (i) there is a feasible
channel leading to iCOM formation through radical-radical combination, (ii)
this channel may possess a barrier, and (iii) there is at least a competitive
reaction where radicals exchange a hydrogen atom, the outcome of which is
somewhat a step backwards in chemical complexity as the products are simple
hydrogenated species (CH4 and NH3 in the current work) and CO.

The present computational data does not allow us to definitively exclude the
presence of a barrier in the radicals combination. Indeed, although common
sense would indicate a lack of barrier, calculations show that the presence of
the ice water molecules introduces an inter-molecular interaction that depends
on where the radicals are placed and on the radical polarity. This interaction
probably necessitates energy to be broken: it is not obvious that this energy
is available in the ISM environments. In fact, according to Garrod et al.
model 16, once the radicals are formed, they remain frozen on the ice and
subsequently more ice layers build up on top. The radicals remain imprisoned
in cavity structures inside the mantle and once the temperature reaches $\sim$
30 K due to the evolution of the central protostar, they diffuse and react.
Among the 3 reactions sites discussed in this article (W18, W33-side and
W33-cav), the one best resembling this picture is given by W33-cav (see Figure
1) due to the larger number of inter-molecular interactions. The Rc and dHa
energy barriers in the cavity have been shown to be larger for HCO + CH3 than
for HCO + NH2, due probably to the different polarity of CH3 and NH2. If
converted to Kelvin, the Rc and dHa barriers for HCO + CH3 and HCO + NH2 on
W33-cav are about 800, 1200, 250 and 190 K, respectively (see Table 3). Thus,
the efficiency of these reactions is not expected to be very high, especially
for HCO + CH3.

It has also to be noted that the starting points from where we study the
reactivity in W33-cav (see reactant geometries from Figure 5) contain both
radicals very close by and in the same cavity site (given the computational
cost of higher quality calculations we cannot simulate much larger clusters).
However, in a more realistic situation each radical would be stored in
different cavities and thus the actual barriers to overcome would not only
involve breaking the radical/ice inter-molecular interactions, but also
surmounting the ice surface diffusional barriers, decreasing in this way the
efficiency of Rc and dHa reactions, even if they were ultimately barrierless.

We conclude this part mentioning that astronomical observations can also bring
useful constraints to the formation routes of iCOMs showing alternative routes
to the ones explored in this study. For example, high spatial resolution
observations of formamide line emission towards the protostellar shock site
L1157-B1 have demonstrated that the formation of formamide is dominated by the
gas-phase reaction NH2 \+ H2CO 65 , a reaction theoretically studied by some
of us 66, 67. On the same line, observations of the deuterated forms of
formamide (namely containing D rather than H atoms) also provide strong
constraints on the formation route of this species in the hot corino of the
solar-type protostar IRAS 16293-2422. The comparison of the measured
NHDHCO/NH2CHO and NH2CDO/NH2CHO abundance ratios33 with those predicted by
theoretical quantum chemical calculations67 strongly favors a gas-phase origin
of formamide also in this source. Therefore, it is very likely that both
grain-surface and gas-phase reactions contribute to the enrichment of iCOMs in
the ISM, playing different roles in different environments.

##  5 5\. Conclusions and perspectives

In this work, we have carried out an accurate study of the chemistry of two
couples of radicals, HCO + CH3 and HCO + NH2, on icy surfaces. Our goal was to
understand the possible reactions between the two radicals on water ice
mantles, and how the results depend on the accuracy of the employed quantum
chemical methods and on the adopted surface models. To this end, we used
different quantum chemistry methods, in particular two hybrid DFT methods,
B3LYP and BHLYP, plus the wave function based CCSD(T) and multi reference-
based CASPT2 ones. In addition, we adopted different cluster models simulating
the water surfaces: we started with the simple cases of one and two water
molecules to identify the basic processes and to test the methodology, and
then using two different, large molecular cluster models for the ASW surfaces,
of 18 and 33 water molecules, respectively.

The conclusions of this work are the following:
1\. If the reaction occurs, two channels are possible: (i) the combination of
radicals into acetaldehyde/formamide and (ii) the formation of CH4/NH3 plus
CO, where the H atom of HCO is passed to CH3/NH2 via H abstraction.
2\. The two reaction channels are either barrierless or have relatively low
energy barriers, from about 2 to 10 kJ/mol, as summarized in Table 3.
3\. Comparison of the results obtained with B3LYP-D3 and BHLYP (the latter in
its pure definition and including both D2 and D3 dispersion corrections) with
those provided by CASPT2 for activation energies and those provided by CCSD(T)
for binding energies, using one and two water molecules plus the radicals as
test systems, indicates that B3LYP-D3 underestimates the energy barriers,
while BHLYP-based methods show a reasonably good performance. For the
computations relative to the 18 and 33 water clusters, we adopted BHLYP-D3 as
it has been found, in the test systems, to properly deal with both the
radical/surface binding and the radical-radical activation energies.
4\. The morphology of the water cluster used for the simulations definitely
affects the results of the computations. In particular, radicals would
interact differently depending on whether they sit on a cavity structure,
where they can establish several weak inter-molecular interactions with the
icy water molecules, in addition to the H-bond.
5\. Taking into account the results described in points 1, 2 and 4, the
mechanism that radical combination necessarily produces iCOMs is still to be
validated, and should be taken with care in astrochemical models.

In order to make progresses, more accurate computations would be needed, but
they are not yet within the reach of the current computational capacities. On
the other hand, dynamical simulations would help to understand the effect of
the relative orientation of radicals upon encounter.

##  6 Supporting information

Structures and errors of the benchmark study, PESs of water assisted
H-transfer reactions on W18 and W33-cav, structures, binding energies and PESs
(of radical-radical coupling, direct H-abstraction and water assisted
H-transfer reactions) on the three analysed surfaces (W18, W33-side and
W33-cav) at BHLYP-D2 and on W18 at B3LYP-D3, and optimized Cartesian
coordinates of all the structures.

JER and CC acknowledge funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation program, for the
Project “the Dawn of Organic Chemistry” (DOC), grant agreement No 741002. AR
is indebted to “Ramón y Cajal” program. MINECO (project CTQ2017-89132-P) and
DIUE (project 2017SGR1323) are acknowledged. PU and NB acknowledge MIUR
(Ministero dell’Istruzione, dell’Università e della Ricerca) and from Scuola
Normale Superiore (project PRIN 2015, STARS in the CAOS - Simulation Tools for
Astrochemical Reactivity and Spectroscopy in the Cyberinfrastructure for
Astrochemical Organic Species, cod. 2015F59J3R). This project has received
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No 811312.
We thank Prof. Gretobape for fruitful and stimulating discussions.
Most of the calculations presented in this paper were performed using the
GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly
supported by the Equip@Meso project (reference ANR-10-EQPX-29-01) of the
programme Investissements d’Avenir supervised by the Agence Nationale pour la
Recherche. Additionally this work was granted access to the HPC resources of
IDRIS under the allocation 2019-A0060810797 attributed by GENCI (Grand
Equipement National de Calcul Intensif).

## References

  * Rubin et al. 1971 Rubin, R.; Swenson Jr, G.; Benson, R.; Tigelaar, H.; Flygare, W. Microwave detection of interstellar formamide. _Astrophys. J._ 1971, _169_ , L39
  * Cazaux et al. 2003 Cazaux, S.; Tielens, A.; Ceccarelli, C.; Castets, A.; Wakelam, V.; Caux, E.; Parise, B.; Teyssier, D. The Hot Core around the Low-Mass Protostar IRAS 16293–2422: Scoundrels Rule! _Astrophys. J. Lett._ 2003, _593_ , L51
  * Herbst and Van Dishoeck 2009 Herbst, E.; Van Dishoeck, E. F. Complex organic interstellar molecules. _Annu. Rev. Astron. Astrophys._ 2009, _47_ , 427–480
  * Caselli and Ceccarelli 2012 Caselli, P.; Ceccarelli, C. Our astrochemical heritage. _Astron. Astrophys. Rev._ 2012, _20_ , 56
  * Belloche et al. 2017 Belloche, A.; Meshcheryakov, A.; Garrod, R.; Ilyushin, V.; Alekseev, E.; Motiyenko, R.; Margulès, L.; Müller, H.; Menten, K. Rotational spectroscopy, tentative interstellar detection, and chemical modeling of N-methylformamide. _Astron. Astrophys._ 2017, _601_ , A49
  * McGuire 2018 McGuire, B. A. 2018 Census of Interstellar, Circumstellar, Extragalactic, Protoplanetary Disk, and Exoplanetary Molecules. _Astrophys. J. Suppl. Ser._ 2018, _239_ , 17
  * Bockelée-Morvan 2011 Bockelée-Morvan, D. An overview of comet composition. _Proc. Int. Astron. Union_ 2011, _7_ , 261–274
  * 8 Ceccarelli, C.; Caselli, P.; Bockelée-Morvan, D.; Mousis, O.; Pizzarello, S.; Robert, F.; Semenov, D. Deuterium Fractionation: The Ariadne’s Thread from the Precollapse Phase to Meteorites and Comets Today. In Protostars and Planets VI, ed. by Beuther, Henrik and Klessen, Ralf S. and Dullemond, Cornelis P. and Henning, Thomas (University of Arizona Press, Tucson), 2014, pp 859-882. DOI: 10.2458/azu_uapress_9780816531240-ch037
  * Bianchi et al. 2018 Bianchi, E.; Codella, C.; Ceccarelli, C.; Vazart, F.; Bachiller, R.; Balucani, N.; Bouvier, M.; De Simone, M.; Enrique-Romero, J.; Kahane, C.; Lefloch, B.; López-Sepulcre, A.; Ospina-Zamudio, J.; Podio, L.; Taquet, V. The census of interstellar complex organic molecules in the Class I hot corino of SVS13-A. _Mon. Not. R. Astron. Soc._ 2018, _483_ , 1850–1861
  * Ceccarelli et al. 2017 Ceccarelli, C. et al. Seeds Of Life In Space (SOLIS): The Organic Composition Diversity at 300–1000 au Scale in Solar-type Star-forming Regions. _Astrophys. J._ 2017, _850_ , 176
  * Charnley et al. 1992 Charnley, S.; Tielens, A.; Millar, T. On the molecular complexity of the hot cores in Orion A-Grain surface chemistry as ‘The last refuge of the scoundrel’. _Astrophys. J._ 1992, _399_ , L71–L74
  * Caselli et al. 1993 Caselli, P.; Hasegawa, T.; Herbst, E. Chemical differentiation between star-forming regions-The Orion hot core and compact ridge. _Astrophys. J._ 1993, _408_ , 548–558
  * Ceccarelli et al. 2001 Ceccarelli, C.; Loinard, L.; Castets, A.; Tielens, A.; Caux, E.; Lefloch, B.; Vastel, C. Extended D2CO emission: The smoking gun of grain surface-chemistry. _Astron. Astrophys._ 2001, _372_ , 998–1004
  * Geppert et al. 2005 Geppert, W. D.; Thomas, R. D.; Ehlerding, A.; Hellberg, F.; Österdahl, F.; Hamberg, M.; Semaniak, J.; Zhaunerchyk, V.; Kaminska, M.; Källberg, A.; Paal, A.; Larsson, M. Dissociative recombination branching ratios and their influence on interstellar clouds. _J. Phys. Conf. Ser._ 2005, _4_ , 26–31
  * Horn et al. 2004 Horn, A.; Møllendal, H.; Sekiguchi, O.; Uggerud, E.; Roberts, H.; Herbst, E.; Viggiano, A.; Fridgen, T. D. The gas-phase formation of methyl formate in hot molecular cores. _Astrophys. J._ 2004, _611_ , 605
  * Garrod and Herbst 2006 Garrod, R. T.; Herbst, E. Formation of methyl formate and other organic species in the warm-up phase of hot molecular cores. _Astron. Astrophys._ 2006, _457_ , 927–936
  * Herbst 2017 Herbst, E. The synthesis of large interstellar molecules. _Int. Rev. Phys. Chem._ 2017, _36_ , 287–331
  * Garrod 2008 Garrod, R. A new modified-rate approach for gas-grain chemical simulations. _Astron. Astrophys._ 2008, _491_ , 239–251
  * Acharyya and Herbst 2015 Acharyya, K.; Herbst, E. Molecular Development in the Large Magellanic Cloud. _Astrophys. J._ 2015, _812_ , 142
  * Ruaud et al. 2015 Ruaud, M.; Loison, J.; Hickson, K.; Gratier, P.; Hersant, F.; Wakelam, V. Modelling complex organic molecules in dense regions: Eley–Rideal and complex induced reaction. _Mon. Not. R. Astron. Soc._ 2015, _447_ , 4004–4017
  * Drozdovskaya et al. 2016 Drozdovskaya, M. N.; Walsh, C.; Van Dishoeck, E. F.; Furuya, K.; Marboeuf, U.; Thiabaud, A.; Harsono, D.; Visser, R. Cometary ices in forming protoplanetary disc midplanes. _Mon. Not. R. Astron. Soc._ 2016, _462_ , 977–993
  * Vasyunin et al. 2017 Vasyunin, A. I.; Caselli, P.; Dulieu, F.; Jiménez-Serra, I. Formation of complex molecules in prestellar cores: a multilayer approach. _Astrophys. J._ 2017, _842_ , 33
  * Rimola et al. 2018 Rimola, A.; Skouteris, D.; Balucani, N.; Ceccarelli, C.; Enrique-Romero, J.; Taquet, V.; Ugliengo, P. Can Formamide Be Formed on Interstellar Ice? An Atomistic Perspective. _ACS Earth Space Chem._ 2018, _2_ , 720–734
  * Enrique-Romero et al. 2016 Enrique-Romero, J.; Rimola, A.; Ceccarelli, C.; Balucani, N. The (impossible?) formation of acetaldehyde on the grain surfaces: insights from quantum chemical calculations. _Mon. Not. R. Astron. Soc._ 2016, _459_ , L6–L10
  * Balucani et al. 2015 Balucani, N.; Ceccarelli, C.; Taquet, V. Formation of complex organic molecules in cold objects: the role of gas-phase reactions. _Mon. Not. R. Astron. Soc._ 2015, _449_ , L16–L20
  * Spezia et al. 2016 Spezia, R.; Jeanvoine, Y.; Hase, W. L.; Song, K.; Largo, A. Synthesis of Formamide and Related Organic Species in the Interstellar Medium via Chemical Dynamics Simulations. _Astrophys. J._ 2016, _826_ , 107
  * Taquet et al. 2016 Taquet, V.; Wirström, E. S.; Charnley, S. B. Formation and Recondensation of Complex Organic Molecules during Protostellar Luminosity Outbursts. _Astrophys. J._ 2016, _821_ , 46
  * Suzuki et al. 2016 Suzuki, T.; Ohishi, M.; Hirota, T.; Saito, M.; Majumdar, L.; Wakelam, V. Survey Observations of a Possible Glycine Precursor, Methanimine (CH2NH). _Astrophys. J._ 2016, _825_ , 79
  * Codella, C. et al. 2017 Codella, C., et al. Seeds of Life in Space (SOLIS) - II. Formamide in protostellar shocks: Evidence for gas-phase formation. _Astron. Astrophys._ 2017, _605_ , L3
  * Skouteris et al. 2018 Skouteris, D.; Balucani, N.; Ceccarelli, C.; Faginas Lago, N.; Codella, C.; Falcinelli, S.; Rosi, M. Interstellar dimethyl ether gas-phase formation: a quantum chemistry and kinetics study. _Mon. Not. R. Astron. Soc._ 2018, _482_ , 3567–3575
  * Rosi et al. 2018 Rosi, M.; Skouteris, D.; Casavecchia, P.; Falcinelli, S.; Ceccarelli, C.; Balucani, N. Formation of Nitrogen-Bearing Organic Molecules in the Reaction NH $+$ C2H5: A Theoretical Investigation and Main Implications for Prebiotic Chemistry in Space. Computational Science and Its Applications $-$ ICCSA 2018 (Melbourne, Australia). Springer International Publishing: Cham (Switzerland), 2018; pp 773–782
  * Skouteris et al. 2018 Skouteris, D.; Balucani, N.; Ceccarelli, C.; Vazart, F.; Puzzarini, C.; Barone, V.; Codella, C.; Lefloch, B. The Genealogical Tree of Ethanol: Gas-phase Formation of Glycolaldehyde, Acetic Acid, and Formic Acid. _Astrophys. J._ 2018, _854_ , 135
  * Coutens et al. 2016 Coutens, A. et al. The ALMA-PILS survey: First detections of deuterated formamide and deuterated isocyanic acid in the interstellar medium. _Astron. Astrophys._ 2016, _590_ , L6
  * Müller et al. 2016 Müller, H. S. P.; Belloche, A.; Xu, L.-H.; Lees, R. M.; Garrod, R. T.; Walters, A.; van Wijngaarden, J.; Lewen, F.; Schlemmer, S.; Menten, K. M. Exploring molecular complexity with ALMA (EMoCA): Alkanethiols and alkanols in Sagittarius B2(N2). _Astron. Astrophys._ 2016, _587_ , A92
  * Ligterink et al. 2018 Ligterink, N. F. W.; Terwisscha van Scheltinga, J.; Taquet, V.; Jørgensen, J. K.; Cazaux, S.; van Dishoeck, E. F.; Linnartz, H. The formation of peptide-like molecules on interstellar dust grains. _Mon. Not. R. Astron. Soc._ 2018, _480_ , 3628–3643
  * Strazzulla and Palumbo 2001 Strazzulla, G.; Palumbo, M. Organics produced by ion irradiation of ices: Some recent results. _Adv. Space Res._ 2001, _27_ , 237 – 243
  * Palumbo et al. 2008 Palumbo, M. E.; Baratta, G. A.; Fulvio, D.; Garozzo, M.; Gomis, O.; Leto, G.; Spinella, F.; Strazzulla, G. Ion irradiation of astrophysical ices. _J. Phys. Conf. Ser._ 2008, _101_ , 012002
  * Woods et al. 2013 Woods, P. M.; Slater, B.; Raza, Z.; Viti, S.; Brown, W. A.; Burke, D. J. Glycolaldehyde formation via the dimerization of the formyl radical. _Astrophys. J._ 2013, _777_ , 90
  * Fedoseev et al. 2014 Fedoseev, G.; Ioppolo, S.; Zhao, D.; Lamberts, T.; Linnartz, H. Low-temperature surface formation of NH3 and HNCO: hydrogenation of nitrogen atoms in CO-rich interstellar ice analogues. _Mon. Not. R. Astron. Soc._ 2014, _446_ , 439–448
  * Fedoseev et al. 2015 Fedoseev, G.; Cuppen, H. M.; Ioppolo, S.; Lamberts, T.; Linnartz, H. Experimental evidence for glycolaldehyde and ethylene glycol formation by surface hydrogenation of CO molecules under dense molecular cloud conditions. _Mon. Not. R. Astron. Soc._ 2015, _448_ , 1288–1297
  * Linnartz et al. 2015 Linnartz, H.; Ioppolo, S.; Fedoseev, G. Atom addition reactions in interstellar ice analogues. _Int. Rev. Phys. Chem._ 2015, _34_ , 205–237
  * de Barros et al. 2016 de Barros, A. L. F.; da Silveira, E. F.; Fulvio, D.; Boduch, P.; Rothard, H. Formation of nitrogen- and oxygen-bearing molecules from radiolysis of nitrous oxide ices – implications for Solar system and interstellar ices. _Mon. Not. R. Astron. Soc._ 2016, _465_ , 3281–3290
  * Butscher et al. 2017 Butscher, T.; Duvernay, F.; Rimola, A.; Segado-Centellas, M.; Chiavassa, T. Radical recombination in interstellar ices, a not so simple mechanism. _Phys. Chem. Chem. Phys._ 2017, _19_ , 2857–2866
  * Fedoseev et al. 2017 Fedoseev, G.; Chuang, K.-J.; Ioppolo, S.; Qasim, D.; van Dishoeck, E. F.; Linnartz, H. Formation of Glycerol through Hydrogenation of CO Ice under Prestellar Core Conditions. _Astrophys. J._ 2017, _842_ , 52
  * Arumainayagam et al. 2019 Arumainayagam, C. R.; Garrod, R. T.; Boyer, M. C.; Hay, A. K.; Bao, S. T.; Campbell, J. S.; Wang, J.; Nowak, C. M.; Arumainayagam, M. R.; Hodge, P. J. Extraterrestrial prebiotic molecules: photochemistry vs. radiation chemistry of interstellar ices. _Chem. Soc. Rev._ 2019, _48_ , 2293–2314
  * Dulieu et al. 2019 Dulieu, F.; Nguyen, T.; Congiu, E.; Baouche, S.; Taquet, V. Efficient formation route of the prebiotic molecule formamide on interstellar dust grains. _Mon. Not. R. Astron. Soc.: Letters_ 2019, _484_ , L119–L123
  * Pereira et al. 2019 Pereira, R.; de Barros, A.; Fulvio, D.; Boduch, P.; Rothard, H.; da Silveira, E. N–O bearing molecules produced by radiolysis of N2O and N2O:CO2 ices. 2019, Nucl. Instrum. Methods Phys. Res., Sect. B. DOI: 10.1016/j.nimb.2019.03.048
  * Lamberts et al. 2019 Lamberts, T.; Markmeyer, M. N.; Kolb, F. J.; Kästner, J. Formation of Acetaldehyde on CO-Rich Ices. _ACS Earth Space Chem._ 2019, _3_ , 958–963
  * Abe 2013 Abe, M. Diradicals. _Chem. Rev._ 2013, _113_ , 7011–7088
  * Frisch et al. 2009 Frisch, M. J. et al. Gaussian09 Revision D.01. 2009; Gaussian Inc. Wallingford CT 2009
  * Karlström et al. 2003 Karlström, G.; Lindh, R.; Malmqvist, P.-Å.; Roos, B. O.; Ryde, U.; Veryazov, V.; Widmark, P.-O.; Cossi, M.; Schimmelpfennig, B.; Neogrady, P.; Seijo, L. MOLCAS: a program package for computational chemistry. _Comput. Mater. Sci._ 2003, _28_ , 222 – 239, Proceedings of the Symposium on Software Development for Process and Materials Design
  * Aquilante, Francesco et al. 2010 Aquilante, Francesco,; De Vico, Luca,; Ferré, Nicolas,; Ghigo, Giovanni,; Malmqvist, Per-åke,; Neogrády, P.; Pedersen, T. B.; Pitoňák, M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. MOLCAS 7: The Next Generation. _J. Comput. Chem._ 2010, _31_ , 224–247
  * Aquilante, Francesco et al. 2016 Aquilante, Francesco, et al. Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table. _J. Comput. Chem._ 2016, _37_ , 506–541
  * Fernandez Galván et al. 2019 Fernandez Galván, I.; Vacher, M.; Alavi, A.; Angeli, C.; Autschbach, J.; Bao, J. J.; Bokarev, S. I.; Bogdanov, N. A.; Carlson, R. K.; Chibotaru, L. F.; et al., OpenMolcas: From Source Code to Insight. 2019; https://chemrxiv.org/articles/OpenMolcas_From_Source_Code_to_Insight/8234021/1
  * Lee et al. 1988 Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. _Phys. Rev. B_ 1988, _37_ , 785
  * Becke 1993 Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. _J. Chem. Phys._ 1993, _98_ , 5648–5652
  * Becke 1993 Becke, A. D. A new mixing of Hartree–Fock and local density-functional theories. _J. Chem. Phys._ 1993, _98_ , 1372–1377
  * Grimme et al. 2010 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _J. Chem. Phys._ 2010, _132_ , 154104
  * Grimme 2006 Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _J. Comput. Chem._ 2006, _27_ , 1787–1799
  * Rimola et al. 2014 Rimola, A.; Taquet, V.; Ugliengo, P.; Balucani, N.; Ceccarelli, C. Combined quantum chemical and modeling study of CO hydrogenation on water ice. _Astron. Astrophys._ 2014, _572_ , A70
  * McQuarrie 1976 McQuarrie, D. _Statistical Mechanics_ ; New York: Harper and Row, 1976
  * Boys and Bernardi 1970 Boys, S. F.; Bernardi, F. d. The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. _Mol. Phys._ 1970, _19_ , 553–566
  * Wakelam et al. 2017 Wakelam, V.; Loison, J.-C.; Mereau, R.; Ruaud, M. Binding energies: new values and impact on the efficiency of chemical desorption. _Mol. Astrophys._ 2017, _6_ , 22–35
  * Sameera et al. 2017 Sameera, W. M. C.; Senevirathne, B.; Andersson, S.; Maseras, F.; Nyman, G. ONIOM(QM:AMOEBA09) Study on Binding Energies and Binding Preference of OH, HCO, and CH3 Radicals on Hexagonal Water Ice (Ih). _J. Phys. Chem. C._ 2017, _121_ , 15223–15232
  * Codella et al. 2017 Codella, C. et al. Seeds of Life in Space (SOLIS). II. Formamide in protostellar shocks: Evidence for gas-phase formation. _Astron. Astrophys._ 2017, _605_ , L3
  * Barone et al. 2015 Barone, V.; Latouche, C.; Skouteris, D.; Vazart, F.; Balucani, N.; Ceccarelli, C.; Lefloch, B. Gas-phase formation of the prebiotic molecule formamide: insights from new quantum computations. _Mon. Not. R. Astron. Soc._ 2015, _453_ , L31–L35
  * Skouteris et al. 2017 Skouteris, D.; Vazart, F.; Ceccarelli, C.; Balucani, N.; Puzzarini, C.; Barone, V. New quantum chemical computations of formamide deuteration support gas-phase formation of this prebiotic molecule. _Mon. Not. R. Astron. Soc._ 2017, _468_ , L1–L5

##  7 For TOC Only