# Fermi-Bose crossover using engineered disorder potential

Madhuparna Karmakar  Department of Physics and Nanotechnology, SRM Institute
of Science and Technology,
Kattankulathur, Chennai 603203, India<EMAIL_ADDRESS>Rajesh
Narayanan  Department of Physics, Indian Institute of Technology, Madras,
Chennai 600036, India.

###### Abstract

We present the first instance of a disorder tuned Fermi-Bose crossover that
could be realized in superconducting systems. More specifically, harnessing a
non perturbative numerical technique we analyze the ground state behavior of a
two-dimensional attractive Hubbard model subjected to spin selective disorder
potential. In particular, using spectroscopic properties we provide
unambiguous evidence of the change in the Fermi surface topology as a function
of the disorder, establishing incontrovertibly a Fermi-Bose crossover.
Interplay of strong correlations and strong disorder brings out the spin
selectivity in the properties of this system giving rise to spin selective
“Bose metal/insulator” phase. We propose an experimental set-up where such
disorder tuned Fermi-Bose crossover could be observed in two-dimensional
electron gas formed at oxide interface. Finally, we speculate on the possible
implication of such spin-selective disorder on unravelling signatures of Bose-
Fermi cross-over on doped iron chalcogenide superconductors.

Introduction: The Bardeen-Cooper-Schrieffer (BCS) superfluid phase at weak
coupling is continuously connected to the strongly coupled phase of the Bose-
Einstein condensate (BEC) Randeria and Taylor (2014) through a cross-over
which can be tuned via various control parameters. Of particular interest is
the cross-over regime at intermediate interactions, where neither a completely
fermionic nor a bosonic description of the system holds good. Intriguing
physics pertaining to the pseudogap phase, preformed Cooper pairs and strange
metals are reported to play out in this regime, attributed largely to the
short range pair correlations Emery and Kivelson (1994); Corson  _et al._ ; Li
_et al._ (2010); Keimer  _et al._ . The advent of ultracold atomic gases made
such cross-overs experimentally accessible via the continuous tunability
provided by the Feshbach resonance W. Ketterle (2008); Stewart  _et al._ ;
Gaebler  _et al._ ; Chin  _et al._ (2004, 2010); Giorgini  _et al._ (2008).

The Fermi-Bose crossover has been recently realized experimentally in solid
state systems such as, iron based superconductors Kasahara  _et al._ (2016);
Rinott  _et al._ (2017); Hashimoto  _et al._ (2020); Lubashevsky  _et al._
(2012). Furthermore, it has been predicted that transition metal
dichalcogenide hetero-structures Crépel  _et al._ (2023), could serve as an
ideal platform to realize the BEC-BCS cross-over in a topological p-wave
superconductor. Apart from the Feshbach resonance, other novel methods to
induce the Fermi-Bose crossover have been proposed which relies on the tuning
of the interlayer hopping of a prototypical bilayer triangular lattice via
strain engineering Loh  _et al._ (2016); Hazra  _et al._ (2020).

Figure 1: Ground state phase diagram in the interaction-disorder
($U-V$)-plane, showing the thermodynamic phases viz. (i) gapped
superconductor, (ii) gapless superconductor, (iii) spin selective Anderson
insulator and (iv) spin selective Bose insulator/metal. The dashed curve
$V_{Bose}$ is the Fermi-Bose crossover scale. The solid curves $V_{g}$ and
$V_{sc}$ illustrate crossover scales between the gapped and gapless
superconductors and superconductor- insulator transition, respectively. See
text for further details.

In this work we propose an alternate protocol wherein the Fermi-Bose crossover
is achieved via the manipulation of the random disorder potential. More
specifically, we tune the fermionic dispersion via the implementation of a
“spin-selective” disorder wherein the different spin species are subjected to
different disorder distributions. Apart from other indicators, we show that
the underlying Fermi surface in the BCS superconductor disappears as one tunes
through the cross-over into the BEC limit as a function of the spin-selective
disorder strength. This change in the topology of the Fermi surface is a
salient feature of the Fermi-Bose cross-over.

The impact of quenched disorder, in particular that of spin independent
disorder potential on superconductors, is a well investigated subject. These
investigations have unravelled, for instance the mechanism behind the
localization phenomena Sacépé  _et al._ (2011, 2010) as well as the nature of
the disorder induced superconductor-insulator transition (SIT), in low
dimensional systems of $d<3$ Fisher (1990); Vojta (2006, 2013, 2019); Hoyos
_et al._ (2007); Vojta  _et al._ (2009); Ghosal  _et al._ (2001); Dubi  _et
al._ (2007). Furthermore, these studies have revealed the existence of an
anomalous metallic phase in certain low dimensional disordered superconductors
(see for example Ref. Kapitulnik  _et al._ (2019), for a recent review), and
also to the observation of exotic disorder induced Griffiths effects in such
systems Kaur  _et al._ (2022); Xing  _et al._ (2015); Shen  _et al._ (2016);
Lewellyn  _et al._ (2019); Reiss  _et al._ (2021). In contrast, barring a few
exceptions Jiang  _et al._ (2013); Nanguneri  _et al._ (2012); Karmakar  _et
al._ (2022), the impact of spin-selective disorder on the superconducting
transition has not been so well investigated. In particular, Ref. Karmakar
_et al._ (2022), have fully elucidated the phases and the ensuing phase
transitions inherent in a spin-selectively disordered two dimensional
superconductor, conclusively demonstrating that apart from the usual
superconducting and normal phases, the system hosts a breached pair phase
(away from the half-filling), wherein superconductivity spatially co-exists
with magnetism.

One such system where a spin-selective disorder tuned Fermi-Bose cross-over
can be observed are ultracold atomic gases hosted on an optical lattice McKay
and DeMarco (2010); Mandel  _et al._ (2003a); Feiguin and Fisher (2009); Liu
_et al._ (2004); Mandel  _et al._ (2003b), wherein the disorder is rendered
spin-selective via the introduction of speckle disorder Jiang  _et al._
(2013); Nanguneri  _et al._ (2012); Skolimowski  _et al._ (2015). Another
avenue that we propose where such cross-overs maybe observed is in the two-
dimension electron gas (2DEG) formed at the interface of oxide
heterostructure: for instance, at the oxide interface formed between non
magnetic band insulators LaAlO3 (LAO) and SrTiO3 (STO) Ohtomo and Hwang
(2004); Zhang  _et al._ (2019); Han  _et al._ (2019); Xu  _et al._ ; Niu  _et
al._ (2017); Chen  _et al._ (2011).

Magneto resistance measurements on these systems firmly established the
existence of ferromagnetism in the 2DEG Li  _et al._ (2011); Dikin  _et al._
(2011). Further, more sensitive techniques such as, magnetic mapping by
scanning SQUID Bert  _et al._ (2011), magnetic force microscopy Bi  _et al._
(2014) etc. established the inhomogeneous spatial coexistence of
ferromagnetism with superconducting correlations. Such inhomogeneous local
magnetic polarization can be modeled via an inherently random local Zeeman
field. An interplay between this random local Zeeman-field and the
conventional chemical potential disorder that obtains due to the presence of
quenched impurities result in an asymmetric renormalization of the disorder
potential seen by the two spin species thus, leading to spin-selectivity of
the disorder. Such asymmetric renormalization of the disorder potential has
also for instance been reported in the material Fe1-yCoySi wherein the spin-
splitting due to the exchange interaction leads to the fact that the depth of
the potential wells at the Co sites is increased for the majority carriers and
decreased for the minority carriers Mena  _et al._ (2006).

Our key results are: (i) We propose a novel protocol for tuning the Fermi-Bose
cross-over in solid state superconductors using a tailored spin-selective
disorder potential. (ii) The cross-over is mapped out in terms of the Fermi
surface reconstruction quantified via the spectroscopic signatures. (iii) We
demonstrate that the interplay of strong interaction and disorder brings forth
a spin-selective Bose metal/insulator phase.

Model and Methods: Our starting point to understand the spin-selective
disorder induced Fermi-Bose cross-over is the attractive Hubbard model defined
on a square lattice, described by the Hamiltonian:

$\displaystyle H$ $\displaystyle=$ $\displaystyle-t\sum_{\langle
ij\rangle,\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+c_{j\sigma}^{\dagger}c_{i\sigma})-\mid
U\mid\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}$ (1)
$\displaystyle+\sum_{i\sigma}(V_{i\sigma}-\mu)\hat{n}_{i\sigma}$

Here, the hopping amplitude is set at $t=1$ for the nearest neighbors, which
is the reference energy scale of the system. The spin-selectivity comes in
making the random $V_{i\sigma}$, spin dependent. In our implementation of spin
selectivity we have chosen our disorder such that $V_{i\uparrow}=0$ and
$V_{i\downarrow}\neq 0=V$, selected randomly from a box distribution of
$[-V/2,V/2]$. $|U|>0$ is the attractive interaction which gives rise to the
$s$-wave superconducting pairing. We work in the limit of the half-filled
lattice ($n=1$), wherein for the translationally invariant system there exists
a degeneracy between superconductivity and charge density wave (CDW) orders.
Inclusion of disorder immediately lifts this degeneracy and makes the
superconducting state energetically favorable over the CDW. We tune the global
chemical potential $\mu$ to stay in the limit of the half filled lattice.

As a first step to study Eq. 1 numerically we decompose the interaction term
via the Hubbard-Stratonovich (HS) decomposition. Specifically, we introduce a
complex $\Delta_{i}(\tau)=|\Delta_{i}(\tau)|e^{i\theta_{i}(\tau)}$ and a real
$\phi_{i}(\tau)$ scalar (bosonic) auxiliary fields at each site, where the
former couples to the pairing and the later to the charge channels Hubbard
(1959); Schulz (1990). We numerically address this problem by employing a path
integral based Monte Carlo technique viz. static path approximated (SPA)
quantum Monte Carlo wherein we drop the time dependence of the auxiliary
fields and treat them as “classical” $\Delta_{i}$ and $\phi_{i}$ Karmakar and
Majumdar (2016); Evenson  _et al._ (1970); Dubi  _et al._ (2007). This is
tantamount to retaining only the $\Omega_{n}=0$, Matsubara frequency mode of
the fluctuations of the auxiliary fields. The spatial dependence is however,
treated exactly by retaining the corresponding spatial fluctuations at all
orders, i. e. $\Delta_{i}(\tau)\rightarrow|\Delta_{i}|e^{i\theta_{i}}$ and
$\phi_{i}(\tau)\rightarrow\phi_{i}$, where, $|\Delta_{i}|$ is the
superconducting pairing field amplitude and $\theta_{i}$ is the corresponding
phase. These approximations lead to a coupled spin-fermion model with the
fermions moving in the spatially fluctuating random background of
$\\{\Delta_{i},\phi_{i}\\}$, which obeys the Boltzmann distribution,
$P\\{\Delta_{i},\phi_{i}\\}=\mathrm{Tr}_{c^{\dagger},c}e^{\beta H_{eff}}$,
where $\beta$ is the inverse temperature and $H_{eff}$ is the effective
Hamiltonian of the spin-fermion model. For large and random background the
trace is computed numerically. The equilibrium configurations of
$\\{\Delta_{i},\phi_{i}\\}$ are generated via classical Monte Carlo
simulation, diagonalizing $H_{eff}$ for each attempted update of
$\\{\Delta_{i},\phi_{i}\\}$. The fermionic correlators are computed on these
equilibrium configurations of the auxiliary fields.

The system is characterized using the following indicators: $(i)$
superconducting pairing field structure factor ($S({\bf
q})=(1/N^{2})\sum_{ij}\langle\Delta_{i}\Delta_{j}^{*}\rangle e^{i{\bf q}.({\bf
r}_{i}-{\bf r}_{j})}$), $(ii)$ spin resolved single particle density of states
(DOS)
($N_{\uparrow}(\omega)=(1/N)\langle\sum_{i,n}|u_{n}^{i}|^{2}\delta(\omega-
E_{n})\rangle$,
$N_{\downarrow}(\omega)=(1/N)\langle\sum_{i,n}|v_{n}^{i}|^{2}\delta(\omega+E_{n})\rangle$)
and the associated energy gap $E^{\sigma}_{g}$, $(iii)$ momentum and spin
resolved spectral function ($A^{\sigma}({\bf k},\omega)=(-1/\pi){\rm
Im}G^{\sigma}({\bf k},\omega)$), where, $u_{n}^{i}$ and $v_{n}^{i}$ are the
Bogoliubov eigenvectors with the eigenvalues $E_{n}$ and $\sigma$ is the spin
label. The results presented here corresponds to a system size of $24\times
24$ and are found to be robust against finite system size effects. The
observables are averaged over $10$ disorder realizations and are found to be
invariant to the increase in disorder realizations.

The phase diagram: Our main results predicated based on the physical
observables described above is encapsulated in the ground state phase diagram
shown in Fig. 1: the translationally invariant or clean system clearly is a
gapped superconductor. There are two cross-over scales labelled as $V_{\rm
g}$, and $V_{\rm Bose}$. In the context of this manuscript, one of the most
relevant scale is $V_{\rm Bose}$, which quantifies the disorder strength where
the Fermi superconductor described by a BCS-like state continuously evolves
into a strongly coupled BEC state. More specifically, across the cross-over
scale set by $V_{\rm Bose}$ the fermionic dispersion spectra, mapped out in
terms of the spectral function (see Fig 2) shows a change in Fermi surface
topology: from being downward convex with an energy gap minima at a finite
wave vector (${\bf k}={\bf k}_{\rm{F}}$) in the BCS regime, the dispersion
branches get modified to being flat, with a gap minima at ${\bf k}=0$ in the
BEC regime Loh  _et al._ (2016). This behavior of the spectral function
showing the BCS-BEC cross-over is accompanied by similar signatures in the
single particle DOS: In line with Ref. Loh  _et al._ (2016), as we tune
through the cross-over, the (DOS) sharp van Hove singularities at the
superconducting gap edges give way to prominent discontinuities, as the system
crosses over from a fermionic to a bosonic state.

The cross-over scale $V_{\rm g}$ is weakly spin-selective and typifies the
disorder strength at which the system evolves from a fully gapped
superconductor to a gapless superconductor. The physics of this cross-over is
very much akin to the effect of magnetic impurities Abrikosov and Gor’kov
(1961); Jiang  _et al._ (2013); Nanguneri  _et al._ (2012); Karmakar  _et al._
(2022) on conventional $s$-wave superconductors. The gapped phase is
characterized with a finite spectral gap $E_{g}^{\sigma}$, and a non-vanishing
pairing field structure factor $S({\bf q})$. However, in the gapless phase the
spectral gap $E^{\sigma}_{g}=0$, and the $S({\bf q})\neq 0$.

This gapless superconductor is fermionic in nature below the cross-over scale
$V_{\rm Bose}$, while for $V\geq V_{\rm Bose}$ it evolves into a Bose
superconductor. Finally, the disorder scale $V_{\rm SC}$ encodes the impurity
strength at which the transition from the superconducting state to the normal
state occurs. Above this scale the pairing field structure factor $S({\bf
q})\rightarrow 0$, indicating the loss of global superconducting phase
coherence. This normal state is a gapless correlated Anderson insulator at
small interactions which adiabatically morphs into a gapped Bose
metal/Insulator with increasing interaction strength. In what follows, we
discuss the various facets of the data which underpins the construction of the
phase diagram shown in Fig. 1.

Figure 2: Fermi surface topology at the ground state, at representative $U-V$
cross sections, corresponding to fermionic ($U=4t,V=0.25t$) and bosonic
($U=7t,V=t$) superconducting states. (a) Single particle DOS for the
$\uparrow$-spin species ($N_{\uparrow}(\omega)$) for fermionic superconducting
phase showing prominent van Hove singularities at the gap edges, (b) Single
particle DOS for the bosonic superconducting phase showing gap edge
discontinuities, (c) corresponding spectral function ($A^{\uparrow}({\bf
k},\omega)$) showing the characteristic back-bending of the Bogoliubov quasi
particle spectra with the gap minima at ${\bf k}={\bf k}_{F}$, signifying the
presence of underlying Fermi surface. (d) corresponding spectral function with
flat dispersive branches having gap minima at ${\bf k}=0$, signifying the
disappearance of the Fermi surface. Figure 3: Disorder tuned Fermi-Bose
crossover as demonstrated via single particle spectral function
($A^{\uparrow}(\omega)$) at the ground state for $U=4t$. The weak ($V=0.25t$)
and strong ($V=2.25t$) limits are characteristic to the fermionic and bosonic
phases, while at intermediate disorder ($V=1.25t$) a (gapless) mixed phase is
realized with spin dependent fermionic and bosonic characteristics. Figure 4:
Thermodynamic and spectroscopic indicators at the ground state. (a) Disorder
dependence of the spectral gap ($E_{g}^{\uparrow}$) as determined from the
corresponding single particle DOS at different interactions. The scale $V_{\rm
g}$ is set by the gap closing condition $E_{g}^{\uparrow}\rightarrow 0$ (b)
Disorder dependent suppression of the superconducting pairing field structure
factor ($S(0,0)$) corresponding to the loss of (quasi) long range phase
coherence ($V_{sc}$). (c)-(d) Disorder dependence of the single particle DOS
($N_{\uparrow}(\omega)$) at the representative interactions of $U=4t$ and
$U=7t$ corresponding to the weak and strong couplings. Note the systematic
crossover of the gap edges from van Hove singularities at weak disorder to gap
edge discontinuities at strong disorder, as the system undergoes disorder
tuned Fermi-Bose crossover. (see text for further details).

Fermi surface reconstruction and Fermi-Bose crossover: Here we highlight the
Fermi-Bose crossover in detail by using the disorder averaged $(i)$ spin
resolved single particle DOS ($N_{\sigma}(\omega)$) and $(ii)$ spin resolved
spectral function ($A^{\sigma}({\bf k},\omega)$), shown in Fig. 2. In
particular, Fig. 2 compares the representative behavior of the spin-resolved
DOS $N_{\uparrow}(\omega)$ and the spin-resolved spectral function
$A^{\uparrow}({\bf k},\omega)$ in the Fermi superconducting phase (for example
at, $U=4t$, $V=0.25t$), (see Fig. 2(a) and 2(c) respectively), with their
behavior in the Bose superconducting phase, (see Fig. 2(b) and Fig. 2(d))
(exemplified by the situation at $U=7t$, $V=t$). The Fermi superconducting
phase is characterized by a robust spectral gap at the Fermi level with sharp
van Hove singularities (inverse square root divergence) at the gap edges, as
shown in Fig. 2(a). The corresponding spectral function shown in Fig. 2(c)
contains the signatures of the usual Bogoliubov spectrum comprising of the
dispersion branches $\pm E_{\bf k}$ with weights $u_{\bf k}^{2}$ and $v_{\bf
k}^{2}$ ($u_{\bf k}$ and $v_{\bf k}$ are the eigenvectors corresponding to the
energy eigenvalues $\pm E_{\bf k}$). The spectral function shows that the
minimum energy gap occurs on a contour at a finite wave vector ${\bf k}={\bf
k}_{\rm F}$, indicating the presence of an underlying Fermi surface.

On the other hand, the spectral function corresponding to the Bose
superconducting phase (Figure 2(d)) comprises of a single point (instead of a
locus) at the wave vector ${\bf k}=0$, where the gap minima occurs, indicating
the absence of a Fermi surface. The corresponding single particle DOS, shown
in Figure 2(b), no longer hosts van Hove singularities, and rather exhibits
discontinuity or kink at the gap edges. This somewhat smoothened (due to the
effect of quenched disorder), discontinuity is very reminiscent of the jump in
the DOS predicted to occur in the parameter window where a bosonic description
obtains Ref. Loh  _et al._ (2016). The system thus undergoes a disorder tuned
Fermi-Bose crossover, across which the Fermi surface undergoes a change in its
topology. The spectroscopic signatures characteristic of this crossover can be
experimentally observed via radio frequency (rf) and angle resolved photo
emission spectroscopy (ARPES) measurements.

The systematic evolution of the Fermi surface topology with the disorder
potential at $U=4t$ is shown in Figure 3, in terms of the fermionic spectral
function $A^{\uparrow}({\bf k},\omega)$. As a function of increasing disorder
the spectrum undergoes cross-over from being downward convex (minima at ${\bf
k}={\bf k}_{\rm F}$) to flat (minima at ${\bf k}=0$) as the underlying Fermi
surface vanishes. At intermediate disorder potential a gapless superconducting
state is realized; for the corresponding dispersion spectra neither a
completely fermionic nor a completely bosonic description holds good. The
system exhibits a spin-selective behavior such that, the $\uparrow$-spin
species show signatures of bosonic correlations in its spectroscopic
properties while the $\downarrow$-spin species behaves akin to a fermionic
system. We observe that the disorder tuned evolution of the Fermi surface
topology is in agreement with the one reported recently in ARPES measurements
across the Fermi-Bose crossover in doped Fe-based superconductor, FeSe1-xSx
Hashimoto  _et al._ (2020).

The energy scales $V_{\rm g}$ and $V_{\rm SC}$: We now use the disorder
dependence of $E^{\uparrow}_{g}$ ($E^{\downarrow}_{g}$ behaves similarly) and
$S({\bf q}=0)$ at different interactions (as shown in Fig. 4(a) and 4(b),
respectively), to capture the scales $V_{\rm g}$, and $V_{\rm SC}$,
respectively. In Fig. 4(a), we plot the behavior of $E^{\uparrow}_{g}$ as a
function of the disorder $V/t$, for a range of interactions. The energy scale
$V_{\rm g}$ is gleaned from Fig. 4(a) as the disorder strength (for a
particular value of $U/t$) at which the $E^{\uparrow}_{g}\rightarrow 0$. The
$E^{\uparrow}_{g}$ is tied to the pairing field amplitude ($|\Delta_{i}|$)
which monotonically increases with increasing $U/t$ and this monotonic
increase is reflected on the scale $V_{\rm g}$, which also shows concomitant
monotonic increase with $U/t$, (see Fig. 1).

To extract the energy scale $V_{\rm SC}$, we lay recourse to the pairing field
structure factor $S({\bf q})$. In Fig. 4(b), we depict the behavior of $S({\bf
q})$ as a function of $V/t$ for various interaction strength $U/t$. The scale
$V_{\rm SC}$ is captured by the value of disorder strength at which the
pairing field structure factor vanishes for each value of $U/t$. From Fig.
4(b), we immediately see that $S({\bf q})$ and consequently $V_{\rm SC}$
display a non-monotonic behavior as a function of interaction strength. This
non-monotonicity of the $S({\bf q})$, and subsequently $V_{\rm SC}$,
(especially in comparison with with increasing $E^{\uparrow}_{g}$), as a
function of $U/t$ is predicated on the different origins that control the loss
of superconducting order in the weak coupling and the intermediate and strong
coupling limits: In the weak-coupling limit, the loss of the superconducting
order is tied to the suppression of the pairing field amplitude, (thus
implying that $E^{\uparrow}_{g}$, is a suitable indicator), whereas at
intermediate and strong interactions it is the phase fluctuations that are
responsible for the suppression of the (quasi) long range order.

Figure 5: Real space maps corresponding to the superconducting pairing field
amplitude ($|\Delta_{i}|$), superconducting phase correlation between a
reference site ($0$) and any other lattice site
($\cos(\theta_{0}-\theta_{i})$) and the charge field ($\phi_{i}$), at selected
interaction-disorder ($U-V$) cross sections. The spatial maps characterize the
excursion of the system from a gapped Fermi superconductor $\rightarrow$
gapless Fermi superconductor $\rightarrow$ gapless Bose superconductor
$\rightarrow$ correlated Anderson insulator, as a function of disorder in the
weak coupling regime ($U=4t$). In the representative strong coupling regime
($U=7t$), the system follows the trajectory of a gapped Fermi superconductor
$\rightarrow$ gapped Bose superconductor $\rightarrow$ spin selective Bose
metal/insulator $\rightarrow$ correlated Anderson insulator, as the function
of disorder.

The behavior of the single particle DOS ($N_{\uparrow}(\omega)$) across the
disorder tuned Fermi-Bose crossover is shown in Figure 4(c) and 4(d), at two
different interactions $U=4t$ and $U=7t$, respectively. This behavior of the
DOS exhibited in Fig. 4(c) should be read in conjunction with Fig. 5 which
depicts the real space signatures of the system across the selected $U-V$
cross sections in terms of the pairing field amplitude ($|\Delta_{i}|$) and
pairing field phase coherence between a reference site and all other sites of
the lattice ($\cos(\theta_{0}-\theta_{i})$). For the sake of completeness the
spatial maps pertaining to the bosonic auxiliary field $\phi_{i}$ coupled to
the charge channel is also presented. At $U=4t$ the system is a gapped Fermi
superconductor upto $V_{g}\sim 1.10t$, as attested by the robust spectral gap
at the Fermi level with inverse square root singularities at the gap edges, as
shown in Figure 4(c). In the gapped Fermi superconductor regime, (as typified
for the case $U=4t$, $V=0.5t$ in Fig. 5), for all disorder strengths below the
scale $V_{Bose}$, the system corresponds to a (quasi) long range phase cohered
uniform Fermi superconductor characterized by a robust and uniform pairing
field amplitude.

Over the regime $1.10t\lesssim V<V_{Bose}\sim 1.25t$ there is finite spectral
weight at the Fermi level accompanied by the suppression of the coherence
peaks at the gap edges, characteristic of the gapless Fermi superconductor.
The real space signature of this gapless Fermi superconductor (at the
prototypical value of $U=4t$, $V=1.15t$), indicate that the corresponding
phase coherence is still (quasi) long ranged but the pairing field amplitude
is spatially fragmented into ”correlated” islands of suppressed
$|\Delta_{i}|$. At $V_{Bose}\sim 1.3t$, discontinuity appears at the gap
edges, indicating that the underlying phase is now a Bose superconductor. The
gapless superconductor is characterized by a finite spectral weight at the
Fermi level, and by a small but finite value for the superconducting structure
factor $S(0,0)$, shown in Fig. 4(b). This phase is characterized by a real
space picture, (see Fig. 5 for $U=4t$, $V=1.65t$) wherein one sees a
pronounced suppression of $|\Delta_{i}|$. Furthermore, the loss of phase
coherence in the gapless Bosonic superconductor is also affirmed by the
spatial maps for the phase correlations shown in Fig. 5 (for $U=4t$,
$V=1.65t$).

For $V\gtrsim V_{sc}$ a correlated Anderson insulator obtains characterized by
a gapless spectra, accompanied by large transfer of spectral weight away from
the Fermi level and absence of any phase coherence peak at the gap edges. The
gapless Bose superconductor and the correlated Anderson insulator are
distinguished from each other by the presence of (quasi) long range
superconducting phase coherence, in the former and its absence in the later
phase. The real space signature of the correlated Anderson Insulator, (
typified by the case at $U=4t,V=2.25t$ in Fig. 5) show the lack
superconducting correlations both in terms of the pairing field amplitude and
phase, as observed via the corresponding spatial maps.

Now, at strong coupling for instance at $U=7t$, the system is a gapped Fermi
superconductor up to $V_{Bose}\lesssim 0.65t$, shows spatial signatures that
corresponds to a robust (quasi)-long range phase cohered uniform Fermi
superconductor characterized by a spatially homogenous pairing field amplitude
as seen for the case $U=7t$, $V=0.25t$, akin to the weak coupling scenario
discussed above. At strong coupling, in the regime $V_{Bose}\leq V<V_{sc}\sim
1.10t$, the system is a gapped Bose superconductor: The DOS in this regime,
(see Fig. 4(d)), clearly exhibits a well developed gap at the Fermi level, and
the appearance of discontinuity at the gap edges. This gapped Bose
superconductor (exemplified by $U=7t,V=0.75t$ in Fig. 5), is characterized
with a real space map that shows a robust $|\Delta_{i}|$ and (quasi) long
range phase coherence. Over the regime $V_{sc}\lesssim V<V_{g}\sim 2t$, the
system is a spin-selective Bose insulator or metal with a robust spectral gap
at the Fermi level and absence of superconducting phase coherence.

Fig. 5 (for the case $U=7t,V=1.75t$), illustrates the prototypical behavior of
the gap and the phase correlations in such a spin-selective Bose
metal/insulator: this state is phase uncorrelated but hosts fairly large
amplitude $|\Delta_{i}|$, which suggests that the system is characterized by a
random distribution of the pairing field with a suppressed magnitude, i. e.
$|\Delta_{i}|$ are large but randomly distributed. The strong coupling regime,
for $V>V_{\rm g}$ as evidenced for the typical case ($U=7t$, $V=2.50t$) is a
correlated Anderson insulator. Pairing field amplitudes are significantly
suppressed in this regime akin to that of the weak coupling, represented by
$U=4t,V=2.25t$.

The fermionic number density remains pinned at the half filling ($n=1$) across
the $U-V$ space as attested by the largely uniform distribution of $\phi_{i}$,
barring the strong disorder regime where weak local fluctuations are observed
in the spatial distribution of $\phi_{i}$. Real space probes for instance, STM
would be able to discern the changes in the DOS as we tune through the Bose-
Fermi cross-over Lin  _et al._ (2023).

Figure 6: Spin-selective behavior of the DOS at $V=1.35t$ as a function of the
interaction strength: In the strong coupling regime $N_{\downarrow}(\omega)$
(see Fig. 6(a)), undergoes large transfer of spectral weight away from the
Fermi level indicating a Bose insulator phase. In comparison, the
$N_{\uparrow}(\omega)$ as seen in Fig. 6(b) shows prominent gap edge
discontinuities. (c)-(d) Disorder dependence of the ratio between the spin
selective energy gap at the Fermi level ($E_{g}^{\sigma}$) and the global
chemical potential ($\mu$) adjusted to the Fermi energy, as obtained from our
numerical calculations at $U=4t$ and $U=7t$ respectively. At weak coupling
both the spin species behaves identically, however the spin selective
dependence on the disorder potential becomes prominent at strong coupling. The
shaded region in panel (d) corresponds to the regime where only one of the
spin species is rendered gapless by disorder.

Spin selectivity: In Fig. 6(a) and Fig. 6(b) we contrast the spin resolved
single particle DOS for different interactions at a selected disorder
potential of $V=1.35t$. At this choice of the disorder potential the system is
a gapless Fermi superconductor for $U=3t$ and a gapless Bose superconductor at
$U=4t$ and $U=5t$ (see Fig. 4(b) for details). Over the regime $U\gtrsim 6t$ a
spin selective Bose metal/insulator state is realized. We understand these
phases based on the $N_{\sigma}(\omega)$ as follows: at $U=3t$ there is finite
spectral weight at the Fermi level in the single particle DOS, along with
prominent phase coherence peaks at the gap edges, as is expected from a
gapless Fermi superconductor. The spectrum continues to be gapless upto $U=5t$
following which the spectrum contains broader peaks with weak signatures of
discontinuity, indicating a cross-over into a gapless Bose superconducting
state. The spectral gap opens up at $U\geq 6t$ and the discontinuity at the
gap edges get progressively pronounced spin-selectively: with increasing
interaction for the $\uparrow$-spin species the discontinuities at the gap
edges become more pronounced. In contrast, for the $\downarrow$-spin species
there is progressively larger transfer of weight away from the Fermi level
with the corresponding gap edge discontinuities becoming more rounded. The
spin dependent dichotomy in the behavior of the DOS can be easily ascribed to
the fact that the $\downarrow$-spin species sees the disorder directly thus
rounding the discontinuities at the gap edges very strongly, (refer Fig.
6(a)). In contrast, the $\uparrow$-spin species sees a renormalized disorder
via residual interactions and thus sees a weaker disorder, as a consequence of
which the discontinuities at the gap edges in the $N_{\uparrow}(\omega)$ are
more prominent, (see Fig. 6(b)). This spin-selectiveness is once again brought
to the fore in Fig. 6(c) and Fig. 6(d) which depicts the disorder dependence
of the ratio of the spin-selective gap to the global chemical potential (set
to the Fermi energy), for two prototypical values at, $U=4t$, and $U=7t$. In
the weak coupling regime, ($U=4t$), as seen in Fig. 6(c) the ratio of the
superconducting gap to the global chemical potential $E_{\rm
g}^{\sigma}/|\mu|$, behaves identically for both the spin species. However, in
the strong coupling regime for instance, at the prototypical value $U=7t$, the
same ratio of $E_{\rm g}^{\sigma}/|\mu|$ shows asymmetric renormalization as a
function of the disorder strength, (see Fig. 6(d)): The superconducting gap
completely collapses for the $\downarrow$-spin fermions as a function of the
disorder strength whilst the superconducting gap is still robust for the
$\uparrow$-spin species.

We now turn our attention to the nature of the non superconducting state that
obtains at large values of $U$ above the scale set by $V_{\rm sc}$: the
anomalous spectrum of the $\uparrow$-spin species in this parameter range
characterized by $(i)$ a spectral gap at the Fermi level, $(ii)$ pronounced
discontinuity at the gap edges and $(iii)$ phase uncorrelated puddles of
superconducting pairing field with large amplitude (shown in Fig.5) Phillips
and Dalidovich (2003); Tsen  _et al._ (2016) leads us to speculate that for
$V\geq V_{\rm sc}$, for large values of interaction strength one might
stabilize a spin-selective Bose metal or insulator. The $\downarrow$-spin
species in contrast loses the bosonic correlations faster owing to the direct
influence of the disorder potential. We note that the exact nature of this
spin selective Bose metal/insulator regime can be ascertained only via the
computation of the relevant transport signatures, a task which is beyond the
scope of this work.

Discussion and conclusions: In conclusion, we have demonstrated that suitably
engineered random disorder potential can serve as a tuning parameter for the
Fermi-Bose crossover in 2D superconductor. Based on the thermodynamic and
spectroscopic signatures as determined via a non perturbative numerical
approach to the problem we have established the Fermi surface reconstruction
as the system undergoes disorder tuned Fermi-Bose crossover. Potential
candidates in solid state systems where such Fermi-Bose cross-over can be
realized are the the iron chalcogenides Shimojima  _et al._ (2017); Miao  _et
al._ (2015). These materials have the ratio $E_{g}/\epsilon_{F}\sim 1$ which
could be renormalized across the cross-over regime by chemical doping: for
example in FeSe1-xSx Hashimoto  _et al._ (2020), it has been observed that
substituting Se by the isovalent S in FeSe results in significant suppression
of the $E_{g}/\epsilon_{F}$ ratio, without altering the carrier concentration
in the material. This suppression of the $E_{g}/\epsilon_{F}$ ratio is very
qualitatively reminiscent of the suppression seen for example in Fig. 6(c),
and Fig. 6(d). However, the single orbital model used in our work does not
take into account the multiorbital band structure of real FeSe
superconductors. It is expected that a multiorbital approach to the problem
will bring in complexities arising out of the interplay between the inter and
intra-band superconducting pairing. It has been recently demonstrated that for
a multiorbital model the pseudogap phase is orbital selective and gives rise
to rich BCS-BEC cross-over physics Tajima  _et al._ (2019). Our work opens up
a tantalizing possibility that the BCS-BEC cross-over physics in iron
chalcogenides could plausibly occur when different orbitals see different
disorder: a sort of orbital selective counterpart of the spin-selective
disorder seen in this manuscript.

We next turn our attention to some open problems with regards to the effect of
spin-selective disorder in superconducting systems: for instance, an open
problem is to calculate the transport signatures so as to completely
characterize the nature of the normal state (i. e. Bose metal/insulator
regime). In a similar vein, the impact of thermal fluctuations in altering the
ground state properties is a problem that is being currently pursued MKR .
This is particularly relevant as thermal fluctuations induced spatial
fragmentation of the superconducting state into isolated phase uncorrelated
islands with large superconducting pairing amplitude could give rise to the
pseudogap phase at $T\neq 0$. Recently it has been argued that the finite
temperature pseudogap phase is in fact a Bose metal Maciej and Nandini (2020).
In a similar spirit the disorder induced spatial fragmentation of the
superconducting state can be thought of as the origin of the low temperature
Bose metal phase.

Acknowledgements: The authors acknowledge the use of high performance
computing cluster Aqua, at IIT, Madras, India. R.N. acknowledge funding from
the Center for Quantum Information Theory in Matter and Spacetime, IIT Madras
and from the Department of Science and Technology, Govt. of India, under Grant
No. DST/ICPS/QuST/Theme-3/2019/Q69. We thank Prof. G. Baskaran (IIT-Madras)
for valuable discussions. We also thank Dr. Prabuddha Chakraborty (ISI-
Bangalore), for a collaboration in the early part of this project.

## References

  * Randeria and Taylor (2014) M. Randeria and E. Taylor, Annual Review of Condensed Matter Physics 5, 209 (2014).
  * Emery and Kivelson (1994) V. J. Emery and S. A. Kivelson, Nature 374, 434 (1994).
  * (3) J. Corson, R. Mallozzi, J. Orenstein, J. N. Eckstein, and I. Bozovic, Nature 1999, 221.
  * Li  _et al._ (2010) L. Li, Y. Wang, S. Komiya, S. Ono, Y. Ando, G. D. Gu, and N. P. Ong, Phys. Rev. B 81, 054510 (2010).
  * (5) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, .
  * W. Ketterle (2008) M. W. Z. W. Ketterle, in  _Proceedings of the International School of Physics ”Enrico Fermi”, Course CLXIV, Varenna, 20â30 June 2006_ (M. Inguscio, W. Ketterle, C. Salomon, Eds. (IOS Press), 2008).
  * (7) J. T. Stewart, J. P. Gaebler, and D. S. Jin, Nature 2008, 744.
  * (8) J. P. Gaebler, J. T. Stewart, T. E. Drake, D. S. Jin, A. Perali, P. Pieri, and G. C. Strinati, Nature Physics 2010, 569.
  * Chin  _et al._ (2004) C. Chin, M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, J. H. Denschlag, and R. Grimm, Science 305, 1128 (2004).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Giorgini  _et al._ (2008) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
  * Kasahara  _et al._ (2016) S. Kasahara, T. Yamashita, A. Shi, R. Kobayashi, Y. Shimoyama, T. Watashige, K. Ishida, T. Terashima, T. Wolf, F. Hardy, C. Meingast, H. v. Lohneysen, A. Levchenko, T. Shibauchi, and Y. Matsuda, Nature Communications 7, 12843 (2016).
  * Rinott  _et al._ (2017) S. Rinott, K. B. Chashka, A. Ribak, E. D. L. Rienks, A. Taleb-Ibrahimi, P. L. Fevre, F. Bertran, M. Randeria, and A. Kanigel, Science Advances 3, e1602372 (2017).
  * Hashimoto  _et al._ (2020) T. Hashimoto, Y. Ota, A. Tsuzuki, T. Nagashima, A. Fukushima, S. Kasahara, Y. Matsuda, K. Matsuura, Y. Mizukami, T. Shibauchi, S. Shin, and K. Okazaki, Science Advances 6, eabb9052 (2020).
  * Lubashevsky  _et al._ (2012) Y. Lubashevsky, E. Lahoud, K. Chashka, D. Podolsky, and A. Kanigel, Nature Physics 8, 309 (2012).
  * Crépel  _et al._ (2023) V. Crépel, D. Guerci, J. Cano, J. H. Pixley, and A. Millis, Phys. Rev. Lett. 131, 056001 (2023).
  * Loh  _et al._ (2016) Y. L. Loh, M. Randeria, N. Trivedi, C.-C. Chang, and R. Scalettar, Phys. Rev. X 6, 021029 (2016).
  * Hazra  _et al._ (2020) T. Hazra, N. Trivedi, and M. Randeria, (2020), arXiv:2011.06598 [cond-mat.supr-con] .
  * Sacépé  _et al._ (2011) B. Sacépé, T. Dubouchet, C. Chapelier, M. Sanquer, M. Ovadia, D. Shahar, M. Feigel’man, and L. Ioffe, Nature Physics 7, 239 (2011).
  * Sacépé  _et al._ (2010) B. Sacépé, C. Chapelier, T. I. Baturina, V. M. Vinokur, M. R. Baklanov, and M. Sanquer,  1, 140 (2010).
  * Fisher (1990) M. P. Fisher, Physical Review Letters 65, 923 (1990).
  * Vojta (2006) T. Vojta, J. Phys. A 39, R143 (2006).
  * Vojta (2013) T. Vojta, AIP Conference Proceedings 1550, 188 (2013).
  * Vojta (2019) T. Vojta, Ann. Rev. Condens. Mat. Phys. 10, 233 (2019).
  * Hoyos  _et al._ (2007) J. A. Hoyos, C. Kotabage, and T. Vojta, Phys. Rev. Lett. 99, 230601 (2007).
  * Vojta  _et al._ (2009) T. Vojta, C. Kotabage, and J. A. Hoyos, Phys. Rev. B 79, 024401 (2009).
  * Ghosal  _et al._ (2001) A. Ghosal, M. Randeria, and N. Trivedi, Phys. Rev. B 65, 014501 (2001).
  * Dubi  _et al._ (2007) Y. Dubi, Y. Meir, and Y. Avishai, Nature 449, 876 (2007).
  * Kapitulnik  _et al._ (2019) A. Kapitulnik, S. A. Kivelson, and B. Spivak, Rev. Mod. Phys. 91, 011002 (2019).
  * Kaur  _et al._ (2022) S. Kaur, H. K. Kundu, S. Kumar, A. Dogra, R. Narayanan, T. Vojta, and A. Bid, (2022), arXiv:2206.10215 [cond-mat.supr-con] .
  * Xing  _et al._ (2015) Y. Xing, H.-M. Zhang, H.-L. Fu, H. Liu, Y. Sun, J.-P. Peng, F. Wang, X. Lin, X.-C. Ma, Q.-K. Xue, _et al._ , Science 350, 542 (2015).
  * Shen  _et al._ (2016) S. Shen, Y. Xing, P. Wang, H. Liu, H.-L. Fu, Y. Zhang, L. He, X. Xie, X. C.and Lin, J. Nie, and J. Wang, Phys. Rev. B 94, 144517 (2016).
  * Lewellyn  _et al._ (2019) N. A. Lewellyn, I. M. Percher, J. Nelson, J. Garcia-Barriocanal, I. Volotsenko, A. Frydman, T. Vojta, and A. M. Goldman, Phys. Rev. B 99, 054515 (2019).
  * Reiss  _et al._ (2021) P. Reiss, D. Graf, A. A. Haghighirad, T. Vojta, and A. I. Coldea, Phys. Rev. Lett. 127, 246402 (2021).
  * Jiang  _et al._ (2013) M. Jiang, R. Nanguneri, N. Trivedi, G. G. Batrouni, and R. T. Scalettar, New Journal of Physics 15, 023023 (2013).
  * Nanguneri  _et al._ (2012) R. Nanguneri, M. Jiang, T. Cary, G. G. Batrouni, and R. T. Scalettar, Phys. Rev. B 85, 134506 (2012).
  * Karmakar  _et al._ (2022) M. Karmakar, S. Roy, S. Mukherjee, and R. Narayanan, Phys. Rev. Res. 4, 043159 (2022).
  * McKay and DeMarco (2010) D. McKay and B. DeMarco, New Journal of Physics 12, 055013 (2010).
  * Mandel  _et al._ (2003a) O. Mandel, M. Greiner, A. Widera, T. Rom, T. W. Hansch, and I. Bloch, Nature 425, 937 (2003a).
  * Feiguin and Fisher (2009) A. E. Feiguin and M. P. A. Fisher, Phys. Rev. Lett. 103, 025303 (2009).
  * Liu  _et al._ (2004) W. V. Liu, F. Wilczek, and P. Zoller, Phys. Rev. A 70, 033603 (2004).
  * Mandel  _et al._ (2003b) O. Mandel, M. Greiner, A. Widera, T. Rom, T. W. Hänsch, and I. Bloch, Phys. Rev. Lett. 91, 010407 (2003b).
  * Skolimowski  _et al._ (2015) J. Skolimowski, D. Vollhardt, and K. Byczuk, Phys. Rev. B 92, 094202 (2015).
  * Ohtomo and Hwang (2004) A. Ohtomo and H. Y. Hwang, Nature 427, 423 (2004).
  * Zhang  _et al._ (2019) M. Zhang, K. Du, T. Ren, H. Tian, Z. Zhang, H. Y. Hwang, and Y. Xie, Nature Communications 10, 4026 (2019).
  * Han  _et al._ (2019) K. Han, K. Hu, X. Li, K. Huang, Z. Huang, S. Zeng, D. Qi, C. Ye, J. Yang, H. Xu, A. Ariando, J. Yi, W. Lae, S. Yan, and X. R. Wang, Science Advances 5, eaaw7286 (2019).
  * (47) P. Xu, W. Han, P. M. Rice, J. Jeong, M. G. Samant, K. Mohseni, H. L. Meyerheim, S. Ostanin, I. V. Maznichenko, I. Mertig, E. K. U. Gross, A. Ernst, and S. S. P. Parkin, Advanced Materials 29, 1604447.
  * Niu  _et al._ (2017) W. Niu, K. Du, S. Wang, M. Zhang, M. Gao, Y. Chen, H. Liu, W. Zhou, F. Song, P. Wang, Y. Xu, X. Wang, J. Shen, and R. Zhang, Nanoscale 9, 34 (2017).
  * Chen  _et al._ (2011) Y. Chen, N. Pryds, J. E. Kleibeuker, G. Koster, J. Sun, E. Stamate, B. Shen, G. Rijnders, and S. Linderoth, Nano Letters 11, 3774 (2011).
  * Li  _et al._ (2011) L. Li, C. Richter, J. Mannhart, and R. C. Ashoori, Nature Physics 7, 762 (2011).
  * Dikin  _et al._ (2011) D. A. Dikin, M. Mehta, C. W. Bark, C. M. Folkman, C. B. Eom, and V. Chandrasekhar, Phys. Rev. Lett. 107, 056802 (2011).
  * Bert  _et al._ (2011) J. A. Bert, B. Kalisky, C. Bell, M. Kim, Y. Hikita, H. Y. Hwang, and K. A. Moler, Nature Physics 7, 767 (2011).
  * Bi  _et al._ (2014) F. Bi, M. Huang, S. Ryu, H. Lee, C.-W. Bark, C.-B. Eom, P. Irvin, and J. Levy, Nature Communications 5, 5019 (2014).
  * Mena  _et al._ (2006) F. P. Mena, J. F. DiTusa, D. van der Marel, G. Aeppli, D. P. Young, A. Damascelli, and J. A. Mydosh, Phys. Rev. B 73, 085205 (2006).
  * Hubbard (1959) J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
  * Schulz (1990) H. J. Schulz, Phys. Rev. Lett. 65, 2462 (1990).
  * Karmakar and Majumdar (2016) M. Karmakar and P. Majumdar, Phys. Rev. A 93, 053609 (2016).
  * Evenson  _et al._ (1970) W. E. Evenson, J. R. Schrieffer, and S. Q. Wang, Journal of Applied Physics 41, 1199 (1970).
  * Abrikosov and Gor’kov (1961) A. A. Abrikosov and L. P. Gor’kov, JETP 12, 1243 (1961).
  * Lin  _et al._ (2023) H. Lin, W. Huang, G. Rai, Y. Yin, L. He, Q.-K. Xue, S. Haas, S. Kettemann, X. Chen, and S.-H. Ji, Phys. Rev. B 107, 104517 (2023).
  * Phillips and Dalidovich (2003) P. Phillips and D. Dalidovich, Science 302, 243 (2003).
  * Tsen  _et al._ (2016) A. W. Tsen, B. Hunt, Y. D. Kim, Z. J. Yuan, S. Jia, R. J. Cava, J. Hone, P. Kim, C. R. Dean, and A. N. Pasupathy, Nature Physics 12, 208 (2016).
  * Shimojima  _et al._ (2017) T. Shimojima, W. Malaeb, A. Nakamura, T. Kondo, K. Kihou, C.-H. Lee, A. Iyo, H. Eisaki, S. Ishida, M. Nakajima, S. ichi Uchida, K. Ohgushi, K. Ishizaka, and S. Shin, Science Advances 3, e1700466 (2017).
  * Miao  _et al._ (2015) H. Miao, T. Qian, X. Shi, P. Richard, T. K. Kim, M. Hoesch, L. Y. Xing, X. C. Wang, C. Q. Jin, J. P. Hu, and H. Ding, Nature Communications 6, 6056 (2015).
  * Tajima  _et al._ (2019) H. Tajima, Y. Yerin, A. Perali, and P. Pieri, Phys. Rev. B 99, 180503 (2019).
  * (66) Madhuparna Karmakar, and Rajesh Narayanan, Unpublished.
  * Maciej and Nandini (2020) M. Maciej and T. Nandini, Phys. Rev. B 102, 144506 (2020).


$\dagger$$\dagger$footnotetext: To whom correspondence should be addressed.
Email<EMAIL_ADDRESS>or<EMAIL_ADDRESS>
# Direct visualization of an elusive molecular reaction: Time-resolved H2
roaming in acetonitrile

Debadarshini Mishra*† Department of Physics, University of Connecticut,
Storrs, Connecticut, 06269, USA  These authors contributed equally. Aaron C.
LaForge*† Department of Physics, University of Connecticut, Storrs,
Connecticut, 06269, USA  Lauren M. Gorman  Department of Physics, University
of Connecticut, Storrs, Connecticut, 06269, USA  Sergio Díaz-Tendero
Departamento de Química, Módulo 13, Universidad Autónoma de Madrid, 28049
Madrid, Spain, EU  Condensed Matter Physics Center (IFIMAC), Universidad
Autónoma de Madrid, 28049 Madrid, Spain, EU  Institute for Advanced Research
in Chemical Sciences (IAdChem), Universidad Autónoma de Madrid, 28049 Madrid,
Spain  Fernando Martín  Departamento de Química, Módulo 13, Universidad
Autónoma de Madrid, 28049 Madrid, Spain, EU  Condensed Matter Physics Center
(IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain, EU  Instituto
Madrileño de Estudios Avanzados en Nanociencia (IMDEA-Nano), Campus de
Cantoblanco, 28049 Madrid, Spain, EU  Nora Berrah  Department of Physics,
University of Connecticut, Storrs, Connecticut, 06269, USA

(August 27, 2024)

###### Abstract

We introduce a novel approach to imaging the complex phenomenon of roaming, an
unconventional mechanism that contributes to molecular dissociation. During
roaming, fragments remain weakly-bound to each other due to long-range
electrostatic interactions. This extended period of interaction time allows
for the formation of new molecular compounds. However, the neutral character
of the roaming fragment and its indeterminate trajectory have made it
difficult to experimentally identify and systematically study. To address this
challenge, we utilize intense, femtosecond IR radiation in combination with
coincident Coulomb explosion imaging to directly reconstruct the momentum
vector of the neutral roaming H2 from ionized acetonitrile, CH3CN. This
technique not only provides a kinematically complete picture of the underlying
molecular dynamics that contribute to H3+ formation, but more significantly,
yields an unambiguous experimental signature of roaming. We corroborate this
with the aid of quantum chemistry calculations, which allow us to fully
determine how this unique dissociative process occurs. Our method enables a
direct visualization of roaming reactions, making the investigation of this
complex process more feasible.

## Introduction

Understanding the multitude of complex molecular processes proceeding through
photoexcitation is of fundamental importance to a wide variety of biological
and chemical systems. Dissociation, one of the more common photochemical
processes, typically proceeds by the fragmenting molecule moving along the
path of minimum energy. However, there exists a more indirect pathway, known
as roaming, where the neutral fragment explores the relatively flat regions of
the potential energy surface far from the minimum energy path due to long-
range, weakly-bound interactions. As the name suggests, the indirect nature of
dissociation through roaming makes this channel more complex than the more
typical fragmentation channels. That said, roaming has been a subject of great
interest in recent years [1, 2, 3]. From its initial indirect observation in
formaldehyde [4], roaming has been proposed to occur in several small
molecules such as acetaldehyde [5, 6], acetone [7], nitrate [8], methyl
formate [9, 10], and propane [11]. Recently, roaming of H2 in monohydric
alcohols, initiated by intense femtosecond laser pulses, was shown to lead to
H3+ formation [12, 13, 14]. In general, photoexcitation of organic molecules,
using high-powered, ultrafast lasers, has proven to be an efficient means to
produce H3+ [15, 16, 17, 18, 19, 20, 21]. Additionally, H3+, an abundant
cation in interstellar media [22, 23], is of great interest due to its
importance in gas-phase chemistry where it plays a critical role in the
production of more complex molecules [24].

Figure 1: Schematic of H2 roaming and H${}_{3}^{+}$ formation in CH3CN. The
pump pulse creates $[\ce{CH3CN}]^{2+}$ by multiphoton absorption followed by
neutral fragmentation of H2. Some H2 remains weakly-bound to the dication in a
roaming state. Afterward, an H+ is abstracted from the dication forming H3+.
If the time-delayed probe pulse arrives before abstraction, the process
prevents H3+ formation.

Although roaming has been proposed to occur in a multitude of small molecules,
a direct visualization of roaming has remained elusive from experimental
observations due to its complex nature. The central problem is the random
nature of the roaming process itself, where the relative position of the
roaming particle with respect to the rest of the system is not deterministic.
For instance, pioneering experimental studies [12, 13, 14] inferred the role
of roaming neutral H2 in the H3+ formation process by relying heavily on
simulations. These studies could not fully exploit the potential of coincident
momentum imaging, which is crucial in distinguishing a true experimental
signature for roaming from other fragmentation processes. More recently, Endo
et al. [25] demonstrated time-resolved signature of roaming in formaldehyde by
investigating a three-ion coincidence channel using Coulomb explosion imaging
(CEI). However, this channel carries imprints of various reactions that are
triggered via UV excitation of formaldehyde, including molecular and radical
dissociation, as well as roaming. As a result, it was necessary to incorporate
high-level simulations to identify the signature of roaming from the
experimental results.

In this paper, we conduct a systematic, time-resolved investigation of the
formation of H3+ in acetonitrile, CH3CN, by directly tracking the dynamics of
the roaming neutral H2, which is ‘invisible’ to charged-particle detection
schemes. We make use of coincident momentum imaging in combination with
femtosecond IR-IR pump-probe spectroscopy. This allows us to access the full,
time-resolved 3D momentum information of each detected fragment, while also
enabling us to reconstruct the momentum vector of the undetected neutral or
‘invisible’ fragment in the reaction, by virtue of momentum conservation. Such
kinematically complete information, even for a coincidence channel involving a
neutral fragment, is crucial for mapping out a direct signature of roaming
reactions, without any interference from other excited state dynamics.
Furthermore, the choice of acetonitrile for our current study makes the
analysis of H3+ formation easier, since it only has three hydrogen atoms that
are bonded to the same carbon atom, thereby avoiding complications that may
arise with molecules that have more than three hydrogen atoms. A schematic of
H2 roaming and H3+ formation in CH3CN is given in Fig. 1. We produce the
parent dicationic state with an intense, femtosecond pump pulse. H3+ is formed
as a result of a roaming H2 neutral which remains in the flat part of the
potential energy surface before abstracting a hydrogen ion from the parent
molecule. The formation of H3+ is probed with a time-delayed fs pulse that
disrupts the ensuing dynamics by ionizing the neutral H2, thereby preventing
the formation of H3+ and depleting its yield. A more in-depth discussion of
our specific experimental setup is given in the Methods section. The momenta
of the resulting ionic fragments from the laser-molecule interaction are
measured in coincidence using the COLTRIMS technique [26, 27].

## Results and Discussion

Figure 2: (a) Time-dependent kinetic energy release of H3+ + C2N+. (b)
Projection of KER signal for H3+ + C2N+ (light green) and D3+ + C2N+ (dark
green) on the pump-probe delay axis. The projected signal intensity in (b) is
normalized to 1 at 0 fs pump-probe delay for both H3+ and D3+ channels.

To elucidate the role of H2 roaming in the formation of H3+ in photoionized
acetonitrile, we begin by investigating its kinematically complete time-
resolved dynamics. Fig. 2 (a) shows the complete KER map of the channel H3+ +
C2N+ as a function of pump-probe delay while (b) shows the projection of the
KER band between 3 eV and 7 eV. We primarily observe a horizontal band
centered around a KER of 5 eV, due to a single pulse exciting the molecule
directly into the dicationic state. Typically, these types of KER maps show a
dynamic, energy-dependent behavior with respect to pump-probe delay. In that
case, the pump pulse excites the system into the ionic state and the time-
delayed probe pulse further ionizes the system into a higher charged state,
resulting in a mapping of the kinetic energy from Coulomb repulsion as a
function of time. In contrast, here we only observe a very weak energy-
dependent behavior in the KER map, since the formation of H3+ is induced
directly through the dicationic state. That said, within the horizontal band,
there is a strong dependence on the pump-probe delay. From the projection,
shown in Fig. 2 (b), we see an enhancement in the total yield of this channel
around time zero of the pump-probe delay due to the overlap of the two pulses.
However, on either side of time zero, we see an immediate decrease in the
yield before the signal steadily increases again. As discussed above, this
behavior is the result of the disruptive nature of the probe pulse that
prevents H3+ formation by ionizing the neutral roamer H2 before it can
abstract an H+. In other words, if the probe pulse arrives early enough, then
it ‘disrupts’ the roaming H2, resulting in a decreased H3+ yield. For longer
time delays, the probe pulse is less likely to disrupt H2 roaming; thus, the
yield of H3+ increases until reaching saturation.

We have additionally performed similar measurements using deuterated
acetonitrile, which allows us to better resolve different masses when using
ion time-of-flight spectroscopy. Overall, we observed no significant
differences in the time-resolved measurements for either sample. This, in
itself, is a rather interesting result since it reveals that the roaming
process relies mostly on electronic properties rather than on nuclear
dynamics. This is clearly illustrated in Fig. 2 (b), where we show the KER
projections of the formation of both H3+ (light green line) and D3+ (dark
green line). The timescales of depletion and enhancement of the H3+ and D3+
projections are almost identical, thus showing that the process is independent
of the mass of the roamer. Additionally, we have fitted the projections with a
rising exponential of the form: $y=A(1-e^{(t-t_{0})/\tau})+y_{0}$ where $A$ is
the normalization constant, $t$ is the pump-probe delay, $t_{0}$ is the time
offset, $\tau$ is the time constant, and $y_{0}$ is the intensity offset. $A$
and $y_{0}$ are arbitrary constants since the intensity is normalized to 1 at
0 fs pump-probe delay. However, $t_{0}$ and $\tau$ have physical significance
since they reveal the timescales of the respective depletion and enhancement
of the H3+ and D3+ formation. The values from the fitting function are given
in Table S1 in the supplementary information (SI) while a comparison of the
fits with the experimental data is given in Fig. S2 of the SI. The time
offset, $t_{0}$, has an average value of 125 fs, which is larger than the
temporal resolution of the experiment (60 fs) and the time constant, $\tau$
has an average value of 105 fs. Overall, the fits reveal that H2 roaming and
H3+ formation occur on ultrafast timescales of a few hundred femtoseconds.

Figure 3: (a) Detailed analysis of a trajectory leading to the formation of
H3+. H3+ is formed as a result of roaming neutral H2 followed by a proton
abstraction. (b) Time-dependent change of the three C-H bond distances for the
same trajectory used in (a), showing H2 roaming and H abstraction leading to
formation of H3+. (c) Time-dependent change in the angle between NCβ and the
three H atoms obtained from the same trajectory used in (a). In this
trajectory, the three angles coincide around 400 fs implying ejection of H3+.
Acetonitrile molecule is shown in the inset identifying the labels for all the
carbon and hydrogen atoms. Figure 4: Photoion-photoion coincidence lines for
(a) D3+ + C2N+ and (b) D+ + C2N+ + D2 from acetonitrile, integrated over all
pump-probe delays. Black dashed lines of slope -1 are overlaid on the PIPICO
channels (a) and (b) to demonstrate the agreement between the expected and
experimentally-obtained slope of these channels. (c) Pump-probe delay-
dependent kinetic energy projection of the roaming D2 neutral. (d) Angular
distribution between the momentum vectors of D+ and D2 (blue line), C2N+ and
D2 (orange line), D+ and C2N+ (green line), and D3+ and C2N+ (crimson line)
within the first 200 fs time-delay window. (e) Newton diagram for the channel
D+ + D2 neutral + C2N+ integrated over the first 200 fs of pump-probe delay.
The momentum vector of C2N+ lies along the x-axis while those of D+ and D2
neutral are plotted in the top and bottom halves, respectively.

To better understand the mechanism of H2 roaming and H3+ formation, we have
performed molecular dynamics simulations, described in detail in the Methods
section. We have found that the primary fragmentation channels producing
cationic hydrogen species and their respective abundances are: ${\rm
H^{+}/H_{2}C_{2}N^{+}}$ \- 79.6%, ${\rm H_{2}^{+}/HC_{2}N^{+}}$ \- 7.6%, and
${\rm H_{3}^{+}/C_{2}N^{+}}$ \- 0.8%. For the current work, we focus the
analysis on those trajectories producing ${\rm H_{3}^{+}}$. Figure 3 shows the
results from one of such trajectories leading to H3+. Overall, several
processes occur within the first few hundred femtoseconds which impact the
${\rm H_{3}^{+}}$ formation. H2 is produced in the first tens of fs and then
roams around up to $\sim$ 150 fs. Subsequently, H exchange occurs wherein the
H atoms of the roaming H2 molecule are bonded back to the carbon atom,
immediately followed by a different pair of H atoms being ejected to form a
new roaming H2 molecule. After another $\sim$ 200 fs, this H2 molecule
abstracts the last remaining H bonded to the C atom, thereby forming H3+ in
coincidence with C2N+. The plots in Fig. 3 (b) and (c) show time-evolution of
the bond distances and angles of the three H atoms relative to the molecular
skeleton of CCN for a typical trajectory. Results for the other trajectories
are given in Figs. S3 and S4, and the corresponding variations of the nuclear
kinetic energy with respect to time are given in Fig. S8. In all cases, H3+ is
formed via neutral H2 roaming, however, due to the small number of
trajectories leading to H3+, we cannot provide an exact value for the roaming
time. That said, all simulated trajectories in the present work predict a
formation timescale between 100 fs and 400 fs for H3+, which overall agrees
very well with the experimental measurements. The calculated nuclear kinetic
energies in the final state are also in good agreement with the experimental
KER (see Fig. S8).

In order to gain a detailed understanding of the roaming mechanism in
acetonitrile, we track the dynamics of the neutral D2 which appears in
coincidence with D+ and C2N+ using photoion–photoion coincidence (PIPICO)
maps. In general, negative-sloping distributions observed in the PIPICO
coincidence maps identify fragments created by multiple ionization.
Specifically, the shape and slope of the distribution give information about
the dissociation process [28]. For instance, the back-to-back emission of two
ions, where momentum is conserved, results in a sharp line of slope -1. For
deuterated acetonitrile, this is observed in Fig. 4 (a) for D3+ in coincidence
with C2N+. Prior to D3+ formation, the PIPICO channel that allows us to track
the behavior of the neutral roamer is the incomplete channel of D+ in
coincidence with C2N+, shown in Fig. 4 (b), which is missing the mass from D2.
Incomplete coincidence channels can have unique shapes and slopes depending on
the missing fragment [28]. When there is deferred charge separation of the
different fragments with mass $m$ such that a neutral $m_{3}$ is ejected
first, followed by $m_{1}m_{2}^{2+}$ breaking apart into two ions, the
resulting coincidence channel has a broad, lozenge shape of slope -1. This is
precisely what we observe in Fig. 4 (b) for the channel D+ + C2N+ + D2. To
show the agreement between the calculated slopes and experimental data, black
dashed lines of slope -1 are overlaid in the PIPICO channels in Figs. 4 (a)
and (b), respectively. To further prove that the missing fragment in this
channel is indeed neutral and not an ion that was simply lost in the detection
process, we have plotted the kinetic energy of the missing D2 in Fig. S5 of
the SI. If it is an ion, one would observe a time-dependent decrease in its
kinetic energy as a result of Coulomb explosion from the remaining charged
moiety. On the other hand, a neutral dissociating from a charged fragment
would result in very low, time-independent kinetic energy. The latter is
observed in the reconstructed kinetic energy of D2, thus proving that the
missing fragment is indeed neutral. Fig. 4 (c) shows the reconstructed kinetic
energy projection of neutral D2 as a function of pump-probe delay where the
time-dependent behavior of this neutral fragment is similar to that of H3+ +
C2N+, shown in Fig. 2 (b).

After isolating the coincidence channel involving the target neutral roamer
D2, we then study its angular correlations with respect to the coincident
ionic fragments. In Fig. 4 (d), we plot the angle between the momentum vectors
of D2 and D+ (blue line), and D2 and C2N+ (orange line). The momentum vector
of D2, which is invisible to our detection scheme, is reconstructed by
implementing momentum conservation among the three fragments. Overall, we see
a very broad distribution implying a highly uncorrelated emission angle
between the neutral D2 and the two ionic fragments, which is a key signature
of roaming. Unlike direct dissociation where the fragments follow a fixed
minimum energy path, roaming fragments traverse relatively flat regions of the
potential energy surface, which results in wide variations in their pathways.
Such an uncorrelated behavior is exactly what makes roaming difficult to
observe since the position of the roaming neutral is not well-defined as
compared to a typical dissociation channel. In that regard, we have a
potentially more direct method for visualizing roaming by directly mapping the
momentum of the reconstructed neutral. For comparison, the ion-ion angular
distribution for D+ and C2N+ (green line) in this incomplete channel is much
sharper and centered at $\sim$ 160∘, indicating strong Coulomb interaction
resulting in a stronger momentum correlation. Such an angular distribution
between D+ and C2N+ is also observed from the calculated roaming trajectories.
For the trajectory shown in Fig. 3 (c), the angle between D+ and C2N+ in the
first 200 fs is approximately 150∘ with a dispersion around this value of the
order of 30∘. Finally, the angular distribution for D3+ and C2N+ (crimson
line), which forms a complete two-body channel, is even narrower and centered
at $\sim$ 175∘ indicating back-to-back emission, as would be expected from a
two-body fragmentation process.

Beyond examining simple angular correlations, we can also construct the full
3D momentum of the reaction, using Newton diagrams. In general, Newton
diagrams are very useful tools for providing structural information about the
reaction dynamics of a system [29, 30, 31]. One of their unique features is
the ability to distinguish different mechanisms of fragmentation. For
instance, concerted, instantaneous fragmentation would appear as sharp,
localized distributions in the Newton diagram [32], whereas sequential
fragmentation would appear as a semicircle [33, 31] due to the rotation of the
intermediate fragment whose lifetime is larger than its rotational period. A
typical Newton diagram for a three-body fragmentation channel in acetonitrile,
D+ + D2+ + C2N+, is shown in Fig. S6 of the SI. In this case, the momentum
distributions are sharp with strong angular correlations.

Figure 4 (e) shows the Newton diagram for the channel D+ + D2 neutral + C2N+.
Here, the pump-probe delay is constrained to values less than 200 fs. The
momentum vector of C2N+ is fixed along the x-axis, while those of D+ and the
reconstructed D2 are plotted on the top and bottom halves, respectively. Due
to momentum conservation amongst the three fragments, the Newton diagram lies
on a plane and shows the relative momentum distributions of two fragments with
respect to C2N+. Overall, there are several features that are markedly
different from the analogous all-ion channel (D+ + D2+ + C2N+), mentioned
above. The distribution of the neutral D2 in Fig. 4 (e) appears quite
different from a typical fragmentation channel since there is no clear angular
dependence in addition to a very broad, albeit lower momentum magnitude
distribution. This shows that the neutral D2 departs with a smaller fraction
of the total momentum involved in this process, as would be expected from
neutral fragmentation. It must be noted that, in sequential processes,
rotation of the intermediate fragment can also result in broad angular
distributions. Therefore, the Newton diagram in Fig. 4 (e) is integrated over
the first 200 fs pump-probe delay window, which is much smaller than the
rotational timescale for C2ND+. Thus, the broad angular correlation that we
observe here cannot be due to the rotation of the intermediate fragment C2ND+,
thereby removing any ambiguity between rotation and roaming. Such a lack of
clear angular correlation in conjunction with a broad distribution of momentum
magnitude of the neutral D2 is a potentially unambiguous signature for roaming
in molecules. Contrary to that, the two ions in this incomplete channel, D+
and C2N+, share a stronger momentum correlation mediated via Coulomb
interaction as evidenced by a narrower angular spread along with a localized
distribution in the magnitude of momentum.

Figure 5: (a) Minimum Energy Path (MEP) followed in the potential energy
surface of the doubly ionized acetonitrile; relative energy in eV referred to
the neutral molecule as a function of the Intrinsic Reaction Coordinate (IRC).
(b) and (c) Electrostatic potential (ESP) of $\rm HCCN^{2+}$ mapped on the
electronic density with isovalue 0.0004 a.u.. Color code for the extreme
values: red = 0.32 a.u. and dark blue = 0.48 a.u.. The geometry used in (b)
corresponds to the channel $\rm H_{2}/HCCN^{2+}$, i.e. after release of
neutral $\rm H_{2}$, and in (c) to the weakly bonded $\rm H_{2}...HCCN^{2+}$,
i.e. last point in the MEP.

To demonstrate how robust this method is for distinguishing neutral
fragmentation channels, in Fig. S1 of the SI, we compare our results to a
different type of neutral fragmentation, which was measured in
2-propanol(CH3CHOHCH3). In this case, the fragmentation channel, CH3+ + H2O
neutral + C2H3+, is initially two primary fragments, CH3+ and C2H5O+, the
latter of which subsequently dissociates into C2H3+ and neutral H2O. In
comparison to neutral roaming, we observe stark differences in the PIPICO
channel, angular distributions, and Netwon diagram, which are discussed in
detail in the SI.

Further insight into H2 roaming has been obtained in a careful exploration of
the potential energy surface (PES). In particular, we computed the minimum
energy path (MEP) that a doubly charged acetonitrile molecule would follow
after vertical ionization, i.e., the maximum gradient on the PES, starting
from the optimized geometry of the neutral system (Fig. 5(a)). We clearly
observe the release of neutral $\rm H_{2}$, followed by the production of a
weakly bonded $\rm H_{2}...HCCN^{2+}$ complex. Furthermore, this shows that
the neutral $\rm H_{2}$ is polarized by the potential exerted by the cationic
fragment. We evaluated such potential by computing the electrostatic potential
(ESP) of $\rm HCCN^{2+}$ in two different configurations: using the geometry
of the optimized structure from the moiety (Fig. 5(b)) and using the geometry
of the last point given in the MEP exploration (Fig. 5(c)). The ESP is plotted
in both cases by projecting the actual value on the electronic density
isosurface, which provides a visualization of the potential that polarizes the
neutral $\rm H_{2}$. Overall, Figs. 5 (b) and (c) reveal that the regions with
the highest electrostatic potential (blue in color) are nearest to the C–H
bond in $\rm HCCN^{2+}$. It is most likely that the roaming H2 remains weakly-
bound in this region. Upon fragmentation of H+ and CCN+, the H2 is
approximately at 90∘ with respect to the two ionic fragments, which
qualitatively agrees with the angular correlations given in Fig. 4(d).

The fact that acetonitrile has only 3 hydrogen atoms that can contribute to
D3+ formation makes it an ideal candidate to unambiguously track such roaming
reactions. However, in order to be exhaustive in terms of all the possible
ways that D3+ can be formed in acetonitrile, we must also consider the
possibility of a roaming neutral D which can abstract D2 from the remaining
moiety to form D3+. For this, we have to consider the incomplete coincidence
channel: D2+ + C2N+ + D. We note that the yield for this channel is nearly 4.5
times lower than D+ + C2N+ + D2, consequently making it a weak contributor to
the D3+ channel. This is not surprising, since a free H or D atom is a
radical, which makes it less likely to be formed than H2 or D2. In fact, this
channel is not observed in our theoretical simulations due to the limited
number of trajectories used in the present work. To our knowledge, it has also
not been observed in previous works. Following the same sequence as before, we
first ensure that the missing D is indeed a neutral by looking at its kinetic
energy as a function of time-delay. Furthermore, in Fig. S7, we show the
angular correlation between all pairs of fragments in the channel D2+ + D
neutral + C2N+. The broad and symmetric angular correlations between the
neutral D and the other two ionic fragments are indications of roaming D.
Unlike in Fig. 4, the correlation between D and C2N+ is quite symmetric. This
could be a result of D being a lighter roaming fragment compared to D2.

## Conclusion

In conclusion, we have developed a more direct means to image neutral roaming
reactions in small molecular systems. Using ultrafast IR pump-IR probe
spectroscopy combined with coincident Coulomb explosion imaging, we have time-
resolved the formation of H3+ from neutral H2 roaming in acetonitrile, which
was measured to occur within a few hundred femtoseconds. This technique
enables us to directly track the ‘invisible’ neutral roamer, the precursor to
H3+ formation, from incomplete fragmentation channels. Additionally, with the
aid of quantum chemistry calculations, we have fully simulated the possible
roaming trajectories, allowing us to follow some of the unique intramolecular
processes which occur in this type of dissociation. In general, our novel
technique gives us a more straightforward means to observe neutral fragments
which can allow us to gain a better understanding of the underlying molecular
dynamics in roaming reactions.

## Methods

The experimental setup has been previously described elsewhere [34, 35, 36].
We used a 5 kHz Ti:Sapphire laser producing 35 fs pulses with a central
wavelength of 800 nm. The laser beam was split into two independent arms which
were time-delayed with respect to one another. Each arm had a variable
intensity that was optimized to produce the highest contrast on the H3+/D3+
signal. For this experiment, the power was in the range of $10^{14}$ W/cm2.
The two beams were directed into the interaction region where they were
overlapped and back-focused to a spot size of about 10 $\mu$m with a temporal
resolution of about 60 fs. Both beams were linearly polarized in the direction
of the time-of-flight axis of the spectrometer. A cold, molecular jet of
acetonitrile, from a room-temperature bubbler, was produced by expansion
through a 30 $\mu$m nozzle seeded with 1 bar of helium gas and propagated into
the spectrometer perpendicular to the two laser pulses. Acetonitrile was
ionized and resulting charged fragments were detected by a Cold Target Recoil
Ion Momentum Spectrometer (COLTRIMS) [26, 27]. Using a weak, homogeneous
electric field, the ions were directed towards a position-sensitive detector,
which is capable of measuring the three-dimensional momentum distributions of
the charged particles. For CEI, we primarily focus on measuring the fragmented
ions, and for the current experimental conditions, we have an ion momentum
resolution of 0.1 a.u. By applying the coincidence technique, we can isolate
specific mass channels of acetonitrile and gain the most relevant information
about the fragmentation dynamics.

_Ab initio_ molecular dynamics simulations were performed, using the Atom
Centered Density Matrix Propagation method (ADMP) [37, 38, 39] as implemented
in the Gaussian16 program [40]. The electronic structure was computed,
employing the density functional theory (DFT), in particular, the B3LYP
functional [41, 42] in combination with the 6-31++G(d,p) basis set [43, 44].
To ensure the adiabaticity of the dynamics, we established a time step of
$\Delta t$ = 0.1 fs and a fictitious electron mass of $\mu$ = 0.1 amu.
Mimicking the experimental conditions, we considered a vertical double
ionization of the acetonitrile molecule in a Franck–Condon manner, introducing
a given amount of internal energy, $E_{exc}$ = 3 eV, which was randomly
redistributed over the nuclear degrees of freedom. This value was chosen to
ensure that the energy available in the system agrees with the experimentally
measured KER. 500 trajectories were propagated up to 1ps. We analyzed the
statistical population of the different channels considering separated
fragments with distances larger than 2.5 Å between atoms forming a moiety.
Charges in the different atomic positions were computed from a Mulliken
population analysis in the last step of each trajectory. We have successfully
used this computational strategy in the past to infer the fragmentation
dynamics in ionized molecules and clusters in the gas phase (see e.g. [45, 46,
47, 48, 34, 35, 49, 50, 36, 51]).

## Acknowledgements

The experimental work was funded by the National Science Foundation under
award No. 1700551. The theory was supported by the MICINN (Spanish Ministry of
Science and Innovation) projects PID2019-105458RB-I00 and
PID2019-110091GB-I00, funded by MCIN/AEI/10.13039/501100011033, the ‘Severo
Ochoa’ Programme for Centres of Excellence in R $\&$ D (CEX2020-001039-S) and
the ‘María de Maeztu’ Programme for Units of Excellence in R $\&$ D
(CEX2018-000805-M). We acknowledge the generous allocation of computer time at
the Centro de Computación Científica at the Universidad Autónoma de Madrid
(CCC-UAM).

## References

  * Bowman and Shepler [2011] J. M. Bowman and B. C. Shepler, Annu. Rev. Phys. Chem. 62, 531 (2011).
  * Bowman and Houston [2017] J. M. Bowman and P. L. Houston, Chem. Soc. Rev. 46, 7615 (2017).
  * Suits [2020] A. G. Suits, Annu. Rev. Phys. Chem. 71 (2020).
  * Townsend  _et al._ [2004] D. Townsend, S. A. Lahankar, S. K. Lee, S. D. Chambreau, A. G. Suits, X. Zhang, J. Rheinecker, L. Harding, and J. M. Bowman, Science 306, 1158 (2004).
  * Houston and Kable [2006] P. Houston and S. Kable, Proc. Natl. Acad. Sci. U.S.A. 103, 16079 (2006).
  * Heazlewood  _et al._ [2008] B. R. Heazlewood, M. J. Jordan, S. H. Kable, T. M. Selby, D. L. Osborn, B. C. Shepler, B. J. Braams, and J. M. Bowman, Proc. Natl. Acad. Sci. U.S.A. 105, 12719 (2008).
  * Maeda  _et al._ [2010] S. Maeda, K. Ohno, and K. Morokuma, J. Phys. Chem. Lett. 1, 1841 (2010).
  * Grubb  _et al._ [2012] M. P. Grubb, M. L. Warter, H. Xiao, S. Maeda, K. Morokuma, and S. W. North, Science 335, 1075 (2012).
  * Nakamura  _et al._ [2015] M. Nakamura, P.-Y. Tsai, T. Kasai, K.-C. Lin, F. Palazzetti, A. Lombardi, and V. Aquilanti, Faraday Discuss. 177, 77 (2015).
  * Lombardi  _et al._ [2016] A. Lombardi, F. Palazzetti, V. Aquilanti, H.-K. Li, P.-Y. Tsai, T. Kasai, and K.-C. Lin, J. Phys. Chem. A 120, 5155 (2016).
  * Rauta and Maiti [2016] A. K. Rauta and B. Maiti, Chem. Phys. Lett. 661, 83 (2016).
  * Ekanayake  _et al._ [2017] N. Ekanayake, M. Nairat, B. Kaderiya, P. Feizollah, B. Jochim, T. Severt, B. Berry, K. R. Pandiri, K. D. Carnes, S. Pathak, _et al._ , Sci. Rep. 7, 1 (2017).
  * Ekanayake  _et al._ [2018] N. Ekanayake, T. Severt, M. Nairat, N. P. Weingartz, B. M. Farris, B. Kaderiya, P. Feizollah, B. Jochim, F. Ziaee, K. Borne, _et al._ , Nat. Commun. 9, 1 (2018).
  * Livshits  _et al._ [2020] E. Livshits, I. Luzon, K. Gope, R. Baer, and D. Strasser, Communications Chemistry 3, 49 (2020).
  * Furukawa  _et al._ [2005] Y. Furukawa, K. Hoshina, K. Yamanouchi, and H. Nakano, Chem. Phys. Lett. 414, 117 (2005).
  * Hoshina  _et al._ [2008] K. Hoshina, Y. Furukawa, T. Okino, and K. Yamanouchi, J. Chem. Phys. 129, 104302 (2008).
  * Hoshina and Tsuge [2010] K. Hoshina and M. Tsuge, Chem. Phys. Lett. 489, 154 (2010).
  * Kraus  _et al._ [2011] P. M. Kraus, M. C. Schwarzer, N. Schirmel, G. Urbasch, G. Frenking, and K.-M. Weitzel, J. Chem. Phys. 134, 114302 (2011).
  * Okino  _et al._ [2012] T. Okino, A. Watanabe, H. Xu, and K. Yamanouchi, Phys. Chem. Chem. Phys. 14, 4230 (2012).
  * Schirmel  _et al._ [2013] N. Schirmel, N. Reusch, P. Horsch, and K.-M. Weitzel, Faraday Discuss. 163, 461 (2013).
  * Kotsina  _et al._ [2015] N. Kotsina, S. Kaziannis, and C. Kosmidis, Int. J. Mass Spectrom. 380, 34 (2015).
  * Geballe and Oka [1996] T. Geballe and T. Oka, Nature 384, 334 (1996).
  * McCall  _et al._ [1998] B. McCall, T. Geballe, K. Hinkle, and T. Oka, Science 279, 1910 (1998).
  * Herbst and Klemperer [1973] E. Herbst and W. Klemperer, Astrophys. J. 185, 505 (1973).
  * Endo  _et al._ [2020] T. Endo, S. P. Neville, V. Wanie, S. Beaulieu, C. Qu, J. Deschamps, P. Lassonde, B. E. Schmidt, H. Fujise, M. Fushitani, _et al._ , Science 370, 1072 (2020).
  * Ullrich  _et al._ [1997] J. Ullrich, R. Moshammer, R. Dörner, O. Jagutzki, V. Mergel, H. Schmidt-Böcking, and L. Spielberger, J. Phys. B: At., Mol. Opt. Phys. 30, 2917 (1997).
  * Ullrich  _et al._ [2003] J. Ullrich, R. Moshammer, A. Dorn, R. Dörner, L. P. H. Schmidt, and H. Schmidt-Böcking, Rep. Prog. Phys. 66, 1463 (2003).
  * Eland [1991] J. Eland, Laser Chem. 11, 259 (1991).
  * Boll  _et al._ [2022] R. Boll, J. M. Schäfer, B. Richard, K. Fehre, G. Kastirke, Z. Jurek, M. S. Schöffler, M. M. Abdullah, N. Anders, T. M. Baumann, _et al._ , Nat. Phys. 18, 423 (2022).
  * Bhattacharyya  _et al._ [2022] S. Bhattacharyya, K. Borne, F. Ziaee, S. Pathak, E. Wang, A. S. Venkatachalam, X. Li, N. Marshall, K. D. Carnes, C. W. Fehrenbach, _et al._ , J. Phys. Chem. Lett. 13, 5845 (2022).
  * Pathak  _et al._ [2020] S. Pathak, R. Obaid, S. Bhattacharyya, J. Bürger, X. Li, J. Tross, T. Severt, B. Davis, R. C. Bilodeau, C. A. Trallero-Herrero, _et al._ , J. Phys. Chem. Lett. 11, 10205 (2020).
  * Xie  _et al._ [2015] X. Xie, C. Wu, Y. Yuan, X.-Z. Li, C. Wu, P. Wang, Y. Deng, Y. Liu, and Q. Gong, AIP Adv. 5, 097213 (2015).
  * Neumann  _et al._ [2010] N. Neumann, D. Hant, L. P. H. Schmidt, J. Titze, T. Jahnke, A. Czasch, M. Schöffler, K. Kreidi, O. Jagutzki, H. Schmidt-Böcking, _et al._ , Phys. Rev. Lett. 104, 103201 (2010).
  * Kling  _et al._ [2019] N. G. Kling, S. Díaz-Tendero, R. Obaid, M. Disla, H. Xiong, M. Sundberg, S. Khosravi, M. Davino, P. Drach, A. Carroll, _et al._ , Nat. Commun. 10, 1 (2019).
  * McDonnell  _et al._ [2020] M. McDonnell, A. C. LaForge, J. Reino-Gonzalez, M. Disla, N. G. Kling, D. Mishra, R. Obaid, M. Sundberg, V. Svoboda, S. Díaz-Tendero, _et al._ , J. Phys. Chem. Lett. 11, 6724 (2020).
  * Mishra  _et al._ [2022] D. Mishra, J. Reino-González, R. Obaid, A. C. LaForge, S. Díaz-Tendero, F. Martín, and N. Berrah, Phys. Chem. Chem. Phys. 24, 433 (2022).
  * Schlegel  _et al._ [2001] H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. Daniels, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 114, 9758 (2001).
  * Iyengar  _et al._ [2001] S. S. Iyengar, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 115, 10291 (2001).
  * Schlegel  _et al._ [2002] H. B. Schlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 117, 8694 (2002).
  * Frisch  _et al._ [2016] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, “Gaussian16 Revision C.01,” (2016), gaussian Inc. Wallingford CT.
  * Becke [1993] A. D. Becke, J. Chem. Phys. 98, 5648 (1993), https://doi.org/10.1063/1.464913 .
  * Lee  _et al._ [1988] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  * Clark  _et al._ [1983] T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. V. R. Schleyer, J. Comput. Chem. 4, 294 (1983).
  * Krishnan  _et al._ [1980] R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, J. Chem. Phys. 72, 650 (1980).
  * Maclot  _et al._ [2013] S. Maclot, D. G. Piekarski, A. Domaracka, A. Meéry, V. Vizcaino, L. Adoui, F. Martín, M. Alcamí, B. A. Huber, P. Rousseau, _et al._ , J. Phys. Chem. Lett. 4, 3903 (2013).
  * Kukk  _et al._ [2015] E. Kukk, D. Ha, Y. Wang, D. G. Piekarski, S. Díaz-Tendero, K. Kooser, E. Itälä, H. Levola, M. Alcamí, E. Rachlew, and F. Martín, Phys. Rev. A 91, 043417 (2015).
  * Piekarski  _et al._ [2015] D. G. Piekarski, R. Delaunay, S. Maclot, L. Adoui, F. Martín, M. Alcamí, B. A. Huber, P. Rousseau, A. Domaracka, and S. Díaz-Tendero, Phys. Chem. Chem. Phys. 17, 16767 (2015).
  * Maclot  _et al._ [2016] S. Maclot, R. Delaunay, D. G. Piekarski, A. Domaracka, B. Huber, L. Adoui, F. Martín, M. Alcamí, L. Avaldi, P. Bolognesi, _et al._ , Phys. Rev. Lett. 117, 073201 (2016).
  * Rousseau  _et al._ [2020] P. Rousseau, D. G. Piekarski, M. Capron, A. Domaracka, L. Adoui, F. Martin, M. Alcami, S. Diaz-Tendero, and B. A. Huber, Nat. Commun. 11, 3818 (2020).
  * Barreiro-Lage  _et al._ [2021] D. Barreiro-Lage, P. Bolognesi, J. Chiarinelli, R. Richter, H. Zettergren, M. H. Stockett, L. Carlini, S. Diaz-Tendero, and L. Avaldi, J. Phys. Chem. Lett. 12, 7379 (2021).
  * Ganguly  _et al._ [2022] S. Ganguly, D. Barreiro-Lage, N. Walsh, B. Oostenrijk, S. t. L. Sorensen, S. Diaz-Tendero, and M. Gisselbrecht, Commun. Chem. 5, 16 (2022).