Thermal fluctuations in crystalline membranes with long-range dipole interactions

TL;DR

This paper investigates how long-range dipole interactions influence the thermal fluctuations of crystalline membranes, showing they mainly renormalize elastic properties without altering universal scaling laws.

Contribution

It demonstrates that dipole-dipole interactions do not change the fundamental scaling behavior of membranes, unlike scalar interactions, but significantly affect elastic coefficients.

Findings

Dipole interactions renormalize elastic coefficients.

Universal scaling laws remain valid with dipole interactions.

Scalar interactions lead to different fluctuation scaling.

Abstract

We study the effects of long-range electrostatic interactions on the thermal fluctuations of free-standing crystalline membranes exhibiting spontaneous electric polarization directed at each point along the local normal to the surface. We show that the leading effect of dipole-dipole interactions in the long-wavelength limit consists in renormalizations of the bending rigidity and the elastic coefficients. A completely different result was obtained in the case of scalar two-point interactions decaying as $R^{-3}$, where $R$ is the distance. In the latter case, which was addressed in previous theoretical research, the energy of long-wavelength bending fluctuations is controlled by power-law interactions and it scales with the wavevector $q$ as $q^{3}$, leading to a modified large-distance behaviour of correlation functions. By contrast, in the case of dipole interactions, the $q^{3}$ dependence of the bending energy vanishes. Non-local terms generated by the expansion of the electrostatic energy are suppressed in the limit of small wavevectors. This suggests that the universal scaling behaviour of elastic membranes holds even in presence of dipole interactions. At the same time, the shift of the Lam\'e coefficients and the bending rigidity induced by electrostatic interactions can be quantitatively important for two-dimensional materials with a permanent out-of-plane polarization.

Equations

$$H_{\text{el}}=\int {\mathrm{d}}^{2}x\left[\frac{1}{2}\kappa {\left({\nabla }^{2}h\right)}^2+\frac{\lambda }{2}{\tilde{u}}_{\alpha \alpha }^2+\mu {\tilde{u}}_{\alpha \beta }{\tilde{u}}_{\alpha \beta }\right],$$

$${\tilde{u}}_{\alpha \beta }=\frac{1}{2}({\partial }_{\alpha }{u}_{\beta }+{\partial }_{\beta }{u}_{\alpha }+{\partial }_{\alpha }h{\partial }_{\beta }h)$$

$$V(\mathbf{R})=A/{\left|\mathbf{R}\right|}^{\sigma }$$

$$\begin{split}H_{\text{sr}}&=\int{\rm d}^{2}x\bigg{[}\frac{1}{2}\kappa\left(\nabla^{2}\bm{r}\right)^{2}\\ &+\frac{\lambda}{2}u_{\alpha\alpha}^{2}+\mu\,u_{\alpha\beta}u_{\alpha\beta}+\tau u_{\alpha\alpha}\bigg{]}~{},\\ \end{split}$$

$$\begin{array}{rl}{u}_{\alpha \beta }& =\frac{1}{2}\left({\partial }_{\alpha }\mathbf{r}\cdot {\partial }_{\beta }\mathbf{r}-{\delta }_{\alpha \beta }\right)\\ & =\frac{1}{2}\left({\partial }_{\alpha }{u}_{\beta }+{\partial }_{\beta }{u}_{\alpha }+{\partial }_{\alpha }h{\partial }_{\beta }h+{\partial }_{\alpha }{u}_{\gamma }{\partial }_{\beta }{u}_{\gamma }\right).\end{array}$$

$$\int {\mathrm{d}}^{2}x\tau u_{\alpha \alpha }\left(\mathbf{x}\right)$$

$$-\sum _{\mathbf{q}}{q}^2{\left|h_{\mathbf{q}}\right|}^2.$$

$$\begin{split}H_{\text{sr}}&=\int{\rm d}^{2}x\bigg{[}\frac{1}{2}\kappa\left(\nabla^{2}h\right)^{2}\\ &+\frac{\lambda}{2}\tilde{u}_{\alpha\alpha}^{2}+\mu\,\tilde{u}_{\alpha\beta}\tilde{u}_{\alpha\beta}+\tau u_{\alpha\alpha}\bigg{]}~{},\\ \end{split}$$

$$\begin{split}&H=H_{\text{sr}}+H_{\text{lr}}~{},\\ &H_{\text{lr}}=\frac{1}{2}P^{2}\int{\rm d}^{2}x\int{\rm d}^{2}x^{\prime}\bigg{[}\frac{\bm{n}(\bm{x})\cdot\bm{n}(\bm{x}^{\prime})}{\left|\bm{r}(\bm{x})-\bm{r}(\bm{x}^{\prime})\right|^{3}}\\ &\;\;\;\;\;\;-\frac{3\left(\bm{n}(\bm{x})\cdot\bm{\nu}\right)\left(\bm{n}(\bm{x}^{\prime})\cdot\bm{\nu}\right)}{\left|\bm{r}(\bm{x})-\bm{r}(\bm{x}^{\prime})\right|^{3}}\bigg{]}~{}.\end{split}$$

$$\mathbf{\nu }\equiv \frac{\mathbf{r}(\mathbf{x})-\mathbf{r}({\mathbf{x}}^{\prime })}{\left|\mathbf{r}(\mathbf{x})-\mathbf{r}({\mathbf{x}}^{\prime })\right|}.$$

$$\begin{split}H_{\text{lr}}&=\frac{1}{2}\sum_{i}\sum_{j\neq i}p_{i}p_{j}\bigg{[}\frac{\bm{n}_{i}\cdot\bm{n}_{j}}{\left|\bm{r}_{i}-\bm{r}_{j}\right|^{3}}\\ &-\frac{3\left(\bm{n}_{i}\cdot\bm{\nu}\right)\left(\bm{n}_{j}\cdot\bm{\nu}\right)}{\left|\bm{r}_{i}-\bm{r}_{j}\right|^{3}}\bigg{]}~{}.\end{split}$$

$$\begin{array}{rl}& ⟨{\left(h\left(\mathbf{x}\right)-h\left({\mathbf{x}}^{\prime }\right)\right)}^2⟩\approx {\left|\mathbf{x}-{\mathbf{x}}^{\prime }\right|}^{2-\eta }\\ & ⟨{\left(\mathbf{u}\left(\mathbf{x}\right)-\mathbf{u}\left({\mathbf{x}}^{\prime }\right)\right)}^2⟩\approx {\left|\mathbf{x}-{\mathbf{x}}^{\prime }\right|}^{2-2\eta },\end{array}$$

$$\begin{split}H_{\text{lr}}&=\frac{1}{2}\sum_{i}\sum_{j\neq i}p_{i}p_{j}\frac{1}{\left|\bm{r}_{i}-\bm{r}_{j}\right|^{3}}\\ &-\frac{1}{2}\sum_{i}\sum_{j\neq i}p_{i}p_{j}\bigg{[}\frac{\frac{1}{2}\left(\bm{n}_{i}-\bm{n}_{j}\right)^{2}+3\left(\bm{n}_{i}\cdot\bm{\nu}\right)\left(\bm{n}_{j}\cdot\bm{\nu}\right)}{\left|\bm{r}_{i}-\bm{r}_{j}\right|^{3}}\bigg{]}\\ &\equiv H_{\text{lr},1}+H_{\text{lr},2}~{}.\end{split}$$

$$\mathbf{n}(\mathbf{x})={\mathbf{e}}_z-\nabla h(\mathbf{x}).$$

$$\begin{split}H_{\text{lr},2}&\simeq H_{\text{lr},2}^{(2,0)}=-\frac{1}{2}P^{2}\int{\rm d}^{2}x\int{\rm d}^{2}x^{\prime}\bigg{[}\\ &\frac{\frac{1}{2}\left(\bm{\nabla}h(\bm{x})-\bm{\nabla}h(\bm{x}^{\prime})\right)^{2}}{\left|\bm{x}-\bm{x}^{\prime}\right|^{3}}\\ &+\frac{3\left(h(\bm{x})-h(\bm{x}^{\prime})-\bm{\nabla}h(\bm{x})\cdot\left(\bm{x}-\bm{x}^{\prime}\right)\right)}{\left|\bm{x}-\bm{x}^{\prime}\right|^{5}}\\ &\times\left(h(\bm{x})-h(\bm{x}^{\prime})-\bm{\nabla}h(\bm{x}^{\prime})\cdot\left(\bm{x}-\bm{x}^{\prime}\right)\right)\bigg{]}~{}.\\ \end{split}$$

$$h(\mathbf{x})=\frac{1}{\sqrt{A}}\sum _{\mathbf{q}}{h}_{\mathbf{q}}{e}^{i\mathbf{q}\cdot \mathbf{x}},$$

$$H_{\text{lr},2}^{(2,0)}=-\frac{\pi }{3}{P}^2\sum _{\mathbf{q}}{\left|\mathbf{q}\right|}^3{\left|h_{\mathbf{q}}\right|}^2.$$

$$H_{\text{lr},2}^{(2,0)}=P^{2}\sum_{\bm{q}}\Big{[}-\frac{\pi}{3}\left|\bm{q}\right|^{3}+\frac{1}{a^{3}}f_{2}\left(a\bm{q}\right)\Big{]}\left|h_{\bm{q}}\right|^{2}~{},$$

$$\begin{split}&H_{\text{lr},1}\simeq H_{\text{lr},1}^{(0,0)}+H_{\text{lr},1}^{(2,0)}+H_{\text{lr},1}^{(0,1)}+H_{\text{lr},1}^{(0,2)}\\ &=\frac{1}{2}\sum_{i}\sum_{j\neq i}p_{i}p_{j}\bigg{[}\frac{1}{\left|\bm{x}_{i}-\bm{x}_{j}\right|^{3}}-\frac{3}{2}\frac{\left(h_{i}-h_{j}\right)^{2}}{\left|\bm{x}_{i}-\bm{x}_{j}\right|^{5}}\\ &-\frac{3}{2}\;\frac{2\left(\bm{u}_{i}-\bm{u}_{j}\right)\cdot\left(\bm{x}_{i}-\bm{x}_{j}\right)+\left(\bm{u}_{i}-\bm{u}_{j}\right)^{2}}{\left|\bm{x}_{i}-\bm{x}_{j}\right|^{5}}\\ &+\frac{15}{2}\frac{\left(\left(\bm{u}_{i}-\bm{u}_{j}\right)\cdot\left(\bm{x}_{i}-\bm{x}_{j}\right)\right)^{2}}{\left|\bm{x}_{i}-\bm{x}_{j}\right|^{7}}\bigg{]}~{}.\\ \end{split}$$

$$\begin{split}&H_{\text{lr},1}\simeq\frac{1}{2}\frac{b_{3}}{a}P^{2}A-\frac{3}{2}\sum_{i}\sum_{j\neq i}p_{i}p_{j}\frac{\left(\bm{u}_{i}-\bm{u}_{j}\right)\cdot\left(\bm{x}_{i}-\bm{x}_{j}\right)}{\left|\bm{x}_{i}-\bm{x}_{j}\right|^{3}}\\ &+P^{2}\sum_{\bm{q}}\bigg{[}-\frac{3}{8}\frac{b_{3}}{a}\,q^{2}+\frac{\pi}{3}\left|\bm{q}\right|^{3}+\frac{1}{a^{3}}f_{1}\left(a\bm{q}\right)\bigg{]}\left|u_{\alpha,\bm{q}}\right|^{2}\\ &+P^{2}\sum_{\bm{q}}\bigg{[}-\frac{3}{8}\frac{b_{3}}{a}\,q^{2}+\frac{\pi}{3}\left|\bm{q}\right|^{3}+\frac{1}{a^{3}}f_{1}\left(a\bm{q}\right)\bigg{]}\left|h_{\bm{q}}\right|^{2}\\ &+P^{2}\sum_{\bm{q}}\bigg{[}\frac{15}{16}\frac{b_{3}}{a}q^{2}\left(\frac{1}{2}\delta_{\alpha\beta}+\frac{q_{\alpha}q_{\beta}}{q^{2}}\right)\\ &-\pi\left|\bm{q}\right|^{3}\Big{(}\frac{1}{3}\delta_{\alpha\beta}+\frac{q_{\alpha}q_{\beta}}{q^{2}}\Big{)}+\frac{1}{a^{3}}g_{\alpha\beta}\Big{(}a\bm{q}\Big{)}\bigg{]}u_{\alpha,\bm{q}}u_{\beta,-\bm{q}}~{}.\end{split}$$

$$b_3\equiv \frac{a}{P^2}\frac{1}{A}\sum _i\sum _{j\ne i}{p}_i{p}_j\frac{1}{{\left|{\mathbf{x}}_i-{\mathbf{x}}_j\right|}^3}.$$

$$\begin{array}{rl}& -\frac{3}{2}\sum _i\sum _{j\ne i}{p}_i{p}_j\frac{\left({\mathbf{u}}_i-{\mathbf{u}}_j\right)\cdot \left({\mathbf{x}}_i-{\mathbf{x}}_j\right)}{{\left|{\mathbf{x}}_i-{\mathbf{x}}_j\right|}^3}\\ & \approx -\frac{3}{4}\frac{b_3}{a}{P}^2\int {\mathrm{d}}^{2}x{\partial }_{\alpha }{u}_{\alpha }\left(\mathbf{x}\right).\end{array}$$

$$\begin{split}H_{\text{lr}}&=H_{\text{lr},1}+H_{\text{lr},2}\\ &\simeq\frac{1}{2}\frac{b_{3}}{a}P^{2}A-\frac{3}{4}\frac{b_{3}}{a}P^{2}\int{\rm d}^{2}x\,u_{\alpha\alpha}\left(\bm{x}\right)\\ &+\int{\rm d}^{2}x\bigg{[}\frac{1}{2}\kappa_{1}\left(\nabla^{2}h\right)^{2}+\frac{\lambda_{1}}{2}\left(\partial_{\alpha}u_{\alpha}\right)^{2}\\ &+\frac{\mu_{1}}{4}\left(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}\right)\left(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}\right)\bigg{]}\\ &-\pi P^{2}\sum_{\bm{q}}\left|\bm{q}\right|q_{\alpha}q_{\beta}\,u_{\alpha,\bm{q}}u_{\beta,-\bm{q}}~{},\end{split}$$

$$\sum _{l,m}\frac{m^2}{{\left(l^2+m^2-lm\right)}^{\frac{5}{2}}}=\frac{2}{3}\sum _{l,m}\frac{1}{{\left(l^2+m^2-lm\right)}^{\frac{3}{2}}},$$

$$-\frac{3}{8}\frac{b_3}{a}{P}^2\sum _{\mathbf{q}}{q}^2\left({\left|h_{\mathbf{q}}\right|}^2+u_{\alpha ,\mathbf{q}}{u}_{\alpha ,-\mathbf{q}}\right)$$

$$\tau =\frac{3}{4}\frac{b_3}{a}{P}^2,$$

$$\begin{split}H&=H_{\text{sr}}+H_{\text{lr}}\\ &\simeq\int{\rm d}^{2}x\,\bigg{[}\frac{1}{2}\kappa^{\prime}\left(\nabla^{2}h\right)^{2}+\frac{\lambda^{\prime}}{2}\left(\partial_{\alpha}u_{\alpha}\right)^{2}\\ &+\frac{\mu^{\prime}}{4}\left(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}\right)\left(\partial_{\alpha}u_{\beta}+\partial_{\beta}u_{\alpha}\right)\bigg{]}\\ &-\pi P^{2}\sum_{\bm{q}}\left|\bm{q}\right|q_{\alpha}q_{\beta}u_{\alpha,\bm{q}}u_{\beta,-\bm{q}}~{},\end{split}$$

$$\sum _{\mathbf{q}}{q}^2{\left|h_{\mathbf{q}}\right|}^2=\int {\mathrm{d}}^{2}x{\left(\nabla h\left(\mathbf{x}\right)\right)}^2,$$

$${\mathcal{O}}_h=\sum _{\mathbf{q}}{\left|\mathbf{q}\right|}^3{\left|h_{\mathbf{q}}\right|}^2$$

$$H_{\text{lr}}=\frac{1}{2}\int \frac{{\mathrm{d}}^{3}Q}{{\left(2\pi \right)}^3}\frac{4\pi }{Q^2}{\rho }_{\mathbf{Q}}{\rho }_{-\mathbf{Q}},$$

Full text

I Introduction

The thermodynamic properties and the elastic response of two-dimensional (2D) crystalline membranes subject to small or vanishing external tension are crucially determined by the thermal fluctuations of their shape. The statistical mechanics of fluctuating elastic membranes has been the subject of extensive research effort Nelson et al. (2004); Katsnelson (2012); Bowick and Travesset (2001); Nelson and Peliti (1987); Aronovitz and Lubensky (1988); Guitter et al. (1989); Le Doussal and Radzihovsky (1992, 2018); Gazit (2009); Kownacki and Mouhanna (2009); Košmrlj and Nelson (2013, 2016); Burmistrov et al. (2016); Gornyi et al. (2015, 2017); Katsnelson and Fasolino (2013); Los et al. (2017, 2016); Burmistrov et al. (2018). The problem is particularly complex because of low dimensionality and the softness of the dispersion relation of bending modes, which leads to strong fluctuations and the breakdown of the harmonic approximation. The energy and the nature of fluctuations depend essentially on the degree of internal order of the membrane Nelson et al. (2004). In the case of fluid membranes, the shear modulus vanishes. It was demonstrated that bending fluctuations then lead to a crumpled state at any finite temperature Nelson et al. (2004). By contrast, membranes with an internal crystalline order and, therefore, a finite shear modulus present a low-temperature flat phase, characterized by a planar average configuration. The stability of the flat phase long-range order to thermal bending fluctuations is made possible by the anharmonic coupling of out-of-plane and in-plane modes Nelson and Peliti (1987). The statistical mechanical properties of the flat phase are controlled by long-wavelength fluctuations of out-of-plane (bending) and in-plane (stretching) modes. In absence of long-range interactions and anisotropies in the elastic coefficients, the most general effective Hamiltonian for slowly varying deformations is Nelson et al. (2004); Nelson and Peliti (1987); Aronovitz and Lubensky (1988); Guitter et al. (1989); Le Doussal and Radzihovsky (1992)

[TABLE]

where $h(\bm{x})$ and $\bm{u}(\bm{x})$ are out-of-plane and in-plane displacements and

[TABLE]

is the strain tensor. $\lambda$ and $\mu$ are Lamé coefficients and $\kappa$ is the bending rigidity. Terms including higher powers of $\partial_{\alpha}h$ or $\partial_{\alpha}u_{\beta}$ or higher-order gradients can be neglected, because they are irrelevant at long wavelengths Nelson and Peliti (1987); Aronovitz and Lubensky (1988); Guitter et al. (1989); not . The field $h(\bm{x})$ plays the role of the Goldstone mode related to broken rotational symmetry Guitter et al. (1989); Le Doussal and Radzihovsky (2018); Bowick and Travesset (2001) and it exhibits a particularly soft dispersion relation. In the harmonic approximation, the energy of bending fluctuations scales as $\kappa q^{4}$, where $q$ is the wavevector. The anharmonic coupling between in-plane and out-of-plane modes constitutes a strongly relevant perturbation to the Gaussian fixed point, making the statistical mechanics very non-trivial Nelson et al. (2004). Correlation functions of $h$ and $\bm{u}$ fields can be calculated perturbatively in the anharmonic coupling strength for wavevectors larger than the ’Ginzburg scale’ $q_{*}\approx\sqrt{3Y/16\pi\kappa^{2}}$, with $Y=4\mu(\lambda+\mu)/(\lambda+2\mu)$. For $q\ll q_{*}$, instead, perturbation theory breaks down Nelson and Peliti (1987); Katsnelson (2012); Katsnelson and Fasolino (2013). The long-wavelength behaviour of correlation functions is determined by a non-trivial fixed point Hamiltonian, characterized by anomalous scaling exponents Nelson et al. (2004); Bowick and Travesset (2001); Nelson and Peliti (1987); Aronovitz and Lubensky (1988); Guitter et al. (1989); Kownacki and Mouhanna (2009). In the small-$q$ regime, the bending rigidity diverges as $q^{-\eta}$ and the elastic Lamé coefficients both vanish as $q^{2-2\eta}$, where $\eta$ is an universal critical exponent. The determination of $\eta$ has been addressed through several field-theoretic methods Nelson et al. (2004); Aronovitz and Lubensky (1988); Guitter et al. (1989); Le Doussal and Radzihovsky (1992, 2018); Gazit (2009); Kownacki and Mouhanna (2009). The self-consistent screening approximation Le Doussal and Radzihovsky (1992, 2018) and the functional renormalization group Kownacki and Mouhanna (2009) lead to $\eta\simeq 0.821$ and $\eta\simeq 0.849$, respectively.

The anomalous scaling of the elastic and bending coefficients has crucial consequences in the physical properties of crystalline membranes. The very stability of the flat phase and the existence of long-range order in the orientation of the normal to the membrane relies on the power-law stiffening of the bending rigidity and would be impossible if anharmonic effects were neglected Nelson and Peliti (1987). In addition, the renormalization of elastic coefficients determines the macroscopic mechanical properties of crystalline membranes in an essential way, leading to anomalous strain and size dependent elastic response Los et al. (2017, 2016); Košmrlj and Nelson (2013, 2016); Katsnelson and Fasolino (2013); Gornyi et al. (2017); Burmistrov et al. (2018).

The statistical mechanics of anharmonic thermal fluctuations of 2D crystals in the flat phase has fundamental importance in the theory of thermodynamic and mechanical properties of free-standing graphene or suspended samples of other atomically-thin materials Katsnelson (2012); Katsnelson and Fasolino (2013); Los et al. (2017). The Ginzburg wavelength $\lambda_{*}=2\pi/q_{*}$ for these materials at room temperature is of the order of 10 Å, which makes the anomalous scaling behaviour important already at short length scales Le Doussal and Radzihovsky (2018). The elasticity of graphene is currently under extensive investigation, both theoretically Los et al. (2016, 2017); Košmrlj and Nelson (2016) and experimentally Nicholl et al. (2015); López-Polín et al. (2017).

It has been understood long time ago that the scaling behaviour of crystalline membranes can be modified in presence of long-range interactions which decay sufficiently slowly with the distance. The effects of power-law interactions of the form

[TABLE]

on $D$-dimensional membranes fluctuating in a $d$-dimensional embedding space were studied in Refs. Toner (1989); Kantor and Kardar (1989); Guitter and Palmeri (1992); Mori and Wadati (1993). In Ref. Toner (1989), in particular, the author focused on effects of long-range interactions in the flat phase and addressed the scaling behaviour for generic $D$, $d$ and $\sigma$. It was shown that, for small fluctuations around the flat configuration, the dominant effect of long-range interactions is to modify the long-wavelength behaviour of the quadratic part of the bending energy, which acquires a wavevector dependence $\left|\bm{q}\right|^{\sigma+2-D}$. When $\sigma+2-D<4-\eta$, where $\eta$ is the short-range critical exponent, long-range interactions dominate the long-wavelength behaviour and the universal anomalous elasticity predicted by conventional membrane theory does not hold Toner (1989). If $\sigma$ is decreased further, $\sigma<\frac{3}{2}D$, the bending stiffness becomes larger and the exponents approach the Gaussian ones. The case of scalar interactions decaying as $\left|\bm{R}\right|^{-3}$ in two dimensions, lies precisely at the boundary of the domain of validity of the Gaussian approximation. The renormalized in-plane elastic coefficients $\lambda(q)$ and $\mu(q)$ then vanish as 1\big{/}\ln\left(1/q\right) Toner (1989).

The fact that long-range interactions can crucially affect the behaviour of systems with massless modes is very general. In two-dimensional spin systems, for example, dipole-dipole interactions lead to the breaking of the conditions of the Mermin-Wagner theorem and have dramatic effects on thermodynamic properties Grechnev et al. (2005). Long-range dipolar interactions also play a crucial role in spin systems at the ferromagnetic phase transitions, leading to a modified critical behaviour Aharony and Fisher (1973); Fisher et al. (1972).

Electrostatic interactions also play an important role in charged fluid membranes immersed in electrolyte solutions. Vast research effort was devoted to the quantitative prediction of the renormalization of bending moduli induced by electrostatic interactions Daicic et al. (1996); Fogden and Daicic (1997); Fogden et al. (1996); Kumaran (2001); Loubet et al. (2013); Netz (2001); Shojaei et al. (2016).

Several theoretical studies also addressed the curvature moduli of membranes subject to van der Waals interactions Netz (2001); Shojaei et al. (2016); Dean and Horgan (2006) and the elasticity of flexoelectric membranes Peliti and Prost (1989); Liu and Sharma (2013).

In this paper, we discuss the effect of dipole-dipole interactions on the fluctuations of free-standing crystalline membranes. We focus, in particular, on the case of two-dimensional crystals with a permanent polarization of fixed magnitude oriented along the local normal, which we describe as arising from point dipole moments located at the lattice sites. We show that dipole-dipole interactions lead to a very different result from the one obtained for scalar interactions decreasing with distance as $\left|\bm{R}\right|^{-3}$. In the quadratic part of the energy, the non-analytic wavevector dependence proportional to $\left|\bm{q}\right|^{\sigma+2-D}=\left|\bm{q}\right|^{3}$ of the inverse $h$ field propagator vanishes exactly. Expanding the electrostatic energy generates, in the leading order in the long-wavelength limit, to a local functional of the strain tensor. Non-local terms, instead, constitute small corrections in the limit of long wavelengths and they are not expected to modify the large-distance behaviour. Electrostatic interactions, on the other hand, lead to a shift of the Lamé coefficients $\lambda$ and $\mu$ and to the bending rigidity $\kappa$, which can be quantitatively important in realistic two-dimensional materials.

Our study is relevant for 2D materials which lack symmetry under inversion and under reflection with respect to the crystal plane and which can exhibit, therefore, a finite out-of-plane polarization. Important examples are the graphene derivatives C2F (fluorinated on one side), C2HF (fluorinated on one side and hydrogenated on the opposite side) and C2LiF Klintenberg et al. (2010). In these materials, polarization arises from the charge transfer between carbon and halogen and alkali atoms. Other examples of two-dimensional materials exhibiting spontaneous out-of-plane polarization are MoSSe monolayer and multilayers Riis-Jensen et al. (2018). Recently, out-of-plane spontaneous polarization and ferroelectric phenomena were observed in bi- and trilayer WTe2, which were identified as rare examples of 2D ferroelectric metals Fei et al. (2018).

Thermal fluctuations of membranes with a spontaneous polarization directed along the normal to the surface were investigated in the previous work Peliti and Prost (1989). However, in Ref. Peliti and Prost (1989), the expansion of the electrostatic energy up to second order in $h(\bm{x})$ was calculated only approximately. This lead to the prediction of a $\left|\bm{q}\right|^{3}$ wavevector dependence of the bending energy, which contrasts with our results. Fluctuations of nearly-planar dipole layers were also investigated in Ref. Miller and Scriven (1974), which reports the calculation of the electrostatic energy up to second order in the out-of-plane fluctuation amplitude and to order $q^{2}$ for a dipole layer at the interface between media with two different dielectric constants $\epsilon_{a}$ and $\epsilon_{b}$. In our work we determined, in the case $\epsilon_{a}=\epsilon_{b}$, the fluctuation energy to higher orders in $\bm{q}$ and in the fluctuation amplitudes, which is crucial for addressing the statistical mechanical properties of the system.

A comparison between our calculations and the results in Refs. Daicic et al. (1996); Fogden et al. (1996); Fogden and Daicic (1997); Kumaran (2001); Loubet et al. (2013); Netz (2001); Shojaei et al. (2016); Dean and Horgan (2006) is not direct because of the differences of the systems under examination. In the case of charged membranes immersed in solutions with added salt, Coulomb interactions are Debye-screened and become effectively short-ranged. In this case, effects of electrostatic interactions can a priori be encoded in a renormalization of local elastic coefficients, in the limit of long wavelengths Daicic et al. (1996); Fogden et al. (1996); Kumaran (2001); Loubet et al. (2013). On the other hand, electrostatic interactions are unscreened in the case of membranes immersed in solutions without added salt, in which the only mobile ions are counterions which ensure an overall charge neutrality Fogden et al. (1996); Fogden and Daicic (1997). In Refs. Fogden et al. (1996); Fogden and Daicic (1997), bending fluctuations of two charged membranes with the same charge and with intervening counterions were addressed. It was found that, for in-phase fluctuation modes of the two layers, the wavevector dependence of the energy is regular at $\bm{q}\rightarrow 0$ despite the long-range nature of the interactions. This finding is similar to our results for the quadratic part of the energy of out-of-plane fluctuations in a membrane with dipole-dipole interactions, although the two systems and the methods to describe them are very different.

II Hamiltonian for crystalline membranes with long-range interactions

Fluctuating configurations of a crystalline membrane in the flat phase are described by specifying the displacement of all points from their position in the minimal energy configuration. In the rest minimum energy state, the membrane lies in a plane, which we can choose as the $xy$ plane. Displacements in the membrane plane are denoted as $\bm{u}\left(\bm{x}\right)$ and displacements in the orthogonal $z$ direction are denoted as $h\left(\bm{x}\right)$. The position of points in the three-dimensional embedding space is denoted as $\bm{r}(\bm{x})=\bm{x}+\bm{u}(\bm{x})+h(\bm{x})\bm{e}_{z}$. When the crystal lattice is sufficiently isotropic and interactions are short-ranged, long-wavelength fluctuations can be described through the coarse-grained Hamiltonian Nelson et al. (2004); Nelson and Peliti (1987); Aronovitz and Lubensky (1988); Guitter et al. (1989); Gornyi et al. (2015):

[TABLE]

where $\kappa$ is the bending rigidity, $\lambda$ is the first Lamé coefficient and $\mu$ the shear modulus. The strain tensor $u_{\alpha\beta}$ is defined as:

[TABLE]

The linear term

[TABLE]

represents an elastic tension. If only the short-range forces described by Eq. (3) were present, $\tau$ should vanish because $\bm{u}(\bm{x})=h(\bm{x})=0$ was defined as the state of minimum energy. However, dipole-dipole interactions will be included. These tend to stretch the membrane. Therefore, a finite value of $\tau$ must be taken into account in order to ensure that the reference configuration will correspond to the state of minimum energy of the total Hamiltonian. The appearance of a finite elastic tension will be important in the following. As shown in Section II.1 (see Eq. (18)), electrostatic interactions lead in the bending energy to a negative, quadratic-in-q term, proportional to

[TABLE]

This term will be exactly cancelled by the positive tension $\tau$, as it was already noticed in Ref. Toner (1989). The exact cancellation is ensured by rotational symmetry and is related to the ’Goldstone’ nature of the out-of-plane mode $h(\bm{x})$.

Eq. (3) can be simplified by neglecting in-plane anharmonicity and the term $\left(\nabla^{2}\bm{u}\right)^{2}$, which can be shown to be irrelevant for the long-wavelength behaviour Nelson et al. (2004); Nelson and Peliti (1987); Aronovitz and Lubensky (1988); Guitter et al. (1989); Gornyi et al. (2015). If these terms are ignored, $H_{\text{sr}}$ becomes:

[TABLE]

where the approximate strain tensor $\tilde{u}_{\alpha\beta}$ is defined as in Eq. (2).

In presence of dipolar interactions, the Hamiltonian for displacement fields is:

[TABLE]

Here $\bm{n}(\bm{x})$ is the normal to the membrane surface at the point $\bm{x}$ and $\bm{\nu}$ is the unit vector:

[TABLE]

$P$ is the two-dimensional polarization density, defined as the dipole moment per unit area. The electrostatic energy defined in Eq. (7) is actually divergent when $\bm{x}^{\prime}$ approaches $\bm{x}$ and a short-distance cutoff must be introduced. Motivated by the study of ordered crystals such as CF, C2HF or C2LiF we take care of the divergence of the energy by the following procedure. We consider the polarization to arise from finite dipole moments located at the sites of a perfect lattice. The energy becomes:

[TABLE]

The discrete indices $i$ and $j$ label the lattice sites at which dipole moments are located, $p_{i}$ denotes the magnitude of the dipole moment at the site $i$. $p_{i}$ is assumed to have the periodicity of the lattice, so that the dipole moment distribution is identical in all unit cells. For simplicity, despite the discretization of the problem, we still describe configurations through functions $\bm{u}(\bm{x})$ and $h(\bm{x})$ defined on the continuum two-dimensional space. For consistency, the functions $\bm{u}(\bm{x})$ and $h(\bm{x})$ must be slowly varying, with finite Fourier components only for wavevectors well within the first Brillouin zone. In Eq. (9), we identify $\bm{r}_{i}$ and $\bm{n}_{i}$ with $\bm{r}(\bm{x}_{i})$ and $\bm{n}(\bm{x}_{i})$. This corresponds to picturing the membrane as a 2D continuum elastic medium in which dipole moments are embedded at discrete lattice positions Tamm (1979); *purcell_electricity_1965. If the displacement fields vary slowly with respect to the atomic scale, the approximate character of the continuum description is inessential. It allows however a simpler definition of the unit vector $\bm{n}$ normal to the surface and it avoids the necessity to derive a discretized version of the short-range elastic Hamiltonian in Eq. (6).

To study how electrostatic interactions perturb the properties of membranes which lie deep in the flat phase and far from the crumpling transition point, it is legitimate to expand Eq. (9) in powers of the displacement fields $h(\bm{x})$, $\bm{u}(\bm{x})$. According to the conventional statistical mechanics of membranes, if only the elastic energy, Eq. (6), is considered, the fluctuating field correlation functions scale at large distances according to:

[TABLE]

with $\eta\simeq 0.8$. Therefore, distances fluctuate much less than their average value. In addition, the mean square fluctuation of the normal around the average $z$-axis is finite and small if $T\ll\kappa\approx T_{\text{cr}}$. This is the case for graphene and most two-dimensional materials at room temperature.

II.1 Electrostatic energy: harmonic approximation

Expanding Eq. (9) up to second order in the fields gives the long-range contribution to the Gaussian fluctuation theory. It is convenient to rewrite the expression for the energy as:

[TABLE]

The second line in Eq. (11), denoted $H_{\text{lr},2}$, is now regular when $\bm{x}_{i}$ approaches $\bm{x}_{j}$ and the lattice sums can be safely transformed into integrals. To first order in the displacement fields, the normal $\bm{n}(\bm{x})$ can be written as:

[TABLE]

Keeping only terms of second order in $H_{\text{lr},2}$ yields:

[TABLE]

We will denote as $H^{(n,m)}$ the term of order $h^{n}\,u^{m}$ in the expansion of the energy. Introducing the Fourier transform:

[TABLE]

gives:

[TABLE]

An explicit calculation of the discrete lattice sum, based on the Ewald method (see the Appendix), actually gives:

[TABLE]

where $a$ is the lattice constant and $f_{2}$ is a regular function of the components of $\bm{q}$ for $a\bm{q}\rightarrow 0$, vanishing as $\left(aq\right)^{4}$ for small wavevectors. The function $f_{2}$ depends on the geometrical details of the lattice. It is in general anisotropic, although it must be invariant under the symmetry operations of the crystal. We will assume that the lattice structure is symmetric enough to ensure that the only invariant quadratic and quartic functions of $\bm{q}$ are $q^{2}$ and $\left(q^{2}\right)^{2}$. Also, we will suppose that the lattice has some inversion or reflection symmetry so that all invariant functions of $\bm{q}$ must be even for $\bm{q}\rightarrow-\bm{q}$. These are actually the same requirements which underlie the isotropy of the elastic and bending coefficients in Eq. (6).

We now turn to the expansion for small fluctuations of $H_{\text{lr},1}$. Expanding up to second order yields:

[TABLE]

Fourier transformation gives:

[TABLE]

Details of the calculations are presented in the Appendix. In Eq. (18), $b_{3}$ is a dimensionless coefficient defined as:

[TABLE]

$f_{1}$ and $g_{\alpha\beta}$ are regular functions of the components of the wavevector at $\bm{q}=0$. $f_{1}$ and $g_{\alpha\beta}$ are both of order $q^{4}$ at small $\bm{q}$. The linear term in the first line of Eq. (18) represents the tendency of the membrane to stretch under the repulsion between dipole moments. In the continuum limit, it can be written as:

[TABLE]

The electrostatic energy can now be obtained by summing Eqs. (16) and (18). All terms with a regular small momentum dependence can be represented in real space as a local functional of $\partial_{\alpha}h(\bm{x})$, $\partial_{\beta}u_{\alpha}(\bm{x})$ and their gradients. The energy can then be cast into the form:

[TABLE]

up to local terms containing higher gradients of the displacement fields 111The term $-\frac{3}{8}\frac{b_{3}}{a}P^{2}\sum_{\bm{q}}q^{2}\left|h_{\bm{q}}\right|^{2}$ in Eqs. (21) and (18) is equivalent to Eq. (16) of Ref. Miller and Scriven (1974) in the case $\epsilon_{a}=\epsilon_{b}$ for the hexagonal lattice geometry considered in Ref. Miller and Scriven (1974). This can be seen by taking advantage of the identity

$\sum_{l,m}\frac{m^{2}}{\left(l^{2}+m^{2}-lm\right)^{\frac{5}{2}}}=\frac{2}{3}\sum_{l,m}\frac{1}{\left(l^{2}+m^{2}-lm\right)^{\frac{3}{2}}}~{},$

where the sum runs over all integer pairs with $(l,m)\neq(0,0)$. This identity follows from the symmetry of the hexagonal lattice.. $\kappa_{1}$, $\lambda_{1}$ and $\mu_{1}$ are non-universal constants depending on the details of the lattice geometry. Notice that the negative tension terms

[TABLE]

in Eq. (16) have been absorbed into the term linear in the strain tensor in Eq. (21).

While summing together the electrostatic and the elastic energy, we must set, in Eq. (6),

[TABLE]

so that terms linear in displacement fields disappear from the total energy. The result for the total energy is therefore, in the harmonic approximation:

[TABLE]

with values of the elastic coefficients and bending rigidity which are shifted by electrostatic interactions. The most striking feature of Eq. (23) is that the wavevector dependence of the energy of out-of-plane fluctuations remains proportional to $q^{4}$ and is therefore, as soft as in the purely elastic case. The cancellation from the Hamiltonian of the tension term Toner (1989), proportional to

[TABLE]

is a consequence of rotational invariance in the ambient space. However, a term proportional to

[TABLE]

is not a priori excluded by rotational symmetry. This is evident from Eq. (11). Both $H_{\text{lr},1}$ and $H_{\text{lr},2}$ are individually rotationally invariant interactions and they both lead to a cubic wavevector dependence of the bending energy. This dependence cancels out only from the sum $H_{\text{lr},1}+H_{\text{lr},2}$.

The cancellation of the non-local term $\mathcal{O}_{h}$ from the energy of bending undulations could also be predicted by using, instead of Eq. (11), the equivalent representation of the energy:

[TABLE]

where

[TABLE]

is the Fourier transform in three-dimensional space of the charge density. This method is similar to the approach in Ref. Peliti and Prost (1989). There, however, approximations were made which lead to the prediction of a $\left|\bm{q}\right|^{3}$ dependence of the quadratic part of the bending energy, in contrast with our results. Eq. (25) differs from Eq. (11) only because the former includes the infinite energy of self-interaction of dipole moments. This energy, however, is independent on the configuration of the membrane and constitutes a shift of the Hamiltonian by a global constant. Non-local terms in the expansion of Eq. (25) up to second order in $h$ and $u$ can only arise from the terms:

[TABLE]

where $\rho_{\bm{Q}}^{(1,0)}$ and $\rho_{\bm{Q}}^{(0,1)}$ are the terms linear in $h$ and $u$ in the expansion of $\rho_{\bm{Q}}$, respectively. Denoting the in-plane and the out-of-plane components of $\bm{Q}$ as $\bm{q}$ and $q_{z}$, we find for $\bm{q}$ well within the first Brillouin zone:

[TABLE]

Plugging Eq. (28) into Eq. (27), yields:

[TABLE]

The fluctuation energy can then be obtained by performing the integral over $q_{z}$ explicitly. In the part quadratic in $h$, $\bar{H}^{(2,0)}$, the Coulomb interaction factor $1/Q^{2}$ cancels against the factor $Q^{4}$. The resulting interaction is local. The divergence of the $q_{z}$ integral in Eq. (29) is not problematic, because in order to get the full expression for the energy at order $h^{2}$ and $u^{2}$ we still need to add the purely local terms:

[TABLE]

where $\rho_{\bm{Q}}^{(0,0)}$, $\rho_{\bm{Q}}^{(2,0)}$ and $\rho_{\bm{Q}}^{(0,2)}$ are the terms of zeroth and second order in the expansion of $\rho_{\bm{Q}}$. Since the energy is finite order by order in the expansion in powers of $u$ and $h$, divergences in the $q_{z}$ integral must cancel out. The only non-local term in Eq. (29), for which the Coulomb interaction factor $1/Q^{2}$ does not cancel out, is

[TABLE]

Notice that $\bar{\mathcal{O}}_{u}$ coincides with the non-local term in Eq. (23).

The physical interpretation of the cancellation of non-local terms in the quadratic part of the $h$-field energy is more transparent if the previous argument is translated from Fourier to real space. Consider a localized out-of-plane deformation $h\left(\bm{x}\right)$, which is nonzero only within a small region of the membrane. The electrostatic potential at a point $\bm{R}$ in three-dimensional space far from the deformed region is given, in a continuum approximation, by:

[TABLE]

The change in electrostatic potential due to the fluctuation to first order in $h$ can be calculated by making use of the expansion of the normal in Eq. (12). The result is:

[TABLE]

where $X_{\alpha}$ are the in-plane components of $\bm{R}$. Since the integrand is a total divergence and the deformation $h\left(\bm{x}\right)$ is assumed to be localized, $\delta V^{(1,0)}$ vanishes. If we now consider a second localized wavepacket, far from the first, we can conclude that there will not be a long-range interaction between the two, at least at the lowest order in the amplitude of the $h$ field.

The cancellation of the non-local term $\mathcal{O}_{h}$, Eq. (24), from the Hamiltonian has crucial physical consequences. If it did not vanish, it would completely modify the critical behaviour, as predicted in Ref. Toner (1989).

II.2 Electrostatic energy: expansion to higher orders

We now discuss the properties of higher-order terms in the expansion of the electrostatic energy. For long-wavelength, slowly-varying fluctuations we can calculate the energy approximately by replacing, in Eq. (11):

[TABLE]

This replacement gives the first order in an expansion in the number of gradients, which translates, in Fourier space, in an expansion in powers of the momentum. The derivative $\partial_{\alpha}\bm{n}$ can be neglected, because it is of higher order. Because the interaction is long-ranged, the gradient series is divergent and presents infinite expansion coefficients. The first order, however, is finite. It gives, therefore, the leading behaviour at long wavelengths of interaction vertices.

When Eq. (33) is substituted into Eq. (11), we see immediately that $H_{\text{lr},2}$ vanishes in the first order of the gradient expansion. This follows from the fact that $\partial_{\alpha}\bm{r}(\bm{x}_{i})\cdot\bm{n}_{i}=0$. Therefore:

[TABLE]

Expanding in powers of the strain tensor and summing over $j$ leads to a finite local functional of $u_{\alpha\beta}$, which represents the leading part at small wavevectors of the electrostatic energy. Therefore, as anticipated in Section II.1, at long-wavelengths, dipolar interactions lead to a shift of anharmonic coupling constants of the short-range Hamiltonian. Non-local behaviour appears as a correction, which is suppressed by a power of the wavevector scale. The first terms of the expansion of Eq. (34) are:

[TABLE]

if the crystal structure is sufficiently symmetric to ensure the isotropy of elastic coefficients. Eq. (35) reproduces the results obtained in Sec. II.1 for the electrostatic tension $\tau$ and for the shifts of the Lamé coefficients $\lambda_{1}$ and $\mu_{1}$, as can be seen by comparison with Eqs. (18) and (21).

We now discuss the leading non-local corrections up to order $h^{2}u$ and $h^{4}$. We show that the correction $\bar{\mathcal{O}}_{u}$ in the harmonic part of the energy is promoted to:

[TABLE]

a rotationally invariant interaction, which can be interpreted as a $q$-dependent correction to the Lamé coefficient $\lambda$. Up to order $uh^{2}$ and $h^{4}$, $\mathcal{O}_{u}$ is the only non-local interaction to be generated.

In order to calculate the non-local part of interaction vertices, it is convenient to expand the electrostatic energy starting from the expression in Eq. (25). Using the expansion for the normal to second order in $h$ and zeroth order in $u$:

[TABLE]

we can calculate $\rho_{\bm{Q}}^{(2,0)}$, the term of order $h^{2}$ in the expansion of the Fourier transform of the charge density $\rho_{\bm{Q}}$:

[TABLE]

As in Eq. (28), $\bm{q}$ and $q_{z}$ denote the in-plane and the out-of plane components of the three-dimensional wavevector $\bm{Q}$. Again Eq. (38) holds under the assumption that $\bm{q}$ lies well within the first Brillouin zone of the lattice. Let us now consider the expansion of the energy, Eq. (25), at order $uh^{2}$:

[TABLE]

The non-local behaviour of the vertex is only generated by the first term:

[TABLE]

The third term in Eq. (39) is clearly local because fields are inserted at only one of the two $\rho_{\bm{Q}}$. In the second term, the Coulomb interaction $1/Q^{2}$ cancels against the factor $Q^{2}$ in $\rho_{\bm{Q}}^{(1,0)}$, see Eq. (28). The latter cancellation has the same origin of the cancellation of the propagator correction $\mathcal{O}_{h}$. Similarly, it can be seen that the non-local behaviour of the quartic vertex $H^{(4,0)}$ is completely encoded in the term:

[TABLE]

Eqs. (40) and (41) severely restrict the possible form of the long-range, non-analytic part of interaction vertices. The momentum transfer $\bm{q}$, the only on which the energy can have a singular dependence, can not be exchanged by a single $h$ field, but only by pairs. In real space this can again be visualized by recalling that the potential generated by a localized wavepacket of out-of-plane deformation vanishes to first order in its amplitude.

Combining Eqs. (28), (38), (40) and (41) and keeping only the non-local part of vertices, the part in which the factor $1/Q^{2}$ does not cancel out from the calculations, we obtain:

[TABLE]

In Eq. (42) all equivalencies are intended up to local functionals of the fields. We can now obtain the complete expression for the long-wavelength Hamiltonian by combining Eqs. (35), (42) and Eqs. (6), (22). The result is:

[TABLE]

This result extends Eq. (23), which was found within the harmonic approximation.

The non-local term $\mathcal{O}_{u}$ is suppressed by the factor $\left|\bm{q}\right|$ in the long-wavelength limit. It constitutes, therefore, an irrelevant perturbation by power counting in a neighbourhood of $D=4$, which is the upper critical dimension for crystalline membranes Aronovitz and Lubensky (1988). Non-local terms of higher order in powers of the fields $u$ and $h$ have a structure which is similar to $\mathcal{O}_{u}$. They are suppressed by a power $\left|\bm{q}\right|$ with respect to local couplings and they also constitute irrelevant operators by power counting.

The fact that $\mathcal{O}_{u}$ is irrelevant even for physical 2D membranes is suggested by the following considerations. In a scaling analysis, we expect that $\mathcal{O}_{u}$ is irrelevant if $\left|\lambda_{\rm R}(q)\right|,\left|\mu_{\rm R}(q)\right|\gg 2\pi P^{2}\left|\bm{q}\right|$ at small $q$. This comparison relies on the assumption that $\mathcal{O}_{u}$ is not renormalized, which is suggested by its non-locality. For $D$-dimensional membranes embedded in a $d$-dimensional space, the renormalized elatic moduli $\lambda_{\rm R}\left(q\right)$ and $\mu_{\rm R}\left(q\right)$ vanish as $q^{4-D-2\eta}$ for $q\rightarrow 0$ Nelson et al. (2004); Aronovitz and Lubensky (1988); Bowick and Travesset (2001). The condition $\left|\lambda_{\rm R}\left(q\right)\right|,\left|\mu_{\rm R}\left(q\right)\right|\gg 2\pi P^{2}\left|\bm{q}\right|$, then, is fulfilled provided that $\eta>(3-D)/2$. Field-theoretical calculations Le Doussal and Radzihovsky (1992); Gazit (2009); Le Doussal and Radzihovsky (2018); Kownacki and Mouhanna (2009); Guitter et al. (1989) indicate that $\eta>1/2$ in the case of clean two-dimensional membranes, thus giving an indication that $\mathcal{O}_{u}$ is irrelevant in the physical case $D=2$, $d=3$. However, scaling with an exponent $\eta<1/2$ was predicted to occur in 2D disordered membranes Le Doussal and Radzihovsky (2018); Gornyi et al. (2015). In this case our considerations open a question about the possible role of non-local terms.

If, on the other hand, the statistical mechanical properties of the membrane are addressed within the self-consistent screening approximation (SCSA) Le Doussal and Radzihovsky (1992, 2018), $\mathcal{O}_{u}$ appears to be irrelevant irrespective of the value of $\eta$. This follows immediately from the fact that $\mathcal{O}_{u}$ is equivalent to a $\bm{q}$-dependent shift of the first Lamé coefficient. Within the SCSA, the local Lamé coefficient $\lambda^{\prime}$ and its non-local correction $-2\pi P^{2}\left|\bm{q}\right|$ enter the self-consistent equations only through their sum. For sufficiently small $\bm{q}$, the non-local correction is negligible. The SCSA equations in the long-wavelength limit Le Doussal and Radzihovsky (1992, 2018), therefore, remain unchanged and the exponent $\eta$ remains necessarily as in the pure short-range case.

II.3 Shift of elastic coefficients

We now discuss the shift of elastic coefficients induced by electrostatic interactions. The relations between the shifted elastic constants $\lambda^{\prime}$ and $\mu^{\prime}$ and the Lamé coefficients $\lambda$ and $\mu$ which appear in the short range Hamiltonian $H_{\text{sr}}$, Eq. (3) can be read from Eq. (18), or equivalently from Eq. (35):

[TABLE]

The shift in the bending rigidity $\kappa^{\prime}-\kappa$ could be determined from Eq. (16) and  (18) by expanding $f_{2}\left(a\bm{q}\right)$ and $f_{1}\left(a\bm{q}\right)$ up to fourth order in $\bm{q}$. This would require an explicit calculation of $q$-dependent lattice sum, which could be performed with the Ewald method (see the Appendix) with an explicit choice of the Ewald cutoff function $\varphi\left(x\right)$ 222For widely used Ewald cutoff functions, see for example Aharony and Fisher (1973); Grechnev et al. (2005); Cohen and Keffer (1955). In this work, we did not address the renormalization of the bending rigidity explicitly.

The renormalization of in-plane elastic coefficients, instead, is determined by the simple expression in Eq. (44). We notice, however, that the ’bare’ coefficients $\kappa$, $\lambda$ and $\mu$ do not coincide exactly with the elastic coefficients $\kappa_{0}$, $\lambda_{0}$ and $\mu_{0}$ of the membrane when dipole interactions are switched off. The reason is that dipole-dipole repulsion induces a finite stretching of the membrane which results in an increase of the elastic coefficients. The elastic constants $\kappa$, $\lambda$ and $\mu$, being defined through Eq. (3), represent elastic coefficients after stretching has already occurred.

If the effects of dipole repulsion are weak and the stretching is small so that linear elasticity theory applies while electrostatic interactions are turned on, we can relate the two sets of coefficients through:

[TABLE]

where $\xi$ is the stretching factor, determined by:

[TABLE]

Here $a$ and $P$ are the ’final’ values of the lattice constant and the polarization density, calculated at the equilibrium value of the stretching factor in presence of electrostatic interactions, the same $a$ and $P$ which were used throughout the rest of this paper. By the previous considerations, Eq. (46) holds under the assumption $\xi-1\ll 1$.

We now give an order of magnitude estimate of the effects of dipolar interactions on the elastic parameters of the membrane. We consider, as an example, the case of graphene derivatives involving Group IA and Group VIIA atoms Klintenberg et al. (2010). In these materials alkali and halogen atoms form partially ionic bonds with carbon atoms in the graphene layer. If alkali atoms bind on one side of the graphene layer and halogen atoms on the opposite side, an out-of-plane permanent polarization results. This occurs, for example, in C2HF and C2LiF in the ’chair’ geometry described in Ref. Klintenberg et al. (2010). For these materials, we expect that the typical scale of the polarization density $P$ is roughly 1 e/Å. From Eq. (44), it follows that corrections to the Lamé coefficients are of the order of 1 eV/Å2. They are, therefore, of the same order of magnitude of the ’bare’ elastic coefficients $\lambda$ and $\mu$. The same considerations suggest that the shift of the bending rigidity is of the order of 1 eV.

More specifically, we expect that among the class of graphene derivatives involving Group IA and Group VIIA atoms, the largest dipole moments and the strongest polarization effects arise in C2LiF Klintenberg et al. (2010). An estimate of the shift of elastic coefficient can be done, within our model, by describing C2LiF (in the ’chair’ geometry) as a lattice with dipole moments $q_{\text{Li}}d_{\text{Li}}$ and $-q_{\text{F}}d_{\text{F}}$ located at the two sublattices of the hexagonal graphene structure. For the charge transfer to lithium and fluorine atoms, we assume in an order of magnitude estimate, $q_{\text{Li}}=+\text{e}$ and $q_{\text{F}}=-\text{e}$. The structural properties of C2LiF were calculated in Ref. Klintenberg et al. (2010). The bond lengths $d_{\text{Li}}$ and $d_{\text{F}}$ and the lattice constant $a$ are 2.16 Å, 1.44 Å and 2.59 Å, respectively. We obtain a polarization density $P\simeq$ 0.6 e/Å. The value of the lattice sum $b_{3}$ is close to $9$ Grechnev et al. (2005) in the case of a square lattice and we expect that it remains of the same order in a hexagonal lattice. Eq. (44) then leads to corrections to the Lamé coefficients of the order of 10 eV/Å2. This result is probably an overestimation. It suggests, however, that electrostatic effects can play an important role in the elastic properties of C2LiF.

III Conclusions

In conclusion, we addressed the effects of dipole-dipole interactions on the long-wavelength thermal fluctuations of membranes exhibiting a permanent out-of-plane polarization directed at each point along the local normal. We focused on crystalline membranes in the flat phase, a case which is particularly important for the thermodynamic properties of two-dimensional materials at room temperature.

We found that, even in presence of long-range dipole-dipole interactions, the dispersion relation of bending fluctuations remains as soft as in the short-range case and can be described at long wavelengths in terms of a bending rigidity. Furthermore, the leading behaviour of the electrostatic energy in the limit of long wavelengths can be represented through a local functional of the strain tensor. These results have important consequences in the statistical mechanical properties of the membrane. They imply that the effective Hamiltonian, in the leading order in the limit of long wavelengths, has the same local form of the Hamiltonian of elasticity theory. Non-local terms are suppressed in the long-wavelength limit and they are irrelevant by power counting. The expected result, therefore, is that the large-distance behaviour and the scaling exponents of elastic membranes are not modified by dipole-dipole interactions, at least in a neighbourhood of $D=4$. A scaling analysis and the self-consistent screening approximation support the expectation that the non-local correction of smallest order in powers of the displacement fields remains irrelevant even for $D=2$, at least in the case of clean membranes, for which the exponent $\eta$ is larger than $1/2$. The condition $\eta>1/2$, however, can be violated in disordered membranes Le Doussal and Radzihovsky (2018); Gornyi et al. (2015). In this case, the role of dipole-dipole interactions in the scaling behaviour raises a new question.

We believe that the model which we studied in this work can help to understand thermal fluctuations of two-dimensional materials, such as the graphene derivatives CF, C2HF and C2LiF Klintenberg et al. (2010), which exhibit permanent out-of-plane polarization. More and more 2D materials with finite out-of-plane polarization have been predicted or discovered in recent years, including, for example, MoSSe Riis-Jensen et al. (2018) and WTe2 Fei et al. (2018).

In addition, we addressed the shift of elastic coefficients due to long-range interactions. As a model for crystalline membranes with out-of-plane polarization, we considered a deformable elastic sheet with an embedded lattice of point dipole moments. We derived explicit expressions for the shift of the Lamé coefficients $\lambda$ and $\mu$ induced by electrostatics and an implicit expression for the renormalization of the bending rigidity. We estimated that the renormalization of elastic moduli can be important in graphene derivatives such as C2LiF Klintenberg et al. (2010), in which polarization arises from charge transfer between carbon, alkali and halogen atoms which form partially ionic bonds. Analytic expressions for the shift of elastic coefficients in our model could be useful, for example, as a basis for comparison with computational results for the structural properties of 2D materials.

Acknowledgements.

This work was supported by the Netherlands Organisation for Scientific Research (NWO) via the Spinoza Prize.

Appendix A Harmonic part of the electrostatic energy: details of the calculations

All calculations in Section II.1 rely on Fourier transforms of power-law interactions on a lattice. It is convenient to begin by studying the sum:

[TABLE]

for $n=3$. As in the main text, $A$ is the total area, $P$ the polarization density and $a$ is the lattice constant. The interesting, long-range, character of power-law interactions results in a non-analytic behaviour of $S_{n}\left(\bm{q}\right)$ as $\bm{q}\rightarrow 0$. The non-analyticity arises from the large distance part of the sum in Eq. 47, where the summation can be accurately approximated with a continuum integral. Replacing the lattice sum with an integral is however impossible at small $\left|\bm{x}_{i}-\bm{x}_{j}\right|$, because of the short-distance divergence. It is convenient, therefore to make use of an Ewald method (See e.g. Aharony and Fisher, 1973; Grechnev et al., 2005; Cohen and Keffer, 1955). We decompose $S_{3}\left(\bm{q}\right)$ into a short-range part $S_{3}^{<}\left(\bm{q}\right)$ and a long-range part $S_{3}^{>}\left(\bm{q}\right)$:

[TABLE]

Here $\varphi$ is a cutoff function, which vanishes faster than any power of $\left|\bm{x}\right|$ for $\left|\bm{x}\right|\rightarrow\infty$ and such that $1-\varphi\left(\left|\bm{x}\right|/a\right)$ vanishes faster than $\left|\bm{x}\right|$ at short distances. The particular choice of $\varphi$ is unimportant in the following discussion. $S_{3}^{<}\left(\bm{q}\right)$, due to the large-distance cutoff, is a regular function of the components of its vector argument. The non-analytic wavevector dependence arises only from $S_{3}^{>}\left(\bm{q}\right)$. In order to calculate $S_{3}^{>}\left(\bm{q}\right)$, it is convenient to transform the sum over the direct lattice into a sum over reciprocal lattice vectors, by making use of the Poisson summation formula. For a generic function, this transformation reads Aharony and Fisher (1973); Grechnev et al. (2005); Cohen and Keffer (1955):

[TABLE]

Here $\Omega$ is the area of the unit cell of the Bravais lattice, $\bm{x}_{\ell}$ are primitive vectors of the direct Bravais lattice and $\bm{g}_{m}$ are reciprocal lattice vectors. In the second and the third line of Eq. (49), the sum over $i$ and $j$ is restricted to the dipoles lying within the first unit cell, denoted as $\Omega_{0}$. The positions of dipoles within the cell $\Omega_{0}$ are the generating vectors $\bm{t}_{i}$. Eq. (49) can be directly applied to the evaluation of $S_{3}^{>}\left(\bm{q}\right)$ by choosing:

[TABLE]

The transformation gives:

[TABLE]

In Eq. (51), all terms with a finite $\bm{g}_{m}$ are regular functions of $\bm{q}$ for $\bm{q}\rightarrow 0$. Singular behaviour can only occur in the $\bm{g}_{m}=0$ term of the sum Aharony and Fisher (1973). This term coincides precisely with the continuum integral approximation to $S_{3}^{>}\left(\bm{q}\right)$. Therefore:

[TABLE]

where $r\left(a\bm{q}\right)$ is a regular function.

In order to calculate the integral in Eq. (52), it is useful to make use of the Fourier transform:

[TABLE]

The integral can be written as:

[TABLE]

The second term in the left-hand side of Eq. (54) is a regular function of $\bm{q}$, due to the large-distance cutoff. In addition, the limit $\varepsilon\rightarrow 0$ must be finite, because the initial expression in Eq. (52), in which $\varepsilon=0$ is well defined. As a result:

[TABLE]

where $r^{\prime}\left(a\bm{q}\right)$ is a new regular function of the wavevector. Summing $S_{3}^{<}(\bm{q})$ and $S_{3}^{>}(\bm{q})$ finally gives:

[TABLE]

$s_{3}$ is regular and vanishes as $q^{2}$ at small $\bm{q}$. The constant term $b_{3}=S_{3}(0)$ coincides with the coefficient defined in Eq. (19) in the main text.

The calculations in Section II.1 require, besides $S_{3}\left(\bm{q}\right)$, also the lattice sums $S_{5}\left(\bm{q}\right)$ and $S_{7}\left(\bm{q}\right)$. These are related to $S_{3}(\bm{q})$ through:

[TABLE]

If the functions $S_{n}\left(\bm{q}\right)$ were perfectly isotropic, Eq. (57) would be sufficient to determine the coefficient of $\left|\bm{q}\right|^{m+2}$ in $S_{n+2}\left(\bm{q}\right)$ from the coefficient of $\left|\bm{q}\right|^{m}$ in $S_{n}\left(\bm{q}\right)$. In general, $S_{n}\left(\bm{q}\right)$ are anisotropic, although they exhibit the same symmetry of the crystal. However, the non-analytic part of any $S_{n}\left(\bm{q}\right)$ must be perfectly isotropic. The steps which led to the calculation of $S_{3}\left(\bm{q}\right)$ could be repeated similarly for arbitrary $n$. The singular dependence on $\bm{q}$ arises only from a term analogue to Eq. (52), a continuum isotropic integral with no traces of the crystal structure anisotropy. For the non-singular part, we will make the assumption that the symmetry of the crystal is large enough to forbid anisotropic functions at least up to fourth order in $\bm{q}$. Under this assumption, using Eq. (57) recursively yields:

[TABLE]

where $s_{5}$ and $s_{7}$ are regular and vanish as $q^{4}$ and $q^{6}$ respectively. Now all necessary ingredients for the calculation of $H_{\text{lr},1}^{(2)}$ and $H_{\text{lr},2}^{(2,0)}$ have been derived. Fourier transformation of $H_{\text{lr},2}^{(2,0)}$ gives:

[TABLE]

where $f_{2}(a\bm{q})$ is an analytic function at $\bm{q}=0$, which vanishes as $q^{4}$ for $\bm{q}\rightarrow 0$. The Fourier transform of $H_{\rm lr,1}$ is calculated similarly. From Eq. (17) in the main text:

[TABLE]

Using Eq. (56) and Eq. (58) leads to the result in the main text, Eq. (18).