Analyzing barium titanate TiO2 surface interactions with tert-butylphosphonic acid using density functional theory

Barium titanate (BTO) is a widely researched ferroelectric useful for energy storage. While BTO’s surface chemistry is commonly studied using density functional theory, little has been published on the TiO2 surface. Here, we determined that BTO’s surface response can be decoupled from the ferroelectric response by using a pre-optimized ferroelectric slab and allowing only the top three atomic z-layers to respond to ligand binding. Multiple favorable binding modes were identified for hydrogen, hydroxyl, water, and tert-butyl phosphonic acid on BTO’s TiO2 surface. Of these ligands, tBuPA dominates surface binding with binding energies as low as − 2.61 eV for its nine configurations.


Introduction
Barium titanate (BTO) is a ceramic material of high interest due to its perovskite atomic structure, yielding properties such as a high dielectric constant, ferroelectricity, piezoelectricity, and pyroelectricity. [1]These properties make both bulk and nanoscale BTO well suited for various applications, ranging from capacitors [2,3] to biomedical engineering [4,5] and to nonlinear optics. [6]However, the dielectric properties vary significantly with nanoparticle size, [7][8][9] making it difficult to reliably produce BTO-based materials with the desired properties and within acceptable tolerances. [10]Thus, research on how to increase the predictability of manufactured BTO properties is valuable as many applications require large, yet precise, dielectric constants.
One manufacturing strategy intended to increase the predictability and strength of BTO's dielectric constant is to coat the BTO nanoparticles in a surfactant that helps prevent agglomeration and increase BTO nanoparticle dispersion in polymer matrices. [11]These surfactants modify the BTO surface, which has been shown to alter the surface chemistry and properties of BTO and other perovskites. [12]To better understand the surface chemistry, recent work has simulated atomic-scale interactions between the barium oxide (BaO)terminated surface of BTO and common ligands present in nanoparticle manufacturing. [13]While the BaO surface is the more common of the two possible (100) termination surfaces, [14] all surface termination types are exposed to ligands in solvent during the manufacturing process, and studies of surfactant ligands interacting with the titanium dioxide (TiO 2 )-terminated surface are rare.
To investigate this, density functional theory (DFT)based simulation is a promising method that has been used to elucidate atomic-scale ligand-slab interactions and investigate ferroelectric material surface chemistry. [15]Past research has used DFT simulations to report on the interactions of the BaO surface with water, [16] ethanol, [17] and phosphonic acids. [18]Each investigation uses one of the many DFT surface chemistry simulation techniques such as constrained optimization, [17] surface mirroring, [19] and external electrical boundary conditions. [20]DFT simulations on BTO can elucidate key parameters, such as BTO atom displacements and ligand binding energy that provide valuable insight into the surface interactions of commonly used ligands.
In this work, we present DFT simulation results for the TiO 2 -terminated surface of BTO interacting with hydrogen, hydroxyl, water, and tert-butylphosphonic acid (tBuPA), one of the possible surfactants used to deagglomerate BTO nanoparticles in the manufacturing process. [21]We describe a systematic approach to constrained optimization that can be used to determine the ideal number of freely moving z-direction atom slab layers to relax when simulating the surface chemistry of a ferroelectric material.In doing so, we elaborate on some complexities associated with simulating ferroelectric materials both with and without ligand interaction, highlighting how the initial slab relaxation affects the binding of ligands and the apparent ferroelectric response in BTO.For each ligand of interest, we present favorable binding modes with their associated reaction energies as well as the ligand's impact on surface BTO atoms.We determine that all ligands are able to bind favorably to the BTO surface when exposed, but the relative binding energies suggest that tBuPA dominates surface binding interactions.A flowchart of our simulation techniques is presented in Supplementary Fig. S1.

Computational methods
All simulations were done using SeqQuest, a periodic-capable, local-orbital DFT code. [22]For all simulations, we used the Perdew, Burke, and Ernzerhof (PBE) generalized gradient functional [23] with spin polarization and optimized norm-conserving pseudopotentials [24] for all atoms.PBE is the standard for surface chemistry interactions and though this method tends to slightly overestimate BTO lattice parameters which results in a slightly softer lattice compared to experimental results, surface relaxation and surface chemistry should be well described. [17]he SeqQuest code performs atomic force minimization to determine a relaxed atomic configuration through one of the standard optimization routines available including steepest descent, accelerated steepest descent, and broyden.Simulations ran until either the relaxed position was found or the system oscillated between states that were close in energy.For similar system configurations, the lower system energy was reported.Simulations were first performed at a looser force convergence tolerance to identify each local energy minima before then tightening the convergence criteria.For most simulations, a maximum force tolerance of 0.05 eV/Å was used to refine the identified local minima configuration.This is near the practical force tolerance threshold for our chemical system.A lower force tolerance would prove significantly more computationally expensive and would be unlikely to yield different results, as shown in Supplementary Fig. S2 and Supplementary Table S1.
Some acid simulations failed to converge to 0.05 eV/Å and were only able to reach a force tolerance of 0.1 eV/Å due to higher forces within the acid molecule, as shown in Supplementary Table S7.These systems demonstrated continual rotation and sliding/bouncing between slightly different local states over many hundreds of optimization steps with negligible change in energy and force.This difficulty did not occur for previous work on the BaO-terminated surface [30] which we attribute to the higher flexibility of the TiO 2 top layer and the availability of several additional local, closely related binding sites.
We performed simulations on various slab types, all of which were constructed from the optimized cubic paraelectric unit cell in which both surfaces terminate with TiO 2 .We choose to begin with a paraelectric configuration for several reasons.First, the goal of this study is to describe relative binding energies for surface interactions separate from polarization response.Second, the experimental nanoparticle system under study is expected to operate in a temperature range that will result in a paraelectric system for the majority of the time.Finally, bare ferroelectric slabs using the boundary conditions of our calculation set tend toward paraelectric in the z-direction, with net zero polarization in z-direction, [30] so a paraelectric starting configuration is the logical first step, with the smallest array of polarization-related complications.Larger slabs were constructed by first replicating the unit cell then allowing all atoms to relax into the force optimized configuration before adding any interacting ligands.To determine the binding location and energy of ligands on the BTO surface, the ligand started on or near the TiO 2 -terminated surface.
All slabs can be specified using z-x-y coordinates with the z-direction depth in terms of layers and the x-y plane in terms of unit cell width.In the z-direction, one layer refers to either of the (100) BTO planes, TiO 2 or BaO.Thus, a structure that is one unit cell deep has three layers, a structure that is two unit cells deep has five layers, and so on.As a further example, a cube of BTO consisting of eight unit cells would be referred to as a five-layer 2 × 2 (or 5 × 2 × 2) slab, as it is five layers deep and two unit cells in width and length.For this paper, results are presented for a five-layer 1 × 1 BTO slab (13 atoms), a ninelayer 1 × 1 BTO slab (23 atoms), a five-layer 4 × 4 slab (208 atoms), and a nine-layer 4 × 4 BTO slab (368 atoms).Atom coloration is consistent across all figures: barium is green, titanium is light gray, oxygen is red, hydrogen is white, carbon is dark gray, and phosphorus is orange.

Successive layer relaxation
BTO nanoparticles can be represented by an infinite periodic slab of ferroelectric material using methods such as surface mirroring, pre-conditioned ferroelectric phase slab, and constrained optimization.We used a constrained optimization technique, wherein simulations are performed on a pre-optimized slab with a specified number of z-layers fixed in place.When simulating ferroelectric materials such as BTO, the ferroelectric response induced when all atoms of a slab are free to move is much stronger than what is observed experimentally and can be same order of magnitude as the surface interactions, or even larger. [25,26]Constrained optimization allows for the free layers at the top of the slab to move in response to an interacting ligand while the lower, fixed layers of the slab help to limit the observed polarization response.In doing so, the ferroelectric response and surface response can be decoupled to isolate the binding energy of surface interactions.
The constrained optimization approach is commonly employed when simulating BTO with most choosing to relax the top three layers of a five-, seven-, or nine-layer slab. [25,26]owever, there are no works we know of that present the rigorous comparative analysis needed to quantitatively justify the optimal number of fixed and free layers.We investigate two slab types to determine the optimal number of fixed layers when simulating the TiO 2 -terminated BTO slab: a pre-optimized five-layer 1 × 1 BTO slab and a pre-optimized nine-layer 1 × 1 BTO slab.In doing so, we are able to determine if there is a relationship between the depth of a slab and the optimal number of layers to freeze for that particular slab.Additionally, results from these small 1 × 1 slabs are extrapolated to wider slabs because we expect the ligand-induced response of the slab for individual ligands to remain similar at larger slab sizes.
First, a fully fixed pre-optimized BTO slab was simulated with a single hydrogen atom on the surface.Additional simulations consisted of successively relaxing layers, starting with the layer closest to the interacting hydrogen.The binding energy associated with the hydrogen interacting with the TiO 2 surface is calculated using Eq.(1).
The simulation set for the five-layer slab consisted of six simulations: one for each choice of relaxed layers, from zero to five.The set for the nine-layer slab also consisted of six simulations: successive relaxation down to four layers (half of the slab) and an additional simulation with all layers relaxed.
The binding orientation of hydrogen changes as more layers of the pre-optimized slab are allowed to move freely, as shown by the simulation set for the five-layer slab in Fig. 1.An equivalent figure for the nine-layer simulation set is shown in Supplementary Fig. S3.As more layers become free to move, the hydrogen bond shifts from an upright position over the oxygen to a side-angled position, moving closer to the adjacent oxygen atoms without binding to them.In all configurations, the hydrogen forms what is likely a covalent bond to a surface oxygen, as indicated by the measured bond length of 0.98 Å.
For both slab sizes, the binding energy becomes sufficiently negative once the top three slab layers are allowed to relax.This indicates a switch from an energetically unfavorable (1) process to an energetically favorable process.Based on the experimental results that indicate hydrogen-BTO surface binding is favorable, we are looking for the number of free slab layers that produces a negative binding energy.Thus, these results show that relaxing the top three slab layers will give a reliable binding energy while minimizing the ferroelectric response of the slab.Additionally, this result validates the common choice to relax the top three layers of a BTO slab. [25,26]

Slab relaxation
One challenge in the field of surface chemistry simulations is constructing a larger slab of a ferroelectric material such that the ferroelectric response is decoupled from the surface response.Wider slabs are often necessary to fit larger ligands on the surface with enough space that the ligand does not interact with neighboring ligands through simulation periodicity.
We constructed and optimized four 4 × 4 slab types and allowed each optimized slab to interact with a single hydrogen atom in order to determine which would yield the most useful simulation results going forward.For our design of the slab constructions, there are two choices for slab height, five-or nine-layers, and two choices for slab geometry, cubic or distorted.First, we constructed a five-layer 4 × 4 slab from the optimized cubic paraelectric unit cell.With all layers free, the slab relaxed into a distorted configuration, as shown in Fig. 2(a).The formation of this non-specific ferroelectricappearing configuration despite maintaining fixed cubic xy parameters during optimization hints at the complexity and difficulty of performing rigorous, useful surface chemistry calculations for this material.
Second, we constructed a nine-layer 4 × 4 slab from the optimized cubic paraelectric unit cell.This slab maintained the non-distorted cubic structure through relaxation, as shown in Fig. 2(d).Third, we constructed a "short" five-layer version of the nine-layer cubic slab by removing the bottom four layers.Then, the bottom two slab layers were held fixed as the slab relaxed to maintain the cubic structure, as shown in Fig. 2(b).Fourth, we constructed a "tall" nine-layer version of the five-layer distorted slab by placing the distorted five-layer slab above the bottom four layers of the cubic nine-layer slab.Then, the whole slab relaxed such that the distorted nature of the top layers propagated throughout the entire slab structure, as shown in Fig. 2(c).The non-interacting relaxed slabs are shown in Supplementary Fig. S4.The slab relaxation energies and hydrogen binding energies, calculated using Eq. ( 1), for all slab types are presented in Supplementary Table S2.
These results provide insight on the optimal slab type to use for future simulations as well as highlight the difficulty of decoupling the ferroelectric response from the surface response.For all interactions, hydrogen formed a 0.98 Å O-H covalent bond with a BTO surface oxygen in a side-angled position.The reaction energy was consistent between slab depths but varied significantly between the cubic and distorted slab geometries.This is likely because the distorted slab exhibits a ferroelectric response and is, as a result, significantly lower in energy compared to the cubic slab.In the cubic slab case, the hydrogen interacts with the cubic structure and induces a ferroelectric response meaning that the energy associated with both the hydrogen-surface interaction and the ferroelectric response are quantified in the calculated binding energy.In contrast, the ferroelectric response is already present in the distorted slab, as quantified in Supplementary Table S3, which helps isolate the surface response from the ferroelectric response when calculating the binding energy for the interaction.
The choice of slab geometry for surface chemistry simulations greatly impacts ligand interaction results, as seen in the significantly lower calculated binding energies for the cubic slab compared to the distorted slab.The difference in calculated binding energy is not enough to justify using the more computationally expensive nine-layer slab rather than the five-layer slab.Therefore, all further simulations presented in this work were performed using the five-layer 4 × 4 distorted slab with the top three slab z-layers free to relax.

Surface interactions with small ligands
Hydrogen (H), a neutral hydroxyl group (OH), and water (H 2 O) were simulated interacting with the TiO 2 surface of BTO.Hydrogen and hydroxyl groups are both common terminations of various acids and organics, thus serving as base cases for how larger molecules with these groups may interact with BTO.The binding energy of these groups on the BTO surface will also factor into the equation to calculate the binding energy if hydrogen atoms or hydroxyl groups dissociate onto the surface.While previous work has examined interactions with radical • OH and identified favorable binding modes, [27] this work examined the interactions with a neutral, nonradical version of OH.This is insightful because larger molecules, such as tert-butylphosphonic acid, are terminated with a neutral, rather than radical, version of OH.Water is present in some steps of the BTO nanoparticle manufacturing process, so we study it to ensure that water binding does not dominate that of the acid.Thus, investigating these small ligands interacting on the BTO surface can either directly or indirectly provide insight into how the surface is modified during the manufacturing process.
The reaction energy associated with OH and H 2 O interacting with the BTO surface is calculated using ( 2) and (3), respectively.Equation ( 2) is slightly different from Eq. ( 1) because it is necessary to compare the bound system to a system that includes nonradical OH.In this scenario, the energy difference between the slab and free water is compared to H and OH bound to the slab at an infinite distance apart.
To determine the ligand-slab binding modes, the ligands were placed near the surface of a five-layer 4 × 4 distorted slab with the top three layers free.Multiple starting locations relative to the TiO 2 surface were explored for all small ligands in ( 2) which all simulations terminated in one of the configurations shown in Fig. 3, indicating that these are the most likely ligand-TiO 2 surface binding modes.
Two local energy minima were identified for hydrogen interacting with the BTO surface.The most favorable hydrogen-surface interaction, shown in Fig. 3(b), occurs when hydrogen binds to one of the surface oxygen atoms.In this configuration, the titanium atoms near the interacting hydrogen are repelled while the oxygen atoms near the interacting hydrogen are attracted, creating a surface rumpling effect that is evidence of electronic reconfiguration within the slab. [28]In contrast, the other identified local energy minima of hydrogen at 3 Å from the surface, shown in Fig. 3(a), induces little impact on the BTO atoms.This configuration would likely result if the hydrogen was initially in a high energy state with the BTO surface, for example, if it was adjacent to a surface titanium.Previous research on the anatase surface has found similar configurations and binding energies for hydrogen. [29]hree local energy minima were identified for hydroxyl interacting with the BTO surface.In all favorable binding modes, the O-H bond in the hydroxyl group remains intact, and the hydroxyl oxygen is bound to a surface Ti.The slightly unfavorable binding mode, shown in Fig. 3(d), indicates that H is not likely to dissociate from an OH group.Additionally, the small difference in binding energy between (e) and (f) is a result of the position of hydroxyl's hydrogen atom relative to surface oxygen atoms, as shown in Supplementary Fig. S8.Further, the large difference in binding energy between Fig. 3(c) and (d) -(f) indicates that a sufficiently close OH will attempt to bind to the BTO surface though the OH may pass through a high energy state before binding to the surface titanium.
Three local energy minima were identified for water interacting with the BTO surface, all of which are favorable.The least favorable configuration, shown in Fig. 3(g), appears to be an unstable local minimum and occurs when one of the water's OH groups is aligned directly above a surface oxygen.The more favorable water binding configurations, shown in Fig. 3(h) and 3(i), occur when water adsorbs to the surface and forms a bond with the surface titanium.Further, the most favorable configuration, shown in Fig. 3(i), occurs when both water hydrogen atoms align themselves above surface oxygen atoms at approximately equal distances from the slab.In contrast, the slightly less favorable configuration, shown in Fig. 3(h), occurs when the two water hydrogen atoms are at unequal distances from the slab.Previous research on the anatase surface has found similar configurations and binding energies for water adsorption as those presented here. [29]or all small ligands, favorable binding modes were identified on the five-layer 4 × 4 distorted BTO slab with the top three layers free.An additional simulation was performed for all ligands on a five-layer 4 × 4 cubic BTO slab, shown in Supplementary Fig. S5.As expected, these results show increased calculated binding energies for the cubic slab compared to the distorted slab.Compared to the smaller five-layer 1 × 1 slab, the five-layer 4 × 4 slab provides a significant decrease in ligand surface coverage such that the ligands no longer interact through the periodic boundary.As a result, binding energies and configurations differ between slab widths, shown in Supplementary Fig. S6 and Supplementary Fig. S7 and quantified in Supplementary Table S4 and Supplementary Table S5.

tert-Butylphosphonic acid interactions
Finally, we simulated tBuPA interacting with the TiO 2 surface of BTO.tBuPA is a common surfactant used in manufacturing and is a polyprotic acid with two acidic protons.One or both protons of the hydroxyl groups on tBuPA may dissociate from the acid and bind to the surface separate from the acid.The binding energy associated with diprotonated tBuPA interacting with the TiO 2 surface is calculated with Eq. ( 4).
The binding energy associated with monoprotonated tBuPA interacting with the TiO 2 surface is calculated using Eq. ( 5).
Equation ( 5) assumes that the missing tBuPA proton binds to the TiO 2 surface of BTO at a theoretically infinite distance away.Similarly, Eq. ( 6) calculates the binding energy associated with deprotonated tBuPA interacting with the TiO 2 surface.
To determine the ligand-slab binding modes, tBuPA was placed near the surface of a five-layer 4 × 4 distorted slab with the top three layers free.We explored myriad starting locations relative to the TiO 2 surface for tBuPA, and all simulations terminated in one of the configurations shown in Fig. 4, indicating that these are the most common tBuPA-TiO 2 surface binding modes.Each configuration represents a local energy minimum and favorable binding occurs when one or more of tBuPA's oxygen atoms bind to BTO surface titanium atoms.Multiple favorable configurations show "near dissociation," in which a hydrogen dissociates from tBuPA and binds to a surface oxygen but remains close to and interacts with tBuPA.Single and double hydrogen dissociation configurations are also identified as local energy minima in which case a dissociated hydrogen is assumed to bind with a surface oxygen sufficiently far from tBuPA such that the two are not interacting.
Nine local energy minima were identified for tBuPA interacting with the BTO surface, providing several insights into binding configuration trends.First, the difference in binding energy between Fig. 4(d) and (e, g) indicates that, in some configurations, it can be more favorable for a single hydrogen to dissociate from tBuPA.Second, the difference in binding energy between Fig. 4(e-g) and (h, i) indicates that a second hydrogen dissociation always produces a higher energy state.Third, the difference in binding energy between Fig. 4(c) and (d) as well as (f) and (e, g) indicates that it is generally neutral or more favorable when a second tBuPA oxygen-surface titanium bond forms.Additionally, surface rumpling is evident in all cases where tBuPA forms a surface bond and causes BTO atoms to displace, as quantified in Supplementary Table S6.
The most favorable configuration, shown in Fig. 4(g), occurs when one hydrogen completely dissociates from tBuPA and a second hydrogen dissociates to bind to a nearby surface oxygen.This configuration allows for two of the three oxygen atoms in tBuPA to bond to surface titanium.Surface rumpling is most apparent in this configuration, as the tBuPA oxygen atoms are not also bound to hydrogen which allows for the oxygen atoms to more strongly attract the surface titanium atoms.
These results can be applied to improve upon current BTO nanoparticle manufacturing processes.There are more binding modes for tBuPA interacting with the TiO 2 surface compared to the BaO surface as presented in previous research, [30] which is expected, as the increased atom density of the TiO 2 provides more locations for tBuPA to interact in reasonable configurations.The increased number of binding configurations might indicate that a higher surface coverage can be ( 6) achieved on the TiO 2 surface of BTO.Additionally, binding on both surface termination types is more favorable for the monoprotonated and diprotonated forms of tBuPA.Thus, a neutral or slightly acidic solution of tBuPA is expected to produce a more reliable surfactant coverage of the BTO nanoparticle than a basic solution.

Summary and conclusions
We have used DFT simulations to determine how various ligands interact with the titanium dioxide-terminated surface of BTO.Constrained optimization analysis showed that the TiO 2 surface of BTO can be simulated using a five-or ninelayer slab in which the top three z-layers move freely.All ligands (H, OH, H 2 O, and tBuPA) were determined to bind favorably to the TiO 2 surface and induce some degree of surface rumpling.Additionally, we identified nine local energy minima for tBuPA interacting with the BTO surface, ranging from + 0.02 eV to -2.61 eV, which indicates that the surfactant forms reliable bonds to a BTO nanoparticle.Additionally, the lower binding energy of tBuPA compared to those of hydroxyl and water indicates that tBuPA can effectively replace these ligands on the BTO surface, an important step in the manufacturing process.
While the presented theoretical results are difficult to compare directly to experimental work due to the inability to isolate binding sites in experiment, our results do fall within the somewhat wide range of expected binding energies. [30]In addition, we were able to validate our simulation results by comparing them to the theoretical work published on similar surfaces, such as BTO's BaO surface and pure TiO 2 , interacting with all ligands of interest. [16,27,29,30]Further, we have highlighted some of the complexities that remain present after 100 years of ferroelectric surface modeling efforts. [5]Our results provide a more comprehensive understanding of ligand interaction with the TiO 2 surface of BTO, adding to the existing literature of ferroelectric surface interaction strategies.There are many routes that future work can take to build upon the existing field and what has been presented in this work.First, examining reaction pathways for each ligand presented here would help inform the most probable binding modes.Next, further analysis on other BTO surface termination types would help fully characterize the nanoparticle surface interactions with relevant ligands.Additional options for future work include examining the ligand and slab behavior using a controlled construction of specific ferroelectric phase slabs, investigating thicker slabs to confirm relative ligand binding effects, simulating extremely deep slabs to illustrate long-range slab relaxation effects, and applying of advanced methods for optimization of rotational modes to improve resolution of difficult to converge binding sites.

Figure 1 .
Figure 1.(a) Five-layer 1 × 1 cubic BTO slab configuration and (b) plot of hydrogen binding energy for each step in successive layer relaxation.The extent of relaxed layers is indicated by the bracket to the left of each slab.As more layers relax, the hydrogen atom shifts from a vertical orientation to sitting just above the surface, the oxygen atoms shift upwards toward the interacting hydrogen, and the barium atoms shift downwards.The calculated binding energy increases as more layers relax, indicating an increase in binding favorability.

Figure 2 .
Figure 2. Comparison of hydrogen binding configuration and calculated energy for each slab type: (a) five-layer 4 × 4 distorted slab (b) fivelayer 4 × 4 cubic slab (c) nine-layer 4 × 4 distorted slab (d) nine-layer 4 × 4 cubic slab.Slab distortion for (a) and (c) is most noticeable along the vertical titanium oxide slab backbone.Calculated hydrogen binding energy varies significantly between the cubic and distorted slab geometries despite similar hydrogen binding interactions.

Figure 3 .
Figure 3. Identified local energy minima for hydrogen (a, b), hydroxyl (c-f), and water (g-i) interacting with the TiO 2 surface of BTO.Other ligand-slab configurations were explored, but all relaxed to states nearly identical to these shown.