{'fraction_non_alphanumeric': 0.04443568001051732, 'fraction_numerical': 0.015469225881373387, 'mean_word_length': 3.965727341964966, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sine-square deformation applied to classical Ising models

Sine-square deformation (SSD) is a treatment proposed in quantum systems, which spatially modifies a Hamiltonian, gradually decreasing the local energy scale from the center of the system toward the edges by a sine-squared envelope function. It is known to serve as a good boundary condition as well as to provide physical quantities reproducing those of the infinite-size systems. We apply the SSD to one- and two-dimensional classical Ising models. Based on the analytical calculations and Monte Carlo simulations, we find that the classical SSD system is regarded as an extended canonical ensemble of a local subsystem each characterized by its own effective temperature. This effective temperature is defined by normalizing the system temperature by the deformed local energy scale. A single calculation for a fixed system temperature provides a set of physical quantities of various temperatures that quantitatively reproduce well those of the uniform system.


I. INTRODUCTION
The Hamiltonian in condensed matter physics is spatially uniform in most cases, and its symmetry determines the physical properties of the system. For this reason, deforming the Hamiltonian may usually mean modifying the physical state itself. However, this turned out not always to be the case for a series of operation called "sine-square deformation" (SSD). The deformed Hamiltonian, H SSD , is generated from the original Hamiltonian, H = drH(r), using an envelope function, f SSD (r), as where r is a coordinate of the system with its origin at the center, and f SSD (r) = 1 2 1 + cos πr R , with R being half the system length. As shown in Fig.1, the sine-square function, f SSD , governs the whole system by a single wavelength of 2R, and f SSD (r)H(r) varies smoothly from the maximum value at the center toward the edges. The SSD was proposed as one of the smooth boundary conditions to remove boundary effects 1,2 , e.g. Friedel oscillations from the open boundaries, or artificial potentials that emerges for a chosen cluster sizes and shapes which may stabilize fictitious orders 3 . The SSD Hamiltonian is empirically known to generate the wave function that recovers the translational symmetry, in perfect coincidence with the wave function under the periodic boundary condition (PBC) 2,4 . This coincidence is proved in an XY model and a transverse-field Ising chain 5,6 . The SSD ground state is also shown to serve as a restricted class of conformal field theories with some applications to a wider class of conformal mappings [7][8][9][10][11][12] , and in that context the quantum dynamics of the SSD system is tested [13][14][15][16][17] . When an external field is applied, the physical quantities as a response function to the field is evaluated from the SSD Hamiltonian, whose accuracy against the exact solution in the thermodynamic limit is O(10 −4 ) 18 . From these studies, it is found that the SSD Hamiltonian loses its translational symmetry but could still or better reproduce the physical properties of the original Hamiltonian. A key to understanding this phenomenon is the idea of real-space renormalization 19 , which stems from Wilson's poor man's scaling 20 . In a system with translational symmetry, a quantum state is characterized by a wave number. When the SSD is imposed, a wavenumber is no longer a quantum number, and a scattering term generated by f SSD (r) mixes this original eigen states 21 . Such mixing generates a series of localized wave packets that serves as another basis set of the Hamiltonian. Since these states are localized, they are no longer influenced by the system size nor by the boundary. Therefore, it allows us to obtain physical quantities that reproduce those in the thermodynamic limit 22 .
In the present paper, we apply this SSD to classical Ising models, as they provide a good platform to test approximate methods based on their exactly solvable structures 23 . The Hamiltonian of the uniform classical Ising model is, with the uniform coupling constant J. The index γ represents a lattice site and the summation is taken over all the neighboring pairs of spins. The exact solutions are known for the one-dimensional(1D) chain and twodimensional(2D) square lattice. A second-order phase transition occurs in the 2D square lattice model at the temperature, k B T c = 2/ ln(1 + √ 2) ∼ 2.2692J.
In constructing the SSD Ising Hamiltonian, we replace the uniform interaction J with f SSD (r)J, where the position vector r is defined at the center of each bond. Suppose that the temperature of this deformed system is k B T . Then, the system can be regarded as an assemblage of Ising spins with "uniform" interactions at a renormalized "effective temperature", k B T eff ≡ k B T /f SSD (r). Figure 1(a) shows an effective temperature as a function of r; the minimum value, k B T eff (r = 0) = k B T , at the center gradually increases toward infinity at the system edge.
The deformation in 1D chain is straightforward. For the site index, i = 1, · · · , L, in Fig.1(c), the bond connecting the i and (i+1)-th sites is located at r i = i−L/2, and we set R = L/2 to generate the values f SSD (r i ) in Eq. (2).
For 2D square lattice, we consider an L × L lattice shown in Fig. 1(b) to keep the aspect ratio as unity 24 , and define a coupling along the row between sites (i, j) and (i + 1, j) as J 1;ij , and its location is defined by r 1;ij . We take a bond along the column between sites (i, j) and (i, j + 1) as J 2;ij located at r 2;ij . Here, the vectors r 1;ij , r 2;ij are not the ordinary position vectors but are introduced to define the deformation function along the axis r parallel to the row and column for Case (i) and (ii), respectively, as shown in Fig. 1(c). We consider two cases: Case (i) deforms the bond interaction only along the row direction, keeping those along the column uniform; we plug in r 1;ij = i − L/2 and r 2;ij = i − (L − 1)/2 to Eq.(2) with R = L/2. The second one, Case (ii) deforms the interaction along the column, while keeping the row direction uniform; we take r 1;ij = j − (L − 1)/2 and r 2;ij = j −L/2 for this case. The two cases formally differ in the analytical treatment as we show in §.III, where we construct the column-to-column transfer matrix.
The aim of this paper is to clarify the role of SSD in a classical Ising model. We show that the energy and related quantities accurately reproduce those of k B T eff , which means that one can obtain a set of data for a wide temperature range in a single calculation. The physical implication is that the SSD system is an assemblage of local subsystems with different temperatures, which form a modified canonical ensemble. The neighboring local subsystems have similar effective temperatures and work to each other as a heat bath. In the final part of the manuscript, we address the possibility of taking other types of deformation.
The paper is organized as follows; In §.II we demonstrate that the transfer matrix method is exactly applied to the SSD Ising model in 1D. In §.III, we analyze the 2D Ising model using a fermionic representation, and obtain an exact form of the partition function for a finite system size for both Case (i) and Case (ii). Then, we numerically evaluate the bond energy using these formulations in §IV. We also perform a classical Monte Carlo simulation for the SSD Hamiltonian in §IV C to test the numerical applicability of SSD. §V gives a summary and discussion.  1. (a) Change in effective temperature, kBT eff , from the center (r = 0) toward the edge (r/R = 1), and the corresponding fSSD(r). (b) Two-dimensional square lattice with interactions J1,ij, J2,ij on bonds running in the positive directions of x-and y-axes from site γ = (i, j). The transfer matrix Vi is constructed in a unit of L sites along the y-direction. (c) SSD along row (Case (i)) and column (Case (ii)).

II. EXACT SOLUTION OF ONE DIMENSIONAL SSD ISING MODEL
Let us first consider a partition function of a 1D periodic lattice consisting of L sites, using the conventional notation, Here, r i = i − L/2 is the location of bond connecting site i and i + 1. The eigenvalues of the transfer matrix V i are, λ ± i = e Ki ± e −Ki , which explicitly depend on index-i, whereas, the corresponding eigen vectors, p ± i = (1, ±1)/ √ 2, are site-independent. In this way, all the transfer matrices are simultaneously diagonalized and the partition function is exactly given as Then, the exact expectation value of the bond energy is obtained as Taking the limit of open boundary, J L → 0, we find λ + L → 2 and λ − L → 0, and the bond energy σ i σ i+1 converges to Reminding that energy per bond in the uniform system takes a form, e bulk = − tanh(J/k B T ), in the thermodynamic limit, one finds that the site-dependent bond energy of the SSD Hamiltonian is e bulk at their "local" effective temperatures, k B T eff (r i ). The partition function at L → ∞ is given as Z = i λ + i which matches exactly the partition function of a system consisting of Ising bond degrees of freedom λ + i with interaction J and i-dependent effective temperature k B T eff (r i ).
For later convenience, we show that Eq.(4) is written using the Pauli matrix, τ ν , ν = x, y, z (we use τ instead of σ to avoid confusion) and unit matrix I and the pa- where we used a relation sinh 2K i sinh 2K * i = 1.

III. EXACT SOLUTION OF TWO DIMENSIONAL SSD ISING MODEL
In this section, we expand an analytical formula to compute the partition function Z of the SSD Ising model at a finite system size in two dimensions. The results shown here are also applied to deformations with functional forms other than f SSD . Among several different approaches 23,25-32 , we build our work based on the analysis given by Schultz, Mattis, and Lieb 32 . It provides a good description of the exact transfer matrix of the 2D Ising model by the 1D free fermionic degrees of freedom. In §.III A, we give a review of this work originally applied to the uniform Hamiltonian and show that it can be applied to our system where the interactions are sitedependent. We introduce the operators for the columnto-column transfer matrix in the fermionic representation and the related trace formula. §.III B is devoted to derivations of Z for Cases (i) and (ii) with the SSD Hamiltonian, which uses the formula obtained in §.III A.
A. Preliminaries

Column-to-column transfer matrix
The explicit form of a transfer matrix, V i , from column-i to column-(i+1) is given for a spatially nonuniform Hamiltonian as in the following. Let σ i ≡ {σ i,j |j = 1, 2, · · · L} be a complete orthonormal basis of dimensions 2 L of the ith column, where σ i,j = ±1 is the Ising degrees of freedom on site (i, j) in Fig.1(b). The following operators,V 1;i andV 2;i , include the interactions along ith column (J 1;ij ) and jth row (J 2;ij ), respectively, and give a partition function Z as Generalization of Eq.(8) to L degrees of freedom along the ith column immediately gives a description of V 1;i and V 2;i as with K 1;i,j = J 1,ij /k B T and K 2;i,j = J 2,ij /k B T , and K * 1;ij is defined as tanh K 1;ij = e −2K * 1;ij . We use a 2 L ×2 L matrix defined as a direct product, τ ν j = I ⊗ I ⊗ · · · ⊗ τ ν ⊗ I · · · ⊗ I, with τ ν operating on the jth Ising degrees of freedom on the column. This form is further transformed first by rotating the axis of the Pauli matrices as τ x j → −τ z j and τ z j → τ x j , and then by using a set of Pauli operators {τ ν j } operating on jth site.V where the trace of the operators over the 2 N -Hilbert space gives Z = Tr L i=1V 1;iV2;i+1 .
By making use of the Jordan-Wigner transformation, we obtain a fermionic representation of the operators aŝ where c † j (c j ) is a creation(annihilation) operator of spinless fermion. In the standard approach, the partition function is given as Z = Tr( (15) being an operator representing the column-(i)-tocolumn-(i + 1) transfer matrix. IfV i does not depend on i, one is able to diagonalize the 2 L × 2 L representation ofV i , and the product of the largest eigenvalues will give Z. However, in Case (i), since the transfer matrices depend on column-i, they cannot be diagonalized simultaneously, and this approach does not straightforwardly apply. In the next subsection, we review the derivation of the exact solution for the uniform Ising model. The formula Eqs. (16)(17)(18)(19)(20) will be adopted to calculate the partition function of Case (i) in the later section. The eigenvalue solution of the transfer matrix of the uniform system is also to be compared with those obtained for the SSD Hamiltonian.

Spatially uniform 2D Ising model
For the interaction parameters, K 1;i,j and K 2;i,j , let us omit the row index-j while keep the column indexi to clarify that they do not depend on j. Since the Hamiltonian is uniform along the column, the transfer matrix on the ith column is block diagonalized by using the Fourier transform of fermionic operators along the j-direction, where we take q = (2l − 1)π/L and 2πl/L when the number of fermions, N = q η † q η q , is even and odd, respectively. The even-and odd-N sectors originate from the constraint imposed on Eq.(12) due to anti-periodic (APBC) and periodic(PBC) boundary conditions, respectively. The operators in Eqs. (13) and (14) are rewritten as, where the last two terms ofĤ 2;i with wave vectors q = 0, π are present only in the odd-N sector. As η q 's with different |q| commute, one can decompose the exponentials of the transfer matrix aŝ We now prepare a matrix representation of Fock operators eĤ 2;iq /2 and eĤ 1;iq by one-body states, .
The formulation in the rest of this subsection holds only when H 2;iq does not depend on column index, i; it is omitted from the interaction parameters as K 1 and K 2 . We multiply e H2;iq/2 = I cosh K 2 + (cos qτ z + sin qτ x ) sinh K 2 and e H1;iq = I cosh 2K * 1 − τ z sinh 2K * 1 using the Pauli matrices τ ν , and find a real symmetric matrix e H 2;iq 2 e H1;iq e where in the final term we diagonalized the matrix by an orthogonal matrix P , and By the same matrix P as in Eq.(21), the operator η q undergoes a Bogoliubov transformation, Using this transformation, the product part of Eq. (19) is rewritten as eĤ i , where the Fock operatorĤ i represent-ing the transfer matrix is given aŝ Here, one can summarize all the q-terms by reading off η 0 = ξ 0 and η π = ξ π , since we see from Eq.(22) that . The partition function is obtained (for reference, see Eq. (32) in the next subsection) as The bond energy is written as Considering only the contributions from the largest term in the product of the second equation in Eq. (25), which is valid for L → ∞, we find which reproduces the result in Ref. [32]. This corresponds to the vacuum state of the Bogoliubov quasi-particle, i.e. having ξ † q ξ q = 0 for all-q. In this limit, Eqs. (26) and (27) are modified to those taking tanh(L 1 q ) → 1. Figure 2 shows pairs of ± (u) q at several different temperature. Here, by introducing a hole creation operator as an anihilation of particle,ξ † q = ξ q , Eq.(24) can be rewritten asĤ where the lowest energy eigen state ofĤ i is obtained by fully occupying a hole band and by keeping particle bands empty. This is a "vacuum"state of Bogoliubov quasiparticle. Exciting a Bogoliubov quasi-particle ξ † q ξ q = 0 corresponds to creating a particle-hole pair with the excitation energy 2 (u) q at wave number q. The gap at q = 0 closes at the transition temperature of the uniform 2D Ising model, T c = 2.2692.

Full counting statistics and the trace formula
To count the full statistics Z consisting of columndependent transfer matrices, we introduce another basic formula proved by Klich in Ref. [33]. Consider a second quantized single-particle operator acting on the Fock space asΓ where A is the M × M matrix representation A lm = l|Â|m of operatorÂ on a single particle Hilbert space of spinless fermions {|l } = {c † l |0 }, with a creation operator, c † l (l = 1 · · · M ), applied on a vacuum |0 . For two operators,Â andB, we find [Γ(Â),Γ(B)] =Γ([Â,B]). We know from the Baker-Campbell-Hausdorff formula that for given matrices A and B, there is a matrix C that fulfills e A e B = e C 35 . Then, the corresponding Fock operators are also given as As we see shortly, this relationship guarantees that one can rewrite Eq.(15) in a single exponential form eĤ i as far as the operatorsĤ 1;i ,Ĥ 2;i ,Ĥ 2;i+1 are written by the common single-particle basis. Any matrix C can be transformed to P −1 CP = diag(ξ 1 , ξ 2 , · · · ξ M ) + D with D being an upper triangular matrix, where we are familiar with D = 0 for the symmet- (1 + e ξm ) This trace formula holds for more than two products of exponentials of the operators. We apply this formula in obtaining Z for Case (i).
B. Exact solutions of the deformed 2D Ising models 1. Case (i): when the Hamiltonian is uniform along the column and non-uniform along the row In this subsection, the parameters K 1;i,j and K 2;i,j depend on i but not on j, where we omit index-j for simplicity. The formula Eqs. (16)-(20) still holds. However, since the transfer matrix V i depends on column-index i, and since V i is no longer symmetric, we cannot straightforwardly prepare an orthogonal matrix P that diagonalizes all V i 's simultaneously.
Instead of dividing Z into columns, we first take the whole product over the columns for each q to obtainV q , and then combine all q sectors as We again find a final single exponential form of the Fock operatorĤ q sinceĤ 1;iq andĤ 2;iq fulfill the condition for Γ(Â) andΓ(B) in Eq. (31). By multiplying the 2 × 2 matrix in Eq. (20), we obtain an explicit form This matrix is diagonalized to diag(e −2Eq , e 2Eq ) by performing a Bogoliubov transformation (η † q , η −q ) → (ξ † q ,ξ −q ) similarly to Eq.(23), and we find a final form Here, E q ≥ 0 is an order-L quantity, which is the energy carried by the Bogoliubov quasi-particleξ q . In obtaining E q , fermions, η q and η −q , mix for 0 < q < π so that the operator in Eq. (35) is represented by a 2 × 2 matrix. Whereas for q = 0 and π, there is no mixing and we The partition function is obtained as (36) Among the contributions from Bogoliubov quasi-particles to Z, namely the last product term in the above equation, a so-called largest eigenvalue is obtained solely from a "vacuum" state of Bogoliubov quasi-particle as and for L → ∞ we obtain Z = Λ 0 . Exciting a single Bogoliubov quasi-particle with a minimum excitation energy, 2min(E q ), yields the next-largest eigenvalue, so that Λ 1 /Λ 0 = e −2min(Eq) . The rest of the eigenvalues are determined by successively exciting quasi-particles The energy densities per bond along the column and row are formally given as In numerically evaluating these quantities, the derivatives are much less accurate than those we obtain for Case (ii) in the next section.

Case (ii): when the Hamiltonian is non-uniform along the column and uniform along the row
In this subsection, we consider Case (ii). We start from Eqs.(13)- (15). From the discussions given in §III A 3, the transfer matrix operator in Eq.(15) is rewritten asV i = eĤ 2;i/2 eĤ 1;i eĤ 2;i/2 = eĤ i , whereĤ i takes a quadratic form of a set of one-body operators {c i }. SinceV i 's do not depend on a column-index i, their representations are separately diagonalized simultaneously for all columns. This time, however, the Hamiltonian is nonuniform along the column. Then, the L×L matrix representation ofĤ i can no longer be block diagonalized into smallest pieces by the Fourier transformation, nor can we apply a simple Pauli matrix representation we used in obtaining Eq. (21). Instead, the form ofĤ i is obtained through the following processes.
We first describeĤ 1;i andĤ 2;i /2 in a quadric form of the one-body operators. Since the number of fermions does not conserve, we need to prepare a set of Lindependent creation and annihilation operators. We reduce the 2L operators {c † j , c j } (j = 1, · · · , L) by half to avoid redundancy 36 . For this purpose, we use the reflection symmetry of f SSD about the center of the system, K 2;ij = K 2,iL−j+1 . The operators that fulfill M −1 a j M = −a j and M −1 b j M = b j about the parity operator M of the mirror reflection are , a 1 , · · · a L 2 ), the following expressions are obtained; where A ± and B ± are the L/4 × L/4 matrices. Next, we find a matrix Q to transform η forĤ 2;i /2 as Here, γ l with l = 1 ∼ L 2 are nonnegative and the other half with L 2 + 1 ∼ L are nonpositive. Since eĤ 2;i/2 = l e γ l η † l η l , one can put this back to the original representation as We thus obtain the Fock operator where from Eq.(40), (e H1;i ) lm = δ lm e −2K * 1;ij for l ≤ L 2 and δ lm e 2K * 1;ij for L 2 + 1 ≤ l. As a third step, we diagonalize Eq.(45) as P −1 T P = diag(ζ l ), by a Bogoliubov transformation ξ † = (ξ † 1 , · · · , ξ † L 2 ,ξ 1 , · · · ,ξ L 2 ) = ΦP . Here, the distribution of eigenvalues is such that half of ln ζ l are nonpositive and the other half are nonnegative. Therefore, by setting ζ l in ascending order and by putting it back to the exponential form with 2 l = | ln ζ l |, we find where for m = l + L 2 we have ln ζ m ≥ 0. We also applied a particle-hole transformationξ mξ † The partition function is obtained from Eq.(32) as with Λ 0 = L j=1 (2 sinh 2K 1;ij ) L/2 L l=1 e l L , being a partition function at L → ∞.
The bond energy along the column and the row are given as

IV. NUMERICAL EXAMINATION
In this section, we numerically demonstrate how SSD works on the Ising model by using the formula in the previous section. We also compare these results with those of classical Monte Carlo simulations.

A. 1D systems
In §II, we found that the maximum eigenvalue of the transfer matrix on the bond at position r i connecting site i and (i + 1) serves as a local partition function on that bond, at its effective temperature, k B T eff (r i ) = k B T /f SSD (r i ).
In the thermodynamic limit, the second largest eigenvalue λ − i is neglected and the total partition function becomes a product of λ + i . This means that the system is an ensemble of (L − 1)-different noninteracting bond degrees of freedom, and unlike a uniform system, each is exposed to its own temperature k B T eff that depends on its location. The form of Eq.(51) indicates that one can obtain a set of equilibrium states with different temperatures k B T eff ranging from k B T to ∞, simultaneously, in a single system.
To examine how accurate the above mentioned description would be at finite L, we numerically evaluate the bond energy of the SSD Hamiltonian in Eq.(6) as a function of k B T eff for several choices of L = 4, · · · , 100, which is shown in Fig. 3(a). The bond energy even at L = 4 shows relatively good agreement with the exact bond energy e bulk = − tanh(βJ) of the L = ∞ uniform Ising model. The inset of Fig. 3(a) shows i (λ − i /λ + i ), which is a rapidly decreasing function of both k B T and L. When this quantity is sufficiently small, Z = i λ + i holds, and λ + i given in Eq.(51) serves as a local partition function for the corresponding local effective temperature, which is fulfilled for most of the temperature range k B T O(0.1J). Therefore, one can realize a canonical ensemble of systems with a variety of temperatures in a single system by properly setting k B T . Figure 3(b) shows the deviation of bond energy against e bulk , which we call an SSD error. For L 40, it is less than 10 −5 . Here, a single partition function for a fixed k B T generates L/2-independent data points with different k B T eff . Therefore, various L/2-sets of data are obtained by varying k B T = 0.1n with n = 1, 2, · · · . The inset shows the SSD error of these series for L = 10; the data obtained near the system center has better accuracy than those near the edges. We also found that the accuracy is improved for higher k B T . We may explain this tendency by a slope of the effective temperature. The slope is gentle at the center and becomes steeper in approaching the edge of the system as shown in Fig. 1(a). The location that gives a certain fixed value of the effective temperature becomes closer to the system center if the system temperature k B T is larger. The local thermal equilibrium is better attained if an additional energy flow caused by the slope of the temperature can be neglected.
The results presented above are in good agreement with the tendency observed in the quantum many-body systems under the SSD. In the quantum cases 18 , the physical quantities at T = 0 measured at the system center reproduce the values in the infinitely large system even when the system is as small as L 20. At finite temperature, it works quite well even at L = 4 in the 1D system 34 . In the same context, the measurements of − σ i σ j in the present classical system is accurately performed near the system center by varying k B T even at small L. Generally, it is easy to increase L by one order of magnitude in classical systems, and the measurements over a wide range of system become accurate enough (see Fig.3(b)). This fact will be of great advantage in utilizing the SSD for Monte Carlo simulations.
In Figs. 3(c) and 3(d), we show the results of a uniform Ising model at finite L to compare with Figs. 3(a) and 3(b), respectively. When λ − /λ + shown in the inset becomes small enough, the bond energy approaches e bulk . In contrast to the case of SSD, the bond energy at L = 4 disagrees with e bulk by about 10 −2 . An advantage of the SSD system over the uniform system is particularly significant at around k B T eff ∼ 1. The finite-size correction of the uniform system remains of order 10 −2 even when increasing L up to 100, where we find λ − /λ + ∼ 1 accordingly.

B. 2D systems
We apply the SSD along the row in Case (i) and along the column in Case (ii) (see Fig.1(c)). The fermionic dispersions along the column, E q /L in Case (i) and l in Case (ii), are presented. For the SSD Hamiltonian, these two cases give different energy dispersions, since for Case (i) q is a good quantum number along the column but for Case (ii) it is not. We checked that when J ij is uniform, Cases (i) and (ii), whose Hamiltonian is the same but the formulation differ, give the same results.

Fermionic dispersions
Let us first compare the energy dispersions of fermions as functions of q in the SSD and the uniform systems. The purpose here is to examine whether the partition function Z → Λ 0 at L → ∞ for the SSD in Case (i) given in Eq.(37) is equivalent to the product of local partition function of each column with different i-dependent k B T eff (r i ). The latter is formally obtained by replacing the uniform temperature k B T of Λ 0 with the i-dependent k B T eff (r i ) in Eq.(28) in §.III A 2. This replacement is equivalent to having the relation where (u) q of the uniform system calculated for each k B T eff (r i ) is summed over different columns on the l.h.s. and E q on the r.h.s. is obtained in Eq. (37). If this equation exactly holds, the classical SSD system in 2D is an extended canonical ensemble of a local subsystems each in an equilibrium of different temperature. Since E q cannot be obtained analytically, we show numerically that this relationship holds almost exactly except for the small deviation at around q ∼ 0. Figure 4(a) shows E q /L in Eq. (35) for Case (i) with L = 50 for several choices of k B T . Since E q is the order-L quantity obtained after multiplying the transfer matrices of all columns, we expect that it can be approximately divided into contributions from different columns if these columns can be regarded as independent subsystems, which is the implication of Eq.(52). In Fig.4(b) we plot together a set of L/2 independent dispersions of fermions representing a single transfer matrix of a uniform system (u) q , each obtained for column(bond-i) dependent k B T eff using Eq. (22). A set of effective temperatures k B T eff for k B T = 0.5 adopted in this calculation is shown in Fig.4(c). At i = 1 ∼ 7 where k B T eff ≤ k B T c in the left panel, the dispersion is a descending function of i, and then for i ≥ 8 in the right panel it ascends with i. By averaging all the dispersions over i = 1 ∼ L, the data points in panel (b) (brown symbols in the l.h.s. panel) are obtained, which is the same data as the one plotted in solid line in panel (a) marked by arrows. We found that except for the very vicinity q/π ∼ 0, the average i (u) q /L and E q /L are in almost perfect agreement. The same comparison holds for other k B T .
In this way, the contributions to E q from each column are well approximated by (u) q under locally-defined effective temperature k B T eff . The result indicates that a picture we proved in 1D also holds in 2D, namely, the system can be regarded as an assemblage of small subsystems having a different canonical temperature k B T eff .

Bond energy
Since we found that the 2D SSD system can simultaneously host L/2 different subsystems with different effective temperature, we can use this fact to evaluate the local physical quantities in each subsystem. Figure 5(a) shows bond energy − σ γ σ γ as a function of k B T eff obtained for Case (ii) using Eqs.(49) and (50), where we set k B T = 0.5 and plot the results for L = 10 ∼ 500. Solid line e bulk is an exact solution of L = ∞-uniform Ising model. The data at L = 10 still deviate from e bulk but when L 50 they almost perfectly overlap with e bulk . The SSD error |e bulk + σ γ σ γ | for the same data set is shown in Fig. 5(b). There are two series of data points following different curves for the same value of L. One is the bond energy evaluated along the columns, and the other one is along the rows. The error is suppressed to less than 10 −3 in a wide range of k B T eff when L 100. Setting k B T to k B T c further suppresses the SSD error as shown in Fig. 5(c). This is because the spatial slope of the effective temperature becomes gentle at the critical temperature k B T c where the finite-size effect is very strong. The local subsystem can more easily attain the  thermodynamic equilibrium at the target effective temperature when the differences of k B T eff with its neighbors are smaller. For comparison, we also calculate the bond energy of a uniform system and plot a finite size error in Fig. 5(d). They are calculated using Eqs. (26) and (27). Near T c ∼ 2.2692J, correlation length diverges and the finite size effect becomes large. This fact is consistent with the peak of SSD error near k B T eff .

Eigen states of fermions
We now examine the spatial distribution of wave functions of fermions when Case (ii) SSD is applied. Fig-ure 6(a) shows the energy levels l of fermions obtained using Eq.(46). We plot the data for several choices of k B T . Since the column direction is no longer uniform, the label l is not a wave number but an index in a descending order of l for l = 1 ∼ 25. The latter half, l = 26 ∼ 50, takes the same value with the former half due to the reflection symmetry.
Following Refs. [21] and [18], we first explain how the SSD term reorganizes the eigen states of fermions. By introducing f SSD (j) = 1−g(j), which fulfills g(j) = g(L+ 1 − j), where δ = 2π/L, one could separate the exponent of Eq.(14) into two parts as, The first term of Eq.(54) is a q-component of Eq. (18). In the uniform system, the eigen state is characterized by a wave number q. By introducing the SSD, one-body states of different values of q mix as in Eq.(55). The particular form of f SSD allows this mixing only between neighboring q's which are discretized in a unit of δ = 2π/L. The amplitude of mixing also depends on q. It takes the largest values at q = 0 and π for the first term, and at q = ±π/2 for the second term. These three are the top or bottom and the middle, of the energy band in Fig. 6(a), respectively. Such moderate mixing generates a wave packet as an eigenstate ofĤ i;2 . Strictly speaking, the final eigenstates are those of Eq.(45) and not ofĤ 2;i , but once we already have a wave-packet state localized in real space, it does not change much by the operation of e H1;i , as it has only a diagonal form in the real-space representation. In Fig. 6(b), a weight of the one-body state (ξ † l |0 ) at a site-j, |(Q i ) jl |, for the dispersion of the fermions at k B T = 0.5 in panel (a) is shown for several energy levels l. Those labeled by l = 1, 2 are the ones providing the largest two l , and are almost completely localized at the edge sites. This is because the coefficients in Eq.(55) is the largest and it efficiently mixes the states over the wide range of wavenumbers, so that the wave packets become a delta function. In all energy levels in Fig. 6(b), the wave packet typically spans over at most three to four lattice spacings, and carries a "bulk" energy l , which gives a rough characteristic energy scale of k B T eff . In the uniform and finite systems, the size effect in quantum state occurs because they are plain waves with discretized wave number. But for SSD, the wave-packet state is localized and does not feel the size of the system. Their local physical quantities behave nearly scale free 22 . This wave-packet-type distribution of fermions in real space supports the picture that the system is an assemblage of local subsystems at different effective temperatures.

C. Monte Carlo simulation
We solved the classical Ising model in a finite system size exactly using the fermionic representation. However, even in the classical systems, the cases with exact solutions are limited. A Monte Carlo (MC) simulation usually serves as a good approximate solver. Applying the SSD to an MC simulation raises a question, whether an SSD system converges to a proper equilibrium state even though the effective temperature of each spin depends on the location? We thus carry out the standard singlespin-flip MC simulations to the 2D Ising model deformed in one direction. The effective temperature differs for each column as shown in Fig.1(c). We collect bond energy for each column separately and plot it against the effective temperature. We discarded the first 5000 MC steps(MCS) and measured the bond energy for 50000 MCS after that. We also performed ten independent MC runs and took an average of the data. The initial spin configuration is the ferromagnetic state with σ γ = 1. Figure 7 shows the MC data obtained at different system sizes. We checked that the MC result of L = 150 is consistent with that of the transfer-matrix method shown in Fig. 5. As shown in Fig. 7(a), the temperature dependence of the bond energy agrees well with the exact results for the whole temperature region. We find a small discrepancy only in the vicinity of the critical temperature, which decreases as the system size increases. The trend is clearly observed in a plot of the SSD error in Fig. 7(b). Both a peak value at T eff = T c and the width were found to scale with 1/ √ L. We also plot in this figure a result of L = 15000 measured only for 500 MCS after discarding 5000 steps. The peak shape at T c is the same as the original measurement with 50000 MCS but the base of the peak is shifted upward roughly by one order of magnitude (∼ 10 times), which is the ratio of a square root of two MCS. Therefore, the SSD error in this off-critical region is controlled by the standard MC statistical error, 1/ √ MCS × 1/ √ L. It suggests that the MC approximation would become exact in the limit of an infinite number of steps in this off-critical region of the SSD system. On the other hand, the SSD error in the critical region near T = T c is due to the systematic one that solely depends on L.
We confirm the validity of the present SSD simulation by examining how precisely we can reproduce the specific heat from our data. The specific heat is usually evaluated as a fluctuation of energy based on the two-point correla-tion of the bond energy over the whole system. However, the correlations between different bonds no longer make sense when the SSD is applied. Instead, we evaluate it from the derivatives of − σ γ σ γ against k B T eff . It is easily performed by a difference between the neighboring effective-temperature data. As shown in Fig. 7(c), the specific heat C is also consistent with the exact result for the off-critical temperature region. The SSD error of the specific heat is consistent with that of the bond energy shown in Fig. 7(b).
Remind that a massive number of temperature data are generated in the SSD simulation only by a small numerical effort. Taking full advantage of this we apply the Gaussian kernel regression (GKR) coupled with the Bayesian inference 37 . The GKR is a machine-learningbased statistical data analysis and with this we can estimate the critical temperature, and can also obtain bond energy as a continuous function of temperature without assuming any analytical function form. Since a larger number of data sets gives better performance of GKR, our SSD MC provides a suitable playground for it. In our case, the specific heat is obtained continuously without taking the numerical derivatives. We only need to take an analytic derivative of the Gaussian distribution function used there. The regression was already proved to reproduce the critical temperature of the classical 2D Ising model within the accuracy of 10 −6 , and the stateof-art temperature-dependent critical exponent that converges to β = 1/8 at T = T c 38 .
We randomly choose 500 data in a range of 1.2 ≤ T ≤ 4.0 and apply the GKR by setting the regression variables (x i , y i ) as where T c and E c (the bond energy at the critical temperature) are parameters to be estimated by the Bayesian inference. Here, we exchange x i and y i from the conventional definition because the bond energy exhibits a singular behavior at T = T c . It is much easier for the regression to model a function with a gentle slope than to model a function with a steep slope. We also know that the specific heat is symmetric and the bond energy is antisymmetric by a mirror reflection of the temperature at the critical temperature in the critical region. To take this prior information into account, we introduce a set of mirror data 39 with respect to the critical point as, x i = −x i and y i = −y i . We mix the data below and above the critical temperature only in the critical temperature region, |y i | < ∆T , where ∆T , which is the width of the critical region, is another parameter to be estimated by the Bayesian inference. The GKR results of L = 1500 and L = 15000 are plotted with lines in Figs (2)) for L = 15000, whereas the exact bulk values are (T c , E c ) = (2.2692 · · · , −0.70710 · · · ). The width of the critical region was estimated as ∆T = 0.150(7) for L = 1500 and ∆T = 0.11(1) for L = 15000. The critical temperature and the bond energy deviate from the exact results only by an order of 10 −4 . Our data in the vicinity of the critical temperature include the SSD error by an order of 10 −3 ∼ −10 −2 as shown in Fig. 7(b). The difference between the neighboring effective temperatures at T c is more than 10 −3 even in the system of L = 15000. The Bayesian inference realizes an accuracy almost ten times better than these SSD errors.
Since the SSD approximation is generally good at the center of the system, we can reduce the SSD error in the critical region by setting the simulation temperature to the critical temperature (in the same context as shown in Figs. 5(b) and 5(c)). Figure 7(d) shows the result of the specific heat. The error bars near T c are much larger than the result of k B T = 0.5. The present number of MCS may not be sufficient because a real simulation temperature is a critical temperature and the critical slowing down may occur. In this model, the specific heat diverges at T c . The exact solution shows that the specific heat reaches C ∼ 2 when the temperature approaches |T − T c | = 4 × 10 −4 . This exact value is reproduced by the Monte Carlo simulation data for L = 15000 within the error bar.

D. Other deformation functions
So far we have studied the effect of SSD on classical Ising models. However, unlike for the quantum models, our results may suggest that the SSD is not the special deformation for classical models. To clarify this point, we performed the same calculation as Case (ii) for other envelope functions. We adopt three functions, whose spatial dependence and the corresponding k B T eff as functions of r/R are shown in Fig. 8(a). One finds that the slope of f sin (r) and f pw (r) are the decreasing and increasing functions of r/R while the slope of f linear (r) is a constant. Whereas, the slope of SSD is first a decreasing function and then becomes an increasing function. Figure 8(b) shows two panels with different system temperatures k B T = 0.5 and 2.2692(T c ) calculating the deformation error (SSD error) of the bond energy along the rows as a function of k B T eff . For k B T = 0.5, f pw gives smaller deformation error than the other three, but for k B T = 2.2692 the same f pw gives the largest SSD error.
The ones for f sin and f linear also have different tendencies depending on k B T . This is because, the deformation error overall tends to increase as the slope d(k B T eff )/dr becomes larger as we discussed previously. The temperature slope at each k B T eff varies depending on the system temperature and the choice of envelope functions. However, since the slope of SSD varies with the moderate tendency compared to the other three cases, and having the optimal zero-slope d(k B T eff )/dr = 0 at both r/R = 0 and 1, it sustains as a moderately stable function, not depending much on the system parameters. Therefore, although one may choose other functions at their purposes, the SSD may be regarded as an optimal function in the sense that it does not require tuning of parameters.

V. SUMMARY AND DISCUSSION
We analyzed the 1D and 2D ferromagnetic Ising model with spatially deformed interactions in the sine-square functional form. We found that this interaction-deformed system is equivalent to the uniform-interaction system with the spatially deformed temperature. To be more precise, we showed by the analytical and numerical analyses that this deformed classical system is an assemblage of small subsystems. Each subsystem locally realizes the equilibrium of a uniform system at its own effective temperature. We propose that this classical SSD at finite temperature gives the approximate extended canonical ensemble with its "state indices" spanning over real space. In the analytical calculation, we first extended the formulation of the conventional transfer matrix method to those of the nonuniform system. We showed that the partition function is exactly obtained even though the interaction is deformed in one direction; this fact is rather trivial in 1D, as the eigenstate of the transfer matrix defined on each bond does not depend on the strength of the interactions. In 2D, the transfer matrices are defined in a unit of a column of the lattice of length L, describing the contributions from 2 L different configurations of the Ising variables. Referring to the previously established approach, the Ising variables are transformed to the noninteracting 1D fermionic operators, and 2 L ensemble average of the Ising variables are mapped to the summation of 2 L different many-body states constructed from the noninteracting one-body states of Bogoliubov fermions. For demonstration, these formulas are numerically evaluated in 1D and 2D Ising models with SSD for system size L = 10 ∼ 500. Notice that it is practically possible to extend it to L 10000 if needed, as it is a one-body problem of fermions on a chain of length L. In the uniform 2D system, the above mentioned Bogoliubov quasi-particles are itinerant plane waves characterized by wavenumbers, but once the SSD turns on, they mix via scattering of the SSD potential and form a set of spatially localized wave packet states. At the same time, according to our picture, "a system with spatially nonuniform interaction bonds J i at the temperature k B T " could be interpreted as "a spatially uniform system with interaction bonds J exposed to the spatially varying effective temperature k B T eff /J = k B T /J i ". Then, the quasi-particle localized on a certain bond feels the effective temperature, and carries the energy corresponding to that of the bulk system at k B T eff . The trace of the product of the exponentials of these quasi-particle energy gives the partition function. The constituent of this product in a unit of the localized wave packet gives the local partition function, carrying the energy typical to that location. The system thus becomes a canonical ensemble of wave-packet states representing the thermal equilibrium at temperature k B T eff . We showed that this picture is valid by evaluating the bond energy numerically exactly for a finite size L. Its deviation from the exact value in the thermodynamic limit is suppressed to less than 10 −3 .
A practical advantage of using the SSD is that it generates a massive number (L/2) of data points with different k B T eff by a single calculation at fixed k B T . Therefore, we can perform Monte Carlo simulations on a large system, and obtain a wide profile of the energy and the specific heat within a sufficient accuracy with a low numerical cost. The SSD error is found to scale roughly with 1/ √ L. When combined with the Gaussian kernel regression, the accuracy improves beyond the SSD error. The applications of the SSD to the MC simulations are very promising.