# From a bistable adsorbate to a switchable interface: tetrachloropyrazine on
Pt(111)

Lukas Hörmann  Andreas Jeindl  Oliver T. Hofmann<EMAIL_ADDRESS>[

###### Abstract

Virtually all organic (opto)electronic devices rely on organic/inorganic
interfaces with specific properties. These properties are, in turn,
inextricably linked to the interface structure. Therefore, a change in
structure can introduce a shift in function. If this change is reversible, it
would allow constructing a switchable interface. We accomplish this with
tetrachloropyrazine on Pt(111), which exhibits a double-well potential with a
chemisorbed and a physisorbed minimum. These minima have significantly
different adsorption geometries allowing the formation of switchable interface
structures. Importantly, these structures facilitate different work function
changes and coherent fractions (X-ray standing wave measurements), which are
ideal properties to readout the interface state.

We perform surface structure search using a modified version of the SAMPLE
approach and account for thermodynamic conditions using ab-initio
thermodynamics. This allows investigating millions of commensurate as well as
higher-order commensurate interface structures. We identify three different
classes of structures exhibiting different work function changes and coherent
fractions. Using temperature and pressure as handles we demonstrate the
possibility of reversible switching between those different classes, creating
a dynamic interface for potential applications in organic electronics.

Graz University of Technology] Institute of Solid State Physics, Graz
University of Technology, Petersgasse 16, 8010 Graz, Austria

##  1 Introduction

Organic/inorganic interfaces are essential for the construction of organic
(opto)electronic devices. To date research has mainly focused on attaining
specific interface properties through controlling the structure and chemistry
of the organic adlayer.1, 2, 3, 4, 5 In this theoretical study we go beyond
these efforts and study organic adsobate layers with switchable properties.

Switchable interfaces in literature rely on several function principles
including modifying interactions at the interface6 and applying external
stimuli such as optical signals7, 8, electric fields9, 10, magnetic fields11,
temperature,12 biochemical processes13 or pH-value14. Here, we use temperature
and pressure to switch adlayers of molecules that exhibit a double well
potential when adsorbing on a substrate. Examples for systems with double well
potentials include benzene derivatives on Pt(111) and anthradithiophene on
Cu(111).15, 16, 17 These systems show great potential as molecular switches15
and have even been suggested as data storage.18 Here, we focus on the system
of tetracholoropyrazine (TCP) on Pt(111), which displays a double-well
potential on the Pt(111) surface.15 The two types of adsorption geometries,
called on- and off-state, exhibit different types of bonding and adsorption
geometries. Crucially, both minima have approximately the same adsorption
energy, allowing the formation of diverse surface structures which are
sufficiently close in energy to be reversibly switched with external stimuli.
Interestingly, in this work we find that individual chemisorbed molecules are
energetically more favorable than physisorbed ones. However, this ordering
reverses in continuous layers (at low temperatures).

Reading out the interface states requires that the switch in structure entails
a shift in interface properties. We focus on two of these properties, namely
the work-function change $\Delta\Phi$ (compared to the clean surface) and the
coherent fraction (as obtained by X-ray standing wave methods), and study how
to modify them using environmental stimuli such as temperature and pressure.
Different possible polymorphs exhibits diverse $\Delta\Phi$s and coherent
fractions. However, the amount by which they can be switched is significantly
smaller the total range. We identify the reasons for this and provide ideas to
overcome it.

Molecules on surfaces may arrange into a large number of possible polymorphs.
Therefore, we investigate the $\Delta\Phi$s and the coherent fractions of
millions of surface structures. To determine these structures and their
adsorption energies we use the SAMPLE approach,19 which combines dispersion-
corrected density functional theory (DFT) with machine learning (see Methods
Section). Further, ab-initio thermodynamics allows modeling the impact of
temperature and pressure on the surface polymorphism. This enables us to
conduct in-silico simulations of experiments.

A number of experimental methods exist to determine work functions (and in
turn $\Delta\Phi$s), such as photo-emission spectroscopy and the Kelvin probe
method. The change in work function resulting from molecules adsorbing on a
surface depends on a number of factors. These include pushback of electron
density into the surface, charge transfer across the interface, formation of
new interface states and bonds between adsorbates and the substrate, image
charge effects, permanent dipoles of the molecules as well as the coverage.20,
21, 22 On- and off-state differ significantly in type of bonding and
adsorption geometry. Therefore, we expect that adsorbate layers consisting of
on- and off-state adsorption geometries yield dissimilar $\Delta\Phi$s,
allowing the construction of a switchable interface.

A useful experimental method to determine the adsorption geometry is the X-ray
standing wave (XSW) technique.23 This method yields two measures, namely the
coherent position and the coherent fraction. The coherent position allows
determining the mean adsorption height of the adsorbates. The coherent
fraction is effectively an order parameter24 containing information about
differences in adsorption heights, which would occur if on- and off-state
molecules coexisted in a particular interface state. If all atoms (of a
particular species) within all molecules adsorb at an identical adsorption
height, the coherent fraction is $1$, while it decreases with variations in
adsorption height. Thus, the coherent fraction will allow differentiating
between on-, off- and mixed-state layers.

The fact that both $\Delta\Phi$ and the coherent fraction are experimentally
readily accessible makes these properties ideal candidates to readout the
state of a switchable interface.

##  2 Results

###  2.1 Individual Molecules on the Surface

When investigating the adsorption of single molecules, we are interested in
local minima of the potential energy surface (PES) which constitute
energetically favorable adsorption geometries. We will hereafter refer to
these as “geometries”. As stated above, TCP on Pt(111) exhibits a double well
potential, with minima occurring at two different adsorption heights. Hence,
two different types of adsorption geometries exist, namely on- and off-state
geometries. While on-state geometries are chemisorbed, off-state geometries
are physisorbed.

We determine the geometries in two steps: Initially, a Gaussian-process
regression (GPR) algorithm identifies the minima geometries of a coarse
grained PES (see Methods Section). The algorithm uses DFT-calculated energies
of a few adsorption geometries as input and interpolates between them. Hereby,
it only considers the most important degrees of freedom, i.e., position,
orientation and bending of the molecule (see Methods Section). If necessary
(see below), we refine the GPR minima using DFT geometry optimizations.

For off-state geometries the coarse grained PES is sufficiently accurate since
the molecules mainly bond via spatially uniform van der Waals interactions.
This leads to a weakly corrugated PES, with an energy range of approximately
$0.2~{}eV$ (see Figure 1). Hence, we can directly use the GPR minima as
adsorption geometries. In the off-state the molecules remains flat and has
five adsorption geometries with (GPR-calculated) bonding energies of
$-0.97~{}eV$ to $-0.95~{}eV$ and adsorption heights of approximately
$3.3~{}\mbox{\AA}$. DFT geometry optimizations of these minima yield only very
small changes in geometry (molecule remains flat) and the gain in adsorption
energy is within the uncertainty of our method.

Conversely, the PES of on-state geometries is strongly corrugated, with minima
and maxima spanning an energy range of approximately $2.0~{}eV$ (see Figure
1). This is due to a pronounced spatial dependence of the covalent
interactions between the molecule and the surface. The molecules are strongly
bent, with the Cl atoms pointing away from the surface. Due to the higher
complexity of the on-state, a PES considering only position, orientation and
bending of the molecule is insufficient to accurately determine adsorption
geometries. Hence, we refine the GPR minima with DFT, allowing the molecule
and the substrate atoms of the two topmost layers to relax. Hereby, the
relaxation of substrate atoms below the molecule contributes up to $-0.5~{}eV$
to the adsorption energies (more negative is more bonding). For the on-state
we find four different adsorption geometries that have adsorption energies
(calculated with DFT) of $-1.07~{}eV$ to $-1.04~{}eV$ and an adsorption height
of approximately $2.1~{}\mbox{\AA}$ above the unrelaxed substrate (This is how
an XSW-experiment would determine the adsorption height.23).

For visualizations of the potential energy surfaces and additional details
regarding adsorption geometries see the Supporting Information.

Figure 1: Cut through the GPR PES of TCP on Pt(111) in [110] direction; (red)
PES of the on-state geometry; (blue) PES of the off-state geometry; energies
are plotted relative to the respective minimum energy $E_{min}$.

###  2.2 Interface Structures

Having discussed the adsorption geometries of individual molecule on the
surface, we will now focus on the close-packed layers they form. In general,
adlayers assume different types of commensurability depending on a delicate
balance between molecule-molecule interactions and the corrugation of
molecule-substrate interactions.25, 26 Strong molecule-molecule interactions
and a comparatively weak corrugation of molecule-substrate interactions allow
maximizing the energy gain from interactions between molecules. This will most
likely lead to incommensurate layers. Conversely, a large corrugation of
molecule-substrate interactions and comparatively weak molecule-molecule
interactions forces the molecules to remain in energetically favorable
adsorption sites. Concurrently, any energy gain from favorable molecule-
molecule interactions would be outweighed by the energy penalty from
unfavorable molecule-substrate interactions. This case leads to commensurate
layers. Using the arguments from the previous paragraph, we can anticipate if
an adlayer will be commensurate or incommensurate.

The two types of local adsorption geometries can form three different classes
of adlayer structures. We will hereafter refer to a close-packed adlayer of
molecules as a “motif”.

First, off-state motifs consist purely of off-state geometries. Here, the PES
of the single molecule is weakly corrugated (see Figure 1). Hence, the
molecule can seek out the most beneficial molecule-substrate interactions,
making off-state motifs incommensurate.

Next, on-state motifs consist only of on-state geometries. The PES of on-state
geometries is strongly corrugated (see Figure 1). Therefore, the molecules
seek the most favorable molecule-substrate interaction, which leads on-state
motifs to be commensurate.

Finally, mixed-state motifs contain both on- and off-state geometries. Here,
the off-state geometries seek out the most beneficial molecule-substrate
interaction, while the on-state geometries remain stuck due to the strongly
corrugated molecule-substrate interactions. Hence, we expect mixed-state
motifs to be commensurate.

To find the energetically most favorable motifs, we first perform commensurate
structure search using the SAMPLE approach (see Methods Section).19 Hereby, we
generate all possible commensurate motifs with different coverages and up to
three molecules per unit cell and predict their adsorption energies. However,
as we state above, we expect that some energetically favorable adlayers are
incommensurate. Truly incommensurate motifs contain an infinite number of
molecules per unit cell, each with a different adsorption site. This makes it
obviously impossible to determine their energies from first principles.
Nevertheless, we can approximate incommensurate motifs using higher-order
commensurability. To consider higher-order commensurate motifs, we use a
generalized SAMPLE approach (SAMPLE-GPR) based on a GPR model (see Methods
Section). This approach also provides energies for geometry optimizations with
limited degrees of freedom.

These additional capabilities are necessary to describe higher-order
commensurate off-state motifs and to optimize molecule positions in mixed-
state motifs, while the SAMPLE approach would be sufficient for on-state
motifs. For a consistent treatment of all three motif classes, we rerank the
$1000$ energetically most favorable motifs of every type and considered
coverage (in total approximately $37000$ motifs) via SAMPLE-GPR and consider
only SAMPLE-GPR energies hereafter.

For the off-state we expect that incommensurate motifs are most energetically
favorable. Hence, we use SAMPLE-GPR to optimize the $20$ most energetically
favorable motifs using the unit cell parameters as well as location and
orientation of the molecules as degrees of freedom (see Methods Section).

In mixed-state structures the physisorbed molecules can move relatively
freely. Therefore, we optimize their positions and orientations in the $40$
energetically most favorable mixed-state motifs of every coverage, while
keeping the chemisorbed molecules fixed.

Figure 2 shows these reranked (and optimized) energies per area plotted
against the coverage for on-, off- and mixed-state motifs. In this chapter, we
consider a temperature of $0~{}K$ where close-packed adsorbate layers seek to
minimize the energy per area, making this the measure of interest.27 At
temperature above $0~{}K$, which we discuss in the next chapter, the Gibbs
free energy of adsorption, which contains the energy per area as well as
temperature-dependent potentials, becomes the relevant measure.28

Figure 2: Adsorption energies per area for the $1000$ energetically most
favorable, as well as the optimized on-, off- and mixed-state motifs at all
considered coverages; the lines show the convex hull of the minimum energy;
the shaded areas around the convex hulls show the prediction uncertainty.

The different classes of motifs have their global minimum energy motif at
different coverages. For the off-state the energetically most favorable motifs
are higher-order commensurate, strongly indicating a preference of
incommensurate layers. This allows off-state motifs to pack more densely than
on- and mixed-state motifs. The energetically most favorable off-state motif
(and overall energetically most favorable motif) has a coverage of
$1.76~{}molecules/nm^{2}$. Mixed-state motifs also profit from closer packing.
Although they are commensurate, the off-state geometries can adjust their
positions relatively freely, allowing for more beneficial molecule-molecule
interactions and tighter packing than in on-state motifs. Mixed-state motifs
have a minimum at a coverage of $1.67~{}molecules/nm^{2}$. On-state motifs
have the lowest packing density with the energetic minimum occurring at a
coverage of $1.58~{}molecules/nm^{2}$. This is due to the strongly corrugated
PES of on-state geometries which prevents them from leaving their potential
well and thus makes it difficult to achieve favorable molecule-molecule
interactions (see Figure 1).

Isolated on-state geometries have an approximately $0.10~{}eV$ more attractive
molecule-substrate interaction than off-state geometries. However, in a
tightly packed layer the molecules exhibit attractive molecule-molecule
interactions. These interactions outweigh the energy penalty accrued from off-
state geometries. This makes densely packed off- and mixed-state motifs
energetically more favorable than on-state motifs.

###  2.3 Switchable Interface States

Having discussed possible interface motifs, we will now investigate how to use
them for a switchable interface. This requires controllably shifting the work-
function change ($\Delta\Phi$) or the coherent fraction by a measurable, or
better, a technically applicable margin. Hereby, the coherent fraction is
technologically less relevant, but allows validating the switching behavior,
since it is not directly related to $\Delta\Phi$. Conversely, $\Delta\Phi$ is
highly relevant for devices29, 30, 31 and will therefore serve as focus of
this discussion. For instance, a change in $\Delta\Phi$ of $200~{}meV$ would
correspond to the typical threshold voltage of a germanium diode. Therefore,
the motifs should provide an equally large or larger variety of $\Delta\Phi$s.
To test this, we will present an overview of the possible $\Delta\Phi$s and
coherent fractions.

In off-state motifs the molecules are flat and assume very similar adsorption
heights. The geometric uniformity of different possible polymorphs leads to
similar $\Delta\Phi$s with differences being largely due to coverage. Motifs
with similar coverage exhibit $\Delta\Phi$s that vary by only $100~{}meV$ and
the largest $\Delta\Phi$ is $-690~{}meV$. Due to the similar adsorption
geometries all off-state motifs exhibit coherent fractions of C-atoms (as well
as N- and Cl-atoms) of close to $1.00$.

In on-state structures the molecules are strongly bent. Moreover, the
different on-state geometries exhibit different adsorption heights, bending
and tilting. These dissimilarities result in motifs with similar coverages
exhibiting $\Delta\Phi$s that differ by as much as $600~{}meV$ from each
other. The maximum $\Delta\Phi$ is $-1020~{}meV$. Furthermore, the distortions
of the adsorption geometries result in different z-positions of the C-atoms,
for which we therefore find coherent fractions as low as $0.76$. We note in
passing that the impact of the internal molecular geometry on the coherent
fractions is known and has been used to elucidate adsorption geometries.24, 32

Mixed-state motifs contain flat lying off-state geometries as well as strongly
distorted on-state geometries. Due to different adsorption geometries, mixed-
state motifs of similar coverage exhibit $\Delta\Phi$s that differ by as much
as $500~{}meV$ from each other. The largest $\Delta\Phi$ is $-930~{}meV$.
Since molecules sit at different adsorption heights, the z-position of C-atoms
in off-state geometries is approximately $1.2~{}\mbox{\AA}$ higher than that
in on-state geometries, which is about half of the Pt lattice distance
($2.35~{}\mbox{\AA}$). Therefore, individual motifs exhibit coherent fractions
for C-atoms that are as low as $0.01$. Coherent fractions smaller than $1$ are
sometimes taken to indicate a disordered structure.33, 34 However, our motifs
are highly ordered and commensurate. While thought experiments have shown the
possibility of such a behavior,24 this is, to our knowledge, the first report
of such low coherent fractions for a system of lying molecules which could, in
principle, be observed in an experiment. This effect is seen most strongly for
C- and N-atoms and is less pronounced for Cl-atoms.

For additional plots and details regarding $\Delta\Phi$ and the coherent
fraction see the Supporting Information.

This discussion shows that on-, off- and mixed-state motifs theoretically have
a large enough diversity in $\Delta\Phi$s and coherent fractions to construct
a switchable interface. Now we must find a way to control this diversity. An
obvious strategy would be using temperature and pressure to shift the
thermodynamic equilibrium and thereby influence which motifs form. In
thermodynamic equilibrium the thermal occupation governs the probability of
finding a particular motif. In fact, multiple motifs may coexist on the
surface and contribute to an average $\Delta\Phi$ or coherent fraction.
Therefore, we must consider a set of motifs, rather than just the
energetically most favorable one. We use a sufficiently large set containing
the $37000$ energetically most favorable motifs as well as the clean Pt(111)
surface (details in the Supporting Information). To include the influence of
temperature and pressure in the adsorption energies of these motifs we use ab-
initio thermodynamics.28 This allows generating the phase diagram which we
show in Figure 3. Panels (a) and (b) show how the expectation values of
$\Delta\Phi$ and the coherent fraction depend on the temperature and the
molecule’s partial pressure. Panels (c) and (d) show the temperature
dependence (at a constant pressure of $10^{-6}~{}Pa$) of both properties
separated into contributions of on-, off- and mixed-state motifs. Regarding
Figure 3c, we note that initially one would expect that the $\Delta\Phi$ of
mixed-state motifs lies between that of on- and off-state motifs. However,
mixed-state motifs have the lowest expectation values for $\Delta\Phi$ below
$200~{}K$ and the highest one above. This results from the fact that here the
adsorption energies do not correlate with $\Delta\Phi$ and on-state motifs
with large $\Delta\Phi$ are not energetically favorable (see Supporting
Information).

Figure 3: (a-b) Thermodynamically populated phase diagram showing (a) the
expectation value of $\Delta\Phi$ and (b) the expectation value of the
coherent fraction of C atoms, hatched areas indicate the thermodynamic range
where adsorption is not energetically favorable; (c-d) cut through the phase
diagram at a constant pressure of $10^{-6}~{}Pa$ showing expectation values
for (c) $\Delta\Phi$ and (d) the coherent fraction separately for on-, off-
and mixed-state motifs; the line width shows how much each type of motif
contributes to the thermal occupation of all states.

Within the phase diagram we define a number of experimentally accessible
temperatures and pressures at which we readout the interface state. Pressure
in our case refers to the partial pressure of the TCP molecules in gas phase.
Since this is difficult to control in practice, we will primarily switch our
interface using temperature. For the discussion we consider temperatures of
common coolants, namely that of liquid helium ($\approx 4~{}K$), that of
liquid nitrogen ($\approx 77~{}K$) and room temperature ($\approx 300~{}K$).

At $4~{}K$, the thermal occupation is dominated by higher-order commensurate
off-state motifs (off-state interface). Here we predict a $\Delta\Phi$ of
approximately $-600~{}meV$, which matches the contribution from off-state
motifs in Figure 3c. Since molecules in off-state motifs adsorb at similar
adsorption heights and are undistorted, the expectation value for the coherent
fraction is $1.00$ (see Figure 3d).

At $77~{}K$, mixed-state motifs dominate the thermal occupation (mixed-state
interface). Here we find a $\Delta\Phi$ of approximately $-650~{}meV$ in
accordance with the expectation value for mixed-state motifs in Figure 3c. The
fact that these motifs contain both on- and off-state geometries (which adsorb
at different heights) leads to a decrease of the coherent fraction to about
$0.67$.

At $300~{}K$, a large number of motifs, including off- and mixed-state motifs,
contribute to the thermal occupation (see Supporting Information). However,
on-state motifs constitute the majority (on-state interface). The $\Delta\Phi$
lies between $-550~{}meV$ and $-500~{}meV$. The coherent fraction at $300~{}K$
is approximately $0.70$. Furthermore, at partial pressures lower than
approximately $10^{-7}~{}Pa$ it is no longer energetically favorable for
molecules to adsorb on the substrate at all (indicated by the hatched areas in
Figure 3).

In passing, we must discuss the uncertainty of our prediction method and the
impact this has on the phase diagram. The uncertainty of our adsorption energy
predictions is approximately $0.04~{}eV$ per molecule. The small energy
differences between off- and mixed-state motifs lie within this uncertainty.
This mostly affects the off-state interface, where only very few motifs
contribute to the thermal occupation. At higher temperatures a larger number
of motifs contribute to the thermal occupation averaging out the prediction
error. A detailed investigation of the (minor) influence of those errors on
our phase diagrams is provided in the Supporting Information.

Based on our results, two ways of switching the interface exist: A switch from
the off-state to the mixed-state interface shifts the coherent fraction by
approximately $0.3$. Switching from the mixed-state to the on-state interface
shifts $\Delta\Phi$ by more than $100~{}meV$. This is in the order of
magnitude of the typical threshold voltage of a germanium diode ($200~{}mV$),
demonstrating the possibility of using TCP on Pt(111) to construct a
switchable interface.

However, the amount by which we can switch $\Delta\Phi$ is an order of
magnitude smaller than the range of possible $\Delta\Phi$s. This is partly due
to the thermal occupation of energetically higher-lying motifs leading to an
averaged interface property. However, the main reason comes to light when
looking closely at the adsorption geometries and their surface dipoles. On-
and off-state geometries exhibit significantly different adsorption geometries
leading to different surface dipoles (see Figure 4).

Figure 4: Surface dipoles and (Mulliken) change transfer for on- and off-state
geometries.

Off-state geometries are very similar and have a surface dipole of $-1.2~{}D$
(data shown for most stable off-state geometry). Since they are physisorbed,
interface charge transfer is small with molecules receiving a Mulliken charge
of $0.1$ electrons.

Conversely, on-state geometries can be separated into two significantly
different groups: Three on-state geometries bind to the surface via N-atoms
(N-bonded) while the fourth binds via C-atoms. The three N-bonded geometries
exhibit a surface dipole of approximately $-1.7~{}D$ and receive a Mulliken
charge of $0.4$ electrons (data shown for most stable N-bonded geometry). The
C-bonded geometry has a surface dipole of only $-0.6~{}D$ while getting a
Mulliken charge of $0.2$ electrons. However, all on-state geometries have
similar adsorption energies.

Therefore, on- and mixed-state motifs containing the C-bonded geometry yield a
small $\Delta\Phi$ (comparable to the off-state). Since such motifs are
energetically favorable, they contribute to the thermal occupation and
decrease the magnitude of the switch in $\Delta\Phi$. A possible solution
would therefore be using a bi-stable system where on- and off-state geometry
have similar surface dipoles within the respective state but differ when
compared to each other. This could be fulfilled by a molecule that exhibits
only one type of bonding chemistry which would lead to more similar on-state
geometries.

##  3 Conclusion

In this contribution we discuss a switchable interface using the different
adsorption states of TCP on Pt(111). We perform surface structure search and
find three different classes of motifs that exhibit different $\Delta\Phi$s
and coherent fractions. Using temperature and pressure as handles we can
switch between three interface states. A change from the off-state to the
mixed-state interface shifts the coherent fraction by approximately $0.3$.
More interestingly, switching from the mixed-state to the on-state interface
shifts $\Delta\Phi$ by more than $100~{}meV$. This is in the order of
magnitude of the typical threshold voltage of a germanium diode ($200~{}mV$),
demonstrating the possibility of using TCP on Pt(111) to construct a
switchable interface.

However, the achieved switch is small compared to the range of possible
$\Delta\Phi$s. Aside form the thermal occupation averaging interface
properties, this has two reasons: (i) on-state geometries exhibit different
bonding chemistry and surface dipoles and (ii) motifs with extreme
$\Delta\Phi$s are energetically unfavorable. This could be overcome by using a
molecule which allows for only one type of bonding chemistry.

##  4 Methods

To perform first principles calculations, we use the FHI-aims code35 with
numerical atom centered orbitals, the PBE exchange correlation functional36
and the TSsurf vdW correction37, 38. We use periodic boundary conditions to
model continuous layers. Since we are dealing with surfaces, we employ the
repeated slab approach and decouple the unit cells vertically by using a
vacuum of approximately $90~{}\mbox{\AA}$ as well as a dipole correction.39
Furthermore, we use k-grids equivalent to the $(48,48,1)$ grid of the
primitive substrate unit-cell.

The property we primarily consider is the adsorption energy (or bonding
energy) $E_{ads}$. We define $E_{ads}$, using the total energy of the combined
system $E_{mol+sub}$, the energy of a tetracholropyrazine molecule in vacuum
$E_{mol}$ and the energy of the clean Pt(111) substrate $E_{sub}$.

$E_{ads}=E_{mol+sub}-E_{mol}-E_{sub}$ (1)

To perform structure search for commensurate interfaces we use the SAMPLE
approach.19 SAMPLE uses coarse-graining and Bayesian linear regression. The
key premise of this approach is that the unit cell of the adsorbate layer is a
super cell of the substrate. This allows us to generate all possible motifs
within a given range of coverages and a limited number of molecules per unit
cell. The first step of building motifs is generating an exhaustive list of
substrate super cells. Then we place molecules into these unit cells. Hereby
we coarse-grain the possible positions and orientations of the molecule on the
surface. Therefore, we determine the local minima geometries of the isolated
molecule on the substrate. This yields a small number of molecular geometries
that we can place in the substrate super cells to assemble motifs. The number
of possible motifs we regularly deal with is in the order of $10^{6}$. Since
first principles calculations cannot be done on all these motifs we use an
energy model (equation (2)) to determine the adsorption energies $E_{ads}$.
The energy model consists of molecule-substrate and molecule-molecule
interactions ($U_{i}$ and $V_{p}$). The sum over all interactions that occur
in a motif then yields its energy $E_{ads}$.

$E_{ads}=\sum_{i}n_{i}U_{i}+\sum_{p}n_{p}V_{p}$ (2)

The interactions $U_{i}$ and $V_{p}$ are initially unknown. To determine them,
we use a type of machine learning called Bayesian linear regression. Hereby,
we compute a small, D-optimally selected, number of motifs with DFT. Using
these calculations, the Bayesian linear regression algorithm then learns the
molecule-substrate and molecule-molecule interactions, which allows predicting
the energies of all motifs.

Since we expect that tetrachloropyrazine on Pt(111) forms incommensurate
layers we need the ability to also predict such motifs. Therefore, we
generalize the SAMPLE approach (SAMPLE-GPR). We replace the energy model with
a Gaussian process, which drops the requirement of discrete molecule-substrate
and molecule-molecule interactions. This algorithm is similar to that
described in a previous publication.40 GPR based algorithms have been used
before to find adsorption geometries of individual molecules on substrates.41
Our algorithm can handle isolated molecule on the surface as well as
continuous layers. Put simply, a GPR algorithm is a sophisticated method to
interpolate adsorption energies and work functions (or other scalar
properties). Hereby the key assumption is that two geometries/motifs that are
geometrically similar have similar properties. The trick is finding a good
measure for this similarity. We use radial distance functions (RDF)
$\boldsymbol{f}$ which contain interatomic distances between molecules and the
substrate as well as interatomic distances between molecules. To determine the
similarity $C_{\alpha\beta}$ of two geometries/motifs, $\alpha$ and $\beta$,
we only need to calculate the overlap integral between the two RDFs
$\boldsymbol{f_{\alpha}}$ and $\boldsymbol{f_{\beta}}$.

$C_{\alpha\beta}=\langle\boldsymbol{f_{\alpha}},\boldsymbol{f_{\alpha}}\rangle=\int\boldsymbol{f_{\alpha}}(x)\cdot\boldsymbol{f_{\beta}^{*}}(x)~{}\mathrm{d}x$
(3)

The RDFs are normed such that
$\langle\boldsymbol{f_{\alpha}},\boldsymbol{f_{\alpha}}\rangle=1$. Therefore,
$C_{\alpha\beta}$ is $1$ for identical geometries/motifs and yields smaller
values for dissimilar geometries/motifs.

Like SAMPLE, SAMPLE-GPR also requires training data. For isolated molecules we
evaluate the model uncertainty and choose the data points with the largest
uncertainty. In case of continuous adlayers we reuse the training sets from
the SAMPLE approach.

To optimize motifs we employ simulated annealing. Hereby SAMPLE-GPR provides
the energy predictions. Our algorithm can optimize the most relevant degrees
of freedom, namely the unit cell parameters, as well as the positions and
orientations of all molecules in the unit cell.

So far we have only discussed surface structure search relying on energies
from DFT calculations. These energies do not account for the effects of
temperature and pressure, which are vital for a comparison with experiment. To
model the impact of different thermodynamic conditions, we use ab-initio
thermodynamics.28 Hereby we consider the thermodynamic equilibrium at a given
temperature and pressure making the Gibbs free energy of adsorption the
measure of interest. When determining the Gibbs free energy we neglect the
contributions of the vibration enthalpy, the configuration entropy and the
mechanical work as is commonly done in literature.28, 42, 43 Using the Gibbs
free energy, we can compute the probability for each motif to occur at a given
temperature and pressure. All probabilities combined yield the thermal
occupation, which we can use to determine expectation values for $\Delta\Phi$s
and coherent fractions. Therefore, we calculate the mean weighted by the
thermal occupation.

For additional details regarding the methods we refer to the Supporting
Information.

Financial support by the FWF START award (Y 1157-N36) is gratefully
acknowledged. Computational results have been achieved using the Vienna
Scientific Cluster (VSC).

The Supporting Information contains convergence tests, additional explanations
of the methods used as well as additional details regarding intermediate and
final results including phase diagrams and uncertainty estimations.

##  5 Data Availability

The data that support the findings of this study are openly available in the
NOMAD repository at https://nomad-lab.eu/prod/rae/gui/
dataset/id/WQfZt1jPQqC8gSIOU-tocA.

## References

  * Casalini et al. 2017 Casalini, S.; Bortolotti, C. A.; Leonardi, F.; Biscarini, F. Self-Assembled Monolayers in Organic Electronics. _Chemical Society Reviews_ 2017, _46_ , 40–71
  * Slowinski et al. 1997 Slowinski, K.; Chamberlain, R. V.; Miller, C. J.; Majda, M. Through-Bond and Chain-to-Chain Coupling. Two Pathways in Electron Tunneling Through Liquid Alkanethiol Monolayers on Mercury Electrodes. _Journal of the American Chemical Society_ 1997, _119_ , 11910–11919
  * Rivnay et al. 2018 Rivnay, J.; Inal, S.; Salleo, A.; Owens, R. M.; Berggren, M.; Malliaras, G. G. Organic Electrochemical Transistors. _Nature Reviews Materials_ 2018, _3_ , 1–14
  * Katz and Huang 2009 Katz, H. E.; Huang, J. Thin-Film Organic Electronic Devices. _Annual Review of Materials Research_ 2009, _39_ , 71–92
  * Wu et al. 2018 Wu, X.; Mao, S.; Chen, J.; Huang, J. Strategies for Improving the Performance of Sensors Based on Organic Field-Effect Transistors. _Advanced Materials_ 2018, _30_ , 1705642
  * Li et al. 2019 Li, J.; Yang, S.; Ren, J.-C.; Su, G.; Li, S.; Butch, C. J.; Ding, Z.; Liu, W. Deep Molecular Orbital Driven High-Temperature Hydrogen Tautomerization Switching. _The Journal of Physical Chemistry Letters_ 2019, _10_ , 6755–6761
  * Blinov et al. 2002 Blinov, L. M.; Kawai, T.; Kozlovsky, M. V.; Kawata, Y.; Ichimura, K.; Seki, T.; Tripathy, S.; Li, L.; Oliveira Jr, O.; Irie, M., et al. _Photoreactive Organic Thin Films_ ; Elsevier, 2002
  * Szymański et al. 2013 Szymański, W.; Wu, B.; Poloni, C.; Janssen, D. B.; Feringa, B. L. Azobenzene Photoswitches for Staudinger–Bertozzi Ligation. _Angewandte Chemie International Edition_ 2013, _52_ , 2068–2072
  * Fujii et al. 2019 Fujii, S.; Koike, M.; Nishino, T.; Shoji, Y.; Suzuki, T.; Fukushima, T.; Kiguchi, M. Electric-Field-Controllable Conductance Switching of an Overcrowded Ethylene Self-Assembled Monolayer. _Journal of the American Chemical Society_ 2019, _141_ , 18544–18550
  * Tam et al. 2010 Tam, T. K.; Pita, M.; Trotsenko, O.; Motornov, M.; Tokarev, I.; Halámek, J.; Minko, S.; Katz, E. Reversible “Closing” of an Electrode Interface Functionalized with a Polymer Brush by an Electrochemical Signal. _Langmuir_ 2010, _26_ , 4506–4513
  * Willner and Katz 2003 Willner, I.; Katz, E. Magnetic Control of Electrocatalytic and Bioelectrocatalytic Processes. _Angewandte Chemie International Edition_ 2003, _42_ , 4576–4588
  * Song and Hu 2010 Song, S.; Hu, N. “On-Off” Switchable Bioelectrocatalysis Synergistically Controlled by Temperature and Sodium Sulfate Concentration Based on Poly(N-isopropylacrylamide) Films. _The Journal of Physical Chemistry B_ 2010, _114_ , 5940–5945
  * Tam et al. 2008 Tam, T. K.; Zhou, J.; Pita, M.; Ornatska, M.; Minko, S.; Katz, E. Biochemically Controlled Bioelectrocatalytic Interface. _Journal of the American Chemical Society_ 2008, _130_ , 10888–10889
  * Bocharova and Katz 2012 Bocharova, V.; Katz, E. Switchable Electrode Interfaces Controlled by Physical, Chemical and Biological Signals. _The Chemical Record_ 2012, _12_ , 114–130
  * Liu et al. 2013 Liu, W.; Filimonov, S. N.; Carrasco, J.; Tkatchenko, A. Molecular Switches from Benzene Derivatives Adsorbed on Metal Surfaces. _Nature Communications_ 2013, _4_ , 1–6
  * Borca et al. 2015 Borca, B.; Schendel, V.; Pétuya, R.; Pentegov, I.; Michnowicz, T.; Kraft, U.; Klauk, H.; Arnau, A.; Wahl, P.; Schlickum, U., et al. Bipolar Conductance Switching of Single Anthradithiophene Molecules. _ACS Nano_ 2015, _9_ , 12506–12512
  * Schendel et al. 2016 Schendel, V.; Borca, B.; Pentegov, I.; Michnowicz, T.; Kraft, U.; Klauk, H.; Wahl, P.; Schlickum, U.; Kern, K. Remotely Controlled Isomer Selective Molecular Switching. _Nano Letters_ 2016, _16_ , 93–97
  * Chen et al. 2020 Chen, C.; Yang, S.; Su, G.; Li, J.; Ren, J.-C.; Liu, W. Conductance Switching in Molecular Self-Assembled Monolayers for Application of Data Storage. _The Journal of Physical Chemistry C_ 2020, _125_ , 1069–1074
  * Hörmann et al. 2019 Hörmann, L.; Jeindl, A.; Egger, A. T.; Scherbela, M.; Hofmann, O. T. SAMPLE: Surface Structure Search Enabled by Coarse Graining and Statistical Learning. _Computer Physics Communications_ 2019, _244_ , 143–155
  * Zojer et al. 2019 Zojer, E.; Taucher, T. C.; Hofmann, O. T. The Impact of Dipolar Layers on the Electronic Properties of Organic/Inorganic Hybrid Interfaces. _Advanced Materials Interfaces_ 2019, _6_ , 1900581
  * Otero et al. 2017 Otero, R.; VÃ¡zquez de Parga, A.; Gallego, J. Electronic, Structural and Chemical Effects of Charge-Transfer at Organic/Inorganic Interfaces. _Surface Science Reports_ 2017, _72_ , 105–145
  * Braun et al. 2009 Braun, S.; Salaneck, W. R.; Fahlman, M. Energy-Level Alignment at Organic/Metal and Organic/Organic Interfaces. _Advanced Materials_ 2009, _21_ , 1450–1472
  * Woodruff et al. 1987 Woodruff, D.; Seymour, D.; McConville, C. F.; Riley, C.; Crapper, M.; Prince, N.; Jones, R. G. Simple X-Ray Standing-Wave Technique and Its Application to the Investigation of the Cu (111)($\sqrt{3}\sqrt{3}$) R30°-Cl Structure. _Physical Review Letters_ 1987, _58_ , 1460
  * Woodruff and Duncan 2020 Woodruff, D.; Duncan, D. X-Ray Standing Wave Studies of Molecular Adsorption: Why Coherent Fractions Matter. _New Journal of Physics_ 2020, _22_ , 113012
  * Hooks et al. 2001 Hooks, D. E.; Fritz, T.; Ward, M. D. Epitaxy and Molecular Organization on Solid Substrates. _Advanced Materials_ 2001, _13_ , 227–241
  * Mannsfeld and Fritz 2005 Mannsfeld, S.; Fritz, T. Understanding Organic–Inorganic Heteroepitaxial Growth of Molecules on Crystalline Substrates: Experiment and Theory. _Physical Review B_ 2005, _71_ , 235405
  * Reuter and Scheffler 2003 Reuter, K.; Scheffler, M. First-Principles Atomistic Thermodynamics for Oxidation Catalysis: Surface Phase Diagrams and Catalytically Interesting Regions. _Physical Review Letters_ 2003, _90_ , 046103
  * Rogal and Reuter 2007 Rogal, J.; Reuter, K. _Experiment, Modeling and Simulation of Gas-Surface Interactions for Reactive Flows in Hypersonic Flights_ ; RTO/NATO, 2007; pp 2–1 – 2–18
  * Fahlman et al. 2019 Fahlman, M.; Fabiano, S.; Gueskine, V.; Simon, D.; Berggren, M.; Crispin, X. Interfaces in Organic Electronics. _Nature Reviews Materials_ 2019, _4_ , 627–650
  * Li and Lu 2021 Li, P.; Lu, Z.-H. Interface Engineering in Organic Electronics: Energy-Level Alignment and Charge Transport. _Small Science_ 2021, _1_ , 2000015
  * Chen et al. 2020 Chen, Q.; Wang, C.; Li, Y.; Chen, L. Interfacial Dipole in Organic and Perovskite Solar Cells. _Journal of the American Chemical Society_ 2020, _142_ , 18281–18292
  * Mercurio et al. 2013 Mercurio, G.; Maurer, R. J.; Liu, W.; Hagen, S.; Leyssner, F.; Tegeder, P.; Meyer, J.; Tkatchenko, A.; Soubatch, S.; Reuter, K., et al. Quantification of Finite-Temperature Effects on Adsorption Geometries of $\pi$-Conjugated Molecules: Azobenzene/Ag (111). _Physical Review B_ 2013, _88_ , 035421
  * Okasinski et al. 2004 Okasinski, J.; Kim, C.-Y.; Walko, D.; Bedzyk, M. J. X-Ray Standing Wave Imaging of the 1/3 Monolayer Sn/Ge-(111) Surface. _Physical Review B_ 2004, _69_ , 041401
  * Bedzyk and Liu 1995 Bedzyk, M.; Liu, C. Order-to-Disorder Phase-Transition Study of Pb on Ge (111). _Physical Review B_ 1995, _51_
  * Blum et al. 2009 Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab Initio Molecular Simulations with Numeric Atom-Centered Orbitals. _Computer Physics Communications_ 2009, _180_ , 2175 – 2196
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Physical Review Letters_ 1996, _77_ , 3865–3868
  * Grimme 2006 Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. _Journal of Computational Chemistry_ 2006, _27_ , 1787–1799
  * Ruiz et al. 2012 Ruiz, V. G.; Liu, W.; Zojer, E.; Scheffler, M.; Tkatchenko, A. Density-Functional Theory with Screened van der Waals Interactions for the Modeling of Hybrid Inorganic-Organic Systems. _Physical Review Letters_ 2012, _108_ , 146103
  * Neugebauer and Scheffler 1992 Neugebauer, J.; Scheffler, M. Adsorbate-Substrate and Adsorbate-Adsorbate Interactions of Na and K Adlayers on Al (111). _Physical Review B_ 1992, _46_ , 16067
  * Hörmann et al. 2020 Hörmann, L.; Jeindl, A.; Hofmann, O. T. Reproducibility of Potential Energy Surfaces of Organic/Metal Interfaces on the Example of PTCDA on Ag (111). _The Journal of Chemical Physics_ 2020, _153_ , 104701
  * Todorović et al. 2019 Todorović, M.; Gutmann, M. U.; Corander, J.; Rinke, P. Bayesian Inference of Atomistic Structure in Functional Materials. _NPI Computational Materials_ 2019, _5_ , 1–7
  * Herrmann and Heimel 2015 Herrmann, P.; Heimel, G. Structure and Stoichiometry Prediction of Surfaces Reacting with Multicomponent Gases. _Advanced Materials_ 2015, _27_ , 255–260
  * Reuter and Scheffler 2001 Reuter, K.; Scheffler, M. Composition, Structure, and Stability of RuO2 (110) as a Function of Oxygen Pressure. _Physical Review B_ 2001, _65_ , 035406


# Spectral Monic Chebyshev Approximation for Higher Order Differential
Equations

M. Abdelhakem111Correspondence, Email<EMAIL_ADDRESS>Department of
Mathematics, Faculty of Science, Helwan University, Cairo, Egypt Canadian
International College, School of Engineering, Cairo, Egypt  Aya Ahmed
Department of Mathematics, Faculty of Science, Helwan University, Cairo, Egypt
M. El-kady  Department of Mathematics, Faculty of Science, Helwan University,
Cairo, Egypt Canadian International College, School of Engineering, Cairo,
Egypt

###### Abstract

This paper is focused on performing a new method for solving linear and
nonlinear higher-order boundary value problems (HBVPs). This direct numerical
method based on spectral method. The trial function of this method is the
Monic Chebyshev polynomials (MCPs). This method was relying on derivative of
MCPs which explicit in the series expansion. The advantage of this method is
solved HBVPs without transforming it to a system of lower-order ordinary
differential equations (ODEs). This method supported by examples of HBVPs in
wide application. The mentioned examples showed that the proposed method is
efficient and accurate.

Keywords: High-order Boundary value problems, Monic Chebyshev polynomial,
Spectral method, Galerkin method

##  1 Introduction

In the last few decades, the popularity of the spectral methods in approximate
computational has been gained. It has proven highly effective for solving many
problems in different branches like physics and chemistry [1, 2]. Spectral
methods solved ODEs by expressing these equations in terms of a series of
known and smooth functions. The basic concept of spectral methods is to use a
set of trial functions, also called basis or expansion approximating
functions. Very smooth, globally and orthogonal are considered to be the
essential properties of this polynomials.

Many applications solved by spectral methods which provided better results [3,
4, 5] time-space with subdiffusion and superdiffusion, Abel’s integral
equations and the multi-dimensional fractional rayleigh-stokes problem in
fluid, respectively, have been solved by spectral methods. Worth mentioning,
spectral method has been displayed in [6], fifth kind Chebyshev [7] and is
still being used and developed in [8].

The superior approximation of the spectral methods depends on the type of the
basis function. It involved three types namely the galerkin, tau and
collocation (pseudospectral) methods. Authors in [9] mentioned galerkin method
by monic gegenbauer polynomial. K. Issa and R. B. Adeniyi [10] presented
solution of non-Linear ordinary differential equations depends on tau method.
A. Napoli and W. M. Abd-Elhameed in [11] solved initial value problems based
on collocation method.

The HBVPs were treated by an enormous number of algorithms using numerical
analysis. Many applications are modeled by HBVPs. It appears in a variety of
different fields of engineering, technology, applied science, physics and
fluid dynamics are described by boundary value problems (BVPs). Authors in
[12] presented application of HBVPs.

Many authors have treated with high odd-order and high even-order BVPs,
particularly [13]. The solutions of second-order BVPs mentioned in [14] by
predictor-corrector method. The fourth-order BVPs solved numerically in many
papers by a different methods [15], [16] and [17]. The sixth-order is
mentioned in [18] by Legendre wavelet collocation method, also solved in [19].
On the other hand, many authors solved high odd-order BVPs individuals [20].
The fifth-order BVPs were solved by [21].

It should be noted that there are many numbers of the literature for solving
HBVPs. In [22], F. A. Costabile and A. Napoli solved HBVPs by birkhoff
Lagrange collocation methods. Also, the ultraspherical method has been
mentioned in [23].

The trial function used in this paper is The MCPs. It is defined by dividing
Chebyshev polynomial by its leading coefficient. In [24], authors constructed
MCPs differentiation matrix. It used in many applications as optimal control
[25].

The outline of this paper is: in section 2, we expose an overview of MCPs and
the applicable characteristics required in what follows. The description of
the MCPs method (the m-th derivatives of MCPs) is proven, in section 3.
Section 4, the algorithm of the method is given. Section 5 exhibited some
numerical results which revealed the accuracy and efficiency of our postulated
numerical algorithms. In The 6th and last section, few concluding remarks were
presented.

##  2 Preliminaries

In this section, we shall present the base formula for MCPs and the relation
between Chebyshev polynomials and MCPs. The Chebyshev polynomials of first
kind [26] are denoted by $"T_{n}(z);\quad n=1,2,3,\dots\quad on\quad
z\in[-1,1]"$ , where

$T_{n}(z)=\cos(n\arccos(z))$ (1)

The fundamental recurrence relation is:

$T_{n}(z)=2zT_{n-1}(z)-T_{n-2}(z);\quad n\geq 2$ (2)

where $T_{0}(z)=1,\quad T_{1}(z)=z.$
The Chebyshev polynomials $"T_{n}(z)"$ can be expanded in power series:

$T_{n}(z)=\frac{n}{2}\sum_{k=0}^{\lfloor
n/2\rfloor}\frac{(-1)^{k}(n-k-1)!(2z)^{n-2k}}{(n-2k)!(k)!}$ (3)

where $\lfloor n/2\rfloor$ is the integer part of $n/2$. This explicit
representation allows us to derive many useful formulas respective to the
Chebyshev polynomials [27]

$2T_{n}(z)=\frac{1}{n+1}T^{{}^{\prime}}_{n+1}(z)-\frac{1}{n-1}T^{{}^{\prime}}_{n-1}(z)$
(4)

The Chebyshev polynomials form a complete orthogonal set on the interval
$[-1,1]$ w.r.t. the weighting function $\varpi(z)=(1-z^{2})^{-1/2}$

$\langle
T_{i},T_{j}\rangle=\int_{-1}^{1}T_{i}(z)T_{j}(z)\varpi(z)dz=\left\\{\begin{array}[]{cc}0,&i\neq
j\\\\[5.69046pt] \frac{\pi}{2},&i=j\neq 0\\\\[5.69046pt]
\pi,&i=j=0\end{array}\right.$ (5)

Consider MCPs of the first kind [28] denoted by $"Q_{n}(z)\quad;\quad
n=0,1,2,....\quad on\quad z\in[-1,1]"$ is given by:

$Q_{n}(z)=\left\\{\begin{array}[]{cc}1,&n=0\\\\[5.69046pt]
2^{1-n}T_{n}(z),&n\geq 1\end{array}\right.$ (6)

Using the previous relation, we have:

$Q_{1}(z)=z,\quad Q_{2}(z)=z^{2}-\frac{1}{2}$ (7)

The recursive formulae of MCPs as follows

$Q_{n}(z)=zQ_{n-1}(z)-\frac{1}{4}Q_{n-2}(z);\quad n\geq 3$ (8)

The MC polynomials constitutes an orthogonal basis w.r.t.$\varpi(z)$:

$\langle
Q_{i},Q_{j}\rangle=\int_{-1}^{1}Q_{i}(z)Q_{j}(z)\varpi(z)dz=\left\\{\begin{array}[]{cc}0,&i\neq
j\\\\[5.69046pt] 2^{1-2i}\pi,&i=j\neq 0\\\\[5.69046pt]
\pi,&i=j=0\end{array}\right.$ (9)

##  3 Monic Chebychev polynomials

In this section, some important relations and properties of the MCPs have
presented and introduced. These relations are essential to setup the Monic
Chebychev differentiation expansions. The next lemma introduced the recurrence
relation for MCPs. This recurrence relation is different from the recurrence
relation Chebychev polynomials. This illustrates MCPs is not a copy of
Chebychev polynomials. Also, this deviation leads us the totally different
differentiation expansion rather those in Chebychev polynomials.

Lemma 3.1

Let $Q_{n}$, be a MC of degree n, then

$Q_{n}(z)=\frac{1}{n+1}Q^{{}^{\prime}}_{n+1}(z)-\frac{1}{4(n-1)}Q^{{}^{\prime}}_{n-1}(z);\quad
n\geq 2$

Proof.

The proof is straight forward from Eq. (4) and Eq. (6).

Theorem 3.1

The m-th derivative of MCPs is expanded as follows:

$Q^{(m)}_{n}(z)=\sum_{k=0}^{\lfloor(n-m)/2\rfloor}b_{nk}c^{(m)}z^{n-2k-m}$
(10)

where

$b_{nk}=\frac{(-1)^{k}2^{-2k}n\varGamma(n-k)}{\varGamma(k+1)\varGamma(n-2k+1)}$
(11) $c^{(m)}=\sum_{l=0}^{m-1}(n-2k-l),\quad m\geq 1$ (12)

Proof.

The following formula is important for proving theorem. From eq. (3), by
replacing $z$ with $2^{1-n}z$ then MCPs $Q_{n}(z)$ can be expanded in power
series:

$Q_{n}(z)=\sum_{k=0}^{\lfloor n/2\rfloor}b_{nk}z^{n-2k},\quad n>0$ (13)

where $\lfloor n/2\rfloor$ is the integer part of $n/2$ and $b_{nk}$ as
equation (11) The first derivative of the formula (13) is the following:

$Q^{(1)}_{n}(z)=\sum_{k=0}^{\lfloor n/2\rfloor}(n-2k)b_{nk}z^{n-2k-1},\quad
n>0$ (14)

Then, using the mathematical induction to prove the theorem: For $m=1$, we
have

$Q^{(1)}_{n}(z)=\sum_{k=0}^{\lfloor(n-1)/2\rfloor}(n-2k)b_{nk}z^{n-2k-1},\quad
n>0$ (15)

which leads to the same result as the direct first differentiation. Now,
assume that theorem is satisfied at $m=j$

$Q^{(j)}_{n}(z)=\sum_{k=0}^{\lfloor(n-j)/2\rfloor}\frac{(-1)^{k}2^{-2k}n\varGamma(n-k)}{\varGamma(k+1)\varGamma(n-2k+1)}z^{n-2k-j},$
(16)

For $m=j+1$, we have

$\displaystyle Q^{(j+1)}_{n}(z)$ $\displaystyle=$
$\displaystyle\frac{d}{dz}Q^{(j)}_{n}(z)=\frac{d}{dz}\sum_{k=0}^{\lfloor(n-j)/2\rfloor}\left(\frac{(-1)^{k}2^{-2k}n\varGamma(n-k)}{\varGamma(k+1)\varGamma(n-2k+1)}z^{n-2k-j}\right)$
(17) $\displaystyle=$
$\displaystyle\sum_{k=0}^{\lfloor\frac{n-(j+1)}{2}\rfloor}\Big{(}\frac{(-1)^{k}2^{-2k}n\varGamma(n-k)}{\varGamma(k+1)\varGamma(n-2k+1)}(n-2k-j)\,c^{(j)}z^{n-2k-j-1}\Big{)}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{\lfloor\frac{n-(j+1)}{2}\rfloor}\big{(}\frac{(-1)^{k}2^{-2k}n\varGamma(n-k)}{\varGamma(k+1)\varGamma(n-2k+1)}(n-2k-j)\,c^{(j+1)}z^{n-2k-j-1}\big{)}$

Hence the theorem is satisfied for every positive integer m.

##  4 Proposed method

For solving the problems in the next section, we need to calculate the
$y(z_{i})$ such that $z_{i}$ is Chebyshev Gauss Lobatto points where
$z_{i}=\cos(\pi i/N);\quad 0\leq i\leq N$
Consider the linear HBVPs:

$\displaystyle f_{m}(z)y^{(m)}(z)$ $\displaystyle+$ $\displaystyle
f_{m-1}(z)y^{(m-1)}(z)+f_{m-2}(z)y^{(m-2)}(z)$ (18) $\displaystyle+$
$\displaystyle\dots+f_{0}(z)y(z)=g(z);a\leq x\leq b$

Subject to the boundary conditions:

$\left\\{\begin{array}[]{l}y(a)=\alpha_{0},\quad
y(b)=\beta_{0},\\\\[5.69046pt] y^{(1)}(a)=\alpha_{1},\quad
y^{(1)}(b)=\beta_{1},\\\\[5.69046pt] \vdots\\\\[5.69046pt]
y^{(q)}(a)=\alpha_{q},\quad y^{(q)}(b)=\beta_{q}.\end{array}\right.$ (19)

such that $f(z)$ and $g(z)$ is a function of $z$ or constant, $\alpha_{0}$ and
$\beta_{0}$ are constants. The number of boundary conditions is equal to the
order of the problem. Consider the non-linear HBVPs as:

$y^{(m)}_{n}=g(z,y,y^{{}^{\prime}});\quad a\leq x\leq b$ (20)

Such that $g(z,y,y^{{}^{\prime}})$ is a non-linear function. For examples,
$e^{y^{(m)}(z)}$ and $(y^{(m)}(z))^{n},n>1\quad etc\dots$, Subject to the
boundary conditions as in (19). Furthermore, by approximating the unknown
function $y(z_{i})$ as

$y(z_{i})=\sum_{i=0}^{N}A_{i}Q_{i}(z_{j}),\quad z_{j}=\cos(\frac{\pi
j}{N});0\leq j\leq N$ (21)

Where $A_{i}$ denotes constants. Realizing that $Q_{i}(z_{j})$ is square
matrix. i.e. the elements of this matrix are the values of Monic Chebyshev
polynomials at $z_{j}$ The derivatives can be expressed in differential
equation as the following form:

$\frac{d^{m}y}{dz^{m}}=\sum_{i=0}^{N}A_{i}Q_{i}(z)$ (22)

Such that m is the order of the derivative. At m=0, the form (22) become
equivalent to (21). Substituting from (22) into (18) yields:

$\displaystyle f_{m}\sum_{i=0}^{N}A_{i}Q^{(m)}_{i}(z)$ $\displaystyle+$
$\displaystyle f_{m-1}\sum_{i=0}^{N}A_{i}Q^{(m-1)}_{i}(z)$ (23)
$\displaystyle+$ $\displaystyle
f_{m-2}\sum_{i=0}^{N}A_{i}Q^{(m-2)}_{i}(z)+\dots$ $\displaystyle+$
$\displaystyle f_{0}\sum_{i=0}^{N}A_{i}Q_{i}(z)=g(z);\quad a\leq x\leq b$

By shifting the interval $[a,b]$ such that, $z$ will be belongs to the
interval from -1 to 1 and shifting the boundary conditions (20) as the
following:

$\left\\{\begin{array}[]{l}\sum_{i=0}^{N}A_{i}Q_{i}(-1)=\alpha_{0},\quad\sum_{i=0}^{N}A_{i}Q_{i}(1)=\beta_{0},\\\\[5.69046pt]
\sum_{i=0}^{N}A_{i}Q^{(1)}_{i}(-1)=\alpha_{1},\quad\sum_{i=0}^{N}A_{i}Q^{(1)}_{i}(1)=\beta_{1}\\\\[5.69046pt]
\vdots\\\\[5.69046pt]
\sum_{i=0}^{N}A_{i}Q^{(q)}_{i}(-1)=\alpha_{q},\quad\sum_{i=0}^{N}A_{i}Q^{(q)}_{i}(1)=\beta_{q}.\end{array}\right.$
(24)

Substituting (21), (22) in (18)

$\displaystyle\sum_{i=0}^{N}A_{i}Q^{(m)}_{i}(z_{j})$ $\displaystyle+$
$\displaystyle f_{m-1}\sum_{i=0}^{N}A_{i}Q^{(m-1)}_{i}(z_{j})$ (25)
$\displaystyle+$ $\displaystyle
f_{m-2}\sum_{i=0}^{N}A_{i}Q^{(m-2)}_{i}(z_{j})+\dots$ $\displaystyle+$
$\displaystyle f_{0}\sum_{i=0}^{N}A_{i}Q_{i}(z_{j})=g(z_{j});\quad 0\leq j\leq
N$

then equations (25) construct a system of algebraic equations. Particularly,
linear or non-linear system of equations depend on formulation of HBVPs ”(18)
and (20)” of the problem. Hint, we shall be replacing some rows of $D_{ij}$ by
the conditions (24). In case of the non-linearity, the non-linear algebraic
system may be treated by numerous methods like newton method, secant method
and SOR- Steffensen-Newton method for finding solution of systems of nonlinear
equations [29]. $A_{i}$ can be approximated by solving the previous system.
Finally we can get the approximate solution expanded depend on MCPs.

##  5 Numerical examples

In this part, the approximation results of linear and nonlinear HBVPs are
solved by applying the previous method. The results emphasized that the method
satisfied high accuracy and efficiency.

Example 1

Non-linear fourth order BVP:

$16y^{(4)}(z)-6e^{-4y(z)}=-12(1.5+0.5z)^{-4};\quad-1\leq z\leq 1$ (26)

With boundary conditions:

$\left\\{\begin{array}[]{cc}y(-1)=0,&y(1)=ln(2),\\\\[5.69046pt]
y^{\prime}(-1)=0.5,&y^{\prime}(1)=0.25.\\\\[5.69046pt] \end{array}\right.$
(27)

with analytical solution:

$y(z)=ln(1.5+0.5z).$

Applying the method, so equation (26) and (27) can be written as:

$\displaystyle
16\sum_{i=0}^{N}A_{i}Q^{4}_{i}(z_{j})-6e^{\sum_{i=0}^{N}}A_{i}Q_{i}(z_{j})$
$\displaystyle=$ $\displaystyle-12(1.5+0.5z_{j})^{-4}$ (28) ; $\displaystyle
0\leq j\leq N$

Subject to:

$\left\\{\begin{array}[]{cc}\sum_{i=0}^{N}A_{i}Q_{ij}(-1)=0,&\sum_{i=0}^{N}A_{i}Q_{ij}(1)=ln(2),\\\\[5.69046pt]
\sum_{i=0}^{N}A_{i}Q^{\prime}_{ij}(-1)=0.5,&\sum_{i=0}^{N}A_{i}Q^{\prime}_{ij}(1)=0.25.\\\\[5.69046pt]
\end{array}\right.$ (29)

Table 1: The approximate solutions for example 1 at N=10. z |  Analytical solution |  Approximate solution [30] |  Approximate solution
---|---|---|---
-1 |  0 |  0 |  0.0002098641
-0.8 |  0.0953101798 |  0.0950147533 |  0.0953558832
-0.6 |  0.1823215568 |  0.1814496227 |  0.1823739721
-0.4 |  0.2623642645 |  0.2609546573 |  0.2624093359
-0.2 |  0.3364722366 |  0.3347370220 |  0.3364999034
0 |  0.4054651081 |  0.4036840381 |  0.4054692075
0.2 |  0.4700036292 |  0.4684459279 |  0.4699817600
0.4 |  0.5306282511 |  0.5294932609 |  0.5305817961
0.6 |  0.5877866649 |  0.5871580370 |  0.5877209165
0.8 |  0.6418538862 |  0.6416636708 |  0.6417780693
1 |  0.6931471806 |  0.6931471806 |  0.6930743766
Table 2: The AE for example 1 at N=10. z |  AE [30] | AE (Present method)
---|---|---
-1 | 0 | 2.098641e-05
-0.8 | 2.954265e-04 | 4.570342e-05
-0.6 | 8.719341e-04 | 5.241534e-05
-0.4 | 1.4096072e-03 | 4.507148e-05
-0.2 | 1.7352146e-03 | 2.766681e-05
0 | 1.7810699e-03 | 4.099373e-06
0.2 | 1.5577013e-03 | 2.186930e-05
0.4 | 1.1349902e-03 | 4.645491e-05
0.6 | 6.286279e-04 | 6.574845e-05
0.8 | 1.902154e-04 | 7.581684e-05
1 | 0 | 7.280400e-05
Table 3: The MAE for example 1 . N | MAE
---|---
10 | 7.5991e-05
12 | 3.5789e-07
14 | 1.2290e-09
16 | 3.5426e-12
18 | 1.5099e-14
20 | 4.4409e-16

As shown in table (2), the MAE for the present method is more efficient
compared by the [30]. In table (3), the present method gets e-16 as the MAE at
N=20 which proved that method is accuracy.

Example 2

The following third order BVP:

$\displaystyle 8y^{(3)}(z)$ $\displaystyle-$ $\displaystyle
0.5(z+1)y(z))=(0.125(z+1)^{3}$ (30) $\displaystyle-$ $\displaystyle
0.5(z+1)^{2}-2.5(z+1)-3)e^{0.5(z+1)}$ ; $\displaystyle\quad-1\leq z\leq 1$

Subject to:

$y(-1)=y(1)=0,\quad y^{(1)}(-1)=0.5.$ (31)

With analytical solution

$y(z)=0.25(1-z^{2})e^{0.5(z+1)}.$

Appling the same algorithm as a previous example, the results are summarized
as table (4):

Table 4: The MAE for example 2 . N |  [30] | present method
---|---|---
3SHLM1 | 3SHLM2 | 3SHLM3
6 | 1.079e-06 | 9.13e-08 | 8.33e-08 | 2.563e-04
9 | 1.210e-7 | 6.80e-09 | 3.30e-09 | 2.809e-08
10 | — | — | — | 9.966e-10
12 | 2.770e-8 | 1.10e-09 | 9.53e-10 | 8.312e-13
14 | — | — | — | 4.441e-16
15 | 8.900e-09 | 2.00e-10 | 1.13e-10 | 1.110e-16

In table (4), we notice that, at N=12, the MAE is 8.312e-13 for our method but
the MAE is 9.53e-10 for [31]. The presented algorithm gets the double
precision ”e-16” as a MAE at N=15. That compression emphases the efficiency
and accuracy of the presented method.

Example 3

The following sixth order BVP:

$y^{(6)}(z)+y(z)=12z\cos(z)+30\sin(z);\quad-1\leq z\leq 1$ (32)

Subject to:

$\left\\{\begin{array}[]{c}y(\pm 1)=0,\quad y^{\prime}(\pm
1)=2\sin(1),\\\\[5.69046pt]
y^{{}^{\prime\prime}}(-1)=-y^{{}^{\prime\prime}}(1)=-4\cos(1)-2\sin(1).\\\\[5.69046pt]
\end{array}\right.$ (33)

With analytical solution:

$y(z)=(z^{2}-1)\sin(z).$

Using the same method as a previous example, the results are summarized as the
following:

Table 5: The MAE for example 3 . N | [32] | [32] | Present method
---|---|---|---
12 | 1.866e-8 | 1.883e-8 | 3.163e-09
15 | — | — | 4.235e-14
16 | 2.384e-13 | 2.394e-13 | 1.920e-14
17 | — | — | 4.996e-16
20 | 2.797e-16 | 2.797e-16 | 3.330e-16
24 | 2.797e-16 | 3.674e-16 | 6.106e-16

In table (5), we notice that, at N=12, the MAE is 3.1631e-09 for our method
but the MAE is 1.866e-8 for [32]. The presented method gets the double
precision ”e-16” as a MAE at N=17, while method in [32] gets it at N=20. That
compression emphases the efficiency and accuracy of the presented method.

Example 4

The ninth order BVP:

$512y^{(9)}(z)-y(z)=-9e^{0.5(z+1)};\quad-1\leq z\leq 1$ (34)

Subject to:

$\left\\{\begin{array}[]{cc}y(-1)=1,&y(1)=0,\\\\[5.69046pt]
y^{\prime}(-1)=0,&y^{\prime}(1)=-0.5e,\\\\[5.69046pt]
y^{{}^{\prime\prime}}(-1)=-0.25,&y^{{}^{\prime\prime}}(1)=-0.5e,\\\\[5.69046pt]
y^{{}^{\prime\prime\prime}}(-1)=-0.25,&y^{{}^{\prime\prime\prime}}(1)=-0.375e\\\\[5.69046pt]
y^{(4)}(-1)=-0.1875.\end{array}\right.$ (35)

with analytical solution:

$y(z)=(0.5-0.5z)e^{0.5(z+1)}.$

Table 6: The Aproximate solution for example 4 at N=10. z |  Analytical solution |  Approximate solution [33] |  Approximate solution
---|---|---|---
-1 |  1.000000 |  1.000000 |  1.000000
-0.8 |  0.994654 |  0.994654 |  0.994654
-0.6 |  0.977122 |  0.977122 |  0.977122
-0.4 |  0.944901 |  0.944901 |  0.944901
-0.2 |  0.895095 |  0.895095 |  0.895095
0 |  0.824361 |  0.824361 |  0.824361
0.2 |  0.728848 |  0.728848 |  0.728848
0.4 |  0.604126 |  0.604126 |  0.604126
0.6 |  0.445108 |  0.445108 |  0.445108
0.8 |  0.245960 |  0.245960 |  0.245960
1 |  0.000000 |  -6.940873e-09 |  -1.765061e-16
Table 7: The AE for example 4 at N=10. z |  AE [33] | AE (present method)
---|---|---
-1 | 0 | 1.110223e-16
-0.8 | 2.289423e-10 | 6.585842e-12
-0.6 | 4.623567e-09 | 1.205764e-10
-0.4 | 2.081388e-08 | 4.923047e-10
-0.2 | 4.783428e-08 | 1.029267e-09
0 | 7.122960e-08 | 1.397438e-09
0.2 | 7.339381e-08 | 1.321320e-09
0.4 | 4.961130e-08 | 8.449000e-10
0.6 | 1.629848e-08 | 3.074025e-10
0.8 | 1.814813e-08 | 3.313241e-11
1 | 6.940873e-08 | 1.765061e-16
Table 8: The MAE for example 4 . N | MAE
---|---
10 | 1.3974e-09
12 | 6.5798e-12
14 | 1.3212e-14
16 | 4.4409e-16
18 | 3.3307e-16
20 | 4.4409e-16

Table (7) shows the MAE for the present method is more efficient compared by
the [33]. Both tables (7) and (8) proved that the presented method is
efficient and accurate.

Example 5

The following linear 12th-order boundary value problem:

$\displaystyle 2096y^{(12)}(z)$ $\displaystyle+$ $\displaystyle
0.5(z+1)y(z)=-(120+11.5(z+1)$ (36) $\displaystyle+$ $\displaystyle
0.125(z+1)^{3})e^{0.5(z+1)};\quad-1\leq z\leq 1$

Subject to:

$\left\\{\begin{array}[]{cc}y(\pm 1)=0,&y^{\prime}(-1)=0.5,\\\\[5.69046pt]
y^{\prime}(1)=-0.5e,&y^{{}^{\prime\prime}}(-1)=0,\\\\[5.69046pt]
y^{{}^{\prime\prime}}(1)=-e,&y^{{}^{\prime\prime\prime}}(-1)=-0.375,\\\\[5.69046pt]
y^{{}^{\prime\prime\prime}}(1)=-1.125e,&y^{(4)}(-1)=-0.5,\\\\[5.69046pt]
y^{(4)}(1)=-e,&y^{(5)}(-1)=-0.46875,\\\\[5.69046pt]
y^{(5)}(1)=-0.78125e.\end{array}\right.$ (37)

With exact solution:

$y(z)=0.25(1-z^{2})e^{0.5(z+1)}.$

Table 9: The MAE for example 5. N | [34] | Present method
---|---|---
12 | 3.122E-11 | 1.711e-12
14 | 3.104E-12 | 9.103e-15
16 | 4.324E-13 | 3.330e-16
18 | 1.279E-13 | 1.110e-16
20 | 2.312E-14 | 1.110e-16
22 | 4.153E-15 | 1.665e-16

Table (9) shows the results of MAE between the present method and the method
in [34]. Since the MAE at N=16 is 3.3307e-16 while the method in [34] gets to
4.153e-15 at N=22 which proven that our present method is more efficient and
accurate.

## Conclusion

We have shown a Galerkin method expanded depend on the MCPs to solve HBVPs.
The developed method has been utilized to solve HBVPs without reducing to
system of lower order differential equations. The essential advantages of the
presented method are the simplicity and easy in apply. Four examples linear
and nonlinear of HBVPs have been solved. The obtained results of numerical
examples showed that the presented method is an accurate and efficient for the
HBVPs.

## Orcid

M. Abdelhakem *  — https://orcid.org/0000-0001-7085-1685

## References

  * [1] E. Celik and M. Bayram, Journal of the Franklin Institute 342, 1-6 (2005).
  * [2] E. Celik, Applied Mathematics and Computation 153, 13-17 (2004).
  * [3] A. Bhrawy, Proceedings of the Romanian Academy - Series A: Mathematics, Physics, Technical Sciences, Information Science 17, 39-47 (2016).
  * [4] M. Abdelkawy, S. Ezz-Eldien, and A. Z. M. Amin, Progress in Fractional Differentiation and Applications, 1, 187-200 (2015).
  * [5] M. A. Zaky, Computers $\&$ Mathematics with Applications 75, 2243-2258 (2018).
  * [6] N. H. Sweilam, A. M. Nagy, I. K. Youssef, and M. M. Mokhtar, Int. J. Appl. Comput. Math, 3, 333-345 (2017).
  * [7] W. M. Abd-Elhameed and Y. H. Youssri, Comp. Appl. Math., 37, 2897-2921 (2018).
  * [8] W. M. Abd-Elhameed and Y. H. Youssri, International Journal of Nonlinear Sciences and Numerical Simulation, (2019).
  * [9] M. Abdelhakem, Doha M. R., and M. El-kady, Journal of Scientific and Engineering Research 5, 317-321, (2018).
  * [10] K. Issa and R. B. Adeniyi, Journal of the Nigerian Mathematical Society 35, 18-24, (2016).
  * [11] A. Napoli and W. M. Abd-Elhameed, Calcolo 54, 57-76 (2017).
  * [12] P. K. Srivastava, Advances in Difference Equations, 9, 7, (2015).
  * [13] A. Napoli and W. M. Abd-Elhameed, Mediterr. J. Math. 14, 170 (2017) .
  * [14] S. Gholamtabar and N. Parandin, Indian Journal of Fundamental and Applied Life Sciences 5, 193-197 (2015).
  * [15] R. K. Lodhi and H. K. Mishra, Journal of the Nigerian Mathematical Society, 35, 257-265 (2016).
  * [16] M. El-Kady, H. Bakheet, and M. Khalil, International Journal of Mathematics $\&$ ComputationTM, 1, 95-105 (2009).
  * [17] G. Akram and I. A. Aslam, Journal of the Association of Arab Universities for Basic and Applied Sciences, 20, 61-67 (2016).
  * [18] M. Sohaib, S. Haq, S. Mukhtar, and I. Khan, Results in Physics 8, 1204-1208 (2018).
  * [19] N. Jha and L. K. Bieniasz, J Sci Comput, 64, 898-913 (2015).
  * [20] W. M. Abd-Elhameed, E. H. Doha, and M. A. Bassuony, ScientificWorldJournal, 2014 (2014).
  * [21] A. Akgul, M. Inc, and A. Kilicman, Applied Mathematics $\&$ Information Sciences 10, 755-764 (2016).
  * [22] F. A. Costabile and A. Napoli, Applied Numerical Mathematics 116, 129-140 (2017).
  * [23] C. Markett, Journal of Mathematical Analysis and Applications 421, 244-259 (2015).
  * [24] M. El-Kady and N. El-Sawy, International Center for Scientific Research and Studies, 9, 21-37 (2012).
  * [25] M. El-Kady and H. Moussa, Center for Scientific Research and Studies, 14, 23-36, (2013).
  * [26] M. Khalil, M. El Kady, and H. Bakheet, Lap Lambert Academic Publishing, (2012).
  * [27] J. Shen, T. Tang, and L.-L. Wang, Berlin Heidelberg: Springer-Verlag, 41, (2011).
  * [28] M. El-Kady and H. Moussa, International Journal of Contemporary Mathematical Sciences 6, 2261-2274 (2011).
  * [29] M. Darvishi and N. Darvishi, Journal Applied Mathematics 2, 21-27 (2012).
  * [30] G. Mustafa, M. Abbas, S. T. Ejaz, A. I. M. Ismail, and F. Khan, Journal of Computational Analysis $\&$ Applications 23, 607-623 (2017).
  * [31] U. Mohammed, R. B. Adeniyi, M. E. Semenov, M. Jiya, and A. I. Ma’ali, Journal of the Nigerian Mathematical Society 37, 1-22 (2018).
  * [32] E. Doha and W. Abd-Elhameed, Bulletin of the Malaysian Mathematical Sciences Society 37, 383-398 (2014).
  * [33] M. Sadaf and G. Akram, Journal of King Saud University - Science (2018).
  * [34] M. El-Gamel, British Journal of Mathematics $\&$ Computer Science 6, 13-23 (2015).

Mohamed Abdelhakem<EMAIL_ADDRESS>graduated from Helwan University,
Cairo, Egypt, in 1999. He received the M. Sc. degree in Computer Science, from
Helwan University, Cairo, Egypt, in 2003, and the Ph. D. degree in pure
mathematics from Helwan University, Cairo, Egypt, in 2011. He is the vice
Chairman and founder of Helwan School of Numerical Analysis in Egypt ”HSNAE”.


ContraBERT: Enhancing Code Pre-trained Models via Contrastive Learning

Shangqing Liu1,
Bozhi Wu1,
Xiaofei Xie222 Corresponding author.,
Guozhu Meng 3,
and Yang Liu1

1Nanyang Technological University, Singapore

2Singapore Management University, Singapore

3SKLOIS, Institute of Information Engineering, Chinese Academy of Sciences, China

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Large-scale pre-trained models such as CodeBERT, GraphCodeBERT have earned widespread attention from both academia and industry.
Attributed to the superior ability in code representation, they have been further applied in multiple downstream tasks such as clone detection, code search and code translation. However, it is also observed that these state-of-the-art pre-trained models are susceptible to adversarial attacks. The performance of these pre-trained models drops significantly with simple perturbations such as renaming variable names. This weakness may be inherited by their downstream models and thereby amplified at an unprecedented scale. To this end, we propose an approach namely ContraBERT that aims to improve the robustness of pre-trained models via contrastive learning. Specifically, we design nine kinds of simple and complex data augmentation operators on the programming language (PL) and natural language (NL) data to construct different variants. Furthermore, we continue to train the existing pre-trained models by masked language modeling (MLM) and contrastive pre-training task on the original samples with their augmented variants to enhance the robustness of the model. The extensive experiments demonstrate that ContraBERT can effectively improve the robustness of the existing pre-trained models. Further study also confirms that these robustness-enhanced models provide improvements as compared to original models over four popular downstream tasks.

Code Pre-trained Models, Contrastive Learning, Model Robustness

§ INTRODUCTION

§ BACKGROUND

§ APPROACH

§ EXPERIMENTAL SETUP

§ EXPERIMENTAL RESULTS

§ DISCUSSION

§ RELATED WORK

§ CONCLUSION

§ ACKNOWLEDGMENTS

We express our sincere gratitude to Mr Daya Guo from Sun Yat-sen University for his assistance. This research is partially supported by the National Research Foundation, Singapore under its the AI Singapore Programme (AISG2-RP-2020-019), the National Research Foundation, Prime Ministers Office, Singapore under its National Cybersecurity R&D Program (Award No. NRF2018NCR-NCR005-0001), NRF Investigatorship NRF-NRFI06-2020-0001, the National Research Foundation through its National Satellite of Excellence in Trustworthy Software Systems (NSOE-TSS) project under the National Cybersecurity R&D (NCR) Grant award no. NRF2018NCR-NSOE003-0001, the Ministry of Education, Singapore under its Academic Research Tier 3 (MOET32020-0004). IIE authors are supported in part by NSFC (61902395), Beijing Nova Program. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not reflect the views of the Ministry of Education, Singapore.

[1]
M. Allamanis, E. T. Barr, P. Devanbu, and C. Sutton, “A survey of machine
learning for big code and naturalness,” ACM Computing Surveys (CSUR),
vol. 51, no. 4, pp. 1–37, 2018.

[2]
“The github blog,” <https://github.blog/2018-11-08-100m-repos/>.

[3]
S. Lu, D. Guo, S. Ren, J. Huang, A. Svyatkovskiy, A. Blanco, C. Clement,
D. Drain, D. Jiang, D. Tang et al., “Codexglue: A machine learning
benchmark dataset for code understanding and generation,” arXiv
preprint arXiv:2102.04664, 2021.

[4]
M. Allamanis, M. Brockschmidt, and M. Khademi, “Learning to represent programs
with graphs,” arXiv preprint arXiv:1711.00740, 2017.

[5]
Y. Zhou, S. Liu, J. Siow, X. Du, and Y. Liu, “Devign: Effective vulnerability
identification by learning comprehensive program semantics via graph neural
networks,” arXiv preprint arXiv:1909.03496, 2019.

[6]
M. Allamanis, H. Jackson-Flux, and M. Brockschmidt, “Self-supervised bug
detection and repair,” arXiv preprint arXiv:2105.12787, 2021.

[7]
S. Liu, Y. Chen, X. Xie, J. K. Siow, and Y. Liu, “Retrieval-augmented
generation for code summarization via hybrid gnn,” in International
Conference on Learning Representations, 2020.

[8]
U. Alon, S. Brody, O. Levy, and E. Yahav, “code2seq: Generating sequences from
structured representations of code,” arXiv preprint arXiv:1808.01400,

[9]
U. Alon, M. Zilberstein, O. Levy, and E. Yahav, “code2vec: Learning
distributed representations of code,” Proceedings of the ACM on
Programming Languages, vol. 3, no. POPL, pp. 1–29, 2019.

[10]
S. Liu, X. Xie, L. Ma, J. Siow, and Y. Liu, “Graphsearchnet: Enhancing gnns
via capturing global dependency for semantic code search,” arXiv
preprint arXiv:2111.02671, 2021.

[11]
X. Gu, H. Zhang, and S. Kim, “Deep code search,” in 2018 IEEE/ACM 40th
International Conference on Software Engineering (ICSE).1em plus
0.5em minus 0.4emIEEE, 2018, pp. 933–944.

[12]
A. Svyatkovskiy, S. Lee, A. Hadjitofi, M. Riechert, J. V. Franco, and
M. Allamanis, “Fast and memory-efficient neural code completion,” in
2021 IEEE/ACM 18th International Conference on Mining Software
Repositories (MSR).1em plus 0.5em minus 0.4emIEEE, 2021, pp.

[13]
U. Alon, R. Sadaka, O. Levy, and E. Yahav, “Structural language models of
code,” in International Conference on Machine Learning.1em
plus 0.5em minus 0.4emPMLR, 2020, pp. 245–256.

[14]
F. Liu, G. Li, Y. Zhao, and Z. Jin, “Multi-task learning based pre-trained
language model for code completion,” in Proceedings of the 35th
IEEE/ACM International Conference on Automated Software Engineering, 2020,
pp. 473–485.

[15]
S. Iyer, I. Konstas, A. Cheung, and L. Zettlemoyer, “Summarizing source code
using a neural attention model,” in Proceedings of the 54th Annual
Meeting of the Association for Computational Linguistics (Volume 1: Long
Papers), 2016, pp. 2073–2083.

[16]
H. Husain, H.-H. Wu, T. Gazit, M. Allamanis, and M. Brockschmidt,
“Codesearchnet challenge: Evaluating the state of semantic code search,”
arXiv preprint arXiv:1909.09436, 2019.

[17]
R. Russell, L. Kim, L. Hamilton, T. Lazovich, J. Harer, O. Ozdemir,
P. Ellingwood, and M. McConley, “Automated vulnerability detection in source
code using deep representation learning,” in 2018 17th IEEE
international conference on machine learning and applications (ICMLA).1em plus 0.5em minus 0.4emIEEE, 2018, pp. 757–762.

[18]
A. V. M. Barone and R. Sennrich, “A parallel corpus of python functions and
documentation strings for automated code documentation and code generation,”
arXiv preprint arXiv:1707.02275, 2017.

[19]
Y. Zhou, J. K. Siow, C. Wang, S. Liu, and Y. Liu, “Spi: Automated
identification of security patches via commits,” ACM Transactions on
Software Engineering and Methodology (TOSEM), vol. 31, no. 1, pp. 1–27,

[20]
S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural
computation, vol. 9, no. 8, pp. 1735–1780, 1997.

[21]
A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with
deep convolutional neural networks,” Advances in neural information
processing systems, vol. 25, pp. 1097–1105, 2012.

[22]
S. Liu, C. Gao, S. Chen, N. L. Yiu, and Y. Liu, “Atom: Commit message
generation based on abstract syntax tree and hybrid ranking,” IEEE
Transactions on Software Engineering, 2020.

[23]
B. Wu, S. Liu, R. Feng, X. Xie, J. Siow, and S.-W. Lin, “Enhancing security
patch identification by capturing structures in commits,” IEEE
Transactions on Dependable and Secure Computing, 2022.

[24]
X. Li, S. Liu, R. Feng, G. Meng, X. Xie, K. Chen, and Y. Liu, “Transrepair:
Context-aware program repair for compilation errors,” arXiv preprint
arXiv:2210.03986, 2022.

[25]
F. Yamaguchi, N. Golde, D. Arp, and K. Rieck, “Modeling and discovering
vulnerabilities with code property graphs,” in 2014 IEEE Symposium on
Security and Privacy.1em plus 0.5em minus 0.4emIEEE, 2014, pp.

[26]
Y. Li, D. Tarlow, M. Brockschmidt, and R. Zemel, “Gated graph sequence neural
networks,” arXiv preprint arXiv:1511.05493, 2015.

[27]
M. Allamanis, “The adverse effects of code duplication in machine learning
models of code,” in Proceedings of the 2019 ACM SIGPLAN International
Symposium on New Ideas, New Paradigms, and Reflections on Programming and
Software, 2019, pp. 143–153.

[28]
A. Kanade, P. Maniatis, G. Balakrishnan, and K. Shi, “Pre-trained contextual
embedding of source code,” 2019.

[29]
Z. Feng, D. Guo, D. Tang, N. Duan, X. Feng, M. Gong, L. Shou, B. Qin, T. Liu,
D. Jiang et al., “Codebert: A pre-trained model for programming and
natural languages,” arXiv preprint arXiv:2002.08155, 2020.

[30]
D. Guo, S. Ren, S. Lu, Z. Feng, D. Tang, S. Liu, L. Zhou, N. Duan,
A. Svyatkovskiy, S. Fu et al., “Graphcodebert: Pre-training code
representations with data flow,” arXiv preprint arXiv:2009.08366,

[31]
A. Svyatkovskiy, S. K. Deng, S. Fu, and N. Sundaresan, “Intellicode compose:
Code generation using transformer,” in Proceedings of the 28th ACM
Joint Meeting on European Software Engineering Conference and Symposium on
the Foundations of Software Engineering, 2020, pp. 1433–1443.

[32]
L. Buratti, S. Pujar, M. Bornea, S. McCarley, Y. Zheng, G. Rossiello,
A. Morari, J. Laredo, V. Thost, Y. Zhuang et al., “Exploring software
naturalness through neural language models,” arXiv preprint
arXiv:2006.12641, 2020.

[33]
R.-M. Karampatsis and C. Sutton, “Scelmo: Source code embeddings from language
models,” arXiv preprint arXiv:2004.13214, 2020.

[34]
Y. Wang, W. Wang, S. Joty, and S. C. Hoi, “Codet5: Identifier-aware unified
pre-trained encoder-decoder models for code understanding and generation,”
arXiv preprint arXiv:2109.00859, 2021.

[35]
W. U. Ahmad, S. Chakraborty, B. Ray, and K.-W. Chang, “Unified pre-training
for program understanding and generation,” arXiv preprint
arXiv:2103.06333, 2021.

[36]
A. Kanade, P. Maniatis, G. Balakrishnan, and K. Shi, “Learning and evaluating
contextual embedding of source code,” in International Conference on
Machine Learning.1em plus 0.5em minus 0.4emPMLR, 2020, pp.

[37]
J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep
bidirectional transformers for language understanding,” arXiv preprint
arXiv:1810.04805, 2018.

[38]
Y. Liu, M. Ott, N. Goyal, J. Du, M. Joshi, D. Chen, O. Levy, M. Lewis,
L. Zettlemoyer, and V. Stoyanov, “Roberta: A robustly optimized bert
pretraining approach,” arXiv preprint arXiv:1907.11692, 2019.

[39]
L. Mou, G. Li, L. Zhang, T. Wang, and Z. Jin, “Convolutional neural networks
over tree structures for programming language processing,” in
Thirtieth AAAI Conference on Artificial Intelligence, 2016.

[40]
Z. Yang, J. Shi, J. He, and D. Lo, “Natural attack for pre-trained models of
code,” arXiv preprint arXiv:2201.08698, 2022.

[41]
Authors, “Contrabert: Enhancing code pre-trained models via contrastive
learning,” <https://sites.google.com/view/contrabert>.

[42]
A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez,
Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in
Advances in neural information processing systems, 2017, pp.

[43]
J. Svajlenko, J. F. Islam, I. Keivanloo, C. K. Roy, and M. M. Mia, “Towards a
big data curated benchmark of inter-project code clones,” in 2014 IEEE
International Conference on Software Maintenance and Evolution.1em
plus 0.5em minus 0.4emIEEE, 2014, pp. 476–480.

[44]
N. D. Bui, Y. Yu, and L. Jiang, “Self-supervised contrastive learning for code
retrieval and summarization via semantic-preserving transformations,” in
Proceedings of the 44th International ACM SIGIR Conference on Research
and Development in Information Retrieval, 2021, pp. 511–521.

[45]
P. Jain, A. Jain, T. Zhang, P. Abbeel, J. E. Gonzalez, and I. Stoica,
“Contrastive code representation learning,” arXiv preprint
arXiv:2007.04973, 2020.

[46]
R. Sennrich, B. Haddow, and A. Birch, “Improving neural machine translation
models with monolingual data,” arXiv preprint arXiv:1511.06709, 2015.

[47]
E. Ma, “Nlp augmentation,” https://github.com/makcedward/nlpaug, 2019.

[48]
A. v. d. Oord, Y. Li, and O. Vinyals, “Representation learning with
contrastive predictive coding,” arXiv preprint arXiv:1807.03748,

[49]
Z. Wu, Y. Xiong, S. X. Yu, and D. Lin, “Unsupervised feature learning via
non-parametric instance discrimination,” in Proceedings of the IEEE
conference on computer vision and pattern recognition, 2018, pp. 3733–3742.

[50]
K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for
unsupervised visual representation learning,” in Proceedings of the
IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp.

[51]
W. Wang, G. Li, B. Ma, X. Xia, and Z. Jin, “Detecting code clones with graph
neural network and flow-augmented abstract syntax tree,” in 2020 IEEE
27th International Conference on Software Analysis, Evolution and
Reengineering (SANER).1em plus 0.5em minus 0.4emIEEE, 2020,
pp. 261–271.

[52]
X. Chen, C. Liu, and D. Song, “Tree-to-tree neural networks for program
translation,” arXiv preprint arXiv:1802.03691, 2018.

[53]
M.-A. Lachaux, B. Roziere, L. Chanussot, and G. Lample, “Unsupervised
translation of programming languages,” arXiv preprint
arXiv:2006.03511, 2020.

[54]
M. Tufano, C. Watson, G. Bavota, M. D. Penta, M. White, and D. Poshyvanyk, “An
empirical study on learning bug-fixing patches in the wild via neural machine
translation,” ACM Transactions on Software Engineering and Methodology
(TOSEM), vol. 28, no. 4, pp. 1–29, 2019.

[55]
Z. Li, D. Zou, S. Xu, X. Ou, H. Jin, S. Wang, Z. Deng, and Y. Zhong,
“Vuldeepecker: A deep learning-based system for vulnerability detection,”
arXiv preprint arXiv:1801.01681, 2018.

[56]
D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by jointly
learning to align and translate,” arXiv preprint arXiv:1409.0473,

[57]
J. Wang and J. Zhu, “Portfolio theory of information retrieval,” in
Proceedings of the 32nd international ACM SIGIR conference on Research
and development in information retrieval, 2009, pp. 115–122.

[58]
M. M. Palatucci, D. A. Pomerleau, G. E. Hinton, and T. Mitchell, “Zero-shot
learning with semantic output codes,” 2009.

[59]
L. Van der Maaten and G. Hinton, “Visualizing data using t-sne.”
Journal of machine learning research, vol. 9, no. 11, 2008.

[60]

[61]
T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for
contrastive learning of visual representations,” in International
conference on machine learning.1em plus 0.5em minus 0.4emPMLR,
2020, pp. 1597–1607.

[62]
M. Kim, J. Tack, and S. J. Hwang, “Adversarial self-supervised contrastive
learning,” arXiv preprint arXiv:2006.07589, 2020.

[63]
A. Radford, J. W. Kim, C. Hallacy, A. Ramesh, G. Goh, S. Agarwal, G. Sastry,
A. Askell, P. Mishkin, J. Clark et al., “Learning transferable visual
models from natural language supervision,” arXiv preprint
arXiv:2103.00020, 2021.

[64]
B. Dai and D. Lin, “Contrastive learning for image captioning,” arXiv
preprint arXiv:1710.02534, 2017.

[65]
Z. Yang, Y. Cheng, Y. Liu, and M. Sun, “Reducing word omission errors in
neural machine translation: A contrastive learning approach,” 2019.

[66]
H. Fang, S. Wang, M. Zhou, J. Ding, and P. Xie, “Cert: Contrastive
self-supervised learning for language understanding,” arXiv preprint
arXiv:2005.12766, 2020.

[67]
T. Gao, X. Yao, and D. Chen, “Simcse: Simple contrastive learning of sentence
embeddings,” arXiv preprint arXiv:2104.08821, 2021.

[68]
D. Shen, M. Zheng, Y. Shen, Y. Qu, and W. Chen, “A simple but tough-to-beat
data augmentation approach for natural language understanding and
generation,” arXiv preprint arXiv:2009.13818, 2020.

[69]
Q. Chen, J. Lacomis, E. J. Schwartz, G. Neubig, B. Vasilescu, and C. L. Goues,
“Varclr: Variable semantic representation pre-training via contrastive
learning,” arXiv preprint arXiv:2112.02650, 2021.

[70]
S. Liu, Y. Li, and Y. Liu, “Commitbart: A large pre-trained model for github
commits,” arXiv preprint arXiv:2208.08100, 2022.

[71]
A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, I. Sutskever et al.,
“Language models are unsupervised multitask learners,” OpenAI blog,
vol. 1, no. 8, p. 9, 2019.

[72]
S. Srikant, S. Liu, T. Mitrovska, S. Chang, Q. Fan, G. Zhang, and U.-M.
O'Reilly, “Generating adversarial computer programs using optimized
obfuscations,” arXiv preprint arXiv:2103.11882, 2021.

[73]
N. Yefet, U. Alon, and E. Yahav, “Adversarial examples for models of code,”
Proceedings of the ACM on Programming Languages, vol. 4, no. OOPSLA,
pp. 1–30, 2020.

[74]
G. Ramakrishnan, J. Henkel, Z. Wang, A. Albarghouthi, S. Jha, and T. Reps,
“Semantic robustness of models of source code,” arXiv preprint
arXiv:2002.03043, 2020.

[75]
P. Bielik and M. Vechev, “Adversarial robustness for code,” in
International Conference on Machine Learning.1em plus 0.5em
minus 0.4emPMLR, 2020, pp. 896–907.

[76]
H. Zhang, Z. Li, G. Li, L. Ma, Y. Liu, and Z. Jin, “Generating adversarial
examples for holding robustness of source code processing models,” in
Proceedings of the AAAI Conference on Artificial Intelligence,
vol. 34, no. 01, 2020, pp. 1169–1176.