Two-dimensional melting via sine-Gordon duality

TL;DR

This paper presents a dual sine-Gordon model to describe two-dimensional crystal melting, capturing the two-stage transition via renormalization-group analysis and connecting it to Coulomb gas results.

Contribution

It introduces a classical dual sine-Gordon framework for 2D melting, linking elasticity, dislocation dynamics, and phase transitions in a unified model.

Findings

Reproduces Coulomb gas flow equations for elastic and defect parameters.

Describes continuous two-stage melting via relevance of cosine operators.

Provides a transparent RG analysis of dislocation and disclination unbinding.

Abstract

Motivated by the recently developed duality between elasticity of a crystal and a symmetric tensor gauge theory by Pretko and Radzihovsky, we explore its classical analog, that is a dual theory of the dislocation-mediated melting of a two-dimensional crystal, formulated in terms of a higher derivative vector sine-Gordon model. It provides a transparent description of the continuous two-stage melting in terms of the renormalization-group relevance of two cosine operators that control the sequential unbinding of dislocations and disclinations, respectively corresponding to the crystal-to-hexatic and hexatic-to-isotropic fluid transitions. This renormalization-group analysis compactly reproduces seminal results of the Coulomb gas description, such as the flows of the elastic couplings and of the dislocation and disclination fugacities, as well the temperature dependence of the associated correlation lengths.

Figures

Equations

$$\begin{array}{rl}\tilde{H}& =\int d^{2}r[\frac{1}{2}{\tilde{C}}_{ij,kl}^{-1}{\partial }_i{\partial }_j\varphi {\partial }_k{\partial }_l\varphi \\ & -g_b\sum _{n=1}^p\cos ({\mathbf{b}}_n\cdot \hat{\mathbf{z}}\times \nabla \varphi )-g_s\cos (s_p\varphi )].\end{array}$$

$$\tilde{H}\sim {\int }_{\mathbf{r}}\left[\frac{1}{2}\tilde{C}({\partial }^2\varphi {)}^2-g_b\cos (\partial \varphi )-g_s\cos (\varphi )\right],$$

$${\tilde{H}}_{\text{cr}}\sim {\int }_{\mathbf{r}}\left[\frac{1}{2}\tilde{C}(\partial \chi {)}^2-g_b\cos (\chi )\right],$$

$${\tilde{H}}_{\text{hex}}\sim {\int }_{\mathbf{r}}\left[\frac{1}{2}{g}_b(\partial \varphi {)}^2-g_s\cos (\varphi )\right],$$

$$\mathcal{H}=\frac{1}{2}{C}_{ij,kl}{u}_{ij}{u}_{kl},$$

$$u_{ij}\approx \frac{1}{2}({\partial }_i{u}_j+{\partial }_j{u}_i).$$

$$C_{ij,kl}=\lambda {\delta }_{ij}{\delta }_{kl}+2\mu {\delta }_{ik}{\delta }_{jl}$$

$${\oint }_{{\mathbf{r}}_0}d\theta =\frac{2\pi }{p}s$$

$$\hat{\mathbf{z}}\cdot \nabla \times \nabla \theta =\frac{2\pi }{p}s{\delta }^2(\mathbf{r}-{\mathbf{r}}_0)\equiv \frac{2\pi }{p}s(\mathbf{r}),$$

$${\oint }_{{\mathbf{r}}_0}d\mathbf{u}={\mathbf{b}}_n,$$

$$\hat{\mathbf{z}}\cdot \nabla \times \nabla u_i=b_{i,n}{\delta }^2(\mathbf{r}-{\mathbf{r}}_0)\equiv b_i(\mathbf{r}),$$

$$u_i$$

$$u_{ij}$$

$$\begin{array}{rl}Z& =\int [d\mathbf{u}]e^{-{\int }_{\mathbf{r}}\mathcal{H}[\mathbf{u}]}\\ & =\int [d\tilde{\mathbf{u}}][d{\mathbf{u}}^s]\int [d{\sigma }_{ij}(\mathbf{r})]e^{-{\int }_{\mathbf{r}}\mathcal{H}[\mathbf{u},{\sigma }_{ij}]},\end{array}$$

$$\begin{array}{rl}\mathcal{H}[\mathbf{u},{\sigma }_{ij}]& =\frac{1}{2}{C}_{ij,kl}^{-1}{\sigma }_{ij}{\sigma }_{kl}+i{\sigma }_{ij}{u}_{ij}\\ & =\frac{1}{2}{C}_{ij,kl}^{-1}{\sigma }_{ij}{\sigma }_{kl}+i{\sigma }_{ij}\left({\partial }_i{\tilde{u}}_j+u_{ij}^s\right).\end{array}$$

$$C_{ij,kl}^{-1}=-\frac{\lambda }{4\mu (\mu +\lambda )}{\delta }_{ij}{\delta }_{kl}+\frac{1}{2\mu }{\delta }_{ik}{\delta }_{jl}$$

$${\partial }_i{\sigma }_{ij}=0,$$

$${\sigma }_{ij}={ϵ}_{ik}{ϵ}_{jl}{\partial }_k{\partial }_l\varphi .$$

$${ϵ}_{ik}{ϵ}_{jl}{\partial }_l{\partial }_k{u}_{ij}^s$$

$$=$$

$$\mathcal{H}[\varphi ]$$

$$H_{\mathbf{b}}=\frac{1}{2}\int \frac{d^{2}q}{(2\pi {)}^2}{b}_i(\mathbf{q}){\tilde{K}}_{ij}(\mathbf{q})b_j(-\mathbf{q}),$$

$${\tilde{K}}_{ij}(\mathbf{q})$$

$$K_{ij}(\mathbf{r})$$

$$\begin{array}{rl}{H}_{\mathbf{b}}& =-\frac{K}{8\pi }{\int }_{{\mathbf{r}}_1,{\mathbf{r}}_2}[\mathbf{b}({\mathbf{r}}_1)\cdot \mathbf{b}({\mathbf{r}}_2)\ln \frac{|{\mathbf{r}}_1-{\mathbf{r}}_2|}{a}\\ & -\frac{\mathbf{b}({\mathbf{r}}_1)\cdot ({\mathbf{r}}_1-{\mathbf{r}}_2)({\mathbf{r}}_1-{\mathbf{r}}_2)\cdot \mathbf{b}({\mathbf{r}}_2)}{|{\mathbf{r}}_1-{\mathbf{r}}_2{|}^2}],\end{array}$$

$$\mathbf{b}(\mathbf{r})$$

$$s(\mathbf{r})$$

$$H=$$

$$+\sum _{{\mathbf{r}}_n}\left[i\frac{2\pi }{6}\varphi ({\mathbf{r}}_n)s_{{\mathbf{r}}_n}-i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_n)\cdot {\mathbf{b}}_{{\mathbf{r}}_n}\right],$$

$$\begin{array}{rl}Z=& \int [d\varphi ]\sum _{\{s_{{\mathbf{r}}_n}\}}\sum _{\{{\mathbf{b}}_{{\mathbf{r}}_n}\}}\prod _{{\mathbf{r}}_n}{e}^{-H[\varphi ,{\mathbf{b}}_{{\mathbf{r}}_n},s_{{\mathbf{r}}_n}]},\\ =& \int [d\varphi ]e^{-\frac{1}{2}{\int }_{\mathbf{r}}{\tilde{C}}_{ij,kl}^{-1}{\partial }_i{\partial }_j\varphi {\partial }_k{\partial }_l\varphi }\sum _{\{s_{{\mathbf{r}}_n}\}}\sum _{\{{\mathbf{b}}_{{\mathbf{r}}_n}\}}\prod _{{\mathbf{r}}_n}\left[e^{-i\frac{2\pi }{6}\varphi ({\mathbf{r}}_n)s_{{\mathbf{r}}_n}-E_s{s}_{{\mathbf{r}}_n}^2}{e}^{i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_n)\cdot {\mathbf{b}}_{{\mathbf{r}}_n}-{\tilde{E}}_b{b}_{{\mathbf{r}}_n}^2}\right],\\ =& \int [d\varphi ]e^{-\frac{1}{2}{\int }_{\mathbf{r}}{\tilde{C}}_{ij,kl}^{-1}{\partial }_i{\partial }_j\varphi {\partial }_k{\partial }_l\varphi }\left[1+e^{-2E_s}\int \frac{d^2{r}_1}{a^2}\frac{d^2{r}_2}{a^2}{e}^{-i\frac{2\pi }{6}\varphi ({\mathbf{r}}_1)}{e}^{i\frac{2\pi }{6}\varphi ({\mathbf{r}}_2)}+\dots \right]\\ & \left[1+e^{-2E_b}\sum _{n=1}^3\int \frac{d^2{r}_1}{a^2}\frac{d^2{r}_2}{a^2}{e}^{-i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_1)\cdot {\mathbf{b}}_n}{e}^{i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_2)\cdot {\mathbf{b}}_n}+\dots \right],\\ =& \int [d\varphi ]e^{-\frac{1}{2}{\int }_{\mathbf{r}}{\tilde{C}}_{ij,kl}^{-1}{\partial }_i{\partial }_j\varphi {\partial }_k{\partial }_l\varphi }\cdot [1+e^{-E_s}\int \frac{d^2{r}_1}{a^2}\left(e^{i\frac{2\pi }{6}\varphi ({\mathbf{r}}_1)}+e^{-i\frac{2\pi }{6}\varphi ({\mathbf{r}}_1)}\right)\\ & +\frac{1}{2!}{e}^{-2E_s}\int \frac{d^2{r}_1}{a^2}\frac{d^2{r}_2}{a^2}\left(e^{i\frac{2\pi }{6}\varphi ({\mathbf{r}}_1)}+e^{-i\frac{2\pi }{6}\varphi ({\mathbf{r}}_1)}\right)\left(e^{i\frac{2\pi }{6}\varphi ({\mathbf{r}}_2)}+e^{-i\frac{2\pi }{6}\varphi ({\mathbf{r}}_2)}\right)+\dots ]\\ & \prod _{n=1,2,3}[1+e^{-E_b}\int \frac{d^2{r}_1}{a^2}\left(e^{i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_1)\cdot {\mathbf{b}}_n}+e^{-i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_1)\cdot {\mathbf{b}}_n}\right)\\ & +\frac{1}{2!}{e}^{-2E_b}\int \frac{d^2{r}_1}{a^2}\frac{d^2{r}_2}{a^2}\left(e^{i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_1)\cdot {\mathbf{b}}_n}+e^{-i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_1)\cdot {\mathbf{b}}_n}\right)\left(e^{i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_2)\cdot {\mathbf{b}}_n}+e^{-i\hat{\mathbf{z}}\times \nabla \varphi ({\mathbf{r}}_2)\cdot {\mathbf{b}}_n}\right)+\dots ],\\ \equiv & \int [d\varphi ]e^{-\tilde{H}},\end{array}$$

$$\tilde{H}=$$

$$-g_b\sum _{n=1}^3\cos ({\mathbf{b}}_n\cdot \hat{\mathbf{z}}\times \nabla \varphi )-g_s\cos \left(\frac{2\pi }{6}\varphi \right)].$$

$$\begin{array}{rl}\mathcal{H}& =\frac{1}{2}{C}_{ij,kl}({\partial }_i{u}_j-\theta {ϵ}_{ij})({\partial }_k{u}_l-\theta {ϵ}_{kl})+\frac{1}{2}K({\partial }_i\theta {)}^2.\end{array}$$

Full text

I Introduction

I.1 Background and motivation

The theory of continuous two-dimensional (2D) melting, developed by Kosterlitz and ThoulessKosterlitz 72 , Halperin and NelsonNelson 79 , and YoungYoung 79 (KTHNY), building on the work of LandauLandau 37 , PeierlsPeierls 35 and BerezinskiiBerezinskii 71 ; Berezinskii 72 , has become one of the pillars of theoretical physics. Mathematically related to simpler normal-to-superfluid and planar paramagnet-to-ferromagnet transitions in films, described by a 2D XY model, it is a striking example of the increased importance of thermal fluctuations in low-dimensional systemsHalperin 79 ; Nelson 83 . In contrast to their bulk three-dimensional analogs, where, typically, fluctuations only lead to quantitative modifications of mean-field predictions (e.g., values of critical exponents), here the effects are qualitative and drastic. Located exactly at the lower-critical dimension, a local-order-parameter distinction between the high- and low-temperature phases that is erased by fluctuations, two-dimensional melting can proceed via a subtle, two-stage, continuous transition, driven by the unbinding of topological dislocations and disclinations defects. This mechanism, made possible by strong thermal fluctuations, thus provides an alternative route to direct first-order melting, argued by Landau’s mean-field analysisLandau 37 to be the exclusive scenario.

As such, the continuous two-dimensional melting (and related disordering of a 2D XY model) is the earliest example of a thermodynamically sharp, topological phase transition between two locally disordered phases, which thus does not admit Landau’s local order-parameter description. It is controlled by a fixed line, which lends itself to an asymptotically exact analysisKosterlitz 72 ; Nelson 79 ; Young 79 .

Although evidence for defects-driven phase transitions has appeared in a number of experiments on liquid crystalsHuang 93 and Langmuir-Blodgett filmsKnobler 07 , finding simple model systems that exhibit these phenomena in experiments or simulations has proven to be more challenging. Most studied systems appear to exhibit discontinuous first-order melting that converts a crystal directly into a liquid. However, it appears that two-stage continuous melting has been experimentally observed by MurrayMurray 88 and ZahnZahn 99 in beautiful melting experiments on two-dimensional colloids confined between smooth glass plates and superparamagnetic colloidal systems, respectively. In these experiments, an orientationally quasi-long-range ordered but translationally disordered hexatic phaseNelson 79 was indeed observed. As was first emphasized by Halperin and NelsonNelson 79 , the hexatic liquid, intermediate but thermodynamically distinct from the 2D crystal and the isotropic liquid, is an important signature of the defect-driven two-stage melting. In these two-dimensional colloids, particle positions and the associated topological defects can be directly imaged via digital video-microscopy, allowing precise quantitative tests of the theory.

I.2 Duality of the two-stage melting transition

The disordering of the simpler 2D XY model (describing, e.g., a superfluid-normal transition in a film) is well known to admit two complementary descriptions, the 2D Coulomb gas of vorticesCoulombGas and its dual sine-Gordon field theorysineGordon ; sineGordonKT ; Nelson 79 ; Halperin 79 ; ChaikinLubensky ; RadzihovskyLubensky . As with other dualities – a subject with long history and of much current interestSenthil 18 – the sine-Gordon duality has been extensively utilized in a variety of physical contexts. Given that elasticity of a crystal can be thought of as a space-spin coupled vector generalization of an XY model (with vector phonon Goldstone modes $u_{x},u_{y}$ replacing the scalar phase angle), it is of interest also to develop an analogous dual sine-Gordon formulation and to use it to study the 2D continuous melting transition.

Indeed, recently, such a complementary description has emerged as a classical limit of the elasticity-to-tensor gauge theory dualityPretkoRL 18/05 ; PretkoRL 18/12 , derived in the context of a new class of topologically-ordered fracton matterNandkishore 18 . As we will detail in the body of the paper the corresponding dual Hamiltonian is given by

[TABLE]

Its key features, which characterize the continuous two-stage melting are the higher-order “Laplacian elasticity”, encoded via elastic constants $\tilde{C}^{-1}_{ij,kl}$, and two sine-Gordon types of operators with couplings $g_{b},g_{s}$, respectively, capturing the importance (fugacities) of dislocation (elementary vector charges ${\bf b}_{n}$) and disclination (elementary scalar charge $s_{p}$) defects.

To flesh out the essence of this dual description, neglecting inessential details, the above Hamiltonian is schematically described by

[TABLE]

where $\int_{\bf r}\equiv\int d^{2}r$. Because of the second-order Laplacian elasticity, standard analysis around the Gaussian fixed line $g_{b}=g_{s}=0$ shows that the mean-squared fluctuations of Airy-stress potential $\phi$ diverge quadratically with system size. This leads to an exponentially (as opposed to power-law in a conventional sine-Gordon model) vanishing Debye-Waller factor, and in turn to a strongly irrelevant disclination cosine, $g_{s}$, that can therefore be neglected, whenever $g_{b}$ is small, i.e., near the Gaussian fixed line.

The schematic Hamiltonian then reduces to

[TABLE]

with $\chi=\partial\phi$. It thus obeys the standard sine-Gordon phenomenology, exhibiting a KT-like “roughenning” transition in $\chi$ with the relevance of $g_{b}$, controlled by the stiffness $\tilde{C}$.sineGordon ; sineGordonKT ; Nelson 79 ; Halperin 79 ; ChaikinLubensky At small $\tilde{C}<\tilde{C}_{c}$, $g_{b}$ is irrelevant, describing the gapless crystal phase, with confined dislocations and disclinations. The melting of the crystal is then captured by the relevance of $g_{b}$ for $\tilde{C}>\tilde{C}_{c}$, corresponding to a transition into a plasma of unbound dislocatons characteristic of a hexatic fluid. Since in this phase $g_{b}$ is relevant, at sufficiently long scales the dislocation cosine in Eq. $(2)$ reduces to a harmonic potential for $\chi$, $-g_{b}\cos(\partial\phi)\sim{\frac{1}{2}}g_{b}(\partial\phi)^{2}$. The effective Hamiltonian is then given by

[TABLE]

where we have neglected the $\tilde{C}$ “curvature” elasticity relative to the gradient one encoded in $g_{b}$, and restored the disclination cosine operator $g_{s}\cos(\phi)$. The resulting conventional sine-Gordon model in $\phi$ can then exhibit the second KT-like “roughenning” transition, capturing the hexatic-to-isotropic fluid transition, asssociated with the unbinding of disclinations. The corresponding RG flow of the dual vector sine-Gordon model is schematically illustrated in Fig. 1. We leave the detailed analysis of this two-stage melting transition to the main body of the paper and the Appendix.

I.3 Outline

The rest of this paper is organized as follows. In Sec. II, after briefly reviewing the elasticity theory of two-dimensional crystal and its topological defects, we give two detailed complementary derivations of the duality transformation to the vector sine-Gordon model. Utilizing the latter we straightforwardly reproduce known results for the crystal-hexatic phase transition in Sec. III, by focusing on the dislocation fugacity cosine operator and neglecting the irrelevant disclinations. Inside the hexatic phase, we derive a scalar sine-Gordon model for the Airy stress potential, that captures the subsequent hexatic-isotropic liquid transition. We conclude in Sec. IV with a summary of our results and a discussion of potential applications of this dual approach.

II Duality of 2D melting

II.1 Two-dimensional elasticity

At low temperatures, the deformations of a crystal do not vary substantially over the atom size, allowing it to be described with a continuum field theory of its phonon Goldstone modes, $u({\bf r})$, with a short-distance cutoff set by the lattice spacing $a$. The underlying translational symmetry (spontaneously broken by the crystal), requires that the elastic energy is an analytic expansion in the strain field $\partial_{i}u_{j}$. Due to the underlying rotational symmetry in the target space (i.e., no substrates and/or external fields), the elastic Hamiltonian is further constrained (in harmonic order) to be independent of the antisymmetric part of $\partial_{i}u_{j}$, i.e., of the local bond angle $\theta({\bf r})=\frac{1}{2}\epsilon_{ij}\partial_{i}u_{j}=\frac{1}{2}\hat{\bf z}\cdot{\bm{\nabla}}\times{\bf u}$. The elastic Hamiltonian density to harmonic order is thus given by

[TABLE]

where $u_{ij}=\frac{1}{2}(\partial_{i}{\bf r}\cdot\partial_{j}{\bf r}-\delta_{ij})$ is the symmetric nonlinear strain tensor (fully rotationally invariant in the target space), which in the harmonic approximation takes the linear symmetrized form

[TABLE]

$C_{ij,kl}$ is the elastic constant tensor, whose number of independent components is restricted by the symmetry of the crystal. For simplicity, we focus on the isotropic hexagonal lattice, where $C_{ij,kl}$ takes the form

[TABLE]

and characterized by two independent elastic constants, namely the Lamé coefficients, $\lambda$ and $\mu$. As we discuss in Appendix B, an external stress $\sigma^{e}_{ij}({\bf r})$ is included through an additional term $-\sigma^{e}_{ij}u_{ij}$, here focusing on the case of a vanishing external stress.

II.2 Topological defects

In addition to the single-valued elastic phonon fields, the crystal also exhibits topological defects – disclinations and dislocations, captured by including a nonsingle-valued part of the phonon distortion field.

Disclinations are topological defects associated with orientational order. A disclination at a point ${\bf r}_{0}$, illustrated in Fig. 2(a), is defined by a nonzero closed line-integral of the gradient of the bond angle around ${\bf r}_{0}$:

[TABLE]

or equivalently in a differential form:

[TABLE]

measuring the deficit/surplus bond angle, $(2\pi/p)s$, with $s$ the integer disclination charge in a $C_{p}$ symmetric crystal. In the case of a hexagonal lattice, $p=6$. In the above equation, $s({\bf r})$ is the disclination charge density.

Dislocations are vector topological defects associated with translational order. A dislocation at ${\bf r}_{0}$ with a Burgers vector-charge ${\bf b}_{n}$ (that is an elementary lattice vector), as illustrated in Fig. 2(b), is defined by a closed line-integral

[TABLE]

or equivalently in the differential form,

[TABLE]

where ${\bf b}({\bf r})$ is the Burgers charge density. As illustrated in Fig. 2(b), a dislocation is a disclination dipole, and it is therefore energetically less costly than a bare disclination charge.

II.3 Vector Coulomb gas formulation

In the presence of topological defects, the distortion field ${\bf u}({\bf r})$ is not single-valued. Along with the associated strain tensor, the distortion field can be decomposed into the single-valued elastic phonon and the singular parts,

[TABLE]

To include these topological and phonon degrees of freedom we focus on the partition function (taking $k_{B}T=1$, i.e., measuring coupling constants in units of thermal energy),

[TABLE]

where the trace over ${\bf u}({\bf r})$ implicitly includes both phonons and topological defects by allowing nonsingle-valued distortions. In the second form, above, we decoupled the elastic energy by introducing a Hubbard-Stratonovich tensor field – the symmetric stress tensor $\sigma_{ij}({\bf r})$ChaikinLubensky , with the resulting Hamiltonian density given by

[TABLE]

Above, for a 2D hexagonal lattice,

[TABLE]

Tracing over the single-valued phonons $\tilde{\bf u}$, enforces the divergenless stress constraint [via $\delta$-function identity $\frac{1}{2\pi}\int_{-\infty}^{\infty}due^{iuf}=\delta(f)$]

[TABLE]

solved with a scalar Airy stress potential, $\phi({\bf r})$,

[TABLE]

Expressing the Hamiltonian density in terms of $\phi({\bf r})$, and integrating by parts in the second linear term, we utilize the defects conditions, Eq. (9) and (11),

[TABLE]

to obtain

[TABLE]

Above $\tilde{C}^{-1}_{ij,kl}=\epsilon_{ia}\epsilon_{jb}\epsilon_{kc}\epsilon_{ld}C^{-1}_{ab,cd}$.

Focussing on dislocations and neglecting the high energy disclination defects, we can straightforwardly integrate out $\phi({\bf r})$ in the partition function, obtaining a dislocations vector Coulomb gas Hamiltonian

[TABLE]

where the tensor Coulomb interaction in Fourier and coordinate spaces is given by

[TABLE]

with $K=\frac{4\mu(\mu+\lambda)}{2\mu+\lambda}$.

Thus, dislocations vector Coulomb gas Hamiltonian in real space reduces to,

[TABLE]

which, in the discrete form and supplemented with core energies (see below) is exactly the vector Coulomb gas model used by Nelson and HalperinNelson 79 and by YoungYoung 79 , as the theory of 2D continuous two-stage melting.

II.4 Dual vector sine-Gordon model

Motivated by the sine-Gordon description of the XY model, we dualize elasticity by transforming the above vector Coulomb gas into a vector sine-Gordon model and re-examine the two-stage continuous 2D melting transition from this complementary approach.

Dislocation and disclination densities on a hexagonal lattice are given as a sum of their discrete charges

[TABLE]

where ${\bf r}_{n}=a(n_{1}\hat{\bf e}_{1}+n_{2}\hat{\bf e}_{2})$, ($n_{1},n_{2}\in\mathbb{Z}$) are triangular lattice vectors spanned by unit vectors $\hat{\bf e}_{1}=\hat{\bf x}$ and $\hat{\bf e}_{2}=\frac{1}{2}\hat{\bf x}+\frac{\sqrt{3}}{2}\hat{\bf y}$, and ${\bf b}_{{\bf r}_{n}}=a(n_{1}\hat{\bf e}_{1}+n_{2}\hat{\bf e}_{2})$, and $s_{{\bf r}_{n}}\in\mathbb{Z}$ are dislocation and disclination charges, respectively.

In terms of these discrete topological defect charges, the Hamiltonian is given by

[TABLE]

where we have added dislocation and disclination core energies $E_{b}=a^{2}\tilde{E}_{b}$ and $E_{s}$ to account for the defects’ energetics at the scales of the lattice constant, not accounted for by the elasticity theoryNelson 79 . The partition function involves an integration over potential $\phi({\bf r}_{n})$ and summation over the dislocation and disclination charges. Following a standard analysisKosterlitz 72 ; Nelson 79 ; CoulombGas ; sineGordon ; ChaikinLubensky ,

[TABLE]

where non-neutral charge configurations vanish automatically after integration over $\phi({\bf r})$. In the last step, above, we have summed over only the positive/negative single charges of dislocation and disclination, and we obtain the dual vector sine-Gordon Hamiltonian,

[TABLE]

Above, the couplings $g_{b}=\frac{2}{a^{2}}e^{-E_{b}},g_{s}=\frac{2}{a^{2}}e^{-E_{s}}$ are proportional to dislocation and disclination fugacities, and the three elementary dislocation Burgers vectors are given by ${\bf b}_{1}=a\hat{\bf x},{\bf b}_{2}=-\frac{a}{2}\hat{\bf x}+\frac{a\sqrt{3}}{2}\hat{\bf y},{\bf b}_{3}=-{\bf b}_{1}-{\bf b}_{2}=-\frac{a}{2}\hat{\bf x}-\frac{a\sqrt{3}}{2}\hat{\bf y}$.

II.5 Vector sine-Gordon duality redux

The above derivation of the dual vector sine-Gordon model departed from the conventional phonon-only elastic model of a 2D crystal, Eq. (5). As discussed in Sec. II.2 target space rotational invariance of the crystal is incorporated by building the theory based on the symmetric tensor part $u_{ij}$, Eq. (6) of the full strain tensor, $\partial_{i}u_{j}$, i.e., forbidding an explicit dependence on the local bond angle $\theta={\frac{1}{2}}\epsilon_{ij}\partial_{i}u_{j}$, which corresponds to an angle of rotation of the crystal.

Alternatively, the rotational invariance of a crystal can be formulated as a gauge-like (minimal) coupling between the full strain tensor $\partial_{i}u_{j}$ and the bond angle $\theta({\bf r})$, encoded in the Hamiltonian density

[TABLE]

It can be straightforwardly verified that an integration over the bond-angle field $\theta$ Higgs’es outRadzihovskyLubensky ; ChaikinLubensky the antisymmetric component of the strain tensor, at long wavelengths recovering the conventional elastic Hamiltonian in (5).

We now decouple the strain and bond elastic terms by introducing two Hubbard-Stratanovich fields – the stress field $\sigma_{ij}$ and the torque “current” $j_{i}$,

[TABLE]

We note that because $\partial_{i}u_{j}$ is not symmetrized, the stress tensor $\sigma_{ij}$ is not symmetric here. In the presence of topological defects, we again decompose the distortion field ${\bf u}$ and the bond angle $\theta$ into the smooth elastic and nonsingle-valued components,

[TABLE]

which allow for dislocation and disclination defects, respectively.

Integrating out the single-valued parts $\tilde{{\bf u}}$ and $\tilde{\theta}$ enforces two constraints

[TABLE]

The first one is solved via a vector gauge field ${\bf A}$ with

[TABLE]

which transforms the second constraint into

[TABLE]

It is then solved by introducing another scalar potential $\phi$, via $j_{i}=\epsilon_{ik}\partial_{k}\phi-A_{i}$. Expressing the Hamiltonian (31) in terms of gauge potentials, ${\bf A}({\bf r})$ and $\phi({\bf r})$, integrating by parts and using the definitions of dislocation ${\bf b}({\bf r})$ and disclination $s({\bf r})$ densities, the Hamiltonian density takes the form

[TABLE]

This model is evidently gauge-covariant under a local transformation,

[TABLE]

Integrating over the vector potential, ${\bf A}({\bf r})$ in the partition function to lowest order yields

[TABLE]

(an effective Higgs mechanismRadzihovskyLubensky ; ChaikinLubensky ) and allows us to eliminate ${\bf A}({\bf r})$ in favor of $\phi({\bf r})$ and to give the effective Hamiltonian density

[TABLE]

which is identical to that found in (II.3), which when supplemented by dislocation and disclination core energies and summed over the defects gives the dual vector sine-Gordon model, Eq. (29).

II.6 Energetics of defects

Inside the crystal state the background defects density vanishes. In terms of the dual defects model, (II.3) and (39), we can simply set the defect charges to zero, $s={\bf b}=0$. In terms of the generalized sine-Gordon model, (29), this corresponds to the irrelevance of both cosines, i.e., vanishing couplings, $g_{s}=g_{b}=0$,

[TABLE]

where $\tilde{C}^{-1}_{ij,kl}=\epsilon_{ia}\epsilon_{jb}\epsilon_{kc}\epsilon_{ld}C^{-1}_{ab,cd}$, and we have specialized it to that of a hexagonal lattice, obtaining (40b) with $K^{-1}=\frac{2\mu+\lambda}{4\mu(\mu+\lambda)}$.

The energy of a single disclination can be obtained by taking $s({\bf r})=2\pi\delta^{2}({\bf r})$. Solving the corresponding Euler-Lagrange equation for $\phi$ gives

[TABLE]

which for the energy of a single disclination in a crystal state gives a well-known result,

[TABLE]

where $L$ is the linear extent of the crystal.

Similarly, for a single dislocation, such as ${\bf b}({\bf r})={\bf b}_{1}\delta^{2}({\bf r})=a\hat{{\bf x}}\delta^{2}({\bf r})$, the corresponding Euler-Lagrange equation for $\phi$ gives

[TABLE]

where $\theta$ is the angle between the direction of ${\bf k}$ and the $\hat{{\bf x}}$ axis. Therefore, the energy of a single dislocation in the crystal state is

[TABLE]

The $C_{3}$ rotational symmetry guarantees that the energies for the other single dislocations, ${\bf b}({\bf r})={\bf b}_{2}\delta^{2}({\bf r})$ and ${\bf b}({\bf r})={\bf b}_{3}\delta^{2}({\bf r})$ are identical.

III Renormalization-group analysis of the melting transition

As discussed in the Introduction, and calculated above, within the crystal state with a vanishing background defect density, the energy of a single dislocation scales as $E_{b}\sim\ln L$, while the energy of a single disclination scales as $E_{s}\sim L^{2}$, where $L$ is the linear extent of the crystal. Thus, as discovered by Nelson and HalperinNelson 79 , above the critical melting temperature $T_{m}$, the dislocations will unbind first, while the disclinations remain confined, leading to the orientationally ordered hexatic liquid, which is stable in a finite temperature range $T_{m}<T<T_{\text{hex}}$.

More formally, within the dual sine-Gordon model this is reflected by the irrelevance of the disclination cosine operator at the Gaussian fixed line. Computing its average in a system of size $L$, we indeed find

[TABLE]

This analysis [that can be more formally elevated to a renormalizaiton group (RG) computation] demonstrates that the disclination cosine operator, $g_{s}$, is strongly irrelevant around the Gaussian fixed line, i.e., when the dislocation cosine, $g_{b}$, is small, corresponding to the absence of screening of disclinations by dislocations.

III.1 Crystal-hexatic melting transition

Thus, within the crystal and near the crystal-to-hexatic transition, we can neglect the disclination cosine, setting $g_{s}=0$, reducing the effective Hamiltonian to

[TABLE]

The RG analysis of this model is more convenient in an equivalent description in terms of a divergenceless vector field ${\bf A}=\hat{\bf z}\times{\bm{\nabla}}\phi$,

[TABLE]

with the constraint ${\bm{\nabla}}\cdot{\bf A}=0$ imposed energetically via a “mass” term $\frac{\alpha}{2}\left({\bm{\nabla}}\cdot{\bf A}\right)^{2}$ added to the Hamiltonian, with $\alpha\to\infty$ taken at the end of the calculation. Interestingly, our model is mathematically closely related to that for the roughenning transition of a crystal pinned by a commensurate substrate, studied by OhtaOhta80 , and by Levin and DawsonLevin 90 .

Specializing $C^{-1}_{ij,kl}$ to a hexagonal lattice, Eq. (16), the Hamiltonian reduces to

[TABLE]

where the couplings are

[TABLE]

In the physical limit $\alpha\to\infty$, the dislocation-free, Gaussian propagator is given by,

[TABLE]

a purely transverse form, with the transverse projection operator, $P_{ij}^{T}({\bf q})=\delta_{ij}-\frac{q_{i}q_{j}}{q^{2}}$ consistent with (23), encoding the target-space rotational invariance of the crystalNelson 79 ; Young 79 .

To describe the melting transition we need to include dislocations, encoded in the $g_{b}$ cosine operator. Although at low temperature, (corresponding to large elastic constants and small $K^{-1}$) a perturbative expansion in $g_{b}$ is convergent, (i.e., the fixed line $g_{b}=g_{s}=0$ is stable), it breaks down for $K$ below a critical value, indicating an entropic proliferation of large dislocation pairs for $T>T_{m}$.

To treat this high-temperature nonperturbative regime requires an RG analysis. Relegating the details to Appendix A, here we present the highlights of the analysis and its results. To control the divergent perturbation theory, we employ the momentum-shell coarse-graining RG by decomposing the vector field ${\bf A}({\bf r})$ into its short-scale and long-scale modes, $A_{i}({\bf r})=A_{i}^{<}({\bf r})+A_{i}^{>}({\bf r})$, with

[TABLE]

where the ultra-violet cutoff $\Lambda=2\pi/a$, and the rescaling factor $b>1$ defines the width of the momentum shell, $\Lambda/b<q<\Lambda$. Following a standard analysis, we integrate short scale modes $A_{i}^{>}({\bf r})$ out of the partition function, obtaining a coarse-grained Hamiltonian for the long-scale modes, $A_{i}^{<}({\bf r})$, with the renormalized coupling $K_{R}^{-1}(b)$, $B_{R}(b)$ and $g_{b_{R}}(b)$ satisfying

[TABLE]

valid to second-order in $g_{b}$. The Greens function appearing above is given by

[TABLE]

and $J_{i}$ factors are defined as,

[TABLE]

where $I_{0}(x)$ and $I_{1}(x)$ are modified Bessel functions.

It is convenient to examine an infinitesimal form of these RG equations by taking $b=e^{\delta l}$ with $\delta l\ll 1$. Near the melting critical point $K_{R}^{-1}(l\to\infty)\equiv K_{R*}^{-1}=\frac{a^{2}}{16\pi},g_{b_{R}}(l\to\infty)\equiv g_{b}^{*}=0$, this then gives RG differential flow equations for the dimensionless coupling constants $\overline{K}^{-1}(l)=K^{-1}/a^{2}$, $\overline{B}(l)=B/a^{2}$ and $\overline{g}_{b}(l)=g_{b}a^{2}$ ,

[TABLE]

Using the definitions in terms of the dimensionless Lamé elastic constants $\overline{\mu}=\mu a^{2}$, $\overline{\lambda}=\lambda a^{2}$, and the fugacity $y$,

[TABLE]

our equations reduce exactly to the seminal RG flows for the inverse shear modulus, $\overline{\mu}^{-1}(l)$, the inverse bulk modulus $[\overline{\mu}(l)+\overline{\lambda}(l)]^{-1}$, and the effective fugacity $y(l)$, respectively,

[TABLE]

first derived by Nelson and HalperinNelson 79 , and YoungYoung 79 .

Following a standard analysisKosterlitz 72 ; Nelson 79 , the characteristic correlation length $\xi_{\text{xtal-hex}}$ near the critical point at $T\to T_{m}^{-}$ can be extracted from the above RG flows, and it is given by

[TABLE]

with the hexagonal lattice exponent given by

[TABLE]

and $c$ a nonuniversal constantNelson 79 .

III.2 Hexatic-isotropic liquid transition

Dislocation-unbinding above the melting temperature, $T_{m}$ destroys the crystal order, restoring continuous translational symmetry. The plasma of unbound dislocations drives the shear modulus to zero, but retains the quasi-long ranged orientational order and the associated bond orientatonal stiffness. Inside this orientationally-ordered hexatic fluid ($T_{m}<T<T_{hex}$) $g_{b}$ is driven to strong coupling, suppressing ${\bf A}$ fluctuations, and allowing us to approximate the dislocation cosine by its harmonic form. With disclinations reinstated, the resulting effective Hamiltonianin takes the standard scalar sine-Gordon form:

[TABLE]

where $J\equiv\frac{3}{2}a^{2}g_{b}$.

Alternatively, we can get to this dual hexatic Hamiltonian by noting that above the critical melting temperature, $T_{m}$, dislocations (dislination dipoles) unbind, leading to an orientationally ordered (a hexatic) fluid. Since the dislocations then appear at finite density, their Burgers charge, ${\bf b}({\bf r})$ can be treated as a continuous (rather than a discrete) vector field. Going back to (27), carrying out a Gaussian integral over a continuous field ${\bf b}_{\bf r}$, and summing over discrete disclination charges, $s_{{\bf r}_{n}}$ we again obtain the hexatic Hamiltonian presented above.

Utilizing $\tilde{H}_{hex}$ we observe that within the hexatic phase, the energy of a disclination, screened by the plasma of proliferated dislocations, is reduced significantly from that of the crystal (where it diverges as $L^{2}$) to $E_{s}^{hex}\sim J_{R}\ln L/a$. Thus, the hexatic-isotropic fluids transition is of the conventional Kosterlitz-Thouless typeKosterlitz 72 ; Nelson 79 ; ChaikinLubensky , taking place at $T_{hex}=(72/\pi)J_{R}(T_{hex})$. Above $T_{hex}$ the fluid is isotropic, characterized by short-ranged translational and orientational correlations.

IV Summary and conclusion

In this paper, starting with a descripton of a crystal in terms of its elasticity and topological defects, we derive a corresponding dual vector sine-Gordon model. In the latter, the disclinations and dislocations are captured by cosine operators of the Airy stress potential and its gradient. The relevance of the latter dislocation cosine signals the continuous Kosterlitz-Thouless-Halperin-Nelson-Young melting transition of a crystal into a hexatic fluidKosterlitz 72 ; Nelson 79 ; Young 79 . The subsequent relevance of the former disclination cosine captures the hexatic-to-isotropic fluid Kosterlitz-Thouless transition, as outlined in Sec. I.2 and illustrated in Fig. 1. Our complementary analysis reproduces straightforwardly the results of Nelson and HalperinNelson 79 and YoungYoung 79 , including the correlation functions, defects energetics, renormalization-group flows, and the correlation length exponent $\overline{\nu}$.

We expect that the simplified vector sine-Gordon formulation, presented here will be useful in further detailed studies of, e.g., the external stress, the defects dynamics, the substrate, and the dynamics of the melting transitiondynamicsMeltingPRB80 .

Acknowledgements.

We acknowledge useful discussions with Michael Pretko, and we thank him for valuable input on the manuscript. This work was supported by the Simons Investigator Award from the James Simons Foundation and by NSF MRSEC Grant No. DMR-1420736. L.R. also thanks the KITP for its hospitality as part of the Fall 2016 Synthetic Matter workshop and sabbatical program, during which this work was initiated and supported by NSF Grant No. PHY-1125915.

Appendix A Derivation of RG equations

In this appendix we present the details of the renormalization-group analysis of the vector sine-Gordon model for the dislocation unbinding transition,

[TABLE]

where the coupling constants are $K^{-1}=\frac{2\mu+\lambda}{4\mu(\mu+\lambda)}$ and $B=\frac{\lambda}{4\mu(\mu+\lambda)}$. The transversality constraint $\nabla\cdot{\bf A}=0$ is imposed energetically by taking $\alpha\to\infty$ at the end of the calculation.

In the physical limit $\alpha\to\infty$, the dislocation-free, Gaussian propagator is given by,

[TABLE]

with the transverse projection operator, $P_{ij}^{T}({\bf q})=\delta_{ij}-\frac{q_{i}q_{j}}{q^{2}}$.

To carry out the momentum-shell RG analysis, we decompose vector field ${\bf A}({\bf r})$ into its short- and long-scale modes, $A_{i}({\bf r})=A_{i}^{<}({\bf r})+A_{i}^{>}({\bf r})$, with

[TABLE]

where the ultra-violet cutoff $\Lambda=2\pi/a$, and the rescaling factor $b>1$ defines the width of the momentum shell, $\Lambda/b<q<\Lambda$.

Integrating out the short-scale modes, $A_{i}^{>}({\bf r})$, the partition function reduces to integration over the long-scale modes, $A^{<}_{i}({\bf r})$ with an effective Hamiltonian,

[TABLE]

where $Z_{0}^{>}=\int[d{\bf A}^{>}]e^{-H_{0}[{\bf A}^{>}]}$ is the harmonic part of the partition function of the short-scale modes with the quadratic Hamiltonian,

[TABLE]

and, the coarse-grained effective Hamiltonian, $H_{<}[{\bf A}^{<}]$ of the long-scale modes given by,

[TABLE]

We drop the last term, $-\ln Z_{0}^{>}$, which is a field independent correction to the free energy that does not affect the flow of the coupling constants. We then compute $H_{<}[{\bf A}^{<}]$ in terms of corrections to elastic constants $\mu,\lambda$ and dislocation fugacity $g_{b}$, arising from coarse-graining $-\ln\langle e^{-H_{\text{int}}[{\bf A}^{<}+{\bf A}^{>}]}\rangle_{0}^{>}$.

We expand $\langle e^{-H_{\text{int}}[{\bf A}^{<}+{\bf A}^{>}]}\rangle_{0}^{>}$ to second order in $g_{b}$,

[TABLE]

finding,

[TABLE]

These are straightforwardly evaluated by Gaussian integration, giving, to first order

[TABLE]

and to second order the connected part,

[TABLE]

where ${\bf v}={\bf r}_{1}-{\bf r}_{2}$, and we have defined

[TABLE]

For $0<v<a$, we approximate the short-scale averaged correlation function by its value at $v=0$,

[TABLE]

For $v>a$, we have the real space correlation functionNelson 79 ,

[TABLE]

Therefore, we can evaluate $G^{>}_{nm}({\bf r}_{1}-{\bf r}_{2})$ explicitly, finding

[TABLE]

This thus gives to second-order in $g_{b}$,

[TABLE]

The above double integral, $\int d^{2}r_{1}\int d^{2}r_{2}(...)$ simplifies using the fact that $G_{nm}^{>}({\bf r}_{1}-{\bf r}_{2})$ is short-ranged, vanishing for $|{\bf r}_{1}-{\bf r}_{2}|$ larger than $b/\Lambda\sim ba$, since $G_{nm}^{>}({\bf r})$ is defined to be composed of Fourier modes only within a thin momentum-shell, $\Lambda/b<q<\Lambda$. Consequently, $\left[e^{\pm G_{nm}^{>}({\bf r}_{1}-{\bf r}_{2})}-1\right]$ is also small everywhere but in the range, $|{\bf r}_{1}-{\bf r}_{2}|\sim b/\Lambda\sim ba$. To utilize these observations we change variables $\{{\bf r}_{1},{\bf r}_{2}\}$ to a their sum and difference,

[TABLE]

obtaining

[TABLE]

and,

[TABLE]

This allows the following simplifications of (75),

[TABLE]

Comparing (79) to the component of the long-scale Hamiltonian $H[{\bf A}^{<}]$, we can extract the corresponding corrections for the coupling constants $K^{-1},B$, and $g_{b}$. To this end, ignoring the field-independent terms, we obtain

[TABLE]

We note that above, we have written $G_{nn}^{>}(0)$ simply as $G^{>}(0)$, since it takes the same value for all elementary Burgers vectors ${\bf b}_{n}$, $n=1,2,3$. Further simplifications lead to the desired form

[TABLE]

where the coefficients $J_{1},J_{2}$ and $J_{3}$ are defined as:

[TABLE]

Above $I_{0}(x)$, $I_{1}(x)$ are modified Bessel functions, we have dropped higher harmonic operators, and have taken, $\delta l\equiv\ln b\ll 1$.

The above analysis now allows us to identify the renormalized couplings $K_{R}^{-1},B_{R}$ and $g_{b_{R}}$,

[TABLE]

obtained to second-order in $g_{b}$. The corresponding RG differential flow equations for the dimensionless couplings $\overline{K}^{-1}(l)=K^{-1}/a^{2}$, $\overline{B}(l)=B/a^{2}$, $\overline{g}_{b}(l)=g_{b}a^{2}$ are then given by

[TABLE]

Near the melting critical point, $\overline{K}^{-1}(l\to\infty)=\frac{1}{16\pi}$, $\overline{g}_{b}(l\to\infty)=0$, we define the reduced temperature, $x(l)=\frac{16\pi}{\overline{K}}-1$, and fugacity, $y(l)=e^{-E_{cs}}=\frac{1}{2}\overline{g}_{b}$. Near the melting point, their flow equations are given by

[TABLE]

where $c_{1}=e^{2}\left(2I_{0}(2)-I_{1}(2)\right)=21.937...$, and $c_{2}=eI_{0}(2)=6.1965...$ are numerical constants, consistent with Halperin and NelsonNelson 79 .

Following a standard analysisKosterlitz 72 ; Nelson 79 , the characteristic correlation length $\xi_{xtal-hex}$ near the critical point at $T\to T_{m}^{-}$ can be extracted from above RG flows, giving

[TABLE]

with the hexagonal lattice exponent given by $\overline{\nu}=0.3696\ldots$ and $c$ a nonuniversal constant.

Using the expressions of $K$ and $B$ in terms of the dimensionless Lamé elastic constants $\overline{\mu}=\mu a^{2}$ and $\overline{\lambda}=\lambda a^{2}$,

[TABLE]

the RG flow equations for the inverse shear modulus, $\mu^{-1}(l)$ and the inverse bulk modulus, $[\mu(l)+\lambda(l)]^{-1}$, near the critical point are then given by,

[TABLE]

Appendix B Elasticity of 2d crystal subject to an external stress field

As the crystal is subject to an external stress tensor field $\sigma^{e}_{ij}({\bf r})$, we need to add an external term into the elastic energy functional

[TABLE]

In the presence of topological defects, the distortion field ${\bf u}({\bf r})$ is not single-valued. It and the associated strain tensor can be decomposed into the single-valued elastic phonon and the singular part,

[TABLE]

To include these topological and phonon degrees of freedom we focus on the partition function (taking $k_{B}T=1$, i.e., measuring coupling constants in units of thermal energy),

[TABLE]

where the trace over ${\bf u}({\bf r})$ in the partition function (as in summing/integrating over the degrees of freedom of the theory) implicitly includes both phonons and topological defects by allowing nonsingle-valued distortions. In the second form, above, we decoupled the elastic energy by introducing a Hubbard-Stratonovich tensor field – the symmetric stress tensor $\sigma_{ij}({\bf r})$, with the resulting Hamiltonian density given by

[TABLE]

and $C^{-1}_{ij,kl}=-\frac{\lambda}{4\mu(\mu+\lambda)}\delta_{ij}\delta_{kl}+\frac{1}{2\mu}\delta_{ik}\delta_{jl}$ for a 2D hexagonal lattice.

Tracing over the single-valued phonons $\tilde{\bf u}$, enforces the divergenceless stress constraint

[TABLE]

solved with a scalar Airy stress potential, $\phi({\bf r})$,

[TABLE]

Expressing the Hamiltonian density in terms of $\phi({\bf r})$, and integrating by parts in the second linear term, we utilize the defects conditions, Eqs. (9) and (11),

[TABLE]

to obtain

[TABLE]

where $\tilde{C}_{ij,kl}=\epsilon_{ia}\epsilon_{jb}\epsilon_{kc}\epsilon_{ld}C_{ab,cd}$, $\tilde{\tilde{C}}_{ij,kl}=\epsilon_{ia}\epsilon_{jb}C_{ab,kl}$, and we have used integration by parts. The total elastic energy functional is therefore composed of an internal part and an external part, with the internal part

[TABLE]

where we have added the dislocation and disclination core energies $E_{cb}$ and $E_{cs}$ to account for their short-scales, and the external part

[TABLE]

Alternatively, we can also start by formulating the elastic energy in terms of both the full strain tensor $\partial_{i}u_{j}$ and the bond angle $\theta({\bf r})$,

[TABLE]

and get the same result following the procedure of Sec. IIC.