Instability of Abrikosov lattice due to nonanalytic core reconstruction of vortices in Bosonic superfluids
IInstability of Abrikosov lattice due to nonanalytic core reconstruction of vortices inBosonic superfluids
Avraham Klein, Oded Agam, and Igor L. Aleiner The Racah Institute of Physics, The Hebrew University of Jerusalem, 91904, Israel Physics Department, Columbia University, New York, NY 10027, USA (Dated: October 23, 2018)We study the impact of the non-analytic reconstruction of vortex cores on static vortex structuresin weakly coupled superfluids. We show that in rotating two-dimensional systems, the Abrikosovvortex lattice is unstable to vortex core deformation: Each zero of the wave function becomes a cutof finite length. The directors characterizing the orientations of the cuts are themselves ordered insuperstructures due either to surface effects or to interaction with shear deformations of the lattice(spiral structure). Similar instability may be also observable in clean superconducting films.
PACS numbers: 67.10.-j, 67.10.Jn,67.85.-d ,67.85.De
Introduction – A rapidly rotating bosonic superfluidforms into a dense lattice of vortices rotating as a rigidbody - the Abrikosov lattice [1, 2]. The favorable latticeis triangular, and its stability and vibrational modes werefirst studied by Tkachenko [3, 4]. Tkachenko waves werefirst detected in superfluid helium [5, 6], and later wereobserved directly in rotating cold atomic condensates [7].The triangular lattice structure minimizes the energyof the currents created by the quantized vortices in therotating frame. Alternatively, using the well-known anal-ogy between quantized vortices in two-dimensional (2D)superfluids and point charges in 2D electrodynamics [8–10], it is the minimum energy configuration of pointcharges in a neutralizing uniform background field. How-ever, the point-charge approach does not take into ac-count the possibility of core deformation in the vortices,which endows each vortex with an additional dipole-likedegree of freedom.In the strong coupling regime, namely when core sizeof the vortex is of order to the interparticle distance, thiscore deformation is negligible. However, Bose-EinsteinCondensates of cold atoms may be realized in the weakcoupling limit, were core size is much larger than theinterparticle distance. Here, it has been recently shownthat vortices in 2D bosonic superfluids experience non-analytic reconstruction of their cores when moving withrespect to the ambient superfluid [11]. In this Letter weshow that this core reconstruction leads to an instabilityof the Abrikosov lattice, and identify the new ordering.
Non-analytic reconstruction of the vortex core – To be-gin with, let us summarize the relevant findings of theanalytic theory [11] of vortices in superfluids at the weakcoupling regime, n ξ (cid:29)
1, where n is the 2D bosonicdensity and ξ is the healing length: i) The low energydynamics are not exhausted by the position of the vor-tex attached to the ambient flow but must include alsothe motion of the vortex with respect to this flow; ii)This motion can be characterized by the kinetic momentaˆ p = (ˆ p x , ˆ p y ); iii) The momentum dependence of the ki- netic energy of the vortex circular motion is non-analytic H ( p ) = p M v ( p ) ; M v ( p ) mnξ ≡ πα ln (cid:18) (cid:126) n ξ p (cid:19) , (1)where m is the boson mass, and α = 0 . . . . . Thenon-analytic dependence of the added mass , M v , on themomentum p is associated with the weak solutions of theGross-Pitaevskii equation inside the vortex core [11], seealso [12] for singular vortex antivortex configurations.An effective theory[11] describing a superfluid systemwith N vortices is conveniently written using Popov’sformalism [9] which maps 2D superfluidity to 2D nonlin-ear electrodynamics. In this mapping, vortices becomecharged particles characterized by the coordinates R k ( t ),and charge 2 π (cid:126) σ k with σ k = ± E is related to the superfluid current j as E = − ˆ (cid:15) j (where ˆ (cid:15) is the antisymmetric tensor of the sec-ond rank acting on the spatial coordinates), and the mag-netic field, B , is the boson density, n . The physical fields E , B are related to the vector A and scalar ϕ potentialsby usual relations E = − ∂ t A − ∇ ϕ , and B = ∇ × A . Inthese variables the Lagrangian L (describing the systemover linear scales much larger than the healing length ξ )of a rotating superfluid with angular frequency Ω is L = N (cid:88) j =1 L v (cid:16) p j ; σ j ; ˙ R j ; { E ; B ; A ; ϕ } r = R j (cid:17) + (cid:90) d r L f ; L v = p · (cid:18) ˙ R − ˆ (cid:15) E B (cid:19) + 2 π (cid:126) σ (cid:16) A · ˙ R − ϕ (cid:17) − H ( p ) ; L f = m (cid:34) ( E + Ω B r ) B − c ( B − n ) n (cid:35) , (2)where c ≡ (cid:126) / ( mξ ) is the sound velocity. Equation (2) isdictated by the Galilean invariance of the system and thetransformation to a rotating frame. Its gauge invarianceis the vorticity conservation. The corresponding action (cid:82) L dt should minimized subject to the condition ϕ =0 at the boundary of the 2D superfluid. (The physical a r X i v : . [ c ond - m a t . o t h e r] A ug meaning of this condition is that the energy of currentgenerated due to placing a vortex in the system is countedfrom the energy of the vortex on the boundary.)The ˆ (cid:15) E B term in L v deserves some discussion. In thefirst place, it highlights the fact that the spectrum ofthe H ( p ) is related to the reconstruction of the core ofthe vortex moving with respect to ambient flow and thusshould vanish for the vortex moving together with theflow, ˙ R = ˆ (cid:15) E /B . Next, unlike charged particles in vac-uum, the vortex can be at rest even for p (cid:54) = 0 providedthat ambient flow is finite. Finally, any reconstructionof the core gives rise to an emergent dipole moment2 πσ j (cid:126) d j , with the displacement vector d j = σ j ˆ (cid:15) p j π (cid:126) B . (3)
Energy of an array of vortices with deformed cores – Inwhat follows we assume that N vortices are sufficientlyfar apart so that one may assume the density to be con-stant, B = n . For simplicity, we assume that the 2 D superfluid is confined in a circle of radius R (cid:29) ξ , andneglect the effect of image vortices, whose main effect isto reduce the effective radius of the system by order of a (cid:112) ln( a/ξ ), where a is the lattice constant [13] . Then,looking for a stationary solution ˙ A = 0; ˙ R = 0, we min-imze the action with respect to ϕ and obtain the energyof the vortices (with σ j = sgnΩ) as a function of theirpositions and displacement vectors: U = π (cid:126) n m (cid:20) N ln Rξ − N m | Ω | R (cid:126) + u (cid:21) ; u = N (cid:88) j =1 (cid:34) ε (cid:32) d j ξ (cid:33) + m | Ω | (cid:126) (cid:0) R j + 2 d j · R j (cid:1)(cid:35) + N (cid:88) i (cid:54) = j (cid:18) d i · ∂∂ R i (cid:19) (cid:18) d j · ∂∂ R j (cid:19) ln ξ | R i − R j | , (4)with the dimensionless function (cid:15) ( x ) given by ε ( x ) = 4 α x ln(1 /x ) . (5)For d j = 0, Eq. (4) reduces to the well known en-ergy of logarithmically repelling particles in a parabolicconfinement potential, yielding the configuration of anAbrikosov lattice with non-deformed cores. Minimiza-tion of the energy (4) with respect to the vortex positions R i gives ∂∂ R i N (cid:88) j (cid:54) = i ln 1 | R i − R j | + m | Ω | (cid:126) N (cid:88) j R j = 0 . (6)The solution of this equation, R j = r j , produces thesites of the lattice corresponding to local energy minima.It is known that the global minimum corresponds to r j comprising a triangular lattice [4]. One useful relationcan be derived independently of any lattice structure.Multiplying both sides of Eq. (6) by r i and summingover i , we obtain m | Ω | (cid:126) N (cid:88) j r j = N ( N − , (7)thus, vortices start to overlap at Ω ≥ Ω c = 2 π (cid:126) / ( √ mξ )for the triangular lattice. For Ω < Ω c , the distance be-tween nearest vortices is a = (Ω c / | Ω | ) / ξ. (8) Instability of the Abrikosov lattice due to core deforma-tion – Consider now a shift of the vortices with respect totheir original position, r j → r j + δ j . It will be convenientto work with the variable ∆ j = δ j + d j . The variable r j + ∆ j should be understood as the virtual position ofthe vortices, defined by the current distribution they cre-ate at distances larger than ξ . Expanding the energy inEq. (4) to second order in d j and ∆ j we obtain that thequadratic contribution to u is u = u ∆2 + u D ; u ∆2 = m | Ω | (cid:126) (cid:88) i | ∆ i | + 12 Re N (cid:88) i (cid:54) = j (∆ i − ∆ j ) ( z i − z j ) ; u D = (cid:88) i (cid:20) ε (cid:18) | D i | ξ (cid:19) − m | Ω | (cid:126) | D i | − Re D i ¯ V (2) i (cid:21) ; (9)Here and henceforth we employ complex vector notation,i.e. z = r x + ir y , d = d x + id y , ∆ = ∆ x + i ∆ y , the directoris defined as D i ≡ d i , and the non-diagonal componentof the electric field gradient ¯ V (2) i tensor is given by¯ V ( n ) i = N (cid:88) j (cid:54) = i z i − z j ) n . (10)The u ∆2 term describes the stability of the triangular lat-tice with respect to small displacements of the vortexcenters and it was studied by Tkachenko [3]. The u D term describes the sensitivity of the Abrikosov lattice tocore deformation and we turn to studying this term.We notice that in the quadratic approximation (9) thecores deformations are independent. This is because eachcore deformation d i is followed by the displacement of thevortex center by δ i = − d i so that the currents created bythe vortex at the distance larger than the healing length ξ remain intact. Moreover, far from the edges of thesample ¯ V (2) i → C rotational symmetry.Minimization with respect to | D i | with the logarithmicaccuracy gives the absolute value of the director | d i | = | D i | = D ∗ = ξ e − γ Ω c / | Ω | , (11)where γ = 2 √ α /π ≈ .
74. The quantization of thevortex motion, [11], sets a lower bound on the value ofthe displacement vector | D | n >
1. Thus, Eq. (11) isapplicable in the frequency intervalΩ c / ln( n ξ ) < Ω < Ω c , (12) i.e. it occurs only in the weak coupling regime, n ξ (cid:29) Boundary effects – Consider a finite size lattice with avortex at z = 0 and | z j | < R . Vortex repulsion leads tolattice deformations near the boundary. Far away fromthe boundary, | z j | (cid:28) R , these deformations may be ne-glected [14], and performing finite sum in Eq. (10) wefind [15] ¯ V (2) i (cid:39) a / R / (cid:16) z i R (cid:17) f (cid:18) Ra (cid:19) , (13a)where f ( x ) is an oscillatory function which may haveeither positive or negative sign [15]. Near the boundary,the value of the field gradient (10) can be estimated usinga straight boundary approximation,¯ V (2) i = 1 a (cid:16) ¯ zz (cid:17) η (cid:16) z ¯ z (cid:17) ; R − | z i | < a (13b)where η ( x ) is of order unity and strongly depends on theboundary facets [15]. It is noteworthy that the angulardependence in Eqs. (13) is consistent with the C v sym-metry of the underlying lattice. Substituting Eqs. (13)into Eq. (9) and minimizing with respect to the anglebetween D i and cystallographic direction we find that D i | D i | = sign (cid:20) f (cid:18) Ra (cid:19)(cid:21) (cid:18) ¯ z i | z i | (cid:19) . (14)This behavior was verified by direct numerical minimiza-tion of u , given by Eq. (4) (see Fig. 1).It is clear from Eq. (13a) that the characteristicanisotropy energy u an = D ∗ (cid:12)(cid:12)(cid:12) ¯ V (2) i (cid:12)(cid:12)(cid:12) diminishes with theincreasing size of the system. Therefore, it is imperativeto consider higher order nonlocal terms in D i , as theyestablish long range order. Higher order corrections to the quadratic energy (9)are obtained from further expansion of (4) in d i , ∆ i . Theresults are most easily written in terms of the Nambuspinors R i = (cid:18) ∆ i ¯∆ i (cid:19) ; D i = (cid:18) D i ¯ D i (cid:19) ; ¯ R i = (cid:0) ¯∆ i , ∆ i (cid:1) ¯ D i = (cid:0) ¯ D i , D i (cid:1) (15a) N =
121 N = FIG. 1. Deformation of a finite Abrikosov lattice due to corereconstruction [14]. The panels show a the configuration ofthe Abrikosov lattice obtained by numerical minimization ofthe energy u for a system with 121(left) and 127 (right) vor-tices, taking into account their emergent dipole moments d i .The lines in the figures represent the directors configuration:The length of each line is 10 (cid:112) | D i | , and each line has an anglecorresponding to (arg D i ) / a = 3. where the complex notation was introduced after Eq. (9).The correction δu = u ∆2 + u + u takes the form u ∆2 = 12 ¯ R i ˆ M (2) ij R j (15b) u = − ¯ R i ˆ M (3) ij D j (15c) u = 34 ¯ D i ˆ M (4) ij D j , (15d)and we used Einstein summation over repeated indices.The matrices ˆ M ( n ) ij , n = 2 , , M ( n ) ij = λδ ij δ n, ˆ − (1 − δ ij ) (cid:32) z i − ¯ z j ) n z i − z j ) n ; 0 (cid:33) , (15e)and λ is found from the condition that the minimal eigen-value of ˆ M (2) ij (understood as 2 N × N matrix) is zero inaccordance with Eq. (7) and rotational invariance of thewhole system.Minimizing the energy u ∆2 + u with respect to theGaussian variables R i we obtain the effective energy ofthe directors˜ u D = 12 ¯ D i M ij D j , M ij ≡ − ˆ M (3) ik (cid:104) ˆ M (2) (cid:105) − kl ˆ M (3) lj + 32 ˆ M (4) ij , (16)where the inverse matrix is defined by (cid:104) ˆ M (2) (cid:105) − ik (cid:104) ˆ M (2) (cid:105) kj = δ ij ˆ λ of M are labeled by quasi-momenta q = q x + iq y , ¯ q = q x − iq y . One easily finds λ ± = − | B || B | | B | − | B | ± (cid:12)(cid:12)(cid:12)(cid:12) ¯ B B | B | − | B | + 3 ¯ B (cid:12)(cid:12)(cid:12)(cid:12) , (17)where B n ( q, ¯ q ) = (cid:88) ω (cid:54) =0 e i ( q ¯ ω +¯ qω ) ω n , B ≡ lim | q |→ | B ( q, ¯ q ) | , (18)and ω labels the positions on the lattice, ω kl ≡ a ( k + le i π/ ). The B n terms can be expressed via Weierstrasselliptic functions. Doing so one finds the minimal energyconfiguration corresponds to the q → q → B , B to third order in q (cid:28) u = (cid:90) d r ∂ z Λ ∂ ¯ z Λ4 π + (cid:90) r 2, i.e. for the states withbroken reflectional symmetry. Fig. 2 illustrates the con-figuration for α = + π/ 2. Near the boundary, the surfaceenergy term (which preserves the reflectional symmetry)becomes important. It gives rise to deformation of thesupervortex configurations in a small boundary layer. In conclusion , we showed that the Abrikosov vortexlattice is unstable with respect to non-analytic deforma-tions of the vortex core. The directors characterizingsuch deformations are themselves ordered in “supervor-tex” structures due either to surface effects or to inter-actions with the shear deformations of the vortex lattice.As a final comment, we point out that the descriptionof static order parameter in superconductors also can berecast in a form Eq. (4). The validity of condition Eq.Eq. (12) should be replaced with constraints on the ap-plied magnetic field H c ln( ξ/L T ) < H < H c , where L T is the characteristic length determining the localityof the Ginzburg-Landau description, H c is the criticalmagnetic field. Since in the vicinity of the critical tem-perature ξ/L T diverges, we expect the instability of theAbrikosov lattice towards the deformation of the vortexcore to be observable in clean superconducting films aswell.We are grateful to N. Katz, A. Kuklov, and E. Zel-dov for discussions of the results. This research wassupported by the United States-Israel Binational ScienceFoundation (BSF) grant No. 2012-134, and the IsraelScience Foundation (ISF) grant No. 302/14 (O.A.), andby the Simons foundation (I.A.). [1] A.A. Abrikosov, J. Phys. Chem. Solids, ,199 (1957).[2] A.A. Abrikosov, Zh. Eksp. i Teor. Fiz. , 1442 (1957) [Sov. Phys. JETP , 1174 (1957)];[3] V.K. Tkachenko, Zh. Eksp. Teor. Fiz. , 1875 (1965) [Sov. Phys. JETP, ,1049 (1966)].[4] V.K. Tkachenko, Zh. Eksp. Teor. Fiz. , 1875 (1965) [Sov. Phys. JETP, , 1282 (1966)].[5] C. D. Andereck, J. Chalups, and W. I. Glaberson,Phys.Rev. Lett. , 33 (1980).[6] C. David Andereck and W.I. Glaberson, J. Low Temp.Phys. , 257 (1982).[7] I. Coddington, P. Engels, V. Schweikhard, and E. A. Cor-nell, Phys. Rev. Lett. , 672 (1972) [Sov.Phys. JETP, ,341 (1973)].[9] V.N. Popov. Functional integrals in quantum field the-ory and statistical physics , Volume 8. (Kluwer Academic,2001).[10] V. Ambegaokar, B.I. Halperin, D.R. Nelson, andE.D. Siggia, Phys. Rev. B, , 1806 (1980).[11] A. Klein, I.L. Aleiner, and O. Agam, Annals of Physics, , 195 (2014).[12] L. Radzihovsky, Phys. Rev. Lett. , 247801 (2015).[13] D. Stauffer and A.L. Fetter, Phys. Rev. , 156 (1968).[14] A.A. Koulakov and B.I. Shklovskii, Phys. Rev. B, ,2352 (1998).[15] See Supplemental Material for details of the calculation.[16] C. Kittel, Phys. Rev. , 965, (1946). Supplementary material for “Instability ofAbrikosov lattice due to nonanalytic core recon-struction of vortices in Bosonic superfluids” In this supplementary material we provide details ofthe calculations of V (2) i , B n ( q, ¯ q ), the solution of theAbrikosov lattice deformation in the continuous limit oflarge systems. To shorten notations, in what follows wewill measure distances in units of the lattice constant a . Calculation of ¯ V (2) i In order to calculate¯ V (2) i = ¯ V (2) ( z i ) = (cid:88) j z i − z j ) , (S.1)for a finite system, | z i | < R , we neglect lattice defor-mation near the boundary and use Poisson’s summationformula:¯ V (2) ( z ) = (cid:88) (cid:126)q ddz (cid:90) ∞ R rdr √ (cid:90) π − π dθ e i(cid:126)q · (cid:126)r z − re iθ , (S.2)where (cid:126)q π = √ ( k (cid:126)b + k (cid:126)b ) where b = ( √ , ), and (cid:126)b = (0 , k and k are integers. Here we haveused the fact that on a lattice point, z = z i , the integralover all space is zero (by symmetry), thus the integralover the interior domain r < R may be replaced by anintegral over the exterior domain r > R .Next we separate the sum over (cid:126)b , to a sum over fami-lies using Miller index notation defined as following. Thefacet directions of triangular Bravias lattice are conve-niently described by Miller 3-index notation. Denotingby (cid:126)a = (1 , (cid:126)a = ( − , √ ), and (cid:126)a = ( − , √ ) thethree principal directions of the triangular lattice, a facetgoing through points (cid:126)a i /l i where l + l + l = 0, is de-noted by the Miller indices ( l , l , l ), where l i are integersand at least two of them are coprime. Equivalent crys-tallographic directions which correspond to π/ (cid:104) l , l , l (cid:105) , all of those directions are ob-tained from each other by the cycling permutations ofindices, or by the changing the sign of all indices.Therefore, each vector of the reciprocal lattice is char-acterized by α = { k , k } with fixed ratio of k /k ,where k and k are coprime integers coinciding with firsttwo Miller indices in 3-indices notation for the facets, k , = l , . Then, each family represents a surface in k -space with fixed angle: φ α = arctan (cid:20) √ (cid:18) k k (cid:19)(cid:21) , (S.3)(0 ≤ φ < π ), and an absolute value | (cid:126)q | = q α n where n isa positive integer, and q α = 2 √ (cid:113) k + k k + k . (S.4)Symmetry implies that each one of these families contains6 members obtained by rotations of π/ 3, or by picking upall the possible combinations from < k , k , − k − k > . Thus¯ V (2) ( z ) = ddz (cid:88) α,n (cid:88) m =0 (cid:90) ∞ R rdr √ a (cid:90) π − π dθ e i πrq α n cos [ φ α − θ + m π ] z − re iθ , (S.5)and performing the sum over m we obtain¯ V (2) ( z ) = 2 √ a (cid:88) α,n ddz (cid:90) ∞ R rdr (cid:90) π − π dθ z e i πrq α n cos( φ α − θ ) z − r e i θ . (S.6)For a < | z | (cid:28) R we may expand to leading order in | z | /R to obtain¯ V (2) ( z ) (cid:39) − z √ (cid:88) α,n (cid:90) ∞ R drr (cid:90) π − π dθe i πq α nr cos( φ α − θ ) − i θ = − √ (cid:16) zR (cid:17) (cid:88) α ∞ (cid:88) n =1 e − i φ α J (2 πq α nR ) q α nR . Noticing that except for α = { , } corresponding to φ α = π and α = { , ¯1 } corresponding to φ α = 0, allother contributions come in pairs, such that for some0 < φ n < π/ α (cid:48) such that φ (cid:48) α = π − φ α , and q α = q α (cid:48) we obtainEq. [13a] with f ( x ) = − √ √ x (cid:88) α,n δ α cos(6 φ α ) J (2 πq α nx ) q α n , (S.7)where the sum is over all vectors α such that 0 ≤ φ α ≤ π/ δ α = 1 if α = { , } or α = { , ¯1 } and δ α = 2otherwise.Using the large argument asymptotic limit of theBessel function one can sum over n and express the resultin the form: f ( x ) (cid:39) − √ √ (cid:88) α δ α cos(6 φ α ) q / α ζ (cid:18) − , q α x (cid:19) , (S.8)where ζ ( β, γ ) = 2Γ(1 − β )(2 π ) − β ∞ (cid:88) n =1 sin (cid:0) πγn + β π (cid:1) n − β (S.9)is a generalized Riemann ζ -function. A sample region of f ( x ) is depicted in Fig. S.1.Consider now the behavior of ¯ V (2) i near the boundaryof the system, i.e. when R −| z i | (cid:46) 1. We focus our atten-tion on the points z i which are close to the main facets ofthe triangular lattice (in real space), and ignore the shiftsin the positions of the vortices that are located near theboundary (which is of order of the lattice constant). Thesmall deviations of the boundary from a straight facetcan be considered as a set of quasi random kinks. Thesekinks create a pseudorandom long-range fields in addition 100 101 102 103 104 105 - - x f ( x ) FIG. S.1. The function f ( x ). to a strong short-range field associated with the straightfacets of the lattice.In what follows we shall focus our attention on theprinciple facets: (cid:104) , ¯1 , (cid:105) and (cid:104) , ¯1 , ¯1 (cid:105) in Miller 3-indexnotation.Consider the function,¯ V (2) ( z ) = − (cid:88) | z j | >R z − z j ) (S.10)which is analytic in | z | < R and has second order poleson the lattice points out side the system. To begin with,let us assume a large enough system and a straight facet.Then ¯ V (2) ( z ) assumes a constant value for z which liesalong a straight line parallel to the facet. We denoteby ¯ V (cid:104) l l l (cid:105) ( k ) the value of this function along lattice rowsparallel the facet associated with Miller indices ( l l l ).Here k denotes the row number counted from outside,i.e. k = 1 corresponds to the boundary row, k = 2 to thenext row in, and so on. The analytic structure of ¯ V (2) ( z )implies that close to the (cid:104) , ¯1 , (cid:105) facet:¯ V (cid:104) (cid:105) ( k ) = − ∞ (cid:88) j = k (cid:32) πe i π sin( πe i π j ) (cid:33) = η (cid:104) (cid:105) k e − iφ , (S.11)where φ is the angle of the facet, e − iφ = ¯ zz , while η (cid:104) (cid:105) k = − ∞ (cid:88) j = k π sinh (cid:16) π √ j + i π j (cid:17) . (S.12)In particular, η (cid:104) (cid:105) (cid:39) . η (cid:104) (cid:105) (cid:39) − . − and η (cid:104) (cid:105) (cid:39) − . A similar calculation for the (cid:104) , ¯1 , ¯1 (cid:105) facet gives ¯ V (cid:104) (cid:105) ( k ) = η (cid:104) (cid:105) k e − iφ , where η (cid:104) (cid:105) k = − ∞ (cid:88) j = k π (cid:16) π √ j + i π j (cid:17) . (S.13)In particular, η (cid:104) (cid:105) (cid:39) . η (cid:104) (cid:105) (cid:39) − . η (cid:104) (cid:105) (cid:39) . η (cid:104) l l l (cid:105) on the type of facetstems from the distance between vortices on each facet, which can vary significantly. (E.g. for the (cid:104) (cid:105) facet theseparation is a , whereas for the (cid:104) (cid:105) facet it is √ a .)In order to generalize this result to the case where theboundary is a straight line with random kinks (generatingthe long-range contribution to the field), let us parame-terize the boundary as: z ( t ) = (cid:26) e − i π t + h ( − t ) near (cid:104) , ¯1 , (cid:105) facet it + h ( t ) near (cid:104) , ¯1 , ¯1 (cid:105) facet , (S.14)where h ( t ) is a real function, and let us denote by ˜ h ( t )its integer part.Consider now the (cid:104) , ¯1 , ¯1 (cid:105) facet. Here t = ( z − ¯ z ) / i andtherefore z + ¯ z = 2 h (cid:0) i z − ¯ z (cid:1) . Thus ˜ h (cid:0) i z − ¯ z (cid:1) is analyticfunction in the z -plane with cuts parallel to the real axis.Namely, ∂∂ ¯ z ˜ h (cid:18) i z − ¯ z (cid:19) = i (cid:88) k s k δ (cid:18) i z − ¯ z − y k (cid:19) , (S.15)where y k are the points where ˜ h ( y ) jumps from value toanother and s k = ˜ h ( y k + 0) − ˜ h ( y k − 0) = ± 1. With thesedefinitions, the analytic function (S.10) can be written inthe form:¯ V (2) ( z ) = − (cid:88) j =1+˜ h ( ¯ z − z i ) iπ/ (cid:104) iπ √ ( z + je i π ) (cid:105) − (cid:90) d z πi z − z ∂ ˜ h (cid:0) ¯ z − z i (cid:1) ∂ ¯ z iπ/ (cid:104) iπ √ ( z − − ˜ h (cid:0) ¯ z − z i (cid:1)(cid:105) . (S.16a)Here the first term accounts for the strong and short-range interaction coming from the local principal facet. Thesecond contribution is an integral along the cuts of ˜ h (cid:0) ¯ z − z i (cid:1) . As 1 / sin is a rapidly decaying function along theline of integration, this term represents the long range contribution from pseudo random charges associated with thesurface kinks.A similar expression is obtained for the (cid:104) , ¯1 , (cid:105) facet:¯ V (2) ( z ) = − (cid:88) j =1+˜ h (cid:16) z − ¯ z √ i (cid:17) πe i π sin (cid:104) π ( z − j ) e i π (cid:105) − (cid:90) d z πi z − z ∂ ˜ h (cid:16) z − ¯ z √ i (cid:17) ∂ ¯ z πe i π sin (cid:104) π (cid:16) z − − ˜ h (cid:16) z − ¯ z √ i (cid:17)(cid:17) e i π (cid:105) . (S.16b)In the previous consideration we assumed that the tri-angular lattice is not deformed at all near the boundary.The effect of shifts of the vortices from the triangular lat-tice to their equilibrium position could be investigatedonly numerically. We found that in the finite systemsthe analytic results are modified quantitatively but notqualitatively.Figure S.2 depicts the vortex lattice structure beforeand after relaxation (ignoring the small effect of core re-construction).As is evident from the figure, the most significant re-sult of relaxation is that boundary vortices form a circle,similarly to the situation for point charges interacting via 1 /r interaction, see Ref. [14] in the main text. Thevortices on this circle are approximately uniformly dis-tributed with a separation of one lattice constant. Therest of lattice remains approximately unperturbed. Thestrong dependence of η (cid:104) l l l (cid:105) on the facet stems from theeffective lattice constant of each facet, and therefore thealmost uniform lattice spacing on the boundary circle re-duces the fluctuations. Fig. S.3 shows the argument andthe absolute value of ¯ V (2) i calculated for the outer row ofthe system depicted in Fig. S.2. One can see, the angular¯ z/z dependence of the field is robust, but the modulationof the absolute value is greatly reduced.Finally, the shift of the vortices due to relaxation gen- FIG. S.2. The vortex lattice configuration: Before relaxation(left) and after relaxation (right). The system contains 235vortices, and the circle is only a guide to the eye. - π - π π π - π - π π π arg ( z i ) a r g ( V i ( ) ) - π - π π π ( z i ) V i ( ) FIG. S.3. The argument (up) and the absolute value (down)of V (2) i at the last row of the system shown in Fig. S.2. Thered and blue disks represent the result for relaxed and nonrelaxed configurations of the vortices. erates effective dipole charges near the edge of the sys-tem. The effect of these dipoles in the bulk is smaller byan order of a/R than V (2) i , but they cause strong devia-tions near the boundary. In upper panel of Fig. S.4 weshow the behavior of the argument of ¯ V (2) i as function ofarg( z i ) within the bulk. It shows that there is anothercontribution in addition to the z dependence obtainedin Eq. (13a). The absolute value of ¯ V (2) i as function of | z i | is presented by a log-log plot in the lower panel ofthe figure. The non-relaxed lattice clearly shows the | z | behavior, while the results for the relaxed lattice showsome deviations. However, we attribute these bulk er-rors to numerical inaccuracy of the relaxation algorithm rather than to the true physical effect. - π π π - π - π π arg ( z i ) a r g ( V i ( ) ) × - × - z i V i ( ) FIG. S.4. The behavior of ¯ V (2) i in the bulk the system shownin Fig. S.2. The upper panel shows the argument of ¯ V (2) i for3 < | z i | a < R = 7 . V (2) i for 1 ≤ | z i | a ≤ 6. Thered and blue disks represent the results for relaxed and nonrelaxed configurations of the vortices. Calculation of B n ( q, ¯ q ) and the eigenvalues ofEq. (17) We begin with the expression (18) given as an infinitesum, where ω = nω + mω , ω = 1 , ω = e iπ/ definethe triangular lattice. B n ( q, ¯ q ) inherits several propertiesfrom the lattice structure, B n ( q, ¯ q ) = B n (¯ q, q ); (S.17a) B n ( e i π q, e − i π ¯ q ) = e − i nπ B n ( q, ¯ q ); (S.17b) B n ( q + ν, ¯ q + ¯ ν ) = B n ( q, ¯ q ); (S.17c)where ν = nν + + mν − , ν ± = 2 π (1 ± i √ ) define the re-ciprocal lattice. To identify the general form of B n wenotice that ∂ n B n ∂ ¯ q n = (cid:18) i (cid:19) n (cid:88) ω (cid:54) =0 e i ( q ¯ ω +¯ qω ) (S.18)= (cid:18) i (cid:19) n (cid:34)(cid:88) ω e i ( q ¯ ω +¯ qω ) − (cid:35) . For any q (cid:54) = ν the infinite sum vanishes, and therefore B n ( q, ¯ q ) should have the form: B n ( q, ¯ q ) = − (cid:18) i (cid:19) n n ! (cid:2) ¯ q n + ¯ q n − f ( q ) + ¯ q n − f ( q ) + · · · + f n ( q ) (cid:3) , (S.19)where f k ( q ) are analytic functions of q which are quasi-periodic on the lattice. They can be expressed in termsof the Weierstrass elliptic functions ζ ( q ) and ℘ ( q ), whichhave the following properties: ζ ( q + ν ) = ζ ( q ) + ¯ νβ ; β = 8 π √ ℘ ( q + ν ) = ℘ ( q ); (S.20b)where ν is any vector of the reciprocal lattice. In thelimit q → ζ ( q ) = 1 q + O ( q ) , (S.21a) ℘ ( q ) = 1 q + O ( q ) . (S.21b)Consider first the case n = 2: B ( q, ¯ q ) = 18 (cid:2) ¯ q + ¯ qf ( q ) + f ( q ) (cid:3) . (S.22)In order to have property(S.17c), the shift in ¯ q mustbe compensated by the quasi-periodic behavior of ζ ( q ),(S.20a), i.e. it can only contain the combination ¯ q − βζ ( q ).Thus B ( q, ¯ q ) = 18 (cid:110) [¯ q − βζ ( q )] + C ( q ) (cid:111) , (S.23)where C ( q ) is some periodic function on the lattice. Thisfunction is uniquely determined from the condition that B ( q, ¯ q ) does not diverge at | q | → 0. Taking into accountthe asymptotic behavior (S.21) we obtain: B = 18 (cid:110) [¯ q − βζ ( q )] − β ℘ ( q ) (cid:111) . (S.24)From here it follows that the small q asymptotic behavioris: B ( q, ¯ q ) = − π √ qq + 18 ¯ q ; | q | (cid:28) . (S.25)The calculation of the B ( q, ¯ q ) and B ( q, ¯ q ) follows asimilar procedure. For instance from (S.19) it followsthat B ( q, ¯ q ) = i (cid:110) [¯ q − βζ ( q )] + C ( q ) (cid:111) , (S.26)where again C ( q ) is periodic function which removes thesingularities of the first term. The result is: B = i (cid:104) ˜ ζ − β ˜ ζ℘ ( q ) + β ℘ (cid:48) ( q ) (cid:105) , (S.27)where ˜ ζ = ¯ q − βζ ( q ) . (S.28) FIG. S.5. A density plot of the eigenvalue λ in q space. Bluecolures denotes small values while red color represents highvalue of λ . The hexagon denotes the first Brillouin zone. Aminimum at the K would imply tripling of the lattice constant,while a minimum at the M -point will generate stipe phase.However, the global minimum, is at the Γ -point, similar toferroelectric systems. The asymptotic behavior in this case is: B ( q, ¯ q ) = − iπ √ q q + i ¯ q 48 ; | q | (cid:28) . (S.29)Finally, for n = 4 we have: B = − (cid:104) ˜ ζ − β ˜ ζ ℘ ( q ) + 4 β ˜ ζ℘ (cid:48) ( q ) − β ℘ ( q ) (cid:105) . (S.30)Using the above expressions for B n and Eq. (17) for theeigenvalue λ − one obtains the energy of the lattice perunit cell, as function of q . A density plot of this energyis shown in Fig. S.5. The energy and configuration of large systems Minimizing the functional (19) with respect to ∆within the bulk we obtain ∂ ∆ ∂z = − ∂ F∂z − ∂ Λ ∂z ; ∂ ¯∆ ∂ ¯ z = − ∂ ¯ F∂ ¯ z − ∂ Λ ∂ ¯ z . (S.31)The solution of the first of these equations is: ∂ ∆ ∂z = − ∂F∂z − 2Λ + ¯ λ (¯ z ) , (S.32)where ¯ λ (¯ z ) is the integration constant. However ∂ ∆ /∂z is a scalar and therefore ¯ λ (¯ z ) can only be a constant.This constant represents a rigid rotation of the systemand may be set to be zero. Thus ∂ ∆ ∂z = − ∂F∂z − ∂ ¯∆ ∂ ¯ z = − ∂ ¯ F∂ ¯ z − . (S.33)0Thus within the bulk ρ = 2 √ (cid:18) ∂F∂z + ∂ ¯ F∂ ¯ z + 4Λ (cid:19) , (S.34)and substituting the second equation of (20) we obtain:12 π ∂ Λ ∂z∂ ¯ z = ∂F∂z + ∂ ¯ F∂ ¯ z + 4Λ . (S.35)This equation can be solved by iterations. In the firstapproximation one neglects the gradient term of Λ toobtain: Λ = − (cid:18) ∂F∂z + ∂ ¯ F∂ ¯ z (cid:19) , (S.36)while next iteration givesΛ = − (cid:18) ∂F∂z + ∂ ¯ F∂ ¯ z (cid:19) + 18 π ∂ Λ ∂z∂ ¯ z . (S.37)Substituting the solution (S.31) in the expression for theenergy u and ignoring the surface term in ρ (as will bejustified later) we have u = (cid:90) d z √ (cid:20)(cid:18) 12 ∆ + F (cid:19) (cid:18) − ∂ Λ ∂z − ∂ F∂z (cid:19) + C.c (cid:21) + 18 π ∂ Λ ∂z ∂ Λ ∂ ¯ z . (S.38)Now integrating by parts and substituting (S.33) we ob-tain u = (cid:90) d z √ (cid:34)(cid:18) ∂F∂z (cid:19) + (cid:18) ∂ ¯ F∂ ¯ z (cid:19) − (cid:35) + 18 π ∂ Λ ∂z ∂ Λ ∂ ¯ z . (S.39)Finally substituting the leading order approximation forΛ (S.37) we obtain: u = (cid:90) d z √ (cid:18) ∂F∂z − ∂ ¯ F∂ ¯ z (cid:19) + 18 π ∂ Λ ∂z ∂ Λ ∂ ¯ z . (S.40)To proceed further we assume that the directors froma cylindrically symmetric configuration (21), where α ( r )is some arbitrary function (to be fixed by minimization)of the radius r = √ z ¯ z . Substituting in the third equationof (20) we obtain F ( z, ¯ z ) = − π √ D ∗ z (cid:90) Rr dr r exp[ iα ( r )] , (S.41)where we took into account that the force F vanishes onthe boundary r = R . From here it follows ∂F ( z, ¯ z ) ∂z = − π √ D ∗ r ddr (cid:40) r (cid:90) Rr dr r exp[ iα ( r )] (cid:41) . (S.42)Now to find α ( r ) one should substitute this expression in(S.40) and vary with respect to α ( r ). To first approxi-mation we can neglect the higher gradient term and takeonly the first term, namely u = − (2 π ) √ D ∗ (cid:90) R drr (cid:18) dY ( r ) dr (cid:19) , (S.43) where Y ( r ) = r (cid:90) Rr dr r sin[ α ( r )] . (S.44)Minima appear either for δYδα = 0, giving α = ± π , or for ddr (cid:18) r dYdr (cid:19) = 0 , (S.45)which gives Y ( r ) = R − r α ( R )] , (S.46)implying that sin[ α ( r )] = sin[ α ( R )] R r . (S.47)This solution clearly holds when r > R ∗ for some R ∗ < R and it should match the other solution at smaller valuesof r . Matching at r = R ∗ we have1 = | sin[ α ( R )] | R R ∗ . (S.48)In order to find R ∗ we should minimize the total energy ofthe system which includes also the surface energy due tothe interaction of the directors with the field V (2) i whichis very strong near the boundary: u s = −(cid:60) (cid:88) i ¯ V (2) i D i . (S.49)The orientation of the directors at the most outer rowsof the system is dictated by V (2) i , however this field atinner rows is much smaller [see Eqs. (S.12) and (S.13)].In order to take it into account we substitute V (2) i ≈ η ¯ z i z i where η (cid:28) D i , and weapproximate the sum (S.49) by an integral: u s = − η π √ D ∗ R cos[ α ( R )] . (S.50)To calculate the bulk energy we first perform the integral(S.44) to obtain Y ( r ) at r < R ∗ and find Y ( r ) = σr ln (cid:18) R ∗ r (cid:19) + r σ − sin[ α ( R )]) , (S.51)where σ = sign[sin[ α ( R )]] and therefore dY ( r ) dr = (cid:40) σr ln (cid:16) R ∗ r (cid:17) − rσ | sin[ α ( R )] | r < R ∗ − r sin[ α ( R )] r > R ∗ . (S.52)Substituting in u we obtain1 u = − (2 π ) √ D ∗ R (cid:20) sin α (cid:18) − (cid:90) dηη ln 1 η (cid:19) + 4 | sin[ α ( R )] | (cid:90) dηη ln η (cid:21) . (S.53)Both integrals in this expression equal 1 / u + u s = − (2 π ) √ D ∗ R (cid:20) | sin[ α ( R )] | − 12 sin [ α ( R )] (cid:21) − η π √ D ∗ R cos[ α ( R )] . (S.54)This expression has two minima within the range | α ( R ) | ≤ π . For large R , they are given by α ( R ) = ± π ∓ (2 γ ) / , where γ = 3 η/ ( πD ∗ R ). Substituting this result in Eq. (S.48) we find: R − R ∗ = 14 (2 γ ) / R = (cid:18) η πD ∗ (cid:19) / R / ..