# THz-assisted microscopy of silica matrix for biological materials
encapsulation: a theoretical and experimental study

Matteo De Tullio  Univ Rouen Normandie, INSA Rouen Normandie, CNRS, Normandie
Univ, GPM UMR 6634, F-76000 Rouen, France  and Giovanni Novi Inverardi
Department of Physics, University of Trento and European Centre for
Theoretical Studies in Nuclear Physics and Related Areas (ECT*-FBK) and Trento
Institute for Fundamental Physics and Applications (TIFPA-INFN), Trento, Italy
Jonathan Houard  Univ Rouen Normandie, INSA Rouen Normandie, CNRS, Normandie
Univ, GPM UMR 6634, F-76000 Rouen, France  Marc Ropitaux  Univ Rouen
Normandie, GLYCOMEV UR4358, SFR Normandie Végétal FED 4277, Innovation Chimie
Carnot, IRIB, F-76000 Rouen, France.  Ivan Blum  Univ Rouen Normandie, INSA
Rouen Normandie, CNRS, Normandie Univ, GPM UMR 6634, F-76000 Rouen, France
Francesco Carnovale  Department of Physics, University of Trento and European
Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT*-FBK)
and Trento Institute for Fundamental Physics and Applications (TIFPA-INFN),
Trento, Italy  Gianluca Lattanzi  Department of Physics, University of Trento
and Trento Institute for Fundamental Physics and Applications (TIFPA-INFN),
Trento, Italy  Simone Taioli  European Centre for Theoretical Studies in
Nuclear Physics and Related Areas (ECT*-FBK) and Trento Institute for
Fundamental Physics and Applications (TIFPA-INFN), Trento, Italy  Gustav
Eriksson  Department of Chemistry and Chemical Engineering, Chalmers
University of Technology Gothenburg 41296, Sweden  Mats Hulander  Department
of Chemistry and Chemical Engineering, Chalmers University of Technology
Gothenburg 41296, Sweden  Martin Andersson  Department of Chemistry and
Chemical Engineering, Chalmers University of Technology Gothenburg 41296,
Sweden  Angela Vella<EMAIL_ADDRESS>Univ Rouen Normandie, INSA
Rouen Normandie, CNRS, Normandie Univ, GPM UMR 6634, F-76000 Rouen, France
Tommaso Morresi<EMAIL_ADDRESS>European Centre for Theoretical Studies in
Nuclear Physics and Related Areas (ECT*-FBK) and Trento Institute for
Fundamental Physics and Applications (TIFPA-INFN), Trento, Italy

###### Abstract

In this study, we use THz-assisted atom probe tomography (APT) to analyse
silica matrices used to encapsulate biomolecules. This technique provides the
chemical composition and 3D structure without significantly heating the
biosample, which is crucial for studying soft organic molecules such as
proteins. Our results show that THz pulses and a positive static field trigger
controlled evaporation of silica matrices, enabling 4D imaging with chemical
sensitivity comparable to UV laser-assisted APT. To support the interpretation
of these experimental results, we devise a computational model based on time-
dependent density functional theory to describe the interaction between silica
matrices and THz radiation. This model captures the nonlinear dynamics driven
by THz-pulses and the interplay between the THz source and the static electric
field in real time. This interdisciplinary approach expands the capabilities
of APT and holds promise for other THz-based analyses offering new insights
into material dynamics in complex biological environments.

Significant progress has been made in the development of powerful tools for
visualising protein structures with sub-nano resolution, such as X-ray
crystallography 1 and cryogenic electron microscopy 2. To gain new insights in
the fields of biology and medicine, these 3D spatial characterisation
techniques need to be further developed into 4D approaches that would provide
information on the local chemical composition as a fourth dimension.

APT is one such 4D imaging technique based on the controlled field evaporation
of atoms from a nanometric needle-shaped sample, achieving sub-nanometre
spatial resolution and high chemical sensitivity for the entire periodic table
and its isotopes3, 4. It has already been used to determine protein and
biomaterial structures making it a promising technique for biological research
5, 6, 7, 8, 9, 10. However, a significant drawback of UV laser-assisted APT is
the heating of specimens due to the use of a laser beam in the ultraviolet
spectral region, which can affect the integrity of soft organic molecules,
such as proteins 11, 12. To address this issue, terahertz (THz) lasers 13, 14
have been proposed as promising candidates for APT-based microscopy, offering
high spatial resolution and chemical sensitivity without damaging the
structure of biological samples. In addition, a breakthrough was achieved in
sample preparation for encapsulating biological samples at room temperature.
In this approach, a physiological solution containing the protein sample and
matrix precursor molecules, typically orthosilicates, is dried to obtain a
glassy phase 7that provides a solid support for measurements with APT without
distorting the native protein structure 15, 16. However, the use of THz-
assisted APT has not yet been explored for materials with larger band gaps,
and a comprehensive understanding of the interaction of THz pulses with the
silica matrix is still lacking.

Therefore, a thorough experimental and theoretical investigation of field
emission dynamics under THz pulse illumination is expected to open up
significant opportunities for the development of APT and other THz-based
analyses, such as time-dependent spectroscopy and THz-driven scanning
tunnelling microscopy 17, 18.

Computational approaches such as molecular dynamics (MD) and ab initio methods
are valuable for understanding the principles of APT applied to soft matter.
For instance, classical MD simulations have been used to study the early
stages of protein embedding processes 19. However, to investigate THz-based
ion field emission through computational modeling, it is essential to
explicitly consider the electronic degrees of freedom. In this context, Time-
Dependent Density Functional Theory (TDDFT) 20 can be combined with MD.
Unfortunately, the high computational cost of TDDFT to describe evaporation
processes limits the size of the simulated system for a given level of
accuracy. Currently, there are few in silico studies applying TDDFT to THz-
induced ion field emission, with only a few focusing on the use of UV
femtosecond laser pulses 21 or static methods 22, 23, 24.

In this work, we present a computational model based on TDDFT to study the
silica matrix under an electrostatic field and THz laser irradiation using
Si(OH)4 as a model. Si(OH)4 is the precursor molecule of the amorphous solid
matrix synthesized in this work using the sol–gel method for encapsulating
biomolecules 16. We emphasize the importance of a time-dependent approach to
capture the nonlinear effects related to the interaction of the sample with
the THz source and static electric field, leading to a dynamical change of the
interaction forces within the sample. This approach allows us to identify the
critical THz laser fields and mechanisms for the ion evaporation process of
the hydroxyl groups of Si(OH)4.

We have also carried out experiments to investigate amorphous silica
structures using a nanometric needle-shaped specimen. We successfully achieved
controlled evaporation of the silica matrix with THz pulses and a positive
static field. Our APT setup is equipped with a detection surface for spatially
resolved mass spectrometry analysis, providing information about the impact
position and time of flight of the ions 25. By controlling the evaporation
layer by layer, we reconstructed the 3D atomic positions of the sample and
determined the chemical nature by calculating the mass-to-charge ratio 4. Our
approach allowed us to generate 4D images of the silica sample with a signal-
to-noise ratio (S/N) and mass resolution comparable to those of UV laser-
assisted APT. The compositional analysis showed a slight decrease in oxygen
loss with increasing electric field, which can be partially explained by the
charge analysis from numerical calculations. Our study shows that THz-based
ion field emission is possible in insulating materials and that the charge
state of evaporated atoms at the THz and static field thresholds obtained in
our first-principles simulations agrees well with the measured values of the
mass spectra.

Silica matrices are amorphous materials with various bonding schemes and
substitutional defects. The studied molecule is a realistic model for the
typical Si-O bond in silica, and its size allows affordable first-principles
calculations. Our numerical simulations evolved the nuclei (ions) using
Ehrenfest dynamics 26, and the electronic wave function was propagated using
the TDDFT formalism (see Supplementary Information for more details). In this
framework, atomic ions and electrons are exposed to tunable external fields,
which can be either static ($\mathrm{E_{dc}}$), time-dependent
($\mathrm{E_{laser}}$) or both. The total external field (static + time-
dependent) is expressed by 21:

$v_{ext}(\mathbf{x},t)=\mathbf{x}\cdot\hat{\mathbf{e}}\left[\mathrm{E_{dc}}+\mathrm{E_{laser}}\
\text{exp}\left(-\frac{(t-t_{0})^{2}}{\sigma^{2}}\right)\text{cos}(\omega
t)\right]$ (1)

where $\hat{\mathbf{e}}$ is the polarisation vector, $\frac{\omega}{2\pi}=5$
THz, $t_{0}$=50 fs, and $\sigma$=20 fs.

The laser field and molecule interaction is described using the dipole
approximation. In our simulations, the box edge is $a=2$ nm, and the laser
wavelength is $\lambda\approx 6\times 10^{4}$ nm, satisfying the dipole
condition $\lambda\gg a$. All calculations use the OCTOPUS code 27, 28, with
an additional internal function for handling static external fields
interacting with ions.

The DFT equilibrium structure of the Si(OH)4 molecule is used as the starting
configuration for the system’s temporal evolution. The molecule has a
tetrahedral shape with the silicon atom coordinated with four hydroxyl groups
29. The molecule’s geometry is shown in Figure 1a. The pair distribution
function in the lower part of the figure indicates that at equilibrium, the
bond lengths for all O-H and O-Si pairs are equal.

Figure 1: Si(OH)4 molecule as a function of the static external field: (a)
$\mathbf{E}$=(0,0,0) [V/nm]; (b) $\mathbf{E}$=(5,0,0) [V/nm]; (c)
$\mathbf{E}$=(10,0,0) [V/nm]. The insets (i) at the top left of panels (a),
(b) and (c) are the geometries, while in the insets at the top right (ii) the
electron localisation function (ELF) at the isosurface level 0.66 has been
added to the structure; the insets at the bottom (iii) show the corresponding
pair distribution functions. The geometries were optimised by fixing the Si
atom, the y- and z-coordinates of the oxygen atom O$1$ and the y-coordinate of
the corresponding hydrogen atom H$1$.

When an external static electric field ($\mathbf{E}_{\text{dc}}$) is applied,
the molecule polarizes in the direction of the field, as shown in Figures 1b
and 1c for fields of $\mathbf{E}_{\text{dc}}=(5,0,0)$ and
$\mathbf{E}_{\text{dc}}=(10,0,0)$ [V/nm]. The more intense the static field,
the more strongly the hydroxyl groups are oriented in the direction of the
field, causing the $\widehat{\mathrm{Si-O1-H1}}$ angle in the Si(OH)4 molecule
to increase compared to the unpolarised case (Figure 1a) and the bond distance
between Si and O$1$ to expand to 2.4 Å.

During the dynamics, the positions of the silicon atom and the O2H2, O3H3, and
O4H4 ions are fixed to prevent the molecule from drifting. We are particularly
interested in the evaporation of the O1 and H1 atoms in response to static and
oscillating fields, resembling the evaporation of individual atoms from a
solid surface 21. For the time-dependent simulations, we focused on two static
fields: $\mathrm{E}_{{\mathrm{dc}},x}=0$ and $\mathrm{E}_{{\mathrm{dc}},x}=10$
[V/nm] (see Figures 2a,b and c,d) using therefore as initial configurations
the geometries in Figures 1a,c respectively. We note that these configurations
are only a small subset of the many possible alignments between the external
fields and the molecule. We specifically chose the configuration where the
O1-Si bond is aligned along the $x$-axis to uniquely evaluate the
contributions of the Coulomb force, static electric field, and time-dependent
laser for each value of the external fields. These results should be
considered as an upper bound for the real physical system, since different
angles between fields and bonds as well as other influencing factors such as
defects may exist in reality. We find a clear difference in the behaviour of
the system, which of course depends on the intensity of the external fields.
Results of the time-dependent simulations are shown in Figure 2. In
particular, the distance as function of time between the hydrogen atom H$1$
and the Si atom is shown in Figures 2a,c and the distance between the oxygen
atom O1 and the Si atom in Figures 2b,d. Different symbols represent different
laser field intensities.

Figure 2: (a,c): Distances between the hydrogen atom H1 and silicon as a
function of time. (b,d): Distances between the oxygen atom O1 and silicon as a
function of time. The labels follow the convention of Figure 1. (e) Plot of
the electronic density’s isosurface equal to $10^{-4}$ and ions for five
different snapshots of the time evolution in the case of
$\mathrm{E}_{{\mathrm{laser}},x}(t=50~{}\mathrm{fs})=61.4$ [V/nm]; the upper
plots report the evolution without static electric field while the lower ones
with $\mathrm{E}_{\mathrm{dc},x}=10$ [V/nm].

A laser field with a peak amplitude of 27.5 [V/nm] can evaporate hydrogen in
the presence of a static external field of $\mathrm{E}_{{\mathrm{dc}},x}=10$
[V/nm] (see blue dots in Figure 2c), but not without the static field (Figure
2a). Predicting this behaviour at the static level by calculating the PES is a
challenge due to the non-linear interplay between the fields leading to
different ionisations of the system. For oxygen, a larger amplitude of 54.5
[V/nm] is required to break the Si-O1 bond, as shown by the red diamonds in
Figures 2b,d. The distance between hydrogen and silicon increases rapidly,
making it easier to determine the evaporation point, while the slower increase
in the distance of oxygen from silicon makes it more difficult to determine
the evaporation point within the time frame of our simulations.

Nevertheless, our ab initio calculations provide information about the forces
and velocities acting on the atoms allowing us to determine the kinetic energy
of the oxygen atom. The analysis of this information, together with the charge
density, reveals the primary role of ionisation in the evaporation of oxygen.
Intense laser fields lead to increased ionisation of the fixed part of the
molecule. As the evaporating H1 and O1 atoms are always positively charged,
this results in a stronger Coulomb repulsion between the evaporating atoms and
the fixed part of the molecule, pushing the atoms far away, even in the
absence of static external fields. Figure 2e shows the electronic density at
the isosurface level $10^{-4}$ as well as the ions for a peak laser amplitude
of 61.4 [V/nm] in the absence of a static field (upper plot) and in the
presence of $\mathrm{E}_{\mathrm{dc},x}=10$ [V/nm] (lower plot).

$E_{\text{dc},x}$ $\left[\frac{V}{\text{nm}}\right]$ |  $E_{\text{las},x}(t=50\ fs)$ $\left[\frac{V}{\text{nm}}\right]$ | $N^{*}_{mol}$ | $N^{*}_{H1}$ | $N^{*}_{O{1}}$ |  $d_{H{1}-Si}$ [Å] |  $d_{O1-Si}$ [Å]
---|---|---|---|---|---|---
0 | 8.7 | 0.0 | - | - | 2.30 | 1.61
10 | 8.7 | 0.1 | - | - | 2.36 | 1.59
0 | 27.5 | 0.8 | - | - | 2.39 | 1.58
10 | 27.5 | 0.6 | 1.0 | - | > 10 | 1.58
0 | 41.1 | 1.2 | 1.0 | - | > 10 | 1.60
10 | 41.1 | 2.0 | 1.0 | - | > 10 | 1.49
0 | 48.8 | 1.6 | 1.0 | - | > 10 | 1.71
10 | 48.8 | 2.6 | 1.0 | - | > 10 | 2.20
0 | 54.5 | 2.4 | 1.0 | 0.5 | > 10 | 3.52
10 | 54.5 | 2.5 | 1.0 | 0.9 | > 10 | 4.37
0 | 61.4 | 2.7 | 1.0 | 0.7 | > 10 | 4.88
10 | 61.4 | 2.8 | 1.0 | 1.0 | > 10 | 6.27
0 | 86.7 | 3.9 | 1.0 | 1.9 | > 10 | > 10
10 | 86.7 | 4.0 | 1.0 | 2.1 | > 10 | > 10
Table 1: Summary of the results from the Ehrenfest dynamics simulations.
$d_{H1-Si}$ and $d_{O1-Si}$ are the distances between the hydrogen atom H$1$
and the Si atom and between the oxygen atom O$1$ and the Si atom;
$N^{*}_{mol}$ is the final total charge of the fixed molecule, while
$N^{*}_{O1}$ and $N^{*}_{H1}$ are the charges of the evaporating oxygen and
hydrogen atoms. The latter two values are only given if the distances between
these atoms and the Si are greater than 3 Å. The number of electrons emitted
is calculated as $N^{*}_{mol}+N^{*}_{O1}+N^{*}_{H1}$.

The snapshots illustrate the density at different time steps of the
simulation. As expected, the static field clearly accelerates the escaping
hydrogen and oxygen atoms. Additionally, we observe a significant difference
in the charge densities of the same oxygen atom with and without a static
field. The results for the charge state and distance reached by the atoms at
the end of the simulation are summarised in Table 1. We observed that the
hydrogen atom always evaporates fully ionised with a charge of +1, while the
charge of the evaporated oxygen atom increases with the laser amplitude. This
theoretical prediction is supported by our analysis of the mass spectra. We
also stress that, despite the analysis of an isolated molecule without
defects, the critical amplitudes of THz-fields for evaporating H1 and O1 are
close to the experimental values 13 as shown below.

On the experimental side, the silica samples were prepared using a lift-out
procedure with a focused ion beam scanning electron microscope (FIB-SEM) 30,
31 to cut the matrix into a thin needle shape with a radius of less than 100
nm (more details in the Supplementary Information).

Figure 3: (a) Schematic representation of a THz-assisted APT analysis chamber.
(b) Left panel: Measured electro-optical signal. Right panel: 3D atomic
distribution of Si (brown dots), oxygen atoms (green dots) and silicon dioxide
(orange dots), reconstructed with a starting radius of 25 nm and a cone angle
of 1.5∘. (c) Mass spectrum of a silica sample (data recorded on 1.1$\times
10^{6}$ atoms).

The spatially resolved mass spectrum is measured with a THz-assisted APT
system, as shown in Figure 3a 32. The nanometric needle-shaped specimen is
polarised with a DC voltage of -0.8 kV to 20 kV, resulting in an apex electric
field of over 10 V/nm. The atoms in the layers of the sample are evaporated
due to the strong electrostatic field and the irradiation with single-cycle
THz pulses, which trigger the evaporation process 13. Single-cycle THz pulses
are generated by a two-color laser field-induced plasma in N2 33, 34. The THz
pulses have a repetition rate of 1 kHz, matching that of the two-color laser
system, which consists of the fundamental harmonic (FH) and the second
harmonic (SH) of the 800 nm wavelength laser. The amplitude of the THz pulses
is proportional to the input power of the FH, which can be adjusted at the
entry of the plasma. The THz fields are predominantly generated by asymmetric
electron currents triggered by photoionisation in N2 gas. The generated THz
field is characterized via electro-optical (EO) scanning in a non-linear
crystal before it is focused on the APT sample. An example of the EO signal is
presented in Figure 3b on the left, where a signal of 3 mV corresponds to an
electric field of 20 MV/m. The THz pulse is then amplified by a factor of
several thousand at the apex of the sample, as shown in previous works on
metallic and non-metallic samples 35, 14. Therefore, the maximum THz field
amplitude ranges between 5-30 V/nm (further details in the Supplementary
Information).

After evaporation, ions are directed onto a detection surface using an
electrostatic field. The detector provides information about the impact
position and time of flight of the ions 25, enabling the reconstruction of the
3D structure of the sample. Each ejected ion is assigned an initial position
through back-projection, as shown in Figure 3b and its chemical nature is
determined by calculating the mass-to-charge state ratio, as shown in the mass
spectrum of Figure 3c 4. The 3D distribution of silicon, silica (SiO and
SiO2), and oxygen atoms is reconstructed by in-house software, with observed
inhomogeneity attributed to the porosity of the amorphous silica sample.

Figure 4: (a) Mass spectrum of the silica sample in the range 0-65 Da. (b)
Focus on the normalized Si2+ peak at 14 Da and Si+ peak at 28 Da. Data set of
250k atoms.

Figure 4 displays the results from the homogeneous part of a second silica
sample. The THz amplitude was maintained at its maximum during the analysis,
and the DC voltage was automatically adjusted to ensure a constant detection
rate during the field evaporation process, as shown in Figure 4b.

The mass spectrum (Figure 4a) exhibits characteristic peaks of silicon in the
2+ charge state at 14, 14.5, and 15 Da and in the 1+ state at 28, 29, and 30
Da. Additionally, peaks of molecular ions of silica, SiO and SiO2, are
observed at 44 and 60 Da for the 1+ charge state and less prominently at 22
and 30 Da for the 2+ state. The peak at 30 Da overlaps with the 30Si+ isotope,
which affects the natural isotopic distribution for the peaks at 28, 29 and 30
Da. The spectrum also indicates the presence of platinum in the 2+ state and
carbon in both the 1+ and 2+ states, likely from the C9H16Pt gas used in the
sample preparation by FIB-SEM milling. A peak at 23 Da suggests the presence
of sodium, likely a residue from the sol-gel preparation.

We see a peak at 17 Da for the OH+ ion, but there are 100 times fewer of these
ions than O+, indicating a low probability of evaporation. This is consistent
with our ab-initio simulation results, which also did not show evaporation of
these ions.

The Si2+ and Si+ peaks (Fig. 4 (b)) have a signal to noise ratio (SNR) larger
than 103 and a mass resolution at 10% of the peak maximum of 140 and 234,
respectively. These values are comparable to those obtained with UV laser-
assisted APT on oxides36, 37. The asymmetry of the Si2+ and Si+peaks, with a
long tail on the right side, suggests that thermal effects contribute to the
evaporation process

Figure 5: (a) Joint plot of the mass-to-charge ratio (as density plot) and the
detection rate (green line) as a function of the temporal progression of the
analysis on a data set of 360k ions in the range 13-15.5 Da to emphasize the
Si2+ species. (b) Focus on the normalized Si2+ peak at 14 Da. A signal-to-
noise ratio of more than 103 was measured. A THz amplitude of 2.43 mV in
combination with a DC bias voltage of 4 kV provided the optimal mass resolving
power at both 14 Da. (c) Red line: DC voltage temporal trend coupled with the
THz signal amplitudes used in each case; black line: temporal progression of
the electric field of the tip surface. (d) Total compositional ratio of
silicon and oxygen as a function of the electric field on the tip surface.

In the second set of analyses we increased the DC voltage and observed a
corresponding decrease in THz amplitude (Figure 5c). This was done to maintain
a constant detection rate and to study how the sample evaporation changes.
Figure 5a shows the mass-to-charge state ratio and detection rate against the
sequence of evaporated atoms, capturing the temporal progression of the
analysis. Slight shifts in the mass spectrum lines were observed, not directly
corresponding to changes in flux but associated with the change in THz
amplitude. The decrease in THz pulse corresponds to a decrease in the average
energy of the ions 14, leading to a shift towards high Da values in the mass
spectrum. The detection rate tends to decrease when the THz field amplitude
decreases by about 20%, despite a 20% increase in DC voltage. We used the Si+
to Si2+ count ratio to calculate the surface field at the top of our sample,
as per Kingham’s post-ionisation theory 38. However, the field calculation
shows a slight but costly increase in the field. This indicates that the
evaporation process has a non-linear behavior depending on the THz field, as
confirmed by our TDDFT simulations. We intentionally avoided increasing the DC
voltage further to prevent tip fracture from mechanical stress 39.

The Si2+ peaks in Figure 5b have a mass resolution at 10% of the peak maximum
of 187, which is consistent across all analysis conditions.

Figure 5d shows the analysis of the composition as a function of the electric
field on the tip surface. The resulting composition deviates from the expected
nominal ratio in silica of 66.6% oxygen and 33.3% silicon, instead showing a
tendency towards an equal distribution of 50% oxygen and 50% silicon,
suggesting an oxygen deficit during the evaporation process. The peak at 16 Da
can be attributed to either O${}_{2}^{2+}$ or O${}_{1}^{+}$. If the peak at 16
Da were attributed to O${}_{2}^{2+}$, the total oxygen content would increase
to 66%, which is the nominal composition. However, APT analysis of
18O-enriched $\alpha$-Fe2O3 ruled out the evaporation of the O${}_{2}^{2+}$
species and confirmed the presence of O${}_{1}^{+}$ 40. DFT simulations
predict the evaporation of the O${}_{1}^{+}$ species and the migration of
O${}_{1}^{-}$ in the cluster, potentially leading to the formation of O2 and
subsequent evaporation as O${}_{2}^{2+}$ 41.

Previous research has shown that a fraction of oxygen is lost in UV laser APT
analysis of oxides, suggesting that it evaporates as a neutral entity rather
than a positively charged ion 42. Furthermore, our results indicate that
changes in electric field had no significant effect on the composition of the
molecular species (see Supplementary Information).

In summary, in our study we used THz-assisted APT to analyze the chemical and
structural properties of silica matrices for encapsulating biomolecules. This
approach does not cause significant heating of the sample and enables high-
resolution 4D imaging with high chemical sensitivity, revealing the emission
of single ionic species over molecular ions. Our results agree with TDDFT
calculations, showing a significant change in charge density with applied
static fields. In particular, first-principles simulations indicate that
oxygen can be evaporated by THz field of 54 V/nm in 50 fs, at odds with
previous predictions of 0.5-3 ps 43, 44, 45. New simulations are planned to
investigate the evaporation of O2 molecular ions and the formation of
negatively charged oxygen atoms. This interdisciplinary approach demonstrates
the potential of THz-assisted APT for advanced microscopy of nanoscale surface
features, including in biological environments.

This action has received funding from the European Union under grant agreement
n. 101046651. MA acknowledges the support of the Swedish Research Council (VR)
and the Knut and Alice Wallenberg Foundation through their Wallenberg Academy
Fellows program.

##  1 Supplementary Information

###  1.1 Sample Synthesis

The silica samples were prepared according to a protocol developed for silica
encapsulation of proteins for APT analysis 7. In all experimental steps, the
water used was filtered through a Milli-Q (Millipore) system. All chemicals
were purchased from Sigma-Aldrich and used without further purification. A
sodium silicate solution (25.5 - 28.5 % SiO2 7.5 - 8.5 % Na2O), also known as
water glass, was used to form silica through a sol-gel process. Gelling was
initiated by adjusting the pH of the sodium silicate solution from alkaline to
physiological (pH 7.4). For this purpose, an acidic ion exchange column was
prepared by adding 1 g of ion exchange resin (Dowex 50WX8 hydrogen form 50-100
mesh, Sigma-Aldrich) to 0.5 g of glass wool in a 10 mL plastic syringe. The
ion exchange resin was activated by washing with 10 mL NaOH (4 M), followed by
water and finally 10 mL HCl (1 M). The activity was adjusted by rinsing with
small amounts of water until the eluate had a neutral pH. The original sodium
silicate solution was then diluted with water in a volume ratio of 1:3 before
1 mL of the solution was passed through the column, collected in a vial and
allowed to harden in an open vial at 37 °C for $\geq$ 48 h to obtain a solid
silicate material.

###  1.2 Experimental Setup

Experiments were conducted in a custom-built ultrahigh vacuum chamber with a
10 cm flight path and a base pressure of 5 $\times$ 10-11 mbar. The sample was
cooled to 50 K using a compression cryostat and polarized with a positive
applied DC bias, resulting in an intense field at the tip apex.

Field ion evaporation was initiated by a THz pulse generated by a two-color
ultra-short laser pulse. The laser system used was a Spitfire model from
Spectra Physics, which produced a 45-fs pulse in the infrared at a wavelength
of 800 nm, with a repetition rate of 1 kHz and a maximum energy of 3.5 mJ. The
output beam was a TEM00 Gaussian beam with a diameter of 8 mm. The beam was
split into two parts (90/10 in energy); the higher-energy beam was used for
terahertz generation and was mechanically chopped at 145 Hz. The lower-energy
beam was used as a probe for the temporal characterization of the terahertz
pulses in an EO sampling setup based on a GaP crystal outside the APT chamber.
The terahertz pulse was generated by a two-color plasma generation process 33.

The NIR laser beam was focused in air by a 40 cm focal lens and passed through
a BBO crystal placed before the focal plane to generate the second harmonic at
400 nm. The focused bi-chromatic femtosecond laser pulses emitted THz
radiation with very high field strengths. The energy of the second harmonic
was fixed at 2.65 mJ, corresponding to an intensity of 3$\times 10^{15}$ W
cm2. The amplitude of the emitted THz pulses could be controlled by changing
the aperture diameter of the pump beam. The produced THz field was then
focused on the sample inside the APT ultra-vacuum chamber.

Surface atoms were ionized by the combined action of the DC field and the THz
pulse and left the surface, flying toward a detector sensitive to the impact
position of the ion and its time of flight. The time of flight of the ions
allowed for the identification of their mass/charge ratio, i.e., (with some
limitations) their chemical identity. Furthermore, the impact position on the
detector made it possible to reconstruct the original position of each ion on
the tip through a back-projection algorithm. The Position Sensitive Detector
(PSD) was a delay line detector with improved multihit capabilities. It
consisted of a conventional 8 cm Roentdek delay line detector following a pair
of microchannel plates (MCPs) in chevron configuration. Signal output from
each end of the two wire pairs was digitized by means of fast digitizer boards
PXI DC 271 (or Keysight U1051A Acquiris). Ion flight time signals were taken
as output from the MCP, and the positioning was calculated using signals
provided by the delay line outputs. A deconvolution algorithm was used to
improve the multihit capabilities of the PSD 25.

###  1.3 Theoretical and computational details

The coupled electron-nuclear dynamics is formulated as follows (a.u. are
adopted) 20:

$i\frac{\partial\phi_{j}(\mathbf{r},t)}{\partial
t}=\Big{[}-\frac{\nabla^{2}}{2}+v_{ext}(\mathbf{r},t)+v_{H,xc}(\mathbf{r},t)+\hat{W}_{en}\left(\\{\mathbf{r}\\},\\{\mathbf{R}(t)\\}\right)\Big{]}\phi_{j}(\mathbf{r},t)$
(2) $M_{J}\frac{\partial^{2}\mathbf{R}_{J}}{\partial
t^{2}}=-\nabla_{\mathbf{R}_{J}}\Big{[}v_{ext}(\mathbf{R}_{J},t)+\hat{W}_{nn}\left(\\{\mathbf{R}(t)\\}\right)+\int
d^{3}r\
n(\mathbf{r},t)\hat{W}_{en}\left(\\{\mathbf{r}\\},\\{\mathbf{R}(t)\\}\right)\Big{]}$
(3)

where $\phi_{j}$ is the electronic wave function of the occupied orbitals,
$v_{ext}$ is external scalar potential acting on both electrons at position
$\\{\mathbf{r}\\}$ and nuclei at position $\\{\mathbf{R}\\}$, $v_{H,xc}$ is
the Hartree plus exchange-correlation potential, $\hat{W}_{nn}$ is the Coulomb
interaction between the nuclei, $M_{J}$ is the mass of the $J$-th nucleus,
$\hat{W}_{en}$ is the Coulomb interaction between the electrons and the
nuclei, and $n(\mathbf{r},t)\equiv\sum_{j}^{occ}|\phi_{j}(\mathbf{r},t)|^{2}$
is the electron density.

The TDDFT simulation box for the Si(OH)4 molecule is a 20 Å cube, which was
chosen to ensure that the TDDFT simulation preserves the energy of the system
without a laser field. The real space mesh spacing is 0.16 Å, resulting in a
total energy convergence below 5 meV. To prevent the electrons emitted by the
laser from being reflected at the edges of the box, we implement absorbing
boundary conditions that suppress the electronic density. For this purpose, we
use a mask function that is equal to one in the inner region and smoothly
decreases to zero in the edge region. For the TDDFT simulations, the time step
integrator is equal to $4\cdot 10^{-4}$ fs and in Eq. 1 we use $t_{0}$=50 fs
(corresponding to the maximum amplitude of the laser signal) and $\sigma$=20
fs 21. The ion trajectories and the electron wave functions are propagated for
100 fs, through 2.5$\cdot$105 total steps.

The initial geometries reported in Figure 1 were optimized using the FIRE
algorithm 46 with a force threshold of 0.03 eV/Å. To prevent global
translations and rotations of the molecule in the presence of the static
external field, the optimizations were carried out by fixing the central Si
atom, the $y$ and $z$ coordinates of the oxygen atom O$1$, and the $y$
coordinates of the hydrogen atom H$1$ (see Figure 1). The static electric
field values analysed in this work are $\mathbf{E}_{\text{dc}}=(0,0,0)$,
$\mathbf{E}_{\text{dc}}=(5,0,0)$ and $\mathbf{E}_{\text{dc}}=(10,0,0)\
\left[\frac{\text{V}}{\text{nm}}\right]$. It is worth noting that the maximum
(and minimum) value of the static fields was chosen so that the electrons
bound for $\mathbf{E}_{\text{dc}}=0$ remain bound to the molecule even in the
presence of the static external field. The explicit treatment of the electric
field is necessary to prevent the global minimum of the potential for the
electrons from moving to the edge of the box when the amplitude of the field
is large enough. Although this can be avoided by treating electric fields in
the framework of the modern polarisation theory 47, we retain the explicit
treatment of $\mathbf{E}_{\text{dc}}$ because it has been shown in previous
investigations 22, 23, 24, 21 to be useful for describing the evaporation
dynamics of atoms from molecules and periodic systems.

The nuclear-electron interaction is modeled using the pseudopotential method
within the local density approximation (LDA) 48 for the static DFT electron
exchange-correlation and the adiabatic LDA (ALDA) for the TDDFT simulations.
In this framework, the Si(OH)4 molecule has 32 valence electrons.

## References

  * Smyth and Martin 2000 Smyth, M.; Martin, J. x Ray crystallography. _Molecular Pathology_ 2000, _53_ , 8
  * Bai et al. 2015 Bai, X.-C.; McMullan, G.; Scheres, S. H. How cryo-EM is revolutionizing structural biology. _Trends in biochemical sciences_ 2015, _40_ , 49–57
  * Gault et al. 2006 Gault, B.; Vurpillot, F.; Vella, A.; Gilbert, M.; Menand, A.; Blavette, D.; Deconihout, B. Design of a femtosecond laser assisted tomographic atom probe. _Review of Scientific Instruments_ 2006, _77_ , 043705
  * Kelly and Miller 2007 Kelly, T. F.; Miller, M. K. Atom probe tomography. _Review of scientific instruments_ 2007, _78_
  * Grandfield et al. 2022 Grandfield, K.; Micheletti, C.; Deering, J.; Arcuri, G.; Tang, T.; Langelier, B. Atom probe tomography for biomaterials and biomineralization. _Acta Biomaterialia_ 2022, _148_ , 44–60
  * Perea et al. 2016 Perea, D. E.; Liu, J.; Bartrand, J.; Dicken, Q.; Thevuthasan, S. T.; Browning, N. D.; Evans, J. E. Atom probe tomographic mapping directly reveals the atomic distribution of phosphorus in resin embedded ferritin. _Scientific reports_ 2016, _6_ , 22321
  * Sundell et al. 2019 Sundell, G.; Hulander, M.; Pihl, A.; Andersson, M. Atom probe tomography for 3D structural and chemical analysis of individual proteins. _Small_ 2019, _15_ , 1900316
  * Qiu et al. 2020 Qiu, S.; Zheng, C.; Garg, V.; Chen, Y.; Gervinskas, G.; Li, J.; Dunstone, M. A.; Marceau, R. K.; Fu, J. Three-dimensional chemical mapping of a single protein in the hydrated state with atom probe tomography. _Analytical chemistry_ 2020, _92_ , 5168–5177
  * Sundell et al. 2017 Sundell, G.; Dahlin, C.; Andersson, M.; Thuvander, M. The bone-implant interface of dental implants in humans on the atomic scale. _Acta biomaterialia_ 2017, _48_ , 445–450
  * Karlsson et al. 2014 Karlsson, J.; Sundell, G.; Thuvander, M.; Andersson, M. Atomically resolved tissue integration. _Nano letters_ 2014, _14_ , 4220–4223
  * Vella et al. 2013 Vella, A.; Silaeva, E. P.; Houard, J.; Itina, T. E.; Deconihout, B. Probing the thermal response of a silicon field emitter by ultra-fast Laser Assisted Atom Probe Tomography. _Annalen der Physik_ 2013, _525_ , L1–L5
  * Houard et al. 2011 Houard, J.; Vella, A.; Vurpillot, F.; Deconihout, B. Three-dimensional thermal response of a metal subwavelength tip under femtosecond laser illumination. _Physical Review B_ 2011, _84_ , 033405
  * Vella et al. 2021 Vella, A.; Houard, J.; Arnoldi, L.; Tang, M.; Boudant, M.; Ayoub, A.; Normand, A.; Da Costa, G.; Hideur, A. High-resolution terahertz-driven atom probe tomography. _Science Advances_ 2021, _7_ , eabd7259
  * Karam et al. 2023 Karam, M.; Houard, J.; Damarla, G.; Rousseau, L.; Bhorade, O.; Vella, A. THz driven field emission: energy and time-of-flight spectra of ions. _New Journal of Physics_ 2023, _25_ , 113017
  * McCarroll et al. 2020 McCarroll, I. E.; Bagot, P.; Devaraj, A.; Perea, D. E.; Cairney, J. New frontiers in atom probe tomography: a review of research enabled by cryo and/or vacuum transfer systems. _Materials Today Advances_ 2020, _7_ , 100090
  * Sundell et al. 2019 Sundell, G.; Hulander, M.; Pihl, A.; Andersson, M. Atom Probe Tomography for 3D Structural and Chemical Analysis of Individual Proteins. _Small_ 2019, _15_ , 1900316
  * Cocker et al. 2013 Cocker, T. L.; Jelic, V.; Gupta, M.; Molesky, S. J.; Burgess, J. A.; Reyes, G. D. L.; Titova, L. V.; Tsui, Y. Y.; Freeman, M. R.; Hegmann, F. A. An ultrafast terahertz scanning tunnelling microscope. _Nature Photonics_ 2013, _7_ , 620–625
  * Wimmer et al. 2014 Wimmer, L.; Herink, G.; Solli, D. R.; Yalunin, S. V.; Echternkamp, K. E.; Ropers, C. Terahertz control of nanotip photoemission. _Nature Physics_ 2014, _10_ , 432–436
  * Novi Inverardi et al. 2023 Novi Inverardi, G.; Carnovale, F.; Petrolli, L.; Taioli, S.; Lattanzi, G. Silica In Silico: A Molecular Dynamics Characterization of the Early Stages of Protein Embedding for Atom Probe Tomography. _Biophysica_ 2023, _3_ , 276–287
  * Ullrich 2012 Ullrich, C. A. _Time-Dependent Density-Functional Theory: Concepts and Applications_ ; Oxford University Press, 2012
  * Silaeva et al. 2015 Silaeva, E. P.; Uchida, K.; Suzuki, Y.; Watanabe, K. Energetics and dynamics of laser-assisted field evaporation: Time-dependent density functional theory simulations. _Phys. Rev. B_ 2015, _92_ , 155401
  * Feibelman 2001 Feibelman, P. J. Surface-diffusion mechanism versus electric field: Pt/Pt(001). _Phys. Rev. B_ 2001, _64_ , 125403
  * Sánchez et al. 2004 Sánchez, C. G.; Lozovoi, A. Y.; Alavi, A. Field-evaporation from first-principles. _Molecular Physics_ 2004, _102_ , 1045–1055
  * Tamura et al. 2012 Tamura, H.; Tsukada, M.; McKenna, K. P.; Shluger, A. L.; Ohkubo, T.; Hono, K. Laser-assisted field evaporation from insulators triggered by photoinduced hole accumulation. _Phys. Rev. B_ 2012, _86_ , 195430
  * Costa et al. 2012 Costa, G. D.; Wang, H.; Duguay, S.; Bostel, A.; Blavette, D.; Deconihout, B. Advance in multi-hit detection and quantization in atom probe tomography. _Review of Scientific Instruments_ 2012, _83_
  * Marx and Hutter 2009 Marx, D.; Hutter, J. _Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods_ ; Cambridge University Press, 2009; p 11–84
  * Andrade et al. 2015 Andrade, X. et al. Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems. _Phys. Chem. Chem. Phys._ 2015, _17_ , 31371–31396
  * Tancogne-Dejean et al. 2020 Tancogne-Dejean, N. et al. Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems. _The Journal of Chemical Physics_ 2020, _152_ , 124119
  * Belton et al. 2012 Belton, D.; Deschaume, O.; Perry, C. An overview of the fundamentals of the chemistry of silica with relevance to biosilicification and technological advances. _FEBS J_ 2012, _279_ , 1710–1720
  * Miller et al. 2005 Miller, M.; Russell, K.; Thompson, G. Strategies for fabricating atom probe specimens with a dual beam FIB. _Ultramicroscopy_ 2005, _102_ , 287–298
  * Thompson et al. 2007 Thompson, K.; Lawrence, D.; Larson, D.; Olson, J.; Kelly, T.; Gorman, B. In situ site-specific specimen preparation for atom probe tomography. _Ultramicroscopy_ 2007, _107_ , 131–139
  * Blavette et al. 1993 Blavette, D.; Bostel, A.; Sarrau, J.-M.; Deconihout, B.; Menand, A. An atom probe for three-dimensional tomography. _Nature_ 1993, _363_ , 432–435
  * Bartel et al. 2005 Bartel, T.; Gaal, P.; Reimann, K.; Woerner, M.; Elsaesser, T. Generation of single-cycle THz transients with high electric-field amplitudes. _Optics Letters_ 2005, _30_ , 2805–2807
  * Kim et al. 2007 Kim, K.-Y.; Glownia, J. H.; Taylor, A. J.; Rodriguez, G. Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields. _Optics express_ 2007, _15_ , 4577–4584
  * Houard et al. 2020 Houard, J.; Arnoldi, L.; Ayoub, A.; Hideur, A.; Vella, A. Nanotip response to monocycle terahertz pulses. _Applied Physics Letters_ 2020, _117_
  * Marquis et al. 2010 Marquis, E. A.; Yahya, N. A.; Larson, D. J.; Miller, M. K.; Todd, R. I. Probing the improbable: Imaging C atoms in alumina. _Materials Today_ 2010, _13_ , 34–36
  * Vella et al. 2011 Vella, A.; Mazumder, B.; Da Costa, G.; Deconihout, B. Field evaporation mechanism of bulk oxides under ultra fast laser illumination. _Journal of Applied Physics_ 2011, _110_
  * Kingham 1982 Kingham, D. R. The post-ionization of field evaporated ions: A theoretical explanation of multiple charge states. _Surface Science_ 1982, _116_ , 273–301
  * Rigutti et al. 2017 Rigutti, L.; Venturi, L.; Houard, J.; Normand, A.; Silaeva, E. P.; Borz, M.; Malykhin, S. A.; Obraztsov, A. N.; Vella, A. Optical Contactless Measurement of Electric Field-Induced Tensile Stress in Diamond Nanoscale Needles. _Nano Letters_ 2017, _17_ , 7401–7409, PMID: 29095635
  * Bachhav et al. 2013 Bachhav, M.; Danoix, F.; Hannoyer, B.; Bassat, J. M.; Danoix, R. Investigation of O-18 enriched hematite ($\alpha$-Fe2O3) by laser assisted atom probe tomography. _International Journal of Mass Spectrometry_ 2013, _335_ , 57–60
  * Karahka and Kreuzer 2013 Karahka, M.; Kreuzer, H. Field evaporation of oxides: A theoretical study. _Ultramicroscopy_ 2013, _132_ , 54–59, IFES 2012
  * Mancini et al. 2014 Mancini, L.; Amirifar, N.; Shinde, D.; Blum, I.; Gilbert, M.; Vella, A.; Vurpillot, F.; Lefebvre, W.; Lardé, R.; Talbot, E.; others Composition of wide bandgap semiconductor materials and nanostructures measured by atom probe tomography and its dependence on the surface electric field. _The Journal of Physical Chemistry C_ 2014, _118_ , 24136–24151
  * Müller 1956 Müller, E. W. Field Desorption. _Phys. Rev._ 1956, _102_ , 618–624
  * Gomer 2004 Gomer, R. Field Desorption. _The Journal of Chemical Physics_ 2004, _31_ , 341–345
  * Kellogg 1984 Kellogg, G. Measurement of activation energies for field evaporation of tungsten ions as a function of electric field. _Physical Review B_ 1984, _29_ , 4304
  * Bitzek et al. 2006 Bitzek, E.; Koskinen, P.; Gähler, F.; Moseler, M.; Gumbsch, P. Structural Relaxation Made Simple. _Phys. Rev. Lett._ 2006, _97_ , 170201
  * Umari and Pasquarello 2002 Umari, P.; Pasquarello, A. Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field. _Phys. Rev. Lett._ 2002, _89_ , 157602
  * Perdew and Wang 1992 Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electron-gas correlation energy. _Phys. Rev. B_ 1992, _45_ , 13244–13249