Universality of the Berezinskii–Kosterlitz–Thouless type of phase transition in the dipolar XY-model

The standard paradigm of phase transitions in a planar system of electrically neutral particles representing an effective XY-model spin dictates that long-range order does not exist at any finite temperature [1, 2]. On the other hand, below a critical temperature ${{T}_{c}}$, spin–spin correlations in two spatial dimensions decay in a power-law fashion. As a consequence, short-range cooperative phenomena, such as superfluidity, can exist at these low temperatures. Below the critical point, thermal excitations predominantly occur in the form of vortex–antivortex pairs. Due to the attraction between vortex and antivortex, the vortex–antivortex pairs remain bound at $T<{{T}_{c}}$, and the gas thus remains superfluid, while the dissociation of vortex–antivortex pairs at $T>{{T}_{c}}$ leads to an exponential decay of correlations. This standard scenario, corresponding to the celebrated Berezinskii–Kosterlitz–Thouless (BKT) phase transition [3–5], however conventionally only applies when the particle interact by short-range (contact) potentials.

Spin ordering phenomena in planar systems with long-range forces mediated by the dipole–dipole interaction between particles carrying static dipole moments μ, have been extensively studied in the past, cf e.g., [6–24]. The dipole–dipole interaction in many cases tends to stabilize the long-range order against thermal fluctuations, and the ground state of the spin system may thus be spontaneously polarized [6–10], or acquire various structures [11–23]. So far, most efforts aimed at understanding a possible phase transition in such systems were a combination of renormalization group (RG) arguments with phenomenological approaches [7–15], or Monte Carlo and Molecular Dynamics simulations [16–21]. However, it is fair to say that a complete picture of the nature of the phase transition in the particular case of the two-dimensional (2D) dipolar XY model, where the dipoles sample the full anisotropy of the dipolar interaction, is elusive. While [19, 20] state that a BKT type transition is observed from their numerical results, a physical mechanism explaining the phase transition is missing. To clarify the nature of this transition is the main aim of our work.

Possible applications, once a thorough understanding of the 2D dipolar XY model has been obtained, span a broad range [16], of which we quote just a few. Besides the commonly studied ferromagnetic and ferroelectric thin films (cf e.g., [11, 23]), these include ultracold dilute gases. There are currently major efforts undertaken to cool heteronuclear molecules with large dipole moments to quantum degeneracy [25–27], which will ultimately lead to studies of the dipolar BKT transition [28–31]. Effective spin models with dipolar interactions in the highly controllable environment of ion traps have received considerable attention as well [32, 33]. In a biological context, dipole–dipole interactions determine, for example, the formation of a 2D hydrogen-bond network and the large-scale polarization of water molecules in the hydration layers of proteins [34, 35]. Finally, we note that confinement in isolated vortex–antivortex pairs by the string tension is one of the rare instances outside the realm of Quantum Chromodynamics, in which linear interaction potentials between a particle and its antiparticle, in Quantum Chromodynamics between quark and antiquark, occur [36, 37].

For contact interactions, vortex and antivortex in a vortex–antivortex pair interact by a logarithmic potential, which in superfluids is due to kinetic energy of the flow, and the attraction force decays with the inverse distance. This logarithmic interaction is the primary requirement for the occurrence of the BKT transition. On the other hand, in the presence of a dipolar interaction term in the Hamiltonian, a spatially constant attraction force ${{K}_{0}}$ between vortex and antivortex in an isolated pair occurs [8, 10] (also see below), which bears an obvious potential importance for the phase transition, which has been overlooked in most previous investigations of the 2D dipolar XY model, with the notable exception of Maier and Schwabl [10].

The latter detailed consideration of the 2D dipolar XY-model within a RG treatment has shown, quoting [10], that the '...flow diagram of the ferromagnetic transition is strikingly similar to the Kosterlitz-Thouless transition.' On the other hand, the existence of the linear interaction between vortices, the vortex confinement, led the authors of [10] to the statement that a novel phase transition, distinct from BKT, takes place in the 2D dipolar XY model. In other words, the dipolar interaction is argued there to be relevant for the nature of the phase transition. In the following, we critically examine the latter conclusion. We demonstrate that the sole effect of the dipole–dipole interaction is that a vortex–antivortex pair dissociation transition of the BKT type occurs at slightly higher temperatures. To this end, we use an analytical model, applicable sufficiently close to the transition point, backed up by numerical simulations for the full range of temperatures. We show that the apparently inconsistent pictures of confinement of isolated vortex–antivortex pairs and occurrence of a phase transition driven by the familiar BKT mechanism of pair dissociation can be made fully consistent with each other if one correctly accounts for shielding of the bare vortex–antivortex pair tension at finite temperatures. The shielding effect of the linear attraction potential within a large pair is mediated by the large number of small vortex–antivortex pairs in which it is immersed. The vortex–antivortex interaction only remains linear at small distances R between the vortices in a pair, $R<{{r}_{0}}$, where the length parameter ${{r}_{0}}$ is defined in equation (14) below and discussed in more detail in appendix B; it however becomes logarithmic at sufficiently large vortex–antivortex pair size, $R>{{r}_{0}}$. Hence, we conclusively demonstrate that the dipolar interaction contribution is irrelevant in the sense of the RG.

The polarization (effective spin) states of the 2D dipolar XY model, which correspond to a vortex–antivortex pair gas, are described by the continuous two-component vector-field $\boldsymbol{s}(\boldsymbol{r})=s\cos \theta {{\boldsymbol{e}}_{x}}+s\sin \theta {{\boldsymbol{e}}_{y}}$, representing the hydrodynamic, coarse-grained average of the spin, where θ is azimuthal angle. The topological charge inside a 2D contour C is defined as usual to be $Q={{(2\pi )}^{-1}}{{\mathop{\int }^{}}_{\mathcal{C}}}\text{d}\theta$. It is the number of rotations of vector $\boldsymbol{s}$ 'winding' around the vortex core along closed contour $\mathcal{C}$; $Q\ne 0$ implies that some vortices are inside $\mathcal{C}$. One can locate the vortex core by contracting the contour to a point. If $Q\ne 0$ is retained in the process, the vortex core resides at this point. For example, the upper vortex in figure 1 (left) has Q = 1 and the lower one $Q=-1$.

Note that the spin vortices are dual to the conventional superfluid vortices in the sense that the spin vector is always oriented perpendicular to the corresponding 'fluid' flow direction, cf figure 1. When the wave function of the superfluid is written in Madelung representation $\Psi =\sqrt{n}\exp [i\theta ]$, for constant density n we have the kinetic energy $\propto \mathop{\int }^{}\text{d}f{{(\nabla \theta )}^{2}}$ (where df is the surface element). For spin vortices, the kinetic energy reads $\propto \mathop{\int }^{}\text{d}f{{\mathop{\sum }^{}}_{a,b}}{{\nabla }_{a}}{{s}_{b}}{{\nabla }_{a}}{{s}_{b}}$, see equation (2) below, which transforms for constant s into an identical expression, $\propto \mathop{\int }^{}\text{d}f{{(\nabla \theta )}^{2}}$.

Zoom In Zoom Out Reset image size

In dimensionless form, the hydrodynamic free energy functional (playing the role of the hydrodynamic 'Hamiltonian' of the system) is

$$\begin{eqnarray}&&{{G}_{S}}\left[ \boldsymbol{s}(\boldsymbol{r}) \right]={{G}_{0}}+{{G}_{dd}}\end{eqnarray} \tag{ 1 }$$

The short-ranged term ${{G}_{0}}$ generally reads ($a,b=x,y$) [38]:

$$\begin{eqnarray}&&{{G}_{0}}=\frac{1}{2}\mathop{\int }^{}\text{d}f[C\underset{a,b}{\mathop{\mathop{\sum }^{}}}\,{{\nabla }_{a}}{{s}_{b}}{{\nabla }_{a}}{{s}_{b}}+C\prime {{(\nabla \cdot \boldsymbol{s})}^{2}}+g{{\left( {{s}^{2}}-1 \right)}^{2}}],\end{eqnarray} \tag{ 2 }$$

where $C,C\prime$ are constants, and $g>0$ is a 'contact interaction' coupling. The 'elastic' constants $C,C\prime$ determine the ground-state spin structure, which for $C>C\prime \geqslant 0$ becomes ferromagnetic. For concreteness, we set $C\equiv 1$ as the unit of energy, and $C\prime =0$, as well as put $g\equiv 1$, choosing as our unit of length the vortex core size (an ultraviolet cutoff). This does not affect our results on the nature of the phase transition (also see the remark after equation (4)). Note that in the hydrodynamic expression equation (4), the constants $C,C\prime$ in equation (2) only affect the ultraviolet cutoff, which is irrelevant for the nature of the phase transition; the same applies to the coupling g in equation (2).

The dipolar interaction energy functional entering the hydrodynamic free energy in equation (1) is given by [6, 7, 10]

$$\begin{eqnarray}&&{{G}_{dd}}=\Lambda \mathop{\int }^{}\text{d}f\text{d}{{f}^{\prime }}\frac{{{\boldsymbol{\rho}}_{P}}(\boldsymbol{r}){{\boldsymbol{\rho}}_{P}}\left( {{\boldsymbol{r}}^{\prime }} \right)}{|\boldsymbol{r}-{{\boldsymbol{r}}^{\prime }}|}.\end{eqnarray} \tag{ 3 }$$

Here ${{\boldsymbol{\rho}}_{P}}(\boldsymbol{r})=-\nabla \cdot \boldsymbol{s}(\boldsymbol{r})$ is the density of polarization charges, and $\Lambda \propto {{\mu }^{2}}$ represents the dipole–dipole interaction coupling constant. In appendix A, we provide some considerations on the microscopic derivation of the above effective Hamiltonian equation (1).

We first undertake a qualitative consideration of the scaling behavior of the contributions to the hydrodynamic free energy functional of the planar dipolar XY model for a single vortex–antivortex pair. Note that as in superfluids (see [39], for example), the energy ${{G}_{0}}$ increases approximately proportional to ${{Q}^{2}}$. Therefore the existence of vortices with $|Q|>1$ at low temperatures is energetically disfavored. Focusing our attention on the regime $T\approx {{T}_{c}}$, we have, typically, that the vortex–antivortex pair size $R\gg 1$. For these vortex–antivortex pairs, using that $\nabla \cdot \boldsymbol{s}\sim 1/R$ and hence ${{G}_{0}}\sim \log R$, ${{G}_{dd}}\sim R$. Therefore, the energy of a vortex–antivortex pair with $R\gg 1$ increases linearly with their size, and vortex and antivortex attract each other with a constant force [8, 10], in sharp distinction to contact interactions, where the force decreases as the inverse of the vortex–antivortex pair size.

For contact interactions, it is well established that the total free energy of the system can be written in the form of a double sum $(i,j=1\ldots 2N$, where N is the number of vortex–antivortex pairs)

$$\begin{eqnarray}&&{{G}_{S}}=-\frac{1}{2}\underset{i,\,j}{\mathop{\mathop{\sum }^{}}}\,{{q}_{i}}{{q}_{j}}u\left( {{r}_{ij}} \right)=\frac{1}{2}\underset{i}{\mathop{\mathop{\sum }^{}}}\,{{q}_{i}}\Phi \left( {{\boldsymbol{r}}_{i}} \right).\end{eqnarray} \tag{ 4 }$$

where, in the purely contact interaction case, $u(r)=2\pi \log \left( 1+\alpha r \right)$, with α is a constant, and ${{q}_{i}}=\pm 1$ are the topological charges of the vortices [4, 5] located at ${{\boldsymbol{r}}_{j}}=\left( {{x}_{j}},{{y}_{j}} \right)$, ${{r}_{ij}}=\left| {{\boldsymbol{r}}_{i}}-{{\boldsymbol{r}}_{j}} \right|$.

On the other hand, when strong dipolar interactions come into play, which lead to long-ranged and anisotropic forces in the planar dipolar XY model, the validity of the pairwise summation formula equation (4) is, in distinction to the contact interaction case, highly nontrivial and needs to be thoroughly justified. Because equation (4) is at the heart of our analytical description in terms of plasma physics (see below), we have therefore set up the corresponding Langevin Dynamics simulations and ran extensive checks to establish the validity of equation (4) (see for a detailed discussion of the numerical procedure appendix D). For an illustration, typical results for the polarization field and energy distributions are presented in figure 1.

We find that the function

$$\begin{eqnarray}&&u(r)=2\pi \log \left( 1+\alpha r \right)+{{K}_{0}}r\end{eqnarray} \tag{ 5 }$$

describes very well the energy of a vortex–antivortex pair system in the dipolar gas ⁶ . The quantity ${{K}_{0}}\equiv {{K}_{0}}\left( \Lambda \right)$ is the 'vortex–antivortex pair tension' coefficient, depending on the coupling strength $\Lambda$. The logarithmic contribution in $u(r)$ arises from ${{G}_{0}}$ in equation (2) and describes the interaction of vortices in the limit of $\Lambda \to 0$ [3–5]. The same numerical calculations also let us establish the functional dependence of ${{K}_{0}}(\Lambda )$, which is presented in the inset to figure 2 below. We observe a saturation in the increase of ${{K}_{0}}$ with $\Lambda$, an effect explained in detail in appendix D.

Zoom In Zoom Out Reset image size

The 'potential' $\Phi (\boldsymbol{r})$ in (4) is introduced in analogy with the electrostatics of charges:

$$\begin{eqnarray}&&\Phi (\boldsymbol{r})=\underset{j}{\mathop{\mathop{\sum }^{}}}\,{{q}_{j}}F\left( \boldsymbol{r}-{{\boldsymbol{r}}_{j}} \right)=\mathop{\int }^{}F\left( \boldsymbol{r} \right){{\boldsymbol{\rho}}_{T}}\left( {{\boldsymbol{r}}^{\prime }} \right){{\text{d}}^{2}}{{\boldsymbol{r}}^{\prime }}.\end{eqnarray} \tag{ 6 }$$

Here, $\boldsymbol{r}=\boldsymbol{r}-{{\boldsymbol{r}}^{\prime }}$, $F\left( \boldsymbol{r} \right)=-2\pi \ln \left( 1+\alpha R \right)-{{K}_{0}}R$, and the vortex topological charge density ${{\boldsymbol{\rho}}_{T}}(\boldsymbol{r})={{\mathop{\sum }^{}}_{j}}{{q}_{j}}{{\delta }^{\left( 2 \right)}}\left( \boldsymbol{r}-{{\boldsymbol{r}}_{j}} \right)$.

By its definition, the potential $\Phi (\boldsymbol{r})$ satisfies a Poisson type equation, ${{\hat{L}\,}_{\boldsymbol{r}}}\Phi (\boldsymbol{r})={{\boldsymbol{\rho}}_{T}}(\boldsymbol{r})$. Here the linear operator ${{\hat{L}\,}_{\boldsymbol{r}}}$ is defined such that ${{\hat{L}\,}_{\boldsymbol{r}}}F\left( \boldsymbol{r}-{{\boldsymbol{r}}^{\prime }} \right)={{\delta }^{\left( 2 \right)}}\left( \boldsymbol{r}-{{\boldsymbol{r}}^{\prime }} \right)$, which in Fourier representation reads ${{L}_{{\boldsymbol{k}}}}=1/{{F}_{{\boldsymbol{k}}}}={{\left\{ 2\pi {{K}_{0}}{{k}^{-3}}+4{{\pi }^{2}}\alpha /\left[ {{k}^{2}}\left( k+\alpha \right) \right] \right\}}^{-1}}$. Let us follow the electrostatic analogy further. Since the energy of a charge (vortex) q placed in the external potential $\Phi$ is $U=q\Phi ,$ the force, acting on the charge (vortex) is $\boldsymbol{F}=q\boldsymbol{e}$, where the quasi-electric field vector $\boldsymbol{e}=-\nabla \Phi$, the mean field at the charge location. The energy of the pair in this field is

$$\begin{eqnarray}&&U={{q}_{+}}\Phi \left( {{\boldsymbol{r}}_{+}} \right)+{{q}_{-}}\Phi \left( {{\boldsymbol{r}}_{-}} \right)\approx -\boldsymbol{d}\cdot \boldsymbol{e}\end{eqnarray} \tag{ 7 }$$

Here the topological dipole moment of a pair is introduced according to $\boldsymbol{d}={{\mathop{\sum }^{}}_{j}}{{q}_{j}}{{\boldsymbol{r}}_{j}}={{q}_{+}}{{\boldsymbol{r}}_{+}}+{{q}_{-}}{{\boldsymbol{r}}_{-}}\equiv \boldsymbol{r}$, where $\boldsymbol{r}={{\boldsymbol{r}}_{+}}-{{\boldsymbol{r}}_{-}}$. The density of the topological polarization charges is

$$\begin{eqnarray}&&{{\boldsymbol{\rho}}_{TP}}=-\nabla \cdot \boldsymbol{P},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&\boldsymbol{P}(\boldsymbol{r})={{n}_{P}}{\left\langle \boldsymbol{d}\right\rangle }\end{eqnarray} \tag{ 9 }$$

is the polarization vector of the vortex–antivortex pair gas, and ${{n}_{P}}$ is the surface density of vortex–antivortex pairs; $\left\langle ...\right\rangle$ denotes statistical averaging. After averaging the polarization topological charge inside a contour C, ${{Q}_{TP}}=\mathop{\int }^{}\text{d}f{{\boldsymbol{\rho}}_{TP}}$, becomes a fractional number in general, while before averaging it must be integer.

The equation for the potential of a point charge Q placed at the origin is ${{\Phi }_{Q}}(\boldsymbol{r})$ is ${{\hat{L}\,}_{\boldsymbol{r}}}{{\Phi }_{Q}}(\boldsymbol{r})={{\boldsymbol{\rho}}_{T}}$, where ${{\boldsymbol{\rho}}_{T}}=Q{{\delta }^{\left( 2 \right)}}(\boldsymbol{r})+{{\boldsymbol{\rho}}_{TP}}$. In the weak-field approximation (whose validity we explain in detail in appendix C), that is keeping only the linear terms in $\nabla \Phi$ of equation (10), we obtain a general formula for $\boldsymbol{P}(\boldsymbol{r})$:

$$\begin{eqnarray}&&\boldsymbol{P}(\boldsymbol{r})=-\mathop{\int }^{}\text{d}{{f}^{\prime }}\psi \left( \left| \boldsymbol{r}-{{\boldsymbol{r}}^{\prime }} \right| \right)\nabla \Phi \left( {{\boldsymbol{r}}^{\prime }} \right).\end{eqnarray} \tag{ 10 }$$

Here $\psi \left( \left| \boldsymbol{r}-{{\boldsymbol{r}}^{\prime }} \right| \right)$ is a general nonlocal kernel in an isotropic and translationally invariant medium, connecting electric field and polarization [38]. Hence ${{\left( {{\boldsymbol{\rho}}_{TP}} \right)}_{{\boldsymbol{k}}}}=-{{k}^{2}}\psi \left( {\boldsymbol{k}} \right){{\Phi }_{Q{\boldsymbol{k}}}}$, where $\psi \left( {\boldsymbol{k}} \right)=\mathop{\int }^{}\text{d}f$ $\exp \left( -i{\boldsymbol{k}}\cdot \boldsymbol{r} \right)\psi (\boldsymbol{r})$. Our aim is to describe large-scale effects in a slowly varying field ${{\Phi }_{Q}}(\boldsymbol{r})$ corresponding to large r and small k. Henceforth, we can thus take ${{\left( {{\boldsymbol{\rho}}_{TP}} \right)}_{{\boldsymbol{k}}}}\approx -\chi {{k}^{2}}{{\Phi }_{Q{\boldsymbol{k}}}}$, where $\chi \equiv \psi \left( {\boldsymbol{k}}=0 \right)=\mathop{\int }^{}\text{d}f\psi (\boldsymbol{r})$. This approximation holds when the kernel $\psi (\boldsymbol{r})$ decays fast enough with r, such that $\mathop{\int }^{}\text{d}f\psi (\boldsymbol{r})$ converges. The physical reason behind the latter assumption is that small pairs are far from dissociation and, thus, the vortex–antivortex pair gas at any moment can be approximately subdivided into separated vortex–antivortex pairs. In this approximation the polarization vector is given by (9), where the mean dipole moment of a pair $\left\langle \boldsymbol{d}\right\rangle$ is calculated with the help of a Boltzmann type formula with potential energy (7). It is clear from this qualitative consideration that the characteristic length of $\psi (\boldsymbol{r})$ should be of order of the typical small pair dimension, ${{R}_{P}}$: the mean size of a typical small pair with unscreened interaction energy $u(r)$ is given by ${{R}_{P}}={{J}_{1}}/{{J}_{0}}$. Here, we defined thermal averages of moments of the radial distance as follows

$$\begin{eqnarray}\begin{matrix} {{J}_{n}} & = & \mathop{\int }^{}{{\text{d}}^{2}}r{{r}^{n}}\exp \left[ -u(r)/T \right] \\ {} & = & 2\pi \underset{0}{\overset{\infty }{\mathop{\mathop{\int }^{}}}}\,\text{d}r{{r}^{n+1}}{{e}^{-\gamma r}}{{\left( 1+\alpha r \right)}^{-\beta }}, \\ &&\end{matrix}\end{eqnarray} \tag{ 11 }$$

where $\beta =2\pi /T$, $\gamma ={{K}_{0}}/T$. At $T\approx {{T}_{c}}$, all parameters approach unity (in our units), therefore ${{R}_{P}}\sim 1$ as well. Note that close to the phase transition this length scale is comparable with the other typical scale of our problem, $\sim n_{P}^{-1/2}$, therefore there is only one characteristic scale at short distances close to ${{T}_{c}}$. So, for the slowly varying field

$$\begin{eqnarray}&&\boldsymbol{P}(\boldsymbol{r})\approx \chi \boldsymbol{e}(\boldsymbol{r})=-\chi \nabla \Phi (\boldsymbol{r}).\end{eqnarray} \tag{ 12 }$$

Therefore χ is the susceptibility of the vortex–antivortex pair gas. Putting everything together, we have for the potential

$$\begin{eqnarray}&&{{\Phi }_{Q}}(\boldsymbol{r})=Q\mathop{\int }^{}\frac{{{d}^{2}}k}{{{\left( 2\pi \right)}^{2}}}\frac{\exp \left( i{\boldsymbol{k}}\cdot \boldsymbol{r} \right)}{{{L}_{{\boldsymbol{k}}}}+\chi {{k}^{2}}}.\end{eqnarray} \tag{ 13 }$$

We define the distance scale

$$\begin{eqnarray}&&{{r}_{0}}=\frac{1}{2\pi {{K}_{0}}\chi },\end{eqnarray} \tag{ 14 }$$

which involves the vortex–antivortex pair tension ${{K}_{0}}$ and the susceptibility χ of the vortex–antivortex pair gas. At large distances, $r\gg {{r}_{0}}$, the dominant contribution to the integral comes from small values of k, for which the ${{L}_{{\boldsymbol{k}}}}$ term in the denominator is negligible. Therefore, since ${{\left( \ln r \right)}_{{\boldsymbol{k}}}}=-2\pi /{{k}^{2}}$, the potential of the 'charge' at large distances is logarithmic: ${{\Phi }_{Q}}(\boldsymbol{r})=-\left( Q/2\pi \chi \right)\ln \left( r/{{C}_{1}} \right)$, where the constant ${{C}_{1}}\sim {{r}_{0}}$. In the opposite limit, $r\ll {{r}_{0}}$, the potential is linear: ${{\Phi }_{Q}}(\boldsymbol{r})\approx -Q{{K}_{0}}r$. We propose a simple interpolating expression: ${{\Phi }_{Q}}(\boldsymbol{r})\approx -\left( Q/2\pi \chi \right)\ln \left( 1+r/{{r}_{0}} \right)$, so that the energy of a sufficiently large pair of size R is given by

$$\begin{eqnarray}&&U\left( R \right)=-{{\Phi }_{Q\,=\,1}}(R)=\frac{1}{2\pi \chi }\ln \left( 1+\frac{R}{{{r}_{0}}} \right).\end{eqnarray} \tag{ 15 }$$

The distance scale ${{r}_{0}}$ in equation (14) represents an analogue of the Debye shielding radius for 2D interactions of topological charges. We refer the reader for further details on the properties of the distance scale ${{r}_{0}}$ to appendix B.

The standard calculation procedure of the transition temperature for a gas of polarizable vortex–antivortex pairs with interaction (15) gives the following implicit equation for the transition temperature in terms of the susceptibility [3–5], ${{T}_{c}}=1/4\pi \chi$. At the transition temperature, $T={{T}_{c}}$, vortex–antivortex pairs begin to dissociate. This implies that at ${{T}_{c}}-T\ll {{T}_{c}}$ only a small fraction of the pairs is large and close to dissociation. For this reason it is possible to neglect the interactions between the largest vortex–antivortex pairs and calculate the energy of a single large pair approaching its dissociation limit, which is permeated by a cloud of comparatively small bound vortex–antivortex pairs. Hereinafter we will subdivide vortex–antivortex pairs into two classes: small pairs and large, close to dissociation, pairs. The shielding effect arises due to the polarization cloud provided by small pairs, influencing the potential energy of a large pair.

Let us calculate first the polarizability ${{\alpha }_{P}}$ of a single small pair. The energy of a small pair in an external field $\boldsymbol{e}$ equals $V(\boldsymbol{r})=u(r)-\boldsymbol{r}\cdot \boldsymbol{e}$. The average dipole moment of the small pair and the susceptibility of the vortex–antivortex pairs gas within the framework of the weak-field approximation (cf appendix C) are given by the relations $\left\langle \boldsymbol{d}\right\rangle =\mathop{\int }^{}\text{d}f\boldsymbol{r}\exp \left( -V/T \right)/$ $\mathop{\int }^{}\text{d}f\exp \left( -V/T \right)\approx {{\alpha }_{P}}\boldsymbol{e},$ $\chi ={{\alpha }_{P}}{{n}_{P}}$. Here, ${{\alpha }_{P}}={{J}_{2}}/\left( 2T{{J}_{0}} \right)$ is the small pair polarizability we are looking for, where we used the definition equation (11). The potential energy of a small pair in the field of the charge Q, $\left\langle \boldsymbol{d}\right\rangle \nabla \Phi \approx \left\langle \boldsymbol{d}\right\rangle \nabla {{\Phi }_{V}}\left( \boldsymbol{\rho } \right)$, approximately does not depend on the position of the pair, $\boldsymbol{\rho}$, so that ${{n}_{P}}$ can be considered as a constant.

The surface density of vortex–antivortex pairs, ${{n}_{P}}$, at ${{T}_{c}}$ can be calculated as follows. At $T\approx {{T}_{c}}$ the vortex–antivortex pairs only start to dissociate and the fraction of large pairs is small. The typical pair is small, and the interaction $\approx u(r)$ between its vortices is still unscreened. Then, we conclude from basic arguments of statistical mechanics that the partition function of N pairs on a surface with area A is given by $z=z_{1}^{N}/N!\approx {{\left( {{z}_{1}}e/N \right)}^{N}}$, using Stirling's formula, and where ${{z}_{1}}=\mathop{\int }^{}{{\text{d}}^{2}}{{r}_{1}}{{\text{d}}^{2}}{{r}_{2}}\exp \left[ -u\left( |{{\boldsymbol{r}}_{1}}-{{\boldsymbol{r}}_{2}}| \right)/T \right]=A{{J}_{0}}$, which again employs equation (11). Minimization of the free energy ${{F}_{P}}\left( N \right)=-T\ln z$ gives $N={{z}_{1}}$, ${{n}_{P}}=N/A={{J}_{0}}$. Using this result and the relation ${{T}_{c}}={{(4\pi \chi )}^{-1}}$, there follows $2\pi {{J}_{2}}=1$, leading to an implicit equation for the critical temperature

$$\begin{eqnarray}&&4{{\pi }^{2}}\mathop{\int }^{}_{0}^{\infty }\text{d}r\frac{{{r}^{3}}\exp \left( -{{K}_{0}}r/{{T}_{c}} \right)}{{{\left( 1+\alpha r \right)}^{2\pi /{{T}_{c}}}}}=1.\end{eqnarray} \tag{ 16 }$$

Next we compare the semianalytical solution for ${{T}_{c}}$ from (16) with the results of the numerical calculation, cf figure 2. As described in detail in appendix D, we used Langevin Dynamics and Binder's method to calculate ${{T}_{c}}$ in a series of simulations with increasingly larger realizations of the model system. In the non-dipolar case, ${{K}_{0}}=0$ (or $\Lambda =0$), the equation (16) yields ${{T}_{c}}\approx 0.7$. From our numerical calculations, we found that ${{T}_{c}}=0.85$. At large $\Lambda$ the transition temperature tends to the constant value ${{T}_{c}}(\infty )\approx 1.4$. From (16) we derive ${{T}_{c}}\left( \infty \right)\approx 1.3$. Note that the RG arguments of [10] gave the much larger prediction ${{T}_{c}}(\infty )=2\pi$. This in turn implies that it is rather difficult to account for all essential Feynman graphs in order to adequately describe the shielding effect in a RG calculation (which in fact applies equally well to the short-range case $\Lambda =0$).

We note that the equation (16) is based on the assumption of the noninteracting small pairs. Rigorously, though, as ${{R}_{P}}\sim n_{P}^{-1/2}$, this fails around $T\simeq {{T}_{c}}$. Hence we expect this approximation to give only an order of magnitude estimate for the susceptibility χ, and the equality sign in equation (16) is to be replaced in that rigorous sense by an ≈. On the other hand, we find from our numerical calculations that the approximation leading to equation (16) still proves to be rather reliable for a sufficiently accurate prediction of the critical temperature, cf figure 2.

We have demonstrated both numerically and analytically using an analogy to plasma physics, that vortex–antivortex pairs in the 2D dipolar XY-model dissociate at a critical temperature in a manner familiar from the BKT transition. This is due to the dipole–interaction-induced linear confinement potential in an isolated large vortex–antivortex pair being shielded by a gas of small pairs, in which the large pair becomes immersed around the transition point. Therefore, the logarithmic attraction between vortices in large pairs is restored. By obtaining a physically transparent scenario, we have therefore provided an unambiguous proof that the BKT mechanism is applicable to a much broader class of systems than hitherto established. Our simulations, combined with the analytical approach presented above, give a rigorous confirmation of the qualitative assumptions discussed in [42]. Shielding of the linear interaction in a vortex–antivortex pair implies that the dipole interaction term ${{G}_{dd}}$ is irrelevant (in the sense of the RG), and that the 2D dipolar XY-model belongs to the same universality class as the contact interaction BKT-model corresponding to ${{G}_{dd}}=0$. This conclusion is confirmed in appendix D by arguments based on the numerical calculation of the Binder cumulant. While the scale hierarchy, discussed in our paper only qualitatively, can more adequately be formulated in the language of the RG [10], we have provided clear evidence that a correct understanding of the physical nature of the phase transition can not be obtained within the RG approach.

We finally stress that the planar dipolar XY model (1) fundamentally differs from the commonly studied 2D system with all dipoles oriented perpendicular to the plane ⁷ , cf e.g., [43], which is a scalar model. Also, it differs from a 'purely' dipolar model (see [24], for example), whose Hamiltonian does not contain the short-range (gradient) term, ${{G}_{0}}$, in equation (1). Therefore, the physics behind these models is fundamentally different from the dipolar XY case. For example, the ground state for square lattice in [24] is antiferroelectric, but we have a ferroelectric ground state.

The work of AYuV, AET, LIM, and POF was supported by Quantum Pharmaceuticals. The research of URF was supported by the NRF of Korea, Grant Nos. 2010-0013103 and 2011-0029541.

We briefly outline in what follows the derivation of the Hamiltonian equation (1), underlying our analysis, from microscopics. The interaction energy of ${{N}_{M}}$ polar molecules with dipole moments ${{\mu }_{A}}$, where $A=1,\ldots ,{{N}_{M}}$, equals

$$\begin{eqnarray}&&V=\frac{1}{2}\underset{\alpha \beta ,A\ne B}{\mathop{\mathop{\sum }^{}}}\,{{\mu }_{A\alpha }}{{\mu }_{B\beta }}{{f}_{\alpha \beta }}(\boldsymbol{\rho})={{V}^{L}}+{{V}^{S}},\end{eqnarray} \tag{ A.1 }$$

where $\alpha ,\beta$ are the spatial indices, $\boldsymbol{\rho}={{\boldsymbol{r}}_{A}}-{{\boldsymbol{r}}_{B}}$, and ${{f}_{\alpha \beta }}(\boldsymbol{\rho})=f_{\alpha \beta }^{L}(\boldsymbol{\rho})+f_{\alpha \beta }^{S}(\boldsymbol{\rho})$. Here $f_{\alpha \beta }^{L}(\boldsymbol{\rho})$ describes the long-range interaction of molecules $f_{\alpha \beta }^{L}(\boldsymbol{\rho})=({{\delta }_{\alpha \beta }}-3{{n}_{\alpha }}{{n}_{\beta }})/{{\boldsymbol{\rho}}^{3}}$, $\boldsymbol{n}=\boldsymbol{\rho}/\boldsymbol{\rho}$. The function $f_{\alpha \beta }^{S}(\boldsymbol{\rho})$ represents the short-range part, for example the hydrogen bond interaction in the case of water molecules. In the continuum approximation of hydrodynamics, when the polarization vector, ${{\boldsymbol{P}}_{d}}(\boldsymbol{r})=\mathop{\sum }^{}_{A=1}^{{{N}_{M}}}{{\boldsymbol{\mu}}_{A}}{{\delta }^{(D)}}(\boldsymbol{r}-{{\boldsymbol{r}}_{A}})$, is considered as a slowly varying function, the latter replaces ${{\mu }_{A}}$, and the summation over molecules $\mathop{\sum }^{}_{A=1}^{{{N}_{M}}}$ is replaced by an integration $\mathop{\int }^{}{{\text{d}}^{D}}{{r}_{A}}{{n}_{M}}$, where ${{n}_{M}}$ is the moment density and D the spatial dimension. The hydrodynamic approximation is capable to describe the long-range effects which are considered in our paper. After an integration by parts, the long-range part ${{V}^{L}}$ takes the form of equation (3).

In the short-range term, we can take $\boldsymbol{\rho}\ll {{\boldsymbol{r}}_{AB}}=|{{\boldsymbol{r}}_{A}}+{{\boldsymbol{r}}_{B}}|/2$, and the energy density is expanded in ρ. Writing

$$\begin{eqnarray}&&{{V}^{S}}=\frac{1}{2}n_{M}^{2}\mathop{\int }^{}{{\text{d}}^{D}}{{r}_{A}}{{\text{d}}^{D}}{{r}_{B}}{{P}_{d,\alpha }}\left( {{\boldsymbol{r}}_{A}} \right){{P}_{d,\beta }}\left( {{\boldsymbol{r}}_{B}} \right)f_{\alpha \beta }^{S}\left( \boldsymbol{\rho} \right),\end{eqnarray} \tag{ A.2 }$$

we use the following definitions and expansions

$$\begin{eqnarray}&&{{\boldsymbol{r}}_{A}}={{\boldsymbol{r}}_{AB}}+\frac{1}{2}\boldsymbol{\rho},\ {{\boldsymbol{r}}_{B}}={{\boldsymbol{r}}_{AB}}-\frac{1}{2}\boldsymbol{\rho},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&{{\text{d}}^{D}}{{r}_{A}}{{\text{d}}^{D}}{{r}_{B}}={{\text{d}}^{D}}{{R}_{AB}}{{\text{d}}^{D}}\boldsymbol{\rho},\ {{\text{d}}^{D}}\boldsymbol{\rho}={{\boldsymbol{\rho}}^{D-1}}\text{d}\boldsymbol{\rho}d{{\Omega }_{n}},\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{{P}_{d,\alpha }}\left( {{\boldsymbol{r}}_{A}} \right)\approx {{P}_{d,\alpha }}\left( {{\boldsymbol{r}}_{AB}} \right)+\frac{1}{2}\left( \boldsymbol{\rho}{{\nabla }_{R}} \right){{P}_{\alpha }}\left( {{\boldsymbol{r}}_{AB}} \right),\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&{{P}_{d,\beta }}\left( {{\boldsymbol{r}}_{B}} \right)\approx {{P}_{d,\beta }}\left( {{\boldsymbol{r}}_{AB}} \right)-\frac{1}{2}\left( \boldsymbol{\rho}{{\nabla }_{R}} \right){{P}_{\beta }}\left( {{\boldsymbol{r}}_{AB}} \right).\end{eqnarray} \tag{ A.6 }$$

The tensor $f_{\alpha \beta }^{S}\left( \boldsymbol{\rho} \right)$ depends only on the distance vector $\boldsymbol{\rho}$. Due to space isotropy it should have the same form in any Cartesian frame, therefore

$$\begin{eqnarray}&&f_{\alpha \beta }^{S}\left( \boldsymbol{\rho} \right)=A\left( \boldsymbol{\rho} \right){{\delta }_{\alpha \beta }}+B\left( \boldsymbol{\rho} \right){{n}_{\alpha }}{{n}_{\beta }}.\end{eqnarray} \tag{ A.7 }$$

After integration over $\text{d}{{\Omega }_{n}}$, the nonvanishing contributions contain even powers of $\boldsymbol{n}$ only. Finally, identifying $\boldsymbol{s}={{\boldsymbol{P}}_{d}}/(\mu {{n}_{M}})$ (assuming that all ${{\mu }_{A}}$ have the magnitude μ), the integration $\mathop{\int }^{}_{0}^{\infty }{{\boldsymbol{\rho}}^{D-1}}\text{d}\boldsymbol{\rho}...$ reproduces the bilinear terms in equation (2).

The quartic 'self-interaction' term, $\propto {{({{s}^{2}}-1)}^{2}}$, stems from a spin saturation effect, $s\to 1$ for large moment densities ${{n}_{M}}$. This saturation can be explained, for example, by the Langevin formula $s=L\left[ \frac{\mu {{E}_{d}}}{{{k}_{B}}T} \right]$, with the Langevin function $L(x)=\coth (x)-1/x$, yielding $s\to 1$ for a large polarizing electric field ${{E}_{d}}$.

The polarization topological charge density of the vortex–antivortex pair gas close to a single topological charge Q equals: ${{\boldsymbol{\rho}}_{TP}}=-\nabla \cdot \boldsymbol{P}\approx -\chi Q{{K}_{0}}/r$. We conclude that the total charge inside a circle of radius r is given by ${{Q}_{t}}(r)=Q+{{Q}_{TP}}\approx Q\left( 1-2\pi {{K}_{0}}\chi r \right)$. From this qualitative consideration we thus come to an important conclusion: the charge is essentially shielded, ${{Q}_{t}}\approx 0$, which occurs at a distance scale $r\sim {{r}_{0}}$. Close to the phase transition ${{r}_{0}}\sim {{R}_{P}}\sim 1$. The polarization of the vortex–antivortex pair gas inhibits the linear attraction within large vortex–antivortex pairs which would prevail with shielding not taken into account. Hence we are led to conclude that the phase transition associated with the dissociation of pairs is qualitatively very similar to the BKT transition in a system with $\Lambda \to 0$.

Starting from the order of magnitude estimate above, we now consider a more rigorous approach. According to (4) two unshielded, probe charges at a distance r interact as $G=-{{q}_{1}}{{q}_{2}}{{K}_{0}}r={{q}_{1}}{{\Phi }_{{{q}_{2}}}}(r)$ (in the presented qualitative consideration we neglect the logarithmic term in $u(r)$), where ${{\Phi }_{{{q}_{2}}}}(r)=-{{q}_{2}}{{K}_{0}}r$ is the topological potential of the charge ${{q}_{2}}$. Similarly, for the point charge Q placed at the origin and the probe charge q: $G=q{{\Phi }_{Q}}(r)$, where ${{\Phi }_{Q}}(\boldsymbol{r})\approx -{{Q}_{t}}(r){{K}_{0}}r$. From here and equations (8), (12), we conclude

$$\begin{eqnarray}&&{{\boldsymbol{\rho}}_{TP}}(r)\approx -\chi {{K}_{0}}\text{ }\!\!\Delta\!\!\text{ }\left( r{{Q}_{t}} \right)=-\frac{\chi {{K}_{0}}}{r}\frac{\text{d}}{\text{d}r}\left[ r\frac{\text{d}}{\text{d}r}\left( r{{Q}_{t}} \right) \right].\end{eqnarray} \tag{ B.1 }$$

This yields a differential equation for ${{Q}_{t}}(r)$:

$$\begin{eqnarray*}&&{{Q}_{t}}(r)=Q+\mathop{\int }^{}_{0}^{r}2\pi r\text{d}r{{\boldsymbol{\rho}}_{TP}}(r)\end{eqnarray*}$$

$$\begin{eqnarray}&&=\;Q-2\pi \chi {{K}_{0}}r\frac{\text{d}}{\text{d}r}\left( r{{Q}_{t}} \right).\end{eqnarray} \tag{ B.2 }$$

Its solution is given by

$$\begin{eqnarray}&&{{Q}_{t}}(r)=\frac{Q{{r}_{0}}}{r}{{e}^{{{r}_{0}}/r}}\mathop{\int }^{}_{{{r}_{0}}/r}^{\infty }dx\frac{{{e}^{-x}}}{x}.\end{eqnarray} \tag{ B.3 }$$

This total charge monotonically diminishes with r from ${{Q}_{t}}=Q$ at r = 0 to Q = 0 at $r=\infty$. From the formula above for the effective charge shielding we again conclude that its characteristic scale is ${{r}_{0}}$.

We verify in this part of the appendix the applicability of the weak-field approximation. Using equation (15), the average interaction between the vortices at finite temperature is $\left\langle U\left( R \right)\right\rangle =z_{P}^{-1}\mathop{\int }^{}_{0}^{\infty }\text{d}R\;U\left( R \right)g\left( R \right)={{T}_{c}}\left( 1-4\tau \right)/\left[ \left( -\tau \right)\left( 1-2\tau \right) \right]$, where ${{z}_{P}}=\mathop{\int }^{}_{0}^{\infty }\text{d}Rg\left( R \right)$, $g(R)=R\exp [-U(R)/T]$, and the relative temperature $\tau =(T-{{T}_{c}})/{{T}_{c}}$. On the other hand, the typical size of a close to dissociation large pair, ${{R}_{DP}}$, can be estimated from the relation $\left\langle U\left( R \right)\right\rangle \equiv {{\left( 2\pi \chi \right)}^{-1}}\ln \left( 1+{{R}_{DP}}/{{r}_{0}} \right)$. Therefore, next to the phase transition, $|\tau |\ll 1$, the dimension of a dissociated pair is large, ${{R}_{DP}}\sim {{r}_{0}}\exp (1/|\tau |)$, and its topological 'electric' field is small, $\boldsymbol{e}=-\nabla U\left( R \right)\propto \exp (-1/|\tau |)\ll 1$. A condition for the applicability of the weak-field approximation is therefore the existence of a small parameter $\exp (-1/|\tau |)$, representing the ratio of the typical topological 'electric' field of large pair, $\sim |\nabla U(R)|$ to the field inside a small pair, $\sim |\nabla u(R)|$.

At $|\tau |\ll 1$, the dissociating pair is so (exponentially) large, that most of the small pairs inside the large pair are far away from the vortex sources of the field such that the field created by these sources is small. Therefore, the weak-field approximation is applicable close to the transition temperature.

We studied the thermodynamics of our model using Langevin Dynamics [44]. Using a discretized representation of the Hamiltonian at every grid point $\alpha =1,2,...,{{N}_{g}}$, we performed a fixed temperature run of sufficient length to get reliable averages. The calculations were performed using periodic boundary conditions on square lattices L×L with the number of independent nodes ${{N}_{g}}={{L}^{2}}$. The Langevin Dynamics dynamical equations are given by:

$$\begin{eqnarray}&&\frac{\text{d}{{s}_{\alpha }}}{\text{d}t}=-\gamma \frac{\partial }{\partial {{s}_{\alpha }}}({{G}_{0}}+{{G}_{dd}})+{{\zeta }_{\alpha }},\end{eqnarray} \tag{ D.1 }$$

where γ is a constant that determines the time scale of relaxation. The stochastic thermal noise terms satisfy $\left\langle {{\zeta }_{\alpha }}(t)\right\rangle =0$ and $\left\langle {{\zeta }_{\alpha }}(t){{\zeta }_{\beta }}(t\prime )\right\rangle =2\;T\gamma {{\delta }_{\alpha \beta }}\delta (t-t\prime )$. In the Langevin Dynamics simulations we use a second-order Runge–Kutta algorithm. The equations of motion above are integrated numerically with the sufficiently small discrete time step $\Delta t=0.005$. To compute the dipole–dipole interaction term in equation (2) in an efficient $O({{N}_{g}}\ln {{N}_{g}})$ way we used a NumPy FFTW realization [45].

The software was first used to check the assumptions leading to equation (4). First we consider the case of zero temperature, T = 0, and checked the additivity rule expressed in (3) numerically by using an imaginary time relaxation method, which is capable to find local minima of the free energy functional on a configuration space of 2D vectors ${{\boldsymbol{s}}_{\alpha }}$ taken in nodes of a dense square lattice. Following [46], the initial approximation was specified by the complex-valued skyrmion type expression ${{s}_{0}}(z)=2W/(1+|W{{|}^{2}})$, with $W=\mathop{\prod }^{}_{j=1}^{2N}{{(z-{{z}_{j}})}^{{{q}_{j}}}},$ and $z=x+iy,$ ${{z}_{j}}={{x}_{j}}+i{{y}_{j}}$. For a gas of pairs the total charge vanishes, ${{\mathop{\sum }^{}}_{j}}{{q}_{j}}=0$. Therefore, $W\to 1$ at $z\to \infty$, which agrees with the physical boundary condition $\boldsymbol{s}\to \left( 1,0 \right)$ at $r\to \infty$. At the first stage of imaginary time propagation the cores of the vortex–antivortex pairs begin very slowly to approach which, finally, leads to their annihilation. Due to this mutual attraction, the vortex–antivortex pair therefore is not a real local minimum. Any initial configuration inevitably transfers at T = 0 to the absolute ground state $\boldsymbol{s}\to \left( 1,0 \right)$. To stabilize vortex–antivortex pairs, 'pinning' of vortices was used by adding to the free energy a term ${{G}_{\text{pin}}}={{\mathop{\sum }^{}}_{j}}{{G}_{j}},{{G}_{j}}=\mathop{\int }^{}\text{d}f{{V}_{j}}(\boldsymbol{r}){{\left[ \boldsymbol{s}(\boldsymbol{r})-{{\boldsymbol{s}}_{0}}(\boldsymbol{r}) \right]}^{2}}$, where ${{V}_{j}}(\boldsymbol{r})={{V}_{0}}\exp \left[ -{{\left( \boldsymbol{r}-{{\boldsymbol{r}}_{j}} \right)}^{2}}/{{a}^{2}} \right]$, $a\sim 1$. We investigated multiple configurations with different numbers of vortex–antivortex pairs using the expression in equation (4), and reproduced numerically the total energy within a small error, less than $\sim 5%$ (examples of typical distributions are presented in figure D1). The string tension constant ${{K}_{0}}$ was calculated by analyzing the mutual forces which prevail in a vortex–antivortex pair by calculating averages of the derivatives of the pinning potential with respect to the positions of the vortices. The accuracy we were able to achieve is limited by the perturbations introduced by the pinning potential, and is sufficient to prove the reliability of the additivity rule equation (4), also cf figure D1. We found that at $\Lambda \gtrsim {{\Lambda }_{CR}}\simeq 0.4$, corresponding to ${{K}_{0}}\left( {{\Lambda }_{CR}} \right)\approx 3.1$, a physical instability of the single vortex–antivortex pair configuration arises: a new small vortex–antivortex pair is spontaneously created in the center of a large vortex–antivortex pair. With further increase of $\Lambda$, the 'parent' vortices are immersed into a cloud of small polarized vortex–antivortex pairs, i.e. dipoles with zero total topological charge. Polarization of vortex–antivortex pairs leads to a net topological charge density (cf the discussion after equation (9)) and, hence, to a reduced increase of the line tension ${{K}_{0}}$ with $\Lambda$. Ultimately, this leads to a saturation effect: ${{K}_{0}}\left( \Lambda \right)\approx {{K}_{0}}\left( {{\Lambda }_{CR}} \right)\approx 3.1$ at $\Lambda \gtrsim {{\Lambda }_{CR}}$ (see the inset of figure 2). The spontaneous creation of pairs follows also from the expression (4): at $\Lambda >{{\Lambda }_{CR}}$ the energy of the large ($R\gg 1$) pair decreases with the emergence of a new small pair in the center. This occurs at $\Lambda ={{\Lambda }_{CR}}$, if we choose $\alpha ={{K}_{0}}\left( {{\Lambda }_{CR}} \right)/2\pi \approx {{K}_{0}}\left( \infty \right)$/$2\pi$.

Zoom In Zoom Out Reset image size

To explore the critical behavior of the model depending on the dipole–dipole interaction coupling constant $\Lambda$ numerically we use Binder's method [47] and calculate the parameter ${{U}_{L}}=1-\left\langle {{\boldsymbol{s}}^{4}}\right\rangle /3/\left\langle {{\boldsymbol{s}}^{2}}\right\rangle ^{2}$ ('Binder's cumulant') versus temperature. Here $\boldsymbol{s}=\boldsymbol{s}/{{N}_{g}}$, $\boldsymbol{s}=\mathop{\sum }^{}_{\alpha =1}^{{{N}_{g}}}{{\boldsymbol{s}}_{\alpha }}$ and the statistical averaging is done in a manner equivalent to the average over polarization configurations obtained with Langevin Dynamics. The intersection point of Binder's cumulants for different values of the system size L gives ${{T}_{c}}$, as shown in figure D2. In the symmetric phase, $T>{{T}_{c}}$, ${{U}_{L}}=0+O\left( 1/A \right)$ as $A\to \infty$, where A is the surface area. In the symmetry-broken phase, $T<{{T}_{c}}$, ${{U}_{L}}=2/3+O\left( 1/A \right)$. At the critical point, ${{U}_{L}}$ tends towards a universal value $0<U_{L}^{\star }<2/3$, which is specific for each model, and is determined by its universality class [47, 48]. To determine the universality class of the model at hand it seems natural to simply compare our value $U_{L}^{\star }\approx 0.621$ with that for the pure BKT-model without dipole–dipole interaction. According to detailed numerical results [49, 50], the magnitude of ${{T}_{c}}$ does not depend on the simulation details, but, in fact, the value $U_{L}^{\star }$ itself '...depends sensitively on boundary conditions, details of the clusters used in calculating the cumulant, and symmetry of the interactions or, here, lattice structure...,' quoting [49]. As a consequence, we have to compare our result with others under the same conditions. We know two such results for $\Lambda =0$: $U_{L}^{\star }\approx 0.61$ for the XY-model [49], and $U_{L}^{\star }\approx 0.62$ for the generalized XY-model [51]. As for the present dipolar XY-model, one has $U_{L}^{\star }\approx 0.62$ in a wide range of $\Lambda$ values [19, 20].

Zoom In Zoom Out Reset image size