Parallel Tempering Monte Carlo Studies of Phase Transition of Free Boundary Planar Surfaces

Abstract

We numerically study surface models defined on hexagonal disks with a free boundary. 2D surface models for planar surfaces have recently attracted interest due to the engineering applications of functional materials such as graphene and its composite with polymers. These 2D composite meta-materials are strongly influenced by external stimuli such as thermal fluctuations if they are sufficiently thin. For this reason, it is very interesting to study the shape stability/instability of thin 2D materials against thermal fluctuations. In this paper, we study three types of surface models including Landau-Ginzburg (LG) and Helfirch-Polyakov models defined on triangulated hexagonal disks using the parallel tempering Monte Carlo simulation technique. We find that the planar surfaces undergo a first-order transition between the smooth and crumpled phases in the LG model and continuous transitions in the other two models. The first-order transition is relatively weak compared to the transition on spherical surfaces already reported. The continuous nature of the transition is consistent with the reported results, although the transitions are stronger than that of the reported ones.

1. Introduction

The two-dimensional surface model proposed by Helfrich is a model for biological membranes composed of lipid molecules, and it shares almost the same mathematical structure with the Polyakov’s rigid string model for elementary particles in subatomic scales [1,2]. In those models, the extrinsic curvature or bending energy plays an essential role in maintaining the smooth shape of surfaces, and because of its mathematical transparency, many studies have been conducted on the phase structure of the model between the smooth and crumpled phases [3,4,5,6,7,8,9]. The discrete surface models defined on triangulated lattices have also been extensively studied by Monte Carlo (MC) simulations [10,11,12,13,14,15,16,17,18,19].

However, the order of the crumpling transition is still controversial. Since the discrete model depends on the discretization of continuous Hamiltonian, we have a variety of discrete models [20,21]. Indeed, the curvature energy itself has a lot of variation such as extrinsic and intrinsic curvatures. For almost all of these discrete models, MC studies predict that the models undergo a first-order crumpling transition if the lattice is of spherical topology and allowed to self-intersect (⇔ self-intersecting) [20]. Here, we should note that only self-intersecting and fixed-connectivity lattice models are studied in this paper, and self-avoiding models and fluid surface models are not considered [15,16]. Intrinsic curvature models also have a first-order crumpling transition even on a disk surface [21], and the intrinsic curvature models are also out of consideration in this paper. The first-order crumpling transition is supported by theoretical studies on the basis of non-perturbative renormalization group techniques [8,9]. In contrast, it is reported that the order of transition is of second-order on the free boundary lattices [22,23]. Therefore, this continuous transition combining the above-mentioned first-order one indicates that the order of transition depends on the surface topology; either spheres or free boundary planar disks.

In this paper, we study three different surface models by the parallel tempering MC (PTMC) technique to check whether or not the transition is of second-order on triangulated disks with a free boundary. The PTMC technique was developed to simulate the spin glass models at very low temperature, where the standard Metropolis MC technique is not effective [24,25]. It is also reported that this PTMC technique can be applied to phenomena which undergo first-order transitions [26,27]. Therefore, we expect that the PTMC technique can also be used to study the phase structure of the surface models in this paper even if these models have first-order transitions [28].

We should here comment on the reason why the crumpling transition is of interest. Indeed, graphite oxide sheets in solvents have a crumpled state [29,30]. Crumpled states can also be observed in graphenes [31,32]. The surface condition of graphenes is altered by corrugations, and therefore ripples, wrinkles and crumples emerge [33,34]. These surface states modify or enhance the material properties such as mechanical, electrical and optical ones [35]. The crumpled graphene, for example, is expected to have enhanced chemical activities and energy storage capacities [36]. In addition to pure graphenes, polymer-graphene nano-composite or graphene-based polymers (or polymer-based graphenes) also has the crumpled states [37,38,39,40]. For the application of crumpled states of these graphene-based materials, it is interesting to study their stability against thermal fluctuation or some other stimuli in environmental conditions [41]. Therefore, the crumpled state is worthwhile studying in terms of phase transition.

This paper is organized as follows: In Section 2, the three different models and the PTMC technique are described in detail. Readers who are familiar with or not interested in these topics can skip Section 2 and go to Section 3, where the simulation results including snapshots of surfaces are presented. In Section 4, we summarize the results and comment on the future study.

2. Models and Monte Carlo Technique

2.1. Triangulated Disk

Discrete surface models are defined on such a triangulated disk shown in Figure 1. The lattice is characterized by the numbers $(N,N_E,N_T)$, which are the total number of vertices, the total number of edges, and the total number of triangles. Using the number L of division of the hexagon edge, we have the expressions for these numbers such that $(N,N_E,N_T)=(3L^2+3L+1,9L^2+3L,6L^2)$. The lattice shown in Figure 1 is obtained by $L=5$, and hence $(N,N_E,N_T)=(91,240,150)$. We chose a sufficiently small L to visualize the lattice structure. The lattices used in the simulations are larger than the lattice in Figure 1. The lattice spacing a for the edge length can be used as the length scale [42]. However, the simulation data are not directly compared to the experimental ones in this paper, and for this reason we fix a to $a=1$ for simplicity.

2.2. Landau-Ginzburg surface Model

The so-called Landau-Ginzburg surface model is first introduced and studied by Paczuski, Kardar and Nelson in [43], and this model is also numerically studied in [44]. Let $\mathbf{r}(\in {\mathbf{\Re }}^3)$ be the surface position. The continuous Hamiltonian is given by

$$\begin{array}{ccc}\hfill S_{\mathrm{LG}}\left(\mathbf{r}\right)& =& \frac{t}{2}\int d^{2}x{\left({\partial }_a\mathbf{r}\right)}^2+\frac{\kappa }{2}\int d^{2}x{\left({\partial }^2\mathbf{r}\right)}^2+u\int d^{2}x{\left({\partial }_a\mathbf{r}·{\partial }_b\mathbf{r}\right)}^2+v\int d^{2}x{\left({\partial }_a\mathbf{r}·{\partial }_a\mathbf{r}\right)}^2,\hfill \end{array}$$

(1)

where ${\partial }_a\mathbf{r}=\partial \mathbf{r}/\partial x_a,(a=1,2)$ is a tangential vector along the local coordinate axis $x_a$ on the surface and plays a role of the order parameter. The real numbers $t,\kappa ,u,v$ are the coefficients of the energy terms, which are square and quadratic with respect to ${\partial }_a\mathbf{r}$. The second term is defined by the square of second-order differentials ${\left({\partial }^2\mathbf{r}\right)}^2$.

The continuous energy $S_{\mathrm{LG}}\left(\mathbf{r}\right)$ in Equation (1) can be written more explicitly such that

$$\begin{array}{c}\hfill \begin{array}{c}{S}_{\mathrm{LG}}=t{S}_1+\kappa S_2+u{S}_3+v{S}_4,\hfill \\ S_1=\frac{1}{2}\int d^{2}x\left[{\left({\partial }_1\mathbf{r}\right)}^2+{\left({\partial }_2\mathbf{r}\right)}^2\right],\hfill \\ S_2=\frac{1}{2}\int d^{2}x{\left({\partial }_1^2\mathbf{r}+{\partial }_2^2\mathbf{r}\right)}^2=\frac{1}{2}\int d^{2}x\left({\partial }_1^2\mathbf{r}·{\partial }_1^2\mathbf{r}+{\partial }_2^2\mathbf{r}·{\partial }_2^2\mathbf{r}+2{\partial }_1^2\mathbf{r}·{\partial }_2^2\mathbf{r}\right),\hfill \\ S_3=\int d^{2}x\left[{\left({\partial }_1\mathbf{r}·{\partial }_1\mathbf{r}\right)}^2+{\left({\partial }_2\mathbf{r}·{\partial }_2\mathbf{r}\right)}^2+2{\left({\partial }_1\mathbf{r}·{\partial }_2\mathbf{r}\right)}^2\right],\hfill \\ S_4=\int d^{2}x\left[{\left({\partial }_1\mathbf{r}·{\partial }_1\mathbf{r}\right)}^2+{\left({\partial }_2\mathbf{r}·{\partial }_2\mathbf{r}\right)}^2+2{\left({\partial }_1\mathbf{r}\right)}^2{\left({\partial }_2\mathbf{r}\right)}^2\right].\hfill \end{array}\end{array}$$

(2)

The detailed information on the role of each term is written in Ref. [44], and we briefly describe the outline of each term. The first term $S_1$ is given by the integration of length squares of the tangential vectors, and it simply plays a role of the tensile energy. The second term $S_2$ plays a role of the bending energy and the in-plane shear energy, and the third term $S_3$ contains both the tensile and the shear energy components. The final term $S_4$ is a quadratic tensile energy term.

The discrete Hamiltonian is as follows [44]:

$$\begin{array}{c}\hfill \begin{array}{c}{S}_{\mathrm{LG}}=t{S}_1+\kappa S_2+u{S}_3+v{S}_4,\hfill \\ S_1=\frac{2}{3}\sum _{ij}{\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)}^2=\frac{2}{3}\sum _i{\mathbf{e}}_i^2,\hfill \\ S_2=\frac{1}{3}\sum _{ij}{\left({\mathbf{e}}_i-{\mathbf{e}}_j\right)}^2+\frac{1}{3}\sum _{\left(ij\right),\left(kl\right)}\left({\mathbf{e}}_i-{\mathbf{e}}_j\right)·\left({\mathbf{e}}_k-{\mathbf{e}}_l\right),\hfill \\ S_3=\frac{2}{3}\sum _{\Delta }\left[{\left({\mathbf{e}}_1^2\right)}^2+{\left({\mathbf{e}}_2^2\right)}^2+{\left({\mathbf{e}}_3^2\right)}^2+{\left({\mathbf{e}}_1·{\mathbf{e}}_2\right)}^2+{\left({\mathbf{e}}_2·{\mathbf{e}}_3\right)}^2+{\left({\mathbf{e}}_3·{\mathbf{e}}_1\right)}^2\right],\hfill \\ S_4=\frac{2}{3}\sum _{\Delta }\left[{\left({\mathbf{e}}_1^2\right)}^2+{\left({\mathbf{e}}_2^2\right)}^2+{\left({\mathbf{e}}_3^2\right)}^2+\left({\mathbf{e}}_1^2\right)\left({\mathbf{e}}_2^2\right)+\left({\mathbf{e}}_2^2\right)\left({\mathbf{e}}_3^2\right)+\left({\mathbf{e}}_3^2\right)\left({\mathbf{e}}_1^2\right)\right],\left(\mathrm{LG}\right).\hfill \end{array}\end{array}$$

(3)

We briefly summarize how to obtain the discrete Hamiltonian in Equation (3) from the continuous one in Equation (2). First of all, we note that the local coordinate origin of the triangle 123 in Figure 2a is at the vertex 1, and hence the differentials ${\partial }_1\mathbf{r}$ and ${\partial }_2\mathbf{r}$ in $S_1$ of Equation (2) are replaced by

$$\begin{array}{c}\hfill {\partial }_1\mathbf{r}\to {\mathbf{e}}_1={\mathbf{r}}_2-{\mathbf{r}}_1,{\partial }_2\mathbf{r}\to {\mathbf{e}}_2={\mathbf{r}}_3-{\mathbf{r}}_1.\end{array}$$

(4)

Recalling that there are two other local coordinate origins on the triangle 123, and including all possible terms for the differentials ${\partial }_1\mathbf{r}$ and ${\partial }_2\mathbf{r}$ and with the factor $1/3$, we have

$$\begin{array}{c}\hfill S_1=\frac{1}{3}\sum _{\Delta }\left({\mathbf{e}}_1^2+{\mathbf{e}}_2^2+{\mathbf{e}}_3^2\right),\end{array}$$

(5)

where ${\mathbf{e}}_3(={\mathbf{r}}_3-{\mathbf{r}}_2)$ is not written in Figure 2a. The integration $\int d^{2}x$ is replaced by the sum over triangles $\Delta$ such that $\int d^{2}x\to {\sum }_{\Delta }$. By replacing the summation convention from the sum over triangles ${\sum }_{\Delta }$ to the sum over bonds ${\sum }_i$, we have $S_1$ in Equation (3). We should note that this energy $S_1$ is also defined on the triangles of the boundary.

The second order derivative such as ${\partial }_1^2\mathbf{r}$ in $S_2$ is replaced by ${\partial }_1^2\mathbf{r}\to {\mathbf{e}}_j-{\mathbf{e}}_i$ using the vectors ${\mathbf{e}}_j$ and ${\mathbf{e}}_i$ on the hexagon in Figure 2b, because ${\mathbf{e}}_j$ and ${\mathbf{e}}_i$ are considered as the discretization of ${\partial }_1\mathbf{r}$ along the coordinate axis $2^{\prime }12$ corresponding to the local coordinate axis $x_1$ at the vertex 1. Another diagonal line $3^{\prime }13$ on this hexagon is considered as the $x_2$ axis, and thereby the square of Laplacian ${\left({\partial }^2\mathbf{r}\right)}^2={({\partial }_1^2\mathbf{r}+{\partial }_2^2\mathbf{r})}^2$ is replaced by ${({\mathbf{e}}_j-{\mathbf{e}}_i)}^2+{({\mathbf{e}}_l-{\mathbf{e}}_k)}^2+2({\mathbf{e}}_j-{\mathbf{e}}_i)·({\mathbf{e}}_l-{\mathbf{e}}_k)$ using the vectors ${\mathbf{e}}_i$, ${\mathbf{e}}_j$,${\mathbf{e}}_k$ and ${\mathbf{e}}_l$ in Figure 2b. Thus, we have $S_2$ in Equation (3) as a discrete bending energy corresponding to the continuous one $(1/2)\int d{x}^2{\left({\partial }^2\mathbf{r}\right)}^2$ in Equation (3). The reason for the factor $1/3$ in the discrete expression is that every vertex inside the boundary is assumed to be the center of hexagon, and therefore the summation is triply duplicated. Strictly speaking, the duplication at the vertices close to the boundary is not triple, however, the coefficient is not so important and is simply fixed to $1/3$, which is the same as in the spherical model [44]. In the discrete $S_2$, ${\sum }_{ij}$ denotes the sum over the three different diagonal lines, and ${\sum }_{\left(ij\right),\left(kl\right)}$ denotes the sum over the corresponding local coordinates on the hexagon.

Please note that on the hexagonal lattice, such as shown in Figure 1, the coordination number at the vertices inside the boundary is given by $q=6$, where the coordination number $q_i$ is the total number of bonds emanating from the vertex i. This is in sharp contrast to spherical lattices, which must include the vertices with $q\ne 6$. Therefore, the vertex with $q=6$ in Figure 2b is sufficient for the discretization of ${\left({\partial }^2\mathbf{r}\right)}^2$ for all internal vertices. On the boundary vertices, the definition of energy $S_2$ is slightly different from that on the internal vertices. On the vertices with $q=4$, the square of Laplacian ${\left({\partial }^2\mathbf{r}\right)}^2={({\partial }_1^2\mathbf{r}+{\partial }_2^2\mathbf{r})}^2$ is simply replaced by ${({\mathbf{e}}_j-{\mathbf{e}}_i)}^2$ instead of ${({\mathbf{e}}_j-{\mathbf{e}}_i)}^2+{({\mathbf{e}}_l-{\mathbf{e}}_k)}^2+2({\mathbf{e}}_j-{\mathbf{e}}_i)·({\mathbf{e}}_l-{\mathbf{e}}_k)$. Moreover, the bending energy $S_2$ is not defined on the vertices with $q=3$ on the boundary (see Figure 1).

In the continuous $S_3$ and $S_4$ of Equation (2), the derivatives ${\partial }_1\mathbf{r}$ and ${\partial }_2\mathbf{r}$ are replaced by ${\mathbf{e}}_1$ and ${\mathbf{e}}_2$ in Equation (4) on the triangle in Figure 2a. The discretization technique is exactly the same as the one assumed in $S_1$, and hence we have the discrete energies $S_3$ and $S_4$ in Equation (3).

The partition function Z is given by the multiple integrations of the vertex positions such that

$$\begin{array}{c}\hfill Z={\int }^{\prime }\prod _{i}d{\mathbf{r}}_{i}exp\left(-S_{\mathrm{LG}}\right),\end{array}$$

(6)

where the prime in ${\int }^{\prime }{\prod }_{i}d{\mathbf{r}}_i$ denotes that the center of the mass of surface is fixed at the origin of ${\mathbf{\Re }}^3$.

From the scale invariance of Z, we have

$$\begin{array}{c}\hfill \langle S_1^{\prime }\rangle /N=3/2,S_1^{\prime }=t{S}_1+\kappa S_2+2u{S}_3+2v{S}_4,\end{array}$$

(7)

where $\langle Q\rangle$ is defined by $\langle Q\rangle ={\int }^{\prime }{\prod }_{i}d{\mathbf{r}}_{i}Qexp\left(-S_{\mathrm{LG}}\right)/Z$ [7,44]. This relation in Equation (7) can be used to check whether the simulation is correctly performed or not.

2.3. Canonical Model

We start with the continuous form of Hamiltonian $S(\mathbf{r},g)$ of the canonical model, which is defined by a mapping $\mathbf{r}$ from a two-dimensional surface M to the three-dimensional Euclidean space ${\mathbf{\Re }}^3$, such that

$$\begin{array}{c}\hfill \mathbf{r}:M\ni (x_1,x_2)⟼\mathbf{r}(x_1,x_2)=(X,Y,Z)\in {\mathbf{\Re }}^3.\end{array}$$

(8)

The variable $\mathbf{r}$, which is originally used to denote the surface position in ${\mathbf{\Re }}^3$, is now used for the symbol of the mapping. Another variable denoted by g in $S(\mathbf{r},g)$ is the metric function $g_{ab}$ on M, where $g_{ab}$ is a $2\times 2$ matrix. The metric $g_{ab}$ originally is not a variable but is determined by a local coordinate, which is fixed arbitrarily by hand from the reparametrization invariance. This invariance is a symmetry of the model under 2D coordinate transformations on the surface $\mathbf{r}\left(M\right)$ in ${\mathbf{\Re }}^3$. Thus, the model is slightly extended from the original one in the sense that $g_{ab}$ is a variable that should be physically determined. As a consequence, we have the possibility to obtain surfaces which cannot be in ${\mathbf{\Re }}^3$ in the extended model, whether this is meaningful or not. For example, let us consider the metric

$$g_{ab}=\left(\begin{array}{cc}{\ell }_{12}^2& {\ell }_{12}{\ell }_{13}cos\mathsf{\Phi }\\ -{\ell }_{12}{\ell }_{13}cos\mathsf{\Phi }& {\ell }_{13}^2\end{array}\right)$$

(9)

on a surface discretized by piece-wise linear triangles such as in Figure 3a, where $cos\mathsf{\Phi }$ is not always identical to $cos\varphi ={\overrightarrow{\ell }}_{12}·{\overrightarrow{\ell }}_{13}/{\ell }_{12}{\ell }_{13}$ between the edge vectors ${\overrightarrow{\ell }}_{12}$ and ${\overrightarrow{\ell }}_{13}$ [45]. If $\mathsf{\Phi }=\varphi$ for all triangles, then this metric is identical with the discrete induced metric $g_{ab}={\mathbf{e}}_a·{\mathbf{e}}_b$. However, if the angles $\left\{\mathsf{\Phi }\right\}$ do not satisfy the triangle equality, i.e., the sum of three internal angles is not always equal to $\pi$, then the surface with such metric in Equation (9) generally cannot be realized in ${\mathbf{\Re }}^3$.

The continuous Hamiltonian is given by

$$\begin{array}{c}\hfill S(\mathbf{r},g)=\gamma {\int }_M\sqrt{g}{d}^{2}x{g}^{ab}{\partial }_a\mathbf{r}·{\partial }_b\mathbf{r}+\frac{\kappa }{2}{\int }_M\sqrt{g}{d}^{2}x{g}^{ab}{\partial }_a\mathbf{n}·{\partial }_b\mathbf{n},(\gamma =1)\end{array}$$

(10)

where the surface tension coefficient $\gamma$ is always fixed to $\gamma =1$. The symbol $g^{ab}$ denotes the inverse of $g_{ab}$, and g is its determinant. We should note that $g_{ab}$ assumed in the expression of $S(\mathbf{r},g)$ in Equation (10) is a variable that should be determined just like the mapping $\mathbf{r}$ as mentioned above [46]. Note also that this Hamiltonian is a two-dimensional extension of the polymer model of Doi-Edwards [47].

Here, we should note that the real surface $\mathbf{r}\left(M\right)$ in ${\mathbf{\Re }}^3$ corresponding to the material under consideration is described by the induced metric $g_{ab}={\partial }_a\mathbf{r}·{\partial }_b\mathbf{r}$. This allows us to consider that the surface with a given metric $g_{ab}$ is different from the real surface, pointing to the possibility for the surface to correspond to this $g_{ab}$. Therefore, from this set of surfaces, a physically meaningful surface should be uniquely determined by the modeling. For this reason, we introduce a two-dimensional surface M in addition to the real surface with ${\partial }_a\mathbf{r}·{\partial }_b\mathbf{r}$, both of which should be physically determined. This is another extension of the surface model, and this is the meaning of the mapping described in Equation (8). Please note that the surface M is not necessarily a manifold [45,48,49].

The problem is then how to determine $g_{ab}$ for M. One possible and simple technique is to fix $g_{ab}$ to the Euclidean metric such that $g_{ab}={\delta }_{ab}$. In this case, M is a simple two-dimensional Euclidean space and plays no role in describing the model, however, we use M to express the domain in the two-dimensional integrations in $S\left(\mathbf{r}\right)$. Thus, the only variable to be determined is $\mathbf{r}$, and $S\left(\mathbf{r}\right)$ is now given by

$$\begin{array}{c}\hfill S\left(\mathbf{r}\right)=\gamma {\int }_M{d}^{2}x{\partial }_a\mathbf{r}·{\partial }_a\mathbf{r}+\frac{\kappa }{2}{\int }_M{d}^{2}x{\partial }_a\mathbf{n}·{\partial }_a\mathbf{n},(\gamma =1).\end{array}$$

(11)

One simple reason for why we assume the Euclidean metric ${\delta }_{ab}$ for $g_{ab}$ is as follows: The Hamiltonian in Equation (10) is invariant under an arbitrary conformal transformation $g_{ab}\to g_{ab}^{\prime }=f\left(x\right)g_{ab}$ with a positive function $f\left(x\right)$ on M. This invariance is described by $S(\mathbf{r},g)=S(\mathbf{r},g^{\prime })$ and implies that the metric $g_{ab}$ can be chosen relatively freely, and therefore it is fixed to the simplest one such that $g_{ab}={\delta }_{ab}$.

The discrete Hamiltonians are obtained from the continuous one in Equation (11) by the replacement of the differentials in Equation (4). Here, we use the symbols ${\overrightarrow{\ell }}_{12}={\mathbf{e}}_1$ and ${\overrightarrow{\ell }}_{13}={\mathbf{e}}_2$ for the edge vectors (Figure 3a), and ${\ell }_{ij}=|{\overrightarrow{\ell }}_{ij}|=|{\mathbf{r}}_j-{\mathbf{r}}_i|$ for the edge (or bond) length. Thus, we obtain

$$\begin{array}{c}S\left(\mathbf{r}\right)=S_1+\kappa S_2,\hfill \\ S_1=\sum _{(i,j)}{\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)}^2=\sum _{(i,j)}{\ell }_{ij}^2,S_2=\sum _{(i,j)}\left(1-{\mathbf{n}}_i·{\mathbf{n}}_j\right),\left(\mathrm{cano}\right),\hfill \end{array}$$

(12)

where the factor $4/3$ is eliminated from both $S_1$ and $S_2$ for simplicity. This $S_1$ is exactly the same as $S_1$ of the LG model in Equation (3) up to the numerical factor. The symbol ${\sum }_{\left(ij\right)}$ in $S_1$ denotes the sum over all bonds $\left(ij\right)$ connecting the vertices i and j. In contrast, $\left(ij\right)$ in the sum ${\sum }_{\left(ij\right)}$ of $S_2$ denotes the triangles sharing a common bond (Figure 3b), and the unit normal vector ${\mathbf{n}}_i$ is defined on the triangle i. The partition function of the canonical model is exactly the same as Z in Equation (6) for the LG model except the Hamiltonian in the Boltzmann factor.

2.4. Modified Canonical Model

The third model, which we call “modified canonical model”, is obtained from the same continuous Hamiltonian in Equation (11) assumed for the canonical model. The only difference is a discretization of the bending energy $S_2$, where the unit normal vector $\mathbf{n}\left(i\right)$ at the vertex i (Figure 4a) is used as well as the normal vector ${\mathbf{n}}_i$ on the triangle i (Figure 3b). The discrete Hamiltonian is given by [20]

$$\begin{array}{c}S\left(\mathbf{r}\right)=S_1+\kappa S_2,\hfill \\ S_1=\sum _{(i,j)}{\ell }_{ij}^2,S_2=\sum _{i=1}^N\sum _{j\left(i\right)}\left[1-\mathbf{n}\left(i\right)·{\mathbf{n}}_{j\left(i\right)}\right],\left(\mathrm{modi}\right),\hfill \end{array}$$

(13)

where ${\mathbf{n}}_{j\left(i\right)}$ in $S_2$ is the unit normal vector of the triangle $j\left(i\right)$ connected to the vertex i. The normal vector $\mathbf{n}\left(i\right)$ at the vertex i is defined by (Figure 4a)

$$\mathbf{n}\left(i\right)={\mathbf{N}}_i/\left|{\mathbf{N}}_i\right|,{\mathbf{N}}_i=\sum _{j\left(i\right)}{\mathbf{n}}_{j\left(i\right)}{A}_{j\left(i\right)},$$

(14)

where $A_{j\left(i\right)}$ is the area of the triangle $j\left(i\right)$. Please note that the interaction range of the normal vectors ${\mathbf{n}}_i$ is slightly larger than that of the canonical model (Figure 4b). In fact, only two nearest neighbor vectors ${\mathbf{n}}_i$ and ${\mathbf{n}}_j$ are directly coupled in $S_2$ in Equation (12) of the canonical model, and as a consequence only three nearest neighbor vectors ${\mathbf{n}}_i(i=1,2,3)$ are coupled to ${\mathbf{n}}_0$ (see Figure 3b). In contrast, as shown in Figure 4b, the non-nearest neighbor ${\mathbf{n}}_i$ and ${\mathbf{n}}_j$, of which the triangles i and j do not directly contact each other, are coupled to ${\mathbf{n}}_0$ in $S_2$ in Equation (13) of the modified model.

2.5. Parallel Tempering Monte Carlo Technique

The so-called parallel tempering Monte Carlo (PTMC) technique developed for the spin glass model at low temperatures is successfully applied to the first-order crumpling transition of the canonical model on spherical lattices [28]. In this subsection, the outline of the PTMC technique applied to the tethered surface model is briefly presented.

Let $(\mathbf{r},\kappa )$ represent a system of configuration $\mathbf{r}(=\{{\mathbf{r}}_1,{\mathbf{r}}_2,\cdots ,{\mathbf{r}}_N\})$ with a given $\kappa$, which is assumed to be changed. In the case of the LG model, the parameter $\kappa$ is also changed while the other three parameters $(t,u,v)$ are fixed in the simulations in this paper. In this PTMC, $N_R$ replicas $\{({\mathbf{r}}_1,{\kappa }_1),({\mathbf{r}}_2,{\kappa }_2),\cdots ,({\mathbf{r}}_{N_R},{\kappa }_{N_R})\}$ are simulated in parallel by the standard Metropolis MC (MMC) simulation technique [50,51], and the systems are exchanged after sufficiently long MMC runs. Because of this exchange, the total number of different combinations $\{({\mathbf{r}}_1,{\kappa }_1),({\mathbf{r}}_2,{\kappa }_2),\cdots ,({\mathbf{r}}_{N_R},{\kappa }_{N_R})\}$ is $N_R!$.

The $N_R$ replicas are updated as follows:

(P₁)

Perform long MMC simulations for $N_R$ replicas

(P₂)

Exchange all nearest neighbor systems $(\mathbf{r},\kappa )$ and $({\mathbf{r}}^{\prime },{\kappa }^{\prime })$ with the probability

$$\begin{array}{c}\hfill W(\mathbf{r},\kappa |{\mathbf{r}}^{\prime },{\kappa }^{\prime })=\mathrm{Min}\left[1,exp(-\Delta )\right],\Delta =({\kappa }^{\prime }-\kappa )\left[S_2\left(\mathbf{r}\right)-S_2\left({\mathbf{r}}^{\prime }\right)\right]\end{array}$$

(15)

(P₃)

Repeat ${\mathrm{P}}_1$ and ${\mathrm{P}}_2$

We should note that the process ${\mathrm{P}}_1$ can be performed in parallel as mentioned above. In the exchange process ${\mathrm{P}}_2$, only bending energy $S_2$ is used, and no information on the other energy $S_1$ is used in the canonical and modified canonical models. This is also true for the LG model, and no information on the energies $S_1$, $S_3$ and $S_4$ is left out in the exchange process.

Figure 5a,b intuitively show the difference of MMC and PTMC simulations for four replica systems. The numbers $1,2,3,4$ denote the bending rigidities such as ${\kappa }_1,\cdots ,{\kappa }_4$, and the color blocks denote the configurations ${\mathbf{r}}_1,\cdots ,{\mathbf{r}}_4$. The combination of $({\mathbf{r}}_i,{\kappa }_j)$ is exchanged as the iterations evolve in the PTMC simulation.

Here, we briefly show that all micro states $\{\mathbf{r},\kappa \}$ satisfy the canonical Boltzmann distribution as a result of PTMC simulations. Let $P\left(\right\{\mathbf{r},\kappa \left\}\right)$ be a probability distribution for all the states $\{\mathbf{r},\kappa \}$ such that

$$\begin{array}{c}\hfill P\left(\{\mathbf{r},\kappa \}\right)=\prod _{m=1}^{N_R}{P}_{\mathrm{eq}}({\mathbf{r}}_m,{\kappa }_m),P_{\mathrm{eq}}(\mathbf{r},\kappa )=Z^{-1}exp\left[-S(\mathbf{r},\kappa )\right],S(\mathbf{r},\kappa )=\overline{S}\left(\mathbf{r}\right)+\kappa S_2\left(\mathbf{r}\right),\end{array}$$

(16)

where $\overline{S}\left(\mathbf{r}\right)$ is independent of $\kappa$ and given by $\overline{S}\left(\mathbf{r}\right)=t{S}_1\left(\mathbf{r}\right)+u{S}_3\left(\mathbf{r}\right)+v{S}_4\left(\mathbf{r}\right)$ for the LG model and $\overline{S}\left(\mathbf{r}\right)=S_1\left(\mathbf{r}\right)$ for the other two models. The $P_{\mathrm{eq}}(\mathbf{r},\kappa )$ in Equation (16) is the Boltzmann distribution function of the state $(\mathbf{r},\kappa )$. If the PTMC exchange satisfies the detailed balance condition described by

$$\begin{array}{c}\hfill P(\dots ,;\mathbf{r},\kappa ;\dots ;{\mathbf{r}}^{\prime },{\kappa }^{\prime };\dots )W(\mathbf{r},\kappa |{\mathbf{r}}^{\prime },{\kappa }^{\prime })=P(\dots ,;{\mathbf{r}}^{\prime },\kappa ;\dots ;\mathbf{r},{\kappa }^{\prime };\dots )W({\mathbf{r}}^{\prime },\kappa |\mathbf{r},{\kappa }^{\prime }),\end{array}$$

(17)

then we understand from the well-known uniqueness theorem that this P is the uniquely determined probability [28]. Please note that the right-hand side of Equation (17) is obtained from the left-hand side by exchanging $\mathbf{r}$ and ${\mathbf{r}}^{\prime }$. Thus, we should prove that $P\left(\right\{\mathbf{r},\kappa \left\}\right)$ in Equation (16) satisfies the condition in Equation (17). This can be accomplished in two steps. The first step is to see that Equation (17) for $P\left(\right\{\mathbf{r},\kappa \left\}\right)$ in Equation (16) is equivalent with the relation

$$\begin{array}{c}\hfill \frac{W(\mathbf{r},\kappa |{\mathbf{r}}^{\prime },{\kappa }^{\prime })}{W({\mathbf{r}}^{\prime },\kappa |\mathbf{r},{\kappa }^{\prime })}=exp\left(-\Delta \right)\end{array}$$

(18)

under the condition given by Equation (15). The second step is to see that the relation in Equation (18) is correct. This second step is almost trivial from the definitions of $W(\mathbf{r},\kappa |{\mathbf{r}}^{\prime },{\kappa }^{\prime })$ and $\Delta$ in Equation (15). The first step is also straightforward to prove.

The assumed parameters for the simulations including the total number of MC sweeps (MCS) are listed in

Table 1. The parameters assumed for the simulations; $n_{P_1}$ denotes the total number of MMC iterations performed in the $P_1$ process per 1 MCS for each replica, ${\kappa }_1$ and ${\kappa }_{N_R}$ (${\kappa }_1<{\kappa }_{N_R}$) are the bending rigidity of the replica 1 and $N_R$, and $\Delta \kappa (=({\kappa }_{N_R}-{\kappa }_1)/N_R)$ is the difference of $\kappa$ between two neighboring replicas.

Model	N	#Total (MCS)	#Therm (MCS)	${\mathit{n}}_{{\mathit{P}}_1}$	${\mathit{N}}_{\mathit{R}}$	${\mathit{\kappa }}_1$	${\mathit{\kappa }}_{{\mathit{N}}_{\mathit{R}}}$	$\mathit{\Delta }\mathit{\kappa }$
LG	7351	$2.5\times {10}^8$	$2.5\times {10}^7$	10	24	0.1835	0.187	$1.46\times {10}^{-4}$
LG	5677	$2.5\times {10}^8$	$2.5\times {10}^7$	10	24	0.1842	0.1878	$1.5\times {10}^{-4}$
LG	4219	$9\times {10}^7$	$3\times {10}^7$	10	24	0.183	0.194	$4.58\times {10}^{-4}$
LG	2611	$1\times {10}^7$	$2\times {10}^6$	10	24	0.18	0.2	$8.33\times {10}^{-4}$
cano	44,287	$3.8\times {10}^7$	$1.8\times {10}^7$	20	24	0.766	0.782	$2.5\times {10}^{-4}$
cano	30,907	$8\times {10}^7$	$3\times {10}^7$	20	24	0.764	0.787	$9.58\times {10}^{-4}$
cano	20,917	$1\times {10}^8$	$1.5\times {10}^7$	20	24	0.762	0.792	$1.25\times {10}^{-3}$
cano	12,097	$7.5\times {10}^7$	$7.5\times {10}^6$	10	24	0.77	0.804	$1.42\times {10}^{-3}$
modi	20,917	$6.4\times {10}^7$	$1\times {10}^7$	20	24	0.451	0.467	$6.67\times {10}^{-4}$
modi	12,097	$1.9\times {10}^8$	$2.5\times {10}^7$	10	24	0.448	0.47	$9.167\times {10}^{-4}$
modi	7351	$5\times {10}^8$	$5\times {10}^7$	10	16	0.444	0.476	$2\times {10}^{-2}$
modi	4219	$5\times {10}^8$	$1.5\times {10}^7$	10	16	0.44	0.49	$6.96\times {10}^{-2}$

, where 1 MCS consists of the processes $P_1$ and $P_2$ of PTMC. In Table 1, the symbol N is the total number of vertices, #total (MCS) and #therm (MCS) are the total number of MCS and the total number of thermalization MCS, and $n_{P_1}$ denotes the number of iterations performed in $P_1$ per 1 MC sweep. This $n_{P_1}$ is fixed to $n_{P_1}=10$ (or $n_{P_1}=20$), and this implies that the total number of MMC (including thermalization) iterations for each replica is 10 (or 20) times larger than the #total (MCS). The $N_R$ is the total number of replicas, and ${\kappa }_1$ and ${\kappa }_{N_R}$ (${\kappa }_1<{\kappa }_{N_R}$) are the bending rigidity of the replica 1 and $N_R$, $\Delta \kappa (=({\kappa }_{N_R}-{\kappa }_1)/N_R)$ is the difference of $\kappa$ between two neighboring replicas, which will be exchanged in the process $P_2$. The total number of iterations #total for the large lattices in the latter two models is not so large compared to those for smaller lattice in the LG model.

3. Simulation Results

3.1. Snapshots

Firstly, in the presentation section, we show snapshots, which are obtained as one of the configurations from the replicas $i=1$ and $i=N_R$ in each model. In the upper (lower) low in Figure 6, the snapshots are obtained from the replica $i=N_R$ ($i=1$) of the (a) LG, (b) canonical, and (c) modified canonical models. The size of lattices are the largest for these snapshots, and hence from Table 1 all $N_R$ are given by $N_R=24$. From Table 1, we find that the values of $\kappa$ corresponding to these replicas are given as follows: (a) $\kappa =0.1835$ (upper), $\kappa =0.187$ (lower), (b) $\kappa =0.766$ (upper), $\kappa =0.782$ (lower), and (c) $\kappa =0.451$ (upper), $\kappa =0.467$ (lower). In the LG model simulations, the parameters $(t,u,v)$ are fixed to

$$\begin{array}{c}\hfill (t,u,v)=(-6,0.2,0.2),\left(\mathrm{LG}\right).\end{array}$$

(19)

We find that the surfaces in the upper low are globally bending but almost flat while those in the lower low shrink to a small ball. From this observation, we understand that the surfaces in the upper and lower lows are in the flat and crumpled phases, respectively, in all of the models.

3.2. Bending Energy and Mean Square Gyration

The bending energy of the LG model is calculated only on the internal vertices, because the definition of $S_2$ on the boundary vertices is slightly different from that on the internal vertices as described in Section 2.2. In Figure 7a–i, we plot the bending energy $S_2/N_B$ of the LG model and the other models, the specific heat

$$\begin{array}{c}\hfill C_{S_2}=\frac{{\kappa }^2}{N}〈{\left(S_2-\langle S_2\rangle \right)}^2〉,\end{array}$$

(20)

and the peak $C_{S_2}^{\max }$ of the specific heat. We should note that $N_B$ is the total number of internal bonds $N_B=3(N-6L)$ for the LG model, where $N-6L$ is the total number of internal vertices and hence $3(N-6L)$ is the total number of the diagonal lines such as $2^{\prime }12$ in Figure 2b. The reason why this $N_B$ is used for $S_2$ in the LG model is because the first term of $S_2$ in Equation (3) is defined on such diagonal lines, and the second term is also defined on a pair of diagonal lines $2^{\prime }12$ and $3^{\prime }13$ in Figure 2b, and there are three different such diagonal lines on each internal vertex. For the other two models, $N_B(=N_E-6L)$ is the total number of internal bonds, on which the bending energy $S_2$ is defined, where $N_E$ is the total number of edges $N_E$ given in Section 2.1. Please note that $C_{S_2}$ in Equation (24) for the LG model is defined also by using N. The curves of $C_{S_2}$ in Figure 7h for the modified model are fluctuating, and this fluctuation may be due to the fact that the total number of iterations is not sufficient for this model as mentioned above, though reasonable peak values $C_{S_2}^{\max }$ are obtained.

From the LG model data plotted in Figure 7a–c, we see that the $S_2/N_B$ abruptly changes against $\kappa$, and correspondingly the $C_{S_2}$ has the anomalous peak showing the existence of the crumpling transition. As mentioned above, the parameters $(t,u,v)$ are fixed to the values in Equation (19) and remain unchanged. The peak heights $C_{S_2}^{\max }$ vs. N are plotted in log-log scale in Figure 7c, where N is used. The data obtained by the other models are also plotted in Figure 7d–i. Fitting the largest three data of $C_{S_2}^{\max }$ to the scaling relation

$$\begin{array}{c}\hfill C_{S_2}^{\max }\sim N^{\sigma },\end{array}$$

(21)

we have

$$\begin{array}{c}\hfill \sigma =1.96\pm 0.07\left(\mathrm{LG}\right),\sigma =0.72\pm 0.01\left(\mathrm{cano}\right),\sigma =0.79\pm 0.05\left(\mathrm{modi}\right).\end{array}$$

(22)

We find from these results that the transition of the LG model is of the first order while in the other two models the transition is of second-order, because $\sigma \ge 1$ in the LG model and $\sigma <1$ in the other models. We note that $\sigma$ should be $\sigma =1$ for the first order transitions, because $C_{S_2}^{\max }$ is expected to scale as $L^2(=N)$, where L is the linear size of the lattice [52]. However, the present result is almost twice of this expectation and indicates that the dimension D of lattice is effectively $D=4$. One possible reason is that the surface is self-intersecting, and hence, the total number or density of possible configurations is larger than that in the case of $D=2$. The first result in Equation (22) of the LG model is qualitatively comparable to $\sigma =1.58\pm 0.08$ of the canonical model on spherical lattice in Ref. [28], because both results indicate that the transition is of first-order. In contrast, the latter two results in Equation (22) of the other models indicate a second-order transition, and hence these two are completely different from the result in Ref. [28]. We consider that this difference simply stems from the difference in the lattice topology or structure; sphere and disk, which are surfaces without a boundary and with a free boundary, or are compact and non-compact.

We calculate the coefficient $\sigma$ in Equation (21) for the LG model using the specific heat $C_{S_2}$ of the bending energy $S_2$ in Equation (12), and we have $\sigma =1.96\pm 0.07$, which is almost identical to the first result $\sigma =1.97\pm 0.07$ in Equation (22). This indicates that the results $\sigma =1.96\pm 0.07$ in Equation (22) for the LG model is reliable, although this result is obtained from the bending energy on the internal vertices only.

The problem here is ascertaining why the result of LG model in Equation (22) is different from the others. Firstly, we should be reminded that it is predicted that the LG model has a continuous crumpling or wrinkling transition from the mean filed analysis [11]. However, as mentioned in Ref. [28], the basic assumption for the mean field analysis that the surface fluctuation is relatively small is not always true even for the models. This is considered to be why the result of the LG model deviates from the mean field prediction and has the first-order transition even on the lattice with a free boundary. Next, we have to consider where the difference comes from in the order of transitions between the LG model and the other two models. One possible origin is the shear stress or resistance to shear deformation expected in the LG model; it is not expected in the other models. The bending energy $S_2$ assumed in the LG models is identical to $S_2$ in the other models up to a numerical factor, and this $S_2$ has no shear resistance to in-plane deformation of triangles; it has only a resistance to the bending or out-of-plane deformation. Thus, the only source for the shear resistance is the Gaussian bond potential $S_1$ in the two models, while in the LG model, $S_3$ and $S_4$ as well as $S_1$ have resistance to shear stresses. Here, we should note that all models are fixed-connectivity or tethered models, and the shear resistance is not negligible compared to the fluid surface models on dynamically triangulated lattices. In the tethered models, any deformations of triangles except simple expansion/shrinkage accompany a shear resistance, because all triangles tend to become regular in the equilibrium configurations. Please note that the simple expansion and shrinkage of triangles are suppressed because the mean bond length remains constant due to the scale invariant-property of the partition function.

The mean square radius of gyration $R_g^2$ defined by

$$\begin{array}{c}\hfill R_g^2=\frac{1}{N}\sum _i{\left({\mathbf{r}}_i-\overline{\mathbf{r}}\right)}^2,\overline{\mathbf{r}}=\frac{1}{N}\sum _i{\mathbf{r}}_i,\end{array}$$

(23)

its variance

$$\begin{array}{c}\hfill C_{R_g}=\frac{1}{N}〈{\left(R_g^2-\langle R_g^2\rangle \right)}^2〉,\end{array}$$

(24)

and the peak values $C_{R_g}^{\max }$ are plotted in Figure 8a–i. The fluctuations of the data $C_{R_g}$ in Figure 8h are relatively large due to the same reason for $C_{S_2}$ mentioned above. We also find that $R_g^2$ rapidly changes, where $C_{R_g}$ has a peak $C_{R_g}^{\max }$. The peak position on the $\kappa$ axis for $C_{R_g}^{\max }$ of the LG model is almost identical with the position for $C_{S_2}^{\max }$, while those for $C_{R_g}^{\max }$ and $C_{S_2}^{\max }$ in the other models are considerably different from each other. This difference is not observed in the canonical model on spherical lattice in Ref. [28], where a first-order transition is expected and where the transition point is very clear and detected uniquely even by numerical simulations. In contrast, continuous transitions are relatively unclear in general, and hence the transition points of the latter two models are relatively hard to observe. Moreover, the total number of iterations for these are not always sufficiently large as mentioned above. These are possible reasons for the deviation of the transition points observed in $C_{R_g}^{\max }$ and $C_{S_2}^{\max }$.

The peak values $C_{R_g}^{\max }$ are expected to scale according to

$$\begin{array}{c}\hfill C_{R_g}^{\max }\sim N^{\mu }.\end{array}$$

(25)

By fitting the largest three data of $C_{R_g}^{\max }$ in Figure 8c,f,i to this relation, we have

$$\begin{array}{c}\hfill \mu =1.73\pm 0.08\left(\mathrm{LG}\right),\mu =0.93\pm 0.12\left(\mathrm{cano}\right),\mu =0.89\pm 0.08\left(\mathrm{modi}\right).\end{array}$$

(26)

These results also support the proposition that the LG model has a first-order transition and the other models a second-order transition, though the coefficient of the canonical model is close to $\mu =1$ and the transition is close to a first-order one. The reason why $\mu$ in Equation (26) of the LG model is $\mu >1$ may be the fact that the surface is not self-avoiding [28]. We should note that the order of transition depends on its definition. If we call a transition first-order only if the coefficient $\mu$ for the variance $C_{S_{*}}$ of Hamiltonian $S_{*}$ is $\mu =1$, then the canonical and modified models clearly have a second-order transition from the results $\sigma$ in Equation (22). In general, the coefficient $\sigma$ for the bending energy $S_2$ is more important than $\mu$ as a coefficient for the determination of the order of transition. Thus, it is reasonable to consider that the crumpling transitions of the canonical and modified canonical models are continuous transitions.

3.3. Binder Quantity and Fractal Dimension

Firstly. in this subsection, we calculate the Binder quantity $B_{S_2}$ defined by [53,54]

$$B_{S_2}=1-\frac{〈{\left(S_2-\langle S_2\rangle \right)}^4〉}{3{〈{\left(S_2-\langle S_2\rangle \right)}^2〉}^2}$$

(27)

for the LG model (Figure 9a). It is expected that $B_{S_2}$ has a peak $B_{S_2}^{\max }$ and $B_{S_2}^{\max }\to 2/3$ at the first-order transition point. To see whether this expectation is satisfied or not, $B_{S_2}^{\max }$ vs. $1/N$ are plotted in Figure 9b, where the Binder quantity $B_{R_g}^{\max }$ for $R_g^2$ is also plotted. The solid lines are also drawn by Mathematica command “Interpolation” for the data of $B_{S_2}^{\max }$ and $B_{R_g}^{\max }$. The value $B_{S_2}^{\max }(N\to \infty )$ on the vertical axis expected from the extrapolations is almost identical with $2/3$, while $B_{R_g}^{\max }(N\to \infty )$ is slightly smaller than $2/3$. This deviation of $B_{R_g}^{\max }(N\to \infty )$ seems due to the size effect, because the lattice size is not so large compared to those used in Ref. [28], where $B_{S_2}^{\max }(N\to \infty )\simeq 0.7$ and $B_{R_g}^{\max }(N\to \infty )\simeq 0.69$ are observed.

A first-order transition is also reflected in the Binder cumulant $V_{S_2}$, which is defined by [53,54]

$$V_{S_2}=1-〈S_2^4〉/\left(3{\langle S_2^2\rangle }^2\right).$$

(28)

The Binder cumulant $V_{R_g}$ for $R_g^2$ is also defined analogously to Equation (28). These quantities are expected to have the minimum $V_{S_2}^{\min }$ and $V_{R_g}^{\min }$ at the first-order transition point, and this expectation is confirmed in Figure 9c, where only $V_{S_2}^{\min }$ is plotted. It is also found that the position of $V_{S_2}^{\min }$ on the $\kappa$ axis is almost the same as that of $C_{S_2}^{\max }$ in Figure 7b. This also implies that $V_{S_2}^{\min }$ reflects the first-order transition.

In the other two models, the quantities $B_{S_2}$ and $V_{S_2}$ (and those for $R_g^2$) are unclear compared to the case of LG model. One of the reasons for this is the continuous nature of the transition; the convergent speed of the simulation is very slow close to the transition point.

Next, we calculate the fractal dimension $D_f$, which is defined by

$$\begin{array}{c}\hfill R_g^2\sim N^{2/D_f}.\end{array}$$

(29)

To calculate this quantity at the transition point, we plot $R_g^2$ vs. N in Figure 9d–f in the log-log scale, where $R_g^2$ are obtained at the peak position of $C_{R_g}$ in Figure 8b,e,h. By fitting the plotted data to Equation (29), we obtain

$$\begin{array}{c}\hfill D_f=2.34\pm 1.01\left(\mathrm{LG}\right),D_f=2.54\pm 0.29\left(\mathrm{cano}\right),D_f=2.36\pm 0.60\left(\mathrm{modi}\right).\end{array}$$

(30)

Since the errors of $R_g^2$ are large in the first and third models (see Figure 9d,f), the errors in $D_f$ are also relatively large. The values of $D_f$ within these errors are comparable to $D_f=2.50\pm 0.30$ at the first-order transition point of the canonical model on the spherical lattice in Ref. [28].

To check the simulations are performed correctly, we can use the relations

$$\begin{array}{c}{S}_1^{\prime }/N=3/2\left(\mathrm{LG}\right),\hfill \\ S_1/N=3/2\left(\mathrm{cano},\mathrm{modi}\right),\hfill \end{array}$$

(31)

where $S_1^{\prime }$ of the first one is given in Equation (7), and $S_1^{\prime }=t{S}_1+\kappa S_2+2u{S}_3+2v{S}_4$. The second equality is also obtained by using almost the same technique for the first one. The symbol $\langle ·\rangle$ is omitted for these $S_1^{\prime }$ and $S_1$ for simplicity. The results plotted in Figure 9g–i are consistent with these predictions in Equation (31).

4. Summary and Conclusions

In this paper, we study the crumpling transition of a planar surface with a free boundary by parallel tempering Monte Carlo (PTMC) simulations using three different tethered lattice models; the Landau-Ginzburg (LG) model, the canonical (cano) model, and the modified canonical (modi) model, defined on triangulated fixed-connectivity lattices of disk topology. The order of the transition is first order in the LG model, while it is of second order in the other models. This second-order nature of the transition is consistent with the result reported in [22,23].

The models studied in this paper are not self-avoiding, however, the transition between the crumpled and smooth phases indicates that these states are stable against thermal fluctuations. This can give insights into studies on real materials such as graphene, where the crumpled states are expected to have many technological applications.

In the presence of impurity, a new phase is expected to appear, and the surface critical exponents are also expected to change from those of the model without impurity [55]. Thus, it is of great interest to study the planar surface model with impurities.