# Quantum entropies of realistic states of a topological insulator

Nicolás Legnazzi<EMAIL_ADDRESS>Omar Osenda
<EMAIL_ADDRESS>Facultad de Matemática, Astronomía, Física y
Comptación, Universidad Nacional de Córdoba, and Instituto de Física Enrique
Gaviola, CONICET, Av. Medina Allende s/n , Ciudad Universitaria, CP:X5000HUA
Córdoba, Argentina

###### Abstract

Nanowires of BiSe show topological states localized near the surface of the
material. The topological nature of these states can be analyzed using well-
known quantities. In this paper, we calculate the topological entropy
suggested by Kitaev and Preskill for these states together with a new entropy
based on a reduced density matrix that we propose as a measure to distinguish
topological one-electron states. Our results show that the topological entropy
is a constant independent of the parameters that characterize a topological
state as its angular momentum, longitudinal wave vector, and radius of the
nanowire. The new entropy is always larger for topological states than for
normal ones, allowing the identification of the topological ones. We show how
the reduced density matrices associated with both entropies are constructed
from the pure state using positive maps and explicitly obtaining the Krauss
operators.

††preprint: APS/123-QED

##  I  Introduction

Topological states appear in different materials and geometries [1]. Their
appearance changes the conducting properties of the material, giving place to
many distinctive phenomena like the quantum spin Hall effect [2, 3, 4],
topological superconductors [5, 6], etc. The compound
$\mbox{Bi}_{2}\mbox{Se}_{3}$ is a 3D topological insulant [1, 7, 8, 9].

The calculation of the spectrum, spin currents, and density of states (DOS) of
BiSe cylinders or quantum dots is the subject of numerous works [10, 11, 12,
13, 14]. Besides first principle calculations, the kp model is a preferred
tool to calculate the spectrum and eigenstates necessary to obtain the spin
currents and the DOS [7]. Once the eigenstates are available, it is possible
to characterize the topological states using quantum entropies to study
different physical regimes.

The cylindrical geometry is well-suited to calculating the Kitaev-Preskill
topological entropy [15]. This quantity should be a constant, independent of
all the quantum numbers that characterize a given topological state. It can
only depend on the topology of the problem. But what are the values of this
entropy for the normal in-band eigenstates that are eigenfunctions of the kp
model? In a situation without a boundary, the topological entropy of a normal
state should be null, but a cylindrical nanowire necessarily has one. So, the
topological entropy for a normal eigenstate obtained using the kp method could
be non-null, but its value should depend on its quantum numbers.

Distinguishing topological from normal states can be done using different
quantum entropies and related quantities, such as the entanglement spectrum
[16, 17, 18, 19, 20]. To this end, we propose a particular reduced-density
matrix whose entropy distinguishes the one-electron topological states of a
$\mbox{Bi}_{2}\mbox{Se}_{3}$ cylinder from the non-topological ones. Our study
differs from those that employ entanglement entropy to study many-electron
wave functions. The entanglement entropy detects the non-local character of
topological states when applied, for instance, to quantum states that are good
approximations of the many-electron ground-state wave function of the
fractional Quantum Hall effect, that is, the Laughlin states [21, 16, 22].

For a given eigenstate of the kp Hamiltonian, which we calculate as a
superposition of a basis set functions as is usual in the Rayleigh-Ritz
variational method, the proposed reduced density matrix depends on the
coefficients of the superposition and integrals of the basis set functions.
The Rayleigh-Ritz variational method accurately provides the band structure
near the gap between the conduction and valence bands in semiconductor
nanostructures when applied to kp Hamiltonians. The application of the
Rayleigh-Ritz variational method to kp Hamiltonians allowed the study of
electronic and optical properties of core-shell nanowires [23, 24], edge
states with and without an external magnetic field applied to a quantum well
[25, 26, 27, 28], the transition between resonance and bound states in quantum
dots embedded in nanowires [29], the entanglement entropy of edge states in
quantum wells [30], spin currents in topological insulators [10], amongst
other physical phenomena.

The paper is organized as follows. In Section 2, we present the kp
Hamiltonian, whose eigenvalues give the band structure of a cylindrical
nanowire made of $\mbox{Bi}_{2}\mbox{Se}_{3}$ and briefly describe how to
obtain a numerical approximation to the spectrum and eigenvectors using the
Rayleigh-Ritz method. Section 3 deals with the calculus of the topological
entropy of an approximate variational eigenvector. The topological entropy
depends on the von Neumann entropy of several reduced-density matrices.
Obtaining each one of these matrices implies tracing out a real-space
partition from a pure density one [15, 31].

We present in Section IV the mode-dependent reduced density matrix (RDM),
$\rho_{MD}$, which contains information about a variational state through the
coefficients of the variational expansion and spatial integrals of the basis
set functions. We show how the von Neumann entropy and the entanglement
spectrum [16] of the mode-dependent RDM allow us to distinguish between non-
topological and topological states.

In Section V, we use quantum state processes [32, 33] as an alternative way to
construct the mode-dependent RDM and some of the RDM necessary to calculate
the topological entropy. Using the quantum process tomography algorithm
described in Reference [34], the quantum process results in a sum of Kraus
operators determined by a gradient-descent algorithm. The quantum process
numerically calculated predicts an RDM, $\tilde{\rho}$, slightly different
from the one intended, $\rho^{\prime}$. We compare both matrices, calculating
their fidelity. Finally, in Section IV, we summarize and discuss our results.

##  II  Model and Hamiltonian

We consider a cylindrical nanowire made of $Bi_{2}Se_{3}$, with a constant
radius on the tens of nanometers and infinitely long in the axial direction.

To obtain the band structure and eigenstates, we employ the
$\mathbf{k}\cdot\mathbf{p}$ Hamiltonian given by

$H_{0}=\varepsilon(\mathbf{k})+\left(\begin{array}[]{c c c
c}{{M(\bf{k})}}&{{B(k_{z})k_{z}}}&{{0}}&{{A\left(k_{\parallel}\right)k_{-}}}\\\
{{B(k_{z})k_{z}}}&{{-M(\bf{k})}}&{{A(k_{\parallel})k_{-}}}&{{0}}\\\
{{0}}&{{A\left(k_{\parallel}\right)k_{+}}}&{{M(\bf{k})}}&{{-B(k_{z})k_{z}}}\\\
{{A\left(k_{\parallel}\right)k_{+}}}&{{0}}&{{-B(k_{z})k_{z}}}&{{-M(k_{z})k_{z}}}\end{array}\right),$
(1)

where

  * •

$\varepsilon(\mathbf{k})=C_{0}+C_{1}k_{z}^{2}+C_{2}k_{\parallel}^{2}$,

  * •

$M(\mathbf{k})=M_{0}+M_{1}k_{z}^{2}+M_{2}k_{\parallel}^{2}$,

  * •

$B(k_{z})=B_{0}$,

  * •

$A(k_{\parallel})=A_{0}$.

The Hamiltonian in Equation 1 was introduced by Zhang [8] and collaborators to
adjust the band structure found for different materials showing topological
states. We consign the parameters that define the Hamiltonian in Equation 1 in
Table 1.

Parameter | Units |
---|---|---
$C_{0}$ | $eV$ | -0.0068
$M_{0}$ | $eV$ | 0.28
$C_{1}$ | $eV\cdot A$ | 1.3
$M_{1}$ | $eV\cdot A$ | -10.0
$C_{2}$ | $eV\cdot A^{2}$ | 19.6
$M_{2}$ | $eV\cdot A^{2}$ | -56.6
$A_{0}$ | $eV\cdot A^{2}$ | 4.1
$B_{0}$ | $eV\cdot A$ | 2.2
Table 1: The Table contains the values of the parameters that enter into the
four band $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian of a
$\mbox{Bi}_{2}\mbox{Se}_{3}$ nanowire.

We calculate approximate eigenvalues and eigenfunctions using the Rayleigh-
Ritz variational method, which reduces the eigenvalue problem

$Hf=ef,$ (2)

where $H$is the $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian, to an algebraic one.
The Rayleigh-Ritz method is suitable for calculating the band structure of
nano-structures near the gap between the conduction and valence bands and for
states lying inside it, for instance, to find the topological states in three-
dimensional topological insulators, to study the transition between localized
and resonance states in quantum dots, to study properties of states in quantum
wells in the Quantum Spin Hall Effect regime, etc. We consign the details
about the implementation of the method to Appendix A.

##  III The topological entropy

The study of topological states has led to different entropic-like quantities
as, for instance, the real-space entropy, the topological entropy of Kitaev-
Preskill, and so on. Also, it is worth mentioning the entanglement spectrum
or, in the case of studies dealing with topological states in spin chains, the
Renyi entropies.

Figure 1: The cartoon depicts the main geometric elements that enter into the
calculus of the topological entropy and other quantities in Sections 2, 3, 4,
and 5. $R$ is the radius of the nanowire, $R_{c}$ is the radius associated to
the circular sectors $A$, $B$, and $C$. The central angle of each sector is
equal to $2\pi/3$. The blue dashed line corresponds to $R$ while the red one
corresponds to $R_{c}$.

Tracing out different subspaces of the whole Hilbert space leads to one or
other entropy. For instance, for a given multiparticle wavefunction $\psi$,
the real-space entropy

$S_{A}(\rho_{A})=-\mbox{Tr}(\rho_{A}\log{(\rho_{A})}),$ (3)

is calculated from the reduced density operator defined by

$\rho_{A}=\frac{1}{Z}\int_{\bar{A}}\,\psi^{\star}(\vec{x})\psi(\vec{x})d\vec{x}^{\prime},$
(4)

where $\bar{A}$ is the spatial region outside region $A$, $Z$ is a
normalization constant, and the integral involves a subset $\vec{x}^{\prime}$
of all the particle coordinates.

From the definition of the reduced density operator above, it is clear that
the value of the real-space entropy could depend on the size of the spatial
region. To avoid this problem, Kitaev and Preskill proposed to trace out over
a set of spatial sectors and combine the corresponding real-space entropies to
single out the topological entropy of the quantum state.

Following the argument in Reference [15], we consider three regions, $A$, $B$,
and $C$, as shown in the cartoon in Figure 1. The triangular sectors are
defined by

$\displaystyle
A:\varphi\in\left[-\pi/6,\pi/2\right]\quad\mbox{and}\quad\rho\leq R_{C},$ (5)
$\displaystyle
B:\varphi\in\left[7\pi/6,11\pi/6\right]\quad\mbox{and}\quad\rho\leq R_{C},$
(6) $\displaystyle
C:\varphi\in\left[\pi/2,7\pi/6\right]\quad\mbox{and}\quad\rho\leq R_{C}.$ (7)

Tracing out the spatial region outside $A$, $B$, $C$, or combinations of them,
the topological entropy is given by

$S_{t}=S_{A}+S_{B}+S_{C}-S_{AB}-S_{BC}-S_{AC}+S_{ABC}.$ (8)

In our case, we calculate $\rho_{A}$ as follows

$\rho_{A}=\int_{D}|\psi(\rho,\varphi,z)\rangle\langle\psi(\rho,\varphi,z)|\;\rho\,d\rho\,d\varphi,$
(9)

where $|\psi\rangle$ is one of the eigenvectors obtained using the Rayleigh-
Ritz method, and the integral is over the surface of the disk of radius $R$
minus the triangular sector $A$, see Figure 1. Proceeding in this way,
$\rho_{A}$, $\rho_{B}$, and so on, are given by $4\times 4$ matrices. The
corresponding von Neumann entropies are given by

$S_{A}(\rho_{A})=-\sum_{i}\lambda_{i}\log(\lambda_{i}),$ (10)

where the $\lambda_{i}$ are the eigenvalues of $\rho_{A}$.

The variational eigenvectors $|\psi\rangle$ are complex column vectors that
depend on the three cylindrical coordinates and are labeled by two quantum
numbers, $k_{z}$ and $L$,

$|\psi_{L,k_{z}}\rangle=\sum_{n}^{N}\left(\\!\\!\\!\begin{array}[]{c}{{b_{L,k_{z},n,\uparrow}A_{L,n}J_{L}\left(\alpha_{n}^{L}\rho/R\right)e^{iL\varphi}}}\\\
{{c_{L,k_{z},n,\uparrow}A_{L,n}J_{L}\left(\alpha_{n}^{L}\rho/R\right)e^{iL\varphi}}}\\\
{{b_{L,k_{z,n,\downarrow}}A_{L+1,n}J_{L+1,n}\left(\alpha_{n}^{L+1}\rho/R\right)e^{i(L+1)\varphi}}}\\\
{{c_{L,k_{z,n,n,\downarrow}}A_{L+1,n}J_{L+1,n}\left(\alpha_{n}^{L+1}\rho/R\right)e^{i(L+1)\varphi}}}\end{array}\\!\\!\\!\right)e^{ik_{z}.z},$
(11)

where $c_{L,k_{z},n,\uparrow}$, $c_{L,k_{z},n,n,\downarrow}$,
$b_{L,k_{z},n,\uparrow}$ and $b_{L,k_{z},n,\downarrow}$ are the linear
variational coefficients, $J_{L}$ is the Bessel function with index $L$, the
$\alpha_{n}^{L}$ are its roots, and $A_{L,n}$ is a normalization constant.
Note that the number of roots employed coincides with the number of basis set
functions used.

Lou and collaborators [10] used the variational basis set in Equation (11) to
calculate the band structure of $\mbox{Si}_{2}\mbox{Be}_{3}$ nanowires. They
focused on nanowires with a diameter of 120 nanometers, and for these
nanowires, they showed that there are two topological states for each pair
$(k_{z},L)$. Of the two topological states, one has energy closer to the
conduction band, while the other lies closer to the valence band. We will
refer to the former state as the upper energy state and to the latter as the
lower energy state, respectively. We will also focus on nanowires with a
diameter of 120 nanometers for comparison purposes.

As is usual with states or eigenvalues calculated using the Rayleigh-Ritz
variational method, the larger the number of the basis set functions, the
greater the accuracy. If $\rho^{m}$ is an RDM obtained from a variational
eigenvector with $m$ basis set functions, we assess the accuracy of the
variational method by comparing two RDM with different values of $m$. The
fidelity between a succession of $\rho_{A}$ matrices obtained for one of the
two topological states with $L=0$, each matrix calculated with a different
number of basis set functions, is shown in Figure 2.

Figure 2: The fidelity between two RDM,
$\mathcal{F}(\rho^{m},\rho^{m^{\prime}})$ vs the basis set size $m$. The
values of $m$ and $m^{\prime}$ are picked as successive values of the set
${20,24,..60,64}$, with $m>m^{\prime}$. The solid blue dots correspond to the
calculated values, while the black line is a guide and corresponds to the
unity value.The reduced density matrix is the $\rho_{ABC}$ obtained for the
topological state with $L=0$, $k_{z}=0.1$, and larger energy.

Note that, while each matrix on the succession is normalized,
$\mbox{Tr}(\rho)=1$, the fidelity approaches the unity for larger basis sets.
This convergence indicates that the successive matrices are more and more
similar. Consequently, in what follows, we will present results obtained with
the basis set with $N=40$. The convergence of the energy values of the
topological states is even better. For given values of $k_{z}$ and $L$, the
figures obtained using the two larger basis sets have a relative difference of
less than $1\times 10^{-4}$.

The topological entropy in Equation (8) should be independent of the quantum
numbers that label the topological quantum state. Figure 3a shows this
property, where we plot the topological entropy of all the variational
eigenstates obtained for different values of $k_{z}$ and $L$ as a function of
the wave number $k_{z}$. The colored lines correspond to topological states
with $L=0,1,2,3$. The topological entropy of the different topological states
is constant as a function of $k_{z}$, except when its energy is close enough
to a band, in which case the entropy value drops to zero.

Figure 3: The panels show the topological entropy and eigenvalues obtained
using the variational method. a) the panel shows the modulus of the
topological entropy $S_{t}$, Equation 8, of topological states with quantum
numbers $L=0,1,2$ and $3$ versus the wave number $k_{z}$. The dashed line
corresponds to $\gamma=\log{2}$. It is clear that $|S_{t}|\sim\log{2}$ for an
ample range of $k_{z}$ values. b) The variational spectrum near the gap. The
black curves correspond to the valence and conduction bands, while the colored
ones correspond to the topological states. Accordingly, these energies lie
inside the gap. For each value of $k_{z}$ and $L$, there are two topological
states, one with an eigenvalue nearer to the conduction band (the upper energy
value) and the other with its eigenvalue closer to the valence band (the lower
energy value). Comparing both panels, it is clear that the sudden drop in the
value of the topological entropy occurs when the energy of the corresponding
state is close enough to the bands. c) The modulus of the topological entropy
calculated for all the eigenstates whose eigenvalue we plotted in b).

In panel b) we show the behavior of the variational eigenvalues near the gap.
The color code used to depict the energy of the topological states is the same
one used in panel a) for the corresponding topological entropies. The black
curves correspond to the valence and conduction bands and contain all the
other eigenvalues obtained from the variational method. These are the
eigenvalues of the normal non-topological states.

Comparing both panels, a) and b), it is clear that the topological entropy
drops to zero when the energy of the topological state is close to the
conduction or the valence bands.

The topological states reduced density matrices are difficult to calculate
accurately close to $k_{z}=0$. Because of this, the topological entropy curves
show a little slump near $k_{z}=0$. The difficulty increases for larger values
of $L$ and, consequently, the deeper one corresponds to the larger $L$.

Finally, in Figure 3c, we show the topological entropy calculated for all the
states whose eigenvalues we plotted in panel b). Again, the colored curves
correspond to the topological states, while the black ones correspond to the
normal non-topological ones. There are a few salient traits worth commenting
on. The topological entropy of the non-topological states depends on $k_{z}$,
and eventually, its value reaches up to $\log{2}$. Since we are dealing with a
finite system, the length scale associated with the different eigenstates is
finite and comparable with $R_{c}$, except for the topological states. So,
changing the value of $R_{c}$ would change the topological entropy of a non-
topological state but not the topological entropy of a topological one.

The topological entropy should also be independent of the length of the
circular curve that defines the triangular sectors $A$, $B$, and $C$ or,
equivalently, of the radius $R_{c}$. Figure 4 shows the topological entropy as
a function of the wave number for three different values of $R_{c}$. As the
Figure clearly shows, $S_{t}$ is a constant with very high precision.

Figure 4: The Figure shows the behavior of the modulus of the topological
entropy for the two topological states with quantum number $L=0$ as a function
of the wave number $k_{z}$. The green curves in each panel correspond to the
topological entropy calculated with a particular $R_{c}$ value. a)
$R_{c}=R/3$, b) $R_{c}=R/2$, and c) $R_{c}=2R/3$. The black curve in each
panel depicts the behavior of the topological entropy calculated for a non-
topological state $\psi^{n}(k_{z},L)$, where $n$ indicates that the state has
an eigenenergy $E_{n}$, with $n$ fixed and such that the eigenenergy lies well
inside the valence/conduction band. It is easy to appreciate that the
topological entropy of the topological states remains constant, irrespective
of $R_{c}$. On the other hand, the entropy corresponding to the non-
topological states changes its value when $R_{c}$ changes. .

##  IV Mode-dependent reduced density matrix, its entropy and entanglement
spectrum

As the results in the previous Section show, the topological entropy is
effectively a constant for the topological states up to the numerical
precision. Nevertheless, this entropy is non-null for the non-topological
ones. Its value depends on both quantum numbers $k_{z}$, $L$, and the radius
$R_{c}$. Note that $\log(2)$ bounds from above the values that achieve
$|S_{t}(\rho_{nt})|$, where $\rho_{nt}$ is the density matrix of a non-
topological state. For some values of $k_{z}$, $L$ and $R_{c}$
$|S_{t}(\rho_{nt})|=log(2)$, which renders the Kitaev-Preskill entropy useless
in this context to distinguish a topological state from normal ones.

In this Section, we propose a reduced-density matrix whose entropy for
topological states is larger than the entropy of the normal ones. The reduced-
density matrix depends on the variational coefficients and spatial integrals
of the variational basis set functions.

Writing the variational eigenfunctions as

$|\psi\rangle=\sum_{n,j}c(j,L,n)f_{j,L,n}(\rho,\varphi)e^{ik_{z}z}|j\rangle,$
(12)

we define the mode-dependent reduced density matrix $\rho_{MD}$, with matrix
elements

$\displaystyle\left[\rho_{MD}\right]_{j,L,n;i,L^{\prime},m}=$
$\displaystyle\frac{1}{\mathcal{N}}c(j,L,n)(c(i,L^{\prime},m))^{*}\times$
$\displaystyle\int_{\Omega}f_{j,L,n}(\rho,\varphi)(f_{i,L^{\prime},m})^{*}(\rho,\varphi)\,\rho\,d\rho\,d\varphi,$
(13)

where $\mathcal{N}$ is a normalization constant such that
$\mbox{Tr}(\rho_{MD})=1$, $\Omega$ is the ring outside $ABC$, and $\rho_{MD}$
is a $4N\times 4N$ matrix.

The mode-dependent RDM is closely related to the pure state given by

$\rho_{p}=|\psi\rangle\langle\psi|,$ (14)

which has matrix elements

$\rho_{p}=c(j,l,n)(c(i,l^{\prime},m))^{*},$ (15)

in the variational basis and is also a $4N\times 4N$ matrix.

In the next Section, we show that $\rho_{MD}$ is also given by

$\rho_{MD}=\sum_{j}K_{j}\rho_{p}K^{\dagger}_{j},$ (16)

i.e., $\rho_{MD}$ results from applying a quantum process to $\rho_{p}$, but
before we want to analyze its spectrum and the behavior of the von Neumann
entropy of $\rho_{MD}$ as a function of $k_{z}$.

Figure 5: The Figure shows the mode-dependent quantum entropy and the spectrum
of the mode-dependent reduced density matrix. a) The mode-dependent entropy
for all the eigenstates shown in Figure 3b versus the wave number $k_{z}$. For
the topological states, depicted with colored lines, the entropy is
independent of the quantum number $L$ and shows a smooth dependence with
$k_{z}$. Again, when the energy of the topological states becomes close enough
to the bands, the value of the entropy drops to zero. b) The spectrum of the
mode-dependent RDM. The eigenvalues are ordered from larger to smaller. The
spectrum of the mode-dependent RDM for topological states and the normal non-
topological ones correspond to the green and black lines, respectively. The
entropies shown in a) were calculated using these spectra. Note that the
vertical scale is logarithmic, and the different behavior exhibited by both
kinds of spectra.

Figure 5 shows a) the von Neumann entropy of the mode-dependent RDM,
$(S\rho_{MD})$, and b) the spectrum of $\rho_{MD}$ as functions of $k_{z}$ and
the eigenvalue number, respectively. The data in panel a) shows that the
entropy of the topological states has larger values than the entropy of the
normal ones, except when the energy of the topological states is too close to
the bands. Moreover, up to the numerical accuracy, its value depends only on
the quantum number $k_{z}$, but not depends on the angular momentum $L$.

On the other hand, Figure 5 b) shows the mode-dependent RDM spectrum obtained
corresponding to all the variational eigenvectors calculated with $N=70$. Note
that the vertical scale is logarithmic. The spectra of the topological states
show a decaying compatible with $\lambda_{k}\sim e^{-\alpha k}$ over an ample
range of values, where $\alpha$ is a constant. This behavior points to the
topological character of quantum states in different physical systems. Figure
6 shows the behavior of the spectrum of the mode-dependent RDM calculated
using three different basis set sizes.

Figure 6: The figure shows the spectrum of the mode-dependent reduced density
matrix calculated for a topological state with quantum numbers $L=0$, and
$k_{z}=0.1$, with three different basis set sizes $N$. The black curve
corresponds to $N=30$, the red to $N=50$, and the blue to $N=70$,
respectively. The black dashed line corresponds to an exponential function and
is a guide to the eye.

The black, red, and blue curves shown in Figure 6 correspond to basis sets
with $N=30$, $50$, and $70$ functions, respectively. The dashed black line is
a guide to the eye. The number of eigenvalues that show an exponential decay
grows with the number of basis functions used to obtain the RDM, although the
increase is slow. This trait is another manifestation of the difficulties
inherent to the obtention of the mode-dependent RDM or, more generally, how
difficult it is to numerically calculate other quantities related to
topological states beyond their spectrum. The entanglement spectrum, which is
given by

$\zeta_{k}=-\log(\lambda_{k}),$ (17)

together with the exponential decay of $\lambda_{k}$ result in

$\zeta_{k}=c+\alpha k.$ (18)

The abrupt decay of the eigenvalues beyond the exponential regime, shown in
Figure 6, marks when the numerical approximation becomes inaccurate, rendering
all the remaining values indistinguishable from zero.

##  V Quantum process tomography

The results in the previous Section suggest different choices for the reduced
density matrices if we wish to distinguish between topological and normal
states using entropic-like quantities, at least when dealing with approximate
ones obtained from phenomenological Hamiltonians.

For a given known quantum state, we could obtain the mode-dependent RDM using
different basis sets, not only the particular one that we employ to implement
the variational method. So, up to a point, some features of the entropy of the
mode-dependent RDM should depend on the basis set chosen.

Instead of implementing our calculations in a different basis set, we prefer
to show that the mode-dependent RDM is equivalent to many others by showing
how to obtain it from the pure state in Equations 13 and 14 using quantum
processes. Doing this has the twofold purpose of not dealing with the
complicated calculations inherent to a change of functions basis set and
showing that RDM obtained through operational processes leads to entropies
that distinguish topological from normal states.

A quantum process takes a given RDM $\rho$ to another one, $\rho^{\prime}$ as
follows

$\rho^{\prime}=\mathcal{E}(\rho).$ (19)

The quantum process $\mathcal{E}$ is a linear superoperator, and $\rho$ and
$\rho^{\prime}$ as operators do not necessarily act on Hilbert spaces of the
same dimension. The particular $\mathcal{E}$ that relates a given pair of
$\rho$ and $\rho^{\prime}$ is determined using quantum process tomography,
which usually is an expensive calculation.

There are numerous algorithms to calculate $\mathcal{E}$ [34, 35, 36, 37, 38,
39, 40, 41, 42, 43], which depend on the representation used for the process
and the dimensionality of the involved RDM. We use the method proposed by
Ahmed et al. [34], which assumes a Kraus representation for the quantum
process

$\mathcal{E}=\sum_{j}K_{j}\rho K_{j}^{\dagger},$ (20)

and that the dimensionality of both Hilbert spaces, where $\rho$ and
$\rho^{\prime}$ act, is a power of two. Ahmed et al. define a cost function
that depends on the two RDM, the Kraus operators, $K_{j}$, and learn them
using a gradient-descent method.

Starting with randomly chosen initial Kraus operators, the non-negative cost
function $\mathcal{L}$ is minimized up to values near zero to learn the
optimal Kraus operators. See Appendix B for details about the calculation. We
studied two cases. For the first case, $\rho$ is the pure state in Equation
(14), and $\rho^{\prime}=\rho_{ABC}$. For the second one, $\rho$ is again the
pure state in Equation (14), and $\rho^{\prime}$ is the mode-dependent RDM.

The cost function is given by

$\mathcal{L}(\mathbb{K})=\sum_{j}\left[d_{j}-\mbox{Tr}\left[\mathcal{M}_{j}\left(\sum_{l}K_{l}\rho_{p}K_{l}^{\dagger}\right)\right]\right]^{2}+\lambda||\mathbb{K}||_{1},$
(21)

where

$d_{j}=\mbox{Tr}(\mathcal{M}_{j}\rho^{\prime}),$ (22)

$\mathcal{M}_{j}$ is a set of measurements, and
$\mathbb{K}=\left[K_{1},K_{2},\ldots,K_{n_{k}}\right]$ is a $n_{k}M\times
M^{\prime}$ matrix formed with all the $M\times M^{\prime}$ Kraus operators .
The matrix norm is given by

$||A||_{1}=max_{j}\sum_{i}|A_{ij}|,$ (23)

and $\lambda\geq 0$ is the strength of the regularization imposed on
$\mathcal{L}$. $\lambda$ is a hyperparameter of the minimization process and
can be fixed or optimized. The Kraus operators must fulfill the condition

$\sum_{l}K_{l}^{\dagger}K_{l}=\mathbb{I}.$ (24)

As we said previously, we consider two cases $\rho^{\prime}=\rho_{ABC}$ and
$\rho^{\prime}=\rho_{MD}$. In the following, we focus on the former case and
will return to the latter near the end of the Section.

###  $\rho^{\prime}=\rho_{ABC}$ case

In Reference [34], Ahmed et al. discussed all the necessary details to
implement the minimization of $\mathcal{L}$ and provided the code to reproduce
their results. To obtain the results found in this Section, we adapted the
code. For the case where $\rho^{\prime}=\rho_{ABC}$, the only adaptation
needed arises from the different dimensions of the Hibert space where $\rho$
and $\rho^{\prime}$ act.

When the dimension of both Hilbert spaces is the same, it is reasonable to use
random unitary square matrices as the initial random Kraus operators needed by
the minimization algorithm. So,

$\frac{1}{n_{k}}\sum U_{l}U_{l}^{\dagger}=1,$ (25)

results in

$K_{l}=\frac{1}{\sqrt{n_{k}}}U_{l}.$ (26)

There are different numerical methods to construct square random unitary
matrices.

Figure 7: The figure shows the typical behavior of the cost function
$\mathcal{L}$ versus the number of iterations of the gradient-descent method.
As the algorithm requires, we instantiate randomly chosen initial Kraus
operators into Equation 21 to produce the initial value of the cost function.
The data corresponds to the case where the target $\rho_{ABC}$ is the RDM of
the topological state with $k_{z}=0.1$, $L=0$, and lower energy see Figure 3.
After a few hundred iterations, the cost function value becomes smaller than
the tolerance value, and the algorithm halts. In this case, the algorithm
employs $n_{k}=20$ Kraus operators.

For the case when $\rho^{\prime}=\rho_{ABC}$ is a $4\times 4$ matrix and
$\rho_{p}$ is an $M\times M$matrix, with $M$ a multiple of four, the quantum
process $\mathcal{E}$ requires $n_{k}$ $4\times M$ Kraus operators. Each Kraus
operator is composed of $M/4$ blocks. We choose, as the initial Kraus
operators, $4\times M$ matrices such that they have only one $4\times 4$ block
different from zero, and this block is a random unitary matrix. If $K_{1}$ has
only the first block different from zero, $K_{2}$ the second one, and so on,
it is clear that the condition $Tr=1$ requires a renormalization of each
random unitary matrix by a factor equal to $\frac{1}{\sqrt{n_{k}}}$.

The measurements $M_{j}$ are given in terms of the eigenvectors of all the
operators of the form $\sigma^{\beta}\otimes\sigma^{\gamma}$, where
$\beta=x,y,z$ and also $\gamma=x,y,z$ that is, all the operators that are the
tensorial product of two Pauli matrices.

Figure 7 shows the typical behavior of the cost function as a function of the
number of iterations performed by the gradient-descent minimization method.
The algorithm that implements the gradient-descent minimization method has a
tolerance parameter. Once the cost function value becomes smaller than the
tolerance value, the algorithm does not further iterate. We used a tolerance
value of $0.01$ for $n_{k}=64$. Imposing lower values for the tolerance does
not necessarily improve the results obtained and, in some cases, leads to
oscillations in the cost function behavior. Since the cost function value does
not become zero, the predicted RDM will differ from the one employed as the
target RDM, that is,
$\tilde{\rho}_{ABC}=\tilde{\mathcal{E}}(\rho_{p})\neq\rho_{ABC}$, where
$\tilde{\mathcal{E}}$ is the quantum process found using a finite non-null
tolerance. Other sources of errors that prevent the predicted RDM from
becoming equal to the target RDM are the number of Kraus operators used,
numerical precision, etc.

Figure 8 shows a comparison between the values obtained for the
$S(\rho_{ABC})$ entropy and the entropy calculated with the RDM resulting from
the quantum process $\tilde{\rho}_{ABC}$. In panel a), the entropy
$S(\rho_{ABC})$ is shown using black solid dots, while for
$S(\tilde{\rho}_{ABC})$, we use red solid dots.

Figure 8: The von Neumann entropy for both the target reduced density matrix
$\rho_{ABC}$ and the predicted one $\tilde{\rho}_{ABC}$ calculated for several
topological states, all with $L=0$ and different wave number $k_{z}$. a) The
von Neumann entropy $S(\rho_{ABC})$ is shown using black solid dots, while
$S(\tilde{\rho}_{ABC})$ is shown using red solid dots. Both entropies show a
plateau for $k_{z}$ values consistent with the one observed in Figure 3. Note
that the topological states used correspond to the upper branch of energies.
b) This panel shows $|S(\tilde{\rho}_{ABC})-S(\rho_{ABC})$ using green solid
dots, versus $k_{z}$. We calculated the modulus of the difference between the
entropies with the data contained in panel a). The number next to each dot is
the fidelity between both sets of matrices,
$\mathcal{F}(\rho_{ABC}(k_{z},L),\tilde{\rho}_{ABC}(k_{z},L))$.

It is easy to appreciate that the quantum process consistently results in RDM
that gives larger entropy values than those corresponding to RDM calculated
with the variational eigenvectors, as in Eq. 9. Panel b) shows the modulus of
the difference between $S(\rho_{ABC})$ and $S(\tilde{\rho}_{ABC})$ using green
points. The differences shown correspond to the data shown in panel a). The
number next to each green dot is the fidelity
$\mathcal{F}(\rho_{ABC},\tilde{\rho}_{ABC})$. Despite the excellent fidelity
between both sets of reduced-density matrices, the calculated, $\rho_{ABC}$,
and the predicted by the quantum process, $\tilde{\rho}_{ABC}$, the
differences in their corresponding entropies can be as large as $0.047$ or a
relative error of 7$\%$. Nevertheless, note that the predicted entropy is also
a constant where the calculated entropy is a constant and that its value also
drops to zero for the values of $k_{z}$ where the energy of the topological
state becomes close enough to the band, see panel a).

###  $\rho^{\prime}=\rho_{MD}^{\prime}$ case

In this case, the quantum process takes $\rho_{p}$, which is $4M\times 4M$
matrix, to a predicted matrix $\rho_{MD}^{\prime}$, which is also a $4M\times
4M$ matrix. In our case, $M=64$ for the largest basis set size, resulting in
$2^{8}\times 2^{8}$ matrices.

Figure 9: The von Neumann entropy for both the target mode-dependent reduced
density matrix $\rho_{MD}$ and the predicted one $\tilde{\rho}_{MD}$
calculated for several topological states, all with $L=0$ and different wave
number $k_{z}$. a) The von Neumann entropy $S(\rho_{MD})$ is shown using black
solid dots, while $S(\tilde{\rho}_{md})$ is shown using red solid dots. Both
entropies show a plateau for $k_{z}$ values consistent with the one observed
in Figure 3. Note that the topological states used correspond to the upper
branch of energies. b) This panel shows $|S(\tilde{\rho}_{MD})-S(\rho_{MD})|$
using green solid dots, versus $k_{z}$. We calculated the modulus of the
difference between the entropies with the data contained in panel a). The
number next to each dot is the fidelity between both sets of matrices,
$\mathcal{F}(\rho_{MD}(k_{z},L),\tilde{\rho}_{MD}(k_{z},L))$.

In Reference [34], the suggested number of measurements employed in Eq. 21 is
$6^{n_{q}}$, where $n_{q}$ is the number of ”qubits” over which the quantum
process acts. Allocating $6^{8}$ vectors, each one of $2^{8}$ components, as
the method requires, becomes impractical. In our case, it is sufficient to
consider only the measurements associated with the eigenvectors of the
operator given by

$\bigotimes_{i=1}^{8}\sigma_{x}^{i}.$ (27)

The operator above has $2^{8}$ eigenvectors, a more manageable number than
$6^{8}$. This adaptation is the only one needed to run the algorithm since the
starting and predicted matrices have the same size, rendering all the Kraus
operators square matrices.

The behavior of the cost function for this case is similar to the previous one
see Figure 7. The number of steps necessary to reach the required tolerance
(0.1) is higher, typically around seven hundred, and $n_{k}=30$.

Figure 9 shows in panel a) both entropies, $S(\rho_{MD})$ and
$S(\tilde{\rho}_{MD})$, for several values of $k_{z}$, while panel b) shows
the modulus of their difference and the fidelity between both matrices. Note
that all the fidelities are better than $0.995$, and the modulus of the
differences are all lower than $0.05$, which results in relative errors of
less than $2\%$.

##  VI Discussion and conclusions

When dealing with one-particle wave functions or spinors, it is possible to
trace out over only one of the coordinates and obtain a coordinate-dependent
RDM, which leads to entropies (or information-like quantities) that detect
transitions in quantum states of one-electron systems [44, 29].

In the case of a many-particle system, there are more possibilities since it
is possible to trace over all the coordinates of a subset of particles
(entanglement between particles) [22, 45, 21] or subsets of coordinates (real-
space entanglement) [31, 46].

Our calculation of the topological entropy exploits the fact that the
eigenvectors of the $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian is a spinor with
four components, $|\psi\rangle$. Using the spinor, we obtain a pure density
matrix $|\psi\rangle\langle\psi|$, a coordinate-dependent $4\times 4$ matrix.
Tracing out a spatial region, as in Equation 9, leads to a coordinate-
independent $4\times 4$ reduced density matrix whose eigenvalues are easy to
calculate.

Interestingly, for the topological states constructed as the spinors in Eq.
11, all the entropies that enter into the calculus of the topological entropy
are constants independent of $k_{z}$ and $L$ except near $k_{z}=0$, that is,
$S(\rho_{sp})=const.$, where $sp$ is any one of the triangular sectors $A,B$
or $C$, or a combination of them. Because all these entropies are constant, it
is enough to look for a quantum process such that
$\rho_{ABC}=\mathcal{E}(\rho_{p})$ to test the method.

A caveat about the ”constant” value of the topological entropy is in order
since the statement is accurate as long as the radius of the triangular
sectors $R_{c}$ is not too close to the radius of the nanowire, $R$, or
becomes too small.

The results in Sections 3, 4, and V correspond to a nanowire 120 nanometers in
diameter and a radius $R_{c}=15$ nanometers. We also ran numerous numerical
tests for other diameters and values of $R_{c}$, and the results were
qualitatively the same.

The entropy and entanglement spectrum obtained with the mode-dependent RDM are
practical tools to distinguish topological from non-topological states, and we
intend to test them in other contexts beyond the topological states in
$\mbox{Bi}_{2}\mbox{Se}_{3}$ nanowires.

That the mode-dependent RDM could result from a quantum process makes us think
that it is possible to define operatively an equivalent RDM, whose entropy
allows us to distinguish between topological and non-topological states
without resorting to a particular basis set of functions.

###### Acknowledgements.

The authors acknowledge partial financial support from CONICET (PIP
11220210100787CO). and SECYT-UNC.

##  Appendix A Variational Method

The simpler Rayleigh-Ritz variational method requires a set of appropriate
basis functions that we denote as $f_{i}$. The expectation value of the
Hamiltonian

$\langle f|H|f\rangle,$ (28)

is obtained using a test function given by

$f=\sum_{i=1}^{N}c_{i}f_{i},$ (29)

and the minimization is performed over the values of the $c_{i}$ coefficients.
This procedure results in an algebraic problem

$\left[H\right]\kappa_{j}=E^{v}_{j}\kappa_{j},$ (30)

where $\left[H\right]$ is an $N\times N$ matrix whose entries are given by

$\left[H\right]_{ij}=\left\langle f_{i}|H|f_{j}\right\rangle,$ (31)

$\kappa_{j}$ is a vector that contains the $c$ coefficients, and the
$E^{v}_{j}$ are the variational eigenvalues. $N$ is the number of functions in
the basis set, and the algebraic problem has $N$ eigenvalues and their
corresponding eigenvectors.

For dealing with a $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian, like the one in
Equation (1), the Rayleigh-Ritz method requires some adaptations. Instead of a
test function, we need a test spinor to calculate the Hamiltonian expectation
value. In particular, we employ the following,

$|\psi_{L,k_{z}}\rangle=\sum_{n}\left(\\!\\!\\!\begin{array}[]{c}{{b_{L,k_{z},n,\uparrow}A_{L,n}J_{L}\left(\alpha_{n}^{L}\rho/R\right)e^{iL\varphi}}}\\\
{{c_{L,k_{z},n,\uparrow}A_{L,n}J_{L}\left(\alpha_{n}^{L}\rho/R\right)e^{iL\varphi}}}\\\
{{b_{L,k_{z,n,\downarrow}}A_{L+1,n}J_{L+1,n}\left(\alpha_{n}^{L+1}\rho/R\right)e^{i(L+1)\varphi}}}\\\
{{c_{L,k_{z,n,n,\downarrow}}A_{L+1,n}J_{L+1,n}\left(\alpha_{n}^{L+1}\rho/R\right)e^{i(L+1)\varphi}}}\end{array}\\!\\!\\!\right)e^{ik_{z}.z},$
(32)

where $J_{L}$ is the Bessel function with index $L$, $\alpha_{n}^{L}$ is its
$n-th$ root, $L$ is an integer number, $R$ is the radius of the cylinder,
$b_{L,k_{z,n,\uparrow}},c_{L,k_{z,n,\uparrow}},b_{L+1,k_{z,n,\downarrow}},c_{L+1,k_{z,n,\downarrow}}$
are the coefficients of the expansion (i.e., the linear variational
parameters), and

$A_{L,n}=\frac{1}{\sqrt{\pi}RJ_{L+1}(\alpha_{n}^{L})},$ (33)

is a normalization constant.

With the provisos mentioned in the paragraph above, the resulting algebraic
problem for the $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian has $4N$ eigenvalues
and their corresponding eigenvectors.

Finally, to calculate the matrix elements necessary for the Rayleigh-Ritz
method, we write the operators $k_{\parallel}$ and $k_{\pm}$ as differential
operators in cylindrical coordinates using that in Cartesian ones

$\displaystyle k_{x}=-i\partial_{x}$ (34) $\displaystyle k_{y}=-i\partial_{y}$
(35) $\displaystyle k_{z}=k_{z}.$ (36)

The change to cylindrical coordinates leads to

$\displaystyle\partial_{x}=cos(\varphi)\partial_{\rho}-\frac{sin(\varphi)}{\rho}\partial_{\varphi},$
(37)
$\displaystyle\partial_{y}=sin(\varphi)\partial_{\rho}-\frac{cos(\varphi)}{\rho}\partial_{\varphi},$
(38)
$\displaystyle\partial^{2}_{x}=cos^{2}(\varphi)\partial^{2}_{\rho}+\frac{sin^{2}(\varphi)}{\rho}\partial_{\rho}+\frac{sin(2\varphi)}{\rho^{2}}\partial_{\varphi}+\frac{sin^{2}(\varphi)}{\rho^{2}}\partial^{2}_{\varphi}$
(39)
$\displaystyle\partial^{2}_{y}=sin^{2}(\varphi)\partial^{2}_{\rho}+\frac{cos^{2}(\varphi)}{\rho}\partial_{\rho}-\frac{sin(2\varphi)}{\rho^{2}}\partial_{\varphi}+\frac{cos^{2}(\varphi)}{\rho^{2}}\partial^{2}_{\varphi},$
(40)

which results in

$\displaystyle k_{\pm}=k_{x}\pm ik_{y}=e^{\pm
i\varphi}\left(-i\left(\partial_{\rho}\pm\frac{i}{\rho}\partial_{\varphi}\right)\right),$
(41) $\displaystyle
k_{\parallel}=k_{x}^{2}+k_{y}^{2}=-\left(\partial^{2}_{\rho}+\frac{1}{\rho}\partial_{\rho}+\frac{1}{\rho^{2}}\partial^{2}_{\varphi}\right).$
(42)

## References

  * [1] Xiao-Liang Qi and Shou-Cheng Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * [2] Kane C L and Mele E J, Phys. Rev. Lett. 95, 226801 (2005).
  * [3] Bernevig B A and Zhang S-C, Phys. Rev. Lett. 96, 106802 (2006).
  * [4] Shuichi Murakami, New J. Phys. 13, 105007 (2011).
  * [5] Masatoshi Sato and Yoichi Ando, Rep. Prog. Phys. 80, 076501 (2017).
  * [6] Kezilebieke, S., Huda, M.N., Vaňo, V. et al.. Nature 588, 424–428 (2020).
  * [7] Chao-Xing Liu, Xiao-Liang Qi, HaiJun Zhang, Xi Dai, Zhong Fang, and Shou-Cheng Zhang Phys. Rev. B 82, 045122 (2010)
  * [8] Zhang, H., Liu, CX., Qi, XL. et al., Nature Phys 5, 438–442 (2009)
  * [9] Liang Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007)
  * [10] Wen-Kai Lou, Fang Cheng, Jun Li, Journal of Applied Physics 110, 093714 (2011).
  * [11] M. Governale et al , New J. Phys. 22, 063042 (2020)
  * [12] Iorio P, Perroni C A and Cataudella V, Eur. Phys. J. B 89, 97 (2016).
  * [13] Linder J, Yokoyama T and Sudbø A, Phys. Rev. B 80, 205401 (2009)
  * [14] L. Gioia, M. G. Christie, U. Zülicke, M. Governale, and A. J. Sneyd, Phys. Rev. B 100, 205417 (2019).
  * [15] A. Kitaev y J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
  * [16] Hui Li and F. D. M. Haldane, Phys. Rev. Lett. 101, 010504 (2008).
  * [17] Lukasz Fidkowski, Phys. Rev. Lett. 104, 130502 (2010).
  * [18] Xiao-Liang Qi, Hosho Katsura, and Andreas W. W. Ludwig Phys. Rev. Lett. 108, 196402 (2012).
  * [19] P. Calabrese and A. Lefevre Phys. Rev. A 78, 032329 (2008).
  * [20] Hong Yao and Xiao-Liang Qi Phys. Rev. Lett. 105, 080501 (2010).
  * [21] M. Haque, O. Zozulya, and K. Schoutens, Phys. Rev. Lett. 98, 060401 (2007).
  * [22] O. S. Zozulya, M. Haque, K. Schoutens, and E. H. Rezayi, Phys. Rev. B 76, 125310 (2007).
  * [23] V. V. Ravi Kishore, B. Partoens and F. M. Peeters, J. Phys.: Condens. Matter 24 135302 (2012).
  * [24] V. V. Ravi Kishore, B. Partoens and F. M. Peeters, J. Phys.: Condens. Matter, 26, 095501 (2014).
  * [25] S. S. Krishtopenko and F. Teppe, Phys. Rev. B 97, 165408 (2018).
  * [26] S. S. Krishtopenko, W. Knap and F. Teppe, Scientific Reports 6, 30755 (2016).
  * [27] R. Skolasinski, D. I. Pikulin, J. Alicea, and M. Wimmer, Phys. Rev. B 98, 201404(R) (2018).
  * [28] Zewei Chen and Tai Kai Ng, Phys. Rev. B 99, 235157 (2019).
  * [29] N. Giovenale and O. Osenda Physica B 627, 413564 (2022)
  * [30] N. Giovenale and O. Osenda Physica E, 144, 115406 (2022)
  * [31] A. Sterdyniak, A. Chandran, N. Regnault, B. A. Bernevig, and Parsa Bonderson, Phys. Rev. B 85, 125308 (2012).
  * [32] E. C. G. Sudarshan, P. M. Mathews, and Jayaseetha Rau, Phys. Rev. 121, 920 (1961).
  * [33] Michael A Nielsen and Isaac L. Chuang, Quantum Computation and Quantum Information (10th ed.). Cambridge: Cambridge University Press. (2010).
  * [34] S. Ahmed, F. Quijandría y A. F. Kockum, Phys. Rev. Lett. 130, 150402 (2023).
  * [35] J. Fiurášek and Z. Hradil, Phys. Rev. A 63, 020101(R) (2001).
  * [36] M. F. Sacchi, Phys. Rev. A 63, 054104 (2001).
  * [37] A. Anis and A. I. Lvovsky, New J. Phys. 14, 105021 (2012).
  * [38] K. Schultz, Phys. Rev. A 100, 062316 (2019).
  * [39] G. C. Knee, E. Bolduc, J. Leach, and E. M. Gauger, Phys. Rev. A 98, 062336 (2018).
  * [40] T. Surawy-Stepney, J. Kahn, R. Kueng, and M. Guta, Quantum 6, 844 (2022).
  * [41] C. H. Baldwin, A. Kalev, and I. H. Deutsch, Phys. Rev. A 90, 012110 (2014).
  * [42] Y. S. Teo, G. I. Struchalin, E. V. Kovlakov, D. Ahn, H. Jeong, S. S. Straupe, S. P. Kulik, G. Leuchs, and L. L. Sánchez-Soto, , Phys. Rev. A 101, 022334 (2020).
  * [43] S. Xue, Y. Liu, Y. Wang, P. Zhu, C. Guo, and J. Wu, Phys. Rev. A 105, 032427 (2022).
  * [44] M. Garagiola and O. Osenda, Physica E 116, 113755 (2020)
  * [45] O. Osenda, F. M Pont, A. Okopińska and P. Serra, J. Phys. A: Math. Theor. 48, 485301 (2015)
  * [46] I. D. Rodríguez and G. Sierra, Phys. Rev. B 80, 153303 (2009)