van der Waals forces control ferroelectric–antiferroelectric ordering in CuInP2S6 and CuBiP2Se6 laminar materials

We show how van der Waals (vdW) forces outcompete covalent and ionic forces to control ferroelectric ordering in CuInP2S6 nanoflakes as well as in CuInP2S6 and CuBiP2Se6 crystals.


Introduction
Of recent recognition is the general principle that van der Waals forces can compete against covalent and ionic forces to control chemical bonding. 1 A critical application of this is the stabilization of gold surfaces and nanoparticles by Au(0)-thiyl bonding rather than the Au(I)-thiolate bonding believed for decades to dominate in these systems. 2 This effect is not restricted to nanotechnology but also leads to Fe(II)-thiyl versus Fe(II)thiolate valence tautomerisation in cytochrome-P450 model compounds. 3 A classic catalytic system is the very strong binding of benzene to copper surfaces, 4 a system traditionally believed to involve Cu-C covalent bonding that instead now forms a paradigm for the understanding of strong van der Waals forces. 5 The van der Waals force is well known as being critical to the control of self assembly, 6,7 hence strongly affecting many conceived devices made from 2D materials. Here, we show that strong specic van der Waals interactions associated with so 8,9 copper ions go beyond this to control internal ferroelectric ordering in polar laminar materials containing sheets of so sulfur or selenium atoms. These forces prove to be sufficiently strong to, in certain cases, overturn the natural electrostatic forces driving ferroelectricity to make antiferroelectric structuresa somewhat rare effect that is only recently becoming understood. 10 Bulk materials of the form ABP 2 X 6 oen take on a laminar shape and are of considerable current interest as possible ultrathin materials for use as ferroelectrics in memory storage and other devices. Bulk systems investigated in this context include: CuInP 2 S 6 , 11-17 CuBiP 2 Se 6 , 18 AgBiP 2 Se 6 , 18 AgBiP 2 S 6 , 18 CuCrP 2 S 6 , 19 and CuVP 2 S 6 , 20 with very closely related ferromagnetic systems including Cr 2 Ge 2 Te 6 . [21][22][23] Specically, materials of interest of the form ABP 2 X 6 usually have X as a chalcogenide S, Se, Te, etc., and A/B as either a monovalent/trivalent or divalent/ divalent combination. Their laminar shape comes about as the materials contain [(PX 3 )-(PX 3 )] 4À anions that assemble with their P-P bonds aligned like a bed of nails that is then covered by two parallel "sheets" of X atoms. Between the sheets, the A and B metal atoms embed into nominally octahedral holes formed between the X atoms of different anions. The property of interest is that the metal atoms may sit at the centre of these octahedral holes, making paraelectric structures with hexagonal symmetry that have no dipole moment, or else they can move to a side of the octahedral hole to become close to just one of the two X planes, creating local polarization orthogonal to the layer. Throughout a layer, these dipole moments can align ferroelectrically, antiferroelectrically, or randomly, giving rise to observed intra-layer phase transitions. 18,19 If a single layer is internally ferroelectric, then adjacent layers can align either ferroelectrically or antiferroelectrically with respect to it, and it is this process with which we are primarily concerned herein. It is relevant to both bulk materials as well as to ultrathin materials, substances oen called "nanoakes". Nanoakes can be made by exfoliating bulk materials 17 and possibly also by self-assembly of individually prepared layers. In both cases, layer stacking will be inuenced by both inter-layer electrostatic forces as well as by van der Waals forces, making nanoakes susceptible to external manipulation. What are the dominant features that control this stacking? If a layer can be thought of in the same way as a laboratory bar magnet, then the electrostatic force will favour ferroelectric alignment. Does the van der Waals force also show preference?
For ABP 2 X 6 , recent experiments reveal results with quite different characteristics: CuInP 2 S 6 bulk 11-13 and nanoakes 17 are ferroelectric whereas the closely related material CuBiP 2 Se 6 forms an antiferroelectric bulk polymorph. 18 In this work, we use rst-principles calculations utilizing density-functional theory (DFT) to show how these features arise.

Methods
All calculations are performed using VASP 5.4.1, 24 where the valence electrons are separated from the core by use of projector-augmented wave pseudopotentials (PAW); 25 the used "POTCAR" les are dated: S, P, Se-06Sep2000, Cu-22Jun2005, and In, Bi-08Apr2002; the number of electrons included per atom are S, Se-6, P-5, Cu-11, In-3, and Bi-5. Test calculations performed including semicore electrons to the calculations resulted in changes to critical energy differences of at most 1 meV. The energy cut-off for the plane-wave basis functions was set at 500 eV for bulk solids but otherwise kept at the default values obtained using "PREC ¼ HIGH". The energy tolerance for the electronic structure determinations was set at 10 À7 eV to ensure accuracy. We use a k-space grids of size: 4 Â 4 Â 1 for bilayers, 4 Â 4 Â 2 for two-layer solids, and 4 Â 4 Â 1 for sixlayer solids. Test calculations involving expansion of these grids to double those sizes were shown to affect calculated energy differences by less than 1 meV. VASP performs periodic imaging in all three dimensions and so to model bilayers we introduced a vacuum region of ca. 15Å between the periodic images in the direction normal to the bilayer. To minimize interactions between these periodic images we applied dipolar correction as implemented in VASP (by setting the "IDIPOL" to "TRUE") in early calculations. However, the effect was found to be negligible (<1 meV) and so this sometimes troublesome procedure was no longer applied. As a second check on the calculation of ferroelectric interaction energies, calculations were performed on systems expanded three times in the z direction, with average energies also varying by less than 1 meV. In the related system CuInP 2 Se 6 , small changes associated with these calculation effects have been shown to have inuence on ferroelectric structure. 26 Geometry optimizations were performed for all structures, terminating when the forces on all atoms fell below 0.01 eVÅ À1 . For nanoakes, in-plane hexagonal unit-cells of dimension 6.5532Å for CuBiP 2 Se 6 (ref. 18) and 6.088Å for CuInP 2 S 6 , 12 are used, whereas all lattice vectors are fully optimized for bulk solids. The starting structures for geometry optimizations were based on the properties of adjacent layers observed in these materials, with test calculations performed for alternate structures using the PBE-D3 method not realizing lower-energy alternatives. Hence, whenever possible, hexagonal symmetry and/or inversion symmetry was enforced in the calculations, with molecular alignments chosen so that the symmetry of the structure paralleled symmetry elements in the plane-wave basis set. To further enhance numerical stability, the "ISYM ¼ 2" ag was used to enforce density-matrix symmetry.

Results and discussion
A. Ferroelectric order in CuInP 2 S 6 crystals and nanoakes Crystal structures of the bulk phase of CuInP 2 S 6 depict internally ferroelectrically ordered layers that ferroelectrically selfassemble with each other, producing ferroelectric crystals that can be exfoliated to produce ferroelectric nanoakes. [11][12][13]17 Its observed crystal structure is depicted in Fig. 1. Nominally hexagonal layers stack above each other in such a way that overall hexagonal symmetry is lost, with six symmetrically equivalent layers stacking vertically to form a unit that is periodic in the direction orthogonal to the layer plane. These six layers are named "a" to "f" in Fig. 1. While the In atoms remain close to the centres of their octahedral holes, the Cu atoms move to take on nearly trigonal-planar arrangements within a single S sheet.
No antiferroelectric arrangement of CuInP 2 S 6 has been observed, so we consider a possible form made simply by moving the Cu atoms from one S sheet in a layer to the other, allowing the In and other atoms to relax accordingly. In Fig. 1, we apply this transformation to alternate layers, converting the original ferroelectric crystal with a structure depicted as [abcdef] to the antiferroelectric structure [a 0 bc 0 de 0 f]. The same procedure could also have been applied to make an alternate antiferroelectric structure [ab 0 cd 0 ef 0 ], but this is symmetrically equivalent and hence not shown. All optimized antiferroelectric structures embody inversion symmetry linking adjacent layers, so these structures have no net dipolar polarization. By sliding adjacent layers with respect to each other, many possible variants of the shown antiferroelectric structure are conceivable. The energy differences between such structures are subtle and will be considered elsewhere. 41 Here we focus on a key feature of all antiferroelectric structures: as Fig. 1 shows, two very different types of inter-layer interactions are apparent depending on whether or not the Cu atoms face each other across the interface.
To highlight this effect, we consider possible bilayers derived from the (observed) ferroelectric and (envisaged) antiferroelectric bulk structures of CuInP 2 S 6 . The ferroelectric bilayer is named "F" and comprised of the layers [bc] extracted from the crystal. Extracting any other pair of layers, or interchanging the layer order to make [cb], produces symmetrically equivalent bilayers, so only one unique bilayer can be produced from the crystal in this manner. However, two unique antiferroelectric structures can be conceived, with one, named "A ii ", made by extracting layers [a 0 c], and the other, named "A oo ", made by extracting layers [bc 0 ] instead. The A ii bilayers are so named as they have both Cu atoms adjacent to each other on the "inside" of the bilayer, whereas the A oo bilayers have them on the "outside", as far away from each other as is possible, see Fig. 1. Fig. 1 indicates the energies DE of the different antiferroelectric bulk structures per cell per layer and bilayer structures per cell relative to that of the ferroelectric structures, obtained using two representative computational methods, (i) the PBE generalized-gradient approximation density functional, 42 and (ii) this combined with the D3 dispersion correction in its original form (PBE-D3). 28 The critical properties of these computational methods of relevance is that they similarly include covalent and ionic effects but only PBE-D3 embodies a realistic description of the van der Waals attraction term.
Both the PBE and PBE-D3 computational approaches predict that the observed ferroelectric bulk polymorph is more stable than its envisaged antiferroelectric variant by ca. 10 meV per cell per layer, a small but signicant quantity. That the two approaches predict very similar outcomes is naively suggestive that electrostatic forces rather than van der Waals forces are most important in determining the electronic structure of laminal materials.
Alternatively, results for the bilayers presents a completely different picture of the nature of the intrinsic processes involved in self-assembly. PBE predicts that the antiferroelectric A ii and A oo bilayers are each about 6 meV less stable per cell per layer than the ferroelectric bilayer F, with all structures barely stable to spontaneous layer separation. As it predicts the bulk ferroelectric phase to be more stable than alternating A oo and A ii bilayer-type arrangements by 9 meV per cell per layer, a longrange favourable ferroelectric electrostatic interaction of 3 meV can be deduced. However, PBE-D3 predicts that A ii is more stable than F by 92 meV per cell whereas A oo is less stable by 88 meV. Hence the van der Waals force is perceived as being very sensitive to the details of inter-layer arrangements, with the overall bulk structure emerging in part as a result of the way these contributions act to cancel each other out. Improved treatments of the electrostatic interactions 26 will not affect this analysis.
Next we consider how robust this bilayer result is likely to be considering possible variations and improvements in the computational methods used. Table 1 shows results for the A ii bilayer energy with respect to the F bilayer energy obtained using the above two computational methods plus 12 others. In one case, the PBE generalize-gradient functional is replaced with the HSE06 (ref. 39) hybrid functional, providing a much enhanced description of electron exchange. The small PBE preference of 5 meV per cell for the ferroelectric bilayer compared to A ii is reversed by a 15 meV preference for A ii , but adding D3 provides an extra 44 meV preference. This change is only half that (97 meV) predicted when D3 is added to PBE, but the general scenario remains unchanged. The 10 other methods utilized in Table 1 all involve variations of the PBE density functional combined with different dispersion approaches. In terms of magnitude of the effect, these approaches show quite a range with the total preference for A ii ranging from 27 to 105 meV per cella signicant range but again the qualitative effect remains preserved. Understanding these method differences is a topic for ongoing research, but a noteworthy results is the net stabilization of 81 meV predicted by revPBE-D3, 29 a method recently highlighted as one of the most robust methods of its type across all types of chemical scenarios. 43 Another is the value of 43 meV predicted by the many-body dispersion 33,34 (MBD) method, an advanced treatment of dispersion that, to treat conductors better, does not assume pairwise additivity; 44 it is also known to be widely robust to chemical scenario. 5 Lastly, we consider the robustness of the main result considering subtle effects not normally considered in energetics analyses of materials. Calculations of the zero-point energy difference between the two structures at the PBE-D3 level indicated a change of less than 1 meV. Also, calculations of the room-temperature Helmholtz free energy change 45 at the PBE-D3 level indicate corrections of just À3 meV per cell. Therefore the best-possible description from rst principles calculations available at the moment is that bilayers of CuInP 2 S 6 should be antiferroelectric.
To demonstrate that this prediction is consistent with the observation of ferroelectric nanoakes of CuInP 2 S 6 , we consider what happens when a third layer is added to an existing A ii bilayer, say by spontaneous self-assembly. As shown in Fig. 2, a third layer can be added in one of two ways, making structure A ii F, in which the new layer adds ferroelectrically to its neighbour, as well as structure A ii A oo , in which it adds antiferroelectrically. PBE-D3 calculations indicate that A ii F is more stable than A ii A oo by 88 meV, paralleling the previous result found for bilayer formation. Making a tetralayer system from  this trilayer then proceeds analogously, with PBE-D3 predicting that A ii FF is 87 meV per cell more stable than A ii FA oo . Taken together, these results suggest that the energies of antiferroelectric multi-layer systems can be usefully described simply in terms of sums of bilayer interactions. The key prediction here is that nanoakes grown under such conditions favouring sequential layer deposition increasingly develop more and more ferroelectric character as the number of layers increases. However, if instead two A ii bilayers self-assemble to make a tetralayer, then the PBE-D3 calculations predict the resultant structure, A ii A oo A ii , to be 8 meV per cell more stable than A ii FF, making the antiferroelectric structure the thermodynamically most stable one. This result is also easily understood in terms of pairwise inter-layer interactions. In Fig. 1 and 2, all results for bilayers, trilayers, and tetralayers were obtained constraining the in-plane lattice vectors to reference values. This minimizes the variables changing between calculations, allowing focus to be placed on key features such as the change in the van der Waals interactions between structures. Alternatively, the shown bulk structures are at fully relaxed lattices obtained individually for each computational method used, allowing focus to be placed on the small energy differences that in the end determine which structure is observed and which is not. Calculations of the bulk materials at the reference lattice vectors indicate energy changes of only 5 meV, however, so lattice relaxation does not have a profound effect. Also, the energies of the ferroelectric crystals at these reference geometries are 26 meV more stable compared to the antiferroelectric structures than one would expect simply by adding bilayer energies, an effect arising from long-range attractive dipole-dipole interactions present in ferroelectric structures. Proper treatment of these effects, 26 plus use of the most accurate computational methods possible, will be critical to making robust predictions for properties of specic nanoakes. For example, the calculated 8 meV energy difference between A ii FF and A ii A oo A ii tetralayers is too small and too susceptible to method variations for this to be regarded as a robust prediction of tetralayer properties.
Nevertheless, robust predictions from the calculations are that CuInP 2 S 6 should be antiferroelectric, that under kinetic control subsequent layers will add ferroelectrically to an initial antiferroelectric bilayer, and that ferroelectric CuInP 2 S 6 structures are intrinsically more stable than antiferroelectric ones for large nanoakes. Given this, then at some point during sequential layer growth, the point will come at which kinetically produced structures transform into their thermodynamically more stable counterparts. In general, the issue of kinetic verses thermodynamics control in molecular self-assembled monolayer production is of considerable currently of interest, with new experimental and computational methods recently developed 46 and reviewed. 1 In the experimentally observed nano-akes, 17 layers were exfoliated from ferroelectric bulk material. Naively, one would expect such a high-energy exfoliation process to provide ample opportunity for the most thermodynamically stable products to be formed.
Through consideration of the details of the structures of the observed 17 nanoakes, the calculated results provide a basis for understanding. These nanoakes were ca. 1 mm in size with ca.
half their area being 3 layers high, a quarter 4 layers high, and the remaining quarter 2 layers high. Trilayer structures, and indeed all structures containing an odd number of internally ferroelectrically ordered layers, must have net dipole polarization, the issue being whether the net polarization corresponds to that of just a single layer or else multiples of this up to the full ferroelectric order of the nanoake. When only a small amount of a fourth layer is added to a trilayer system, the energetics indicate that it can only add ferroelectrically. The fact that antiferroelectric structures come in two fundamental types and that these must be combined in multi-layer systems explains the experimental observations despite predictions that bilayers should be antiferroelectric.

B. Antiferroelectric order in CuBiP 2 Se 6 crystal
This material is observed at low temperature to comprise an antiferroelectric combination of ferroelectric monolayers, see Fig. 3. 18 It differs from CuInP 2 S 6 primarily through the horizontal alignment of the layers, with the alignment in CuBiP 2 -Se 6 crystals such that the hexagonal symmetry of each layer is maintained throughout the crystal. Again, six layers are present in the smallest unit cell that shows periodicity in the direction normal to the layer plane, but for CuBiP 2 Se 6 these consists of layer pairs that are replicated each three times. The layers in the observed crystal are labelled [ab 0 cd 0 ef 0 ] using the same notation as that already applied in Fig. 1 and 2. By simply transferring Cu atoms within octahedral holes, a ferroelectric alternative is envisaged as [abcdef] shown in Fig. 3; this is symmetrically equivalent to the other structure makeable in this fashion, [a 0 b 0 c 0 d 0 e 0 f]. From the observed antiferroelectric structure, two fundamentally different bilayer types can be extracted as A ii (e.g., [b 0 c] shown in Fig. 3) and as A oo (e.g., [ab 0 ] shown in Fig. 3). In contrast to bilayers built based on CuInP 2 S 6 , [bc 0 ] is inequivalent to [ab 0 ] while [a 0 b] is inequivalent to [b 0 c]; hence these other possibilities are also shown in Fig. 3 and labelled A 0 oo and A 0 ii , respectively. Also, two different types of ferroelectric bilayers can be simply constructed from the observed antiferroelectric crystal structure, with [ab] and [bc] shown as examples in Fig. 3 and named F and F 0 , respectively. While PBE-D3 calculations indicate that these ferroelectric variants have very similar energies, large differences are again predicted between antiferroelectric forms, and we focus on this effect.
The results shown in Fig. 3 for PBE and PBE-D3, supported by those obtained using the other methods shown in Table 1, indicate that variations in covalent and ionic bonding poorly discriminate between bilayer structures, whereas dispersion again has a profound effect. As for CuInP 2 S 6 , the attractive arrangement of dipoles in a ferroelectric bulk material stabilizes this form by ca. 3 meV per cell per layer, a small effect compared to the magnitude of the differential van der Waals forces. However, the signicant difference to CuInP 2 S 6 is that for it the stabilization and destabilization of adjacent bilayer interfaces present in the bulk material cancel, the (PBE-D3) stabilization of A ii by 86 meV per cell overpowers the destabilization of A oo by 61 meV per cell to provide the driving force for the formation of the observed antiferroelectric crystal of 13 meV per cell per layer.
C. The van der Waals forces that control ferroelectric ordering Table 2 explores the origins of the differential van der Waals bonding between the F, A ii , and A oo structures, listing the van der Waals contributions coming from the D3 correction to CuBiP 2 Se 6 and CuInP 2 S 6 bilayer energies. The D3 correction takes on a "pairwise additive" form involving summation of London dispersion energies between every pair of atoms in the system. 5,28,44 Hence, in principle it is possible to decompose this contribution into terms arising from constituent atoms, and here we do this considering contributions from Cu to Cu interactions, Cu to non-Cu }Cu} interactions, and non-Cu to non-Cu interactions using:  To do this, we evaluate the D3 energy for subsystems containing either only Cu atoms or else no Cu atoms and compare the results to the full D3 energy of the system. This procedure is somewhat approximate as D3 is strictly not purely pairwise additive, with the van der Waals properties of each atom determined in the context of its environment, but the results will be qualitatively indicative of the role of copper in controlling inter-layer van der Waals interactions. Table 2 shows that the total D3 contributions to bilayer energies range from near À3.6 eV per cell for CuBiP 2 Se 6 to near À4.0 eV per cell for CuBiP 2 Se 6 , the difference arising owing to the higher polarizabilities of Se and Bi compared to S and In. For CuBiP 2 Se 6 , ca. 45% of the dispersion energy arises through interactions with Cu atoms, reducing to ca. 30% for CuBiP 2 Se 6 . However, in terms of the D3 energy differences between the A ii or A oo structures and F, DE D3 , in both cases the Cu-involving contribution is ca. 125% of the whole, with the Cu À Cu contribution opposing it but having only one h of the magnitude. The direct Cu-Cu van der Waals interaction contributes only one sixth to one third of this total, with the remainder coming from Cu À Cu interactions. However, both types of Cu contributions are much more attractive across A ii interfaces than they are across A oo .
Why the copper van der Waals forces are strongest over A ii junctions and weakest over A oo ones is sketched in Fig. 4. The most important van der Waals interactions involve the Cu of one layer interacting with the Cu and the electron-rich part of its neighbour. The magnitude of such interactions clearly scale in the order A ii > F > A oo . Also in Fig. 4 the basic ionic inter-layer interactions are also sketched, clearly favouring F over A ii and A oo . It is the balance achieved between the weak ionic forces and the strong but opposing van der Waals forces that controls the ferroelectric ordering in ABP 2 X 6 laminar materials.

Conclusions
What is revealed for both ABP 2 X 6 laminar materials considered is that antiferroelectric structures intrinsically involve two different types of inter-layer interactions with signicantly different van der Waals forces, one type being much more attractive than the van der Waals force acting between ferroelectric structures while the other is much less attractive. The difference between these two van der Waals bonding effects compared to the weak dipole-dipole stabilization energy operative in ferroelectric structures then controls the properties of crystals and nanoakes.
Realizing the nature of the balance of forces controlling the structure of ABP 2 X 6 laminar materials will facilitate new synthetic strategies for designing desirable materials, as well as new measurements to exploit bilayer properties. How this works out for the bulk properties of CuInP 2 S 6 and CuBiP 2 Se has been developed elsewhere, 41 with a focus added to consider also how these factors contribute to ferroelectric, paraelectric, and antiferroeelectric ordering within a single layer (as has been observed in these and related materials 18,19 ). This occurs as external electric elds can manipulate the dipole-dipole interaction, with controlling eld strengths being determined by the differential van der Waals interaction.
In general, the role of van der Waals forces and the manifold ways they can compete with traditional ionic and covalent bonding scenarios is a topic of great current interest, 1 with the results found here for ABP 2 X 6 laminar materials providing another example. van der Waals forces not only effect selfassembly but also can control localized chemical structure. Which computational methods provide the best description of the van der Waals forces in these systems is also the subject of future work, 41 as is understanding how different methods perceive the competition with chemical forces and, as a consequence, how better methods can be developed in the future. 47

Conflicts of interest
There are no conicts to declare. Fig. 4 The primary ionic and van der Waals (vdW) interactions between layers of ABP 2 X 6 laminar materials. The ionic forces are weak but provide net stabilization of ferroelectric polymorphs compared to antiferroelectric ones, whereas the differences in inter-layer van der Waals interactions focus on the interaction of the Cu atoms with the electron-rich parts of the neighbouring bilayer, strongly stabilizing A ii bilayers while destabilizing A oo bilayers compared to ferroelectric ones.