Comment on"Kosterlitz-Thouless-type caging-uncaging transition in a quasi-one-dimensional hard disk system"[Phys. Rev. Research 2, 033351 (2020)]

Huerta et al. [Phys. Rev. Research 2, 033351 (2020)] report a power-law decay of positional order in numerical simulations of hard disks confined within hard parallel walls, which they interpret as a Kosterlitz-Thouless-type caging-uncaging transition. The proposed existence of such a transition in a quasi-one-dimensional (q1D) system, however, contradicts long-held physical expectations. To clarify if the proposed ordering persists in the thermodynamic limit, we introduce an exact transfer matrix approach to expeditiously generate equilibrium configurations for systems of arbitrary size. The power-law decay of positional order is found to extend only over finite distances. We conclude that the numerical simulation results reported are associated with a crossover, and not a proper thermodynamic phase transition.


I. INTRODUCTION
The work of Huerta et al. identifies a Kosterlitz-Thouless (KT)-type caging-uncaging transition in a system of hard disks confined between parallel walls [1]. That identification is based on the power-law decay of positional order in numerical simulations near close packing, and is supported by the detection of a narrow subpeak in the pair distribution function and of transverse excitation modes in the caging and uncaging regimes. The proposal is intriguing because the presence of a phase transition in such a system is physically unexpected. One-dimensional (1D) and quasi-one-dimensional (q1D) systems with short-ranged interactions have indeed long been considered incapable of exhibiting genuine phase transitions [2][3][4]. (Singularities with respect to changes in structural quantities nonetheless remain possible [5].) KT transitions, however, differ from conventional phase transitions in many respects [6,7]. They leave no thermal feature in the partition function and its derivatives, such as the specific heat [8], and are thus "infinite-order" in nature. The critical phase below the KT transition temperature also differs from conventional ordered phases in that it only exhibits quasi-long-range order. That critical phase is thus only identified from the power-law decay of spatial correlations, with a critical exponent value that changes with system conditions [7]. Do these features exempt KT transitions from traditional expectations for q1D systems? If not, how can one explain the power-law decay of the positional order observed in the simulations of Ref. [1]?
These questions motivate our consideration of an exact transfer matrix scheme that provides equilibrium observables and correlation lengths in the thermodynamic limit, and thus sidesteps hurdles associated with thermalization and finite-size corrections in numerical simulations. Although the memory and computational complexities of * yi.hu@duke.edu transfer matrix treatments increase exponentially with the number of possible pair interactions [9], the approach has already been demonstrated for q1D hard disks with up to next-nearest-neighbor interactions [10,11]. The model studied by Huerta et al. hence lies comfortably within the computationally accessible regime. However, because KT transitions leave no thermal signature and because the correlation length given by the transfer matrix is not directly related the decay of the longitudinal pair distribution function g(x), standard transfer matrix schemes do not suffice. We thus introduce an approach for planting equilibrium configurations at minimal computational cost. By broadening the range of g(x) compared to what Ref. [1] reports, we find that its power-law decay is truncated at large distances, and hence that the KT-like scaling observed in numerical simulations results from a smooth crossover rather than a genuine thermodynamic phase transition.

II. PLANTING METHOD
Because the planting scheme proposed is generic and could be applied to any system solvable by transfer matrices, we first describe it in general terms, and then apply it to the specific q1D system of interest.
A. Transfer matrix setup

Consider a system described by the Hamiltonian
for interactions between a unit a i and its subsequent unit, a i = a i+1 . A unit could be, for instance, m subsequent spins (in a lattice model) or particles (in an off-lattice model) with at most m-th nearest neighbor interactions.
The key step consists of writing the partition function as the trace of a product of transfer matrices, such as Z = tr K N for N identical units. For lattice models at inverse temperature β, matrix entries are with a and a indexing rows and columns, respectively, for the n possible states taken by a. For continuumspace (off-lattice) models, a discretization of space in n segments of size δa similarly gives In both cases, the resulting n × n transfer matrix can be used to write is a matrix of eigenvectors. In the thermodynamic, N → ∞, limit, the leading eigenvalue λ N 0 asymptotically dominates the partition function, and hence the free energy per site is f = − log λ 0 /β. The i-th subleading eigenvalue can then also be used to obtain the i-th correlation length, ξ i = 1/ log(λ 0 /|λ i |).
For 1D and q1D systems with finite-range interactions between units, the transfer matrix is finite with nonnegative entries. According to the Perron-Frobenius theorem, the leading eigenvalue λ 0 is non-degenerate, i.e., λ 0 > |λ 1 |, and the entries of leading left and right leading eigenvectors, u 0 and u −1 0 , are real and non-negative. It then follows that a system described by transfer matrices always presents a finite correlation length for finite β and thus cannot undergo a finite-temperature phase transition. This reasoning should also apply to finitepressure KT-type transitions in q1D hard disks, as we consider below.

B. Generating equilibrium configuration
Unlike molecular simulations, which provide configurations in real space, transfer matrices are probabilistic objects. Structural observables commonly accessible in the former may thus not as easily be obtained from the latter. For instance, although the pair correlation, g(x), can be computed by (inverse) Fourier transforming the structure factor computed as in Ref. [12], this approach is not straightforwardly generalizable to many other structural observables. To sidestep this difficulty equilibrium states can be planted for subsequent analysis from the eigenvectors of these matrices. The marginal probability P (a) (or P (a)δa in off-lattice models) of a state a in equilibrium configurations is indeed given by and the conditional probability that, given a state of a, the subsequent state is a , P (a |a) (or P (a |a)δa in off-lattice models) is Once the leading eigenvector of the transfer matrix is known, we can propagate an equilibrium configuration of arbitrary size by the following algorithm.
1. Generate a "starting" configuration according to Eq. (5) by initializing a cumulative probability distribution array Q(a) based on the marginal probability P (a) such that for indexed states a = 1, ..., n, and Q(0) = 0.
2. Generate a uniformly distributed random variable γ ∈ [0, 1) and choose a state with index a such that 3. Initialize another cumulative probability distribution array Q(a |a) for the conditional probability P (a |a), such that for indexed states a = 1, ..., n, and Q(0|a) = 0.
4. Generate a random variable γ again and choose the subsequent state a from Q(a − 1|a) ≤ γQ(n|a) < Q(a |a).
5. Propagate subsequent states by setting a ← a and repeating steps 3 and 4.
Because the state index a denotes a continuous variable in off-lattice models, the sampling must then be interpolated. Specifically, a and a are found by a = f (γQ(n)) and a = f (γQ(n|a)|a), where f (·) and f (·|a) are interpolation functions that map Q(a) → a and Q(a |a) → a , respectively. For sufficiently large n, the choice of interpolation method does not significantly affect the result, and for the q1D hard disk model considered a simple cubic interpolation scheme suffices.

C. q1D hard disk scenario
We now specialize to the case of hard disks of radius of d in a channel infinite in the x (longitudinal) direction, and defined by hard walls a distance D apart in the y (transverse) direction. The scaled transverse length free to disk centers is thus h = (D − d)/d. As in Ref. [1], we further specialize to the case D/d = 3/2. For 1 < D/d < 1 + √ 3/2 = 1.866 . . ., at most nearestneighbor interactions between disks are geometrically possible, and thus a i in Eq. (3) is naturally taken as the vertical coordinates of disk centers, y i . For such a q1D system, the partition function for the isothermal-isobaric (constant N P T ) ensemble can be written using a transfer matrix approach [13][14][15]. (Although this scheme differs from the canonical, constant N V T ensemble simulated in Ref. [1], in the thermodynamic limit the two are rigorously equivalent.) An entry of the transfer matrix then reads where σ(y, y ) = d 2 − (y − y ) 2 is the contact distance between two neighboring disks and F is the force associated with the longitudinal pressure. Note that we have here replaced δa in Eq. (3) with √ δyδy . The associated change of variable does not affect the eigenvalues of the transfer matrix, but keeps K symmetric even for uneven discretization of y, and hence u = u −1 .
At high pressures, i.e., near close packing, disks are mostly confined near the walls. To improve the numerical accuracy of the discretization scheme, Ref. [10] suggested a change of variable y → t, such that and c is adjustable. In particular, increasing c refines the grid around the walls. The discretization of the new variable t ∈ [−1, 1] is then uniformly spaced by δt and the transfer matrix becomes K tt = exp −βF σ(y(t), y(t )) dy dt t · dy dt t · δt , (10) where dy/dt = a + b c sech 2 (ct). Results for this regime are fairly insensitive to different choices of c ∼ O(1) and n ≥ 100. Without loss of generality, we thus set c = 3 and n = 1000. Note that eigenvalue evaluation is then nearly instantaneous on a standard desktop computer. The leading and subleading correlation lengths given by the transfer matrix for this systems correspond to the correlation length of the zig-zag order, ξ y = ξ 1 , and of the longitudinal spacing, ξ δx = ξ 2 [10]. These lengths hence describe the decay at large distances, |i − j| → ∞, of respectively. Note that the oscillatory nature of Eq. (11) for zig-zag order is associated with λ 2 being negative [10]. From the leading eigenvalues and eigenvectors of K we can also generate y-coordinates of disks according to the scheme described in Sec. II B. Longitudinal spacings between neighboring obstacles, δx i , are then generated according to the distribution rule for the constant N P T ensemble For each pressure considered, we generate configurations of longitudinal size L ≥ 500, and compute the pair distribution function by averaging over 400 independent realizations. Because the transfer matrix is evaluated for the constant N P T ensemble, we first determine the longitudinal force F under the longitudinal number density ρ = lim N,L→∞ N/(L/d) (Fig. 1), and then compute g(x) from planted configurations (Fig. 2). At short distances, the transfer matrix and the simulation results are fully consistent, including the apparition of a small sub-peak at high densities (see [1, Fig. 2(a)]). At long distances, however, while Ref. [1] reports that g(x)−1 decays with a characteristic power-law beyond a certain density, we find that for sufficiently large x, g(x) − 1 decays exponentially for all densities considered. For ρ = 1.1111, in particular, we observe that the power-law-like decay terminates at x ∼ 100, even though it persists at least up to x ∼ 200 in simulations. Because the leading correlation length from the transfer matrix, ξ y = 364 (Fig. 3), is then very close to the simulated system size, N = 400, this discrepancy is most likely a finite-size correction. The resulting self-interactions in x through the periodic boundary condition obfuscates the exponential decay of g(x) − 1. As additional evidence, we note that g(x) computed from the inverse discrete Fourier transform (DFT) of the structure factor S(q) (obtained as in Ref. [12]) can be made to look like the simulation results by choosing a discretization spacing δq that corresponds to a finite system size L = 2π/δq (Fig. 2(c)). Is the intermediate-range algebraic decay then an echo of a two-dimensional transition? The consideration of a purely 1D model suggests not. Recall that for 1D rigid rods of length d [16,17], wherex = x/d andp = βF d are the reduced distance and pressure, respectively. Note that at highp and smallx, the height of the k-th peak is approximately the maximal value of the k-th summand, g(x k max ) = (1 +p) Note also that at high pressures disks in the q1D system studied can be viewed as 1D hard rods of effective Because the variance of y is non-zero, however, the peak height of g(x) is also reduced by a constant multiplicative factor, c < 1. By replacingp = βF d and k = xmax/d +1/p 1+1/p , we thus obtain an approximate form for the algebraic x −1/2 max decay of the peak height, Setting c = 0.7 nicely fits both ρ = 1.1111 and 1.1400 for g(x max ) > O(1) and x max > O(1) (Fig. 2(b)). For g(x max ) − 1 O(1), however, the single summand approximation to g(x) breaks down. The algebraic decay thus only persist up to √ k ∼p, i.e., x/d ∼ (βF d) 2 . Interestingly, the exponential decay of g(x) − 1 as x → ∞, which truncates the algebraic decay, is also a 1D feature [18]. This decay is controlled by the characteristic length ξ g(x) of a Tonk gas [11,Eq. 24]. As h → 0, the rows of the transfer matrix converge, and hence the associated correlation lengths vanish, but ξ g(x) persists. In other words, as long as h is small, the algebraic decay reported in Ref. [1] is essentially a 1D feature, and thus not an echo of KT-type transition. Beyond the nearestneighbor interaction regime, h > √ 3/2, however, the decay of g(x) does genuinely become quite rich [19].
Finally, Ref. [1] investigates the distribution of longitudinal spacing δx, over which window-like defects become substantial, and identified the onset of caging-uncaging transformation around ρ ost = 1.0526. Interestingly, this signature can also be found in ξ y and ξ δx . These correlation lengths, which have been previously studied for q1D systems [10,11,14]-including for D/d = 3/2-are shown in Fig. 3 for reference. The reported ρ ost coincides with the onset of exponential growth of ξ y with pressure (longitudinal force), as analyzed in Ref. [14]. It further coincides with the maximum of ξ δx reported in Ref. [11]. We thus conclude that the phenomenon reported in Ref. [1] is also associated with that crossover.

IV. SUMMARY
Based on direct evidence gleaned from an exact planting scheme derived from transfer matrices, we have demonstrated that the power-law-like decay of positional order reported Ref. [1] only persists over finite distances. The pre-asymptotic power-law decay of the pair correlation is rooted in 1D physics, and not in any KT-type transition or its echo. The suggested uncaging transformation in Ref. [1] instead coincides with anomalies identified in earlier studies, and may thus be considered as a part of that same crossover.