Anomalous Elastic Evolution Induced by Copper Hopping in van der Waals Ferroelectric CuInP2S6

Van der Waals (vdW) ferroelectric CuInP2S6 (CIPS) has great potential in post‐Moore's law electronics owing to their advantages of weak interlayer interaction, stable surface with free dangling bonds, and robust switchable spontaneous polarization. The flexoelectric effect is demonstrated as an alternative switching method for the design of ferroelectric domains in layered vdW CIPS. However, the investigation of the correlation between the polarization and elastic properties remains indistinct. Here, an elastic evolution is explored experimentally and theoretically in layered vdW CIPS with temperature, of which Young's modulus (EYoung) is determined by analyzing the force‐indentation responses of vdW ferroelectric CIPS. Interestingly, an anomalous leap in EYoung of CIPS from 35 to 65 GPa occurs when the temperature rises across the Curie temperature TC of ≈315 K. Deep potential molecular dynamic (DPMD) simulations identify that this abnormal behavior can be attributed to temperature‐dependent change of copper distribution and the local copper dynamic hopping, intertwined with the order‐disorder ferroelectric phase transition, which is quite different from the typical decrease of EYoung in the lattice constant due to temperature increasing. The exploration provides an important reference for the analysis of coupled mechanical properties and ferroelectricity in CIPS and its applications in flexible electronics.


Introduction
Layered van der Waals (vdW) CuInP 2 S 6 (CIPS) has received much attention owing to its robust room temperature ferroelectricity (Curie temperature T C of ≈315 K), atypical switching DOI: 10.1002/aelm.2023003523][4][5][6] Its dangling-bond-free stackable structure makes it ready to physically integrate onto various substrates without interfacial issues.In particular, the diverse electronic and bandgap structures endow the layered vdW ferroelectric CIPS with more advantages in electronic and optoelectronic applications. [7]n addition to the electrical and photoelectric properties, elastic properties of materials are quite basic and essential in the design and integration of new materials into devices as well, especially in the area of structural strength and reliability. [10,11]For ferroelectric materials, the existence of polarization switching may bring magical effects on mechanical properties, for example, the enhanced super elasticity in freestanding BiFeO 3 . [12,13]For layered vdW ferroelectric materials, the complex relationship between mechanical behaviors and ferroelectricity contains more abundant physical mechanisms, especially in the field of low-dimensional physics, because the polarization switching and phase transition are restricted by the separated layered structures.Taking vdW CIPS as an example, with topography variations of CIPS caused by giant strain gradients, the out-of-plane Figure 1.a) A schematic diagram reflecting the lattice structure of CIPS, the influence of flexoelectric effects due to strain gradients on the polarization, and the influence of the redistribution of copper occupation density due to temperature on elasticity.b) Raman spectra of the transferred CIPS under 600 g mm −1 grating under a 532 nm laser.c) The "field-off" and corresponding "field-on" PFM amplitude and phase hysteresis loops at the red cross area in Figure S1 (Supporting Information).d,e) PFM phase images d) before and e) after applying ±10 V DC bias on inner/outer rectangular areas.20] polarization has been demonstrated to be mechanically switchable due to the flexoelectric effect. [14]These strain engineering results suggest the existence of coupling between strain and ferroelectric polarization, as shown in Figure 1a.However, although substantial progress has been made in the significant fundamental mechanism and practical device application in vdW CIPS, the elastic properties of this material and their analytical correlation between elastic and polar ferroelectricity have remained unexplored so far.
In this work, we measured the force-deflection response curves of vdW ferroelectric CIPS membranes with different thicknesses by atomic force microscopy (AFM) nanoindentation methods to reveal the macroscopic mechanical properties of CIPS multilayers.The elastic modulus (Young's modulus) of vdW CIPS was extracted by fitting the measured load forceindentation displacement relationship with an analytical model for the first time.Changes in ferroelectric domains of CIPS as a function of temperature were observed by piezoresponse force microscopy (PFM), and the Curie temperature (T C ) was determined.An anomalous behavior of Young's modulus was observed around transition temperature.Through the combination of elasticity observations below and above T c with the mechanical loading and atomic deep potential molecular dynamic (DPMD) simulations, we establish the atomic correlation between Young's modulus and the copper redistribution and dynamic hopping intertwined with the order-disorder ferroelectric phase transition behaviors, as shown in Figure 1a.

Results and Discussion
CIPS thin membranes were mechanically exfoliated from bulk materials and then transferred to heavily doped silicon wafers (resistivity ≈1 Ω • cm) with pre-etched cylindrical holes of 2 μm in diameter and 600 nm in depth.Prior to the CIPS transfer, the silicon wafers were pre-treated by a 15 min O 2 plasma process to ensure CIPS films bind tightly with the substrates via vdW interactions.Figure S1 (Supporting Information) shows the surface morphology and height profile of a CIPS membrane sample measured by AFM.To verify the quality of CIPS, the Raman spectra under 600 g mm −1 line density grating were measured with a 532 nm laser, as shown in Figure 1b.The characteristic peaks at 264, 315, and 375 cm −1 were caused by S-P-S mode, Cu + ions vibration, and P-P stretching, respectively.The locations of the characteristic peaks are consistent with the previously reported Raman spectra of CIPS. [15,16]To eliminate the interference of electrostatic contribution, the "field-off" PFM hysteresis loops were measured. [5,17,18]Figure 1c presents the robust roomtemperature ferroelectricity in the vdW CIPS membranes based on the typical butterfly-like ferroelectric hysteresis loop in the "field-off" PFM amplitude spectra and the apparent 180°switching of the phase signals conjointly further confirm.Besides, a stiffer probe was also applied to reduce electrostatics, and the PFM hysteresis loops were shown in Figure S2 (Supporting Information).By applying reversed ±10 V DC tip bias on box-inbox patterns (Figure 1d,e), down/up polarization directions were written separately, together with the redistribution of ferroelectric domains.
We conducted AFM nanoindentation experiments (Figure 2a) to measure the intrinsic mechanical properties of vdW CIPS on the cylindrical holes.According to the previous report, with decreasing thickness, vdW CIPS undergoes a phase transition between 90 nm and 100 nm, corresponding to the disappearance of in-plane polarization and an abrupt stiffness. [21]Therefore, the thickness of all the measured samples was no more than 100 nm to ensure that the mechanical properties are not affected by the phase transition caused by the thickness.A vertical load force was applied to suspended CIPS membranes through an AFM probe.The CIPS membranes spontaneously presented a depression with a depth of  0 in the area of the hole due to gravity (Figure 2b).As the AFM probe descended in height, the maximum depth would increase to  0 + , where  represented the indentation depth of the suspended CIPS membranes.In Consideration of the parameters of the probe and the known stiffness of the cantilever beam before the test, the force-indentation curves could be obtained for the vdW CIPS membranes with different thicknesses.The radius of the cylindrical hole was much larger than that of the AFM tip, and the suspended membrane was firmly clamped at the edge of the hole during the nanoindentation process, leading to the following relationship between the load force F and indentation displacement .The relationship can be expressed by the fixed-boundary model: [22][23][24][25] where  2D denotes the pretension of the suspended CIPS membrane, E 2D represents the 2D in-plane elastic modulus, and r = 1 μm is the radius of the hole.The q is a constant related to the Poisson ratio  = 0.185, which was calculated by density function theory (DFT) as q = (1.049− 0.15 − 0.16 2 ) −1 .
Figure 2c shows the experimental and fitting results of the forceindentation curves of CIPS membranes with different thicknesses were shown in Figure 2c.The simulation results fitted by Equation ( 1) are in great agreement with the experimental data with R 2 > 0.999.[24][25] Our results in Figure 2c also confirmed the theory that the linear interval is approximately in the range of 0-20 nm. [26]Both  2D and E 2D could be extracted during the simulation process but only the parameter E 2D is of great interest to this study.Note that the influence of the indenter tip radius was neglected during the process.We collated 30 valid data points of E 2D (Figure 2d) and found that the overall trend increased gradually with increasing thicknesses.The relationship between Young's modulus (E Young ) and E 2D [22][23][24][25] is below where t is the thickness of the CIPS membranes.Figure 2e shows the E Young of CIPS with different thicknesses, along with those of other typical 2D materials such as graphene, [24,27] InSe, [28,29] black phosphorus (BP), [30] BN, [27] and MoS 2 . [23]The E Young of vdW CIPS was calculated statistically to be 47 ± 12 GPa, which was close to E Young of an InSe thin membrane with a thickness of 79 nm. [28]We found that, when the thickness of CIPS was larger than 30 nm, its E Young showed no significant correlation with thickness, likely because the interlayer sliding was negligible. [30,31]The E Young of vdW CIPS obtained is larger than that of crystal CIPS of several microns thickness (27 GPa), [32] which is also consistent with the law that CIPS become softer above 100 nm. [21]Our results are also close to the simulated E Young of CIPS in previous reports. [4,32,33]he little difference between our experimental results and the calculated one may stem from the different occupancy of interlayer Cu in different samples. [4,32,33]n order to explore the relationship between the elasticity and ferroelectricity, we heated the CIPS membrane from room temperature 293 to 333 K, which was higher than its ferroelectric T C of ≈315 K. [1,18] The evolution of ferroelectric domains of 40 nm CIPS membrane was measured by PFM during the heating process (Figure 3; Figure S3, Supporting Information).Since the thickness is much <90 nm, ferroelectric CIPS only exhibit out-plane polarization and in-plane polarization will not be discussed.[21] During the test, the stiffness of the cantilever beam at different temperatures was calibrated by testing the standard sapphire sample to eliminate the interference of probe deformation and ensure the accuracy of the measurement data.We found that the ferroelectric domains in the red dotted ellipse almost vanished during the ferroelectric transition (see Figure 3c,d), and this is consistent with the previous reports.[1,18] Besides, almost all ferroelectric domains disappeared when the temperature was equal to or >328 K, which was also consistent with the previous reports.[18] Meanwhile, AFM nanoindentation tests on CIPS were executed at different temperatures, and the corresponding results stimulated by Equation (1) were shown in Figure 4a.Interestingly, a significant variation can be obviously observed below and above T C . Figure 4b shows the temperaturedependent E Young calculated using Equation (2).A sharp rise of E Young from 35 to 65 GPa was observed between 313 and 323 K, corresponding exactly to the transition through T C .Similar re-sults were replicated for 50 nm CIPS membrane, as shown in Figure S4 (Supporting Information).Significantly, for most materials, the E Young decreases with increasing temperature, [34][35][36][37] which is completely different from the temperature-dependent behavior of E Young we have observed in CIPS.
The coincidence of temperature results suggests a possible correlation between ferroelectric phase transition and abnormal increase of E Young .To reveal the microscopic origin of such abnormal behavior, deep potential molecular dynamic (DPMD) [38][39][40][41] simulations were performed to study the thermally dynamic evolution of the mechanical properties with variation in structure parameters and the behavior of ferroelectric phase transition below and above T c . Figure 5a shows the brief schematic illustration of a deep neural network-based potential of CIPS learned from first-principles DFT calculations data (for details see methods).A deep potential model trained from DFT dataset allows us to directly compare deep potential results with DFT results.We compared the energies and atomic forces calculated using the   We explored the thermally dynamic evolution of bulk-like CIPS by running NPT (constant-pressure, constant-temperature) DPMD simulations from 50 to 500 K.To avoid the supercell size effect, we simulated the temperature-driven phase transition starting with a periodic 80 000-atom ferroelectric phase supercell consisting of 20 × 20 × 10 unit cell.Below T C , the polarization of CIPS was simulated to be 3.3 μC cm −2 , which is close to previous reports. [3,5,6,15,18,42,43]Above T C , CIPS underwent a ferroelectric to paraelectric phase transition, resulting in changes in both lattice constant and distribution of copper ions, [42,44] together with the remnant polarization gradually decreasing to zero, as the DPMD simulations predicted in Figure 5c.Since the remnant polarization has been approximately reduced to zero, it is unlikely that the sudden increase in E Young of CIPS is caused by the phase transition from low-polarization phase to high-polarization phase as mentioned in the previous report. [32]Figure 5d shows the evolution of the lattice constant, which is consistent well with previous experimental observations. [44]In general, an increase in lattice constant (see Figure 5d, around T C ) means an increase in the spacing between atoms and a volume expansion, resulting in a decrease in E Young . [45]However, the situation in our experiment is quite the reverse.The behavior of E Young in CIPS presents an abnormally abrupt enhancement along with the sudden increase in lattice constant around T C .Besides, the interlayer gap width and layer thickness of CIPS have also been proven to increase with increasing temperature, [42] which rules out the possibility of an increase in E Young due to the monoclinic to tetragonal phase transition. [21]Thus, after excluding the possibility of lattice structure change increasing E Young , the sharp stiffening at T C was most likely associated with the change of copper ions.
DPMD simulations showed that the occupancy density of Cu ions varies with the temperature and has an obvious qualitative change near the T C (Figure 5e), which is consistent with the previous report. [42]Figure 5f shows that the simulated temperaturedependent E Young increases abruptly from 34 to 68 GPa around T c , corresponding to the transition from ferroelectric to paraelectric phase, which was in good agreement with our experimental results.The simulated results of 34 and 68 GPa are close to previously reported for both low-polarization and highpolarization phases of CIPS. [32]However, it has been clearly demonstrated that a larger c-lattice parameter will decrease the energy of low-polarization phase and increase the energy of highpolarization phase, resulting in a transition from the high-to the low-polarization phase. [6]Figure 5d shows that the c-lattice parameter increases with increasing temperature, which eliminates the appearance of more high-polarization phase, so the anomalously large increase of E Young of CIPS is independent of highpolarization phase.The schematic diagrams of simulated lattice structure of ordered (ferroelectric) and disordered (paraelectric) CIPS were also presented in Figure 5f.In addition, we calculated the temporal-spatial distribution of copper ions, as shown in Figure S6 (Supporting Information), simulating the trajectory of copper ions under a temperature of 360 K within 50 ps.It is found that copper ions hopped up and down locally within the sulfur octahedron in each unit cell, rather than being locally stationary.Figure 5f compared two kinds of configurations, one with a disordered distribution of copper ions in space coordinates but without dynamical hopping at 0 K (dotted gray line), and the other with spatial-temporal disordered distribution (solid red line).The results showed that pure spatial disorder does not cause any abnormal behavior of mechanical property.Whereas, once the dynamical hopping of copper ions is introduced, the E Young arises abruptly, suggesting that the stiffness of layered CIPS above T c stems from the stiffening effect of the local dynamic hopping copper ions.Note that due to the setting value of initial energy, there was a deviation between the values of T C in the simulation and experiment, but this did not affect our analysis of mechanical and electrical properties in different phases.

Conclusion
In summary, we reported an anomalous evolution of E Young in vdW CIPS with temperature.We measured the force-indentation responses of vdW CIPS under AFM by applying a point force to the suspended membranes.Through fitting the forceindentation responses based on the fixed-boundary model, the E Young of CIPS at room temperature was calculated statistically.Furthermore, an anomalously large increase of E Young of CIPS from 35 to 65 GPa was observed around T C, which could be possibly attributed to the local dynamic hopping of copper ions and intertwined with the order-disorder ferroelectric phase transition.The correlation between elasticity and ferroelectric polarization in CIPS opens a direct window into mechanically engineering polarization states and enriches the functionalities of vdW ferroelectric materials in flexible electronics in future.

Experimental Section
Sample Preparation: High-quality CIPS single crystals were purchased from 6Carbon Technology.2D vdW layered CIPS was prepared by the mechanical exfoliation method from bulk samples and transferred onto heavily-doped p-type silicon wafers with pre-etched holes of 2 μm in diameter and 600 nm in depth.
Raman Spectroscopy and AFM&PFM Measurements: Raman spectra were performed using an excitation laser of 532 nm under 600 g mm −1 line density grating (Labram HR Evolution).All the Raman tests were carried out in the darkroom to shield ambient light interferences.The AFM&PFM measurements were performed on a Bruker Dimension Icon.The SCANANSYST-AIR probes were used for surface topography and height tests.The MESP-RC-V2 probes were used for PFM tests.The TESPA-V2 probes were used for nanoindentation tests.The "Field-off" hysteresis loops and the corresponding "Field-on" hysteresis loops were measured by Asylum Research MFP-3D with the tip HQ:NSC18/Pt.

Deep Potential Molecular Dynamic (DPMD) Simulations:
The DPMD simulations on the 2D vdW CIPS were developed based on a deep learning method trained on the data from first-principles DFT calculations.The temperature-driven ferroelectric-to-paraelectric phase transition of CIPS could be well captured by DPMD simulations, and the nature of the phase transition was characterized at high atomic-level accuracy.Using the recently proposed deep potential method, the study trained a deep neural network-based potential of CIPS learned from first-principles density DFT calculations data.Choosing the configurations for training dataset was crucial for an accurate deep potential model. [46]As the focus was on the study of the structure and its dynamical response at finite temperature.Manually constructing all the relevant configurations was tedious, especially the ones that were far from the equilibrium state.Here, a concurrent learning procedure was used to automatically generate training data that cover all the relevant configuration space. [47]The workflow of each iteration in concurrent learning procedure includes three main steps: To obtain an accurate deep potential model for layered CIPS, the vdW correction with optB86b functional was employed in structurerelated calculations. [48]An energy cutoff of 600 eV and 3 × 3 × 1 k-point mesh are sufficient to converge the energy and atomic force.

Figure 2 .
Figure 2. a) Schematic of the nanoindentation test on a suspended CIPS membrane.b) 3D and corresponding 2D AFM topographies of the CIPS membrane on the hole.c) Force-indentation curves of suspended CIPS membranes with different thicknesses.d) Experimental E 2D and e) E Young of CIPS as a function of the thickness.

Figure 3 .
Figure 3. Evolution of PFM phase of the vdW CIPS with temperatures ranging from 293 to 333 K.

Figure 4 .
Figure 4. a) Force-indentation curves of a 40 nm thick suspended CIPS membrane at different temperatures.The hollow circles represented the experimental data, and the solid lines represented the simulation results.b) Corresponding E Young of CIPS as a function of temperature.

Figure 5 .
Figure 5. a) Schematic illustration of the DPMD simulations.b) Comparison of energies calculated by DPMD and DFT.c) Simulated polarization of CIPS as a function of temperature.d) Simulated lattice constants along a-, b-, and c-axes as a function of temperature.e) Cu ion occupancy density as a function of temperature.f) Simulated E Young of CIPS as a function of the temperature and the diagrams of ferroelectric and paraelectric CIPS.
deep potential model and DFT calculations all configurations in the final training dataset.As shown in Figure 5b and Figure S5 (Supporting Information), the mean absolute error of energy (|ΔE DP-DFT |) and atomic force (|ΔF DP-DFT |) between DP and DFT is 1.1 meV/atoms and 0.7 eV Å −1 , demonstrating the good agreement between trained deep potential and DFT results.Hence, the DP model can accurately predict the thermodynamic and kinetic properties of CIPS at finite temperature with DFT accuracy.

1 )
Training the deep potential model, 2) exploring configurations by running NPT DPMD simulations at various temperatures, 3) labeling configurations according to certain criteria and adding them to the training dataset, and then repeat the step 1 again.The initial training dataset contained 1000 randomly perturbed 4 × 4 supercells containing 160 atoms.Iterating the above procedure 36 times generates ≈22 000 training configurations dataset.The DP models were trained from the above configurations and corresponding PBE-based DFT energies with (240, 240, and 240) fitting networks.