Anharmonic stabilization of ferrielectricity in CuInP$_2$Se$_6$

Using first-principles calculations and group-theory based models, we study the stabilization of ferrielectricity (FiE) in CuInP$_2$Se$_6$. We find that the FiE ground state is stabilized by a large anharmonic coupling between the polar mode and a fully symmetric Raman-active mode. Our results open possibilities for controlling the single-step switching barrier for polarization by tuning the Raman-active mode. We discuss the implications of our findings in the context of designing next-generation optoelectronic devices that can overcome the voltage-time dilemma.

CIPSe was reported to undergo a broad phase transition (200 to 240 K) from the paraelectric (PE) phase (space group P31c) to the FiE phase (space group P31c) via an incommensurate phase where there is coexistence of FiE and antiferroelectric order [18,19]. Similar to other materials in this family [20], the phase transition is expected to have a large order-disorder character in addition to a displacive character [21]. A coupled oscillator model with an on-site potential that describes the local distortions leading to the double-well potential and an inter-oscillator coupling term that described the ordering of dipoles between the different sublattices can together reconcile with both the displacive character and the order-disorder character in these materials [20,22]. While considerable order-disorder character is expected for the phase transition near the critical temperature, the switching under a finite field from -FiE to +FiE at low temperatures should necessarily be described by only the on-site potential. Here, we assume that the energy barrier for switching corresponds to the energy difference between the FiE phase and the PE phase.
The prior works discussed the stability of the FiE phase for various thicknesses under different boundary conditions [10]. They also obtained the on-site potential using a Landau expansion about the polar off-centering of only the Cu atoms. While the distortions of the Cu sites relative to the PE phase are large (1.2 Å) [23] it would be surprising for such a large distortion to stabilize the FiE phase without the participation from the other atoms. In fact, to understand the microscopic origin of the polarization one needs to include the distortions of the In atoms in addition to the Cu atoms [23]. As such, a better model for the onsite potential in terms of all the relevant atomic distortions is needed to understand the stabilization of the FiE phase.
In this paper, we address these questions and discuss the different structural distortions leading to the stabilization of the FiE phase. By choosing an ordered PE parent structure we discuss the structural distortions in the FiE phase using density functional theory (DFT) calculations and group-theoretical methods. We find that the polar displacements of the Cu and In atoms create most of the polarization in the FiE phase. However, this polar distortion alone does not lead to a gain in the energy, and even the full polar mode leads to a very shallow double-well potential (2 meV/f.u.). We report that a strong anharmonic coupling between the polar mode and a fully symmetric Raman-active mode is necessary to stabilize the polar phase. While such a coupling is symmetry allowed even in conventional FEs, we find that its magnitude is large in CIPSe. By analyzing the various contributions to the total energy we discuss the microscopic origin of this large coupling. We demonstrate that strain can be an effective knob to both enhance and suppress the FiE phase. Finally, we discuss the possible implications of this large nonlinear coupling in device applications.

II. METHODS
We computed the first-principles total energies using Vienna ab initio simulation package (VASP) [24], with the PBE functional including van der Waals correction as implemented by Grimme (DFT-D2) [25]. We used projector augmented wave pseudopotentials with valance configurations Cu (3d 10 4s 1 ), In (5s 2 5p 1 ), P (3s 2 3p 3 ), and Se (4s 2 4p 4 ). Structural relaxation was done with a force convergence tolerance of 0.1 meV/Å using a conjugate-gradient algorithm. The convergence criterion for the electronic self-consistent calculations was set to 10 −8 eV. A regular 8 × 8 × 4 -centered k-point grid was used to sample the Brillouin zone with a plane-wave cutoff energy of 600 eV. The soft-modes were computed using the density-functional-perturbation theory. Polarization was calculated using the Berry-phase approach [26]. While reporting the polarization and the distortion amplitudes the center of mass of the P atoms is used as the origin. This choice does not affect our results. The presented values of the Kohn-Sham eignevalues, the Ewald energy, and the 1 The additional three out-of-plane polar modes ( − 2 [A 3 ]) involving only one of the In, P, or Se atoms, and the out-of-plane breathing modes involving only the P atoms ( + 1 [A 3 ]: P) are not shown for the sake of clarify.
Hartree energy were from the values reported in the VASP output. We used the ISOTROPY software suite to aid with the group-theoretic analysis [27].

A. Crystal structure
We chose the bulk ordered high-symmetry structure of CIPSe with space group P31c [see Figs. 1(a) and 1(b)] as the reference PE structure for studying the FiE phase [29]. Apart from serving as a zero-point reference for polarization this phase also allows us to study the one-step barrier for the switching of polarization at 0 K. The computed lattice parameters and the occupied Wyckoff positions agree well with the average experimental parameters in the P31c phase (see Table I) [28]. The crystal structure is similar to the transition metal trichalcogenides with a trigonal lattice [30,31]. However, two different atoms (Cu and In) occupy the metal site. This naturally creates an inversion-asymmetry within each layer. The Cu and In ions have a nominal charge of +1 and +3, respectively, resulting in an outer-shell electronic configuration of d 10 for both the atoms. They both are octahedrally coordinated by six Se atoms. The layers are stacked such that the P 2 Se 6 ligand sits on top of the Cu atom from the other layer, thereby recovering inversion symmetry. This interlayer inversion center is labeled ("×") in Fig. 1. The PE phase is unstable with both zone-center and zoneboundary instabilities. As we are interested in understanding the FiE phase which is reported to have the same cell size as the PE phase [11], we only consider the zone-center instability here. We verified that the relaxed structures from the other instabilities have higher energy than the FiE phase.
The relaxed crystal structure of the FiE phase with space group P31c is shown in Figs. 1(c) and 1(d). The computed lattice parameters and Wyckoff positions compare well with the corresponding experimental values (see Table II) [11]. The Cu and In atoms exhibit an antiparallel out-of-plane distortion relative to each other. Further, the P and Se Wyckoff sites are split into two sites each, corresponding to the top and bottom sub-layers within each monolayer. We label the top (bottom) trigonal Se plane within each CIPSe layer as Se t (Se b ).

B. Mode decomposition of the FiE phase
The polar mode that transforms as the − 2 irreducible representation (irrep) of P31c is the primary order parameter for the transition from the PE phase to the FiE phase. There are six symmetry-adapted modes that transform as the   In addition, the polar mode induces a Raman-active mode which preserves the full symmetry of the parent space group (labeled + 1 ). This mode in turn is made up of a combination of four SAMs that involve the displacement of only the P and Se atoms. Two of these SAMs are labeled + 1 [A 3 ] and correspond to the out-of-plane expansion of the individual P and Se layers. Figure 1( Fig. 1(j)]. Table III shows the amplitude of the corresponding SAMs in the FiE phase. Notably, the largest distortion involves the out-of-plane motion of the Cu atoms towards the Se t trigonal face (labeled − 2 [A 3 ]: Cu), and has an amplitude of 1.15 Å. Such a large distortion is atypical even for FEs. It is reasonable to assume that this will lead to large anharmonicity in the phonon modes, as we demonstrate later.  In (2.24) is larger than that of Cu (0.90), a relatively small displacement of In atoms (0.18 Å) can carry a large dipole moment. This results in the reported FiE phase with partial cancellation of polarization between the two sites, and a net polarization along the direction of the Cu displacement. The fully relaxed FiE phase has a total polarization of 2.58 μC/cm 2 , which compares well with prior first-principles predictions [10]. Even though the polarization is smaller than what is typically observed in conventional FEs, it is larger than other reported ferrielectrics such as ammonium sulphate (0.62 μC/cm 2 ) [32]. Figure 2 shows the total polarization as a function of the fractional amplitude of the fully symmetric mode (Q + 1 ) and the polar mode (Q − 2 ) in the FiE phase with respect to the PE phase. These modes induce local distortions proportional to that listed in Table III. The lattice parameters are fixed to that of the fully relaxed PE lattice parameters. This does not affect our central results. It is evident that the Q + 1 on its own does not carry any polarization, as expected from symmetry. But a condensation of Q + 1 on top of the polar mode Q − 2 modifies the electronic polarization. Still, the polar mode carries most of the polarization in the system. As the metal atom distortions dominate the polar mode, we can conclude that the microscopic origin of the polarization in FiE CIPSe is from the distortion of the metal atoms. This is also corroborated by Table IV which shows that the distortions of the Cu and In atoms create most of the polarization in FiE CIPSe.

D. Energy profile of the polar mode
While the Cu and In atoms lead to most of the polarization in the FiE phase, from Fig. 3 we find that neither the individual distortion of the Cu atoms (blue lines) or In atoms (orange line), nor the collective ferrielectric distortion of both the atoms (red lines) creates a double-well potential with respect to the PE phase [33]. This highlights the inadequacy of models based on a Landau expansion of only the Cu polar offcentering to explain the FiE phase [10]. The energy profile for the full polar-instability in the PE phase (black curve in Fig. 3) shows a shallow double-well potential (2 meV). Here, while we present the total energy as a function of the polar mode along the relaxed PE-FiE direction we verified that similar results are obtained if we choose instead the polar soft-mode eigenvectors from the force-constant matrix (see Fig. S1 in the Supplemental Material [34]). This highlights the importance of the collective distortion of all atoms that transforms as the − 2 irrep of P31c rather than just the displacement of Cu (or In) atom as discussed in prior research [10,18,19], as the primary order parameter for understanding the symmetry lowering phase transition in CIPSe.
Interestingly, both the well-depth and amplitude of the polar mode at the minimum are much smaller than that of the fully relaxed FiE phase. The fully relaxed FiE phase has an energy gain of 98 meV/f.u. with respect to the PE phase. This is in sharp contrast to other proper FEs where the dominant polar distortion on its own accounts for a large portion of the energy gain in the polar phase characterized by a deep doublewell potential [35]. This large energy difference in CIPSe immediately suggests that to understand the microscopic   . The scale of the latter is the same as that in Fig. 3. All the reported energies are relative to the fully relaxed PE phase (P31c). mechanism that stabilizes the polar phase requires not only the polar mode but also the nonpolar modes.

E. Total energy surface
To explain the above mentioned discrepancy we compute the energy surface [see Fig. 4(a)] by varying fractional amplitudes of the fully symmetric mode (Q + 1 ) and the polar mode (Q − 2 ). We find that the total energy goes up with Q + 1 . This is expected as the initial configuration is completely relaxed. We also recover the shallow double-well as we condense Q − 2 mode. The lowest energy structures corresponding to the fractional coordinates of (1, ±1) have 61 meV/f.u. lower in energy than the PE phase. While further stabilization (37 meV/f.u.) of the FiE phase is provided by a full structural relaxation, Fig. 4(a) shows that the dominant factor leading to the stabilization of the FiE phase is the strong anharmonic coupling between the polar mode and the fully symmetric mode. Such anharmonicity has also been reported in other members of the thiophosphate family [36].
To quantify this anharmonic coupling, we fit the energy surface about the PE phase. The total energy (E ) can be written as where E 0 is the energy of the PE phase that is set to zero, a i j and b i j are the real-valued coefficients of the expansion about Q + 1 and Q − 2 , respectively. The free-energy expansion considered in Eq. (1) is only for the on-site interactions. Additional intersite coupling coming from the zone-boundary modes needed to considered to discuss order-disorder transition [20]. As Q s are in fractional units, the coefficients in the fit have units of energy allowing us to compare the contribution from each term. c i j s are the real-valued coefficients corresponding to the anharmonic coupling between the two modes. Anharmonicity due to coupling between multiple modes has attracted a lot of attention in the context of nonproper FEs [37][38][39][40][41] and antiferroelectrics [42,43]. However, in the case of proper FEs even though such a coupling between the polar mode and the fully symmetric mode is allowed by symmetry they are considered to be small and are typically ignored [44].
Our approach of including this coupling between the full polar mode and other structural distortions goes beyond the existing models which looks at only the Cu distortions [10]. We will show that such an anharmonic coupling is crucial to the stabilization of the polar phase, and is key to controlling the polarization in CIPSe. Table V shows the coefficients obtained from fitting the energy surface in Fig. 4(a). For the total energy to be bounded the coefficients of the largest-order polynomial for each mode should be positive. As the sixth-order coefficient of the polar mode (b 06 ) was negative we included up to the eighth-order term. The asymptotic standard error in the fit was less than 8.2%, with most coefficients within 4.5%. The fit can be further improved by adding additional terms (see Table S1 in the Supplemental Material) [34]. As the overall conclusions are independent of these additional terms, we omit them in the discussion of the minimal model.
We find b 02 to be negative and b 04 to be positive as expected of a characteristic double-well potential due to a polar-instability [45]. However, from Table V it is clear that the anharmonic coupling is anomalously large with c 12 about eight times larger than the harmonic term b 02 . For perspective, this is comparable to the improper FE coupling between the polar mode and the zone-boundary mode in YMnO 3 [38]. The effective quadratic coefficient (B eff 02 ), including the anharmonic terms from Eq. (1), is The strong anharmonic coupling renormalizes B eff 02 to become more negative. Similar renormalization of the optimal polarization and the Curie-Weiss temperature in conventional oxide FEs was discussed due to electrostriction [45]. But an anomalously large contribution of the fully symmetric mode leading to this renormalization as we find in CIPSe is new. TABLE V. The value of the coefficients (in meV) in Eq. (1) from fitting the energy surface shown in Fig. 4(a).

F. The microscopic mechanism of the coupling
Prior works discussed the possibility that a filled d 10 outershell electronic configuration is second-order Jahn-Teller active, leading to the off-centering of both the Cu and In atoms in this class of materials [9,[46][47][48]. So to elucidate the physical origin of the large nonlinear coupling between the polar mode and the fully symmetric mode, we analyzed the various contributions to the total energy. Specifically, we focus on the Kohn-Sham eigenvalues and the Ewald energyi.e., the electrostatic ion-ion interaction terms for the FiE phase relative to the PE phase, as shown in Figs. 4(b) and 4(c), respectively. The former corresponds to the electronic component of the second-order Jahn-Teller effect, and should lead to the lowering of the Kohn-Sham energy. The latter can be compared to the effect of distortion without an electronic redistribution which should lead to an increase in energy. Together they model a second-order Jahn-Teller effect [45]. This is shown as a function of Q + 1 and Q − 2 in Fig. 4. We find that the polar mode lowers the total energy by reducing the Kohn-Sham electronic contribution to the energy [see Fig. 4(b)]. But the condensation of the polar mode results in a large Ewald energy penalty [see Fig. 4(c)]. Physically, as the Cu ions move closer to the Se ions [see Fig. 1(e)], thereby lowering the Kohn-Sham energy, the Ewald part of the total energy goes up due to an increased ion-ion repulsion. However, this shorter Cu-Se t bond can be stabilized by an expansion (contraction) of the trigonal Se t (Se b ) plane. As shown in Table III Fig. 1(f) for the − 2 A 1 mode, in addition to this mode the + 1 A 1 mode [see Fig. 1(i)] is needed to stabilize the displaced Cu atoms. This naturally leading to a coupling between the two modes. The expansion and contraction of the Se sublayers we find is consistent with prior x-ray diffraction studies in this class of compounds [9]. Further, Fig. 4(c) shows that this coupling between the polar mode and the fully symmetric mode overcomes this energy penalty. So, the stabilization of the FiE phase in CIPSe is mediated by the coupling between the metal distortions and the chalcogen atoms that also contributes to the fully symmetric Raman active mode. Such a competition between the electronic degrees of freedom and the repulsive forces between the ions is in line with a second-order Jahn-Teller effect [9,[45][46][47][48].
Our analysis gives a clear microscopic mechanism of the unusual coupling between the polar mode and the fully symmetric mode that leads to the stabilization of the FiE phase. We find that this anharmonic stabilization persists even in the monolayer limit (see Fig. S5 in the Supplemental Material [34]).

G. Strain-control of FiE
We also explored the effect of other distortions that transform as the + 1 irrep of the space group P31c. Previous works on oxides showed that polar-distortions naturally couple to strain [35,49]. So, we discuss the case of an in-plane biaxial strain ( x 2 +y 2 ), which uniformly expands the 2D-area for positive values. The effect of this strain is similar to that of the FIG. 5. The total energy as a function of the fractional amplitude of the distortion from the PE phase to the FiE phase, for different value of biaxial strain ( x 2 +y 2 ). Results for tensile strains of 2% (blue), 1% (black), the unstrained case (red), and compressive strain of 1% (green), and 2% (magenta) are reported. The energy difference can be enhanced or suppressed by the application of an in-plane strain.
displacive + 1 mode [see Fig. 1(i)]. While only the Se atoms participate in the displacive mode, the strain mode includes all the atoms. Naturally, we expect this component of the strain to have a larger effect on the anharmonic coupling. We also check the contribution from other fully symmetric strain mode and found that the dominant effect comes from this biaxial strain mode. Figure 5 shows the total energy profile for the FiE phase relative to the PE phase for various value of x 2 +y 2 . For the unstrained case (red solid line) we recover the FiE phase as the low-energy state. We find that for a compressive strain of 1% (green line) the FiE phase, although a local minimum, has higher energy (4 meV/f.u.) than the PE phase with a small energy barrier (12 meV/f.u.) between the states. For intermediate values of compressive strain where the FiE remains the ground state but where the energy barrier is comparable to the quantum fluctuations a quantum PE state could be stabilized by strain where quantum fluctuations suppress the ordering of dipole moments [50,51]. Incidentally, a quantum paraelectric phase was proposed in this family of materials driven by chemical substitution [52].
For a compressive strain of 2% (magenta line) the polar instability is fully suppressed and the PE phase becomes the lowest energy structure. On the other hand, for a tensile strain (1% shown in black lines and 2% shown in blue line) the well-depth around the FiE phase becomes deeper. Note that the FiE phase has a x 2 +y 2 strain of ∼ 1% relative to the PE phase (see Tables I and II), and largely accounts for the additional 37 meV gain from relaxing the lattice constants in the FiE phase. As biaxial strain can be easily controlled experimentally instead of directly tuning the displacive + 1 mode, mechanical strain can be an effective tool to both suppress and enhance the FiE phase, providing a way to tune the energy barrier corresponding to a single-step switching of the polarization.

IV. CONCLUSION AND OUTLOOK
While our analysis is done at 0 K in the absence of an external electric field and domains, it gives some preliminary insights into the factors that govern the switching of polarization under an applied electric field. Within a singlestep switching process we assume that the energy barrier for switching is determined by the energy difference between the FiE phase and the PE phase. Our analysis, shows that this barrier is determined primarily by a large nonlinear coupling between the polar mode and the fully symmetric Ramanactive mode. This opens up possibilities to both enhance and suppress the energy barrier for polarization switching by modifying the fully symmetric mode, instead of the polar mode. Since the fully symmetric mode is Raman active, in addition to mechanical strain, it should also couple to ultra-fast optical probes allowing for the optical control of the switching barrier. [53][54][55]. The application of such an optical field should be faster than applying a mechanical strain.
When the barrier energy is reduced sufficiently, a faster electrical switching of the polar state can be achieved within the transition state theory, assuming a one-step switching process. On removal of the two fields the higher barrier around the FiE phase, for example in the unstrained case, can be recovered. This will lead to a long-lifetime of the polar phase, which can be even enhanced for instance by engineering a tensile strain as shown in Fig. 5. This understanding can provide a first step to addressing the voltage-time dilemma [56].
We note that this approach for overcoming the voltagetime dilemma is not unique to CIPSe. Other materials where the structural stabilization of the ordered phase involves significant contributions from modes other than the modes describing the memory-order parameter has the potential to overcome this inherent bottleneck. In the case of CIPSe, while the metal polar-distortions created the memory orderparameter the stabilization of the polar phase was due to the anharmonic coupling between the polar mode and the fully symmetric mode. Improper FEs, pseudoproper FEs, hybrid improper FEs and triggered phases are all potential candidates for this as the polarization itself is not the primary order parameter and the stabilization of the polar phase comes from the coupling between polarization and other structural distortions. We note that a strain-induced suppression of the barrier was recently reported in hybrid improper FE [57].
Finally, we comment on the nature of the FiE phase. Within a displacive picture CIPSe is a proper FiE as the polar mode is the primary order parameter and fully determines the symmetries of the FiE phase. However, unlike other proper FEs where the energy lowering term is P 2 [35,[58][59][60], the driving term that stabilizes the polar phase is dominated by the anharmonic coupling between the polar mode and the fully symmetric Raman active mode. This distinguishes the ferrielectricity in CIPSe and other layered thiophosphates from other known FEs.
While our discussion is centered around a on-site model for the total energy, our results should also be important in a multi-step switching process. Future works need to include the effect of zone-boundary instabilities also to understand the order-disorder nature of the phase transition and the intermediate incommensurate phase reported at the phase boundary. Similarly, the impact of our findings in the context of the unconventional electromechanical response found in the domain walls in this family of materials, and in the actual path of polarization switching under an applied field also needs to be explored.
In summary, we discuss the microscopic origin of polarization and the structural-stabilization of ferrielectricity in CIPSe. The distortions of the Cu and In atoms create most of the polarization in FiE CIPSe. On the other hand, the anharmonic coupling between the polar mode and the fully symmetric distortions involving the Se atoms stabilizes the FiE phase. The energy barrier for the single-step switching of the polarization in CIPSe can be controlled by changing the fully symmetric mode for instance, by the application of a uniform in-plane strain. This can help overcome the voltagetime dilemma in microelectronic device applications. Further, CIPSe could potentially be driven into a quantum paraelectric phase by the application of a nominal strain.
The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan [61]. In addition, simulation inputs and outputs for the DFT calculations performed in this work are available via the Materials Data Facility [62].