title: Probing the chirality of 1D Majorana edge states around a 2D nanoflake in a superconductor
authors: Andrzej Ptok David J. Alspaugh Szczepan Głodzik Aksel Kobiałka Andrzej M. Oleś Pascal Simon Przemysław Piekarz
text: 
  
  
  
  
  
  
  [e-mail: ]aptok@mmj.pl
  
  
  
  
  [e-mail: ]dalspa1@lsu.edu
  Department of Physics and Astronomy, Louisiana State University, 
  Baton Rouge, Louisiana 70803-4001, USA
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Max Planck Institute for Solid State Research,
  Heisenbergstrasse 1, D-70569 Stuttgart, Germany
  
  
  
  [e-mail: ]pascal.simon@u-psud.fr
  Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, F-91405, Orsay, France
  
  
  
  [e-mail: ]piekarz@wolf.ifj.edu.pl
  
  
  
  
  
  
  
  The interplay between superconductivity, magnetic field and spin–orbit 
  coupling can lead to the realization of non–trivial topological phases.
  Recent experiments 
  
  
  
  have found signatures of such phases in magnetic nanoflakes formed by nanostructures coupled to a superconducting substrate.
  These heterostructures comprise  a topologically non-trivial region  surrounded by a trivial one due to the finite magnetic exchange field induced by the magnetic nanoflake. 
  The analysis of the topological phase diagram of such a system shows that  a similar phase separation occurs by tuning the chemical potential of the nanoflake.
  In this paper, we study such a possibility in detail, analyzing the spatial extent of the edge modes circulating around the nanoflake and discussing some practical implementations. 
  We also show how the chirality of Majorana edge states can be probed using scanning tunneling spectroscopy with a double tip setup.
  
  
  
  Probing the chirality of 1D Majorana edge states around 
   a 2D nanoflake in a superconductor
      Przemysław Piekarz
      December 30, 2023
  ============================================================================================
  
  
  
  
  
  
  § INTRODUCTION
  
  The quest for the realization of Majorana zero
  modes (MZMs), driven by the pursuit of both fundamental physics and their potential application to fault-tolerant topological quantum computation <cit.>, is steering an active research in engineering p-wave superconductivity.
  Non-Abelian braiding is an essential step towards topological
  quantum computing, though it has not yet been experimentally achieved with MZMs. 
  Due to the localized nature of MZMs, their braiding will necessarily involve both coupling and manipulation processes.
  
  However, it has been suggested that non-Abelian braiding is not only restricted to MZMs, but can also be implemented  with one-dimensional (1D) chiral Majorana fermions <cit.>.
  Chiral Majorana fermions can manifest themselves as quasiparticle
  edge states of a two-dimensional (2D) topological p-wave superconductor <cit.>. 
  Signatures of 1D chiral Majorana quasiparticles have been recently observed in 2D heterostructures consisting of a quantum anomalous Hall insulator bar in contact with a superconductor <cit.>.
  Additionally, recent progress in atomic scale engineering <cit.>  is opening up new perspectives for the practical implementation of chiral Majorana fermions by spatially building  non–trivial topological phases separated from trivial ones.
  
  Recent progress includes Co islands grown on a Si substrate and covered by a monolayer  of Pb <cit.> and nano-scale Fe islands of  monoatomic height on a Re surface <cit.>.
  Due to the non–trivial topological phase transition resulting from a gap closure, in-gap edge states surrounding the topological superconducting (SC) domain are observed. 
  In both experiments, these in-gap states are strongly delocalized around the islands and have been interpreted as signatures of chiral Majorana fermions.
  
  
  
  
  
  In 2D superconductors with Rashba spin–orbit coupling (SOC), the transition to a non–trivial phase can be induced by an external Zeeman magnetic field  
  The boundary between the trivial and non–trivial topological phases is given by h_c^2=μ^2+Δ^2,
  see Fig. <ref>(a), where h_c stands for the critical Zeeman field for given values of the doping μ and the SC gap Δ. 
  In the aforementioned experimental results, the Zeeman magnetic energy arises from the presence of magnetic dopants interacting with the substrate, while the SOC and the SC gap are intrinsic to the subsystem. 
  Looking at the phase diagram presented in Fig. <ref>(a), one observes that a line which connects the points A and B in the (μ,h) plane could correspond to an inhomogeneous system in real space, where a non-magnetic trivial domain (point A) surrounds or borders with a topological magnetic domain (point B).
  
  
  
  In this paper, we instead choose to explore  an alternative route. 
  We consider an inhomogeneous system but in a constant magnetic Zeeman energy, which would reside on the C–D line in Fig. <ref>(a).
  The transition to the topological domain occurs due to a change in the chemical potential μ.
  Such a system could be constructed experimentally in many different ways: one option would be to substitute the magnetic atoms in the Co island <cit.> by non-magnetic ones and add a magnetic field parallel to the SC Pb monolayer [Fig. <ref>(b)].
  Another promising approach is to use the versatility offered by 2D van-der-Waals heterostructures <cit.>.
  A possible way to generate an homogeneous Zeeman exchange energy in the proximity of a superconductor could be engineered by stacking recently synthesized 2D magnetic materials <cit.> with a transition metal dichalcogenide (TMD) superconductor such as NbSe_2.
  A non-magnetic island can be obtained by evaporating some alkaline adatoms to enforce charge transfer.
  
  The paper is organized as follows:
  In Sec. <ref>, we first start with a circular geometry for the nanoflake and derive the dispersive chiral Majorana edge states analytically in the continuum limit.
  We also discuss the spatial extent of the chiral Majorana modes in the tranverse direction.
  In Sec. <ref> we compare our results obtained in the continuum limit to exact diagonalization of a tight binding model on a lattice. 
  In Sec. <ref>, we propose and study a setup in order to measure the chirality of the Majorana edge states using two ferromagnetic tips.
  Finally, we present a summary of our results and give conclusions in Sec. <ref>.
  
  
  
  
  
  
  
  
  
  § NANOFLAKE WITH CIRCULAR GEOMETRY
  
  
  
  We begin by considering a nanoflake with a circular symmetry as depicted 
  in Fig. <ref>. For simplification and without a loss 
  of generality, we can also assume a smooth boundary of the nanoflake. 
  Keeping in mind the circular symmetry, the Hamiltonian will commute with 
  the zth component of the total angular momentum operator J_z≡ L_z+S_z. 
  We may then find the energies of the bound state wave functions localized 
  at the edge of the nanoflake in terms of the m_J quantum numbers. 
  This method allows us to determine the existence of chiral subgap states 
  within our system, and has successfully been used in the studies of other 
  2D systems with circular symmetry such as graphene <cit.>.
  
  
  
  Thus, the real-space normal state Hamiltonian has the form,
  
      H(r,∇)=
      (-∇^22m-μ(r))σ_0 
      +α(σ× -i∇)_z + hσ_z ,
  
  where σ_i (for i={0,x,y,z}) are the Pauli matrices acting 
  in spin space. Here, the system is 2D with r≡(x,y), and 
  the chemical potential is given by,
  
      μ(r) = 
      μ_1    r < R_0 , 
      μ_2    r ≥ R_0 .
  
  We may also define the discontinuity δμ=μ_2-μ_1 at the boundary,
  see Fig. <ref>. 
  In other words, δμ corresponds to the spatial variation of the chemical 
  potential induced by the nanoflake. The BdG Hamiltonian may then be expressed as,
  
      ℋ = 12∫ dr Ψ^†(r)
      [     H(r,∇)      iσ_yΔ;     -iσ_yΔ -H^T(r,-∇) ]Ψ(r) ,
  
  where Ψ(r)≡ (ψ_↑(r),ψ_↓(r),
  ψ_↑^†(r),ψ_↓^†(r))^T 
  is the Nambu spinor, with ψ_σ(r) being electron field 
  operators which destroy an electron with spin σ at location r.
  
  Due to the circular symmetry of the nanoflake (or more precisely, the 
  scalar chemical potential μ(r)=μ(r)), the BdG Hamiltonian commutes 
  with the zth component of the total angular momentum operator J_z=L_z+S_z. 
  It follows that the Hamiltonian and J_z share the same eigenstates. 
  The eigenstates of J_z, with the half-integer eigenvalues m_J, are given by
  
      φ_m_J = 
      [ u_m_J↑(r)   e^i(m_J - 1/2)θ;   u_m_J↓(r) e^i(m_J + 1/2)θ; v_m_J↑(r)   e^i(m_J + 1/2)θ;   v_m_J↓(r) e^i(m_J - 1/2)θ ].
  
  To focus on states with small total angular momenta, we take a 
  low-energy approximation and neglect the kinetic energy term in the 
  Hamiltonian <cit.>. By writing the u_m_Jσ(r) 
  and v_m_Jσ(r) functions in terms of the Modified Bessel 
  Functions of the First and Second Kind, we may then solve for the bound state 
  wave functions localized at r=R_0. The details of this approach can be found 
  in the Supplemental Material (SM) [See Supplemental Material at [URL will be 
  inserted by the publisher] for details of the analytical calculations and additional 
  numerical results for the nanoflake with irregular shape. In the second case, 
  the numerical results for different sets of values λ and V_0 are presented.], 
  and the resulting wave functions have the form:
  
      φ_m_J = 
      φ_m_J1    r < R_0 , 
      φ_m_J2    r ≥ R_0 ,
  
  where
  
  
      φ_m_J1   =   ∑_η=±N_m_J1η[ a_m_J1η I_m_J - 12(k_m_J1ηr) e^i(m_J - 12)θ; b_m_J1η I_m_J + 12(k_m_J1ηr) e^i(m_J + 12)θ; c_m_J1η I_m_J + 12(k_m_J1ηr) e^i(m_J + 12)θ;         I_m_J - 12(k_m_J1ηr) e^i(m_J - 12)θ ],
  
  and
  
      φ_m_J2   =   ∑_η=±N_m_J2η[ a_m_J2η e^i(m_J-12)π  K_m_J - 12(k_m_J2ηr) e^i(m_J - 12)θ; b_m_J2η e^i(m_J+12)π  K_m_J + 12(k_m_J2ηr) e^i(m_J + 12)θ; c_m_J2η e^i(m_J+12)π  K_m_J + 12(k_m_J2ηr) e^i(m_J + 12)θ;         e^i(m_J-12)π  K_m_J - 12(k_m_J2ηr) e^i(m_J - 12)θ ].
  
  Here I(z) and K(z) are the Modified Bessel Functions of the first and 
  second kind, respectively. The a, b, and c parameters along with the 
  normalizations are derived within the SM <cit.>, while the radial 
  momenta which control the spatial extent of the bound states wave functions are given by
  
      k_m_Jjη   =   1α√( h^2 - E_m_J^2 + Δ^2 - μ_j^2 + 2η√(E_m_J^2μ_j^2 + Δ^2(h - μ_j)(h + μ_j))),
  
  
  where j = 1,2 and η=±. For each of the φ_m_J bound states 
  localized at r = R_0, the total radial spatial extent ξ_m_J of the wave 
  functions are thus determined by 
  
      ξ_m_J=max{k_m_J1+^-1,k_m_J1-^-1,k_m_J2+^-1,k_m_J2-^-1}.
  
  
  
  
  
  
  The energy spectrum of these bound states vs the m_J quantum numbers is 
  presented in Fig. <ref>, while the spatial profile of the wave functions 
  is given in Fig. <ref> for two different characteristic sets of 
  parameters. The number of in-gap states is quantized due to the finite perimeter of the nanoflake and their energy spacing depends on intrinsic parameters. In the first set of parameters, we find 8 in-gap states while we 
  have 4 in-gap states in the second set.
  We choose the strength of the Zeeman field such that the topological gap, 
  
      Δ_ top=| h - √(Δ^2+μ^2)|,
  
  is equal both inside and outside the nanoflake. From our first set of data 
  we obtain Δ_ top=0.070 meV, while for our second set we obtain 
  Δ_ top=0.0296 meV. The existence of only one branch of states 
  in the topological gap signals that they are chiral (here left handed). 
  This in turn supports the hypothesis of a phase separation in space with a chiral Majorana edge state circulating around the nanoflake.
  
  We also note that these dispersive states are not perfectly linear, but instead 
  exhibit a slight cubic behavior. Comparing both set of parameters, we see 
  that the larger the discontinuity |δμ|=|μ_1-μ_2|, the more 
  localized the edge states are. This behavior is expected; indeed, the in-gap 
  Majorana states are pinned at the domain wall which is controlled by the 
  spatial variation of the chemical potential.
  
  
  
  
  
  
  
  
  § TIGHT BINDING FORMULATION
  
  
  
  In order to confirm our calculations performed in the continuum limit and to go 
  beyond the circular symmetry assumption for a nanoflake, we have also performed
  numerical calculations on a lattice based on a tight binding description.
  
  Our system can be described by the following tight binding Hamiltonian,
  
      ℋ=ℋ_ kin+ℋ_ SO
      +ℋ_ prox+ℋ_ flake.
  
  The first term corresponds to a free particle on a 2D square lattice,
  
      ℋ_ kin = ∑_ijσ [-t+(4t-μ-σ h)δ_ij] c_iσ^†c_jσ^.
  
  Here t is the hopping integral between nearest neighbor sites
  [In the case of an homogeneous 2D system, this gives a 
  dispersion relation ℰ_ k=-2t(cos(k_x)+cos(k_y)). 
  Then, the minimal available energy is equal to -4t.], 
  μ is the chemical potential (calculated from the bottom of the band), 
  while h can be regarded either as a genuine Zeeman energy or as a magnetic 
  exchange energy depending on the situation under consideration. 
  In all cases, we treat it as an effective magnetic field in what follows.
  The second term describes the in-plane SOC,
  
      ℋ_ SO= 
      -iα∑_ijσσ' c_i+ d_jσ^†[(
       d_j×σ̂)·ẑ]_σσ'^
      c_iσ'^,
  
  where vectors d_i∈{±x̂,±ŷ} stand for the 
  locations of the neighbors of the i-th site, while 
  σ̂=(σ_x,σ_y,σ_z) is the vector with the 
  Pauli matrices being its components. A SC gap can be induced in the 
  layer through the proximity effect—this process is described through 
  the third term by the BCS-like form,
  
      ℋ_ prox=Δ∑_i(
      c_i↑^†c_i↓^†+ H.c.).
  
  
  
  
  
  
  The last term in Eq. (<ref>) denotes the influence of the nanoflake 
  at the particle distribution,
  
      ℋ_flake=-∑_i V_i^ c_iσ^† c_iσ^.
  
  We assume that every atom comprising the nanoflake changes the energy levels 
  of the rest of the sites, i.e., we assume a long-range impurity potential 
  given by V_i=V_0∑_mexp(-R_mi/λ), where the summation 
  is carried out over all adatoms in a given configuration 
  𝒱 <cit.>. Here, 
  λ denotes the characteristic length of decay of the impurity potential.
  
  The Hamiltonian ℋ can be diagonalized by the unitary transformation,
  
      c_iσ=∑_n (u_inσ^γ_n^ 
      -σ v_inσ^∗γ_n^†) ,
  
  which leads to the Bogoliubov–de Gennes (BdG) equations <cit.> 
  of the form
  
      ℰ_nΦ_in=∑_jℍ_ijΦ_jn,
  
  with eigenvectors Φ_in=(u_in↑,v_in↓,
  u_in↓,v_in↑)^T. Here
  
      ℍ_ij = ( [    H_ij↑     D_ij  S_ij^↑↓        0;   D_ij^∗ -H_ij↓^∗        0  S_ij^↓↑;  S_ij^↓↑        0    H_ij↓     D_ij;        0  S_ij^↑↓   D_ij^∗ -H_ij↑^∗ ])
  
  is the Hamiltonian in the matrix form, with matrix elements 
  H_ijσ= -t∑_jδ_⟨ i,j⟩ 
  + ( 4t -μ -σ h -V_i)δ_ij as the kinetic term, 
  D_ij=Δδ_ij describing the SC correlations, and 
  S_ij^σσ'=-iα∑_j[ 
  ( d_j×σ̂)·ẑ]_σσ'δ_⟨ i,j⟩ standing for the matrix representation 
  of the spin–orbit coupling. 
  More details of this method can be found e.g. in Ref. <cit.>.
  
  
  
  
   §.§ Numerical results
  
  
  
  We report the calculations on a N_x× N_y=59× 59 square 
  lattice with periodic boundary conditions (PBC). Omitting generality, 
  in the calculations presented in this subsection we use a nanoflake 
  with a circular shape and a radius R_0=15.1, which covers about 
  20% of the total area. We present results for a nanoflake characterized 
  by λ=1 and V_0/t=-0.06. In subsequent calculations, we take 
  α/t=0.15, Δ/t=0.3, and μ/t=0.4. 
  If not mentioned explicitly in the text, the value of the magnetic Zeeman 
  field h was chosen so that the boundary between the trivial and non–trivial 
  phases remains in the center of the artificial domain wall (typically h/t≃ 0.4).
  
  
  
  
  
  
    §.§.§ Density of states
  
  
  
  From the solutions of the BdG equations, we first calculate the local density of 
  states (LDOS) <cit.>,
  
      ρ_i(ω)   =   ∑_σ n[|u_inσ|^2δ(ω-ℰ_n) 
      +|v_inσ|^2δ(ω+ℰ_n)],
  
  where we replace the Dirac function δ(ω) by a Lorentzian, 
  δ(ω)=ζ/[π(ω^2+ζ^2)], 
  with a small broadening 
  
  Figure <ref> shows an example of the LDOS for a chosen path 
  (along the nanoflake with y=30). The non–trivial domain is separated 
  from the trivial phase by in-gap states strongly localized along the 
  edge of the nanoflake (cf. Fig. <ref>). 
  The domain wall is visible in the form of two sets of LDOS peaks with oscillating intensity near the ω/t ≃ 0, which are a result of the discrete nature of the in-gap state's spectrum.
  The LDOS in our system does not exhibit an “X”–shape crossing through 
  the energy gap around the nanoflake edge, in contrast to the results 
  presented in Refs. <cit.>.
  Nevertheless, our results are in agreement with experimental results 
  presented in Ref. <cit.>. In consequence, we 
  do not observe the two-ring localization like in Refs.  <cit.> but a 
  one-ring structure along the nanoflake boundary like in 
  Ref. <cit.>.
  
  
  
  
  
  The localization of the in-gap states is shown Fig. <ref>.
  As we can see, the edge states are localized around the nanoflake at the border between the trivial and topological phases.
  All states create localized states with a circular shape of the same radius (equal approximately to 17), while the radius of the nanoflake is equal R_0≃ 15.1 (shown by doted-dashed line).
  The difference of these results with respect  to the ones discussed in the continuum limit is a consequence of the smearing of the domain wall given by V_i. 
  
  
  
  
  
  The total density of states (DOS) of the system per site is given by a
  summation of the LDOS over the whole 2D space, i.e., 
  ρ(ω)= 1/N ∑_iρ_i(ω).
  Here, we can separate the sum into three different terms
  
      ∑_i→∑_i∈Nanoflake+∑_i∈DW+∑_i∈Bulk,
  
  a contribution from the sites belonging to nanoflake, a contribution from the domain wall (DW), and a contribution from the 
  bulk states  [cf. inset in Fig. <ref>(b)]. 
  As in Ref. <cit.>, we can define the functions 𝒞_i
  in order to classify the states in real space.
  These functions are equal to 1 when a site i belongs to a given region, and 0 otherwise.
  We assume that the NF (bulk) region is located in sites where the nanoflake 
  changes (does not change) the chemical potential significantly, i.e.,
  | V_i | < 0.95 max{| V_i |} (| V_i | > 0.05 max{| V_i |}).
  Otherwise, we treat the site as a part of the DW region, i.e., when 
  0.95max{|V_i|}≤ |V_i|≤ 0.05max{|V_i|}. 
  These conditions can be smoothed arbitrarily without changing the results qualitatively (see Fig. S6 and Fig. S7 in the SM <cit.>).
  
  
  We use this recipe to present the partial density of states (PDOS), 
  which is defined as follows
  
      ρ̃(ω)=∑_i𝒞_i ρ_i(ω).
  
  The contribution of the bulk (NF) presented in Fig. <ref>(a) [Fig. <ref>(b)], looks like the familiar DOS of a pure 2D 
  Rashba spin-orbit coupled superconductor, in which the gap closing is 
  followed by a reopening of the topological gap. Dashed blue lines 
  serve as a guide to the eye, showing the linear continuation of the 
  last negative and first positive eigenvalue.
  The value of the magnetic field h for which these dashed lines cross  zero 
  energy indicate the phase transition from trivial 
  to non–trivial phase.
  
  
  For effective fields larger than this crossing point, in both cases we 
  observe a topological gap opening around 
  ω/t≃± 0.1. Contrary to this, the DW PDOS  
  [Fig. <ref>(c)] ressembles a non–trivial 2D system with egdes 
  <cit.>. 
  Between critical fields for bulk and NF 
  regions (dashed blue lines crossing the Fermi level), we observe the 
  in-gap states associated only with DW, which confirms the existence of
  bound states localized along the DW. Additionally, with the increase of 
  h, the contribution of these states to the total DOS is shifted from the DW 
  to bulk region (cf. Fig. S8 in the SM <cit.>).
  
  
  
  
  
  
  
    §.§.§ Non–trivial topological domains and topological phase diagram
  
  
  
  
  A magnetic field, h, leads to a closing 
  of the trivial SC gap and reopening of a new non–trivial topological 
  gap. This occurs at the critical energy h_c^2=μ^2+Δ^2
   <cit.> 
  (h_c/t≃ 0.5 for our choice of parameters). In our system, the 
  value of the chemical potential varies from site to site i.e., 
  μ_i=μ+V_i. This non-homogeneity can lead to a situation in which 
  the above condition is met only locally 
  <cit.>. We can therefore 
  construct a space dependent indicator <cit.> that 
  describes the spatial distribution of the non–trivial topological 
  phase,
  
      χ_i=√( ( μ + V_i ) ^2 + Δ^2) - h.
  
  Thus, positive (negative) sign of χ_i indicates the topologically 
  trivial (non–trivial) phase.
  Indeed, from the analysis of χ_i under an increase of the  magnetic field h (see Fig. S9 in the SM <cit.>), we  find that the non-trivial phase exists in the system when χ_i < 0. 
  
  
  
  From the above analysis, it follows that the spatial dependence of the 
  χ_i indicator gives a correct information about the emergence of the 
  non–trivial topological domain inside the nanoflake. 
  
  Using this condition, we  construct a topological phase diagram in the two-parameter space defined by the chemical potential 
  μ and the effective magnetic field h shown in
  Fig. <ref>(a).
  
  
  First, we recall that the trivial and non–trivial topological phases are 
  separated by a parabolic boundary  in the 
  homogeneous system. In our system, the nanoflake introduces a 
  non-homogeneity in μ to the system, thus the boundary splits due to the 
  existence of two regions in space where 
  topological phases can 
  emerge. As a consequence, the boundary in the phase diagram evolves 
  into a stripe, whose width is given by max |V_i|. Thus, in the limit 
  V_0→ 0 (homogeneous system without any nanoflake), the “stripe” would narrow 
  down into a line h_c(μ), as seen in Fig. <ref>.
  
  
  Second, the region of the stripe separating the phases coincides with 
  the value of the gap δ E, calculated as the energy difference between 
  the eigenvalues closest to Fermi level [Fig. <ref>(b)].
  The existence of a non–trivial domain in the system is a result of 
  the occurrence of in-gap states with exponentially small eigenenergies 
  (described by δ E).
  
  
  
  
  
  
  
  
    §.§.§ Bond current
  
  
  
  
  When the system has boundaries 
  <cit.>, or if artificial 
  barriers are introduced <cit.>, in-gap states localize on the edges of the system. The in-gap states 
  are localized in a collection of preferable locations, independently of the broadening of domain wall.
  These well–localized in-gap states provide a contribution to the bond current  which can be expressed as the local charge flow and  is obtained from the Heisenberg equation 
  <cit.>,
  
      ħ∂⟨ n_i⟩/∂ t = 
      ⟨ [ ℋ,n_i ] ⟩.
  
  The current vector field can be represented as a sum of the 
  spin-dependent currents, I_i=∑_iσ I_iσ, where 
   can be 
  expressed by the BdG 
  eigenvectors <cit.>.
  
  
  
  In Fig. <ref>, we present the real space map of the bond  
  current. Color of the arrows denotes their magnitude ∝ |I_i|.
  Once again, the width of the domain wall (controlled by λ) is 
  reflected by an observable, however, this time it is through the area 
  in which there is a significant flow of charge 
  (cf. Fig. S3 in the SM <cit.>). 
  These low energy modes bear a close resemblance to the surface states 
  of three-dimensional (3D) topological insulators (TIs)
  <cit.>. 
  Due to the SOC induced band inversion and 
  bulk-boundary correspondence, the 2D surface of a TI hosts metallic 
  states which disperse through the band gap.
  Generally, band inversion 
  is a result of the non–trivial phase transition, but here
  “surface states” are limited to the boundary of the nanoflake, which 
  serves as an edge of the system. The current vector field presented in 
  Fig. <ref> is the sum of the contributions of the 
  spin-↑ and spin-↓ current. 
  
  Due to the breaking of time reversal symmetry, the spin-↑ 
  component dominates and thus resembles the “quasi-helical” situation 
  described in Ref. <cit.>. 
  It is worth noting that a 
  magnetic impurity or a ferromagnetic island proximity–coupled to a SOC 
  superconductor will exhibit a finite spin polarization, thus giving 
  rise to persistent currents 
  <cit.> as a 
  result of the magnetoelectric effect. This type of bond current along the 
  nanoflake can be observed experimentally, e.g. in the differential 
  conductance  
  
  
  
  
   §.§ Numerical results for irregular nanoflake
  
  
  
  
  In order to go beyond the circular 
  limit discussed in previous paragraphs, we locate the substituted atoms at random sites of the lattice as nearest neighbors of an initial atom, located at the center of the surface.
  
  This method results in a nanoflake with rugged boundary (cf. Fig. <ref>) between the substituted (blue) and substrate (gray) atoms, which is similar to experimental setups <cit.>.
  However, artificial construction of nanoflakes could be characterized by a more regular shape too.
  In practice, in the SM <cit.>, we show that the main properties of the system do not depend qualitatively on the shape of the nanoflake. Therefore, all properties described before remain when the flake becomes irregular.
  
  
  
  
  
  
  
  
  
  
  § PROPOSAL TO EXPERIMENTALLY MEASURE THE CHIRALITY USING SCANNING TUNNELING MICROSCOPY 
  
  
  
  
  
  
  
  
  
  Our experimental proposal is based on the double-tip measurement technique <cit.>.
  The in-gap edge states localized around the nanoflake can give a nonlocal response between two spatially separated tips, which is schematically shown in Fig. <ref>.
  Similar transconductance technique was successfully used to measure an in-gap surface band <cit.>.
  
  
  We performed the calculation of the local and non-local differential conductance using the Kwant <cit.> code to numerically obtain the scattering matrix <cit.>:
  
      S = ( [ S_11 S_12; S_21 S_22 ]) ,    S_ij = ( [ S_ij^ee S_ij^eh; S_ij^he S_ij^hh ]) .
  
  Here, S_ij^αβ is the block of the scattering amplitudes of incident particles of type β in tip j to particles of type α in tip i.
  Then, the differential conductance matrix is given as <cit.>:
  
      G_ij ( E ) ≡∂ I_i/∂ V_j = e^2/h( T_ij^ee - T_ij^he - δ_ij N_i^e) ,
  
  where I_i is the current entering terminal i from the scattering region and V_j is the voltage applied to the terminal j, and N_i^e is the number of electron modes at energy E in terminal i.
  Finally, the energy transmission is:
  
      T_ij^αβ = ( [ S_ij^αβ]^† S_ij^αβ) .
  
  
  
  
  The nonlocal response is constituted by two processes:
  
  (i) a
   direct electron transfer between the leads and (ii) the crossed Andreev reflection (CAR) of an electron from one tip into a hole in the second tip <cit.>. 
  In typical cases, the CAR contribution dominates the electron transfer <cit.> and such processes is responsible for the  Cooper pair splitter <cit.>.
  Here, we show that this technique has a potential application in  measuring the chirality of the edge state.
  In order to achieve this goal, we suggest the use of two ferromagnetic tips described by the Hamiltonian:
  
      ℋ_tip^i = ∑_ kσ( ε_ kσ - σ M_i) c_ kσ^† c_ kσ ,
  
  where ε_ kσ is the dispersion relation of the free electrons in the tips, while M_i is the magnetization of the i^ th tip.
  We assume that the tips have opposite magnetization, i.e. M_1 = - M_2.
  The tips are separated from the plane of the system by the barrier potential.
  We assume a system like that shown schematically in Fig. <ref>, i.e. tips are located exactly above the nanoflake edge, in a non-symmetric position.
  As a result, a distance between 1st and 2nd tip differs when measured clockwise or counterclockwise (along the edge state channel at the border of the nanoflake).
  
  
  
  In our system, we can find the local (G_11 and G_22) as well as non-local (G_12 and G_21) differential conductance (Fig. <ref> and Fig. <ref>, respectively).
  The local conductance G_ii can be treated as a probe of the existence of states in the system <cit.>. 
  In this sense, each state gives a positive signal in G_ii, independently of the direction of the magnetic field h (see top and bottom panels in Fig. <ref>).
  As we can see, in our case we observe several in-gap states.
  Due to the use of ferromagnetic tips we observe non-equality of G_11 and G_22 (blue and red lines, respectively).
  
  
  
  
  
  
  
  
  
   
  
  The non-local conductances G_12 and G_21 (Fig. <ref>) describe  different situations.
  In order for CAR processes to occur, electrons from both tips need to have opposite spins to create a Cooper pair and simultaneously eject a hole  from the other tip.
  Due to the applied external magnetic field the incident electrons from the edge bond current would have the same spin as the electrons from the first tip they encounter. 
  The electron with opposite spin  should come the other tip (which has opposite magnetization) and emit a hole with the same spin  to constitute a CAR process.
  Thanks to this nonlocal phenomena, we can easily determine the direction of the edge state propagation as the nonlocal conductance coincides with the direction of chiral bond current.
  Thus, if the sign of magnetic field  changes,  the direction of the edge state propagation reverts too, G_12 becomes G_21 due to the spatial symmetry of the system.
  However, tips are not perfectly magnetized, therefore nonlocal conductance measurement in the direction opposite to the chiral bond current remains nonzero (e.g. G_21 for ⊙ direction of magnetic field).
  As the edge state bond current is constituted of both spins, particles with spins not aligned with the tip magnetization scatter off the tip (cf. wavefunction localizations in Fig. <ref>).
  If the magnetic field is not present in the system, in-gap nonlocal conductance vanishes 
  (not shown), and peaks at the edge of superconducting gap appear <cit.>.
  Differences between the absolute value of non-local conductances are the consequence of the non-symmetric position of the tips.
  The change from negative to positive slope near the gap can be interpreted as a crossover from subgap transport dominated by crossed Andreev reflection to charge imbalance above the gap <cit.>.
  Additionally, the non-local conductance strongly depends on a distance between the tips <cit.>.
  
  
  To explain the above results, we analyzed the propagating modes in the system (Fig. <ref>).
  The nonlocal transport corresponds to situation where incident electron from one tip is transmitted through the edge modes as the chiral mode and is scattered into the second tip. 
  Such propagating modes can be represented in form of wavefunction ψ( r) = ( ũ_↓ , ũ_↑ , ṽ_↓ , ṽ_↑ )^T, where ũ and ṽ correspond to its electron- and hole-like component, respectively. 
  Here,  ψ is composed mostly of ↓-electron and ↑-hole component <cit.>.
  Regions with the non–zero  probability of localization correspond to the particle remnants of scattering processes, which distribution coincides with the direction of edge state propagation.
  As we mentioned above, with relatively small voltage bias, the non-local transport is dominated by CAR processes (which was envisioned by the antiparallel magnetization of G_12 and G_21).
  
  In the case of the “positive” ⊙ magnetic field, we observe propagation from 1st to 2nd tip in clockwise direction [Fig. <ref>, panels (a) and (b)], for both electron and hole  component of the wavefunction.
  Then, if we check the mode propagation from 2nd to 1st tip [Fig. <ref>, panels (c) and (d)], we can see that the chirality of the propagating modes is preserved and clockwise.
  
  In the case of the magnetic field with opposite direction, we observe modes propagating in opposite direction as previously mentioned (bottom panels in Fig. <ref>).
  Without any surprise, mode propagation direction in this case is also preserved.
  
  
  Summarizing, the non-local conductance is not just a fingerprint of the existence  of the chiral mode <cit.>, but also a tool to measure its chirality.
  
  
  
  
  
  
  
  
  
  
  
  
  § SUMMARY AND CONCLUSIONS
  
  
  
  Recent experimental results have presented the possibility of the emergence of non–trivial topological phases in magnetic nanostructures coupled to superconducting substrates <cit.>.
  
  In this paper, we have explored the artificial implementation of topological phase transitions induced by the local modification of the chemical potential.
  In this respect, we performed analytic calculations valid in the continuum limit in the case of a nanoflake with a circular geometry and found the spectrum of the system as a function of the total angular momentum. We have also studied how the transverse spatial extent of the wave function of the chiral Majorana state localized around the nanoflake depends upon the system parameters. 
  We then performed similar calculations for a finite size geometry using a  tight binding formulation.
  In-gap states correspond to prominent peaks in the LDOS only in a distinct region of space, identified as the domain wall, which should be observed relatively simply through STM experiments. 
  We have introduced a real space indicator, which locally characterizes  the topological phase. 
  
  Indeed, for a few sets of parameters, results obtained from this indicator were in agreement with those obtained in the continuum limit.
  In the case analyzed here, the effective magnetic field leads to the realization of a 
  non–trivial phase only in distinct regions of the system, creating a non–trivial superconducting dome surrounded by a trivial superconducting phase. 
  This phase separation could be observed through the measurement of a bond chiral current, which is connected to the existence of strongly localized in-gap states.
  We have shown that this current circulates around the nanoflake.
  Additionally, with the help of the indicator, we have found the topological phase diagram of the system. 
  We have shown that an artificial phase separation can be induced for a finite range of effective magnetic field. 
  The boundary of this phase separation is strongly related to the effective magnetic field, which determines the transition between
  trivial to non–trivial phases 
  
  
  In the last part of our work, we have proposed an experimental method to measure the chiraliy of the edge states, based on a double-tip measurement of the non-local differential conductivity. We have found that
  the non-local transport properties between the two tips allow one to determine  the chirality of the edge state. Although challenging, this type of experiment should be accessible with the present technology.
  
  
  
  
  
  We would like to thank Tristan Cren for interesting discussions.
  This work was supported by the National Science Centre (NCN, Poland) 
  under grants Nos.:
  2017/24/C/ST3/00276 (A.P.),   
  2017/27/N/ST3/01762 (S.G.),   
  2018/31/N/ST3/01746 (A.K.),   
  2016/23/B/ST3/00839 (A.M.O.), 
  and
  2017/25/B/ST3/02586 (P.P.).   
  A. M. Oleś is grateful for the Alexander von Humboldt
  Foundation Fellowship .
  AP appreciate also founding in the frame of scholarships of the Minister of Science and Higher Education (Poland) for outstanding young scientists (2019 edition, no. 818/STYP/14/2019).
  
  
  
  
  
  
  
  
  
  