# Effects of surface chemical modifications on the adhesion of metallic
interfaces. An high-throughput analysis.

E. Poli  Department of Physics and Astronomy, University of Bologna, 40127
Bologna, Italy<EMAIL_ADDRESS>M. Cutini  Department of Physics and
Astronomy, University of Bologna, 40127 Bologna, Italy  M. A. M. A. Nosir
Department of Physics and Astronomy, University of Bologna, 40127 Bologna,
Italy  O. Chehaimi  Department of Physics and Astronomy, University of
Bologna, 40127 Bologna, Italy  M. C. Righi  Department of Physics and
Astronomy, University of Bologna, 40127 Bologna, Italy<EMAIL_ADDRESS>
###### Abstract

Chemical interactions between two surfaces in contact play a crucial role in
determining the mechanical and tribological behavior of solid interfaces.
These interactions can be quantified via the adhesion energy, that is a
measure of the strength by which two surfaces bind together. Several works in
literature report how the presence of chemisorbed atoms at homo- and
heterogeneous solid-solid interfaces drastically change their proprieties. A
precise evaluation of how different species at solid contacts modulates their
adhesion would be extremely beneficial for a range of different technological
fields: from metallurgy to nuclear fusion. In this work we have used and high-
throughput approach to systematically explore the effects of the presence of
non-metallic elements, at different concentrations, on the adsorption and
adhesion energies of different homogeneous metallic interfaces. Together with
the databases for the adsorption and the adhesion energies, we calculated
several other properties such as the charge transferred at the interface, the
d-band edge shift for the substrate the Bond order and the interfacial density
redistribution for the hundreds of systems analyzed. These values were used to
define different trends with respect to chemical and concentration parameters
that could be useful for the development of engineered interfaces with
selected properties. In particular we noticed how the substrate with low
filling of d-band are the most prone to adsorb ad-atoms and how the adsorption
of almost all non-metallic elements decreases the adhesion energy of solid
interfaces, particularly in the case of Fluorine. Carbon and Boron were the
only two ad-atoms species that showed an opposite trend increasing the
adhesion energy instead.

##  1 Introduction

Solid-Solid interfaces are ubiquitous in natural and technological processes1,
2, 3. The mechanical behaviour of such interfaces is dictated by an interplay
between the macro4 and microscopic properties5 of the surfaces comprising the
contact. In fact, while the physical processes that govern friction are not
fully understood 6 it is known that the origins of frictional interactions are
the result of the shearing of microscopic junctions due to nano-asperities
present at contacting surfaces7, 8. In this context adhesive interactions are
extremely relevant given the high surface to volume ratio at nano-asperity
contacts. Adhesion is defined by the atomistic interactions of the mated
surfaces of which some relevant aspects are chemical composition, crystal
orientation and stacking, presence of adsorbates, impurities, defects and
segregates 9, 10, 11, 12, 13. In particular the presence of atomic adsorbates
at metallic interfaces can modify the tribological and mechanical properties
of interfaces with relevant implications in many technological fields such as
metallurgy14, nuclear fusion generator15, 16, 17, automotive, catalysis18,
clean fuel deployment19. The phenomena influenced by the presence of
adsorbates ranges from hydrogen embrittlement of steels19, 20 to
superlubricity21 to electrode adhesion22 to different substrates and many
more. A deeper knowledge of the effects of different adsorbates should allow
the prediction and optimization of the mechanical and the tribological
behaviour of metallic interfaces, opening the possibility for applications
where selected levels of adhesion and friction required is obtained through
chemical engineering. The interplay between theoretical analysis and
experimental evidence is pivotal in this direction. From a theoretical
standpoint the advent of high-throughput software has dramatically enhanced
the input that can be given to the experimental community to speed up/aide the
discovery and testing of novel compounds with selected properties23, 24, 25,
26, 27. This approach hence has become a valuable tool for material
advancements. In this paper we present an high-throughput study, based on an
extension of the TribChem28 program previously developed by our group, on how
the absorption of non-metal atomic species at homogeneous metal-metal
interfaces defines their adhesive and chemical properties. In particular we
focus on how the presence of difference adsorbates influences the Adsorption
and Adhesion energies given their relevance in defining the tribological
properties of an interface29. We also analyse how these central quantities are
correlated to several other chemical descriptors/properties. The study of
these fundamental systems gives us a perspective of the fundamental
interactions at plays at the boundary between metallic interfaces and is a
first step in order to move towards the comprehension of the physics of more
heterogeneus interfaces.

##  2 Computational Methods

###  2.1 Tribchem

The high-throughput screening of the influence on adhesion of adatoms at
homogeneous metal interfaces has been carried out using the TribChem software
developed by our group. The technical details of the code are outlined in a
separate paper28. TribChem role is to automate the creation of input files and
geometries, submit and run the theoretical calculations, and finally collect
and analyze the generated data in external databases. The code is mainly
written in python and leverages Fireworks30 and Atomate31 as workflow
managers. In addition TribChem is interfaced with Pymatgen32 and
MPInterfaces33 that are used to create, manipulate and process the
Input/Output operations. The DFT calculations are run by Tribchem using the
Vienna Ab initio Simulation Package (VASP34). We used for all the simultions
plane wave projector augmented wave (PAW35) basis and the Perdew-Burke-
Ernzerhof36 (PBE) -generalized gradient approximation of exchange correlation.

###  2.2 Workflow outline

  * •

Tribchem retrieves the structure for the material of interest from the
Material Project database.

  * •

The bulk structures are optimised converging their lattice parameters via a
fit with the Birch-Murnaghan equation of state.

  * •

The DFT energy convergence w.r.t the k-point mesh and the plane-wave cutoff is
then evaluated via an automated set of calculations.

  * •

The slabs of the materials selected are the obtained cutting the bulk along a
input/user specified direction. In our specific case the slabs were cut along
the most stable orientation of each crystalline structure: (111) for FCC
crystal, (110) for BCC crystal and (0001) for HPC ones.

  * •

The optimal slab thickness is selected via convergence of the surface energy.

  * •

The adatom is added to the slab in the different high symmetry sites depending
on the orientation of the surface (as an example: Long-bridge, short-bridge,
or top for the (110) BCC crystals. The structures obtained are then relaxed
until the convergence of the forces is minimized below 0.01 eV/Å.

  * •

The lowest energy configuration is then selected and the adsorption energy
with respect to the clean slab is calculated and saved in the Adsorption
Energies (EADS) database.

  * •

The interface is then built reflecting the slab (and not the adatom) through a
mirror plane parallel to the slab surface. The mirrored slab is then placed at
a distance along the z-axis equal to the interlayer spacing of the bulk. The
total supercell is finally formed adding 15 Å of vacuum to decouple the
Periodic Boundary Conditions (PBC) images along the z-direction. The structure
is then relaxed along the z-axis until the convergence of the forces is
minimized below 0.01 eV/Å.

  * •

The mirrored slab is then shifted laterally to match all the high symmetry
positions of the lower slab and. Exploring the different shifts for which the
high-symmetry points of the two surface match, and collecting the
corresponding total energies, a 2-dimensional potential energy surface (PES)
is built. The structure corresponding to the minima of the PES is then
selected. Finally, this structure is fully relaxed until the convergence of
the forces is minimized below 0.01 eV/Å.

The Adhesion energy (Eadh), charge density redistribution ($\rho_{dist}$),
Bader charges 37, Density Derived Eelectrostatic Charges (DDEC) 38, Bond Order
and d-band edge 39 are then calculated and added to the respective databases.
In particular in order to calculate the DDEC charges, Bond Order and d-band
edges, we used the codes VASPKIT40, and CHARGEMOL38. The capabilities of
running post-processing analysis has not been fully implemented in the
workflow yet and a more technical paper about these capabilities will follow.
This full procedure is run for 1x1, 2x1, 2x2 cell dimensions (along the a and
b cell axis respectively) considering in each case only one adatom per system
and simulating in this way three different coverages, $\theta=$ 1.0 monolayer
(ML), $\theta=$ 0.5 and $\theta=$ 0.25, respectively.

###  2.3 Subtrates and adatoms choices

The substrates choice was driven by our previous work (and hence the
availability of the Eadh and $\rho_{dist}$ data) on the clean homogeneous
interfaces of these metals. The series of substrates analysed comprises of:
Ag, Al, Au, Cr, Cu, Fe, Ir, Mg, Mo, Ni, Pt, Ti, W, and Zn. These metals were
chosen given their technological relevance and the fact that they constitute a
representative set of transition metals with different d-orbital occupations.
In addition Mg and Al were also considered as limit cases of metal which
reactivity is dictated by s and p-orbitals. The adatom selection comprises of
all reactive non-metals from the first second and third period (Cl excluded).
Boron was also included given its relevance in previous tribological studies.

##  3 Results and Discussion

The database of the adsorption energies for the different substrates with the
different adatoms adsorbed on top and the one of the adhesion energies of the
interfaces with the adatoms intercalated are shown in Figure1 A and B
respectively. The data are reported for all the three different coverages
considered. A more detailed picture of each individual system (i.e. each
adatom on every substrate) can be found in Fig S1 for the adsorption energies
and Fig S2 for the adhesion energies. The quantities used to compile these
databases are defined as follows:

$\ \ E_{ADS}=\frac{(E_{SLAB+nAd}-nE_{Ad}-E_{SLAB})}{n}$ (1)

Where ESLAB+nAd is the total energy of the substrate covered by n adatoms, EAd
is the energy of the isolated atom and ESLAB is the energy of the clean
substrate).

$\ \ E_{ADH}=\frac{1}{A}(E_{12}-E^{1}_{Ad}-E^{2})$ (2)

Where A is the supercell area, E12 si the total energy of the two slabs and
the adatom in contact, E${}^{1}_{Ad}$ is the energy of the isolated lower slab
with adsorbed the adatom and E2 is the energy of the isolated upper slab.

$\ \
\rho_{dist}=\frac{1}{2z_{0}}\int_{-z0}^{z0}\left|\frac{\rho_{12}-\rho_{1}^{Ad}-\rho_{2}}{A}\right|dz$
(3)

Where 2z0 is the interface distance, $\rho_{12}$ is the planar average of the
interface charge density $\rho_{1}^{Ad}$ is the planar average of the charge
density of the bottom slab with the adsorbed adatom and $\rho_{2}$ is the
planar average of the charge density of the bottom slab for each specific
system. $\rho_{dist}$ is a metric to quantify the accumulation of charge at
the interface 29.

$\ \ \epsilon_{d}=\frac{\int_{-\infty}^{\infty}{n_{d}(\epsilon})\epsilon
d\epsilon}{\int_{-\infty}^{\infty}{n_{d}(\epsilon})d\epsilon}$ (4)

Where $\epsilon_{d}$ is the center of d-projected DOS $\epsilon$ represents
the Kohn–Sham eigenvalues and n($\epsilon$) is the calculated total d-PDOS.
We also calculated for each system the bader37 and DDEC charges and the DDEC
bond order. DDEC is an atom in molecules charge partitioning where the total
electron density (n(r)) is partitioned into overlapping atomic densities
(ni(r)).

$\ \ n_{i}(r)=\frac{w_{i}(r)}{\sum_{k}(w_{k}(r))}n(r)$ (5)

The atomic partial charges are then computed by integrating the atomic
electron densities over all space

$\ \ q_{i}=z_{i}-N_{i}=Z_{i}-\sum{n_{i}(r)d^{3}r}$ (6)

where Ni is the number of electrons assigned to atom i and zi is its effective
nuclear charge. Using the same approach, higher-order atomic multipoles may be
computed as first-order, second-order, nth-order moments of the atomic
electron densities.The weighting function in the DDEC method is described so
that the atomic weights are simultaneously optimized to resemble the spherical
average of ni(r) and the density of a reference ion of the same element with
the same atomic population Ni . Following this approach, the assigned atomic
densities rapidly reach a converging multipole expansion of the QM
electrostatic potential and the resulting atomic populations are chemically
reasonable. Building on this charge partitioning the Bond Order of the
different atoms can be calculated.

A more detailed explanation of these method is beyond the scope of this paper
an the interested reader is remanded to Ref.38, 41

###  3.1 Adsorption Energy Analysis

Several general trends are observable regarding the adsorption energies:
i) The substrates with low filling of d-band are the most prone to adsorb
adatoms (see Figure1, top panel). This trend is consistent with d-band centre
theory for adsorption on metals for which the higher the d states are in
energy relative to the Fermi level, the higher in energy the antibonding
states are and the stronger the bond. This explain the lower propensity of
noble metal and Zinc to adsorb different species 42, 43 with energies in the
order of 4 - 2 eV for the 25% coverage, 3 - 2 eV for 50% coverage and 2 - 1 eV
for 100% coverage. Conversely, the adsorption energies for the adatoms on the
other substrate increase with the decrement of the d electrons. These trends
were all confirmed by the d-band edge values calculated for all system (Table
1 Supporting information). Aluminum (that has only one p electron) shows
adsorption energies that are comparable in magnitude with the ones of low
filled d-bands substrates. It also shows (and Magnesium too) similar energetic
trends to the low filled d-bands substrates with a high propensity to adsorb
C, O and N. The case of Magnesium is difficult to evaluate since for several
systems the adsorbed species percolates into the surface structure, in
particular C, N and O. This is probably due to the high reactivity of the
alkaline-earth metals and the large lattice parameter along z of Mg that
allows for easier percolation.
ii) The most favorable adatom to adsorb depends both by the substrate and the
coverage of the system(see Figure 2a). The only three systems for which the
coverage does not influence the most favorable species adsorbed are: Pt/C,
Ag/F and Al/O. In general the species with highest Adsorption Energy are the
ones with smaller atomic radii exception made for hydrogen (i.e. C,N ,O and
F). It is notable how at higher coverages the adsorption of Fluorine becomes
more advantageous with respect the other adatoms.
iii) H is the least energetically favourite species to adsorb on the different
substrates except for high concentrations where S and P are the least
favourite especially on coniage metals. Figure 2b in fact highlights how at
$\Theta=0.5$ and $0.25$ H is generally the least favoured adatom on the
different surfaces.
iv) Figure 2b show the average change in Adsorption energy passing from
$\Theta=1.0$ to $0.5$ and from $\Theta=0.5$ to $0.25$. It can be seen that for
atoms with smaller atomic radii the changes in Eads are smaller (i.e. the gain
in energy due to adsorption is less coverage dependent). In particular for
Hydrogen the adsorption energies remain on average the same, in the
concentration range analysed. Boron, Carbon, Sulfur and Phosphorous Eads show
instead a marked dependence by the coverage. Decreasing the coverage in these
cases in fact leads to a drastic increase of the adsorption energies that
seems to be correlated to a widening of the distance between the ad-atoms from
the surface (see Figure 1, bottom panel). Figure 2d shows in particular how
the B, C, P and S are on average 0.3-0.2 Å further away from the substrate at
$\Theta=1.0$. For the other species the distance from the substrate remains
almost the same exception made for F. These trends can in part be explained
calculating the Bond Order between the adatoms and the different species
present in the system. (see Figure 2e,f) For B, P, S (and in a more limited
way C) the increment in coverage leads to an increment in the Bond Order of
the adatoms with their periodic images (Figure 2e) and a decrement of the
total Bond Order with the atoms comprising the substrate (Figure 2f). This
implies that when the concentration of the substrate at the interface is
increased the species with larger atomic radii tend to form stabilising bonds
between their adatoms and decrease the interaction (form less bonds) with the
substrate. To further corroborate this observation we calculated the Pearson
correlation coefficients between the total Bond Order considering only
adatoms-adatoms bonds and the different adatom atomic radii (i.e.
B,C,F,H,N,O,P,S). The Paerson correlation coefficients are positive and above
0.2 for B(0.298), C(0.203), P(0.808) and S(0.625), indicating a strong
positive correlation between these quantities, hence confirming our
observation. For N the value is 0.115 indicating a somehow weaker correlation,
while for H(0.015),F(0.038) and O(0.043) the coefficients are all lower that
0.05 indicating almost no correlation between the strength of interactions
between adatoms (hence their coverage) and their Atomic radii.
v) Regarding the Bader and DDEC charge analysis (see Table 2 and 3 in the
Supporting Information, showing the net charge transferred to the adatom) it
can be observed that systems presenting higher adsorption energies also
present bigger charge transfers between the substrates and the adatoms. To
verify this trend for the whole database we calculated the Pearson correlation
coefficients between EAds and these two quantities (i.e. Bader and DDEC
cahrges). We obtain a coefficient of 0.64 for the Bader charges and 0.70 for
the DDEC charges, confirming the strong positive correlation between these
quantities. The charge transfers seems to be especially significant in the two
non-transition metal case and Titanium. It can be also observed, from the
Bader charge analysis, that for Boron and Phosporous in some limited cases the
electronic density is transferred from the adatom to the metal surface. This
behaviour reflect the fact that these two species are the ones with the lowest
electronegativity of the set we considered for the adatoms. This behaviour is
observed also using the DDEC charge analysis but only for Phosphorous
(although the Boron net charge values calculated are very near the zero). This
discrepancy is probably due to the different charge partitioning of the two
algorithms.

###  3.2 Adhesion Energy Analysis

Regarding the adhesion energy database we notice that:
i)The presence of adatoms leads in almost all cases to a reduction of the
Adhesion energy (see Figures 1B ). The important exceptions to this trend are
Boron and Carbon that at low $\Theta$ slightly increase the Adhesion Energy
and at higher coverage highly increase it. From Figure3a it can be seen that
for almost all cases (Mg, Zn and Al excluded) at high concentrations Boron and
Carbon are in fact the adatoms that result in the highest adhesion energy.
Both these species are known for their flexibility in forming bonds 44, 45,
46, 47, hence we can hypothesize that at a certain concentrations such
behavior could lead to the formation of stronger interactions at the interface
w.r.t the other adatoms. However, given the recent literature highlighting the
lubricant properties of graphene 48, 49 is important to specify how the
systems described in our simulations present a very different chemistry from
the graphene case. In fact in our system the carbon atoms are adsorbed at
distances where they cannot form stable bonds hence they are in a less stable
configuration and more prone to interact with the substrate and transfer
charge/electronic density. Our systems geometries (at least using the current
workflow setup) are more suited to be compared to the addition of additive in
alloys and metalworking. Several studies can be found where the addition of
Boron and Carbon lead to an increment of the material strength, toughness, and
high wear resistance50, 51, 52, 53, 54, 55. Many of these works focused on
Iron based systems50, 51 but even the other examples show that as long as the
presence of these elements is interstitial all the previously listed
properties are enhanced. In particular boron is nowadays commonly added to
many structural alloys primarily as a grain-boundary strengthener52, 53.
ii) Figure3b shows instead that the lowest Adhesion energies are obtained for
systems with Fluorine or Sulphur adsorbed at the interface. This trend reflect
the fact that sulfur56, 57 and fluorine compounds are at the basis of several
industrial lubricant additives. Several reports in literature also suggest F
as promising adhesion-reducing element at interfaces 58, 59, 60 and could be
also connected to the hydrophobicity of fluorinated interfaces (such as
Teflon). Our findings also confirm previous results for S that is know for its
lubricating properties, especially at Fe interfaces 61.
iii) It is interesting to observe how the difference in coverage plays a less
clear role for adhesion interactions w.r.t. the adsorption energies case (see
Figure3c). In fact for F, N, O and P there is on average an increment in
adhesion energy going from $\Theta=0.5$ to $\Theta=1.0$ even if the overall
EAdh value remains lower than in the clean interface case. On the contrary,
the adhesion energy is lowered if we increase the coverage passing from
$\Theta=0.25$ to $\Theta=0.5$ (exception made for B, C and N). The only two
adatom species for which the Adhesion decrease constantly with and increment
of the coverage are Sulfur and Hydrogen. Such trend has relevant technological
implications since hydrogen embitterment is a important challenge to resolve
in order to enable hydrogen to become a wide-spread fuel option 62.

iv) A possible explanation/hint of the reasons behind these multivariate
trends of the EAdh with respect to the coverage can be deduced by the charge
analysis of the different components of the systems simulated. Figure 3d-f
show the average change in the net Bader charge of the adatom (Panel D), down
slab (Panel E) and up slab (Panel F) w.r.t to the different change in surface
decoration. Regarding the adatom net charge it can be seen that, excluded
Hydrogen and Nitrogen with coverage passing from $\Theta=0.25$ to
$\Theta=0.5$, for all cases the decrement in coverage leads (on average) to a
gain of negative charge on the adatom. Consequentially, this accumulation of
electronic density on the adatoms generate a loss in negative charge on the up
and down slabs (i.e. the slab on which the adatoms are adsorbed). However, the
redistribution of the charge between the two slab is not symmetrical and leads
to non trivial electrostatic pattern at the interface. A characteristic
example is given by Fluorine where we always observe a loss in Adhesion energy
w.r.t. to the clean interface. However in many of these systems passing from
$\Theta=0.5$ to the monolayer leads to and increment in EAdh where going from
$\Theta=0.25$ to $\Theta=0.5$ results in the opposite effect. From Figure 3d-f
it can be seen that the F adatom loses on average 0.1 e- passing from coverage
0.25 to 0.5 and from 0.5 to 1. In response (on average) the down slab gains
electron density going from $\Theta=0.25$ to $\Theta=0.5$ but not going from
half coverage to the monolayer ( 0 e- gained). The up slab gains electron
density in both the increments in concentration but for the one passing from
$\Theta=0.5$ to $\Theta=1$ the average negative net charged obtained is bigger
(i.e. 0.07 e-). This results in a less negatively charged adatom monolayer/
less positively charged up-slab and hence in a diminished repulsive
electrostatic interaction. This smaller interaction probably leads to an
increment EAdh. The F case is just an example to highlight how the adhesion
energy is modulated case by case by the microscopic environment at the
interface. This observation further underlines the necessity forcareful
atomistic simulations of the interface environment/morphology.

Figure 1: Adsorption energy in eV (Top panel), adhesion energy in (J/m2)
(Middle Panel) and vertical distance along z between the adatom and the xy-
plane of the substrate slab in Å(Bottom Panel) for 25%, 50% and 100%
coverages. Empty squares indicates when the atomic adsorption is not stable.
In this cases, the adsorbate percolates into the surface structure. For the
adhesion energies the clean interfaces values calculated in Ref. 63 are
reported in the Column ”Clean”.In both cases darker shades are related to
higher adsorption and adhesion energies. for the adsorption systems at
coverage $\Theta$=0.25,0.5,1.0. Squares with 0 values indicates when the
atomic adsorption is not stable. In this cases, the adsorbate percolates into
the surface structure. Figure 2: Panel A: Maximum Adsorption Energy
(E${}_{A}bs$) for each substrate considered with respect to the coverage. The
adatom for which the Maximum Adsorption Energy is obtained is specified in the
label above the data-point. Panel B: Minimum Adsorption Energy (E${}_{A}bs$)
for each substrate considered with respect to the coverage. The adatom for
which the Maximum Adsorption Energy is obtained is specified in the label
above the data-point. Panel C: Average difference in Adsorption Energy between
different coverage for all the adatoms considered. Panel D: Average change in
distance from the substrate between different coverage for all the adatoms
considered. Panel E: Average change in Bond Order between adatoms considering
different coverage. Panel F: Average change in Bond Order between adatoms and
the substrate considering different coverage. Figure 3: Panel A: Maximum
Adhesion Energy (E${}_{A}dh$) for each substrate considered with respect to
the coverage. The adatom for which the Maximum Adsorption Energy is obtained
is specified in the label above the data-point. Panel B: Minimum Adhesion
Energy (E${}_{A}dh$) for each substrate considered with respect to the
coverage. The adatom for which the Maximum Adsorption Energy is obtained is
specified in the label above the data-point. Panel C: Average difference in
Adhesion Energy between different coverage for all the adatoms considered. The
dotted line represent the same analysis executed without considering the data
obtained for Mg and Al (i.e. non transition metal). Panel D: Average change in
the adatom charge (e-) between different coverage for all the adatoms
considered. Panel E: Average change in the down slab (i.e. the one with the
adsorbed adatom) charge (e-) between different coverage for all the adatoms
considered. Panel F: Average change in the up slab charge (e-) between
different coverage for all the adatoms considered. All the charge analysis
reported in this graphs are derived using Bader charges. The graphs obtained
using the DDEC charges can be found in the Supporting Info.

###  3.3 Interfacial density redistribution analysis

The Interfacial density redistribution should aid us in further explore the
electronic density redistribution at the interface. $\rho_{redist}$ shows
similar trends (see Image 4) to the adhesion energy one w.r.t to the Adatom
species in fact, for all the adatoms where the adhesion energy was drastically
reduced (P,S and F) we can observe how the $\rho_{redist}$ values are smaller.
This indicates a smaller redistribution of the charge density at boundary that
in previous work has been correlated to lower adhesion and hence friction 29.
$\rho_{redist}$ values generally remain similar at lower concentrations
($\Theta=0.25and0.5$) and than increase at full coverage. This behaviour could
be dictated by the fact that the higher the concentration of the adatoms the
bigger (on average) is the charge transferred with the two slabs (see
Figure3). The only two species that don’t follow this trend are H, N and to a
lesser extent O. In particular the incremental presence of hydrogen decreases
the electronic redistribution between metallic slabs. Almost all the values of
$\rho_{redist}$ calculated for the 25-50% concentrations are lower with
respect to their Clean reference. This in principle respect the correlation
observed between EAdh and $\rho_{redist}$ in previous work 29, 63. In fact the
lower the adhesion the lower the charge density redistribution. However for B
and C at $\Theta=0.5$ decoration the $\rho_{redist}$ value is smaller that the
clean interface one but Eadh is bigger. This deviation from the trend is even
more evident at full concentration where especially for d-metal with low
filled orbitals. These deviations support our previous observations for the
Adhesion Energies and further highlights how the presence, concentration
geometry of the adatoms has to be factored in the analysis interfacial charge
density redistribution. In fact it seems that the concentration of the adatoms
together with the charge transfer they induce and the consequent electrostatic
interactions at the nanoscale result in a non-trivial energetic picture with
respect to the Adhesion interaction.

Figure 4: Charge redistribution factor $\rho_{redist}$ in $e^{-}Ang^{-3}$ for
the adhesion systems at coverage $\Theta$=0.25,0.5,1.0. The clean interfaces
values calculated are reported in the Column ”Clean”.

##  4 Conclusions

In summary, we used an high-throughput approach to the study of the effects of
the presence of different species at several homogeneous metallic interfaces.
This systematic approach allowed us to uncover/confirm several interesting
trends:

  * •

The substrates with low filling of d-band are the most prone to adsorb ad-
atoms. This trend is consistent with d-band centre theory for adsorption on
metals.

  * •

For P, S and B increasing the coverage leads to a drastic reduction of the
adsorption energies that seems to be correlated to a widening of the distance
between the ad-atoms from the surface and the formation of bonds between
adatoms.

  * •

For H Eads varies much less with the change in coverage w.r.t all the other
species analysed.

  * •

Exception made for B and C the presence of non-metal adatoms always reduces
the adhesion. Our adhesion results were confirmed/supported by several
experimental results. Boron and Carbon are well known to be used in metal
casting and alloying to improve material strength, toughness, and high wear
resistance and Boron more specifically is used as a grain-boundary
strengthener. Such correspondence between experimental literature and our
results underlines the potential usefulness of our adhesion database to
engineer material with modular toughness, wear resistant properties.

  * •

P, S and in particular F are strong adhesion reducer. This trend was also
supported by several works we found in literature, however especially for P
and S the experimental studies were more focused on lubricant additives hence
their effects on the sliding behaviour at the interface should be further
explored.

  * •

Hydrogen present a unique behaviour respect to the other species probably due
to its small atomic radius. Its increment in concentration at the interface in
fact promote a clear reduction in adhesion.

  * •

The increment in coverage of the different adatoms lead to a non-trivial non-
symmetric charge transfer picture at the interface that can give rise to
different electrostatic regimes at the nanoscale.

  * •

This non-trivial picture is also reflected in the $\rho_{redist}$ values. The
relation between $\rho_{redist}$ and Eadh observed in previous works seems to
still be present. However the effects of the nature and concentration of the
adatoms give rise to several exception that will need to be further analysed.

This work represents a first step in the methodical study of more realistic
models of solid interfaces and their mechanical properties. Obviously the
systems simulated here are still far from realistic conditions, temperature,
pressure, heterogeneity, disorder and several other factors are missing.
Nonetheless our workflow approach can be extended (with relative ease) to
include some of these parameters and support further experimental
applications. Nonetheless we found several experimental work upholding our
findings. This aims to be the starting point to inform the experimental
community about functional interfaces in adhesion applications. In fact our
future goal is to further expand our workflow in order to consider a wider
range of concentrations materials and heterogeneous contacts so that a well
defined interface with pre-determined adhesion properties can be
engineered/obtained staring from from ab-initio data instead that just
empirical work.

##  5 Supporting Information Available

The Supporting Information is available free of charge. This contains the
databases of the d-band edges calculated for the absorption systems, the
database for the Bader charges of the Adatoms calculated for the absorption
systems and the database for the DDEC charges of the Adatoms calculated for
the absorption systems.

##  6 Competing interests

The authors declare no competing interests.

##  7 Acknowledgements

These results are part of the ”Advancing Solid Interface and Lubricants by
First Principles Material Design (SLIDE)” project that has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation program (Grant agreement No. 865633).

## References

  * Duine et al. 2018 Duine, R. A.; Lee, K.-J.; Parkin, S. S. P.; Stiles, M. D. Synthetic antiferromagnetic spintronics. _Nature Physics_ 2018, _14_ , 217
  * Tu et al. 2018 Tu, Z.; Choudhury, S.; Zachman, M. J.; Wei, S.; Zhang, K.; Kourkoutis,; F., L.; Archer, L. A. Fast ion transport at solid–solid interfaces in hybrid battery anodes. _Nature Energy_ 2018, _3_ , 310
  * Bhushan 2008 Bhushan, B. Nanotribology and nanomechanics in nano/biotechnology. _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 2008, _366_ , 1499–1537
  * Kendall 1988 Kendall, K. Theoretical aspects of solid-solid adhesion. _Science Progress (1933- )_ 1988, _72_ , 155–171
  * Krylov and Frenken 2014 Krylov, S. Y.; Frenken, J. W. M. The physics of atomic-scale friction: Basic considerations and open questions. _physica status solidi (b)_ 2014, _251_ , 711–736
  * Cihan et al. 2016 Cihan, E.; İpek, S.; Durgun, E.; Baykara, M. Z. Structural lubricity under ambient conditions. _Nature Communications_ 2016, _7_ , 12055
  * Medina and Dini 2014 Medina, S.; Dini, D. A numerical model for the deterministic analysis of adhesive rough contacts down to the nano-scale. _International Journal of Solids and Structures_ 2014, _51_ , 2620–2632
  * Müser et al. 2017 Müser, M. H. et al. Meeting the Contact-Mechanics Challenge. _Tribology Letters_ 2017, _65_ , 118
  * Siegel et al. 2002 Siegel, D. J.; Hector, L. G.; Adams, J. B. Adhesion, Atomic Structure, and Bonding at the $\mathrm{Al}(111)/\alpha-{\mathrm{Al}}_{2}{\mathrm{O}}_{3}(0001)$ Interface: A First Principles Study. _Phys. Rev. B_ 2002, _65_ , 085415
  * Feldbauer et al. 2018 Feldbauer, G.; Wolloch, M.; Bedolla, P. O.; Redinger, J.; Vernes, A.; Mohn, P. Suppression of Material Transfer at Contacting Surfaces: The Effect of Adsorbates on Al/TiN and Cu/Diamond Interfaces From First-Principles Calculations. _Journal of Physics: Condensed Matter_ 2018, _30_ , 105001
  * Tripathi et al. 2018 Tripathi, M.; Awaja, F.; Bizao, R. A.; Signetti, S.; Iacob, E.; Paolicelli, G.; Valeri, S.; Dalton, A.; Pugno, N. M. Friction and Adhesion of Different Structural Defects of Graphene. _ACS Applied Materials & Interfaces_ 2018, _10_ , 44614–44623
  * Cornil et al. 2020 Cornil, D.; Rivolta, N.; Mercier, V.; Wiame, H.; Beljonne, D.; Cornil, J. Enhanced Adhesion Energy at Oxide/Ag Interfaces for Low-Emissivity Glasses: Theoretical Insight Into Doping and Vacancy Effects. _ACS Applied Materials & Interfaces_ 2020, _12_ , 40838–40849
  * Wang et al. 2021 Wang, Q.; Li, Y.; Chen, Z.; Wang, M.; Zhu, H.; Wang, H. Understanding Alloying Behaviors of Sc, Ni and Zn Additions on Al/TiB2 Interfaces Based on Interfacial Characteristics and Solute Properties. _Surfaces and Interfaces_ 2021, _26_ , 101427
  * Rabinowicz 1971 Rabinowicz, E. The Determination of the Compatibility of Metals Through Static Friction Tests. _A S L E Transactions_ 1971, _14_ , 198–205
  * Matějíček et al. 2013 Matějíček, J.; Vilémová, M.; Mušálek, R.; Sachr, P.; Horník, J. The Influence of Interface Characteristics on the Adhesion/Cohesion of Plasma Sprayed Tungsten Coatings. _Coatings_ 2013, _3_ , 108–125
  * Heuer et al. 2019 Heuer, S.; Matějíček, J.; Vilémová, M.; Koller, M.; Illkova, K.; Veverka, J.; Weber, T.; Pintsuk, G.; Coenen, J.; Linsmeier, C. Atmospheric Plasma Spraying of Functionally Graded Steel/Tungsten Layers for the First Wall of Future Fusion Reactors. _Surface and Coatings Technology_ 2019, _366_ , 170–178
  * Heuer et al. 2019 Heuer, S.; Lienig, T.; Mohr, A.; Weber, T.; Pintsuk, G.; Coenen, J.; Gormann, F.; Theisen, W.; Linsmeier, C. Ultra-Fast Sintered Functionally Graded Fe/W Composites for the First Wall of Future Fusion Reactors. _Composites Part B: Engineering_ 2019, _164_ , 205–214
  * Ding et al. 2008 Ding, F.; Larsson, P.; Larsson, J. A.; Ahuja, R.; Duan, H.; Rosén, A.; Bolton, K. The Importance of Strong Carbon Metal Adhesion for Catalytic Nucleation of Single-Walled Carbon Nanotubes. _Nano Letters_ 2008, _8_ , 463–468, PMID: 18162001
  * Hafsi et al. 2018 Hafsi, Z.; Mishra, M.; Elaoud, S. Hydrogen embrittlement of steel pipelines during transients. _Procedia Structural Integrity_ 2018, _13_ , 210–217, ECF22 - Loading and Environmental effects on Structural Integrity
  * Louthan and Caskey 1976 Louthan, M.; Caskey, G. Hydrogen transport and embrittlement in structural metals. _International Journal of Hydrogen Energy_ 1976, _1_ , 291–305
  * Ayyagari et al. 2022 Ayyagari, A.; Alam, K. I.; Berman, D.; Erdemir, A. Progress in Superlubricity Across Different Media and Material Systems—A Review. _Frontiers in Mechanical Engineering_ 2022, _8_
  * Pang et al. 2020 Pang, L.; Barras, A.; Mishyn, V.; Sandu, G.; Melinte, S.; Subramanian, P.; Boukherroub, R.; Szunerits, S. Enhanced electrocatalytic hydrogen evolution on a plasmonic electrode: the importance of the Ti/TiO2 adhesion layer. _J. Mater. Chem. A_ 2020, _8_ , 13980–13986
  * Curtarolo et al. 2013 Curtarolo, S.; Hart, G. L. W.; Nardelli, M. B.; Mingo, N.; Sanvito, S.; Levy, O. The high-throughput highway to computational materials design. _Nature Materials_ 2013, _12_ , 191–201
  * Haastrup et al. 2018 Haastrup, S.; Strange, M.; Pandey, M.; Deilmann, T.; Schmidt, P. S.; Hinsche, N. F.; Gjerding, M. N.; Torelli, D.; Larsen, P. M.; Riis-Jensen, A. C.; Gath, J.; Jacobsen, K. W.; Mortensen, J. J.; Olsen, T.; Thygesen, K. S. The Computational 2D Materials Database: high-throughput modeling and discovery of atomically thin crystals. _2D Materials_ 2018, _5_ , 042002
  * Li et al. 2022 Li, Z.; Yoon, J.; Zhang, R.; Rajabipour, F.; Srubar III, W. V.; Dabo, I.; Radlińska, A. Machine learning in concrete science: applications, challenges, and best practices. _npj Computational Materials_ 2022, _8_ , 127
  * Rosen et al. 2022 Rosen, A. S.; Fung, V.; Huck, P.; O’Donnell, C. T.; Horton, M. K.; Truhlar, D. G.; Persson, K. A.; Notestein, J. M.; Snurr, R. Q. High-throughput predictions of metal–organic framework electronic properties: theoretical challenges, graph neural networks, and data exploration. _npj Computational Materials_ 2022, _8_ , 112
  * Hebnes et al. 2022 Hebnes, O. L.; Bathen, M. E.; Schøyen, Ø. S.; Winther-Larsen, S. G.; Vines, L.; Hjorth-Jensen, M. Predicting solid state material platforms for quantum technologies. _npj Computational Materials_ 2022, _8_ , 207
  * Losi et al. 2023 Losi, G.; Chehaimi, O.; Righi, M. C. To be published
  * Wolloch et al. 2018 Wolloch, M.; Levita, G.; Restuccia, P.; Righi, M. C. Interfacial Charge Density and Its Connection to Adhesion and Frictional Forces. _Phys. Rev. Lett._ 2018, _121_ , 026804
  * Jain et al. 2015 Jain, A.; Ong, S. P.; Chen, W.; Medasani, B.; Qu, X.; Kocher, M.; Brafman, M.; Petretto, G.; Rignanese, G.-M.; Hautier, G.; Gunter, D.; Persson, K. A. FireWorks: a dynamic workflow system designed for high-throughput applications. _Concurrency and Computation: Practice and Experience_ 2015, _27_ , 5037–5059
  * Mathew et al. 2017 Mathew, K. et al. Atomate: A high-level interface to generate, execute, and analyze computational materials science workflows. _Computational Materials Science_ 2017, _139_ , 140–152
  * Ong et al. 2013 Ong, S. P.; Richards, W. D.; Jain, A.; Hautier, G.; Kocher, M.; Cholia, S.; Gunter, D.; Chevrier, V. L.; Persson, K. A.; Ceder, G. Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis. _Computational Materials Science_ 2013, _68_ , 314–319
  * Mathew et al. 2016 Mathew, K.; Singh, A. K.; Gabriel, J. J.; Choudhary, K.; Sinnott, S. B.; Davydov, A. V.; Tavazza, F.; Hennig, R. G. MPInterfaces: A Materials Project based Python tool for high-throughput computational screening of interfacial systems. _Computational Materials Science_ 2016, _122_ , 183–190
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 1996, _54_ , 11169–11186
  * Kresse and Joubert 1999 Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. _Phys. Rev. B_ 1999, _59_ , 1758–1775
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Yu and Trinkle 2011 Yu, M.; Trinkle, D. R. Accurate and efficient algorithm for Bader charge integration. _The Journal of Chemical Physics_ 2011, _134_ , 064111
  * Limas and Manz 2018 Limas, N. G.; Manz, T. A. Introducing DDEC6 atomic population analysis: part 4. Efficient parallel computation of net atomic charges, atomic spin moments, bond orders, and more. _RSC Adv._ 2018, _8_ , 2678–2707
  * Nørskov et al. 2011 Nørskov, J. K.; Abild-Pedersen, F.; Studt, F.; Bligaard, T. Density functional theory in surface chemistry and catalysis. _Proceedings of the National Academy of Sciences_ 2011, _108_ , 937–943
  * Wang et al. 2021 Wang, V.; Xu, N.; Liu, J.-C.; Tang, G.; Geng, W.-T. VASPKIT: A user-friendly interface facilitating high-throughput computing and analysis using VASP code. _Computer Physics Communications_ 2021, _267_ , 108033
  * Manz 2017 Manz, T. A. Introducing DDEC6 atomic population analysis: part 3. Comprehensive method to compute bond orders. _RSC Adv._ 2017, _7_ , 45552–45581
  * Hammer and Nørskov 1995 Hammer, B.; Nørskov, J. Electronic factors determining the reactivity of metal surfaces. _Surface Science_ 1995, _343_ , 211–220
  * Hammer and Nørskov 1995 Hammer, B.; Nørskov, J. K. Why gold is the noblest of all the metals. _Nature_ 1995, _376_ , 238–240
  * Bhunia 2020 Bhunia, A. Boron Has More Surprises. _Chem_ 2020, _6_ , 320–321
  * Qiu et al. 2023 Qiu, L.; Zhang, X.; Kong, X.; Mitchell, I.; Yan, T.; Kim, S. Y.; Yakobson, B. I.; Ding, F. Theory of sigma bond resonance in flat boron materials. _Nature Communications_ 2023, _14_ , 1804
  * Nyholm and Espallargas 2023 Nyholm, N.; Espallargas, N. Functionalized carbon nanostructures as lubricant additives – A review. _Carbon_ 2023, _201_ , 1200–1228
  * Zhang et al. 2020 Zhang, K.; Zhang, Y.; Shi, L. A review of linear carbon chains. _Chinese Chemical Letters_ 2020, _31_ , 1746–1756
  * Li et al. 2018 Li, Z.; Xu, C.; Xiao, G.; Zhang, J.; Chen, Z.; Yi, M. Lubrication Performance of Graphene as Lubricant Additive in 4 n pentyl 4 cyanobiphyl Liquid Crystal (5CB) for Steel/Steel Contacts. _Materials_ 2018, _11_
  * Berman et al. 2014 Berman, D.; Erdemir, A.; Sumant, A. V. Graphene: a new emerging lubricant. _Materials Today_ 2014, _17_ , 31–42
  * Izaki et al. 2002 Izaki, M.; Miyamoto, N.; Nakae, A.; Hasegawa, T.; Watase, S.; Chigane, M.; Fujiwara, Y.; Enomoto, M. I. H. Martensitic Iron-Carbon-Boron Alloy Electrodeposit with Improved Mechanical Properties. _Journal of The Electrochemical Society_ 2002, _149_ , C370
  * li LIU et al. 2009 li LIU, Z.; CHEN, X.; xiang LI, Y.; hua HU, K. High Boron Iron-Based Alloy and Its Modification. _Journal of Iron and Steel Research, International_ 2009, _16_ , 37–54
  * Sharma et al. 2022 Sharma, A.; Dixit, S.; Baler, N.; Agrawal, P.; Makineni, S. K.; Chattopadhyay, K. Impact of boron as an alloying addition on the microstructure, thermo-physical properties and creep resistance of a tungsten-free Co-base superalloy. _Materials Science and Engineering: A_ 2022, _855_ , 143899
  * Ma and Zhang 2016 Ma, S.; Zhang, J. WEAR RESISTANT HIGH BORON CAST ALLOY - A REVIEW. 2016
  * Garcia-Cabezon et al. 2020 Garcia-Cabezon, C.; Garcia-Hernandez, C.; Rodriguez-Mendez, M. L.; Herranz, G.; Martin-Pedrosa, F. Role of carbon and nitrogen in the improvement of corrosion resistanceof new powder metallurgy Co-Cr-Mo alloys. _Corrosion Reviews_ 2020, _38_ , 273–286
  * Szkliniarz and Szkliniarz 2023 Szkliniarz, A.; Szkliniarz, W. Carbon in Commercially Pure Titanium. _Materials_ 2023, _16_
  * CHEN et al. 2023 CHEN, Y.; RENNER, P.; LIANG, H. A review of current understanding in tribochemical reactions involving lubricant additives. _Friction_ 2023, _11_ , 489–512
  * Mohamed Faruck et al. 2020 Mohamed Faruck, A. A.; Hsu, C.-J.; Doerr, N.; Weigand, M.; Gachot, C. How lubricant formulations and properties influence the performance of rotorcraft transmissions under loss of lubrication conditions. _Tribology International_ 2020, _151_ , 106390
  * Somers et al. 2012 Somers, A. E.; Biddulph, S. M.; Howlett, P. C.; Sun, J.; MacFarlane, D. R.; Forsyth, M. A comparison of phosphorus and fluorine containing IL lubricants for steel on aluminium. _Phys. Chem. Chem. Phys._ 2012, _14_ , 8224–8231
  * Veldhuis et al. 2007 Veldhuis, S.; Dosbaeva, G.; Benga, G. Application of ultra-thin fluorine-content lubricating films to reduce tool/workpiece adhesive interaction during thread-cutting operations. _International Journal of Machine Tools and Manufacture_ 2007, _47_ , 521–528
  * Ware et al. 2018 Ware, C. S.; Smith-Palmer, T.; Peppou-Chapman, S.; Scarratt, L. R. J.; Humphries, E. M.; Balzer, D.; Neto, C. Marine Antifouling Behavior of Lubricant-Infused Nanowrinkled Polymeric Surfaces. _ACS Applied Materials & Interfaces_ 2018, _10_ , 4173–4182, PMID: 29250952
  * Fatti and Righi 2019 Fatti, G.; Righi, M. C. Selenium Chemisorption Makes Iron Surfaces Slippery. _Tribology Letters_ 2019, _7_ , 125
  * Dwivedi and Vishwakarma 2018 Dwivedi, S. K.; Vishwakarma, M. Hydrogen embrittlement in different materials: A review. _International Journal of Hydrogen Energy_ 2018, _43_ , 21603–21616
  * Restuccia et al. 2023 Restuccia, P.; Losi, G.; Chehaimi, O.; Marsili, M.; Righi, M. C. High-Throughput First-Principles Prediction of Interfacial Adhesion Energies in Metal-on-Metal Contacts. _ACS Applied Materials & Interfaces_ 2023, _15_ , 19624–19633, PMID: 37015021


[orcid=0000-0003-0430-7256]

[1]

[1]

url]www.ariunbold.physics.msstate.edu

1]organization=Department of Physics and Astronomy, addressline=Mississippi
State University, city=Starkville, postcode=39762, state=MS, country=USA

[1]Corresponding author

# Quantum dynamics via a hidden Liouville space

Gombojav O. Ariunbold<EMAIL_ADDRESS>[

###### Abstract

Quantum dynamics for arbitrary system are traditionally realized by time
evolutions of wave functions in Hilbert space and/or density operators in
Liouville space. However, the traditional simulations may occasionally turn
out to be challenging for the quantum dynamics, particularly those governed by
the nonlinear Hamiltonians. In this letter, we introduce a nonstandard
iterative technique where time interval is divided into a large number of
discrete subintervals with an ultrashort duration; and the Liouville space is
briefly expanded with an additional (virtual) space only within these
subintervals. We choose two-state spin raising and lowering operators for
virtual space operators because of their simple algebra. This tremendously
reduces the cost of time-consuming calculations. We implement our technique
for an example of a charged particle in both harmonic and anharmonic
potentials. The temporal evolutions of the probability for the particle being
in the ground state are obtained numerically and compared to the analytical
solutions. We further discuss the physics insight of this technique based on a
thought-experiment. Successive processes intrinsically ’hitchhiking’ via
virtual space in discrete ultrashort time duration, are the hallmark of our
simple iterative technique. We believe that this novel technique has potential
for solving numerous problems which often pose a challenge when using the
traditional approach based on time-ordered exponentials.

###### keywords:

Hilbert space Liouville space harmonic oscillator anharmonic oscillator
intensity-dependent oscillator Schr${\rm\ddot{o}}$dinger representation von-
Neumann equation S-operator superradiance sub-radiance time-ordered
exponentials

A new iterative density-operator technique is introduced.

Quantum dynamics for the driven generalized quantum oscillators are studied
both numerically and analytically.

##  1 Introduction

Although, the standard approach based on time-ordered exponentials is
extremely useful [1, 2, 3], it may occasionally turn out to be challenging,
particularly, in the case of revealing nonlinear quantum dynamics [4, 5] that
requires rigorous numerical simulations [6, 7, 8]. Quantum dynamics for
arbitrary system are traditionally realized by time evolutions of wave
functions in Hilbert space, which can also be expressed in terms of density
operators in the Liouville space [2, 3]. In this letter, we introduce a new
nonstandard iterative technique formulated as follows. (i) Finite time
interval is divided into a large number of discrete subintervals with an
ultrashort width. (ii) The Liouville space is expanded with an additional
(i.e., virtual) space for this ultrashort time duration. The system’s original
Hamiltonian is, then, modified for the system’s space plus virtual space,
where the force terms are replaced with the virtual quantum operators. (iii)
The density operator for the system is extracted by tracing over the virtual
operator space. In principle, various virtual operators can be chosen
depending on the specific quantum system. Here we choose two-state spin
raising and lowering operators because of their simple algebra. In the next
section, we present the standard approach using S-operator defined as time-
ordered exponentials in Hilbert, and then, in the Liouville space. In section
3, we introduce our technique and implement it to the well-known example of a
charged particle in a harmonic potential. The temporal evolutions of the
probability for the particle being in the ground state are obtained by our
technique and compared to the analytical solutions obtained using the standard
S-operator. By extending this example, we perform numerical simulations for
temporal evolutions for the ground state probability for the generalized
systems governed by time-dependent nonlinear Hamiltonians. We further discuss
the physics insight of this technique based on a thought-experiment, in which
a large number of polarized atoms successively interact with a lossless cavity
field. The last section is a conclusion.

##  2 Standard approach

In this section, the standard approach for quantum dynamics both in Hilbert
space and the Liouville space is presented. We consider the system with the
Hamiltonian given by

$\hat{H}=\hat{H}_{0}+\hat{V}$ (1)

here $\hat{H}_{0}$ is the unperturbed (free) and $\hat{V}$ interaction
Hamiltonians and we set $\hbar\equiv 1$.

###  2.1 Quantum dynamics in Hilbert space

We begin with the approach for the Hilbert space. In the interaction
representation, the rapid state evolution due to $\hat{H}_{0}$ is removed as
$|\psi_{I}(t)\rangle={\rm exp}{(i\hat{H}_{0}t)}|\psi_{S}(t)\rangle$, where
$|\psi_{I}(t)\rangle$ and $|\psi_{S}(t)\rangle$ are wave functions in the
interaction and Schr${\rm\ddot{o}}$dinger representations, respectively.
Unitary transformation of initial state in the interaction picture is given as
$|\psi_{I}(t)\rangle=\hat{U}(t)|\psi_{I}(0)\rangle$, here unitary operator
$\hat{U}(t)$ satisfies $\hat{U}(t)^{\dagger}\hat{U}(t)=\hat{1}$ and is
expressed as $\hat{U}(t)={\rm exp}{(i\hat{H}_{0}t}){\rm exp}{(-i\hat{H}t)}$.
Time evolution of $\hat{U}(t)$ can be derived from
$i{\partial\hat{U}}/{\partial t}=\hat{V}_{I}(t)\hat{U}(t)$. The Hamiltonian is
in the interaction representation as $\hat{V}_{I}=\hat{V}_{I}^{\dagger}={\rm
exp}{(i\hat{H}_{0}t)}\hat{V}{\rm exp}{(-i\hat{H}_{0}t)}$. Choosing time
interval between $t_{2}$ and $t_{1}$ $(t_{1}<t_{2})$, unitary transformation
is expressed as [1]

$|\psi_{I}(t_{2})\rangle=\hat{S}(t_{2},t_{1})|\psi_{I}(t_{1})\rangle$ (2)

with S-operator $\hat{S}(t_{2},t_{1})$. We divide time interval $t_{2}-t_{1}$
into $N$ sub-intervals with a width of $\Delta t$. At mid-time
$\tau_{j}=t_{1}+(j-1/2)\Delta t$ in the $j$th interval, the S-operator is
written as [1]

$\hat{S}\left(\tau_{j}+{\Delta t}/{2},\tau_{j}-{\Delta t}/{2}\right)={\rm
e}^{-i\hat{V}_{I}(\tau_{j})\Delta t}$ (3)

where $N\rightarrow\infty$ and $\Delta t\rightarrow 0$ but $t_{2}-t_{1}$ is
finite. Eq.(3) leads to the traditional time-ordered exponential given as [1]

$\hat{S}\left(t_{2},t_{1}\right)={\hat{T}}\left[\prod_{j=1}^{N}{\rm
e}^{-i\hat{V}_{I}(\tau_{j})\Delta t}\right]={\hat{T}}\left[{\rm
e}^{-i\int_{t_{1}}^{t_{2}}\hat{V}_{I}(t)dt}\right]$ (4)

where time ordering for boson operators is defined as
${\hat{T}}[\hat{O}_{1}(\tau_{1})$ $\hat{O}_{2}(\tau_{2})$ $\cdots$
$\hat{O}_{N}(t_{N})]$ $=$ $\hat{O}_{p_{1}}(\tau_{p_{1}})$
$\hat{O}_{p_{2}}(\tau_{p_{2}})$ $\cdots$ $\hat{O}_{p_{N}}(t_{p_{N}})$ with
$\tau_{p_{1}}>\tau_{p_{2}}\cdots>\tau_{p_{N}}$.

###  2.2 Examples

As an example, we consider a driven harmonic oscillator. For that temporal
evolutions of the ground state using S-operator are given in Eq.(4). Let a
particle of a charge $q\equiv 1$, mass $m\equiv 1/2$ be in a harmonic
potential $(\hbar\equiv 1)$. The driving electric field is $E(t)\equiv 1$, if
$T>t>0$, and otherwise, it is zero and $\omega$ is the frequency of the
oscillator. In the interaction picture, the Hamiltonian is written as

$\hat{V}_{I}(t)=\hat{R}(t)\epsilon^{\ast}(t)+\hat{R}^{\dagger}(t)\epsilon(t)$
(5)

where time-dependent operators are $\hat{R}(t)=\hat{b}(t)$ $=$ $\hat{b}{\rm
e}^{-i\omega t}$ and $\hat{R}^{\dagger}(t)$ $=$ $\hat{b}^{\dagger}(t)$ $=$
$\hat{b}^{\dagger}{\rm e}^{i\omega t}$ and the force terms are $\epsilon(t)$
$=$ $\epsilon^{\ast}(t)$ $=$ $-\sqrt{{1}/{\omega}}E(t)$. The probability
$p(T)$ for the particle to remain in the ground state $|\psi_{0}\rangle$ after
time $T$ is written as

$p(T)=|\langle\psi_{0}|\hat{S}(T,0)|\psi_{0}\rangle|^{2}$ (6)

The probability amplitude is given by S-operator from Eq.(4) as

$\langle\psi_{0}|\hat{S}(T,0)|\psi_{0}\rangle=\langle\psi_{0}|{\hat{T}}{\rm
e}^{-i\int_{0}^{T}\hat{V}_{I}(t^{\prime})dt^{\prime}}|\psi_{0}\rangle={\rm
e}^{-iB(T)}$ (7)

where
$B(T)=\int_{0}^{T}dtdt^{\prime}\epsilon(t)^{\ast}G(t-t^{\prime})\epsilon(t^{\prime})$
and $G(t)$ is the Green’s function. For this simple example, the Green’s
function is well known

$G(t)=-i{\rm e}^{-i\omega t}\theta(t).$ (8)

Therefore, for the particle, its probability to remain in the ground state
after time $T$ is analytically found to be as [1]

$p(T)={\rm exp}\left[-\frac{4}{\omega^{2}}{\rm sin}^{2}\left(\frac{\omega
T}{2}\right)\right]$ (9)

with pulse area $\omega T$.

This example is the simplest case when a linear Hamiltonian is considered.
That conveniently ensures to use the well known Green function in Eq.(8).
However, in general, the Green’s functions are mostly unknown and a laborious
numerical method is often needed. Next, we consider two more examples that use
nonlinear Hamiltonians. The first example for the nonlinear Hamiltonian is a
driven anharmonic oscillator. The Hamiltonian is given in the form in Eq.(5)
[2] where degenerate two-boson nonlinear operators $\hat{R}(t)$ $=$
$\hat{b}^{2}(t)$ $=$ $\hat{b}^{2}{\rm e}^{-2i\omega t}$ and
$\hat{R}^{\dagger}(t)$ $=$ $\hat{b}^{2}{}^{\dagger}(t)$ $=$
$\hat{b}^{2}{}^{\dagger}{\rm e}^{2i\omega t}$ and the force terms
$\epsilon(t)$ $=$ $\epsilon^{\ast}(t)$ $=$ $-\sqrt{{1}/{2\omega}}E(t)$ are
assumed to have a similar form as in the harmonic oscillator case. The second
example for the nonlinear Hamiltonian is a driven intensity-dependent
oscillator. The Hamiltonian is given in the form in Eq.(5) [9] with the
intensity-dependent nonlinear boson operators $\hat{R}(t)$ $=$
$\hat{b}(t)\sqrt{\hat{b}^{\dagger}\hat{b}}$ and $\hat{R}^{\dagger}(t)$ $=$
$\sqrt{\hat{b}^{\dagger}\hat{b}}\hat{b}^{\dagger}(t)$ and $\epsilon(t)$ is
assumed to be the same as before. In section 3, we numerically solve for the
probability time evolutions for these nonlinear systems and compare with the
approximate analytical results.

##  3 Nonstandard approach

###  3.1 Quantum dynamics in Liouville space

Before introducing our technique, let us first replace the wave functions in
Hilbert space with density operators in Liouville space [2, 3]. We recall that
$\rho_{N}\equiv\rho(t_{2})=|\psi_{I}(t_{2})\rangle\langle\psi_{I}(t_{2})|$ and
$\rho_{0}\equiv\rho(t_{1})=|\psi_{I}(t_{1})\rangle\langle\psi_{I}(t_{1})|$
from Eq.(2). Using S-operator in the $j$th interval from Eq.(3), we rewrite
Eq.(2) in terms of density operators rather than wave functions as

$\hat{\rho}_{j}={\rm e}^{-i\hat{V}_{I}(\tau_{j})\Delta t}\hat{\rho}_{j-1}{\rm
e}^{i\hat{V}_{I}(\tau_{j})\Delta t}$ (10)

here $\hat{\rho}_{j}=\hat{\rho}(\tau_{j}+{\Delta t}/{2})$ and
$\hat{\rho}_{j-1}=\hat{\rho}(\tau_{j}-{\Delta t}/{2})$. In the traditional
approach, to obtain $\hat{\rho}(t_{2})$ at later time $t_{2}$ for any given
initial state $\hat{\rho}(t_{1})$ at $t_{1}$ Eq.(10) is repeatedly evaluated,
where $\Delta t\ll 1$ and $N\gg 1$ but $t_{2}-t_{1}=N\Delta t$ is finite. It
is also important to note that Eq.(10) is the formal solution of the Liouville
- von Neumann equation [2].

###  3.2 Quantum dynamics in Liouville space restructured with a virtual space

From this point, we implement our new nonstandard approach, rather directly
evaluating Eq.(10). As before, finite time interval $t_{2}-t_{1}$ is divided
into $N$ discrete subintervals with an ultrashort width of $\Delta t$. The
Liouville space is expanded with a two-state spin operator space for duration
of $\Delta t$. The system’s original Hamiltonian is, then, modified for the
system’s space plus spin space, where the force terms are replaced with the
spin operators. The density operator for the system is extracted by tracing
over the spin operator space. In the $j$th interval with an infinitesimally
short width of $\Delta t$, it is an acceptable ansatz where we replace the
original Hamiltonian $\hat{V}_{I}(\tau_{j})$ by new Hamiltonian expanding it
with an additional virtual space $\hat{A}_{j}$ as

$\hat{V}_{I}(\tau_{j})\rightarrow\hat{V}_{I}(\tau_{j})\otimes\hat{A}_{j}$ (11)

where $\left[\hat{V}_{I}(\tau_{j}),\hat{A}_{j}\right]=\hat{0}$. For the sake
of simplicity, $\hat{A}_{j}$ can be chosen to be a two-state spin operator
defined as

$\displaystyle\hat{A}_{j}$ $\displaystyle=$
$\displaystyle|\alpha_{j}|^{2}|\uparrow\rangle\langle\uparrow|+\alpha_{j}^{\ast}\beta_{j}|\downarrow\rangle\langle\uparrow|$
(12) $\displaystyle+$
$\displaystyle\alpha_{j}\beta_{j}^{\ast}|\uparrow\rangle\langle\downarrow|+|\beta_{j}|^{2}|\downarrow\rangle\langle\downarrow|$

with $|\alpha_{j}|^{2}+|\beta_{j}|^{2}=1$. For this choice, the force terms
are replaced with the raising $|\uparrow\rangle\langle\downarrow|$ and
lowering $|\downarrow\rangle\langle\uparrow|$ operators specifically as

$\displaystyle\epsilon(\tau_{j}){\rm
e}^{i\omega\tau_{j}}\rightarrow\eta_{j}\alpha_{j}\beta_{j}^{\ast}|\downarrow\rangle\langle\uparrow|$
$\displaystyle\epsilon^{\ast}(\tau_{j}){\rm
e}^{-i\omega\tau_{j}}\rightarrow\eta_{j}^{\ast}\alpha_{j}^{\ast}\beta_{j}|\uparrow\rangle\langle\downarrow|$
(13)

Thus, in the $j$th interval with $\Delta t$ width, this original Hamiltonian
Eq.(5) can be replaced with a new Hamiltonian $\hat{V}_{A}(\tau_{j})$, also
known as the unified Jaynes-Cummings Hamiltonian [11, 12] as

$\hat{V}_{A}(\tau_{j})=\eta_{j}^{\ast}\hat{R}(\tau_{j})|\uparrow\rangle\langle\downarrow|+\eta_{j}\hat{R}^{\dagger}(\tau_{j})|\downarrow\rangle\langle\uparrow|.$
(14)

Instead of the original approach given by Eq.(10) for $\hat{\rho}_{j}$, we
introduce an iterative relation for new density operator $\varrho_{j}$ using
the modified Hamiltonian given in Eq.(14) as

$\hat{\varrho}_{j}={\rm Tr}_{A}\left[{\rm e}^{-i\hat{V}_{A}(\tau_{j})\Delta
t}{\hat{A}_{j}}\otimes\hat{\varrho}_{j-1}{\rm e}^{i\hat{V}_{A}(\tau_{j})\Delta
t}\right]$ (15)

Therefore, our goal is to demonstrate that the two density operators converge

$\hat{\varrho}_{N}\simeq\hat{\rho}_{N}$ (16)

for the same pure initial state.

###  3.3 A thought-experiment

The essentials of our iterative technique are explained by the following
thought-experiment.

Figure 1: A monokinetic beam of $N$ two-level atoms passing through a lossless
cavity. Each atom interacts with the cavity field for a duration of $\Delta
t$. The cavity field builds up to its final state $\hat{\varrho}_{N}$ (or
$|\psi_{I}(t)\rangle$) from initial state $\hat{\varrho}_{0}$ (or
$|\psi_{I}(0)\rangle$) after a finite time $t=N\Delta t$.

As sketched in Fig. 1, let us consider a monokinetic beam consisting of
individual two-level atoms. Each atom is prepared in arbitrary coherent
superposition of the upper and lower states [13, 14, 8]. The atoms are then
injected into a lossless cavity in a well controlled rate where only one atom
at a time is present inside the cavity for duration $\Delta t$. At the exit
from the cavity the individual atoms are not intended to be measured. Total
number of atoms is $N$ and the $j$th atom-field coupling constant is
$\eta_{j}$. Although the present model can be generalized to multi-level atoms
[10, 9, 8], for the sake of simplicity, we consider only two-level atoms,
where $|\uparrow\rangle$ and $|\downarrow\rangle$ are upper and lower atomic
states, respectively. Correspondingly, $\alpha_{j}$ and $\beta_{j}$ are
probability amplitudes for the $j$th atomic upper and lower states. Thus, as a
result of numerically solving Eq.(15), the final cavity field state is
evaluated from the existing initial quantum state in the cavity after time
$t=\Delta tN$. For example, when atoms are prepared in the same phase then the
cavity field evolves to the so-called superradiant state [5, 8, 17, 18, 15,
16]. The mean number of photons created in the cavity (i.e., field intensity)
is proportional to $N^{2}$ rather than $N$. On the other hand, when each
successive pair of atoms are prepared in perfectly out-of-phase, then the
cavity field evolves to the sub-radiant state [5, 8]. Moreover, we justify
that time evolutions involve pure states after tracing over the virtual space
operator. As demonstrated in our earlier work [8], an initial coherent state
given as $|\gamma_{0}\rangle$ evolves into
$|\gamma(t)\rangle\simeq|\gamma_{A}+\gamma_{0}\rangle$, with
$\gamma_{A}=-i\eta\alpha\beta^{\ast}t$ at later time $t$. Therefore, the above
statement that our technique maintains time evolutions for pure states is
justified not only for infinitesimally short $\Delta t$ interval, but also for
finite time $t$.

###  3.4 Examples

Here we employ our technique for the previous example for a charged particle
in harmonic potential. We evaluate Eq.(15) using the Hamiltonian given in
Eq.(14) both numerically and analytically, however, for the sake of
simplicity, only for initial vacuum state. Thus, the modified Hamiltonian
$\hat{V}_{A}(\tau_{j})$ in Eq.(14) is written in terms of
$\hat{R}(t)\equiv\hat{b}(t)$. Comparing the Hamiltonian in Eq.(5) with the
ansatz in Eq.(13), we obtain $\epsilon(\tau_{j}){\rm e}^{i\omega\tau_{j}}$ $=$
$-\eta_{j}\alpha_{j}\beta_{j}^{\ast}$ $=$ $-\eta_{j}\zeta_{j}^{\ast}$. The
parameter $\zeta_{j}=\alpha_{j}^{\ast}\beta_{j}$ stands for a coherence
between spin states. For example, for parameters chosen to be as $\eta_{j}=1$,
$|\alpha_{j}|=|\beta_{j}|=1/\sqrt{2}$, it is given by $-{\rm
e}^{i\omega\tau_{j}}/2$ with $|\zeta_{j}|=1/2$. In the Fock state
representation, the $j$th density matrix elements are
$\rho_{j}(n,n^{\prime})=\langle n|\hat{\rho}_{j}|n^{\prime}\rangle$. We
numerically evaluate $\varrho_{N}(n,n^{\prime})$ from Eq.(15), to obtain
$\varrho_{N}(0,0)$, at $t_{2}-t_{1}=T$ with $t_{1}=0$ to compare the
probability $p(T)$ given in Eq.(9). Eq.(9) is rewritten in terms of time-
independent parameters $\eta_{j}=1$, $|\zeta_{j}|=|\zeta|$ associated to the
virtual operator space as

$p(T)={\rm exp}\left[-\frac{4|\zeta|^{2}}{\omega^{2}}{\rm
sin}^{2}\left(\frac{\omega T}{2}\right)\right]$ (17)

Figure 2: Analytical (black curves) and numerical (red circles) results for
temporal evolutions of the probability $p(t)$ for a charged particle driving
by the external field being in the ground state as functions of $\omega
t/\pi$. Left column: The system with a linear Hamiltonian. Middle column: The
system with the degenerate two-boson Hamiltonian. Right column: The system
with the intensity-dependent Hamiltonian.

In Fig. 2, the density matrix elements for $\varrho_{N}(0,0)$ (red circles)
and $p(t)$ (black curves) are plotted as functions of pulse area per $\pi$,
$\omega t/\pi$. For the plots in Figs. 2(a) and 3(a), the parameters include
coherence $\zeta=1/2$, total number $N=3750$, width of the subintervals
$\Delta t=0.001$, time $T=N\Delta t=3.75$ frequency $\omega=2\pi/5$. For the
plot in Fig. 2(a’), except for the larger total number $N=10000$ and later
time $T=10$, the rest of parameters remain the same as that given in Fg. 1(a).

Figure 3: Realizations of dynamics of quantum states from initial vacuum.
Numerical calculations of the temporal evolutions displayed by the contour
plots of the Husimi Q-functions accompanied with the trajectories (red curves)
of the centers of these contours. Left column: The system with a linear
Hamiltonian. Middle column: The system with the degenerate two-boson
Hamiltonian. Right column: The system with the intensity-dependent
Hamiltonian. All parameters are the same as used in Fig. 2.

In Fig. 3, the quasi-distributions given by the Husimi Q-functions are
plotted. The Husimi Q-function [2, 4, 5] is defined as
$Q(x,y)=\langle\alpha|\hat{\varrho}_{N}|\alpha\rangle/\pi$, here
$\alpha=\sqrt{x^{2}+y^{2}}$ ${\rm exp}[i\ {\rm atan}(y/x)]$. Because of
coherent state representations, the Q-functions conveniently illustrate the
coherent state as a displaced vacuum state with a perfect ring shape [8],
preserved for entire time. In Fig. 3, the red curves indicate the trajectories
of displacements of the initial coherent state over time. These trajectories
are the centers of single selected contour plots with the fixed value of the
Q-functions at any given time $t$. For example, in Fig. 3(a), this trajectory
follows a circle but is not yet complete circle opposite to that case in Fig.
3(a’). The parameters used for the plots in Fig. 3 (a,b,c,a’,b’,c’) are the
same as those used in Fig. 2(a,b,c,a’,b’,c’), respectively. In Fig.
2(b,b’,c,c’), the realizations of quantum dynamics for nonlinear Hamiltonians
with $\hat{R}(t)=\hat{b}^{2}(t)$ in (b,b’) representing two-boson processes
and $\hat{R}(t)=\hat{b}(t)\sqrt{\hat{b}^{\dagger}\hat{b}}$ in (c,c’)
representing intensity-dependent processes are demonstrated. Similar to Fig.
2, the numerical results for $\varrho_{N}(0,0)$ are compared to approximate
analytical expressions for time evolutions for the probabilities being in the
ground state $p(T)$ after time $T$ in Fig. 3(b,c). In the case of two-boson
transition processes, the approximate analytical expressions are obtained to
be

$p(T)\approx{\rm exp}\left[-\frac{8|\zeta|^{2}}{\omega^{2}}{\rm
sin}^{2}\left(\frac{\omega T}{2}\right)\right]$ (18)

For Figs. 2(b,c), the parameters are given as $N=8000$, $\Delta t=0.0001$,
$T=0.8$, $\omega=\pi$ and $\omega T/\pi=0.8$, while for Figs. 2(b’,c’), the
parameters are the same as in (b) except for frequency $\omega=4\pi$ and,
thus, $\omega T/\pi=3.2$. In Fig. 3(b), the Q-functions display how the
initial vacuum state with a ring shape is transformed to the significantly
squeezed states with its signature oval shape [19, 20, 8, 2] for a slower
process with a frequency of $\omega=\pi$. However, for the fast process with
$\omega=4\pi$, the state remains merely in vacuum state without observable
squeezing. Lastly, Figs. 2(c,c’) and 3(c,c’) represent the temporal evolutions
for the Holstein-Primakoff ${\rm SU}(1,1)$ transformed states [21, 9].
Similarly, the probabilities for slow (Fig. 2(c)) versus fast (Fig. 2(c’))
processes are compared. It is important to note that the analytical formula
for the probability for these processes is identical to Eq.(17). However, the
deviation (i.e., displacement) is not as much pronounced as for coherent
states (see, Fig. 3(c,c’)).

##  4 Conclusions

In the standard approach, quantum dynamics for arbitrary system are realized
by the time evolutions of wave functions in Hilbert space, which can also be
expressed in terms of density operators in Liouville space. However, the
standard quantum simulations may occasionally turn out to be challenging,
particularly, for nonlinear dynamical systems.

In this letter, we introduce a new nonstandard iterative technique, formulated
as follows. (i) A finite time interval is divided into a large number of
discrete subintervals with an ultrashort width. (ii) The Liouville space is
synthesized with an additional virtual space for ultrashort time duration and
the quantum system’s original Hamiltonian is modified accordingly. In
particular, the force terms are replaced with virtual quantum operators. (iii)
The density operator for the system is extracted by tracing over the virtual
operator space. In principle, various virtual operators can be chosen
depending on specific quantum system. For example, the simple algebra of using
two-state spin raising and lowering operators reduces the cost of time-
consuming calculations. After introducing our technique, we implement it to
the well-known example of a charged particle in a harmonic potential. Temporal
evolutions of the probability for the particle being in the ground state are
obtained by the present technique and compared to the analytical solutions
given by the standard approach. We further discuss the physics insight of this
technique based on a thought-experiment. Lastly, we perform numerical
simulations for temporal evolutions for the ground state probability for
generalized systems governed by the time-dependent nonlinear Hamiltonians. The
quantum dynamics are realized by using the quasi-distributions.

Successive processes implicitly ’hitchhiking’ via virtual space for discrete
ultrashort time duration, are the hallmark of our technique. We believe that
this novel technique has potential for solving numerous problems otherwise
challenging to address using the standard approach based on time-ordered
exponentials.

## References

  * [1] P. Coleman, Introduction to many body physics, (Cambridge University Press, 2015)
  * [2] M. Scully and S. Zubairy, Quantum optics, (Cambridge University Press, 1997)
  * [3] S. Mukamel, Principles of nonlinear optical spectroscopy, (Oxford University Press, 1995)
  * [4] J. Perina, Quantum Statistics of Linear and Nonlinear Optical Phenomena 2nd edn (Dordrecht: Kluwer 1991)
  * [5] L. Mandel, and E. Wolf, Optical Coherence and Quantum Optics (Cambridge: Cambridge University Press, 1995)
  * [6] Fisher R A, Nieto M M and Sandberg V D 1984 Phys. Rev. D 29 110\.
  * [7] S. L. Braunstein, and R. I. McLachlan, Generalized squeezing 1987 Phys. Rev. A 35 1659\.
  * [8] G. Ariunbold, J. Perina, and Ts. Gantsog, Nonclassical states in cavity with injected atoms 1999 J. Opt. B: Quantum Semiclass. Opt. 1 219\.
  * [9] G. Ariunbold, and J. Perina, Holsttein-Primakoff $SU(1,1)$ coherent state in micromaser under intensity dependent Janynes-Cummings interaction 1998 Acta Phys. Slov. 48 315\.
  * [10] G. Ariunbold, and J. Perina, Two-mode correlated states in cavity with injected atoms, 1999 Acta Phys. Slov. 49 627\.
  * [11] D. Bonatsos, C. Daskaloyannis, and G. A. Lalassisis, Unification of Jaynes-Cummings models 1993, Phys. Rev. A 47 3448
  * [12] P. Shanta, S. Chaturvedi, V. Srinivasan, and R. Jagannathan, Unified approach to the analogues of single-photon and multiphoton cohernet states for generalized bosonic oscillators 1994, J. Phys. A: Math. Gen. 27 6433\.
  * [13] D. Meschede, H. Walther and G. M${\rm\ddot{u}}$ller, One-atom maser 1985 Phys. Rev. Lett. 54 551\.
  * [14] F. L. Kien, M. O. Scully, and H. Walther, Generation of a coherent state of the micromaser field, 1993 Found. Phys. 23 177\.
  * [15] R. H. Dicke, Coherence in Spontaneous Radiation Processes, Physical Review, 1954, 93, 99
  * [16] M. G. Benedict, Super-radiance: Multiatomic coherent emission, (CRC Press, 1996)
  * [17] G. O. Ariunbold, A cascade superradiance model, 2022, Phys. Lett. A, 452 128468
  * [18] G. O. Ariunbold, V. A. Sautenkov, H. Li, R. K. Murawski, X. Wang, M. Zhi, T. Begzjav, A. V. Sokolov, M. O. Scully, and Yu. V. Rostovtsev, Observations of Ultrafast Superfluorescent Beatings in a Cesium Vapor Excited by Femtosecond Laser Pulses, 2022, Phys. Lett. A, 428, 127945.
  * [19] H. P. Yuen, Two-photon cohernet states of the radiation field 1976 Phys. Rev. A 13 2226
  * [20] D. Stoler, Equivalence Classes of minimum uncertainty packets, 1970 Phys. Rev. D 1 3217
  * [21] T. Holstein, and H. Primakoff, Field dependence of the intrinsic domain magnetization of a ferromagnet 1940, Phys. Rev. 58 1098\.


# Enormous variation in homogeneity and other anomalous features of room
temperature superconductor samples: a Comment on Nature 615, 244 (2023)

J. E. Hirsch  Department of Physics, University of California, San Diego, La
Jolla, CA 92093-0319

###### Abstract

The resistive transition width of a recently discovered room temperature near-
ambient-pressure hydride superconductor roomtlh  changes by more than three
orders of magnitude between different samples, with the transition temperature
nearly unchanged. For the narrowest transitions, the transition width relative
to $T_{c}$ is only $0.014\%$. The voltage-current characteristics indicate
vanishing critical current, and the normal state resistance is unusually
small. These anomalous behaviors and other issues indicate that this system is
not a superconductor. Implications for other hydrides are discussed.

Room temperature at near ambient pressure has recently been claimed for a Lu-
H-N compound roomtlh . In this paper we question the validity ot that claim on
the basis of analysis of some of the reported measurement results that show
anomalous features, and discuss the implications of this analysis.

In Extended Data (hereafter ED) Fig. 15 of Ref. roomtlh , the authors show
resistance versus temperature in the absence and presence of a magnetic field.
The relative width of the resistive transition $\Delta T/T_{c}$ shown in the
inset of the figure is $0.13$ for zero field. The authors explain the
considerable width by stating “The large transition width at zero field
indicates sample inhomogeneities, which is typical for high-pressure
experiments.”, which is not implausible.

Figure 1: Resistance versus temperature for N-doped lutetium hydride at
pressure $\sim 10kbar$, as reported in Ref. roomtlh  ED Fig. 13a and
associated reported raw data. The blue and red curves were measured under cool
down and warm up conditions respectively roomtlh .

However, in the data for resistance versus temperature in ED Fig. 13a of Ref.
roomtlh , shown here in Fig. 1, the width of the resistive transitions is only
$0.04K$, so the relative width is $\Delta T/T_{c}=0.00014$, at comparable
pressures (10 kbar vs 15 kbar respectively). For the resistance curves shown
in Fig. 2 of Ref. roomtlh , at pressures 10 kbar and 16 kbar, the relative
width is in-between those two extremes, $\Delta T/T_{c}\sim 0.008$, 60 times
larger than for ED Fig. 13a, 16 times smaller than for ED Fig. 15.

Following the logic of the authors, the sample used for ED Fig. 13a is 1,000
times more homogeneous than the sample used for ED Fig. 15, the samples used
for their Fig. 2 are in-between. The protocol used in preparing these samples
was presumably similar, as described in the “Methods” section of Ref. roomtlh
. It is not understandable why superconducting samples prepared similarly
would exhibit a degree of inhomogeneity that differs by three orders of
magnitude. The fact that the transitions shown in Fig. 1 show hysteresis is
also in conflict with what is expected for a superconductor and is not
explained in the paper.

Figure 2: Resistance versus temperature for N-doped lutetium hydride at
pressure $15kbar$. Left panel: with background subtraction, as reported in
Ref. roomtlh  ED Fig. 15. Right panel: same without background subtraction,
obtained from the raw data given in Ref. roomtlh . The black, blue and red
curves are for applied magnetic fields $0T$, $1T$, $3T$ respectively roomtlh .

It should also be noted that the authors say “In some cases, small residual
resistance from the instrument offsets was subtracted from the measured
voltage.”, but don’t specify whether “some cases” include the resistance data
shown in their Fig. 2 and ED Fig.13a. For ED Fig. 15 it is explicitly stated
that a background resistance is subtracted out roomtlh . When plotting the raw
data without background subtraction, the curves shown in Fig. 2 right panel
result bar . There is no hint of a superconducting transition in Fig. 2 right
panel.

Figure 3: Top panel: voltage versus current from source data for Fig. 2b of
Ref. roomtlh . Bottom panel: voltage versus current data for a known
superconductor, from Ref. visuper .

In addition, the width of the transitions shown in Fig. 1 is unreasonably
small. No other known superconductor exhibits such sharp transitions even at
ambient pressure, and under pressure additional broadening of the transition
results from pressure gradients. In ref. hmnature , we pointed out that the
narrow width of the transitions reported in Ref. roomt  for another room
temperature superconductor under pressure, CSH, was unreasonably small,
$\Delta T/T_{c}=0.005$, and that was 35 times $larger$ than the width seen in
Fig. 1. Other anomalies in the resistance curves of Ref. roomt  were noted in
Ref. hamlin . Ref. roomt , which has six coauthors in common with Ref. roomtlh
, was recently retracted roomtr .

Furthermore, the voltage-current characteristic shown in Fig. 2b of Ref.
roomtlh , reproduced here in Fig. 3 top panel, is not consistent with what is
expected for a superconductor at temperature well below its critical
temperature. There is no evidence of any region of zero resistance in Fig. 3
top panel, in contrast with the typical behavior shown in Fig. 3 bottom panel
for a known superconductor visuper , showing zero voltage up to a critical
current that increases as the temperature decreases. Other such measurements
for known superconductors showing how this behavior varies with temperature
and magnetic field are shown in Refs. visuper2 ; visuper3 ; visuper4 ;
visuper5 . If one nevertheless insisted to infer a non-zero critical current
from Fig. 3 top panel, it would be certainly smaller than $0.5mA$, at
temperature $T/T_{c}\sim 0.1$. According to the authors, “On average, sample
sizes are on the order of 70-100 $\mu m$ in diameter and 10-20 $\mu m$ thick”.
Assuming the smallest cross-sectional area in that range, $A\sim 10\mu m\times
70\mu m$, yields a critical current density $J_{c}<72A/cm^{2}$ (or smaller if
larger dimensions are assumed). That would be five orders of magnitude smaller
than critical current densities reported for $LaH_{10}$ and $H_{3}S$ at
temperature $T/T_{c}\sim 0.5$ e2021p , and inconsistent with the magnetization
measurements shown in Fig. 3 of Ref. roomtlh .

Figure 4: Voltage versus current for Lu-H-N. Left panel: Fig. 2b of Ref.
roomtlh . Right panel: another example reported by R. P. Dias, Ref. apstalk .

Furthermore, the values of the normal state resistance at room temperature
inferred from Fig. 2b of Ref. roomtlh , and for another example of voltage
current-characteristics presented by the lead author of Ref. roomtlh  in Ref.
apstalk , shown on the left and right panels of Fig. 4, are anomalously small:
$R\sim 1.12m\Omega$ for the left panel, $R\sim 0.05m\Omega$ for the right
panel. From the van der Pauw formula $\rho\sim\pi dR/ln(2)$ (d=thickness,
R=resistance), assuming the largest $d=20\mu m$ yields $\rho\sim 10\mu\Omega-
cm$ for the left panel of Fig. 4 and $\rho\sim 0.4\mu\Omega-cm$ for the right
panel. These values are six times and 130 times smaller that the room
temperature resistivity of lutetium metal. These discrepancies point to the
behavior reported roomtlh ; apstalk  as voltage versus current for Lu-H-N both
at room temperature and at low temperatures being due to experimental
artifacts rather than superconductivity. It is notable that for other samples
that showed a resistance drop from much larger room temperature resistance
values (Fig. 2a and ED Fig. 13 a of Ref. roomtlh ), consistent with what could
be expected, no voltage-current characteristics were shown in Refs. roomtlh ;
apstalk .

Figure 5: Three holy grails. The top three panels show resistance versus
temperature for three different compounds at pressure ranges $\sim 200GPa$
(Refs. roomt ; roomt2 ; roomt3 ), $\sim 20GPa$ (Ref. one ) and $\sim 2GPa$
(Ref. roomtlh ) respectively, all showing room temperature superconductivity.
The bottom three panels show $T_{c}$ versus pressure for the three different
compounds as reported in Refs. roomt , one , roomtlh . Note that even though
Ref. roomt  was retracted roomtr , all the authors disagreed with the
retraction.

We also point out that the ac susceptibility data shown in ED Fig. 5 of Ref.
roomtlh  before background subtraction show a background dependence on
temperature that has positive slope, negative slope, and zero slope, for the
same or comparable pressures. The background ac susceptibility is expected to
reflect the physical properties of the environment of the sample, which should
not drastically change for different measurements.

Figure 6: Resistance versus temperature for various hydrides under pressure
claimed to be high temperature superconductors. The name of the last author in
the reference is shown in each panel. From upper left to lower right: $LuHN$
roomtlh ,$CSH$ roomt , $YH_{9}$ ydias , $LaH_{10}$ hemley , $LaH_{10}$
laeremets , $YH_{6}$ yeremets , $CeH$ ceh , $H_{3}S$ e2019 , $CaH_{6}$ cah6 .

We also point out that the lead author of Ref. roomtlh  R. P. Dias and some of
its coauthors have previously reported room temperature superconductivity in
other compounds, at pressures one one  and two roomt ; roomt2 ; roomt3  orders
of magnitude larger than reported in Ref. roomtlh , as shown in Fig. 5. In the
112 years since superconductivity was discovered, no room temperature
superconductivity has been conclusively established by other researchers in
any compound at any pressure despite intensive searches. The probability that
the same research group would hit this holy grail three separate times is
insignificant.

Finally, we point out that experimental attempts to reproduce the results
reported in Ref. roomtlh  have shown no indication of superconductivity in
samples prepared by following the sample preparation method described in Ref.
roomtlh  rep1 ; rep2 ; rep3 ; rep4 ; rep5 ; rep6 , and theoretical attempts to
calculate $T_{c}$ in this system within the conventional theory of
superconductivity have found values of $T_{c}$ two orders of magnitude smaller
repth ; repth2 ; repth3 ; repth4  than reported in Ref. roomtlh .

In conclusion, the extreme sharpness of the resistive transition curves shown
in Fig. 1, the fact that the width of the resistive transition changes by
three orders of magnitude between different samples, the fact that resistance
data versus temperature plotted without background subtraction show no hint of
superconductivity, the fact that voltage-current characteristics do not show
evidence for a finite critical current, the fact that normal state resistances
measured are anomalously small, the fact that the background ac susceptibility
changes drastically in different measurements, and the fact that several
experimental and theoretical studies have not been able to reproduce the
results reported in Ref. roomtlh , indicate that the behavior observed
reported in Ref. roomtlh  is not due to superconductivity.

If the resistance drops seen in Fig. 5 and other figs. of Refs. roomtlh ;
roomt ; roomt2 ; roomt3  are not due to superconductivity they must be due to
other reasons unrelated to superconductivity, either physical phenomena or/and
experimental artifacts associated with performing resistance measurements on
very small hydrogen-rich samples under high pressure in diamond-anvil cells,
with the expectation grounded in the conventional theory of superconductivity
that superconductivity will be found pickett . The resistance versus
temperature curves shown in Fig. 5 above and in other figs. of Refs. roomtlh ;
roomt ; roomt2 ; roomt3  look similar to resistance versus temperature curves
for other hydrides under high pressure that have been claimed to be high
temperature superconductors in recent years troyan22 , as shown in Fig. 6.
This raises the possibility that those same other reasons unrelated to
superconductivity that may account for the resistance drops in Refs. roomtlh ;
roomt ; roomt2 ; roomt3  reported by Dias and coworkers could account for the
resistance drops seen in all hydrides under high pressure claimed to be high
temperature superconductors troyan22 , in contradiction with the predictions
of Ref. pickett  and in agreement with other theoretical predictions xor .
Magnetic evidence claimed to support high temperature superconductivity in
sulfur and lanthanum hydrides under pressure has been called into question
elsewhere mag1 ; hmmre ; huangmine .

###### Acknowledgements.

The author is grateful to multiple colleagues for discussions on these issues,
and in particular to Erik van Heumen for stimulating discussions on transport
measurements. Some of the anomalies reported here were also noted in comments
to Ref. roomtlh  submitted by readers that are posted at the website of Ref.
roomtlh , at pubpeer.com pubpeer , in the blog post “nanoscale views” nano ,
and in various postings in reddit.com reddit .

Competing interests: the author declares no competing interests. Data
availability statement: The data that support the findings of this study are
available from the author upon reasonable request.

## References

  * (1) Nathan Dasenbrock-Gammon, Elliot Snider, Raymond McBride, Hiranya Pasan, Dylan Durkee, Nugzari Khalvashi-Sutter, Sasanka Munasinghe, Sachith E. Dissanayake, Keith V. Lawler, Ashkan Salamat and Ranga P. Dias, “Evidence of near-ambient superconductivity in a N-doped lutetium hydride”, Nature 615, 244 (2023).
  * (2) J. E. Hirsch and F. Marsiglio, “Unusual width of the superconducting transition in a hydride”, Nature 596, E9 (2021).
  * (3) Elliot Snider, Nathan Dasenbrock-Gammon, Raymond McBride, Mathew Debessai, Hiranya Vindana, Kevin Vencatasamy, Keith V. Lawler, Ashkan Salamat and Ranga P. Dias, “Room-temperature superconductivity in a carbonaceous sulfur hydride”, Nature 586, 373 (2020).
  * (4) See also comment by Bartolomeo Osanza at website of Ref. roomtlh .
  * (5) J. J. Hamlin, arXiv:2210.10766 (2022).
  * (6) E. Snider et al, “Retraction Note: Room-temperature superconductivity in a carbonaceous sulfur hydride”, Nature 610, 804 (2022).
  * (7) T. Freltoft, H. J. Jensen and P. Minnhagen, “Evidence for intrinsic critical current density in high Tc superconductors”, Solid State Communications 78, Issue 7, 635 (1991).
  * (8) J. Lopez and P. Mune, “Anomalous inversion of voltage versus current curves in ZFC and FC conditions in BSCCO ceramic superconductors”, Physica C 261, 173 (1996).
  * (9) L. Qiao et al, “Dynamical instability of the electric transport in superconductors”, Scientific Reports 8, 14104 (2018).
  * (10) G. C. Han, J. Phys: Condens. Matter 7, 8175 (1995).
  * (11) J. W. Ekin et al, “Evidence for weak link and anisotropy limitations on the transport critical current in bulk polycrystalline $Y_{1}Ba_{2}Cu_{3}O_{x}$, Journal of Applied Physics 62, 4821 (1987).
  * (12) R. P. Dias, talk presented at APS March Meeting, 2023, “Observation of Room Temperature Superconductivity in Hydride at Near Ambient Pressure”, Session K20, K20.00002, 3:00 PM–5:48 PM, Tuesday, March 7, 2023. See youtube.com/watch?v=XhhvOMuLF94, 9:03/13:43 or https://jorge.physics.ucsd.edu/K20Dias.mov.
  * (13) V. Minkov, et al, “Magnetic field screening in hydrogen-rich high-temperature superconductors”, Nat Commun 13, 3194 (2022).
  * (14) R. P. Dias, reported at colloquium given at University of California San Diego, April 22, 2021. Same data were presented by R. P. Dias at the 10th Asian Conference on High Pressure Research (ACHPR), November 21-26, 2021, and at CDRTS 2022, L’Aquila, Italy, July 26-29, 2022.
  * (15) G. Alexander Smith, Ines E. Collings, Elliot Snider, Dean Smith, Sylvain Petitgirard, Jesse S. Smith, Melanie White, Elyse Jones, Paul Ellison, Keith V. Lawler, Ranga P. Dias and Ashkan Salamat, “Carbon content drives high temperature superconductivity in a carbonaceous sulfur hydride below 100 GPa”, Chem. Commun. 58, 9064 (2022).
  * (16) Hiranya Pasan, Elliot Snider, Sasanka Munasinghe, Sachith E. Dissanayake, Nilesh P. Salke, Muhtar Ahart, Nugzari Khalvashi-Sutter, Nathan Dasenbrock-Gammon, Raymond McBride, G. Alexander Smith, Faraz Mostafaeipour, Dean Smith, Sergio Villa Cortés, Yuming Xiao, Curtis Kenney-Benson, Changyong Park, Vitali Prakapenka, Stella Chariton, Keith V. Lawler, Maddury Somayazulu, Zhenxian Liu, Russell J. Hemley, Ashkan Salamat and Ranga P. Dias, “Observation of Conventional Near Room Temperature Superconductivity in Carbonaceous Sulfur Hydride”, arXiv:2302.08622 (2023).
  * (17) X. Ming et al, “Absence of near-ambient superconductivity in $LuH_{2\pm x}N_{y}$”, arXiv:2303.08759 (2023).
  * (18) P. Shan et al, “Pressure-induced color change in the lutetium dihydride LuH2”, arXiv:2303.06718 (2023).
  * (19) X. Xing et al, “Observation of non-superconducting phase changes in $LuH_{\pm x}N_{y}$”, arXiv:2303.17587 (2023).
  * (20) S. Cai et al, “No evidence of superconductivity in the compressed sample prepared from the lutetium foil and H2/N2 gas mixture”, arXiv:2304.03100 (2023).
  * (21) N. Wang et al, “Percolation-induced resistivity drop in cold-pressed LuH2”, arXiv:2304.00558 (2023).
  * (22) Ying-Jie Zhang et al, “Pressure induced color change and evolution of metallic behavior in nitrogen-doped lutetium hydride”, arXiv:2303.17453 (2023).
  * (23) K. P. Hilleke et al, “Structure, Stability and Superconductivity of N-doped Lutetium Hydrides at kbar Pressures”, arXiv:2303.15622 (2023).
  * (24) Y. Sun et al, “Effect of nitrogen doping and pressure on the stability of cubic $LuH_{3}$”, arXiv:2303.14034 (2023).
  * (25) Z. Huo et al, “First-principles study on the superconductivity of N-doped fcc-LuH3”, arXiv:2303.12575 (2023).
  * (26) P. P. Ferreira et al, “Search for ambient superconductivity in the Lu-N-H system”, arXiv:2304.04447 (2023).
  * (27) I. A. Troyan et al, “High-temperature superconductivity in hydrides”, and references therein.
  * (28) W. E. Pickett, “Room temperature superconductivity: The roles of theory and materials design”, Rev. Mod. Phys. 95, 021001 (2023) and references therein.
  * (29) E. Snider et al, “Synthesis of Yttrium Superhydride Superconductor with a Transition Temperature up to 262 K by Catalytic Hydrogenation at High Pressures”, Phys. Rev. Lett. 126, 117003 (2021).
  * (30) M. Somayazulu et al, “Evidence for Superconductivity above 260 K in Lanthanum Superhydride at Megabar Pressures”, Phys. Rev. Lett. 122, 027001 2019).
  * (31) A. P. Drozdov et al, “Superconductivity at 250 K in lanthanum hydride under high pressures”, Nature volume 569, 528 (2019).
  * (32)
  * (33) W. Chen et al, “High-Temperature Superconducting Phases in Cerium Superhydride with a $T_{c}$ up to 115 K below a Pressure of 1 Megabar”, Phys. Rev. Lett. 127, 117001 (2021).
  * (34) S. Mozaffari et al, “Superconducting phase diagram of $H_{3}S$ under high magnetic fields”, Nature Comm. 10, 2522 (2019) .
  * (35) L. Ma et al, “High-Temperature Superconducting Phase in Clathrate Calcium Hydride $CaH_{6}$ up to 215 K at a Pressure of 172 GPa”, Phys. Rev. Lett. 128, 167001 (2022).
  * (36) J. E. Hirsch, “Hole superconductivity xOr hot hydride superconductivity”, J. Appl. Phys. 130, 181102 (2021) and references therein.
  * (37) J. E. Hirsch and F. Marsiglio, “On magnetic field screening and expulsion in hydride superconductors”, arXiv:2207.01541 (2023).
  * (38) J. E. Hirsch and F. Marsiglio, “Clear evidence against superconductivity in hydrides under high pressure”, Matter and Radiation at Extremes 7, 058401 (2022).
  * (39) J. E. Hirsch, “Faulty evidence for superconductivity in ac magnetic susceptibility of sulfur hydride under pressure”, National Science Review, 9, nwac086 (2022).
  * (40) pubpeer (2023).
  * (41) nanoscale views (2023)
  * (42) reddit.com (2023).