Measurements of magnetization on the Sierpiński carpet
MMeasurements of magnetization on the Sierpi´nski carpet
Jozef
Genzor , Andrej G endiar , and Tomotoshi N ishino Department of Physics, Graduate School of Science, Kobe University, Kobe 657-8501, Japan and Institute of Physics, Slovak Academy of Sciences,D´ubravsk´a cesta 9, SK-845 11, Bratislava, Slovakia (Dated: April 25, 2019)Phase transition of the classical Ising model on the Sierpi´nski carpet, which has the fractal di-mension log ≈ . T c = 1 . β associatedwith the local magnetization varies by two orders of magnitude, depending on lattice locations,whereas T c is not affected. I. INTRODUCTION
Understanding of phase transitions and critical phe-nomena plays an important role in the condensed matterphysics . Systems on regular lattices are of the majortarget of such studies, where elementary models exhibittranslationally invariant states, which are scale invariantat criticality. It has been known that critical behavioris controlled by global properties, such as dimensionalityand symmetries. This is the concept of the universality .If we focus our attention on inhomogeneous lattices,there is a group of fractal lattices, which are self-similarand exhibit non-integer Hausdorff dimensions. Geomet-rical details, such as lacunarity and connectivity, couldthus modify the properties of their critical phenomena.An important aspect of the fractal lattices is the ramifi-cation, which is the smallest number of bonds that haveto be cut in order to isolate an arbitrarily large boundedsubset surrounding a point. In the early studies by Gefen et al. , it was shown that the short-range classical spinmodels on finitely ramified lattices exhibit no phase tran-sition at nonzero temperature.The ferromagnetic Ising model on the fractal latticethat corresponds to the Sierpi´nski carpet is one of themost extensively studied models with fractal lattice ge-ometry. Monte Carlo studies combined with the finite-size scaling method have been performed , includingMonte Carlo renormalization group method . The criti-cal temperature T c is relatively well estimated within thenarrow range 1 . < ∼ T c < ∼ .
50, where one of the mostrecent estimate is T c = 1 . et al. . Onthe other hand, estimates of critical exponents are stillfluctuating, since it is rather hard to collect sufficientnumerical data for a precise finite-size scaling analysis .This is partially so because an elementary lattice unit cancontain too many sites, and there is a variety of choiceswith respect to boundary conditions. This situation per-sists even in a recent study by means of a path-countingapproach . Yet, a number of issues remains unresolvedconcerning uniformity of fractal systems in the thermo-dynamic limit .Recently, we show that the higher-order tensor renor-malization group (HOTRG) method can be used as an appropriate numerical tool for studies of certain types offractal systems . The method is based on the real-space renormalization group, and, therefore, the self-similar property of fractal lattices can be treated in anatural manner. In this article, we apply the HOTRGmethod to the Ising model on the fractal lattice that cor-responds to the Sierpi´nski carpet. The method enablesus to estimate T c from the temperature dependence ofthe entanglement entropy s ( T ). In order to check theuniformity in the thermodynamic functions, we choosethree distinct locations on the lattice, and calculate thelocal magnetization m ( T ) and the bond energy u ( T ). Asit is trivially expected, these local functions, m ( T ) and u ( T ), yield the identical T c . Contrary to the naive in-tuition, the critical exponent β , which is associated withthe local magnetization m ( T ) ∝ ( T c − T ) β , strongly de-pends on the location of measurement, and the estimatedexponent β can vary within two orders of magnitude withrespect to the three different locations on the fractal lat-tice, where the local functions are calculated.Structure of this article is as follows. In the next sec-tion, we explain the recursive construction of the fractallattice, and express the partition function of the systemin terms of contractions among tensors. In Sec. III weintroduce HOTRG method for the purpose of keepingthe numerical cost realistic. The way of measuring thelocal functions m ( T ) and u ( T ) is explained. Numericalresults are shown in Sec. IV. Position dependence on thelocal functions is observed. In the last section, we sum-marize the obtained results, and discuss the reason forthe pathological behavior of the fractal system. II. MODEL REPRESENTATION
There are several different types of discrete latticesthat can be identified as the Sierpi´nski carpet. Amongthem, we choose the one constructed by the extensionprocess shown in Fig. 1. In the first step ( n = 1), thereare eight spins in the unit, as it is shown on the left.The Ising spins σ = ± a r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r n = = = C (1) C (2) C (3) X (1) X (2) FIG. 1: Build-up process of a discrete analog of the Sierpi´nskicarpet. The circles represent the lattice points, where theIsing spins are located. The vertical and horizontal links de-note the interacting pairs. The first three units n = 1, 2, and3 are shown. For each unit n , we draw the corners C ( n ) bythe thick lines. We label the shaded regions X (1) and X (2) . ond step ( n = 2), the eight units are grouped to form anew extended unit, as shown in the middle. Now, thereare 64 spins on the 9 × n = 3). Generally,in the n -th step, an extended unit contains 8 n spins onthe 3 n × n lattice. The Hausdorff dimension of this lat-tice is d H = log ≈ . n → ∞ .In the series of the extended units we have thus con-structed, there is another type of the recursive structure.In Fig. 1 at the bottom of each unit, we have drawn apyramid-like area by the thick lines. One can identifyfour such pyramid-like areas within each unit (enumer-ated by n ), and each area can be called the corner C ( n ) .The corners are labeled C (1) , C (2) , and C (3) from left toright therein. It should be noted that there are only 2 n − spin sites in common, where two adjacent corners meet.In the case n = 2 drawn in the middle, we shaded aregion on the left, which contains six sites, and label theregion X (1) . Having observed the corner C (2) at the bot-tom, we found out that the corner consists of two rotatedpieces of X (1) and the four pieces of C (1) . In n = 3, weshaded a larger region X (2) (in the similar manner as X (1) ), which now contains 36 sites. We can recognizethat X (2) consists of seven pieces of X (1) and the twopieces of C (1) . We have thus identified the following re-cursive relations, which build up the fractal: • Each n -th unit contains 4 pieces of C ( n ) , • C ( n +1) contains 2 pieces of X ( n ) and 4 pieces of C ( n ) , • X ( n +1) contains 7 pieces of X ( n ) and 2 pieces of C ( n ) .The Hamiltonian of the Ising model, which is con-structed on the series of finite-size systems n = 1 , , , · · · ,has the form H ( n ) = − J (cid:88) (cid:104) ab (cid:105) σ a σ b . (1) ij C ij (1) σ σξ ij (1) ijkl σσσσ ηξ Xkl a C ad X cb b FIG. 2: Structure of the initial corner matrix C (1) ij in Eq. (3)and the 4-leg tensor X (1) ijkl in Eq. (6). The summation runs over all pairs of the nearest-neighbor Ising spins, as shown by the circles in Fig. 1.The spin positions are labeled by the lattice indices a and b . They are connected by the lines, which corre-spond to the ferromagnetic interaction J >
0, and noexternal magnetic field is imposed. First we calculatethe partition function (expressed in arbitrary step n ) Z ( n ) = (cid:88) exp (cid:20) − H ( n ) k B T (cid:21) (2)as a function of temperature T , where the summation istaken over all spin configurations, and where k B denotesthe Boltzmann constant. At initial step n = 1, we definethe corner matrix C (1) ij = (cid:88) ξ = ± exp (cid:2) Kξ ( σ a + σ b ) (cid:3) , (3)where K = J/k B T , and the matrix indices i = ( σ a + 1) / j = ( σ b + 1) / ξ isdenoted by the filled circle. We have chosen the orderingof the indices i and j , which is opposite if comparing C (1) ij with the corresponding graph. The partition function ofthe smallest unit ( n = 1), which contains 8-spins, is thenexpressed as Z (1) = (cid:88) ijkl C (1) ij C (1) jk C (1) kl C (1) li , (4)and can be abbreviated to Tr (cid:2) C (1) (cid:3) . We will express Z ( n ) for arbitrary n > Z ( n ) = Tr (cid:2) C ( n ) (cid:3) (5)by means of the corner matrix C ( n ) ij , where each one un-dergoes extensions , as we define in the following.Let us notice that the region X (1) appears from thestep n = 2. The Boltzmann weight corresponding to thisregion X (1) can be expressed by the 4-leg (order-4) tensor X (1) ijkl = (cid:88) ξη exp (cid:2) K ( σ a σ b + σ c σ d + ξη ) (cid:3) × exp (cid:2) Kξ ( σ a + σ d ) + Kη ( σ b + σ c ) (cid:3) , (6) i j j i j j i k l k a ec f hq s e p ri fd gacblb d CC CC X X CCX X XXXX X FIG. 3: Extension of the local matrix C ( n ) in Eq. (7) (onthe left) and the tensor X ( n ) in Eq. (8) (on the right). where the spin locations are depicted in Fig. 2 (bot-tom). We have additionally introduced new indices k = ( σ c + 1) / l = ( σ d + 1) /
2. Now we can math-ematically represent the recursive relations in terms ofcontractions among the matrices C ( n ) and tensors X ( n ) .Figure 3 shows the graphical representation of the exten-sion processes. Taking the contraction among the twotensors X ( n ) and the four matrices C ( n ) , as shown inFig. 3 (left), we obtain the extended corner matrix C ( n +1) through the corresponding formula C ( n +1) ij = C ( n +1)( i i )( j j ) = (cid:88) abcdef C ( n ) aj X ( n ) abcj C ( n ) fc C ( n ) db X ( n ) dei f C ( n ) i e , (7)where the new indices i and j , respectively, representthe grouped indices ( i i ) and ( j j ). Apparently, thediagram in Fig. 3 (left) is more convenient than Eq. (7)for the better understanding of the contraction geometry.This relation can be easily checked for the case n = 1after comparing Figs. 1, 2, and 3.Similarly, the extension process from X ( n ) to X ( n +1) shown in Fig. 3 (right) can be expressed by the formula X ( n +1) ijkl = X ( n +1)( i i )( j j )( k k )( l l ) = (cid:88) abcdefghprqs X ( n ) abl p X ( n ) bck l X ( n ) cdqk X ( n ) fgda X ( n ) efri X ( n ) ghj s X ( n ) i j he C ( n ) rp C ( n ) sq , (8)where we have again abbreviated the grouped indices to i = ( i i ), j = ( j j ), k = ( k k ), and l = ( l l ). Thisrelation can be checked for the case n = 1 by comparingthe area X (1) and X (2) in Fig. 1.Through the iterative extension of the tensors, we can formally obtain the corner matrix C ( n ) ij for arbitrary n ,and express Z ( n ) by Eq. (5). The free energy per spin isthen f ( n ) = − n k B T ln Z ( n ) (9)since the n -th unit contains 8 n spins. This functionconverges to a value f ( ∞ ) in the thermodynamic limit n → ∞ , where convergence with respect to n is rapid,and n = 35 is sufficient in the numerical analyses. Thespecific heat per site can be evaluated by taking the sec-ond derivative of the free energy c f ( T ) = − T ∂ ∂T f ( ∞ ) . k i j j i j j i k l k a ec f ijhq s jel p r i ifd gacbl b d CC CC X XCCX X XXXX X
FIG. 4: (Color online) The renormalization group transfor-mations in Eq. (11) (on the left) and Eq. (12) (on the right)applied, respectively, to Eq. (8) and (7) (cf. Fig. 3).
III. RENORMALIZATION GROUPTRANSFORMATION
The matrix dimension of C ( n ) is 2 n − by definition.Therefore, it is impossible to keep all the matrix elementsfaithfully in numerical analysis, when n is large. The sit-uation is severer for X ( n ) , which has four indices. Bymeans of the HOTRG method , it is possible to reducethe tensor-leg dimension, the degree of freedom, down toa realistic number. The reduction process is performedby the renormalization group transformation U , which iscreated from the higher-order singular value decomposi-tion (SVD) applied to the extended tensor X ( n +1) ijkl .Suppose that the tensor-leg dimension in X ( n ) ijkl is D for each index, i.e., i, j, k, l = 0 , , . . . , D −
1. As wehave shown in Eq. (8), the dimension of the grouped in-dex i = ( i i ) in X ( n +1)( i i )( j j )( k k )( l l ) is equal to D .We reshape the four tensor indices to form a rectangularmatrix with the grouped index ( i i ) and the remaininggrouped index ( j j k k l l ) with the dimension D . Ap-plying the singular value decomposition to the reshapedtensor, we obtain X ( n +1)( i i )( j j k k l l ) = (cid:88) ξ U ( i i ) ξ ω ξ V ( j j k k l l ) ξ , (10)where U and V are generalized unitary, i.e. orthonormal,matrices U T U = V T V = . We assume the decreasingorder for the singular values ω ξ by convention. Keeping D dominant degrees of dominant freedom for the index ξ at most, we regard the matrix U ( i i ) ξ as the renor-malization group (RG) transformation from ( i i ) to therenormalized index ξ . For the purpose of clarifying therelation between the original pair of indices ( i i ) andthe renormalized index ξ , we rename ξ to i and write theRG transformation as U ( i i ) i . In the same manner, weobtain U ( j j ) j , U ( k k ) k , and U ( l l ) l , where we have dis-tinguished the transformation matrices by their indices.The RG transformation is then performed as X ( n +1) ijkl ← (cid:88) i i j j k k l l U ( i i ) i U ( j j ) j U ( k k ) k U ( l l ) l X ( n +1)( i i )( j j )( k k )( l l ) , (11) ABY
FIG. 5: (Color online) Three positions A (on inner boundary),B (on outer boundary), and Y (in innermost position) chosenfor the observation of local functions. The lower half of theunit n = 4 is drawn only. where the sum is taken over the indices on the con-nected lines in Fig. 4 (left). The left arrow used inEq. (11) represents the replacement of the expandedtensor X ( n +1)( i i )( j j )( k k )( l l ) for the renormalized one X ( n +1) ijkl . Since the RG transformation matrices U are ob-tained from SVD applied to X ( n +1)( i i )( j j )( k k )( l l ) , thereis no guarantee that the RG transformation can bestraightforwardly applied to C ( n +1)( i i )( j j ) , as we have de-fined in Eq. (7). It has been confirmed that the transfor-mation C ( n +1) ij ← (cid:88) i i j j U ( i i ) i U ( j j ) j C ( n +1)( i i )( j j ) (12)is of use in the actual numerical calculation. The corre-sponding diagram is shown in Fig. 4 (right).We add a remark on the choice of the transformationmatrix U . In a trial calculation, once we tried to cre-ate U from the corner matrix C ( n +1) ij by both SVD anddiagonalization. However, we encountered numerical in-stabilities, in which the singular values (or eigenvalues)decayed to zero too rapidly, especially, when n was large.Thus, we always create U from SVD that is applied to X ( n +1) ijkl only.With the use of these RG transformations, it is possi-ble to repeat the extension processes in Eq. (7) and (8),and to obtain a good numerical estimate for Z ( n ) and f ( n ) in Eq. (9). The actual numerical calculations in thiswork were performed by a slightly modified procedure.We split X ( n ) ijkl into two halves and represent each partby 3-leg tensor. This computational trick allowed us toincrease the leg-dimension up to D = 28, or even larger. A. Impurity tensors
In the framework of the HOTRG method, thermody-namic functions, such as the magnetization per site m ( T )and the internal energy per bond u ( T ), can be calcu-lated from the free energy per site f ( ∞ ) . Alternatively, i i j j i j j k k l i l i i j j a b d ec f ijc hq s jk d ga fb el p r ia b d ec f ij C C CXX CX C B CC XXX X XXCXA Y
FIG. 6: (Color online) Extension processes of the impuritytensor in Eq. (15) (top left), Eq. (16) (top right), and Eq. (17)(bottom). The RG transformations U are expressed by theexternal three green lines meeting at the open circles. these functions are obtained by inserting impurity ten-sors (separately derived from C ( n ) and X ( n ) ) into thetensor network of the entire system. Since the fractallattice under consideration is inhomogeneous, these ther-modynamic functions can depend on the position theyare placed. In order to check the dependence, we choosethree typical locations A , B , and Y , as shown in Fig. 5on the fractal lattice.As an example of such a single site function, let us con-sider a tensor representation of the local magnetization.Looking at the position of site A in Fig. 5, one finds thatit is located on the corner matrix C (1) . Thus, the initialimpurity tensor on that location is expressed as A (1) ij = (cid:88) ξ = ± ξ exp (cid:2) K ξ (cid:0) σ i + σ j (cid:1)(cid:3) , (13)similar to Eq. (3). It is also easy to check that the ini-tial impurity tensor B (1) , which is placed on a positiondifferent from A , is expressed by the identical equation,so that we have A (1) ij = B (1) ij . The site Y lies inside thearea X (1) and we define the corresponding initial tensorfor local magnetization as Y (1) ijkl = (cid:88) ξη ξ + η (cid:2) K ( σ i σ j + σ k σ l + ξη ) (cid:3) × exp (cid:2) K ξ (cid:0) σ j + σ k (cid:1) + K η ( σ i + σ l ) (cid:3) , (14)similarly to Eq. (6).We can thus build up analogous extension processesof tensors, each of which contains an impurity tensor wehave defined. The extension process of the impurity cor-ner matrix that contains A (1) ij is then written as A ( n +1)( i i )( j j ) = (cid:88) abcdef C ( n ) aj X ( n ) abcj A ( n ) fc C ( n ) db X ( n ) dei f C ( n ) i e , (15)which is graphically shown in Fig. 6 (top left). Therein,the RG transformation A ( n +1)( i i )( j j ) → A ( n +1) ij is depicted T s T c FIG. 7: The entanglement entropy s ( T ) in Eq. (21). by the green lines with the open circles, which stand for U in accord with Eq. (12). The impurity tensor placedaround the site B obeys the extension procedure B ( n +1)( i i )( j j ) = (cid:88) abcdef C ( n ) aj X ( n ) abcj C ( n ) fc B ( n ) db X ( n ) dei f C ( n ) i e , (16)as shown on the top right of Fig. 6 (top right). For thelocation Y shown in Fig. 5, we take the contraction Y ( n +1)( i i )( j j )( k k )( l l ) = (cid:88) abcdefghprqs X ( n ) abl p X ( n ) bck l X ( n ) cdqk Y ( n ) fgda X ( n ) efri X ( n ) ghj s X ( n ) i j he C ( n ) rp C ( n ) sq , (17)which is depicted in Fig. 6 (bottom), where the graph isrotated by the right angle for book keeping.In the calculation of the local bond energy u ( T ), theinitial tensors satisfy the equations A (1) ij = − J (cid:0) σ i + σ j (cid:1) (cid:88) ξ ξ exp (cid:2) Kξ (cid:0) σ i + σ j (cid:1)(cid:3) , (18) Y (1) ijkl = (cid:88) ξη − Jξη exp (cid:2) K ( σ i σ j + σ k σ l + ξη ) (cid:3) × exp (cid:2) Kξ (cid:0) σ j + σ k (cid:1) + η ( σ i + σ l ) (cid:3) , (19)recalling that B (1) ij = A (1) ij . Starting the extension pro-cesses with these initial tensors, we can calculate the ex-pectation value of the bond energy around the site A bymeans of the ratio u A ( T ) = lim n →∞ Tr (cid:16) A ( n ) (cid:2) C ( n ) (cid:3) (cid:17) Tr (cid:16)(cid:2) C ( n ) (cid:3) (cid:17) . (20)The convergence with respect to n is fast because of thefractal geometry. It is straightforward to obtain the localenergy u B ( T ) and u Y ( T ), as well as the local magnetiza-tion m A ( T ), m B ( T ), and m Y ( T ) in the same manner. T c c Y = ∂ u Y / ∂ Tc A = ∂ u A / ∂ Tc B = ∂ u B / ∂ Tc f = − T ∂ f / ∂ T T -60-40-20020 ∂ c / ∂ T ≡ ∂ T c c f c Y c B c A ∂ T c A ∂ T c B T c T c FIG. 8: (Color online) Specific heats c A ( T ), c B ( T ), c Y ( T ),and c f ( T ). The inset shows the derivative of the specific heatwith respect to temperature ∂ T c . IV. NUMERICAL RESULTS
For simplicity, we use the temperature scale with k B =1 and fix the ferromagnetic interaction strength to J = 1.All the shown data are obtained after taking a sufficientlylarge number of system extensions, provided that theconvergence with respect to n has been reached. Thedegrees of freedom D for each leg-dimension is D = 28at most. Apart from the critical (phase transition) re-gion, where D needs to be the largest, we used D = 18,which sufficed to obtain precise and converged data wehave used for drawing all the graphs.An analogous kind of the entanglement entropy s ( T )can be calculated by the HOTRG method. After apply-ing SVD to the extended tensor, s ( T ) can be naturallyobtained from the singular values ω ξ in Eq. (10) throughthe formula s ( T ) = − (cid:88) ξ ω ξ Ω ln ω ξ Ω , (21)where Ω = (cid:80) ξ ω ξ normalizes the probability. The entan-glement entropy s ( T ) always exhibits stable convergencewith respect to n . Figure 7 shows the temperature de-pendence of s ( T ), which is obtained with D = 18. Thereis a sharp peak at the critical temperature T c , which canbe roughly determined as 1 .
48 from the data shown.Taking the numerical derivative with respect to T forthe calculated local energies u A ( T ), u B ( T ), and u Y ( T ),respectively, we obtain the specific heats c A ( T ), c B ( T ),and c Y ( T ), as shown Fig. 8. We observe a sharp peak in c Y ( T ) at T c , whereas there is only a rounded maximumin c A ( T ) and c B ( T ), and their peak positions do not coin-cide with T c associated with the position Y . The specificheat per site c f ( T ) defined in Sec. II as well as c A ( T ) and c B ( T ) demonstrate a weak singularity at T c . This factcan be confirmed by taking their derivative with respect T m m Y m A m B T c FIG. 9: (Color online) The local magnetization m Y ( T ), m A ( T ), and m B ( T ). to T , i.e., ∂c∂T , which leads to the identical singularity at T c , as shown in the inset of Fig. 8. The result clearlymanifests that the critical behavior strongly depends onthe location, where the measurements of the bond energyis carried out.Figure 9 shows the local magnetizations m A ( T ), m B ( T ), and m Y ( T ) with respect to temperature T , underthe condition that D = 18. They fall to zero simultane-ously at the identical T c , while the critical exponent β in m ( T ) ∝ ( T c − T ) β is significantly different for each case.From the plotted m A ( T ) we obtain β = 0 .
52, and from m B ( T ) we obtain β = 0 .
78. In both cases we use therough estimate T c = 1 . | T c − T | < .
015 are considered for numerical fitting.Since the variation in m Y ( T ) is too rapid to capture β under the condition D = 18, we increase the tensor-legfreedom up to D = 28. Figure 10 shows m Y ( T ) zoomed-in around T ∼ . | T − T c | < ∼ − . Therefore, the data points inthis narrow region were excluded from the fitting anal-ysis. Then, we obtain T c = 1 . β = 0 . β obtained from m A ( T ) and m B ( T ). In the similar manner, as we have observed forthe specific heat, critical behavior of the model stronglydepends on the location of the impurity tensors A , B ,and Y on the Sierpi´nski carpet. V. CONCLUSIONS AND DISCUSSIONS
We have investigated the phase transition of the ferro-magnetic Ising model on the Sierpi´nski carpet. The nu-merical procedures in the HOTRG method are modified,so that they fit the recursive structure in the fractal lat-tice. We have confirmed the presence of the second orderphase transition, which is located around T c = 1 . T m Y D=20D=22D=24D=26D=28 × -4 × -4 T c − T × -11 × -11 × -11 × -11 m Y β − FIG. 10: (Color online) Detailed view of m Y ( T ) when D =20 , , ,
26, and 28. Inset: the power-law behavior below T c = 1 . β = 0 . in accordance with the previous studies . The globalbehavior of the entire system captured by the free en-ergy per site f ( ∞ ) exhibits the presence of a very weaksingularity at T c , as we observed in Ref. 14.What is characteristic of this fractal lattice is the po-sition dependence in the local magnetization m ( T ) andlocal energy u ( T ). For example, we find that the criticalexponent β differs by a couple of orders of magnitude,which corresponds to the fact that the measured magne-tization depends on position, where the impurity tensoris placed on the fractal-lattice. A key feature appears inthe local energy u Y ( T ), where we deduce a sharp peak inits temperature derivative, c Y ( T ); contrary to the smoothbehavior in c f ( T ), being the averaged specific heat. In-tuitively, such position dependence would be explainedby the density of sites around the pinpointed location.Around the site Y, the spins are interconnected moredensely than those around the boundary sites A and Bin Fig. 5. One might find a similarity with the criticalbehavior on the Bethe lattice , where the singular be-havior is only visible deep inside the system, whereas thefree energy is represented by an analytic function of T for the entire lattice.The current study can be extended to other fractal lat-tices, e.g., variants of the Sierpi´nski carpet or to a fractallattice, we had already studied earlier , where the posi-tional dependence of the impurities has not been exam-ined yet. Another point to be considered is to investigatemore variations of the locations on the fractal lattice inorder to analyze the mechanism of the non-trivial behav-ior observed in the current position dependence. Acknowledgments
This work was partially funded by Agent´ura pre Pod-poru V´yskumu a V´yvoja (APVV-16-0186 EXSES) andVedeck´a Grantov´a Agent´ura MˇSVVaˇS SR a SAV (VEGAGrant No. 2/0123/19). T.N. and A.G. acknowledgethe support of Ministry of Education, Culture, Sports,Science and Technology (Grant-in-Aid for Scientific Re-search JSPS KAKENHI 17K05578 and 17F17750). J.G.was supported by Japan Society for the Promotion of Science (P17750). T.N. thanks the funding of the Grant-in-Aid for Scientific Research MEXT Exploratory Chal-lenge on Post-K computer (Frontiers of Basic Science:Challenging the Limits). Phase transitions and critical phenomena , vol. 1-20, ed.C. Domb, M.S. Green, and J. Lebowitz (Academic Press,1972-2001). Y. Gefen, B.B. Mandelbrot, and A. Aharony, Phys. Rev.Lett. , 855-858 (1980). Y. Gefen, Y. Meir, B.B. Mandelbrot, and A. Aharony,Phys. Rev. Lett. , 145-148 (1983). Y. Gefen, A. Aharony, and B.B. Mandelbrot, J. Phys. A:Math. Gen. , 1267-1278 (1983). Y. Gefen, A. Aharony, and B.B. Mandelbrot, J. Phys. A:Math. Gen. , 1277-1289 (1984). J.M. Carmona, U.M.B. Marconi, J.J. Ruiz-Lorenzo,A. Taranc´on, Phys. Rev. B , 14387 (1998). P. Monceau, M. Perreau, F. H´ebert, Phys. Rev. B , 6386(1998). P. Monceau, M. Perreau, Phys. Rev. B , 184420 (2001). G. Pruessner, D. Loison, and K.D. Schotte, Phys. Rev. B , 134414 (2001). M.A. Bab, G. Fabricus, and E. Albano, Phys. Rev. E ,036139 (2005). Pai-Yi Hsiao, P. Monceau, Phys. Rev. B , 064411 (2003). T.W. Burkhardt and J.M.J. van Leeuwen, Real-spacerenormalization, Topics in Current Physics 30 (Springer,Berlin, 1982), and references therein. M. Perreau, Phys. Rev. B , 174407 (2017), and refer-ences therein. J. Genzor, A. Gendiar, and T. Nishino, Phys. Rev. E ,012141 (2016). R. Krcmar, J. Genzor, Y. Lee, H. ˇCenˇcarikovˇa, T. Nishino,and A. Gendiar, Phys. Rev. E , 062114 (2018). Z.Y. Xie, J. Chen, M.P. Qin, J.W. Zhu, L.P. Yang, andT. Xiang, Phys. Rev. B , 045139 (2012). L. de Lathauwer, B. de Moor, J. Vandewalle, SIAM J.Matrix Anal. Appl. , 1324 (2000). R. J. Baxter, Exactly Solved Models in Statistical Mechan-ics (Academic Press, London, 1982). R. Krcmar, A. Gendiar, K. Ueda, and T. Nishino, J. Phys.A: Math. Theor.41