Melting of the two-dimensional solid phase in the Gaussian-core model

A general theory for the melting of two dimensional solids explaining the universal and non-universal properties is an open problem up to date. Although the celebrated KTHNY theory have been able to predict the critical properties of the melting transition in a variety cases, it is already known that it is not able to capture the occurrence of first order transitions observed in certain systems as well as it doesn't provide a clear way to calculate the melting temperature for a specific model. In the present work we have developed an analytical method that combines Self Consistent Variational Approximation with the Renormalization Group in order to deal simultaneously with the phonon fluctuations and the topological defects present in the melting process of two dimensional crystals. The method was applied with impressive success to the study of the phase diagram of the Gaussian-core model, capturing not only the reentrant feature of its 2D solid phase, but also the related critical temperatures as a function of the density in quantitative detail. The developed method can be directly applied to study the melting of any hexagonal simple crystal formed by particles interacting through any finite pairwise interaction potential. Additionally, it has the potential to explain the occurrence of first order transitions in the melting process of two dimensional crystals.

A general theory for the melting of two dimensional solids explaining the universal and non-universal properties is an open problem up to date.Although the celebrated KTHNY theory have been able to predict the critical properties of the melting transition in a variety cases, it is already known that it is not able to capture the occurrence of first order transitions observed in certain systems as well as it doesn't provide a clear way to calculate the melting temperature for a specific model.In the present work we have developed an analytical method that combines Self Consistent Variational Approximation with the Renormalization Group in order to deal simultaneously with the phonon fluctuations and the topological defects present in the melting process of two dimensional crystals.The method was applied with impressive success to the study of the phase diagram of the Gaussian-core model, capturing not only the reentrant feature of its 2D solid phase, but also the related critical temperatures as a function of the density in quantitative detail.The developed method can be directly applied to study the melting of any hexagonal simple crystal formed by particles interacting through any finite pairwise interaction potential.Additionally, it has the potential to explain the occurrence of first order transitions in the melting process of two dimensional crystals.

I. INTRODUCTION
Two dimensional systems described by continuous microscopic variables present unusual phase transitions between the low and high temperature phases.This is related with the fact that in two dimensions thermal fluctuations are quite strong, preventing the stabilization of a long range ordered phase [1][2][3][4] .In this scenario, a vestigial quasi long-range order is often observed at low temperatures with slowly power law decaying correlations functions 5,6 .Such is the case of the nematic, 2D-solid and hexatic phases [7][8][9][10][11][12][13][14][15] in a vast class of different physical systems like ultra-thin magnetic films [16][17][18] , quasi two dimensional copolymer systems 19,20 ), vortex matter in two dimensional superconductors [21][22][23] and many others 13,14 .
The study of these kind of phases has its origin in a series of pioneering works by Berezinskii, Kosterlitz, Thouless (BKT), Halperin, Nelson and Young 5,6,[24][25][26][27][28][29] .In these works it was established the central role of the thermal induced process of pair unbinding and proliferation of topological defects in the disruption of the quasi long-range order.From the technical point of view the development of this research area is historically tied with the development of a proper Renormalization Group (RG) theory capable to describe the proliferation of the topological excitations at high-temperature.The success of this RG scheme granted the 2016 Nobel prize in physics 30 .
Despite this huge success, important questions remain open today in the physics of 2D transitions.Most of these questions originate from the difficulty to properly estimate non-universal properties, such as the transition temperature, in the wide range of diverse physical system, where topological defects unbinding occurs.Indeed, while the BKT RG scheme yields a phenomenological description of the coarse grained topological variables, the connection between their properties and the actual microscopic theory at hand is often unfeasible.Thus, a natural question arises "How can we maintain a close connection between the microscopic details of the system and its BKT description?".In the context of the melting of two-dimensional crystals, one can also ask "How can we construct a RG scheme which accounts for the variety of scenarios observed in simulations and experiments that include first order phase transitions [31][32][33][34][35] ?".The answers to these questions are the topic of the present work.
In order to access non-universal properties, first, it is necessary to build RG equations maintaining full connection with the microscopic model.Second, good estimates of the relevant energy of the defects are needed to capture the temperature scale of the corresponding phase transition 36,37 .And lastly, the effects of smooth fluctuations on the ground state should be taken into account, since close to the melting temperature they could produce a significant deviation of the effective microscopic rigidity from its bare or zero temperature value 6,17,38 .Following this route it has been possible to estimate the properties of the topological phase transition in diverse XY model configurations 37,[39][40][41][42] .In few specific cases the Self-Consistent Harmonic Approximation (SCHA) has been used to calculate the effective rigidity as a function of temperature and simultaneously this information used as an input of the RG equations in order to estimate the melting temperature better than RG or variational Mean Field alone 36,41 .
In this context, the construction of a similar calculation scheme for the study of the melting process of two dimensional crystals still constitutes an unexplored route.The construction of such method shall provide not only good estimates for the melting temperature of two dimensional crystals but shall also yield a self-consistent approach providing the correct qualitative behavior of the phases below and above the transition.This last property constitutes a substantial advantage with respect to the several implementations of Density Functional Theory (DFT), which have been used to tackle the melting problem in two and three dimensional crystals, [43][44][45] .These traditional approaches in the vast majority of cases consider that the 2D solid phase can be treated as a periodic phase which breaks translational symmetry, yielding a major shortcoming that in general produces a wrong description of the critical properties of the melting transition, or even predict incorrectly the nature of the phase transition itself.In light of the previous discussion, a theory of 2D melting, that we intend to construct, represents a major advancement towards the understanding of melting in two-dimensional crystals, as it will be capable to explain the diversity of melting scenarios observed in numerical simulations [31][32][33][34][46][47][48][49][50] .
In this work we implement a Self-Consistent Harmonic technique that retains total connection with the microscopic model.This allows us to determine the effective elastic Lamé's coefficients as a function of the density of particles and the temperature.A first mean field estimate of the melting temperature then can be found as the moment when the effective transversal elastic rigidity goes to zero.Contrasting with all other DFT techniques in two dimensions, this kind of mean field calculation has the distinctive virtue of being able to describe properly the qualitative behavior of the phases below and above the melting transition.To improve the mean field results, we use the obtained values for the effective elastic coefficients as an input into the RG equations for the melting transition of a two-dimensional solid.This produces a strong correction to the phase boundary of the 2D solid phase.To test the quality of this method that incorporates the SCHA into the RG theory we performed extensive overdamped Langevin simulations to determine an accurate melting curve for the two dimensional solid phase.The obtained analytical results shows an impressive agreement with the simulational results.

II. SELF-CONSISTENT HARMONIC APPROXIMATION FOR THE 2D SOLID MELTING
We consider a classical system of particles in two dimensions interacting through a Gaussian pairwise potential of the form V (r) = V 0 exp(−r 2 /r 2 0 ), known as the Gaussian-Core Model (GCM) [51][52][53] .The parameters r 0 and V 0 represent the range and the intensity of our soft-core potential, respectively.The partition function of this model can be written directly in terms of the configurational integral, which is obtained after integrating over all momenta of the particles: where Λ ≡ h/ √ 2πmk B T is the de Broglie thermal wave length and N the number of particles.Having in mind that our goal is to study the melting of the two-dimensional solid phase it is natural to consider that each particle will effectively explore a region of finite "volume" around its equilibrium position.This consideration allows us to write the position of each particle r i as R i + u i , where R i represent the equilibrium lattice position of the i-particle and u i its relative position respect R i .
For the pair interaction potential we are considering, the ground-state configuration is given by a triangular lattice with spacing a = (2/ √ 3ρ) 1/2 , where ρ stand for the density of particles in two dimensions.In this way the equilibrium positions of the particles can be represented as R j = a (ne 1 + me 2 ), with n and m integers and the basis vectors of the lattice taken as e 1 = (1, 0) and e 2 = (−1/2, √ 3/2).Once we include the combinatorial factor of distributing the N particles of the system on the N sites of the solid, the partition function can be rewritten as: This last expression allows us to define our effective Hamiltonian as Before we proceed to implement the Self-Consistent Harmonic Approximation (SCHA) in order to obtain the melting curve of the 2D solid phase, it is worth mentioning that the usual mean-field approaches devised as a variational theory on the local density profile fails to identify the 2D solid phase.Instead, depending on the pair interaction potential, such theory can predict at most a phase with broken translational symmetry in two dimensions -which is forbidden by the Mermim-Wagner theorem.Moreover, in our case, since the Fourier transform of our pair interaction potential is positive-definite, even such a theory would fail to identify the existence of the solid phase.Furthermore, the GCM have a density versus temperature reentrant hexagonal solid phase in its low temperature regime 52 and the melting of the solid phase occurs through a two-step process with an intermediate hexatic phase existing in a very narrow temperature region 53 .Now we proceed with the implementation of the SCHA.The Fourier transform of the field u(R j ) on the triangular lattice is defined as: where the momentum integral is performed over the first Brillouin zone (BZ) of the triangular lattice.We choose our test Hamiltonian in the usual harmonic form for an hexagonal two dimensional solid (see Appendix A): where ϵ 0 (ρ) represent the ground state energy per particle of the GCM at a given density ρ.The dispersion relations f 1 (q) and f 2 (q) represent the longitudinal and transversal elastic test response functions to be found, in principle, through the minimization of the variational free energy functional.Additionally, the longitudinal and transversal elastic fields are given as: û∥ (q) = −iq • û(q)/q and û⊥ (q) = −iq ⊥ • û(q)/q, with q ⊥ = (−q y , q x ).
Although the form of the test Hamiltonian may look highly elaborated it is in fact quite intuitive since it coincides with the expansion of H up to second order in powers of û(q), considering generic dispersion relations for the longitudinal and transversal elastic modes.In the long wave limit (q → 0), the harmonic theory predicts f 1 (q) = (2µ + λ)q 2 and f 2 (q) = µq 2 , where µ and λ represent the so-called Lamé's elastic coefficients -these and other results of harmonic theory for the solid elasticity are reviewed in the Appendix A.
Now we can proceed with the construction of the variational free energy.As it is well known, the actual free energy of the system is bounded from above by the minimum of the functional: where Given the Gaussian character of our test model it is not hard to conclude that: and consequently: Now we focus on the determination of the free energy of the test model F 0 , which can be obtained from the corresponding partition function Z 0 .In this case the quadratic form of the test Hamiltonian (5) allows for a direct integration of the Gaussian degrees of freedom, obtaining: which lead us to a free energy of the test Hamiltonian of the form: To proceed with the construction of our variational free energy, we should now determine ⟨H⟩ 0 .Considering the translational symmetry of the solid phase, it is possible to write: where {R i } corresponds to the equilibrium positions of the particles in a hexagonal solid at a fixed density.At the same time, the average interaction energy between two particles linked to lattice sites separated by a vector R is given by: where we have used the standard form of the Fourier transform in two dimensions to decouple the average process of the potential V (r).Here g q (R) ≡ e iq•(u(R)−u(0)) 0 , represent the positional correlation function calculated with the Boltzmann measure of H 0 .Considering now the Gaussian nature of the stochastic variables u(R), we can express g q (R) as: At this point we can verify that in the limit of zero temperature (β → ∞) the positional correlation function converges to one, as expected for the perfectly ordered crystal.
The expression in Eq. ( 13), once inserted in Eq. ( 12), complete the construction of the variational free energy functional in terms of f 1 (q) and f 2 (q).In this way, the variational free energy per particle in units of k B T can be written as: Demanding δf var /δf 1 (q 0 ) = 0 and δf var /δf 2 (q 0 ) = 0, we obtain the following set of integral equations for f 1 (q) and f 2 (q): where q 0 = (q 0,x , q 0,y ), q ⊥ 0 = (−q 0,y , q 0,x ) and g q (R j ) is defined in Eq. ( 13).The numerical solution of this set of integral equations is a quite difficult task.However, it is expected that the long distance elastic properties of the system is captured by low momentum behavior of f 1 (q) and f 2 (q).This can be used to follow a simpler approach to determine the boundary in the phase diagram of the 2D solid phase.Considering the symmetries of the solid phase under study, we already know that in the low momentum regime (q → 0) the leading order term of the functions f 1 (q) and f 2 (q) is proportional to q 2 .This feature can be used to derive a system of equations for the effective Lamé coefficients of the solid, which we define as: In this way, considering the form of Eq. ( 15), we can conclude that the elastic coefficients satisfy the following set of equations: where e can be taken as (1, 0) and e ⊥ as (0, 1).Once we have obtained an exact system of equations for r 1 and r 2 it is natural to approximate f 1 and f 2 by its low momentum form (f ) in the calculus of g q (R).This consideration lead us, after a simple but lengthy calculation, to: where q ∥ and q ⊥ represent the components of q parallel and perpendicular to R, and the coefficients ω ∥ (R) and ω ⊥ (R) depend on r 1 and r 2 in the following way: and the functions A ∥ (R) and A ⊥ (R) are defined as the following integrals over the first Brillouin zone of the hexagonal crystal: Here it is important to mention that, in the determination of g q (R), we have neglected the cross term proportional to q ∥ q ⊥ in the argument of the exponential function.This approximation is very well justified since, depending on the direction of R, the corresponding coefficient is either zero, or a small and rapidly decreasing function of R whose highest value is of the order of 1% of the other coefficients in the quadratic form.
In this way the obtained expression for g q (R) in Eqs. ( 18) and ( 19) allows us to close the system of Eqs. ( 17) for r 1 (T ) and r 2 (T ).Such system can now be used to determine the effective Lamé coefficients as they are affected by phonon fluctuations.We can then identify the 2D solid phase as the region in the phase diagram in which the effective shear modulus r 2 (T ) > 0.
A. Application to the Gaussian-Core Model (GCM) As mentioned earlier in this work we focus on the analytical study of the GCM for which V (q) = πV 0 e −q 2 r 2 0 /4 .The Gaussian form of this potential allows the analytical integration of the right-hand side of equations (17), leading us to following system of equations for r 1 (T ) and r 2 (T ): where the sum over R, as in eq. ( 17), is performed over the hexagonal lattice with spacing a(ρ), and ω ∥ and ω ⊥ represents the functions of R given by eq.(19).Finally, the coefficients A ∥ (R) and A ⊥ (R) entering in the system (21) are determine numerically for each value of R. Now we can proceed with the numerical solution of the system for r 1 and r 2 at any density and temperature.
In the first column of Fig. 1 we show the numerical solution for µ and 2µ + λ as a function of temperature at several fixed densities of particles.As we can see, the transversal elastic coefficient, or shear modulus, µ(T ) displays the typical behavior of the order parameter in a BKT-like transition with an abrupt decay when temperature approaches the melting temperature (T m ).In terms of densities, we are able to observe the re-entrant behavior of µ(T ) already expected from previous numerical simulations results 53 .In the other hand, the longitudinal elastic coefficient 2µ(T ) + λ(T ) grows steadily with the density, a expected behaviour for any stable solid.Additionally, we observe that this quantity is an increasing function of temperature in the region of densities where µ takes its maximum value, as can be seen in Fig. 1-b.This unusual behavior is related with the fact that in this region of densities the fluctuations in the position of particles tends to increase the effective repulsion between them.
In the second column of Fig.
(1), we show a comparison between the values of µ and 2µ + λ, at zero and at the meting temperature, as the density is increased.We observe that the value of µ decreases with temperatures at all densities, while (2µ + λ) displays a weak non-monotonical behavior with temperature.
Finally, the melting curve in the temperature versus density plane can be observed in Fig. (2) in green.A detailed analysis of this phase diagram for the GCM and its comparison with equivalent results using different techniques will be presented in the following sections.However, we can anticipate that direct comparison with simulation results lead us to conclude that this meanfield technique overestimates significantly the maximum melting temperature of the model.Nevertheless, the SCHA developed in this work has the merit of being a novel mean field calculation describing properly the ground-state properties of the system as well as the qualitative behavior of the 2D solid phase and its melting.

III. DEFECT MEDIATED PHASE TRANSITION AND RG RELATIONS
In order to improve the agreement between the mean field results and the computational results is necessary to take into account the effects of the relevant topological defects, known as dislocations.The proliferation of dislocations is quite effective in disrupting the periodic order of the solid phase as the melting process occurs increasing temperature.The theory for describing the defect mediated melting transition in two dimensions, also known as KTHNY theory, was provided in a series of foundational works by Toner, Halperin, Nelson and Young 26,28,29,54 in which RG equations are obtained for the renormalized Young's modulus and fugacity of the dislocations, respectively.
The RG system of equation for the Young's modulus in units of k B T (K(l)) and the defect's fugacity (y(l)) is given by: with the initial conditions given by the bare values of K and y, i.e.K(0) = 4µ(µ+λ) (2µ+λ)k B T and y(0) = exp(−E c /k B T ).Here µ and λ represent the Lamé's coefficients and E c represent the energy of an isolated defect.Within RG theory the melting temperature corresponds to the lowest temperature at which K(l → +∞) = 0, signaling that the long distance effective rigidity of the system goes to zero.
To use the RG equations for detecting the melting transition the energy of the relevant defects in the melting process should be provided.In our case we estimate such quantity using the harmonic elastic theory to calculate half of the energy corresponding to a pair of conjugate dislocations.The obtained result depends on the orientation of the Burguer's vector of the dislocations (e) with respect to the underlying lattice, and weakly on the orientation of the distance vector between the dislocations of the pair d.Since there is a continuum of possible configurations for the dislocation pair, the energy of the dislocation is estimated as the average energy between all configurations maintaining the minimum possible length of the dislocation pair.In this way the energy of a dislocation in units of k B T is estimated to be: An important consideration to reach this result is related to the minimum distance between a pair of stable dislocations.Numerical simulations performed for the GCM model lead us to the conclusion that such distance is higher than the lattice spacing.One useful way to obtain a pair of stable dislocations in this system is to subtract a particle from a perfect crystal configuration and leave the system to relax to a new stable configuration exhibiting a pair of dislocations separated by a distance that can be roughly estimated to d = √ 3a 55 .More details on how to reach this conclusion and calculate the energy of a dislocation pair can be found in Appendix B.
We realize now that, since the initial condition of the RG flow equations are completely determined by the value of the bare Young's modulus in units of thermal energy K(0), the melting transition will occur at a certain specific value of the parameter K(0).In this case the direct numerical solution of the system of Eqs.(22) allow us to conclude that the melting transition occurs approximately at K(0) c ≈ 23.922π.It is important to mention that within the KTHNY theory, it is well established that the effective Young's modulus at the critical point K(l → ∞) is equal to 16π.However the corresponding value of the bare Young's modulus K(0) it is not an universal quantity and its value at the critical point depends on the specific energy cost of the relevant defects.The value of K(0) c can now be used to build the melting curve considering the calculated Lamé's coefficients, µ(T ) and λ(T ), dressed by phonon fluctuations.The melting temperature (T m ) of the 2D solid phase at a given density is determined from the self-consistency relation 4µ(ρ, The result of this procedure is shown in Fig. 2, which can be compared with our Molecular Dynamics (MD) simulations results for the 2D solid melting temperature, discussed further in the next section.The agreement between the analytical and computational results is quite impressive, indicating that the proposed method not only captures well the phenomenology of the described phase transition but also produce a precise estimation of the melting transition.

IV. NUMERICAL SIMULATIONS
In order to confront the theoretical phase diagrams found in the previous sections, we have performed molecular dynamics (MD) simulations of the GCM in a NVT or canonical ensemble.It is worth mentioning that the analytical results of this work are derived in the same ensemble, while previous simulation results for this model are available in the NPT ensemble 53 .We are interested in sampling the equilibrium configurations, so we have employed the Langevin equation in the overdamped limit to simulate the behavior of the system in contact with a heat bath: where T is the temperature of the thermal bath and γ is the viscosity.We have chosen to measure the timescales in units of γ and the temperature scales in units of V 0 /k B .The random force ξ i (t) is a white noise with ⟨ξ i (t)⟩ = 0 and ⟨ξ i (t)ξ j (t ′ )⟩ = δ ij δ(t−t ′ ).
The simulations have been performed using the Heun algorithm to integrate numerically the N stochastic differential equations of motion 56 , for which we used a time step dt/γ = 0.1.In order to accelerate the relaxation to thermal equilibrium at all temperatures, we have used the parallel tempering technique 57,58 , considering sets of two types of initial conditions, one corresponding to a liquid (disordered) and another to a crystalline hexagonal lattice.This allows to define a criterion to identify thermal equilibrium at each temperature as the stage at which all sets of configurations attain the same stationary state.We have used a number of particles ranging from 1024 up to 8100 particles.
In order to characterize the phase diagram of the system, we measure two relevant and well established order parameters associated to the translational (ψ T ) and bond-orientational orders (ψ 6 ) 9 .They are defined as: where k 0 represents the characteristic wave vector of the reciprocal lattice, estimated from the position of the peaks in the structure factor.Additionally, N nn (j) stand for the number of nearest neighbours of the particle j and θ jl is the bond angle between particles j and l, both determined from a Delaunay triangulation.
The goal of our computational study is to address the construction of the density versus temperature phase diagram for the GCM.In particular, we focus in the estimation of the 2D solid melting temperature.Furthermore, it is well established in previous works that the solid-liquid melting in the GCM occurs through an intermediate hexatic phase, present in a quite slim region next to the 2D solid phase boundary.Also, in this work we estimate the phase boundaries of the hexatic phase.The susceptibility associated to the translational order parameter (χT ) in presented in magenta and the corresponding to the orientational order parameter (χ6) in green.The width of the filled curves represent the statistical uncertainty of the data.(c) Phase diagram of the GCM from MD simulations in NVT ensemble.The points in green are the estimates for the hexatic-liquid phase boundary and points in magenta are the analog for the 2D solid-hexatic phase boundary.The lines are guides to the eyes.The narrow region in grey corresponds to the hexatic phase.
The signatures of the transitions are more pronounced if we study the corresponding order parameter susceptibilities, defined as: In order to estimate the 2D solid melting temperature we observed the maximum of the translational order parameter suscep-tibility (χ T ), whereas for the hexatic melting temperature, the maximum of orientational order parameter susceptibility (χ 6 ).These quantities are shown as a function of temperature in the first row of Fig. 3 for an increasing number of particles and for two values of the density.The corresponding melting temperatures can be estimated by extrapolating linearly the location of the maxima as a function of N −1/2 .It is important to mention that considering the temperature corresponding to the maximum of the susceptibilities for the largest system size would produce a visually indistinguishable phase diagram in comparison with the one presented above.The resulting phase diagram is shown in Fig. 3.c, where we can notice that the extension of the hexatic phase is quite narrow in temperature and increases progressively with density from ρr 2 0 > 0.4.For the comparison between simulational and analytical results in Fig. 2 we have considered the magenta curve from Fig. 3.c, which corresponds to the boundary of the 2D solid phase.

A. Melting criteria from Young modulus measurements
One of the conclusions of our analytical results in the previous section is that the bare Young modulus K(0) in units of k B T at the melting temperature takes approximately the value of 23.9π for the GCM.This value is a consequence of the microscopic model constructed to estimate the core energy of the dislocations, shown in detail in Appendix B. Conversely, the K(l → ∞) = 16π criterion for the renormalized Young modulus is an universal prediction from the KTHNY theory valid for any model in principle.In this scenario, a microscopic calculation of the bare and renormalized Young modulus would be an important and independent validation of our analytical predictions.
The Young modulus can be measured in numerical simulation by well established techniques, such as the one that measures the elastic properties from unperturbed equilibrium distances between particles, shown in early works from Squire et al 59 .This method uses a quadratic expansion in the deformation field for the Helmholtz free energy density f , in order to extract the elastic tensor components C ijkl = ∂ 2 f /∂η ij ∂η kl | η=0 , where η ij represent a small strain tensor perturbation in the equilibrium positions of the particles.These elastic constants can be expressed as a sum of different contributions Here, the Born term is the same for both of these constants and can be written as: while the fluctuation terms are given by: where the stress tensor is defined as: where α, β are the Cartesian components of the distance r ij between pairs of particles.In order to estimate the Young's modulus, we define the shear modulus as G = C 44 − P and the bulk modulus as B = C 44 + C 12 , where the pressure P can be calculated using the virial expression.In this way, the Young's modulus can be found as: where a = (2/ √ 3ρ) 1/2 is the triangular lattice spacing.This quantity has been measured in two-dimensional melting simulation studies using the same technique depicted above [60][61][62] , as well as using other methods [63][64][65] .
The fluctuation contributions to the elastic constants given by Eqs. ( 29) and ( 30) are variances of the intensive variables σ αβ .These variances are expected to behave as 1/N , in order to produce finite values of C (F ) 44 and C (F ) 12 in the thermodynamic limit.Within the SCHA the variables r α ij in Eq. ( 31) have a Gaussian probability distribution.In general, such equation can be rewritten in terms of the Fourier transform of the pair potential V (r), which allow us to translate the nontrivial r ij -dependence of the kernel of σ αβ to an expression that contains r ij only in the argument of complex exponentials.The stochastic average of these kind of quantities decay exponentially with the quadratic fluctuation of r ij .In this way, the fluctuation contributions to the elastic constants given by Eqs. ( 29) and ( 30) can be seen as 'variance of variances' that, within the SCHA, are expect to be particularly small.Within this hypothesis we can identify the Eq. ( 32) calculated without the contributions from Eqs. ( 29) and (30) as K(0), the bare value of the Young's modulus in units of thermal energy as calculated within the SCHA.Consequently, we can think of this quantity measured in simulations as the effective Young's modulus affected only by phonon fluctuations.From this estimation of K(0) obtained within numerical simulations we can independently verify the melting criterion for the GCM obtained from analytical calculation (K(0) c ≃ 23.9π).At the same time, taking into account all contributions to C 44 and C 12 in Eq. ( 32) lead us to the macroscopic effective value of the Young's modulus K ≡ K(l → ∞) within our simulations.This provides a more standard path for the determination of the melting temperature of the 2D solid phase in numerical simulations, resulting from the criterion K(T c ) = 16π.
Finally, we present in Fig. 4.(a) the behavior of the bare Young's modulus K(0) together with the reference value of 23.9π for the melting temperature, and in Fig. 4.(b) the macroscopic Young's modulus K with its reference value of 16π.By selecting, at each density, the corresponding melting temperature we can draw the estimated 2D solid phase boundary corresponding to each melting criterion.The comparison of these phase boundaries with the one shown in Fig. 3.(c) is presented in Fig. 4.(c).As can be observed, the criterion K(0) c ≃ 23.9π for the GCM produces a 2D solid phase boundary that agrees well with the results obtained from more general methods.

V. CONCLUSIONS
In the present work we have developed a novel technique to study the 2D solid melting transition of the GCM.Although the boundary of this phase is already known from previous computational studies, up to our knowledge, this is the first time that it is reproduced analytically with such level of accuracy.In contrast with other MF predictions, we have observed that the SCHA alone is able to predict correctly the qualitative properties of the melting transition, such as the reentrant melting curve and the density for which the melting temperature is higher.In spite of this, only the mechanism of topological defects proliferation related to BKT like transitions is able to explain the significant reduction of the melting temperature when compared with MF predictions.This mechanism is captured in our approach by performing a two-step calculation where the Young's modulus and the defects fugacity calculated by MF are used as initial condition of the RG flow describing dislocation proliferation and unbinding.
Our approach demonstrates how the relative low value of the energy of dislocation pairs causes a drastic reduction of the melting temperature, due to the early proliferation of topological defects.The effect is correctly captured by the RG system in Eq. ( 22), once proper values of the defect parameters are introduced.Our estimate for the energy of the defects can be seen as a leading order approximation, since in general it is expected that higher order contributions will depend on the particle's density and on the Young's modulus itself in a nonlinear way.This is a possible explanation for the small difference observed in Fig. 2 between the analytical and numerical melting curves at higher densities.A better agreement could be reached by improving the estimation of the energy of the defects, however such involved calculation is beyond the scope of the present work.Additionally, it is important to mention that in the low density regime where the effects of the phonon are weak we verified that, as expected, the melting curves calculated using the zero temperature Young's modulus and using the effective one dressed by phonons produces similar melting curves.
We would like to stress that the presented method is an important step not only in the characterization of the melting transition of two 2D crystals but also to the study of this process in other modulated two dimensional systems, like for magnetic 2D textures [66][67][68][69][70] .Finally, it is worth noticing that the developed method has the potential to produce melting scenarios originally not contained in the KTHNY theory.A first order transition melting scenario is in principle a possibility in those cases where the SCHA produces a discontinuous melting transition and the defect energy is high enough for the mean field transition to takes place at lower temperature than the one predicted by RG equations.In this sense, this kind of self-consistent variational plus RG methods can possibly provide a more general framework to explain the variety of melting scenarios observed in different models.where A uc is taken as 1 2 √ 3a 2 , considering that the dislocations "center" have to be positioned on the underlying triangular lattice of the solid.The previous relation lead us to the conclusion that b c (r) = 1 2 √ 3a 2 eδ(r).In this way the field corresponding to a conjugate dislocation pair -as the one shown in Fig. 5  The modulus of the vector e is take it as one, which corresponds to the minimum possible value of the Burguer's vector of a dislocation in units of the lattice spacing.Now we can write our dislocation field in Fourier space as: It is not difficult to realize now that the energy of a dislocation pair will depend on the orientation of the vectors d, which is parameterized by the angle α, and e, parameterized by the angle θ, while both angles are measured respect to the x-axis.Consequently, we observe that the x and y components of the dislocation density field are given by: bx (q) = √ 3a 2 2 (−2i) cos(θ) sin q • d 2 , by (q) = √ 3a 2 2 (−2i) sin(θ) sin q • d 2 . (B6) In this way, the energy of a single dislocation can be calculated as one half of the total energy of a pair of conjugate dislocations, wich in units of k B T can be written as: 4π q 2 bx (q) bx (−q) 1 − q 2 x q 2 + 2 bx (q) by (−q) − q x q y q 2 + by (q) by (−q) 1 − q 2 y q 2 , (B7) which after some simplifications can be recast in the following form 12πa 4 (q ⊥ • e) 2 q 4 sin 2 (q • d/2), (q ⊥ • e) 2 q 4 sin 2 (q • d/2).

(B8)
To finally determine the actual value of the energy of a dislocation, we take the arithmetic average over all possible configurations in θ and α, considering the do not have in fact one single type of defect but a continuous class.This procedure allow us to conclude that: which lead us to the conclusion: In order to reach a final mean value for the energy of a dislocation, we need as an input the distance between the dislocations in the pair, d.We can consider this distance as the minimum distance for which a dislocation pair is stable, since such pairs are the most relevant.
To understand better how to estimate such a quantity we performed M.D. simulations at k B T ∼ 0 (see main text) for the specific model under consideration (GCM).We used the ordered triangular lattice as initial condition for our simulation, removed or added a single particle and study the stationary configuration of the system given those specific initial conditions.Here is important to remark that vacancy or interstitial defects can be both considered as dislocation pairs, nonetheless we often regard them as different kind of defects since they have a higher symmetry 55 than an arbitrary dislocation pair.In our simulations we have noticed that vacancies are stable for ρ ≤ 0.3, and naturally small fluctuations continuously deform its structure into pairs of dislocations with distance d ∼ √ 3a for ρ > 0.3 -both cases are represented in Fig. 5 (b).For higher temperatures, we have also noticed the presence of dislocation pairs more tightly bounded.Different from the defects found by relaxing the system from a vacancy initial condition, these pairs are not stable.They can be thus thought as "virtual" dislocation pairs, as they vanish quickly if we let the system relax at k B T ∼ 0. As an example of this line of reasoning, in ref. 55 , defects are considered as dislocation pairs only if the separation is greater than a certain threshold value necessarily higher than the lattice spacing of the underlying lattice.
Finally, if we consider d = √ 3a in Eq. (B10) and perform the numeric integration over the Brillouin zone, we can conclude our estimate for the energy of a dislocation, E c /k B T = bK, with b ≃ 0.072.

FIG. 1 :
FIG.1: Behaviour of the Lamé coeficients µ and 2µ + λ, resulting from the variational approach as function of the temperature, for some specific density values.The temperature ranges between zero and the critical value Tm, at which the µ(T ) function has a infinite slope.Comparison of the Lamé coefficients at T = 0 and T = Tm.The solid curves represent the ground state values and the dashed lines represent the value at T = Tm.

FIG. 2 :
FIG.2: Phase diagram for simulations are represented by black circles and dashed line (cubic splines used as guide to the eyes).The result for variational mean field is presented by a green line and renormalization group + variational approach in orange.

FIG. 3 :
FIG.3: Susceptibilities of the order parameters for increasing number of particles, for (a) ρr 2 0 = 0.4 and (b) ρr 2 0 = 0.7.The susceptibility associated to the translational order parameter (χT ) in presented in magenta and the corresponding to the orientational order parameter (χ6) in green.The width of the filled curves represent the statistical uncertainty of the data.(c) Phase diagram of the GCM from MD simulations in NVT ensemble.The points in green are the estimates for the hexatic-liquid phase boundary and points in magenta are the analog for the 2D solid-hexatic phase boundary.The lines are guides to the eyes.The narrow region in grey corresponds to the hexatic phase.

FIG. 4 :
FIG. 4: Simulation results for a system with N =1024 particles: (a) bare Young's modulus K(0) and (b) Young's modulus K measurements (see main text) for different densities as a function of the temperature.(c) Phase boundary for the 2D solid phase resulting from criteria K(Tc) = 16π and K(0)c ≃ 23.9π.For comparison the full line correspond to the phase boundary as shown in Fig. 3.(c).

FIG. 5 :
FIG.5: Vacancy and conjugate dislocation pair in a triangular lattice.The particles in green has five neighbours and in orange has seven.The arrows represent the burgers vector of each dislocation.