Quantum paraelastic two-dimensional materials

TL;DR

This study investigates the elastic energy landscape and phase transition of 2D tin oxide monolayers, revealing a quantum paraelastic phase that can be tuned via charge doping, with implications for 2D quantum materials.

Contribution

It provides the first detailed analysis of quantum paraelasticity in 2D SnO monolayers using ab-initio MD and phonon mode analysis.

Findings

Transition temperature of ~8.5 K identified.

Soft phonon modes linked to structural transformation.

Quantum paraelastic phase potentially tunable by charge doping.

Abstract

We study the elastic energy landscape of two-dimensional tin oxide (SnO) monolayers and demonstrate a transition temperature of $T_c=8.5\pm 1.8$ K using ab-initio molecular dynamics (MD), that is close to the value of the elastic energy barrier $J$ derived from $T=0$ K density functional theory calculations. The power spectra of the velocity autocorrelation throughout the MD evolution permits identifying soft phonon modes likely responsible for the structural transformation. The mean atomic displacements obtained from a Bose-Einstein occupation of the phonon modes suggest the existence of a quantum paraelastic phase that could be tuned with charge doping: SnO monolayers could be 2D quantum paraelastic materials with a charge-tunable quantum phase transition.

Tables

Table 1. Table 1: Orbital character of added (negative sign) or removed (positive sign) electronic charge per unit cell, ratio among lattice parameters a 1 / a 2 subscript 𝑎 1 subscript 𝑎 2 a_{1}/a_{2} , and J 𝐽 J ( K 𝐾 K /u.c.).

$q$ ($\left|e\right|$)	$\Delta s$	$\Delta p_x$	$\Delta p_y$	$\Delta p_z$	$a_1/a_2$	$J$
$-0.05$	$-$0.019	$-$0.001	$-$0.001	$-$0.029	1.06	5
$0.05$	0.013	0.000	0.003	0.034	1.12	50
$0.15$	0.041	0.000	0.010	0.099	1.13	54
$0.25$	0.068	0.001	0.013	0.168	1.10	23
$0.35$	0.095	0.002	0.010	0.242	1.05	3

Equations

$$X(a_1,a_2)=\frac{a_1-a_2}{\sqrt{2}};Y(a_1,a_2)=\frac{a_1+a_2-2a_C}{\sqrt{2}},$$

$$U(X,Y)=a{X}^4-b{X}^2+c{Y}^2+\frac{b^2}{4a},$$

$$⟨f(X,Y)⟩\equiv \frac{{\int }_{X_{min}(T)}^{X_{Max}(T)}{\int }_{Y_{-}(X,T)}^{Y_{+}(X,T)}{e}^{-U/k_{B}T}f(X,Y)dXdY}{{\int }_{X_{min}(T)}^{X_{Max}(T)}{\int }_{Y_{-}(X,T)}^{Y_{+}(X,T)}{e}^{-U/k_{B}T}dXdY}.$$

Full text

Structural transformations in two-dimensional (2D) materials were first observed in transition metal dichalcogenides Zhang et al. (1996); Kim et al. (1997); Duerloo et al. (2014) exposed to large gate and bias voltages Zhang et al. (1996); Kim et al. (1997); Wang et al. (2017a). These phase transitions induce two-dimensional ferroic and multiferroic behavior Mehboudi et al. (2016a); Seixas et al. (2016); Chang et al. (2016). Three conditions underpin ferroic behavior: (i) the presence of structural degeneracies leading to a high-to-low-symmetry structural configuration as temperature ($T$) is decreased, (ii) a low-energy barrier $J$ along the path joining the degenerate ground states, and (iii) the existence of macroscopic monodomains Naumis et al. (2017). The existence of macroscopic domains may not be granted when $J$ is too small.

With a predicted multiferroic behavior Seixas et al. (2016), SnO monolayers offer an ideal testbed for this observation: Similar to 2D iron sulfide (FeS), SnO monolayers have a rectangular unit cell in their low-temperature ground state. Neglecting quantum fluctuations, a transition temperature $T_{c}$ can be determined via molecular dynamics, beyond which these compounds turn into a traditional litharge structure; i.e., a square unit cell that otherwise maintains its atomistic coordination. But while the material behavior of FeS is driven by at least three competing energy scales of similar magnitude (the superconducting gap, spin-splitting/magnetic anisotropy barrier, and the elastic barrier leading to the rectangular-to-square unit cell transformation), the structural transformation of a pristine charge-neutral SnO monolayer can be characterized without concern for these competing energy scales.

In this Letter we study the ferro- to para-elastic transition in SnO monolayers, starting from the parametrization of their elastic energy landscape, a discussion of their structural degeneracies, and of their elastic energy barrier $J$. Using ab initio calculations of the unit cell, ab initio molecular dynamics (MD) calculations on the NPT ensemble, and model calculations, we demonstrate that the structural transition temperature $T_{c}$ is very similar to $J$. The power spectra of the MD evolution permits identifying soft phonon modes likely responsible for the observed structural transformation. The Bose-Einstein occupations of the phonon modes lead to increased atomistic fluctuations that may stabilize a quantum paraelastic phase at low temperature that is tunable via electrostatic gating.

Figure 1(a), obtained with the VASP code Kresse and Furthmüller (1996); Kresse and Joubert (1999) that implements the plane-wave pseudopotential approach within DFT with PBE exchange-correlation pseudopotentials Perdew et al. (1996); Kresse and Joubert (1999), shows the energy landscape for a SnO monolayer. Optimizations were performed with forces smaller than $10^{-3}$ eV/Å and using up to 90$\times 90\times 1$ k-points. There, structure $A$ has a rectangular unit cell with lattice constants $a_{10}=4.01$ Å and $a_{20}=3.68$ Å, as highlighted by a star. O atoms have coordinates $\mathbf{b}_{1}$=(0,0,0) and $\mathbf{b}_{2}$=($a_{10}/2$,$a_{20}/2$,$-\delta_{z}$), while Sn atoms are located at $\mathbf{b}_{3}$=($a_{10}/2$,[math],$\Delta_{z}$) and $\mathbf{b}_{4}$=([math],$a_{20}/2$,$-\delta z-\Delta z$), with $\delta z=0.18$ Å, and $\Delta_{z}=1.08$ Å, respectively. This way, $d_{10}=2.28$ Å and $d_{20}=2.23$ Å. Structure $B$ is obtained by swapping $x-$ and $y-$coordinates, making structures $A$ and $B$ degenerate Mehboudi et al. (2016a); Barraza-Lopez et al. (2018); Naumis et al. (2017) and the SnO monolayer ferroelastic Seixas et al. (2016).

Point $C$ in Fig. 1(a) has lattice constants $a_{1}=a_{2}=a_{C}=(a_{10}+a_{20})/2$, and it sits midway through the path that joins points $A$ and $B$ (black dashed line). In structure $C$, O atoms are located at $\mathbf{b}_{1}$=(0,0,0) and $\mathbf{b}_{2}$=($a_{C}/2$,$a_{C}/2$,0), and Sn atoms are at $\mathbf{b}_{3}$=($a_{C}/2$,[math],$\Delta_{C}$) and $\mathbf{b}_{4}$=([math],$a_{20}/2$,$-\Delta_{C}$) with $\Delta_{C}=1.17$ Å, for an interatomic distance $d_{C}$=2.25 Å. Figure 1(b) shows side and top views of rectangular structures $A$ and $B$, and the (litharge) structure $C$.

Defining length variables $X$ and $Y$ as linear combinations of lattice parameters $a_{1}$ and $a_{2}$:

[TABLE]

one creates the energy cuts shown in Fig. 1(c) on the energy landscape along the dashed black line at the $(X,0)-$direction, and the dashed red lines at the $(0,Y)-$ and $(\pm X_{0},Y)-$directions indicated in Fig. 1(a).

The dependence of the landscape on $Y$ was almost identical along the red diagonal lines in Fig. 1(a) cutting through points $A$, $B$ and $C$: it is $J+c(X=0)Y^{2}$ around point $C$, and $c(X=X_{0})Y^{2}$ around point $A$ (or $B$) with $c(X=0)$=23,042 KÅ*-2*/u.c. and $c(X=X_{0})$=22,364 KÅ*-2*/u.c., for a negligible dependency of the coefficient $c$ on $X$; $c$ is approximated by $(c(0)+c(X_{0}))/2$, to simplify the analytic description of the 2D energy landscape:

[TABLE]

with $a=4,252$ K/Å4, $b=387$ K/Å2, and $c=22,703$ K/Å2 per unit cell.

Setting $U(X_{0},0)=0$, $X_{0}=\pm\sqrt{\frac{b}{2a}}=0.21$ Å, and inverting Eq. 1 at the analytic value of $X_{0}=0.21$ Å, $Y=0$ Å, one gets $a_{10}=3.99$ and $a_{20}=3.70$ Å, which agrees with DFT results in Fig. 1 reasonably well. $J\equiv{}U(0,0)-U(X_{0},0)=\frac{b^{2}}{4a}=8.8$ K/u.c.

When expressed in units of K/u.c. $J$ turns out to be almost equal to $T_{c}$, as shown via ab initio MD performed at ambient pressure and at 12 predefined temperatures. These calculations, performed using the NVT ensemble (i.e., with a constant number of electrons and target pressure and temperature), require having two walls parallel to (not in contact with) the 2D material fixed throughout the MD evolution.

The Qbox MD code Gygi (2008); Schlipf and Gygi (2015) implements the plane-wave pseudopotential approach within DFT with highly-tuned PBE exchange-correlation pseudopotentials that have tested against multiple other plane wave codes for numerical consistency Schlipf and Gygi (2015). It permits fixing the two walls not in direct contact with the SnO monolayer. MD data for charge-neutral SnO monolayers is presented in Figs. 2 and 3 and discussed next.

$T_{c}$ is signaled by a sudden change of lattice parameters $a_{1}$ and $a_{2}$ onto $a_{C}$, highlighting the importance of performing MD on an NPT ensemble to capture the transition temperature ($T$) Mehboudi et al. (2016a, b); Barraza-Lopez et al. (2018). MD calculations were performed on a $4\times 4$ SnO supercell originally set at configuration $A$, which assumes the existence of a macroscopic domain with such atomistic configuration (condition (iii) in Ref. Naumis et al. (2017)). Using a 1.54 fs time resolution, we retrieved 16,000 frames, and used the last $2^{13}$ (=8192) frames to obtain ensemble averages and power densities in the frequency domain for any given temperature ($T$).

Figures 2(a), 2(d) and 2(e) show the thermal evolution of lattice parameters ($\langle a_{1}\rangle$, $\langle a_{2}\rangle$), interatomic distances ($\langle d_{1}\rangle$, $\langle d_{2}\rangle$), and angles ($\langle\alpha_{1}\rangle$, $\langle\alpha_{2}\rangle$) (c.f., Fig. 1(b)), where $\langle\rangle$ indicates ensemble averages. There is a sudden collapse of these order parameters at $T_{c}=8.5\pm 1.8$ K, which ratifies the ansatz $T_{c}=\beta J$, with $\beta\simeq 1$ suggested early on. The thickness in Fig. 2(f) is the height difference among the lowest and the highest atom in the supercell, and it captures out-of-plane undulations Barraza-Lopez et al. (2018).

Together with earlier NPT MD results Mehboudi et al. (2016a, b), the value of $T_{c}$ obtained here points to the central role of the energy scale $J$ in understanding multiferroicity in 2D materials. The entirety of our MD data points to a relation $T_{c}\simeq J$.

It is illustrative to describe the structural transition right from the energy landscape in Fig. 1. Consider a particle there, and assign a mean kinetic energy equal to $k_{B}T/2$ to each of the two degrees of freedom ($k_{B}$ is Boltzmann’s constant). In this classical construct, the particle has zero potential energy and a mean kinetic energy of $k_{B}T$ at point ($X_{0},0$). The particle reaches the largest values of $X$ and $Y$ when its mean kinetic energy is zero, and the condition $k_{B}T=U(X,Y)$ determines the contours in the $(X,Y)$ plane the particle is confined to, as a function of $T$. Depending on whether $k_{B}T$ is smaller or larger than $U(0,0)=J$, the particle is either confined to a single well (monodomain $A$), or it oscillates among the two wells (i.e., the material has transitioned).

Equation $k_{B}T=U(X,Y)$ has four roots: $X_{P\pm}(T)=\sqrt{(b\pm\sqrt{4ak_{B}T})/2a}$, and $X_{N\pm}(T)=-X_{P\pm}(T)$, where $N$ ($P$) stands for negative (positive). Being a classical construct, the elastic energy profile forbids direct tunneling among the two wells, so one is constrained to ($X_{min}(T)=X_{P-}(T)\leq X\leq X_{P+}(T)=X_{Max}(T)$) when $k_{B}T\leq J$ at monodomain $A$. Both wells are accessible when $T>J$, and $X_{min}(T)=X_{N+}(T)\leq X\leq X_{P+}(T)=X_{Max}(T)$. $X_{min}(T)$ takes on two different values, depending on whether $T\leq J$ or $T>J$.

Contours are given by: $Y_{\pm}(X,T)=\pm\sqrt{(k_{B}T-aX^{4}+bX^{2}-J)/c}$, and ensemble averages within the model for any function $f(X,Y)$ are obtained from Kittel (2004):

[TABLE]

The solid lines in Fig. 2(a) are $\langle a_{1}\rangle$ and $\langle a_{2}\rangle$; shaded areas are standard deviations. These values were obtained from Figs. 2(b) and 2(c) through an inversion of Eqn. 1. The standard deviations on $a_{1}$ and $a_{2}$ are equal to 0.048 Å, right below $T_{c}$. The agreement between the model and MD –including fluctuations– is remarkable.

Despite of the small supercell size (inset in Fig. 2(f)), the power spectral density (PSD) –the Fourier transform of the MD data onto reciprocal space and frequency Koukaras et al. (2015)– shown in Figs. 3(a,b) displays sharp peaks that uncover the phonon frequencies of the SnO monolayer, and include anharmonic contributions by construction. The PSD is resolved over chemical species (with slow modes mostly due to Sn, and overall faster modes due to O, which was vertically shifted in these plots) and orthogonal displacements $L$, $T$, and $Z$ Kittel (2004). The supercell provides a coarse sampling of reciprocal space that includes $k-$points $i\mathbf{b}_{1}/4+j\mathbf{b}_{2}/4$, where $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ are reciprocal lattice vectors and $i,j$ range in between $-2$ and $2$.

The rectangular unit cell imposes an asymmetry on the PSD shown in Fig. 3(a) at $T=4.0\pm 0.8$ K $<T_{c}$, that is best brought about by placing the spectral density at point $\mathbf{b}_{2}/4$ over that acquired from point $\mathbf{b}_{1}/4$, and that of $\mathbf{b}_{2}/2$ over $\mathbf{b}_{1}/2$, and aided by vertical dashed lines originating from plots $\mathbf{b}_{1}/4$ and $\mathbf{b}_{1}/2$. There, the Z-mode in orange at 14 THz at $\mathbf{b}_{1}/4$ appears at a lower frequency at $\mathbf{b}_{2}/4$. In a similar manner, the O $T$ mode seen in gray at about 8.5 THz at $\mathbf{b}_{1}/4$ is downshifted onto the $L$ mode in light blue at $\mathbf{b}_{2}/4$. A similar downshift is seen at the maximum frequency mode at $\mathbf{b}_{2}/2$ with respect to its value at $\mathbf{b}_{1}/2$. There is a significant spectral density at frequencies $\gtrsim$ 0 THz and at least two closely spaced peaks at the $\Gamma-$point in Fig. 3(a). As showcased at $T=16.7\pm 3.2$ K in Fig. 3(b), the PSD shows aligned peaks for $T>T_{c}$, thus reflecting an enhanced symmetry.

The power spectrum is rearranged into a phonon dispersion in Fig. 3(c) for $T<T_{c}$ (left, mode resolved) and $T>T_{c}$ (right), while the overlaid continuous trend is a phonon dispersion at 0 K, without anharmonic contributions, that was obtained from the Hessian matrix at $T=0$. Direct comparison of the continuous trend lines and the PSD data shows at least one missing phonon line above 10 THz at the $\Gamma-$ point on the left subplot (highlighted by a vertical arrow), which appears to have gone soft to accumulate near 1 THz frequency (as suggested by a dashed overlaid curve). The possibility of the mode at 10 THz being absent due to the small supercell employed still does not answer the existence of the large spectral density near 1 THz at the $\Gamma-$point seen in Figs. 3(a) and 3(c), that does not match any phonon line suggested by the Hessian. Diagonal arrows highlight the reduced symmetry when $T<T_{c}$.

Ferroic behavior necessitates the feasibility of creating macroscopic monodomains Naumis et al. (2017); i.e., that the two wells on the elastic energy landscape are separated by a classical wall that forbids quantum tunneling. This classical limit is not guaranteed here due to the small value of $J$ and the presence of a light element (oxygen) prone to significant zero-point, quantum fluctuations.

Using the phonon spectrum and eigenvectors obtained from the Hessian, one can contrast the phonon energy and the mean displacements obtained using occupations consistent with Bosonic quantum statistics, with the occupations resulting from a classical occupation that is linear in temperature, and which underlies all the results presented up to this point. The drastic enhancement in occupation when Bosonic statistics are considered may drive a quantum paraelastic phase reminiscent of the quantum paraelectric phase in SrTiO3 Cowley (1964); Schneider et al. (1976); Müller and Burkard (1979).

The energy due to phonons with a Bosonic occupation (BE) shown by a red solid line in Fig. 4(a) displays the usual uptick at zero temperature; the classical occupation is linear in temperature, and it agrees with the quantum occupation only at large $T$. The Bosonic occupation in turn produces the uptick on the rms atomic displacements at low $T$ Masri and Dobrzynski (1971) depicted in Figs. 4(b-d), when compared with results obtained using occupations linear in $T$. The solid lines coalesce onto the dashed lines at large $T$. The solid vertical line in Figs. 4(b-d) is $T_{c}$, while the dashed area indicates its standard deviation.

Figure 4(d) indicates that out-of-plane displacements $\langle\Delta z^{2}\rangle$ are not too dissimilar from both statistics. But in-plane components $\langle\Delta x^{2}\rangle$ and $\langle\Delta y^{2}\rangle$ obtained from the BE occupation show large excursions at $T=0$ already, only attainable with a classical occupation at higher $T$. Figure 4(e) suggests that $J$ can be increased by a slight hole doping Seixas et al. (2016) and further decreased by either electron doping or a substantial hole doping: charge doping is a handle to get the SnO monolayer in and out of the quantum paraelastic phase.

$J$ depends on whether the charge added is symmetric or asymmetric on the $p_{x}$ and $p_{y}$ orbitals; c.f. Table 1. The apparent correlation among $a_{1}/a_{2}$ and $J$ in Table 1 is reminiscent of that given in Ref. Mehboudi et al. (2016a). In turn, the location in reciprocal space where electron charge is added or removed seen in Fig. 5 determines the orbital character.

In the right side of Figs. 4(b,c), the in-plane rms displacement with Boson occupations may lead to a lower $T_{c}$, as suggested by the horizontal orange arrows in Figs. 4(b-c). The renormalization of $T_{c}$ is increased as $J$ is brought to smaller values, as displayed by larger horizontal green lines in Figs. 4(b-c). For values of $J$ smaller than 21 K/u.c., fluctuations of the order of 0.048 Å–that drove the classical transition in Fig. 2(a)– will occur at 0 K already, to set the quantum paraelastic phase.

There are experimental and engineering challenges related to sample quality and exfoliation down to monolayers that must be met to observe these phenomena. The SnO monolayer limit is yet to be achieved, though encouraging results demonstrating exfoliation down to 4 monolayers and 2 monolayers have been recently reported Zhang et al. (2017); Singh et al. (2017). SnO samples that show ambipolar doping Nomura et al. ; Varley et al. (2013) will be ideal to test the present predictions, and additional control of electron doping, necessary to further explore the paraelastic phase, may be achieved via electrostatic doping with gates (e.g., Ref. Wang et al. (2017b)).

In conclusion, we have identified charge-neutral SnO monolayers as incipient quantum paraelastic materials; i.e., materials in which a ferroic structural configuration may be forbidden due to quantum fluctuations. It has been shown that the elastic energy barrier, expressed in K/u.c., is a reasonable estimate for the transition temperature of charge-neutral SnO monolayers. The power spectrum obtained directly from the MD evolution show missing modes, that go to lower energies prior to $T_{c}$ and herald the structural transition. The magnitude of $T_{c}$ may be tuned via charge doping, for quantum statistics to dictate the ferroic behavior of these materials.

T.B. was funded by the NSF (Grant No. DMR-1610126), A.P.S. by FONDECYT, project No 1171600 (Chile), P.D. by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. DOE under contract DE-AC02-06CH11357, and S.B.L. by the U.S. DOE, Office of Basic Energy Sciences (Early Career Award DE-SC0016139). Part of this work was performed at the Center for Nanoscale Materials, a U.S. Department of Energy Office of Science User Facility, and supported by the U.S. DOE, Office of Science, under Contract No. DE-AC02-06CH11357. Calculations were performed on Cori at NERSC, a U.S. DOE Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. Conversations with P. Kumar, B. Fregoso and G. Naumis are gratefully acknowledged.