Topological Phase Diagram of an Interacting Kitaev Chain: Mean Field versus DMRG Study

: In this work, we study the topological phase transitions of a Kitaev chain generalized by the addition of nearest-neighbor Coulomb interaction. We show the presence of a robust topological phase as a function of the interaction strength and of the on-site energy with associated non-zero energy Majorana states localized at the chain edges. We provide an effective mean-field model that allows for the self-consistent computation of the mean value of the local particle number operator, and we also perform Density Matrix Renormalization Group numerical simulations based on a tensor network approach. We find that the two methods show a good agreement in reporting the phase transition between trivial and topological superconductivity. Temperature robustness within a physically relevant threshold has also been demonstrated. These findings shed light on an entire class of topological interacting one-dimensional systems in which the effects of residual Coulomb interactions play a relevant role.


Introduction
In the last two decades, topological superconductivity has attracted growing interest, not least because of its potential role in building decoherence-protected, fault-tolerant quantum electronics that are more stable [1][2][3][4][5][6][7][8][9].The simplest model of topological superconductivity was proposed by Kitaev in 2001 [10].It consists of a p-wave superconducting wire, known as the Kitaev chain, in which two zero energy excitations, called Majorana zero modes (MZMs), nucleate at both edges of the wire.This model reproduces essential features of a Symmetry-Protected Topological State [11][12][13][14][15][16]; in particular, the particle-hole symmetry protects the two MZMs from external perturbations until the bulk gap closes.
MZMs are promising candidates as qubits in topological quantum computers [17][18][19][20][21][22][23] due to their non-abelian anyonic nature, characterized by degenerate energy states and distinct fermion parities.Compared to fermions and bosons, anyons exhibit different statistics, the so-called braid group statistics.When we exchange two non-abelian anyons, we obtain a matrix transformation that encompasses the entire degenerate subspace [24].This unique property of non-abelian anyons allows for the implementation of quantum logic gates through anyon braiding.
Considerable efforts have been made to study how Coulomb interactions influence the properties of one-dimensional systems exhibiting topological phases, including Kitaev chain models.Understanding the robustness of topological phases to interactions [41][42][43][44] has remained an elusive challenge up to now.In particular, previous works have studied the effects of Coulomb interaction beyond the mean-field approximation.Since the Kitaev chain model can be mapped to a spin chain model through the Jordan-Wigner transformation, some works have explored the effects of Coulomb interaction using numerical methods, like the Density Matrix Renormalization Group (DMRG) [45][46][47][48][49][50][51][52].Exact solutions, however, are obtainable only within certain parameter regimes [53,54].Investigations leveraging field theoretical/renormalization group methods [55,56] and bosonization [57], as well as combinations of these methods [58,59], have shed light on critical aspects of these systems.In particular the presence of strong electron-electron interactions destroys the topological phase by suppressing the superconducting gap, and repulsive interactions actually increase the parameter range for which the topological phase exists.
In this paper, we demonstrate that a self-consistent approach to the mean value of the local particle number operator of the interaction is able to provide an effective meanfield model that accurately describes the physical properties of the system.We present a comprehensive analysis of the topological phase diagram of the model, in both repulsive (U > 0) and attractive (U < 0) regimes.Moreover, we also perform DMRG numerical simulations [60] for comparison and to test the bona fides of our protocol.Our results demonstrate that the self-consistent mean-field approach is accurate in describing the topological trivial phase transitions.We also comment on the potential role and limitations of such an approach.This paper is organized as follows: In Section 2, we introduce the Kitaev chain model with the nearest-neighbor Coulomb interaction.This section explains the application of the self-consistent method at both zero and finite temperatures.Section 3 shows numerical results, obtained by the self-consistent method and the DMRG technique.This section also explores how temperature affects the topological phase.Conclusions are drawn in Section 4.

Theoretical Model
We consider the Kitaev chain with nearest-neighbor Coulomb interaction [52,53], a 1D model comprising N interacting spinless fermions with a p-wave superconducting pairing, treated in the mean field.It is described by the following Hamiltonian: where ĉ † j and ĉj are the fermionic creation and annihilation operators, respectively; nj = ĉ † j ĉj is the local number operator that describes the occupation of fermions on each site; µ is the on-site energy; w is the hopping amplitude; ∆ = |∆| e i θ is the superconducting pairing constant; and U is the strength of the Coulomb interaction that can be either attractive (U < 0) or repulsive (U > 0).The superconducting pairing, which creates and annihilates pairs of electrons, leads to the formation of Cooper pairs inside the wire.As a consequence, the particle number is not a conserved quantity of the system, whereas the fermion number parity is conserved.

Mean-Field Approximation
We can treat the nearest-neighbor interaction in the mean-field approximation by assuming that the mean value of the local particle number nj is non-zero and neglecting the fluctuations, δ nj ≡ nj − nj , around their respective mean values up to the first order.Then, the Hamiltonian (1) becomes the following: The mean-field approximation redefines the on-site energy term, resulting in a Dsymmetry class for the new Hamiltonian (2) equivalent to the non-interacting case [10].The additional mean-field terms coming from the decoupling of the interaction correspond to bond operators whose expectation values slightly renormalize the hopping and the superconducting gap and are taken into account in the bare values of t and ∆ [61].
The diagonalization of the Hamiltonian ( 2) is achieved through the Bogoliubov-de Gennes (BdG) formalism [62].We define the Nambu spinor as Consequently, the Hamiltonian is transformed into Ĥ Further, the Hamiltonian is reformulated in the diagonal basis as Ĥ = 1 2 ξ † H BdG ξ. Here represents the Nambu spinor in the diagonal basis, with the indices corresponding to positive energy (particle-like) excitations, and U, the unitary matrix of eigenvectors of H BdG , connects ψ and ξ as ψ = U ξ. (3) The quasiparticle fermionic operators αm and α † m follow canonical anticommutation relations.The diagonal matrix H BdG = U † H BdG U is of size 2 N × 2 N with excitation energies on its diagonal.In this diagonal basis, the Hamiltonian (1) represents free quasiparticles.The Hamiltonian dependence on the mean values of the local particle number nj necessitates a self-consistent approach.In Figure 1, we report a sketch of the self-consistent algorithm in which the convergence is assured by monitoring, for two consecutive steps of the algorithm, "old" and "new", both the self-consistently calculated mean values of the local number operator and the variance: This condition guarantees that the energy spectrum of the excitations converges.The application of the BdG formalism makes it possible to examine in detail the distinctive features of MZMs, in particular their robustness and localization at the chain edges.These aspects are crucial for differentiating the topological from the trivial phase.We use the lowest particle-like excitation energy (minimum absolute value of the excitation spectrum) as an indicator: a null value indicates a possible topological phase, while a nonzero value signals the trivial phase.Furthermore, through the self-consistent method, we can investigate the localization of the MZMs at the edges of the system.This is performed by detecting the peaks of mean value of the local number operator at the extremes of the chain.
In this work, we use a mean-field approach at both zero and finite temperatures, as discussed in the following subsections.

Zero-Temperature Case
At zero temperature, the mean value of an operator coincides with that in respect to the system ground state.When the ground state exhibits a degeneracy m dgs , the mean value of the local particle number is given by nj = 1 m dgs where the sum extends over the degenerate ground states |GS⟩.In our scenario, the system has a single ground state |0⟩ in the trivial phase (where m dgs = 1) and a doubly degenerate ground state |0⟩, |1⟩ in the topological phase (with m dgs = 2).The two states differ in the occupation of a non-local zero energy fermionic mode.The mean values of the local number operator with respect to the ground state ⟨m| nj |m⟩ = ⟨m| ĉ † j ĉj |m⟩ = ĉ † ĉ m (m = 0, 1) are derived using the following relationship between the 2 N × 2 N correlation matrices in the no-diagonal basis ψ and diagonal basis ξ: The elements in Equation ( 6) are N × N submatrices.The correlation matrix in the rhs of Equation ( 6) is known.Since the quasiparticles are free and |1⟩ = α † 0 |0⟩ is the occupied degenerate state in the topological phase, the following relations hold: Then, the mean values of the local number operator are computed self-consistently by selecting the diagonal values of the submatrix ⟨ ĉ † ĉ⟩ m .

Finite-Temperature Case
At a finite temperature, the mean value of the local particle number operator is computed using the component-wise form of Equation (3): Moreover, the mean thermal values of the local quasiparticle number operators are expressed as follows: As a consequence, the mean thermal value of the particle number operator reads as follows: where ) −1 is the Fermi-Dirac distribution function and E n represents the energy of excitations.It is noteworthy that the sum in Equation ( 10) runs over states m with positive energy.The coefficients u j,m and v j,m are determined by solving the BdG equations:

Numerical Results
From the self-consistent method outlined in Figure 1, we can obtain the phase diagram in the (U/w, µ/w) parameter space by using the lowest particle-like excitation energy as an indicator of the existence of a topological phase, as shown in Figure 2a.The method well describes the topological and trivial phases of the system.Indeed, when the system is in the topological phase, zero-energy excitations are present, as depicted in Figure 2b.From Figure 2c, it can be observed that the self-consistent method converges for U/w < 1, though it does not capture the density wave (DW) phase, which occurs when U/w > 1 [52].In the strong coupling limit, in fact, one expects a charge density wave that can be either commensurate or incommensurate with the lattice spacing.These phases are characterized by a modulation of the mean value of the local particle number in the many-body ground state.To detect them, we need either to generalize our mean-field algorithm or rely on the DMRG [Figure 3b] that captures them well.The results showed in Figure 2 are obtained by setting a limit of 500 maximum iterations for the self-consistent loop at a specific point in the phase diagram.We observed that the algorithm converges quickly, often within a hundred iterations, at points within the trivial and topological phases.However, along the transition line between the topological and trivial phases, the convergence of the algorithm is slower or not reached.
Additionally, we noticed that for some points it is necessary to increase the number of iterations to achieve convergence.To verify the validity of the self-consistent method used, we compare the previous results with those of the DMRG.
In Figure 3a we derive the phase diagram by the DMRG technique.To construct the full phase diagram, it is essential to look beyond the energy difference between the many-body ground state and the first excited state.A notable feature of the DW phase that distinguishes it from the trivial and topological phases is the equal fermionic parity between the ground state and the first excited state [52].Thus, if the difference in fermionic parity between them is zero, the system is in the DW phase.If not, it is in the trivial or topological phase.To distinguish between these last two phases, we use as an indicator the energy difference between the many-body ground state and the first excited state: if this is zero, the system is in the topological phase, whereas if it is non-zero, it is in the trivial phase.
Figure 4 shows a comparison between the self-consistent method at zero temperature and the DMRG.Fixing the on-site energy µ/w (Columb interaction U/w), as the Coloumb interaction U/w (on-site energy µ/w) varies, both the lowest particle-like energy excitation in the zero-temperature self-consistent method and the energy difference between the first excited state and the ground state calculated by the DMRG are in significant agreement.

Conclusions
We have presented an analysis of the topological phase diagram of a Kitaev chain generalized by the addition of Coulomb interactions.We have based our analysis on a mean-field self-consistent approach where the mean value of the local particle number operator is self-consistently fixed.We have found that topological phases are stable in a certain range of attractive interactions and on-site energies.A comparison with DMRG numerical simulations is also reported to complement the phase diagram of the model, showing the emergence of a DW phase.Our results indicate that the self-consistent meanfield approach is a reliable method to investigate the topological phase diagram of the interacting Kitaev chain, requiring less numerical effort, and thus could be usefully applied to other one-dimensional topological interacting models.

Figure 1 .
Figure1.Flowchart of the self-consistent method used to calculate the excitation spectrum and mean value of the local number operator ⟨ n⟩ new .The way ⟨ n⟩ new is computed differentiates the zero-temperature case from the finite-temperature case.

Figure 2 .Figure 3 .
Figure 2. Results obtained from the self-consistent method at zero temperature: (a) Phase diagram in the (U/w, µ/w) parameter space. (b) Mean value of the local number operator at varying chain sites in both trivial and topological phases, indicated by the red and purple crosses in the phase diagram in panel (a).(c) Maximum value of the convergence condition regarding the mean values of the local number operator, between the last step n and the penultimate step n − 1. System parameters are N = 50, w = |∆| = 1, and ε = 10 −10 .The comparison in Figure4bhighlights a subtle shift between the two methods, observed for the specified fixed on-site energy µ/w = 1.8 and Coulomb interaction U/w = −0.1Finally,we have investigated the effect of the temperature on the topological phase transition and verified that the transition line is not affected by temperatures up to k B T of order w by verifying that the lowest eigenvalues are independent from T.

Figure 4 .
Figure 4. Comparison between the self-consistent zero-temperature method and DMRG: (a) Phase diagram in the (U/w, µ/w) parameter space. (b) Lowest particle-like energy excitation in the self- consistent (SC) method in green and energy difference between the first excited state and ground state calculated by DMRG in red.In the top (bottom) panel, we fix the on-site energy (Coloumb interaction) and vary the Coulomb interaction (on-site energy); see the purple (orange) dashed line in the phase diagram in panel (a).We fix for the tolerance of the self-consistent method the value ε = 10 −10 .System parameters are N = 50 and w = |∆| = 1.