Article pubs.acs.org/cm
Influence of Structural Fluctuations, Proton Transfer, and Electric Field on Polarization Switching of Supported Two-Dimensional Hydrogen-Bonded Oxocarbon Monolayers Shuang Chen,†,§ Axel Enders,‡,§ and Xiao Cheng Zeng*,†,§ †
Department of Chemistry, ‡Department of Physics and Astronomy, and §Nebraska Center of Materials and Nanoscience, University of Nebraska−Lincoln, Lincoln, Nebraska 68588, United States S Supporting Information *
ABSTRACT: The structural alignment, proton transfer, and molecular dipole under an electric field and as a function of simulation time have been investigated computationally for experimentally observed twodimensional sheets of croconic acid (CA) on Ag(111) surface and rhodizonic acid (RA) molecules on Au(111) surface at room temperature. Depending on their local environment, some of the OH···O bonds in the CA monolayer exhibit spontaneous proton transfer especially for those bonds that are part of a trimer unit within the hydrogen-bonding network. In stark contrast, the RA molecules exhibit little proton transfer. It is found that thermal structural fluctuations of the molecular layers translate into considerable fluctuations of the polarization vector within the film plane, and even polarization reversal, at room temperature, which even can mask additional contributions to the polarization from the spontaneous and electric field induced proton transfer in CA monolayer. A common feature for both supported CA and RA monolayers is their constant polarization normal to the film plane.
■
INTRODUCTION From conventional industrial production to current nanotechnology application,1,2 the demands for commercial devices of inorganic ferroelectric materials are rising. A key property of the ferroelectric materials is their electric polarization that can be reversed under the influence of an external electric field. Ferroelectricity typically occurs together with pyroelectric and piezoelectric behaviors, which makes inorganic ferroelectric materials attractive for applications in a variety of devices, including random access memories (RAMs), capacitors, energy harvesters, sensors, and gate insulators.1−5 By comparison to inorganic ferroelectrics, their organic counterparts have additional benefits. They are light, bendable, do not require structural ordering, and can thus be fabricated cheap and on any surface, such as in printing processes. Popular examples of organic ferroelectrics include organic solids with polar components (e.g., thiourea6 and vinylidene fluoride7,8) and nonpolar organic materials such as charge-transfer (CT) complexes (e.g., mixed-stack tetrathiafulvalene (TTF) and halogenated quinones9−11), hydrogen-bonded cocrystals (e.g., phenazine (Phz) and halogenanilic acids (H2xa)12−14), singlecomponent low-molecular-mass molecular crystals,15−17 and hydrogen-bonded CT complexes.18,19 Presently, most applications of organic ferroelectrics are based on poly(vinylidene fluorid) (PVDF) and its related copolymers for ferroelectric field-effect transistors,20,21 nanolithography,22 and solar cells.23 Recent research interest is shifting toward surface-supported single-component hydrogen-bonded molecular crystals with the goal to achieve two-dimensional (2D) ferroelectric monolayers.24,25 Croconic acid (CA, H2C5O5) is particularly suited to © 2015 American Chemical Society
this goal, due to its layered, sheet-like structure in the bulk, where the electric polarization is within the plane of those layers. The molecular dipole of CA itself is within the molecular plane as well (see Supporting Information Figure S1(a)). The ferroelectric switching of CA crystals is realized by changing chemical-structure polarity through the π-bond switching and intermolecular proton-transfer processes under an applied electric field.15 CA belongs to one of a series of planar monocyclic oxocarbon acids, H2CnOn, where n = 3−6. The oxocarbon acids undergo the keto−enol transformation (−C O···HO−C ↔ C−OH···OC−) through resonanceassisted proton transfer. Remarkably, CA crystals exhibit roomtemperature ferroelectricity with a spontaneous polarization of ∼20 μC/cm2 under a coercive field of 11−29 kV/cm along the c axis.15 2D hydrogen-bonded sheets of CA have been prepared in an earlier study by molecular beam epitaxy on Ag(111) surface under ultrahigh vacuum at room temperature.24 Scanning tunneling microscopy (STM) measurements combined with the first-principles calculations demonstrated the possibility of cooperative proton tautomerism in these 2D sheets, where the substrate played a deterministic role for reducing the energy barrier for the proton transfer along the Hbonds.24 2D sheets of the related rhodizonic acid (RA, H2C6O6) were investigated in another study.25 While RA is hygroscopic and thus typically found in the form of RA dihydrate,26 2D crystalline layers of RA have been successfully Received: May 7, 2015 Revised: June 17, 2015 Published: June 17, 2015 4839
DOI: 10.1021/acs.chemmater.5b01717 Chem. Mater. 2015, 27, 4839−4847
Article
Chemistry of Materials
Figure 1. Initial configurations of self-assembled (a) CA monolayer on Ag(111) surface and (b) RA monolayer on Au(111) surface for AIMD simulations. In the CA monolayer, CA molecules exhibit spontaneous proton transfer, while no proton transfer occurs within the RA monolayer. In the AIMD simulations, an in-plane electric field of 1010 V/m is added along the x direction in CA slab models, along the y direction in RA slab models, or along the z or −z direction of both models. The CA and RA molecules are numbered, and their molecular axes are highlighted on the basis of their molecular dipole directions in Supporting Information Figure S1 for further statistical analyses. such, the influence of external electric field on the self-assembled CA monolayer on Ag(111) surface or RA monolayer on Au(111) surface can be investigated. The AIMD simulations are carried out using periodic boundary conditions. The slab models for CA monolayer on Ag(111) and RA monolayer on Au(111) are the same as those used by Kunkel et al.24,25 (see Figure 1). The optimized CA and RA monolayers on metal surfaces are taken as the initial configurations for the AIMD simulations. Because of the high computing cost of AIMD simulations, all atoms of metal substrates are fixed. The vacuum layer is about 20 Å in each system. The Brillouin zones for the slab models are only sampled at the Γ point for the AIMD simulations. The Tkatchenko−Scheffler (TS) dispersion correction30 is adopted to account for weak van der Waals (vdW) interactions. The dipole slab correction is added along the direction perpendicular to the substrate. The double numerical basis set plus polarization (DNP) is selected along with the effective-core potential (ECP) for treating core electrons. The global real space cutoff is set as 4.5 Å. The selfconsistent field (SCF) convergence is set to 10−6 au, and the thermal smearing with 0.005 Ha is applied to the orbital occupation to speed up convergence. Each AIMD simulation is performed in the canonical (NVT) ensemble with the temperature controlled at 298.15 K by using the Nosé−Hoover chain thermostat (Q ratio = 2.0 and chain length = 3). The time step is set as 1 fs. For CA monolayer on Ag(111), the in-plane electric field is applied along the x direction, while for RA monolayer on Au(111), the inplane electric field is along the y direction. Also, the electric fields perpendicular to the metal surfaces (both z and −z) are applied to the two model systems. We take the 3D CA crystal as a benchmark model system first to test a range of magnitude for the applied electric field in our simulations. In the AIMD simulation of the 3D CA crystal, we observe the π-bond switching and intermolecular proton-transfer processes along the c axis when the electric field reaches the magnitude of 1010 V/m, as can be seen in Supporting Information movie S1. Hence, for the computational study of the field effect on the CA or RA monolayer, an electric field of the same magnitude is likely needed. Note that this field is much stronger than the coercive field (e.g., ∼106 V/m) seen from the polarization hysteresis loops of crystals.
synthesized on Au(111) under ultrahigh vacuum.25 The molecular dipole of RA is in the molecular plane as well, as shown in Supporting Information Figure S1(b). The mechanism of electric-field induced switching of molecular dipoles in epitaxial 2D layers of CA24 and RA,25 and the role of the supporting substrates, are poorly understood at present. Thus, in this study, we investigate the field-induced polarization switching on two model systems: (1) a CA monolayer on Ag(111) surface and (2) a RA monolayer on Au(111) surface (see Figure 1). We focus on the 2D polarized states of CA and RA monolayers on metal surfaces, without or with an external electric field. In our ab initio molecular dynamics (AIMD) simulations, the external electric field is applied in the lateral (x−y) dimensions (Ex or Ey) or along the direction perpendicular to the metal surface (Ez or −Ez), as shown in Figure 1, to mimic the orientation-dependent ferroelectric switching of 2D CA sheet or RA sheet. We demonstrate that proton transfer occurs spontaneously within the CA monolayer on Ag(111) surface at room temperature.24 Our AIMD simulations further reveal a field effect on proton transfer in CA monolayer, which has not been addressed in a previous study.24 In contrast, no such proton transfer was observed in RA monolayers on Au(111), regardless of the magnitude and direction of the electric field. We will discuss specifically the role of thermally induced structural fluctuations and how they affect the switching behavior of the adsorbate layers under the electric field.
■
COMPUTATIONAL DETAILS
The AIMD simulations are performed in the framework of density functional theory (DFT) with generalized gradient approximation (GGA) for the exchange-correlation functional in the Perdew−Burke− Ernzenhof (PBE) form,27 as implemented in the DMol3 module of Materials Studio 7.0. The external electric field is treated via adding a static potential as an extra term to the Hamiltonian (in Dmol3).28,29 As 4840
DOI: 10.1021/acs.chemmater.5b01717 Chem. Mater. 2015, 27, 4839−4847
Article
Chemistry of Materials The AIMD trajectories, i.e., variations of structural configurations of the model systems with the simulation time, are recorded, from which various structural configurations of CA monolayer or RA monolayer with specific self-assembled structures can be selected. Subsequently, the Berry phase method,31−33 as implemented in VASP 5.3.5,34,35 is used to compute the polarization of selected CA monolayer or RA monolayer along the x and y directions without the metal substrate. The total spontaneous polarization is decomposed into ionic and electronic contributions, PBerry = Pion + Pe. For the polarization along the z direction, the Berry phase method is not applicable due to the presence of a vacuum layer. Hence, the charge density integral method is employed to estimate the z component of polarization and to further correct the x and y components of polarization from the Berry phase method based on the two-time charge density calculations of the supercells with and without the metal substrate, respectively. The Berry phase method, combined with charge density integral, gives us the three components of the polarization vector. Analysis of this polarization vector for selected CA and RA monolayers (taken from snapshots of AIMD trajectories) will provide new insights into the influence of thermal fluctuations and external electric field on the polarization switching of CA and RA monolayers on metal substrates. To compute the charge density integral of organic monolayer, the integral interval depends on two cutoffs in physical scale: one, organic monolayer/metal substrate interface, which is set to the height of metal atoms in the first layer plus vdW radius of metal atom, and the other, vacuum/organic monolayer interface, which is set to the z coordinate of the highest atom within the organic monolayer plus vdW radius of carbon atom. For estimation of the z component of the polarization, the effect of charge transfer between metal substrate and organic monolayer has been included. Note that the charge density integral cannot be used to estimate the polarizations along x and y directions, because the charges at the periodic boundary are difficult to locate. However, the polarization difference between one with metal substrate and one without from the two-time charge density integral calculations, ΔP = Pwith − Pwithout, can be used here to approximately correct the Berry phase results for x and y components as Px(y) = PBerry + ΔP from the contribution of charge-transfer interaction. The computed polarizations are summarized in Supporting Information Table S1. Because polarization differences are conceptually more physically meaningful than the “absolute” polarization,36 the final corrected polarization changes for different configurations with respect to the initial configurations of CA and RA monolayers are presented in Table 1. For both the polarization and the charge density calculations in VASP, the convergence condition in SCF iteration is set to be less than 10−6 eV. The PBE-D3 method is used, where the D3 refers to the Grimme’s correction37 to account for vdW interactions. The electron− ion interaction is described by the projector augmented wave (PAW) potentials38,39 with an energy cutoff of 600 eV. The Brillouin zones are sampled using a 8 × 8 × 1 k-point mesh within the Monkhorst−Pack scheme.40 To further understand the mechanism that leads to spontaneous proton transfer within supported H-bonded monolayers, smaller supercells with either CA dimer, CA trimer, or RA dimer placed on a metal substrate with three metal layers (see Supporting Information Figure S2) are also used to estimate activation energies for different proton-transfer processes. To this end, the metal substrates are fixed during the geometry optimization, and the vacuum layer is set to about 20 Å. For geometry optimization of reactants and products, 4 × 4 × 1 and 3 × 3 × 1 k-point meshes in the Monkhorst−Pack scheme were used for the CA dimer or CA trimer on Ag(111) and for the RA dimer on Au(111), respectively. The convergence conditions in energy, force, and displacement are set to be 10−5 au, 0.002 Ha/Å, and 0.005 Å, respectively. The transition state connecting the reactant to product is located via synchronous transit methods (complete LST/QST), which is further confirmed by the frequency analysis. For the latter, the sampling is at the Γ point to lower the computing costs. Last, singlepoint energies of the located stationary points are calculated at the denser 4 × 4 × 1 and 3 × 3 × 1 k-point meshes for CA and RA supercells, respectively, using the higher SCF convergence of 10−6 au.
Table 1. Polarization Changes in Various Configurations of CA and RA Monolayers on Metal Substrates with Respect to the Initial Configurationsa component of polarization (μC/cm2) AIMD time (fs) CA@Ag(111) 136 360 670 1689 CA@Ag(111) 382 637 872 1131 2029 2204 CA@Ag(111) 650 1465 1658 2113 CA@Ag(111) 655 1500 RA@Au(111) 378 777 1118 1511 1980 2218 RA@Au(111) 1755 1998 2206 2460 2772 2842 2987 3050 RA@Au(111) 1755 1998 2206 2460 2772 2842 2987 3050
Px without E 18.9 −11.5 1.5 −14.6 with Ex −8.9 −25.6 −2.5 3.8 −23.9 1.3 with Ez −16.3 13.7 1.4 −8.9 with −Ez −14.9 −23.8 without E 8.6 9.2 18.1 8.9 8.7 0.1 with Ez 9.5 −0.2 26.9 9.6 −8.7 −17.9 −8.9 −9.0 with −Ez 9.3 −0.1 −9.0 27.4 0.2 −8.9 −26.9 −0.1
Py
Pz
0.7 −0.1 0.7 −1.3
0.2 0.5 −0.1 0.2
−0.6 −30.6 0.7 0.9 −4.0 −3.4
0.4 −0.5 0.3 −0.3 −0.3 0.4
1.1 3.3 −0.2 −18.3
−0.6 0.2 −0.08 0.2
3.7 −26.2
−0.2 0.2
9.9 31.3 −4.9 −20.4 −36.2 −30.9
0.2 −0.3 −0.3 −0.04 0.00 −0.07
−15.3 −15.4 −20.8 −25.6 −10.1 −15.7 −15.4 −15.5
−0.05 0.1 −0.02 −0.1 −0.1 −0.08 0.03 0.04
−15.3 −20.3 −20.8 −25.7 −10.0 −10.2 −5.0 −10.2
−0.03 0.05 0.1 −0.3 −0.1 0.01 0.00 −0.2
a
These configurations are taken from snapshots in the AIMD simulations. The polarization values are computed using both the Berry phase method and the charge density integral.
■
RESULTS AND DISCUSSION Dynamic Structural Fluctuations in CA and RA Monolayers at Room Temperature. Models of the experimentally observed structures of CA monolayer on Ag(111)24 and RA monolayer on Au(111)25 are the starting point for our AIMD simulations. The variations of O···H distances with the simulation time for selected H-bonds in the absence of applied electric field show that the H atoms move 4841
DOI: 10.1021/acs.chemmater.5b01717 Chem. Mater. 2015, 27, 4839−4847
Article
Chemistry of Materials
Figure 2. (a) Definition of O···H distances r1−r12 in the slab model of CA monolayer on Ag(111) surface where the Ag atoms are omitted to show different hydrogen bonds more clearly. (b−d) Variations of O···H distances (r1−r12) of hydrogen bonds in CA monolayers on Ag(111) surface, in the absence or presence of electric field, as a function of the simulation time. Those peaks with height >1.5 Å can be viewed as the proton-transferred configurations.
corresponding to r2 and r3 curves within the 0.5−0.6 ps time domain). On the other hand, as shown in Figure 2b, the protons belonging to the H-bonds within a CA dimer, such as those characterized by the r1, r7 pair, r4, r10 pair, or r6, r12 pair in Figure 2a, remain static during the time period of simulation, in the absence of an external electric field. The observed proton-transfer behavior can be explained from the computed activation energies of the corresponding protontransfer processes (see Supporting Information Figure S2). A previous study has estimated the reaction barrier to be rather high, ∼130 meV/molecule (or 6 kcal/mol) for concerted proton transfer in CA dimer on Ag(111) surface.24 Here, we find that if the proton transfer occurs in two separated steps, the reaction barrier can be lowered by a factor of 2, to 2.9 kcal/ mol as shown in Supporting Information Figure S2(a) for the first step of proton transfer (ΔEa1). The two-step proton transfer in the CA dimer is suggested, for instance, by the variations of the pair of r1 and r7 as a function of the simulation time, where only r1 signifies a proton transfer and r7 does not during the simulation under the electric field of Ex in Figure 2b. The proton transfer is most likely not occurring concertedly, contrary to that suggested earlier.24 Other pairs, such as r4, r10 and r6, r12, remain nearly unchanged as well. In the CA trimer,
spontaneously along some of the bonds within the CA monolayer (see Figure 2), but not in the RA monolayer (for the unchanged O···H distance during the AIMD simulation as indicated by Supporting Information Figure S3). We distinguish hydrogen bonds in two local environments in the CA monolayer. Within the network in Figure 2a, CA molecules can be found to form either trimers (which incorporate three intermolecular bonds as separately characterized by r8, r9, and r11, for example) or dimers (which incorporate two intermolecular bonds as separately characterized by r4 and r10, for example). Without the electric field, the H atoms belonging to a trimer can easily transfer between two O atoms, as shown by the change of O···H distance as a function of the simulation time (see black curves corresponding to r2, r3, and r5 of trimer I in Figure 2c or r9 and r11 of trimer II in Figure 2d). However, r8 exhibits little change within the 2.4 ps simulation time (see Figure 2d). Note that those peaks (in Figure 2b−d) with height >1.5 Å can be viewed as the protontransferred configurations. Interestingly, multiple proton-transfer events, that is, at least two proton-transfer events, can occur within the same time interval. These events can be described by the corising peaks (of black curves in three panels of Figure 2c or d) within the same time domain (e.g., corising peaks 4842
DOI: 10.1021/acs.chemmater.5b01717 Chem. Mater. 2015, 27, 4839−4847
Article
Chemistry of Materials
Figure 3. Variations of lateral displacements (Δx and Δy) for (a) CA molecules on Ag(111) surface and (b) RA molecules on Au(111) surface, in the absence or presence of electric field, as a function of the simulation time. The definition of lateral displacements is highlighted in the insets.
Figure 4. Variations of molecular heights (h) for (a) CA molecules on Ag(111) surface and (b) RA molecules on Au(111) surface, in the absence or presence of electric field, as a function of the simulation time. The molecular height is defined as the vertical height of the center of mass of each molecule to the underlying substrate.
To further understand the difference in the proton transfer between CA dimers and CA trimers, we also compute the activation energies requested in one representative case of three-step proton transfer in CA trimer (see Supporting Information Figure S2(b)). For the first step of the proton transfer, the activation energy (ΔEa1) is about 2.1 kcal/mol, slightly lower than ΔEa1 involved for the CA dimer (Supporting Information Figure S2(a)). For the second and third steps of
three H atoms of hydrogen bonds also likely follow a stepwise, but not a concerted, transfer mechanism, as indicated by the peaks arising at different time for r2, r3, and r5 or r8, r9, and r11 (Figure 2c and d) even without the electric field. With the presence of the field, interestingly, it can be seen that three proton-transfer events can occur within the same time interval such as the corising peaks (blue curves in three panels of Figure 2c) after 0.6 ps. 4843
DOI: 10.1021/acs.chemmater.5b01717 Chem. Mater. 2015, 27, 4839−4847
Article
Chemistry of Materials
Figure 5. Variations of rotation angles (φ) for (a) CA molecules on Ag(111) surface and (b) RA molecules on Au(111) surface, in the absence or presence of electric field, as a function of the simulation time. The rotation angle is defined as the azimuthal angle of molecular axis to the x direction of the surface.
Figure 6. Variations of orientation angles (θ) for (a) CA molecules on Ag(111) surface and (b) RA molecules on Au(111) surface, in the absence or presence of electric field, as a function of the simulation time. The orientation angle is defined as the azimuthal angle of molecular axis to the normal direction of the surface (z).
6) at room temperature. As shown in Figure 3a, the CA molecules exhibit random movement from their equilibrium position in the lateral (x−y) dimensions with lateral displacements
