Parafermion excitations in superfluid of quasi-molecular chains

We study a quantum phase transition in a system of dipoles confined in a stack of $N$ identical one-dimensional lattices (tubes) polarized perpendicularly to the lattices. In this arrangement the intra-lattice interaction is purely repulsive preventing the system collapse and the inter-lattice one is attractive. The dipoles may represent polar molecules or indirect excitons. The transition separates two phases; in one of them superfluidity (understood as algebraic decay of the corresponding correlation functions) takes place in each individual lattice, in the other (chain superfluid) the order parameter is the product of bosonic operators from all lattices. We argue that in the presence of finite inter-lattice tunneling the transition belongs to the universality class of the $q=N$ two-dimensional classical Potts model. For $N=2,3,4$ the corresponding low energy field theory is the model of Z$_N$ parafermions perturbed by the thermal operator. Results of Monte Carlo simulations are consistent with these predictions. The detection scheme for the chain superfluid of indirect excitons is outlined.


I. INTRODUCTION
Emergence of Majorana fermions (see in [1]) in topological insulators [2] has inspired a search for such fermions in other condensed matter systems. In this article we show that parafermions [3], of which Majorana fermions represent a particular case, describe excitation spectra of quantum chains (strings) of polarized dipoles. Material realization of such systems has became possible due to the recent breakthroughs in creating and trapping high density samples of (polar) molecules [4]. As proposed in Ref. [5], multi-layered structures of indirect excitons [6] may also form similar systems in the form of excitonic chains. Each indirect exciton (not to be confused with the excitons formed at non-Γ point) has static dipole moment due to a spatial separation of electron and hole. Interaction between the dipoles in the N -layered structure can encourage a formation of the excitonic chains similar to chains of polar molecules. Since light field E and excitons are coupled linearly a state of excitonic field ψ is imprinted directly on the emitted light. As a consequence, properties of excitonic chains can be explored through light emission providing a new powerful experimental tool to study strongly correlated systems.
So far, quantum chains have been studied in various analytical approximations which neglect tunneling of particles along the chains. In Ref. [7] it has been proposed that stiff dipolar non-interacting quantum chains may form Bose-Einstein condensate. Inter-layer pairing in bilayered 2D dipolar fermionic systems has been studied in the BCS approximation in Refs. [8]. The dimerization transition in the 2D multilayered geometry of dipolar fermions was analyzed in Ref. [9], and it has also been proposed that for strong dipolar interactions long chains can form by the N-clock phase transition [9]. Fermionic dipolar molecules forming a mixture of single fermions, dimers and fermionic trimers in 1D N = 3-layered system has been discussed in the ideal gas approximation in Ref. [10].
In low dimensions quantum fluctuations are enhanced and therefore quantum particles from different 1d tubes may develop strong correlations. This can lead to interesting physics. The difficulty is that such system, in general, is not amenable to the standard mean field or perturbation expansion methods and one has to resort to a combination of non-perturbative techniques and numerics. A numerical study of quantum chains which takes into account the partial-chain exchanges as well as the intra-chain dynamics has been performed in Ref. [11] for the case of zero inter-layer tunneling. It has been shown that polar molecules in the N -layered geometry can form flexible (quantum rough) chains, and these chains can undergo a quantum phase transition to a superfluid phase characterized by off-diagonal long-range order (ODLRO) in the N -body density matrix, while all M -body density matrices with M < N show insulating behavior regardless of the filling factor (provided it is the same in each layer). If the inter-layer (dipolar) Z X FIG. 1: (Color online) Schematic picture of N = 3 parallel 1d lattices (also called "tubes") stretched along the X-axis. The dipole particles (represented by arrows) occupy lattices sites (full circles). The arrows show the polarization of the particles dipole moments (along the Z-axis). The dipoles may tunnel between nearest sites along the X and the Z directions.
interactions are weak, a stack of N layers features a N -component superfluid (N-SF). Once the interaction becomes stronger, the non-dissipative drag between layers will eventually convert the N-SF phase to the flexible chain superfluid (CSF) characterized by ODLRO only in the N -body density matrix [11]. The corresponding transition is continuous in the 1D-geometry and discontinuous in the 2D-geometry for N > 2 [11].
In the present work we study the transition between N-SF and CSF states in the presence of finite inter-tube tunneling. Since for 2D-layers (N > 2) the transition is of the first order and the inter-layer tunneling cannot change this, we concentrate on the 1D-geometry depicted in Fig. 1, that is, when layers may be considered as tubes.
Our main findings are the following. A quantum phase transition into the N -chain superfluid (in 1+1 dimensions) is in the universality of the classical q = N 2D Potts model. That is, for N = 2, 3, 4 the transition is a continuous one and for N > 4 it is of the first order. For N = 2 we develop a microscopic low energy description in terms of the field theory of two species of Majorana fermions and one gapless bosonic field. For N = 3, 4 instead of a detailed derivation we present arguments based on symmetry of the problem and on results of our numerical calculations.
The paper is organized as follows. In Sec. II we discuss possible phases. Then, in Sec. III we will describe the microscopic lattice model accounting for a system of N coupled tubes. Then we formulate field theoretical description of the system valid in the continuum limit. This description is rigorously derived for the case of two tubes (N = 2) where the low energy limit leads to a model of relativistic Majorana fermions. The correlation functions characterizing light emission from the N = 2 excitonic system are derived. We also present plausible arguments concerning possible field theory for the cases N = 3, 4. These arguments are supported by the Monte Carlo calculations presented in Sec. IV. The Monte Carlo procedure is performed for the coarse-grained dual version of the Hamiltonian in the discretized time approximation. Finally, in the Conclusion we will give a summary of the main results and perspectives for detecting excitonic CSF by N -photon correlation spectroscopy.

II. ORDER PARAMETERS AND UNIVERSALITY OF THE TRANSITION TO THE CSF PHASE
Without inter-tube tunneling in the N-SF state [11] each tube is characterized by its own phase. In the case of finite inter-tube tunneling (along Z in Fig. 1) these phases lock in into a single phase field ϕ so that the superfluid (SF) is characterized by the bosonic operator ψ ∼ e iϕ . As is well known, in 1D the real long range order is substituted by quasi long range order characterized by nonzero stiffness and algebraic decay of certain correlation functions at zero temperature. The SF phase of our system is characterized by an algebraic decay of the bosonic field. Meanwhile in the Chain Superfluid Phase (CSF) correlators of individual Bose operators ψ † (x, z)ψ(x ′ , z ′ ) decay exponentially with respect to |x − x ′ | and an algebraic order pertains to the product of operators of all tubes describing the order of quasi-molecular complexes each consisting of N bosons. It is important that Ψ N , Eq.(1), is invariant with respect to the transformation ψ(x, z) → exp(2πim(z, x)/N )ψ(x, z) where m(z, x) is defined modulo N and obeys the constraint z m(z, x) = pN, p = 0, 1, 2, .... Thus, m(z, x) can be broken as m = m ′ +m into the discrete global part m = p, p = 1, 2, ..., N − 1 and the local gauge-type m ′ (z, x) obeying z m ′ = 0.
Setting aside the discussion of a possible role of the local-gauge symmetry, we note that the global transformation forms a discrete symmetry group which determines the universality of the SF-CSF transition. Among the possible candidates one can consider the p-clock model and the (standard) Potts model (see in Ref. [12]). While the case N = 2 should be assigned to the Ising universality class [13], the nature of the transition for N > 2 is not obvious at all. Naïvely, one may anticipate the p-clock universality because of the nearest neighbor tunneling (between tubes). In what follows we will show that such expectation is not correct, and the criticality is controlled by the standard Potts model (also called as Ashkin-Teller-Potts model). Accordingly, for 1D tubes (that is, D = 1 + 1 membranes) it should be continuous for N = 2, 3, 4 and discontinuous for N > 4.

III. MICROSCOPIC HAMILTONAN AND THE EFFECTIVE MODEL FOR N = 2 IN TERMS OF MAJORANA FERMIONS
Each tube represents an optical lattice occupied by particles with Bose statistics (polar molecules or indirect excitons). The microscopic Hamiltonian H describing SF and CSF has the following form: Here a † xz , a xz are creation (annihilation) operators creating (destroying) a boson at site x belonging to zth tube; n xz = a † xz a xz denotes onsite density operator obeying the hard-core constraint; V xz;x ′ z ′ describes the matrix element for dipole-dipole interaction between sites (xz) and (x ′ z ′ ). It is characterized by the strength V d = d 2 z /b 3 z , where d z stands for the induced dipole moment and b z denotes distance between two nearest tubes. This interaction is mainly attractive along the z-direction and repulsive along the x-direction. In this paper we will be studying a simplified version of the model (2). Specifically, we will reduce the dipole-dipole interaction to just a nearest-neighbor attraction V 1 along the Z-direction and nearest-neighbor repulsion V 0 along the X-direction. Clearly, such approximation cannot change neither the low energy physics nor universality class of the transition.
A. Two tubes (N = 2). Low energy decsription In the low energy limit the microscopic Hamiltonian (2) can be replaced by the effective model describing the SF to CSF transition (in the chosen approximation). Taking into account the single occupancy constraint we can rewrite Hamiltonian (2) in terms of the Pauli matrix operators. Restricting ourselves to the simplest case of two tubes, we have: where σ ± operators stand, respectively, for the bosonic creation and annihilation ones a † , a in Eq.(2) and n = σ z + 1. We assume that here as everywhere throughout the paper the density fluctuations (the total 4 one and the difference between the tubes) are incommensurate with the lattice. In the context of model (3) it is achieved by a proper choice of the chemical potential µ.
Following Schulz [14] we will treat this model at low energies using bosonization technique (see also Ref. [15]). The model describing each tube is the spin S=1/2 XXZ model; in the continuous limit it is equivalent to the Gaussian model: where a = 1, 2 labels the tubes and Θ a is the field dual to Φ a : Luttinger parameter K and velocity v are determined by the intra-chain interactions and the chemical potential. For model (3) at µ = 0 we have [16]: and we assume that t > 0. Then the continuum limit of the spin operators is given by the following bosonization formulae: where dots stand for less relevant operators and C, C z are amplitudes determined by the short range physics and a 0 is the short range cut-off. The Fermi momentum k F for each chain is determined by its chemical potential µ with the convention that k F (µ = 0) = 0. As we have mentioned above, we will always assume that µ = 0 so that the spin fluctuations are incommensurate with the lattice. Substituting (6) into (3) and defining the fields we obtain the Hamiltonian H = H + + H − , where H + describes the symmetric mode (+): and H − contains only the anti-symmetric fields: At this juncture we note that the Hamiltonian (8) describes a mode which is always critical. The corresponding order parameter is Ψ + = σ + 1 σ + 2 ; according to (6) 5 Although in one dimension such order parameters with continuous symmetry do not have vacuum averages, at T = 0 their correlation functions exhibit slow (algebraic) decay. In the present case However, there may be operators with correlators decaying faster than that of Ψ + . They are different in different phases of our model. Depending on which of the cosines in (9) takes over, the ground state of this model describes either quasi long range superfluid order in each tube or pair density wave state. The latter state has a singularity in the density-density correlation function at the finite wave vector 2k F . When both cosines are relevant (that is at 1/2 < K − < 2) these states are separated by a quantum critical point (QCP), the location of which is approximately determined by the relation where Λ is the ultraviolet cut-off determined by the lattice. The vicinity of the QCP can be studied analytically when K − ≈ 1 (this will be our assumption throughout the rest of the paper). The result for K − = 1 is that the transition belongs to the Ising model universality class. By continuity this continue to hold throughout the entire region of existence of QCP 1/2 < K − < 2.
For K − ≈ 1 it is convenient to refermionize (9) with the result Where m ± =t ⊥ ±Ṽ 1 and ρ L,R and η L,R are left-and right-moving components of Majorana (real) fermions: We comment that the fermionization of the system in terms of real (Majorana) fermions is consistent with the Ising type symmetry breaking. This model is equivalent to the continuum limit of two quantum Ising (QI) models coupled by the energy density operators [17]. The transition occurs when one of the Majorana masses becomes zero. To access the correlation functions we need to express the original spin operators in terms of the Ising model fields: where Q = π/a 0 + 2k F , s ± (µ ± ) are Ising order (disorder) parameters for the Ising models represented by ρ and η fermions respectively. These operators are nonlocal in terms of fermioms. Their explicit expressions are not needed here, all we need to know is that in the part of the phase diagram m > 0 we have σ = 0, µ = 0 and for m < 0 we have σ = 0, µ = 0. Therefore att ⊥ >Ṽ 1 > 0 when both masses have the same sign both σ ± have vacuum averages. Replacing these operators in (16) by their vacuum averages we get Since Φ + has a gapless spectrum, the corresponding correlation function decays algebraically with the exponent K + /2 which is four times smaller than the exponent for Ψ + (11). In the other phase m + m − < 0 and a similar replacement can be done for the oscillatory part of the density operator yielding The latter situation corresponds to Pair Density Wave (PDW). The density oscillations in this phase exist alongside the superfluidity of coupled pairs, though in the presence of disorder PDW is pinned [18] and the superfluidity is not destroyed. The QCP separating the two phases occurs when one of the masses (if

B. Correlation Functions
The SF-CSF transition can be determined from measurements of various correlation functions. Here we will specifically address the situation when the dipoles are indirect excitons confined in bi-layered structures. Then, since such excitons can be converted directly into light, the excitation spectrum can be extracted from measurements of light emission. Such emission is described by the linear coupling of the electric field E to the excitonic operatorsâ z,x ,â + z,x : Therefore in the first order of perturbation theory in the excitonic transition matrix element d the emission or absorbtion probability of a single photon is related to the imaginary part of the < σ + σ − > correlation function. According to (16) we have where G σ and G µ are two-point correlation functions of s and µ operators. These correlators of the Ising model are well known. In the ordered state of the Ising model the Lehmann expansion for G s contains matrix elements between the vacuum and the states with even number of Majorana fermions (including zero) and for G µ it contains matrix elements with the odd number of fermions [19]. In the disordered state s and µ are interchanged. Keeping this in mind in the SF phase we obtain the following expansion for (21): where θ 12 = θ 1 − θ 2 , γ = K + /4 and the dots stand for the matrix elements between the vacuum and states with energies higher than 2m − .
For simplicity in what follows we will set v + = v − . Performing the Fourier transformation of (22) and the analytic continuation from Matsubara frequency to real one iω → ω + i0 and then taking the imaginary part, we obtain where with s 2 = ω 2 − (vk) 2 . The thin green line features the dependence (23) in Fig. 2. In the vicinity of the 2-particle threshold s 2 = 4m − 2 (24) behaves as indicating that the 2-particle continuum is very weak. Thus in the superfluid phase the absoption is dominated by the first term in (23) originating from the gapless excitations of the condensate.
To observe clear signs of the Majorana mode one has to make measurements in Pair Density Wave phase where the gapless excitations can be emitted only together with one Majorana fermion. The dominant contribution comes from G s+ G s− term in (21); the function G s+ is replaced by constant as before, but G s− now contains the emission of one Majorana fermion. As a result we get where the dots stand for the terms containing emissions of at least three massive particles. Therefore in the interval (3m 2 − > ω 2 − (vk) 2 > m 2 − ) we have: As we see, the spectral function here has a strong singularity at the one-particle threshold if 2γ − 1 < 0 (that is, K + < 2). Such feature is shown by the thick red line in Fig. 2.

C. N tubes symmetrically coupled
Now we consider a system of N > 2 tubes coupled to each other in such a way that each tube interacts with all others. The treatment in this case is more complicated since we have to resort to 8 non-Abelian bosonization (see, for instance, [15], [20]). As a starting point we take non-interacting tubes with SU(2) symmetry (V 0 = t || ). Then an invidual tube is equivalent to spin S=1/2 isotropic Heisenberg antiferromagnet and in the continuum limit is described by the SU 1 (2) Wess-Zumino-Novikov-Witten (WZNW) model. The sum of N SU 1 (N) WZNW Hamiltonians can be decomposed as (see, for instance [21]) This decomposition should be understood in the sense that operators (primary fields) of the critical theory on the left hand side of the identity can be written as products of operators belonging to the critical field theories on the right hand side. Decompositions of that kind can be very helful outside criticality if the perturbation happens to be such that it does not act in at least one of the sectors. Decomposition (28) is consistent with the fact that central charges of the theories on the left-and right hand side of (28) are equal: The U(1) subsector of (28) corresponds to the symmetric bosonic phase (N-tube generalization of Φ + from the previous subsection). Since the inter-tube interaction does not contain this field, it remains gapless. As far as the other sectors are concerned, we will leave a detailed analysis to future publications and only formulate some conjectures. Information extracted from our numerical calculations suggests the following scenario. Close to the critical point the SU 2 (N )/U N −1 (1) -sector (Gepner's parafermions [22]) is weakly coupled to the rest. This coset sector remains massive throughout the entire phase diagram, at least in the part where t ⊥ > 0,

IV. NUMERICAL RESULTS
The QPT transition discussed above is characterized by disappearance of the algebraic off-diagonal order in all M-body density matrices, where M = 1, 2, ..., N − 1. Specifically, as the inter-tube interaction is increasing the order existing in all M-body density matrices must eventually vanish up to the order M = N − 1. At the same time, the order remains, practically, unaffected in the N-body density matrix.

A. M-body density matrix
The M -body density matrix D M can be written explicitly as where ... stands for the quantum-thermal averaging. In 1D SF D 1 (x, z; x ′ , z ′ ) ∼ 1/|x − x ′ | b , b < 1, exhibits algebraic order at large |x − x ′ |. In the CSF, D 1 (x, z; x ′ , z ′ ) ∼ exp(−|x − x ′ |/ξ 0 ), ξ 0 ∼ 1, that is, it becomes short ranged at T = 0 regardless of the filling factor. Thus, in the CSF phase the N -body density matrix is characterized by the exponential decay D N (x 1 , ..., x m ; x ′ 1 , ..., x ′ m ) ∼ exp(−|x m1 − x m2 |/ξ 0 ) with respect to any pair x m1 , x m2 of coordinates from the set x 1 , ..., x m (or x ′ 1 , ..., x ′ m ). In the CSF there is also the algebraic order The transition SF to CSF can be detected by the critical behavior of any density matrix. In particular, the same criticality controls the long-distance behavior of . That is, in SF phase D N is trivially long-ranged with respect to |R cm − x m | because D N can simply be factorized into a product of D 1 . In contrast, in the CSF-phase, while exhibiting ODLRO with respect to R cm − R ′ cm , D N is short-ranged with respect to |R cm − x m | (or |R ′ cm − x ′ m |). Thus, the criticality can also be detected by measuring the behavior of the relative distances x m (or x ′ m ). Specifically, we have considered the square of so called gyration radius [11] as the mean of with respect to the first set of the coordinates of D N defined in Eq.(30), provided the coordinates from the second set are kept within some distance ∼ ξ 0 from R ′ cm [24]. More specifically, x k , where k = 1, 2, .., N , represents the x-coordinate in the k−th tube.
In the SF of a length L, In what follows we will be calculating the mean of the ratio G N = R 2 g /R 2 0 , so that it is changing from G N ≈ 1 in the SF state to G N ∼ 1/L 2 ≈ 0 in the CSF phase. It is worth mentioning that R g can be viewed as a typical width of a chain. For strongly bound case this width is ∼ ξ 0 ≈ 1, and in the SF phase it is ∼ L, and, thus, it exhibits critical behavior typical for correlation length.

B. J-current formulation
Hamiltonian (2) and the gyration radius (31) have been used for ab initio simulations of a single chain (with exactly one polar particle per layer (in d = 2) or tube (in d = 1) for the case t ⊥ = 0 [11]. It was found that the chain can undergo quantum roughening transition with the tuning parameter being the interaction strength V d . The transition is, practically, insensitive to the interaction range.
The Monte Carlo simulations at finite densities n < 1 ( incommensurate with the lattice along the tubes) in each layer have been conducted in the discrete-time J-current-type formulation [25] of the Hamiltonian (2) [11]. For the purpose of analyzing the universality of the transition this approach turns  out to be much more efficient than the ab initio one. Here we will be using the model where the intertube tunneling is allowed. The actual dipole-dipole interaction will be replaced by onsite attraction between neighboring tubes, with the intra-tube dipole-dipole repulsion ignored. This approximation becomes essentially exact when n << 1: while the inter-tube attraction is not affected, the intra-tube dipole-dipole repulsion between atoms scales as ∼ n 3 → 0 in 1d and, thus, becomes irrelevant. The corresponding space-time action, then, becomes where J b is the integer bond current obeying Kirchhoff's conservation law [25]. In some sense, these conserved currents represent world-lines of particles in imaginary discrete space-time, with H J being the action in Feynman's path integral. The summation in (32) is performed over all space-time bonds b (coming out from a space-time site (x, τ, z) either along ±x or along imaginary time ±τ or along ±ẑ directions); ∇ z J b ≡ J b (x, τ, z + 1) − J b (x, τ, z); µ denotes chemical potential. [Here we tuned µ to have 1/2 filling of bosons per site in each tube]. Periodic boundary conditions along space 0 < x < L − 1, 0 ≤ z ≤ N − 1 and along imaginary time 0 ≤ τ ≤ β, with β = L, where L = 2, 3, ...., have been used. The coefficients K, U can be related to t z,z ′ and V xz; The case K(ẑ) = ∞ corresponds to zero inter-tube tunneling (studied in Ref. [11]). Here we will focus on K(ẑ) = K(x) = K(τ ) situation as the one which naturally represents the whole universality class.
It is worth emphasizing that the inter-tube attraction between two particles (in neigboring tubes) located, respectively, at the space-time points (x, z, τ ) and (x, z ± 1, τ ) is described by the terms ∼ −U J(x, z, τ ) J (x, z ± 1, τ ). Accordingly, when these two particles form a bound state, their world-lines stay close to each other to gain the binding energy ∼ U . Similarly, a chain of several particles is represented by a bundle of several world-lines forming a membrane in the space-time.
The action (32) can be viewed as a coarse grained dual representation of the Hamiltonian (2). While being not precise for quantifying finite energy (non-universal) properties of the system, the J-current model [25] belongs to the same universality class as the original model (2). Thus, for the purpose of this work and for sake of numerical practicality, it will be sufficient to study the model (32).
We also note that, while being formally defined on the lattice, the model (2) and its dual formulation (32) account well for the continuous space-time situations at low energies as long as the filling factor n remains incommensurate with the lattice. Long-range intra-tube repulsion may complicate the situation by inducing crystalization at, say, n = 0.5 and, thus, shifting the CSF phase to lower densities. Such feature, however, does not affect the universality of the N-SF to CSF transition, and, in order to establish it in a most efficient way we simply turn off the intra-tube repulsion and study the case n = 0.5.
Monte-Carlo simulations of the model (32) have been performed within the Worm Algorithm approach [26]. Green's function in imaginary time (as well as the density matrix (30)) is given by the statistics D N ({x m , τ m , z m }; {x ′ m , τ ′ m , z ′ m }) of "sources" and "sinks" of the bond currents located, respectively, at (x m , τ m , z m ), m = 1, 2, ..., N , and (x ′ m , τ ′ m , z ′ m ), m = 1, 2, ..., N , lattice points. In order to insure the condition |R ′ cm − x ′ m | ≤ ξ 0 , while (x m , τ m , z m ), m = 1, 2, ..., N , are free to take any value, we , and, accordingly have evaluated the means of the normalized gyration radius G N and of the center of mass distance For sake of numerical efficiency we have symmetrized the model (32) by choosing U (b) independent of the type of a bond, that is, U (b) = U . The CSF phase has been identified by the condition |R cm − R ′ cm | /L = const and G N ∼ o(L −2 ) for U > U c , where U c corresponds to the QCP. In the SF phase (that is, U < U c ), while the first condition remained, practically, unchanged, G N ≈ 1 with high accuracy. The criticality of the SF-CSF transition has been analyzed through evaluating the divergent behavior of d G N /dU in the vicinity of U = U c .
C. Finite size scaling of the gyration radius As discussed above, d G N /dU exhibits singularity in the limit L → ∞. The change from G N ≈ 1 to G N ≈ 0 occurs in a narrow range δU = |U − U c | around the QCP U c ∼ 1. Such a behavior is clearly seen in Fig. 3: the range of the transition δU narrows as L increases. It is important to emphasize that in the following analysis of the criticality the knowledge of the exact value of U c in the thermodynamical limit is not required. All we need is an approximate range where the derivative exhibits a clear sign of divergence.
The range δU is controlled by the diverging correlation length ξ(U ) ∼ |U − U c | −ν , ν > 0. According to the finite size scaling approach, G N can be represented as some regular function F (y), y = L/ξ(U ) varying from F (y = 0) = 1 to F (y = ∞) = 0 over the range y ∼ 1. Thus, d G N /dU ≈ F ′ y/δU ∼ L 1/ν . Loosely speaking, one can view this relation as d G N /dU ≈ 1/δU, δU ≈ L −1/ν → 0. It should be noted that at any finite L the actual divergence contains the so called subleading terms -powers of L smaller than 1/ν. These terms are the main source of systematic errors in our analysis.
We have evaluated d G N /dU numerically by Monte Carlo [26] and constructed the graphs d G N /dU versus G N by scanning over U around the critical point U c for sizes L = β = 10, 20, ...300. These graphs turn out to be self-similar so that d G N /dU for all sizes > 10 collapsed on a single master curve by simple rescaling of d G N /dU for size L 1 to another size L 2 as d G N /dU → λ(L)d G N /dU . Then, the rescaling coefficient λ, which represents the inverse width δU as λ ∝ 1/δU , has been plotted in the log-logL axes in order to determine the critical exponent ν. The results of this procedure are presented on Figs. 4,5 for the case N = 3. The same procedure has been used in the cases N = 2, 4 as well. The found exponents are: ν = 0.972 ± 0.02 for N = 2, ν = 0.835 ± 0.015 for N = 3, ν = 0.735 ± 0.015 for N = 4. The errors include statistical errors as well as the systematic errors due to the subleading contributions. We note that the value of ν for N = 2 is consistent with the d = 2 Ising (or q = 2 Potts) universality. We also note that the values of ν for N = 3 and N = 4 are consistent with the corresponding ones ν = 0.837 and ν = 0.756 obtained by the Renormalization Group calculations for the q = 3, 4 2d Potts model [27].
The above data collapse fails for N > 4 after reaching some size L ∼ 50 − 100 for N = 5 and much smaller sizes for N > 5. Furthermore, the energy histogram develops bimodality typical for I-st order transitions, Fig. 6. While for N = 5 the bimodality develops on sizes L ≥ 400, β = L, for N = 8 it is already well developed at L = 160, β = L. Such features are consistent with the I-st order transition in 2d Potts model for N = q > 4.
The finite size analysis has been applied to the case of zero inter-tube tunneling as well, when the transition is expected to be in the BKT universality. That is, ξ ∼ exp(...|U − U c| −1/2 ). The variation of the gyration radius G XY in this case can also be represented by some regular function F (y) characterized by the range y ∼ 1 with y = L/ξ . Thus, d G XY /dU ≈ F ′ y ∼ (δU ) −3/2 ∼ (ln(L/L o )) 3 (where L o stands for some microscopic scale) at its maximum. The maximum value of this derivative has been plotted as a function of L = 30, ..., 600 in Fig.7. As can be seen the fit of the data is consistent with the ln 3 L dependence with high accuracy.

V. SUMMARY AND DISCUSSION
As shown above, the Majorana fermions description of the Ising-type transition to CSF state in the twochain system is consistent with the numerical evaluation of the correlation length exponent ν = 0.972±0.02 (versus its exact value ν = 1). The quantum transitions in the cases N = 3, 4 are characterized by emerging enlarged symmetries: Z 3 for N = 3 and Z 4 for N = 4. The predicted values ν = 5/6 ≈ 0.833, N = 3, and ν = 3/4 = 0.75, N = 4, are matched well by the corresponding numerical ones ν = 0.835 ± 0.015 for N = 3 and ν = 0.735 ± 0.015 for N = 4. The symmetry enlargement occurs despite the short-range nature of the tunneling between the tubes. In other words, the critical behavior proceeds as though the tunneling between all tubes is the same. Such feature -identical interaction between all elements -is typical for the standard Potts model [12] and should be contrasted with the p-clock model.
The detection of the CSF order as well as the criticality to SF state can be based on measuring fieldcorrelators. In the case of the layered structures supporting indirect excitons [6] this means analyzing the excitonic emission of light. In the SF phase the emission intensity tests directly sound-like excitations of excitonic Luttinger liquid similarly to Eq. (23) for N = 2 case. In the CSF phase, the one-photon emission is controlled by the gapped parafermionic modes and, therefore, acquires the threshold similar to the case described by Eq. (27) for N = 2. Both features exhibit threshold singularities as shown in Fig. 2. Thus, the light emission can become a crucial tool for detecting parafermionic excitations.
It is also worth mentioning that, as the system enters the CSF phase, there should appear a special feature in the correlated N -photon emission. Since light is linearly coupled to the excitonic operator and in the CSF phase the algebraic order exists only in the product of N excitonic operators, Eq.(1), such order (entanglement) will be imprinted on N emitted photons. We will consider specific proposal for detecting such N -photon entanglement in greater detail elsewhere.
When this work was prepared for publication we learned about the preprint by Lecheminant and Nonne [28] which results have a substantial overlap with ours.