Vortex confinement transitions in the modified Goldstone model

Michikazu Kobayashi, ∗ Gergely Fejős, 3, 4, † Chandrasekhar Chatterjee, ‡ and Muneto Nitta 5, § Department of Physics, Kyoto University, Oiwake-cho, Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan Research and Education Center for Natural Sciences, Keio University, Hiyoshi 4-1-1, Yokohama, Kanagawa 223-8521, Japan Department of Physics, Keio University, Hiyoshi 3-14-1, Yokohama, Kanagawa, 223-8522, Japan Institute of Physics, Eötvös University, 1117 Budapest, Hungary Department of Physics, Keio University, Hiyoshi 4-1-1, Yokohama, Kanagawa 223-8521, Japan


I. INTRODUCTION
The Berezinskii-Kosterlitz-Thouless (BKT) transition [1][2][3][4] is a topological phase transition of two-dimensional systems, which divides a low temperature phase with bound vortex-antivortex pairs from a high temperature phase with free vortices. The phenomenon was first analyzed in terms of the XY model, and one of its most important impacts was that it showed that superfluidity and superconductivity can be realized even in twodimensions; Even though in two-dimensions long-range order with continuous symmetry is forbidden by the Coleman-Mermin-Wagner (CMW) theorem [5][6][7], there is a possibility of quasi-long-range order, which shows algebraically decaying correlations. The BKT transition realizes this scenario, and it also has the unique feature of being a continuous phase transition without breaking any symmetry. It has been experimentally confirmed in various condensed matter systems such as 4 He films [8], thin superconductors [9][10][11][12][13], Josephson-junction arrays [14,15], colloidal crystals [16][17][18], and ultracold atomic Bose gases [19]. The XY model shares common properties including the BKT transition with the twodimensional O(2), or Goldstone model at large distances or low energies, which is a regular version of the XY model described by one complex scalar field, in which the U(1) Goldstone mode for the XY model is complemented by a massive amplitude (Higgs) mode. One of the merits of the latter is to allow vortices as regular so- * michikaz@scphys.kyoto-u.ac.jp † fejos@keio.jp ‡ chandra.chttrj@gmail.com § nitta@phys-h.keio.ac.jp lutions in contrast to the XY model in which vortices are singular configurations.
XY-like models do not necessarily show the BKT transition. For example, for sharply increasing spin-spin potential, the phase transition between the paramagnetic and ferromagnetic phases can be of first order [20]. It is not of any surprise that the so-called modified XY model, where on a square lattice the Hamiltonian of the rotor is extended with a π periodic term, also shows a different scenario. It was predicted long ago that for large enough J coupling, there exists a nematic phase separated from the ferromagnet, and the transition between them is of Ising type [21,22]. This was also confirmed by numerical calculations [23]. The Ising-type transition is related to the presence of domain walls in this model. Moreover, it was conjectured that molecules and anti-molecules of half-quantized vortices play a crucial role for phase transitions, in contrast to a pair of vortices and anti-vortices in the XY model. As of today, the model (1) and its various modifications [24][25][26][27][28][29][30][31][32] are of great importance and interest, especially due to their relevance in condensed matter physics applications, e.g. superfluidity in atomic Bose gases [33], arrays of unconventional Josephson junctions [34], or high temperature superconductivity [35].
The BKT transition of the XY model was originally analyzed via a real-space renormalization group (RG) approach [4], which is rather unconventional and not easily linkable to the Wilsonian picture of the RG [36]. In the past years, the functional RG (FRG) approach, which adopts the Wilsonian idea of mode elimination and averaging to the level of the effective action [37], has also been applied and developed in regard to the BKT transition both in continuum [38][39][40][41][42] and lattice formulations [43,44]. It turned out that the method was capable of displaying in the two-dimensional O(2), or Goldstone model, a line of fixed points that is responsible for the topological nature of the phase transition. This is remarkable in the sense that no vortices need to be introduced explicitly, as opposed to the conventional real space RG description [4]. One of the shortcomings of the FRG treatment, however, is that because one is typically solving the RG flow equation of the scale dependent effective average action via a derivative expansion, as an artifact, only a line of quasi-fixed points are found. That is, the RG flow does not stop along this line, but only slows down significantly compared to other regions of the parameter space. There have been attempts to improve on obtaining a true line of fixed points, but in the continuum theory no scheme has been developed yet, which could have been successful.
The goal of this study is twofold. On one hand, we aim to show an approximation scheme of the FRG flow equations that can show significant improvement on the possibility of reaching a true line of fixed points in the continuum version of the XY model, and which can also be applied naturally to the modified XY model, i.e. the continuum version of (1). To our knowledge, for the first time in the framework of a momentum space RG, we describe the two-step transition in the latter model, and we will also predict that fluctuations may completely make the topological transition disappear. On the other hand, we also aim to provide full numerical simulation of the system, and show that depending on the value of the self-coupling of the scalar field, the structure of the transitions is even richer than it is predicted by the RG.
The paper is organized as follows. In Sec. II, we introduce the modified Goldstone model and construct classical solutions, an integer vortex, a soliton, and a vortex molecule of two half-integer vortices connected by a soliton in that model. In Sec. III, after giving a brief review of the FRG, we reproduce some earlier results of the BKT transition via the FRG, and also show the improvement announced above. Then, this scheme is applied to the modified XY model, and we show how a two-step transition can emerge in the system. In Sec. IV, we confirm this scenario via full numerical simulations, and reveal what nature of the corresponding transitions have. Sec. V is devoted for a summary. In Appendix A, we show how to derive the Hamiltonian of the modified Goldstone model from the microscopic lattice model of the modified XY model, while in Appendix B, we derive some of the corresponding flow equations of the FRG.

II. MODEL AND SOLUTIONS
A. Modified Goldstone model In this study we are interested in the continuum version of the XY model, i.e. the Goldstone model and its modification [for its derivation from the microscopic Hamiltonian (1) see Appendix A]: where ψ is a complex scalar field, and λ, a and b are positive coupling constants. The continuum version of the standard XY model refers to b = 0 and in the modified XY model we have b > 0. The field equation can be obtained from the Hamiltonian (2) as that we call the modified Gross-Piteavskii equation.

B. Classical solutions
Field equations (3) of the modified Goldstone model admit superfluid (or global) vortex solutions. Here, we show how such a vortex solution transforms into a halfquantized vortex molecule, when the second term of Eq. (2) becomes large enough. We work in a simplified parameter space, where a 2 + b 2 = 1, and thus the a = cos θ, b = sin θ parametrization can be used. The transformation of the solution can be seen in Fig. 1. One observes that around θ ≈ 78 • , a clear picture of a vortex molecule emerges, where two half-quantized vortices are connected by a one-dimensional soliton. One expects that at finite temperature, as a function of θ, somewhere close to the aforementioned value, emergence of the molecules will have an effect on the phase structure of the system.
In a vortex molecule shown in Fig. 1, each of the two vortices has a half-quantized circulation d l · (∇Arg[ψ]) = π, and the soliton connecting them has a π-phase jump. To analyze the stability of the soliton, we determine the following one-dimensional stable solution of the modified Gross-Piteavskii equation (3) in 1D with the boundary condition ψ(x → −∞) = √ 2ρ 0 , and ψ(x → ∞) = √ 2ρ 0 e iϕ . Fig. 2 (a) shows profiles of soliton solutions. Fig. 2 (b) shows the total energy H 1D as a function of ϕ and θ. It is clear that if H 1D takes the maximum value at some ϕ < π, then the soliton solution with ϕ = π becomes locally stable (metastable) by having a positive energy barrier ∆H 1D ≡ H 1D (ϕ = (a) where the maximal angle ϕ max is the value of ϕ at which H 1D takes the maximum. Fig. 2 (c) shows the maximal angle ϕ max and the energy barrier ∆H 1D . The former starts to take the nonzero value, ∆H 1D > 0, with ϕ max < π at around θ ≈ 15 • , above which the soliton is, therefore, energetically stable.
That is to say that the appearance of vortex molecules and the stability of the soliton are not related, thus it is not the (de)stabilization of the domain wall that lets molecules emerge.
It is worth to note that these configurations become singular in the limit of λ → ∞, in which the model reduces to the modified XY model. Therefore, the modified XY model does not allow these configurations as solutions to the field equations while the modified Goldstone model does.

C. Type of symmetry and (quasi) breaking of symmetry
Here, we discuss the symmetry properties of the Hamiltonian [Eq. (2)] and show the possible (quasi-)breaking patterns of symmetries. The symmetry of the Hamiltonian with generic parameters is of U(1) as a phase shift of the field, ψ → ψe iα for the arbitrary α ∈ [0, 2π). In the case of a = 0 and b > 0, the two fields ψ and ψ e iπ are identifiable, because the Hamiltonian [Eq. (2)] is the functional of ψ 2 , rather than ψ. Therefore, the symmetry of the Hamiltonian is only U(1)/Z 2 , where the Z 2 symmetry comes from the identification of ψ ∼ ψ e iπ . This Z 2 factor is essential for the presence of (deconfined) halfquantized vortices.
Depending on the parameter regions, the U(1) or U(1)/Z 2 symmetry is spontaneously broken in the ground state in different patterns summarized as follows: for a > 0 and b = 0, (4a) Here, arrows , −→, and =⇒ denote quasi-breaking of symmetry via a BKT transition, ordinary symmetry breaking with a thermodynamic phase transition, and simultaneous (quasi-)breakings of symmetry, respectively. Here, "quasi" breaking means that the symmetry is not exactly broken due to the CMW theorem in the thermodynamic limit but is locally broken at semi-macroscopic scales with an algebraically decaying correlation function. Now let us explain each breaking pattern. In the simplest case, i.e., for a > 0 and b = 0 [Eq. (4a)], the standard BKT transition occurs with the quasi breaking of the U(1) symmetry. In the opposite case, i.e., for a = 0 and b > 0 [Eq. (4b)], the BKT transition occurs with the quasi breaking of the U(1)/Z 2 symmetry, for which half-quantized and anti half-quantized vortices start to form in pairs. In the case of b a > 0 [Eq. (4c)], two successive spontaneous (quasi-)breakings occurs. At the first stage (at higher temperature) the U(1) symmetry is quasi broken to a Z 2 subgroup accompanied with the BKT transition. At the second stage, in the temperature lower than the BKT transition temperature, the remaining Z 2 symmetry is further spontaneously broken due to the thermodynamic transition. In this case, halfquantized and anti half-quantized vortices start to form pairs at the BKT transition, and domain walls appear at the thermodynamic transition. Some domain walls have no endpoint forming loops as well as those in the Ising model, but some others appear between two halfquantized or two anti half-quantized vortices forming vortex or anti-vortex molecules as shown in Fig. 1. In the remaining case of a ≈ b [Eq. (4d)], rather than a conventional BKT transition, the BKT transition occurs with the quasi breaking of U(1)/Z 2 symmetry, and the thermodynamic transition with breaking of Z 2 symmetry simultaneously. All vortices are integer and domain walls do not have endpoints.
In the following sections, we study the modified Goldstone model by the FRG and Monte-Carlo simulation.

III. FUNCTIONAL RENORMALIZATION GROUP CALCULATIONS
In this section, after giving a brief review of FRG, we apply it to the modified Goldstone model approximately, at the leading order of the derivative expansion, and obtain the phase structure.

A. Flow equation: a review
Here we review the basics of the FRG. In the core of the formalism lies the Γ k average effective action, in which fluctuations of the dynamical fields are incorporated up to a momentum scale k. The Γ k function obeys the following flow equation: where Γ (2) k is the second derivative matrix of Γ k with respect to the dynamical variables, R k is a regulator function, which is defined (in Fourier space) through a momentum dependent mass term: added to the classical Hamiltonian (or Euclidean action). We denoted the set of fluctuating field variables by ψ. R k is supposed to give a large mass to modes that have momenta q k, and leave the ones with q k untouched. The classical Hamiltonian is by definition does not contain any fluctuations, therefore, it serves as an initial condition for the RG flow of Γ k=Λ at some microscopic scale Λ. The flow equation (5) is then need to be integrated down to k = 0, where one obtains the full free energy (or quantum effective action). One is free to choose the R k function such that it fulfills the requirement of suppressing low momentum modes, and in this paper we employ the so-called optimal version: where Θ(x) is the Heaviside step function, and Z k is the wavefunction renormalization factor.

B. Local potential approximation'
Here, we solve flow equation (5) for the modified Goldstone model approximately, at the leading order of the derivative expansion. It is sometimes called the Local Potential Approximation' (LPA') ansatz, with the prime referring to nontrivial wavefunction renormalizations. In our approximation, Γ k takes the following form: where instead of a complex variable, the ψ i field is considered as a two-component real vector: ψ i = (ψ 1 , ψ 2 ), while ρ = ψ i ψ i /2, and we have only kept the original couplings in the effective potential. Namely, Eq. (8) is compatible with form of Eq. (2), but it comes with kdependent couplings and field dependent wavefunction renormalization factor [Z k (ρ)]. Projecting the flow equation (5) onto a subspace spanned by homogeneous field configurations, one gets (see also Appendix B) is the anomalous dimension at this order of the approximation, where the wavefunction renormalization is evaluated at the minimum point of the effective potential, where the index t refers to the transverse direction, we arrive at the flow equation for Z k (ρ) (see, again, Appendix B for details): and M 2 t,k are the longitudinal and transverse components of the momentum independent part of the Γ (2) k matrix, respectively: k , we evaluate Eq. (10) atρ =ρ 0,k and get Now one can search for fixed points of Eqs. (9) and (13). The flow diagram in terms ofλ k andρ 0,k can be seen in the left side of Fig. 3. One observes the line of quasi-fixed points, and notes that the flow, even though significantly slows down, is clearly nonzero in the aforementioned region.

C. Wavefunction renormalization improvement
The key of the improvement to be described here is to realize how crucial the role of the wavefunction renormalization factor (Z k ) is in the previous description. In order to escape from the CMW theorem, in the low temperature phase Z k has to diverge so that the renormalized field can condense (the expectation value of the bare field is always zero). Since any rescaling of the field should lead to the same description of the system, one expects that any field derivative of the wavefunction renormalization factor is proportional to Z k itself, Z (n) k ∼ Z k , which means that they also diverge, and in principle none of them should be neglected. Here we take into account the first derivative of Z k , which will indeed lead to a significant improvement, and it also makes possible to treat the modified version of the XY model in the FRG.
If we keep track of the field derivative of Z k , then it is possible to take into account in η k the implicit kdependence coming from the change of the minimum of the effective potential when the RG scale is varied. To our knowledge this effect has always been neglected in earlier computations. In principle, one should have where d k refers to total differentiation, and in the right hand side both Z k and Z k are evaluated atρ =ρ 0,k . We have also introduced the notation ∆η k = −w k k∂ kρ0,k with w k = Z k /Z k . Since Z k has appeared in our formula, we also need to derive a flow equation for it. This can be obtained from Eq. (10) after applying d/dρ to the both sides (note that ∂ k does not commute with d/dρ). Using that dM 2 l,k /dρ = 3λ k , dM 2 t,k /dρ =λ k , we get which leads to At this point, it is important to mention that Eq. (15) is not exact, as deriving Eq. (10) one lets the field operators act only on the potential part of the two-point correlation function and not on Z k (ρ). This would have introduced further Z k (ρ) dependence in the right hand side of Eq. (10), which is neglected here. Now, going back to Eqs. (9), we notice that the flow equations change due to the new approximation of the anomalous dimension, η This does not make much of a difference in the flow ofλ k , but completely changes the flow ofρ 0,k . The reason is that Eq. (9b) becomes an implicit equation, since k∂ kρ0,k also appears in the right hand side through ∆η k ≡ −w k k∂ kρ0,k . After some algebra we arrive at The flow ofλ k is analogous to Eq. (9a), but η (0) k is replaced by η k : Now one solves the coupled equations (13), (16), (17) and (18). The corresponding flow diagram can be seen in the left side of Fig. 3. The comparison shows that taking into account the derivative of the wavefunction renormalization factor in the anomalous dimension significantly stabilizes the flow along the line of (quasi-)fixed points, as in the improved case the freezing of the flow holds on ∼ 20 times longer in RG time t = − log(k/Λ).

D. Phase structure
Now we are in a position to show that in the modified XY model fluctuations can dramatically change the structure of the line of fixed points, as seen in Fig. 3. First, note that, the ansatz of Eq. (8) and the approximation Z k (ρ) ≈ Z k (ρ 0,k ) + Z k (ρ 0,k )(ρ −ρ 0,k ) is compatible with the microscopic Hamiltonian of the modified XY model, since from Eq. (8) we have which is equivalent to where a k = (Z k − Z 2 k w k ρ 0,k )/2, b k = Z 2 k w k /16. Eq. (20) is now of the form of the original Hamiltonian in Eq. (2) using the ψ i vector notation.
The reason why the RG flows of the ordinary XY model can change dramatically is that depending on the initial value w Λ (or b Λ , equivalently) at the UV scale,ρ 0,k can approach a singularity, which sends the flows in thē λ k −ρ k plane away from the line of fixed points. What essentially happens is that the line of quasi-fixed points terminates also at another endpoint, see Fig. 4. The endpoint on the left corresponds to a BKT transition at higher temperature, and the new one on the right signals another transition at lower temperature. Even though the method does not make a definite prediction, this should correspond to the Ising transition already reported in the earlier papers [21][22][23].
Analyzing the flow ofρ k , one notes that already in the ordinary XY model (i.e. for w Λ = 0) at first sight it might seem to be possible that the denominator in the right hand side of Eq. (17) becomes zero, but it turns out that it never happens. The flow equation always makes w k decrease as fluctuations are integrated out, therefore, the flows are regular. Note that, however, if, at the microscopic scale, w Λ > 0, then k∂ kρ0,k can indeed blow up.
The condition that needs to be met for a diverging flow is which shows that for positive w Λ values the line of fixed points can also terminate on the right (see Fig. 4), leading to a two-step transition. For later reference, just as in Sec. II, we restrict ourselves to the following case: i.e. we may use the parametrization a Λ = cos θ, b Λ = sin θ (θ ∈ [0, π/2]), which leads to the following constraints: Solving them for w Λ and Z Λ we get Dropping the last term in the bracket of the right hand side of Eq. (21) (we are interested in a rough estimate) we can get the following condition for the critical value of ρ 0,Λ belonging to the second endpoint of the line of fixed points: 0 = − (cos θ + 8ρ 0,Λ sin θ) 2 4 sin θ + 2(cos θ + 8ρ 0,Λ sin θ)ρ 0,Λ + 1 16π .
For a given θ, this is an equation for the second endpoint of the line of fixed points, in terms of ρ 0,Λ , see Fig. 4. Surprisingly, if θ = 0 is small (i.e. we are close to the XY model), the solution ρ 0,Λ | sol is always negative. This means that since the flows blow up for initial values ρ 0,Λ > ρ 0,Λ | sol , unlessρ 0,Λ | sol ≡ Z Λ ρ 0,Λ | sol ≈ 0.5 (which is the location of the original endpoint of the BKT transition), the line of fixed points completely disappear.
The critical angle at which this happens is That is to say, for 0 = θ < θ c , if there is a transition in the system, it cannot be of topological type, no matter how close we are to the XY model (still, at θ = 0 we have one, and only one BKT transition). However, once θ > θ c , the line of fixed points starts to come back to the picture, now equipped with another endpoint, which indicates that there exist two transitions. A higher temperature transition has to be of BKT type, and a lower temperature transition, presumably an Ising transition [23], is expected to be of second order (though the corresponding fixed point is not seen in this approximation). Note that, the aforementioned structure heavily relies on the assumption a 2 Λ + b 2 Λ = 1. Had we not had this constraint, and just set e.g. a Λ ≡ 1, we would have found a two step transition for 0 < b < b c (the higher temperature one being topological), and no topological transition for b > b c (here b c > 0 is some positive constant).

IV. NUMERICAL SIMULATIONS
In this section, we numerically investigate the equilibrium properties of the modified Goldstone model defined in Eq. (2).

A. Preparation
The discretized Hamiltonian, H ∆x , from Eq. (2), becomes where ψ i is the field ψ at the discretized point x = x i and ∆x is the lattice spacing (which serves as an ultraviolet cutoff scale). In the limit of λ → ∞ and rewriting ψ = √ 2ρ 0 e iθi , the discretized Hamiltonian H ∆x becomes equivalent to the Hamiltonian H mXY in Eq. (1) for the modified XY model with J = 4aρ 0 and J = 8bρ 2 0 . Now we numerically calculate equilibrium ensemble averages

B. Correlation function and transition temperature
We first show our results for the following two correlation functions: where L is the system size and N (r) is the number of points, x i , that satisfy r ≤ |x i | < r + ∆x. When θ = 0 (θ = π/2), we expect the standard BKT transition triggered by integer vortices (half-quantized vortices) for ψ i (ψ 2 i ), and the algebraic decay G 1 (r) ∝ r −η (G 2 (r) ∝ r −η ) below the BKT transition temperature. At the BKT transition temperature, the critical exponent satisfies η = 1/4 [3,4]. To obtain the BKT transition temperature, therefore, we can use the finite-size scaling of the correlation functions, in which G (1,2) /r −1/4 is expected to be a universal function of r/L. shows the dependence of G 1 (r)/L −1/4 with θ = 0 and λ = 8 as a function of r/L at T = 0.6T * where T * is the BKT transition temperature for the standard XYmodel with θ = 0 and λ → ∞. The expected universality of G 1 (r) is sufficiently satisfied at large r, which, therefore, predicts that the BKT transition temperature is T BKT 1 0.6T * . In the same way, we can estimate the temperature T BKT 2 0.21T * with θ = π/2 and λ = 8 from the finite-size scaling of G 2 (r).
We further expect the appearance of a second order, Ising-type phase transition [23], where the domain of definition for the phase of the ψ field is spontaneously broken from [0,2π] to [0,π], which can be thought of as a spontaneous breaking of a discrete Z 2 symmetry. At the critical temperature for this phase transition, the correlation function also shows algebraic decay. Since the critical exponent η takes the same value as that of the BKT transition temperature, i.e., η = 1/4 for the twodimensional Ising-type transition, we can use the same finite-size scaling analysis as shown in Fig. 5. We here define the temperature T 1 (T 2 ) at which G 1 (r) (G 2 (r)) shows the algebraic decay G 1 ∝ r −1/4 (G 2 (r) ∝ r −1/4 ). Then, by definition, T 1 = T BKT 1 at θ = 0 and T 2 = T BKT 2 at θ = π/2. Denoting by θ 1 and θ 2 critical angles, we have found the following results for T 1 and T 2 : 1. When θ is small, i.e., θ ≤ θ 1 , then T 1 > T 2 .
4. When λ → ∞ for the modified XY-model, then θ 1 = θ 2 , i.e., both T 1 and T 2 always exist at any θ. The specific values of θ 1 and θ 2 are shown in TABLE I.

C. Superfluid density and specific heat
To determine the type of the transitions, we calculate the superfluid density ρ s defined as [46,47] and the specific heat C = d H /dT , where F (δ) = −T log e −H/T is the free energy under the argumenttwisted boundary condition ψ(x + L) = e iδ·L ψ(x). When a BKT transition occurs at the transition temperature T BKT , the universal jump, ∆ρ s , of the superfluid density is On the other hand, for second order transitions we expect close to the corresponding critical temperature (T 2nd ) that, the superfluid density obeys ρ s ∝ (T 2nd − T ) ζ . The critical exponent ζ is obtained by the Josephson relation ζ = 2β − νη, where β, ν, and η are the critical exponents of the order parameter, the correlation length, and the correlation function, respectively. By inserting β = 1/8, ν = 1, and η = 1/4 for the Ising-type transition, we obtain ζ = 0, i.e., the superfluid density also jumps at the transition temperature, similarly to the BKT transition. However, no universal relation holds, which allows for a distinction between the two. Fig. 6 shows the dependence of the superfluid density with respect to the temperature for θ = 0 • [the panel (a)] and θ = 10 • [the panel (b)]. The solid line shows the relation ρ s = T /π. In the panel (a), this line intersects ρ s with a good accuracy at T 1 suggesting the standard universal relation related to the BKT transition temperature, i.e., we indeed observe a topological transition. In the panel (b), however, ρ s deviates from the aforementioned line at T 1 , and therefore, we expect that the transition is of second order, with a nonuniversal jump at the transition temperature. Here, we relabel T 1 ≡ T 2nd 1 . In neither of the panels do we find any characteristic structure in ρ s at T = T 2 . We, therefore, conclude that the property of the correlation function G 2 ∝ r −1/4 is just the crossover, and we relabel T 2 as the crossover temperature T 2 ≡ T CO 2 .
(a)  ]. As shown in TABLE I, the value θ = 60 • is between θ 1 and θ 2 for λ = 8, and we find neither a BKT, nor a second order phase transition. Instead, what we see is a first order phase transition due to the sharp jump of the superfluid density ρ s , see Fig. 7 (a). Because the temperature at which the superfluid density ρ s jumps does not really depend on the system size L, its estimation is fairly simple. We denote this transition temperature by T 1st * . In Fig. 7 (b), i.e., for θ = 87 • , θ is larger than θ 2 , and the superfluid density ρ s does show the universal relation (34) at the corresponding temperature, T = T 2 . Therefore, we find, again, a BKT transition with the aforementioned transition temperature, relabeling it as T 2 ≡ T BKT 2 . Fig. 8 shows the jump of the superfluid density ∆ρ s at the phase transition as a function of θ, normalized by ∆ρ s0 , which is the value for the universal jump (34) for (a) the BKT transition. It is specifically defined as (note that T 1 , T 1st * and T 2 depend on θ) We estimate the value of the jump ∆ρ s by fitting the superfluid density ρ s at the transition temperature (i.e. T 1 for 0 ≤ θ < θ 1 , T 1st * for θ 1 ≤ θ < θ 2 , and T 2 for θ 2 ≤ θ ≤ π/2) via the function where a is a θ dependent constant. For θ = 0 and θ > θ 2 , the relation ∆ρ s ∆ρ s0 is satisfied, therefore, we find BKT transitions with the transition temperature T BKT 1 for θ = 0 and T BKT 2 for θ 1 ≤ θ ≤ π/2. For other values, the universal relation does not hold and the transition becomes of second order for 0 < θ < θ 1 , and of first order for θ 1 < θ ≤ θ 2 . Fig. 9 shows the specific heat C. Whereas the specific heat has a single peak near the transition temperature for θ < θ 2 , i.e., in the panels (a)-(c), it has double peaks for θ ≥ θ 2 , suggesting two-step transitions. In the latter case, the first and second peaks of the specific heat correspond to the temperatures T 1 and T 2 , respectively. Because the correlation function G 1 becomes G 1 ∝ r −1/4 at T = T 1 and the phase at T < T 1 should be continuously connected from the phase with θ < θ 2 (see Fig. 10), the transition at T 1 should indeed be of second order. The absence of the peak at T = T 2 for θ < θ 1 consolidates our conclusion about that here T 2 gives not the transition but only a crossover as T CO 2 .
(a)  , integer vortex pairs are bounded to show a quasi long-range ordered phase. For 0 < θ < θ 1 , this BKT transition changes to a second order phase transition with the transition temperature T 1 ≡ T 2nd 1 , implying a true long-range ordered phase for T < T 2nd 1 with the breaking of the Z 2 symmetry. For θ 1 ≤ θ < θ 2 , the two temperatures T 2nd 1 and T CO 2 defined for 0 < θ < θ 1 merges to one first order transition temperature, T 1st * . For θ 2 ≤ θ ≤ π/2, this transition temperature, T 1st * , splits again into two transition temperatures, T 2nd 1 and T BKT 2 . The second order phase transition ultimately disappears, as while θ → π/2, T 2nd 1 → 0. Unlike the BKT transition for θ = 0, the BKT transition for θ 2 ≤ θ ≤ π/2 is triggered by the correlation function G 2 (not G 1 ), and therefore, we expect the quasi long-range ordered phase by the bounding of half-quantized vortex pairs at T 2nd 1 < T < T BKT 2 . Because the low temperature phases (i.e., T < T 1st * for θ 1 ≤ θ < θ 2 and T < T 2nd 1 for θ 2 ≤ θ ≤ π/2) should be continuously connected from the long-range ordered phase at 0 < θ < θ 1 , these phases should also be of true long-range order.
Here, we wish to establish the relationship between the phase diagram and the (quasi-)breaking patterns of symmetry summarized in Eqs.  Fig. 8 and Fig. 10, respectively.
Finally, in Fig. 11, we show the jump of the superfluid density ∆ρ s and the phase diagram in the λ = ∞ limit, in which the modified Goldstone model reduces to the modified XY model. As the coupling λ increases, the region of the first order phase transition for θ 1 < θ ≤ θ 2 shrinks and ultimately disappears.

E. Vortex configurations
Here we discuss the relationship between topological defects (such as integer and half-integer vortices, and one-dimensional solitons considered in Sec. II) and the corresponding phase transitions. At the BKT transition temperature T BKT 1 with θ = 0 • , the number of integer vortex-antivortex pairs is changing rapidly due to their bounding. At the second and first order transition temperatures, T 2nd 2 and T 1st 1 , the Z 2 symmetry breaking causes the rapid decrease of one-dimensional solitons. At the BKT transition temperature T BKT 2 with θ > θ 2 , the number of half-integer vortex-antivortex pairs changes rapidly. The vortex-molecules, which contain two halfquantized vortices should be stable in order for the BKT transition to exist at the temperature T BKT 2 . On the other hand, the stability of one-dimensional solitons is enough for the existence of the Z 2 symmetry breaking. The stability of vortex molecules for θ 78 • and onedimensional solitons for θ 15 • in the case of λ = 8 are consistent with the existence of T BKT 2 for θ > θ 2 ≈ 84.5 • , and the Z 2 symmetry breaking at T 2nd 1 or T 1st * for θ > 0. We next show snapshots of vortex configurations and the phase profile at the transition temperatures in Fig. 12. In all the panels, most vortices and antivortices form paired states with short distances. Furthermore, most of them lie on the solitons that appear as boundaries between the two phases Arg[ψ] ∼ 0 and Arg[ψ] ∼ π. Fig. 13 shows the distribution function P (Arg[ψ]) corresponding to the snapshot of the phase profile. In Fig. 13 (b) for θ = 87 • , the stability of one-dimensional solitons can be clearly seen from the double peaked structure of P (Arg[ψ]) at Arg[ψ] = 0 and Arg[ψ] = π. At T = T BKT 2 , the Z 2 symmetry is not broken and the height of two peaks are the same. On the other hand, breaking of the Z 2 symmetry at T = T 2nd 1 can be confirmed via the existence of imbalanced peaks, P (0) > P (π). This imbalanced distribution can also be seen in Fig. 12 (d), where the region with Arg[ψ] ∼ 0 is apparently larger than that with Arg[ψ] = π and Arg[ψ] = −π. Note that, in Fig. 13 (a) for θ = 10 • and θ = 60 • , however, the double peaked structure is absent and there is only one single peak at Arg[ψ] = 0. We believe that this absence comes from finite-size effects and it is expected that the doublepeaked structure is restored with larger system size. We note that all the peaked structures shown in Figs. 13(a) and 13(b) come from finite-size effects, and they become completely flat in the thermodynamic limit due to the CMW theorem.

V. SUMMARY
In this paper, we first have defined the modified Goldstone model in Eq. (2) as a regular and continuum version of the modified XY model, and constructed a soliton, an integer vortex and a molecule of half-quantized vortices connected by a soliton. Then, we have analyzed the phase structure of the modified Goldstone model in two dimensions via two different approaches. First, by using the functional renormalization group technique, we have shown how to describe BKT transitions by calculating the scale evolution of the effective Hamiltonian. Based on earlier works we have constructed a new approximation scheme of the RG flow equations, where the field dependence of the wavefunction renormalization is taken into account. In the standard Goldstone model it has led to a more accurate description of the underlying structure of line of fixed points, and it has also turned out to be of particular importance when one is interested in the role of the modified kinetic term ∼ |∇ψ 2 | 2 , by revealing a second endpoint of the line of fixed points. The FRG method predicts that in the modified model there can exist a two-step phase transition, depending on the ratio between the coefficients of the standard and modified kinetic terms. It has also been shown that even if the coefficient of the modified kinetic term is not large enough to split the phase transition into two, it is capable of completely destroying its topological nature.
Second, this scenario has been verified to great accuracy via full numerical simulation of the system by the Monte-Carlo method. Through predicting critical temperatures and calculating the superfluid density with the specific heat numerically, we have confirmed the following properties of the phase structure. If only the standard or modified kinetic terms are present, the system undergoes one, and only one phase transition, which is of BKT type, corresponding to vortex and half-vortex unbinding, respectively. If both terms are present, depending on the ratio between their coefficients, and by assuming that their square sum equals unity (a 2 + b 2 = 1), there exists either one or two transitions. If there is only one transition, it is never topological, and can be of both first and second order. If there are two transitions, then the one corresponding to the higher temperature is of BKT type, presumably related to half-vortex unbinding, while the other transition is of Ising type.
It would be interesting to improve upon the present renormalization group approximation scheme. Since higher field derivatives of the wavefunction renormalization factor could also play an important role for BKT-like transitions, one is interested in deriving a tower of equations for the aforementioned factors, and solve them simultaneously. Furthermore, the present scheme has only predicted the existence of a new endpoint of the line of fixed points, which indicated the existence of a second transition, but not any new fixed point that would correspond to the Ising transition at low temperature. It would be particularly important to find a scheme, which can overcome this shortcoming.
The results of this paper can be contrasted to another model admitting a vortex molecule solution of halfquantized vortices connected by a soliton, that is, coherently coupled Bose-Einstein condensates or two-gap superconductors [48] and spin-1 spinor Bose-Einstein condensates under the quadratic Zeeman field [49]. In this case, a two-step phase transition does not occur when two components are coupled by a Josephson interaction or a quadratic Zeeman field, while it can occur when they are decoupled. Essential differences between this case and that of the modified Goldstone model discussed in this paper are yet to be clarified.
Our study of the modified Goldstone model in two Euclidean dimensions has revealed that there exist twostep phase transitions related to half-quantized vortex molecules connected by domain walls. It is an open question whether there is any higher dimensional model allowing a two-step phase transition. For instance in three dimensions, a pair of a monopole and an anti-monopole connected by a string may play crucial role.