Majorana edge states in Kitaev chains of the BDI symmetry class

I read the revised manuscript from the author and I acknowledge that he made a serious effort to address mine and the other referee concerns. Unfortunately, part of these modifications are bound to generate confusion rather than clarify the open issues. In my previous report I requested the author to clarify the issue of the decay of the Majorana edge states for power-law interactions with \alpha,\beta>2 due to the mismatch between the result of the author (and of Ref. [20]) with the ones of Ref. [13]. The author, in agreement with Ref.[20], obtains the power-law decay in Eq.(28)

s_{l}\propto l^{-\delta} with \delta=min(\alpha,\beta).

Ref.[13] finds that the above result is only applicable for \alpha<3 or \beta<2, while at larger decays the approximation of Ref.[13] yields that long-range interactions are in-influential to the edge modes decay.

In the revised version of the manuscript the author claims that the result in Eq.(28) is the leading one and that all additional terms are subleading. In order to justify the discrepancy with the finding of Ref.[13] the author also makes some general claims in a footnote (Ref. [29]), where he states that Ref.[13] is "most likely erroneous". In order to sustain this latter claim the author states that assuming $u_{l},v_{l}\propto exp(-\lambda l)$ with $Re(\lambda)>0$ yields a diverging result in Eq.(2) .

Actually, it is the claim of the author to be erroneous. In fact, if one assumes the exponential decay of the BdG amplitudes with positive decay exponent the summation on the l.h.s. of Eq.(2) is convergent and, accordingly, cannot serve as an argument against the scattering approach. I advice the author to remove his fallacious statement and amend the footnote accordingly. Also, I would advise the author to be more cautious when implicating possible mistakes in other published works.

More in general, it should be noted that all approaches of Ref. [13], of Ref.[20] and of the current manuscript, albeit using different methods, produce similar equations. The main discrepancy between the approaches of Ref. [13] and Ref.[20] depends on the way the polylogarithm singularities in the complex plane are treated. Indeed, Ref.[20] always keeps the non-analytic term in the expansion of the PolyLog. On the contrary, Ref.[13] only focuses on the leading order contribution, which becomes analytic for $\alpha>3$ and $\beta>2$.

The author's result is specular to the one of Ref.[20], but the author states that Eq. (28) can only be applied to the case of non-integer $\alpha,\beta$. This is in agreement with the aforementioned picture in which the non-analyticity of the polylog function plays a crucial role in the Wiener-Hopf method. However, this result leaves open two alternative scenarios:

a) For integer power-laws e.g. $\alpha=6$ and $\beta=4$, the decay of the modes is still power-law.

b) The power-law behaviour only occurs at non-integer values and the integer case displays exponential behaviour.

The author stated that his analysis unambiguously identifies Eq. (28) as the leading behaviour for the decay of the edge mode, but what scenario is to be expected? In the case (a) the author shall clarify how the power-law behaviour is recovered in absence of the polylog branch cut. In the case (b) the result will be quite unphysical since all integer power-law decays will behave differently from the non-analytic case. These question is not physically irrelevant since most of the naturally occurring power-law potentials display integer power-laws (gravity, coulomb interaction, dipolar, van-der-Wals,...).

In the revised version of the manuscript the author also includes some finite-size numerical study of the edge mode decay. This is certainly commendable. Nonetheless, the discussion of these results is incomplete and it generates several further issues rather than clarifying the aforementioned questions.

Issue 1: The author states: "When $\alpha< 1$ the slopes exhibit
significant variations with site number $l$, however, the
variations become less pronounced with increasing chain
length (compare the blue dotted and solid curves, corresponding to L = 10000 and L = 20000, respectively)". The poor convergence of the numerical data does not come as a surprise to me, since strong long-range couplings $\alpha,\beta<1$ are known to generate large finite size corrections. Yet, the numerical curves in Fig.1 clearly do not converge to a constant slope and the author only presents the study of a single size (apart from 1 set of values where both $L=10^{4}$ and $L=2*10^{4}$). This is far from being enough to argue convergence of the numerics to the analytics. I suggest the author to present some form of finite size scaling for various (exponentially increasing) system sizes 2^10, 2^{12}, 2^{14}, and so on, to show that the central slope is converging to the expected value of Eq. (28).

Issue 2: While I can live with some finite size corrections at $\alpha,\beta<1$, the black curve in Fig.1 with $\alpha=1.5$ and $\beta=3$ also does not reproduce the result in Eq. (28). Worst of all, the curve seem to converge to a fixed slope, but the slope does not coincide with the theoretical prediction. On this matter the author only comments: "The curve corresponding to $\alpha = 1.5$ exhibts a relatively stable slope, which is somewhat smaller than the one predicted analytically [Eq. (28)]", but this is far from being sufficient. Eq. (28) was supposed to be exact (at least at $1<\alpha,\beta<2$) how can the numerics not reproduce it? Why is the result smaller than the prediction? Are those very weak finite size effects? I urge the author to present a finite size scaling for the mode decay also in this region and present convincing evidence that Eq. (28) holds in the thermodynamic limit.

Issue 3: Since the author is performing numerical analysis, I would reccommend doing so also in the region $\alpha>3$ and $\beta>2$ where the aforementioned controversy between Eq. (28) and the exponential decay took place. The author may also try some large integer power laws, e.g. \alpha=6 and \beta=4 to investigate what one may expect in the limit of the applicability of Eq. (28).

Unfortunately, in the present form the paper is not suitable to be published, but a more extensive numerical study and a softening of the current harsh claims may make it a valuable addition to the literature.

-) Amend/soften incorrect/harsh statements in the footnote Ref.[29].

-) Comment on the destiny of the edge mode power-law tails for large integer $\alpha$ and $\beta$.

-) Numerical finite size scaling of the edge mode decay at $\alpha,\beta<1$ to verify its convergence to Eq. (32).

-) Numerical finite size scaling of the edge mode decay at $\alpha,\beta>1$ to verify its convergence to Eq. (28).

-) Numerical finite size scaling of the edge mode decay at $\alpha=6,\beta=4$ or large values to verify the possible limitations of Eq. (28).

1. The analytical results are strong
2. Good presentation and comparison with previous results

1. Numerical checks are not properly explained and interpreted

The authors has carefully considered the Referees' suggestions and changed the paper accordingly.
In particular, now recent literature is correctly mentioned and acknowledged.
I think the theoretical and analytical calculations are now complete.

However I have some questions regarding the added part on numerical evaluations of s_l and E_0(l).
Finite size effects can be very strong (especially when alpha<1) so that it is necessary to consider very long chains and values of l far from the (right) boundary of the finite chain considered in numerics. Also the theoretical predictions holds for l large. i.e. far from the left boundary l=0. This can be clearly seen in Fig 1(a), from which one can reasonably argue that the correct range to perform any fit lies between l=40 and l=5000, at least for L=20000.
This is indeed the range considered for the evaluation of the energy E_0(l) (Fig.2), which exhibits results agreeing with the theoretical predictions. Instead, the author has not used the same procedure for s_l and I wonder why a similar fit (in the same range) is not reported for curves such as the one reported in Fig 1(a), for different values of alpha.
Also, it would be appropriate to have an estimate of the error of the fit.
Finally, I have doubts on the relevance of the behaviour of d(ln s_l^2)/d(ln l) as function of l itself -shown in Fig 1(b). This quantity is calculated locally from a quantity such as s_l that on a chain of length L is for sure affected strongly by both the finiteness of the chain itself and of the lattice spacing.

Analysis of the validity of the numerics