Electronic, dielectric and optical properties of two dimensional and bulk ice: a multi-scale simulation study

We have calculated the electronic band structure of 2D-ices using density functional theory (DFT) as implemented in the Quantum-ESPRESSO (QE) package ¹⁷. We have used ultrasoft pseudopotentials to treat the interaction between the ion cores and valence electrons and applied the generalized gradient approximation (GGA) for the exchange-correlation interactions. We also have studied the effect of nonlocal correlations using the van der Waals density functional (vdW-DF) of Dion et al ¹⁸. Accuracy of forces on each atom has been considered about 0.1 mRy/bohr for the variable-cell optimization, relaxing the cell parameters and atomic positions. in order to accurate calculation of structural and electronic properties, the kinetic energy cut-offs of $150~{}Ry$ and $1500~{}Ry$ were found to be sufficient for the wavefunctions and the charge densities, respectively, where the k-point grid for 2D- (bulk) ices was set to $6\times 6\times 1$ ($6\times 6\times 6$). The k-point grid for 2D- (bulk) ices was chose $24\times 24\times 1$ ($24\times 24\times 24$) for non-self-consistent calculation in partial density of state (PDOS) analysis. The smearing parameter of 0.01 Ry has been used for PDOS analysis. In addition, total energy convergence threshold was set to $10^{-12}eV$. The Coulomb cutoff technique ^{19; 20} was used to reduce interactions between periodic images and cost of the ab-initio calculations for 2D-ice structures.

The dielectric properties of the 2D-ices were determined using the DFPT approach. In DFPT, the dielectric constant tensor is defined as a linear response to the perturbative electric field ^{21; 22} and the ionic displacement are considered as a perturbation to the equilibrium system. The response of the electronic charge density to the perturbative electric field in the linear response regime, employed to determine the electronic contribution of the dielectric tensor, i.e. the high-frequency dielectric constant ($\varepsilon_{el}$). Subsequently, the static dielectric constant ($\varepsilon_{0}^{\mu\nu}$) is the summation of the electronic and the ionic parts of the system to the applied electric field:

where

Here $S_{m,μν}$ is the mode oscillator strength tensor, defined in terms of the Born effective charges $Z^{*}$ , the atomic masses $M_{i}$, and the normalized eigenvectors $u_{i,m\mu}$ of the i^th ion along a given direction $\mu$ for a particular mode $m$. Also $\omega_{m}$ is the phonon mode frequency and $\Omega$ is the unitcell volume. Thus, in order to compute $\varepsilon_{0}$, the knowledge of all the phonon frequencies at the zone center of Brillouin zone is needed. The latter requires the solution of the dynamical matrix at the zone center.

By means of equilibrium molecular dynamics simulations employing the large scale atomic/molecular massively parallel simulator LAMMPS ²³, we computed the molecular (dipolar) dielectric constant of 2D-ice at 80 K. Our simulated systems contain 2D-ices, confined between two walls, separated by a distance ² $h=2\,t$ (these walls are used to produce the confining potential). The number of water molecules are 400, 100 and 128 for f-SQ, R-SQ, and HEX, respectively²⁴. The TIP4P model was employed to describe the water molecules ²⁵. The NVT ensemble (Nose-Hoover thermostat) is used to keep the temperature at 80 K for 2D- and bulk ices excluding ice Ih (200 K). We modified TIP4P model to incorporate the lattice structure of the optimized structures by using DFT: we changed bond lengths and bond angles of water molecules relevant to the DFT outputs ². Also, the charges of O and H atoms were modified such that the large dipole moment ($\sim$ 3D) for single water molecule is reached ⁹. In order to verify this modification, we calculated the variation of dipolar dielectric constant of bulk water with time using the modified TIP4P model, see Fig. 3(a). It is seen that the dielectric constant of bulk water lies in acceptable range, i.e. it is between dielectric constant of bulk water using TIP4P and TIP4P/2005 models.

Periodic boundary conditions are applied along $x$, $y$ directions and the confinement was along the $z$-direction. The particle-particle particle-mesh method was used to compute the long-range Coulomb interaction with a relative accuracy of 10^-4. Water bond lengths and bond angles were fixed by the SHAKE algorithm ²⁶. In all MDS a time step of 1 fs was chosen. Following the system relaxation (1 ns), the thermodynamical sampling was done up to 8 ns to ensure the smoothness of the converged dielectric constant. The temperature in MDS for bulk ices (water) set to be 80 K (300 K); because all of the studied bulk ices are stable at this temperature.

A microscopic picture of dielectric properties of 2D-ices could be found by calculating the fluctuations of the total polarization of a system, $\vec{M}$ at finite temperature. By calculating different components of the total dipole moment i.e. $M_{x}$, $M_{y}$, and $M_{z}$ after equilibration, one obtains different components of the molecular dielectric constant tensor as ²⁷

where $\varepsilon_{\infty}$ is the optical dielectric constant and taken to be 1. Also, $\sigma_{\mu\nu}^{2}$ $=$ $\langle M_{\mu}M_{\nu}\rangle-\langle M_{\mu}\rangle\langle M_{\nu}\rangle$ while $\mu,\nu=x,y,z$ and $V$ is the volume of the system. Here, the time averaging was taken for more than 5 ns when in-plane dielectric constant is converged. Note that Eq. (4) can only be used for a homogeneous systems ¹¹.

The optical dielectric function $\varepsilon_{el}(\omega)$ was calculated using norm-conserving pseudopotentials, in the energy range of 0 to 30 eV. In order to increase accuracy for the dielectric functrion calculations, we used the k-point grid for 2D- (bulk) ices as $12\times 12\times 1$ ($12\times 12\times 12$). The optical dielectric constant was calculated (the frequency dependent of the electronic dielectric constant or equivalently the electronic part of the dielectric function) within the framework of the random-phase approximation (RPA) based on DFT ground-state calculations. The mentioned dielectric function consists of frequency dependent real ($\varepsilon_{el}^{r}(\omega)$) and imaginary part ($\varepsilon_{el}^{i}(\omega)$). It is represented as:

The imaginary part of dielectric function ($\varepsilon^{i}(\omega)$) can be calculated using Kubo–Greenwood formalism ²⁸. Once we know the imaginary part, the real part can be obtained using the Kramers–Kronig relations (for more details see Ref. [29]). Note that the optical dielectric constant extracted from RPA (at zero frequency) is compatible to that obtained from DFPT.

The electronic band structure and partial density of states (PDOS) of 2D- and bulk ices are shown in Fig. 2. We found the direct band gap of energy about 5.49, 5.55 and 5.12 $eV$ for f-SQ, b-SQ and HEX 2D-ices and 5.37 and 4.99 $eV$ for ice VIII and ice Ic, respectively. The electronic band structure for b-RH 2D-ice, ice XI and ice Ih confirm the indirect band gap. The corresponding energy gap for b-RH, ice XI and ice Ih are 5.55, 5.06 and 5.08 $eV$, respectively. The PDOS are shown in the right side of the panels of Fig. 2. It is seen that in contrast to the conduction bands, the O (H) atoms contribution in the valance band is larger (negligible). Notice that the difference between the energy gap of 2D-ices (and 2D-h-BN) and bulk water (bulk h-BN) is small, i.e. 2$\%$ (4$\%)$, as compared to the difference between bulk TMDs and 2D-TMDs, see Fig. 4. This is due to the weak hydrogen bond between water molecules in bulk and 2D-ices as well as insulating feature in all phases of water and h-BN. The TMDs have large reductions of energy gap while transiting from 3D to 2D. The obtained energy gap, the electronic band structure and PDOS for ice VIII, ice XI,and ice Ih are in agreement with the results of Refs. [32-35].

Independent of water phase (including ice), the nearest neighbor distance between oxygen atoms of water molecules are almost the same, i.e. 2.8-3Å  i.e. all the ice phases formed under very high pressures satisfy the so-called Bernal-Fowler ice rules where each water molecule has four hydrogen-bonded neighbors with a quasi-tetrahedral configuration with two short O-H distances (the donated protons) and two long ones (the accepted protons). Transiting from bulk into 2D phase of ice, the crystalline structure with larger density can be formed where the nearest neighbor distance are more or less the same but, the preferential tetrahedral bonding geometry is different. In Figs. 3(b,c) we depict the in-plane radial distribution function for O-O distance in f-SQ and 3D-radial distribution function of three bulk ice (Ice XI, ice Ic and ice Ih). As can be seen, as expected the nearest neighbor O-O distance (first peaks) for all structures are equal. However, the second nearest neighbors (shown by arrows in Fig. 3(b)) are different. The latter causes a larger density of f-SQ ⁶ (1.36gcm^-3) as compared to bulk crystalline ices i.e.   0.92 gcm^-3 ice Ih and   0.92 gcm^-3 ice XI.
On the other hand, the most well-known electrical conductivity mechanism for all phases of water/ice is the Grotthuss mechanism (GM) also known as proton jumping. In GM an excess proton or proton defect diffuses through the hydrogen bond network of water molecules and a covalent bond between neighboring molecules are formed and broken continuously. Because of the same nearest neighbor distance in bulk ice (bulk water) and those of 2D-ice, the GM mechanism should be valid for 2D-ices. Subsequently only an infinitesimal increase in the band gap 2D-ices as compared to the bulk ices is seen. This was confirmed by our ab-initio simulations.
Moreover, h-BN is a wide band gap semiconductor with high thermal and chemical stability. In both bulk and monolayer h-BN, N atoms and B atoms are hybridized with sp² at the interlayer forming a honeycomb structure. In bulk h-BN there is weak interactions between each layer of h-BN, such as electrostatic interactions and Van der Waals interactions, which causes the electrical band gap of both bulk h-BN and 2D-h-BN become more or less equal ³⁶. Transition from bulk into 2D-TMDs also obeys the same general role, i.e. the band gap of a 2D-TMD is larger than that of bulk.

Here, we turn our attention to the main focus of this paper. The dielectric constants extracted from QE represent the combined dielectric constant of 2D-ice and surroundings large vacuum with a height “$c$”. In order to distill the dielectric constant of 2D-ices, we eliminate the contribution of the vacuum using a capacitance model ¹⁴. In fact, in the out-of-plane (in-plane) direction, the capacitance of the supercell extracted directly from QE code ($\epsilon_{SC}$) is the series (parallel) combination of vacuum capacitance and the 2D-ice capacitance. This helps us to find the out-of-plane and in-plane relative dielectric constant of 2D-ices using below equations ¹⁴:

“t” is the effective thickness of 2D-ice.

In order to make sure that the used “c” values are sufficiently large and the obtained $\varepsilon^{xx,yy,zz}$ were converged well, we performed several additional calculations and found corresponding $\varepsilon_{SC}$. The results show that when “c” values change from 12Å  to 40Å,  the $\varepsilon^{xx,yy,zz}$ values only increase about 5$\%$. Thus, the used vacuum size (see Table I) are sufficiently large.

Using the aforementioned correction, the results for electronic ($\varepsilon_{el}$), ionic ($\varepsilon_{ion}$), and dipolar ($\varepsilon_{dip}$) dielectric constant for 2D-ices are listed in Table I. The main findings are illustrated in Fig. 5 and are listed as:

• •

  $\varepsilon_{ion}^{xx}\simeq\varepsilon_{ion}^{yy}$ for 2D-ices except for b-RH ice in $\varepsilon_{ion}^{xx}\simeq 10\varepsilon_{ion}^{yy}$, which is due to the remarkable in-plane anisotropic in the lattice structure of b-RH.

• •

  $\varepsilon_{el}^{xx}\simeq\varepsilon_{el}^{yy}$ for 2D- (bulk) ices and larger (equal to) than $\varepsilon_{el}^{zz}$. The larger the $\varepsilon_{el}^{xx}$ and $\varepsilon_{el}^{yy}$, the more flatness of the structures. This is due to the confined electron clouds in quasi 2D-space.

• •

  $\varepsilon_{el}^{xx}$ and $\varepsilon_{el}^{yy}$ of f-SQ ice are larger than those for other 2D-ices. This might be due to the confining electron in 2D-plane and stronger response of system to an in-plane electric field, where for other 2D-ices there is a small buckling in their geometry.

• •

  Because of the crystalline structure of all the studied crystalline 2D-ices and bulk ices their dipolar contributions are remarkably smaller than that of bulk water.

• •

  The electronic dielectric constants of 2D- and bulk ices are smaller than for other 2D-materials such as MoS₂, WSe₂, and WS₂ ¹⁴. The latter is attributed to the semiconducting nature of the TMDs rather than insulating nature of 2D-ices.

The items aforementioned are represented in Figs. 7(a,b). Consequently the total dielectric constant (Fig. 6) of the studied ices (except for b-RH and ice VIII) are smaller than 10.

In Fig. 7, the real ($\varepsilon^{r}$) and imaginary ($\varepsilon^{i}$) parts of the optical dielectric function are shown for 2D- and bulk ices and bulk water. The results show that:

• •

  For 2D-ices, $\varepsilon^{xx}\neq\varepsilon^{zz}$ in both real and imaginary parts which is expectable.

• •

  Because of isotropic (anisotropic) lattice of the f-SQ and b-SQ ices (b-RH and HEX) for both real and imaginary parts satisfy. $\varepsilon_{el}^{xx}(\omega)=\varepsilon_{el}^{yy}(\omega)(\varepsilon_{el}^{xx}(\omega)\neq\varepsilon_{el}^{yy}(\omega)$).

• •

  Despite the bulk cubic ice structures (ice VIII and ice Ic), the bulk hexagonal ice structures (ice XI and ice Ih), the dielectric functions are almost equal, i.e. $\varepsilon_{el}^{xx}(\omega)\simeq\varepsilon_{el}^{yy}(\omega)\simeq\varepsilon_{el}^{zz}(\omega)$, for both real and imaginary parts.

The absorption behavior of ices can be understood by analysing the imaginary part of dielectric function. Also, the absorption edge gives the optical gap ($\Delta_{o}$). Moreover, the peaks in the $\varepsilon^{i}(\omega)$ function is related to interband transitions in electronic band structure. For instance, the first peak corresponds to the energy gap ($\Delta$). Indeed, the first critical point calculated in $\varepsilon_{el}^{i}(\omega)$ is related to the transition from valence band maximum to the conduction band minimum, i.e. the energy band gap. The green-solid and pink-dashed vertical arrows in Fig. 7 denote the optical gap and energy gap. In general, the optical gap is equal to the electronic band gap minus the exciton binding energy. In other words, the optical gap and electronic gap should be equal if the many-body perturbation theory is not taken into account. Therefore, our obtained difference between optical gap and electronic band gap is due to the smearing applied in the implemented Kubo-Greenwood approach ³⁸. Moreover, note that $\varepsilon^{r}(\omega=0)$ in Figs. 7 and $\varepsilon_{el}^{0}$ listed in Table I, apparently, are different. However, the difference originates in the two different methods that we employed to calculate them. The first one is the microscopic dielectric function that obtained by using RPA and the second one is the macroscopic dielectric function which was obtained by using DFPT. Thus, they coincide if the same method are used to obtain them.

To compare the optical properties of 2D- and bulk ices, we demonstrated in Fig. 8 (Fig. 9) the real (imaginary) dielectric functions. The results show that:

• •

  The real part of dielectric function has larger values in bulk ices as compared to 2D-ices in low energy region. This is due to the fact that electrons are distributed in 3D space in bulk ice. Moreover, the peaks are larger in 2D-ices as compared to the bulk ices.

• •

  The negative values of real dielectric function of bulk ices correspond to the largest electromagnetic wave reflection. In fact the inductive properties will dominate in this range of energies with the negative values of real dielectric function. In simple words, at energy ranges (e.g. in Ice VIII the energy range 15-20 eV has the negative dielectric function) the electric displacement vector and the electric field vector have opposite directions.

• •

  There is a redshift (move toward lower energy region of the major peaks of $\varepsilon_{el}^{i}(\omega)$ for bulk ices in comparison to 2D-ices. In all studied ice, the prominent peaks in $\varepsilon_{el}^{i}(\omega)$ correspond to optical transmission which are mainly due to the interband transition from the $p$ valence bands to $s$ conduction bands (i.e. the O-$p$ orbital below Fermi level -valence band- and H-s orbital above Fermi level -conduction band-). PDOS of each system is shown in the RHS panels of Fig. 2.

• •

  The absorption energy ranges for 2D- and bulk ices are in the ultraviolet spectra ($>3.2eV$) and visible spectra (between 2 and 3.2 eV), respectively.

The obtained dielectric functions for ice VIII, ice XI,ice Ic, ice Ih, and bulk water are consistent with the results of Refs. [32,39], Ref. [35], and Refs. [40,41], respectively. Note that $\varepsilon^{r}(\omega=0)$ in Figs. 7 and $\varepsilon_{el}^{0}$ listed in Table I, apparently, are different. However, the difference originates in the two different methods that we employed to calculate each of them. The first one corresponds to the microscopic dielectric function which was obtained by using RPA and the second one corresponds to the macroscopic dielectric function which was obtained by using DFPT.

It is insightful to calculate mechanical stiffness of a typical 2D-ice (the results are selectively presented for f-SQ). We performed an additional calculation for determining the Young’s modulus and Poisson’s ratio of f-SQ ice. We applied both uniaxial and biaxial strains found the corresponding energies $E_{u}$ and $E_{b}$, for uniaxial and biaxial strained systems. Next, the best fits on two equations $E_{b}=2bu^{2}_{b}$ and $E_{u}=\frac{1}{2}(b+\mu)u^{2}_{u}$, helped us to obtain 2D-bulk modulus “b” and 2D-shear modulus “$\mu$”. Here $u_{u}$ and $u_{b}$ are the applied uniaxial and biaxial strains, respectively. Consequently, the Young’s modulus ($Y_{2d}$) and Poisson’s ratio ($\nu_{2d}$) can be computed:  ⁴²

The results are $Y_{f-SQ}$=3.6 N/m and $\nu=0.6$. The Young’s modulus of 2D-ice is two orders of magnitude smaller than the one for graphene, i.e. 340 N/m ⁴³. This is due to the weak hydrogen bonds in 2D-ices compared to the strong planar covalent bonds in graphene. The corresponding Poisson’s ratio of 2D-ice is larger than graphene ( 0.3). The obtained Young’s modulus (Poisson’s ratio) for 2D-ice is more or less equal to bulk modulus of bulk ice. In Table II, we listed Young’s modulus/bulk modulus and Poisson’s ratio of f-SQ/bulk crystalline ice. To convert N/m unit in $Y_{2d}$ to Pas unit, one may use the simple relation $Y_{3d}=Y_{2d}/t_{0}$ where $t_{0}$=3Å is the effective thickness of f-SQ ice (see Table I).