# Prabhupadavani: A Code-mixed Speech Translation Data for 25 Languages

Jivnesh Sandhan1, Ayush Daksh2, Om Adideva Paranjay3,
Laxmidhar Behera1,4 and Pawan Goyal2
1IIT Kanpur, 2IIT Kharagpur, 3University of Pennsylvania, 4IIT Mandi
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Nowadays, the interest in code-mixing has become ubiquitous in Natural
Language Processing (NLP); however, not much attention has been given to
address this phenomenon for Speech Translation (ST) task. This can be solely
attributed to the lack of code-mixed ST task labelled data. Thus, we introduce
Prabhupadavani, which is a multilingual code-mixed ST dataset for 25
languages. It is multi-domain, covers ten language families, containing 94
hours of speech by 130+ speakers, manually aligned with corresponding text in
the target language. The Prabhupadavani is about Vedic culture and heritage
from Indic literature, where code-switching in the case of quotation from
literature is important in the context of humanities teaching. To the best of
our knowledge, Prabhupadvani is the first multi-lingual code-mixed ST dataset
available in the ST literature. This data also can be used for a code-mixed
machine translation task. All the dataset can be accessed at:
https://github.com/frozentoad9/CMST.

##  1 Introduction

Speech Translation (ST) is a task in which speech is simultaneously translated
from source language to a different target language.111We refer to speech
translation as a speech-to-text task. It aids to overcome the language
barriers across different communities for various applications such as social
media, education, tourism, medical etc. Earlier attempts to build a robust
speech translation system mainly focused on a cascaded approach where two
separate architectures for Automatic Speech Recognition (ASR) and Machine
Translation (MT) are used in pipeline mode Cho et al. (2013); Post et al.
(2013); Tsvetkov et al. (2014); Ruiz et al. (2015); Sperber et al. (2017).
However, these approaches mainly suffer from cascading effect of error
propagation. Thus, attention shifted to end-to-end approaches Bérard et al.
(2018); Duong et al. (2016); Weiss et al. (2017); Bansal et al. (2017) due to
their ability to obliviate error propagation and ease of maintaining a single
architecture. However, these end-to-end approaches could not match the
performance of cascaded systems due to the lack of sufficiently large data
Niehues et al. (2021). Notably, with the recent upsurge in ST datasets Di
Gangi et al. (2019); Zanon Boito* et al. (2020); Iranzo-Sánchez et al. (2020);
Wang et al. (2020), this gap has been closed Niehues et al. (2021); Ansari et
al. (2020).

Nowadays, most users prefer to communicate using a mixture of two or many
languages on platforms such as social media, online blogs, chatbots, etc.
Thus, code-mixing has become ubiquitous in all kinds of Natural Language
Processing (NLP) resources/tasks Khanuja et al. (2020); Chakravarthi et al.
(2020); Singh et al. (2018a, b); Dhar et al. (2018). However, the existing NLP
tools may not be robust enough to address this phenomenon of code-mixing for
various downstream NLP applications Srivastava and Singh (2021). Therefore,
there has been a surge in creating code-mixed datasets: (1) to understand
reasons for the failure of existing models, and (2) to empower existing models
for overcoming this phenomenon. Nevertheless, it is challenging to find
natural resources that essentially capture different aspects of code-mixing
for creating datasets for a wide range of NLP tasks. Although there has been
considerable research in generating copious data and novel architectures for
the ST task, we find that not much attention has been given to address the
code-mixing phenomenon on the ST task. Possibly, this can be justified due to
the lack of a code-mixed ST dataset. To the best of our knowledge, no such
sufficiently large, multi-lingual, naturally occurring, code-mixed dataset is
available for the ST.

Thus, in this work, we introduce Prabhupadavani, a multi-lingual, multi-
domain, speech translation dataset for 25 languages containing 94 hours of
speech by 130 speakers. The Prabhupadavani is about Vedic culture and heritage
from Indic literature, where code-switching in the case of quotation from
literature is important in the context of humanities teaching. The multiple
domains cover utterances from public lectures, conversations, debates, and
interviews on various social issues. This is the first code-mixed data for
speech translation to the best of our knowledge. It is code-mixed with
English, Bengali and Sanskrit. From the typological point of view, the
languages covered vary over ten language families. All the audios files have
been manually aligned and translated. We believe that our work will ignite
research in this direction to understand- (1) How to make existing systems
robust for handling this phenomenon effectively? (2) Can multi-lingual
training help to improve performance on code-mixed speech translation? (3)
Will the gap between the cascade and end-to-end systems be closed? (4) Can we
train a single model for all languages using parallel nature of the
dataset?222Prabhupadavani has parallel translations available in all the 25
languages for all the utterances.

Languages | # Types | # Tokens | Types per line | Tokens per line | Avg. token length
---|---|---|---|---|---
English | 40,324 | 601,889 | 10.58 | 11.27 | 4.92
French (France) | 50,510 | 645,651 | 11.38 | 12.09 | 5.08
German | 50,748 | 584,575 | 10.44 | 10.95 | 5.57
Gujarati | 41,959 | 584,989 | 10.37 | 10.95 | 4.46
Hindi | 29,744 | 716,800 | 12.36 | 13.42 | 3.74
Hungarian | 84,872 | 506,608 | 9.13 | 9.49 | 5.89
Indonesian | 39,365 | 653,374 | 11.54 | 12.23 | 6.14
Italian | 52,372 | 512,061 | 9.23 | 9.59 | 5.37
Latvian | 70,040 | 477,106 | 8.69 | 8.93 | 5.72
Lithuanian | 75,222 | 491,558 | 8.92 | 9.20 | 6.04
Nepali | 52,630 | 570,268 | 10.03 | 10.68 | 4.88
Persian (Farsi) | 51,722 | 598,096 | 10.61 | 11.20 | 4.10
Polish | 71,662 | 494,263 | 8.99 | 9.25 | 5.86
Portuguese(Brazil) | 50,087 | 608,432 | 10.80 | 11.39 | 5.12
Russian | 72,162 | 490,908 | 8.96 | 9.19 | 5.79
Slovak | 73,789 | 520,465 | 9.39 | 9.75 | 5.37
Slovenian | 68,619 | 516,649 | 9.35 | 9.67 | 5.30
Spanish | 49,806 | 608,868 | 10.75 | 11.40 | 5.07
Swedish | 48,233 | 581,751 | 10.31 | 10.89 | 5.00
Tamil | 84,183 | 460,678 | 8.37 | 8.63 | 7.65
Telugu | 72,006 | 464,665 | 8.34 | 8.70 | 6.56
Turkish | 78,957 | 453,521 | 8.27 | 8.49 | 6.35
Bulgarian | 60,712 | 564,150 | 10.10 | 10.56 | 5.24
Croatian | 73,075 | 531,326 | 9.58 | 9.95 | 5.28
Danish | 50,170 | 587,253 | 10.40 | 11.00 | 4.98
Dutch | 42,716 | 595,464 | 10.52 | 11.15 | 5.05

Table 1: Statistics of the Prabhupadavani dataset

##  2 Related Work

### Speech Translation:

Recently, there have been increased efforts for creating large speech
translation data sets for many languages. Nevertheless, the available datasets
are limited to certain languages or only underpaid licenses for non-English
languages. Most of the relatively larger datasets are English-centric Di Gangi
et al. (2019); Bérard et al. (2018), domain-specific Zanon Boito* et al.
(2020); Iranzo-Sánchez et al. (2020) or with limited speech hours Zanon Boito*
et al. (2020); Zhang et al. (2021). Recently, there has been upsurge in ST
datasets in the literature Di Gangi et al. (2019); Iranzo-Sánchez et al.
(2020); Wang et al. (2020); Salesky et al. (2021). To the best of our
knowledge, there is no such naturally occurring, sufficiently large and multi-
lingual ST dataset that contains a code-switching phenomenon. We fill this gap
by contributing a code-mixed speech translation dataset for 25 languages.

### Code-mixing:

Code-mixing is ubiquitous and well addressed on variety of downstream NLP
tasks. However, majority of code-mixed datasets are synthetically generated
Gonen and Goldberg (2018); Khanuja et al. (2020); therefore, they may not be
able to capture the different aspects of code-mixing. The code-mixed datasets
for an ASR task are either limited to only one/two languages or contain only a
few hours of speech data Nakayama et al. (2019); Lyu et al. (2015). Synthetic
code-mixed datasets may not capture the different aspects of code-mixing. To
the best of our knowledge, Prabhupadavani is the first code-mixed dataset
available for 25 languages on speech translation task.

##  3 Data Description

### Resource:

Vanimedia’s Multi-language Subtitle
Project333https://vanimedia.org/wiki/Multi-language_Subtitle_Project has
created 1,080 audio mini-clips of Śrīla Prabhupāda’s lectures, conversations,
debates, interviews and is now transcribing them in multiple
languages.444https://vanimedia.org/wiki/Table:_Clips_to_subtitle 700+
translators participate in creating subtitles for all 1,080 mini-clips for
108+ languages. Currently, this work has been completed for 25 languages. To
procure subtitles for each clip in multiple languages, translators are
provided with mini-clips and English subtitles that are manually aligned with
each utterance. They use a third-party software named
Dotsub.com555https://dotsub.com and the task of translators is to provide
translation for the corresponding utterance with the help of given
transcription. On average, there are 3-4 translators for each language, and
each clip takes more or less one hour for translation. Each translator has
invested an average of 6 hours every day in translating these clips.
Collectively, the time taken to translate 1,080 clips into 25 languages is
over 45 weeks. Current release of the dataset contains 25 languages (including
transcription) for which transcription/translations are available. For these
languages, we have over 53K utterances of transcription and translations.

Language | Tokens | Types | Percentage
---|---|---|---
English | 500,136 | 6,312 | 83.6
Bengali | 46,933 | 3,907 | 7.84
Sanskrit | 51,246 | 7,202 | 8.56
Total | 598,315 | 17,421 | 100
Table 2: Statistics of code-mixing in Prabhupadavani

### Preprocessing:

In this section, we describe pre-processing steps followed to arrive at the
final version of the data. First, we scrape transcription and their
translations in 25 languages from Dotsub and extracted the corresponding
Youtube links of all audio clips. We use Selenium666https://www.selenium.dev/
web crawler to automatize the process of downloading. We convert those videos
into MP3 audio clips using a third-party
application777https://ytmp3.cc/uu100cc/. We chop the converted audio files
based on the timestamps provided in the subtitle (.srt)
files.888https://pypi.org/project/audioclipextractor/ This process boils down
to 53,000 utterances with their transcription and translation in 25 languages.
Table 4 illustrates the example from Prabhupadavani. In order to obliviate
time and efforts needed for data pre-processing, we provide train, dev and
test set splits using stratified sampling. There are 51,000, 1,000 and 1,000
utterances in train, dev and test set, respectively. We consider the following
dimensions for stratified sampling: (1) different speakers (2) proportion of
intra-sentential and inter-sentential code-mixing.

Language | Type | Examples
---|---|---
English-Sanskrit | Inter | Kṛṣṇa is assuring. ahaṁ tvāṁ sarva-pāpebhyo mokṣayiṣyāmi
| Intra | Sense gratification means udara-upastha-jihvā
Sanskrit-English | Inter |  Īs̀āvāsyam idaṁ sarvam. Everything belongs to God
| Intra |  andhā yathāndhair upanīyamānāḥ and people, leaders.
English-Bengali | Inter | Give up their. Asat-saṅga-tyāga ei vaiṣṇava ācāra
| Intra | Therefore Caitanya Mahāprabhu said guru-kṛṣṇa-kṛpāya
Bengali-English | Inter |  Guru-kṛṣṇa-kṛpāya pāya bhakti-latā-bīja. Then our devotional service is perfect
| Intra |  tāṅhāra nāhika doṣa means he is not faulty.
Table 3: Examples of code-mixing in Prabhupadavani. English, Sanskrit and Bengali are indicated by black, red and violet color, respectively. | Language | Translations
---|---|---
Source | English | You can become Brahman. Brahma-bhūyāya kalpate
Target | Bulgarian |
| Hindi |
| Russian |
| Tamil |
| Gujarati |
Table 4: Sample data point from Prabhupadavani. For a code-mixed utterance, we
show its English transcription (source) and the corresponding translations for
5 languages.

### Code-mixing:

Prabhupadavani is code-mixed across three languages: English, Sanskrit, and
Bengali. Table 2 reports the overall statistics of the code-mixing present in
Prabhupadavani. If we consider the number of tokens in utterances, then it is
mainly dominated by English (83.0%) tokens; however, it is not the case in
terms of a number of types. This attributes to the contrasting nature of
morphology (English vs Sanskrit/Bengali). Code-mixing is categorized into two
classes- (1) inter-sentential: speaker chooses to switch the language after
completion of utterance (2) intra-sentential: speaker switches the language
within an utterance. Table 3 illustrates examples of code-mixing from
Prabhupadavani. Table 5 shows the statistics of both these types of code-
mixing. Mainly speaker explains Sanskrit verses to the English audience;
therefore, transitions between the Sanskrit-English pair is more. However,
sometimes speaker also use Bengali literature to illustrate the points. Thus,
we observe Bengali-English code-switching. Notably, there is no code-switching
between Bengali-Sanskrit because the audience is English speaking.

| Inter-Sentential | Intra-Sentential
---|---|---
English-Sanskrit | 2,356 | 2,338
Sanskrit-English | 2366 | 851
English-Bengali | 339 | 124
Bengali-English | 339 | 0
Table 5: Code-Mixing type for our dataset

### Language diversity:

From typological point of view, Prabhupadavani covers 25 languages (inlcluding
ASR) from 10 language families. They are listed as follows: (1) Indo-
European:- (a) Romance: Italian, Portuguese, Spanish, French (b) Germanic:
German, Swedish, Danish, Dutch, English (c) Baltic: Lithuanian, Latvian (d)
Slavic: Croatian, Polish, Slovak, Russian, Slovenian, Bulgarian, (e) Indo-
Aryan: Bengali, Hindi, Sanskrit, Nepali (f) Indo-Iranian: Persian (Farsi) (2)
Uralic:- (a) Finno-Ugric: Hungarian (3) Austronesian:- (a) Malayo-Polynesian:
Indonesian (4) Dravidian:- Tamil, Telugu. Prabhupadavani contains fusional
languages: Indo-European, Uralic and Agglutinative languages: Austronesian
(Indonesian), Dravidian and Turkic (Turkish). In the former languages,
grammatical markers bear several meanings and for the latter ones, they
exhibit only one meaning at the same time. Except Hindi (Indo-Aryan), all the
languages in our dataset use nominative-accusative marking. Hindi uses
ergative-absolutive marking upto a limited extent. We can also categorize
languages based on the number of grammatical genders. Some languages pose (1)
three genders: German, Russian, Swedish, etc. (2) two genders: French, Spain,
Hindi, etc. (3) no genders: English, Nepali, Persian, Turkish, etc. Based on a
syntactic construct, we can categorize the languages present in
Prabhupadavani- (1) SVO word order: English, Italian, French, Indonesian, etc.
(2) SOV word order: Indo-Aryan, Dravidian and Turkic (3) flexible word order:
Russian, Hungarian.

Thus, the diversified language coverage of Prabhupadavani along with its code-
mixed nature makes it a suitable dataset to investigate various linguistic
phenomena for the speech translation task, ASR and machine translation.

### Applications of dataset:

In this section, we throw some light on possible applications of
Prabhupadavani: (1) The full dataset of Prabhupadavani contains 2,400 hours of
speech. The current release of Prabhupadavani dataset can be utilized to
facilitate automatic subtitling of the remaining part of data in the different
languages. In this way, it will also be helpful to generate relatively larger
dataset. (2) This dataset can provide a fertile soil to investigate on- How to
make existing systems robust to code-mixing phenomenon? Will the gap between
cascade and end-to-end to approaches be closed? Can multi-lingual training
help to address code-mixing phenomenon? Can we train single ST system for all
languages? How robust will be these models on another domain?

##  4 Conclusion and Discussion

In this work, we focused on a code-mixed speech translation dataset. Although
code-mixing is a spoken language phenomenon, not much attention has been given
to address this phenomenon due to the unavailability of such a dataset. Thus,
we released Prabhupadavani, a high-quality multilingual multi-domain code-
mixed ST dataset containing 94 hours of speech data, 130+ speakers, for 25
languages covering 10 language families. The dataset is code-mixed with three
languages: English, Bengali, and Sanskrit. In order to reduce the efforts
needed for pre-processing, we provide stratified data splits for the dataset.
The same dataset can be utilized for the code-mixed machine translation task.
We believe that these efforts will (1) set a fertile soil for investigating
the applicability of existing solutions, (2) help to analyze the kind of
errors existing systems are making, and (3) facilitate researchers to propose
a novel solution to make existing systems robust for code-mixing. We plan to
extend this dataset for 108+ languages.

In code-switching conversations, speakers prefer a specific communication in a
specific language of choice. In that context, the interesting factor is often
when the boundary between languages becomes fuzzy, as in cases where an
English verb stem is used with Spanish morphology. In the case of
Prabhupadavani, the audience does not necessarily speak Sanskrit or Bengali,
and the code-switching is primarily quotation or explanation. That is
interesting as a phenomenon, but it is not the same as dual bilingual code-
switching. We believe future work may help to get deep insights into this
phenomenon.

### Ethics Statement:

We do not foresee any ethical concerns with the work presented in this
manuscript. We have taken the consent of the Vanipedia team and Bhaktivedanta
Book Trust International to use translations and audio in our dataset.

## Acknowledgements

We would like to thank the Vanipedia team (https://vanipedia.org/) of 700+
translators for establishing this multi-lingual database for us to develop. We
thank the Bhaktivedanta Book Trust International for permitting us to use
Prabhupadavani audio in our dataset. We are grateful to Manish Gupta,
Microsoft, for helping us with insightful discussions. We would like to thank
the anonymous reviewers for their constructive feedback towards improving this
work. The TCS Fellowship supports the first author’s work under Project
TCS/EE/2011191P.

## References

  * Ansari et al. (2020) Ebrahim Ansari, Amittai Axelrod, Nguyen Bach, Ondřej Bojar, Roldano Cattoni, Fahim Dalvi, Nadir Durrani, Marcello Federico, Christian Federmann, Jiatao Gu, et al. 2020.  Findings of the iwslt 2020 evaluation campaign.  In _Proceedings of the 17th International Conference on Spoken Language Translation_ , pages 1–34.
  * Bansal et al. (2017) Sameer Bansal, Herman Kamper, Adam Lopez, and Sharon Goldwater. 2017.  Towards speech-to-text translation without speech recognition.  In _Proceedings of the 15th Conference of the European Chapter of the Association for Computational Linguistics: Volume 2, Short Papers_ , pages 474–479, Valencia, Spain. Association for Computational Linguistics.
  * Bérard et al. (2018) Alexandre Bérard, Laurent Besacier, Ali Can Kocabiyikoglu, and Olivier Pietquin. 2018.  End-to-end automatic speech translation of audiobooks.  In _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 6224–6228. IEEE.
  * Chakravarthi et al. (2020) Bharathi Raja Chakravarthi, Vigneshwaran Muralidaran, Ruba Priyadharshini, and John P McCrae. 2020.  Corpus creation for sentiment analysis in code-mixed tamil-english text.  _arXiv preprint arXiv:2006.00206_.
  * Cho et al. (2013) Eunah Cho, Christian Fügen, Teresa Herrmann, Kevin Kilgour, Mohammed Mediani, Christian Mohr, Jan Niehues, Kay Rottmann, Christian Saam, Sebastian Stüker, et al. 2013.  A real-world system for simultaneous translation of german lectures.  In _INTERSPEECH_ , pages 3473–3477.
  * Dhar et al. (2018) Mrinal Dhar, Vaibhav Kumar, and Manish Shrivastava. 2018.  Enabling code-mixed translation: Parallel corpus creation and mt augmentation approach.  In _Proceedings of the First Workshop on Linguistic Resources for Natural Language Processing_ , pages 131–140.
  * Di Gangi et al. (2019) Mattia A. Di Gangi, Roldano Cattoni, Luisa Bentivogli, Matteo Negri, and Marco Turchi. 2019.  MuST-C: a Multilingual Speech Translation Corpus.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 2012–2017, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Duong et al. (2016) Long Duong, Antonios Anastasopoulos, David Chiang, Steven Bird, and Trevor Cohn. 2016.  An attentional model for speech translation without transcription.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 949–959, San Diego, California. Association for Computational Linguistics.
  * Gonen and Goldberg (2018) Hila Gonen and Yoav Goldberg. 2018.  Language modeling for code-switching: Evaluation, integration of monolingual data, and discriminative training.  _arXiv preprint arXiv:1810.11895_.
  * Iranzo-Sánchez et al. (2020) J. Iranzo-Sánchez, J. A. Silvestre-Cerdà, J. Jorge, N. Roselló, A. Giménez, A. Sanchis, J. Civera, and A. Juan. 2020.  Europarl-st: A multilingual corpus for speech translation of parliamentary debates.  In _ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 8229–8233.
  * Khanuja et al. (2020) Simran Khanuja, Sandipan Dandapat, Sunayana Sitaram, and Monojit Choudhury. 2020\.  A new dataset for natural language inference from code-mixed conversations.  _arXiv preprint arXiv:2004.05051_.
  * Lyu et al. (2015) Dau-Cheng Lyu, Tien-Ping Tan, Eng-Siong Chng, and Haizhou Li. 2015.  Mandarin–english code-switching speech corpus in south-east asia: Seame.  _Language Resources and Evaluation_ , 49(3):581–600.
  * Nakayama et al. (2019) Sahoko Nakayama, Takatomo Kano, Andros Tjandra, Sakriani Sakti, and Satoshi Nakamura. 2019.  Recognition and translation of code-switching speech utterances.  In _2019 22nd Conference of the Oriental COCOSDA International Committee for the Co-ordination and Standardisation of Speech Databases and Assessment Techniques (O-COCOSDA)_ , pages 1–6.
  * Niehues et al. (2021) Jan Niehues, Elizabeth Salesky, Marco Turchi, and Matteo Negri. 2021.  Tutorial proposal: End-to-end speech translation.  In _Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Tutorial Abstracts_ , pages 10–13, online. Association for Computational Linguistics.
  * Post et al. (2013) Matt Post, Gaurav Kumar, Adam Lopez, Damianos Karakos, Chris Callison-Burch, and Sanjeev Khudanpur. 2013.  Improved speech-to-text translation with the fisher and callhome spanish–english speech translation corpus.  In _Proc. IWSLT_.
  * Ruiz et al. (2015) Nicholas Ruiz, Qin Gao, Will Lewis, and Marcello Federico. 2015.  Adapting machine translation models toward misrecognized speech with text-to-speech pronunciation rules and acoustic confusability.  In _Interspeech_. ISCA-International Speech Communication Association.
  * Salesky et al. (2021) Elizabeth Salesky, Matthew Wiesner, Jacob Bremerman, Roldano Cattoni, Matteo Negri, Marco Turchi, Douglas W. Oard, and Matt Post. 2021.  Multilingual tedx corpus for speech recognition and translation.  In _Proceedings of Interspeech_.
  * Singh et al. (2018a) Kushagra Singh, Indira Sen, and Ponnurangam Kumaraguru. 2018a.  A twitter corpus for hindi-english code mixed pos tagging.  In _Proceedings of the Sixth International Workshop on Natural Language Processing for Social Media_ , pages 12–17.
  * Singh et al. (2018b) Vinay Singh, Deepanshu Vijay, Syed Sarfaraz Akhtar, and Manish Shrivastava. 2018b.  Named entity recognition for hindi-english code-mixed social media text.  In _Proceedings of the seventh named entities workshop_ , pages 27–35.
  * Sperber et al. (2017) Matthias Sperber, Jan Niehues, and Alex Waibel. 2017.  Toward robust neural machine translation for noisy input sequences.  In _International Workshop on Spoken Language Translation (IWSLT)_ , page 18.
  * Srivastava and Singh (2021) Vivek Srivastava and Mayank Singh. 2021.  Challenges and considerations with code-mixed nlp for multilingual societies.
  * Tsvetkov et al. (2014) Yulia Tsvetkov, Florian Metze, and Chris Dyer. 2014.  Augmenting translation models with simulated acoustic confusions for improved spoken language translation.  In _Proceedings of the 14th Conference of the European Chapter of the Association for Computational Linguistics_ , pages 616–625, Gothenburg, Sweden. Association for Computational Linguistics.
  * Wang et al. (2020) Changhan Wang, Juan Pino, Anne Wu, and Jiatao Gu. 2020.  CoVoST: A diverse multilingual speech-to-text translation corpus.  In _Proceedings of The 12th Language Resources and Evaluation Conference_ , pages 4197–4203, Marseille, France. European Language Resources Association.
  * Weiss et al. (2017) Ron J Weiss, Jan Chorowski, Navdeep Jaitly, Yonghui Wu, and Zhifeng Chen. 2017.  Sequence-to-sequence models can directly translate foreign speech.  _arXiv preprint arXiv:1703.08581_.
  * Zanon Boito* et al. (2020) Marcely Zanon Boito*, William Havard*, Mahault Garnerin, Eric Le Ferrand, and Laurent Besacier. 2020.  MaSS: A Large and Clean Multilingual Corpus of Sentence-aligned Spoken Utterances Extracted from the Bible.  In _Proceedings of the 12th Language Resources and Evaluation Conference_ , pages 6486–6493, Marseille, France. European Language Resources Association.
  * Zhang et al. (2021) Ruiqing Zhang, Xiyang Wang, Chuanqiang Zhang, Zhongjun He, Hua Wu, Zhi Li, Haifeng Wang, Ying Chen, and Qinfei Li. 2021.  Bstc: A large-scale chinese-english speech translation dataset.  _arXiv preprint arXiv:2104.03575_.


# In-gap states induced by magnetic impurities on wide-band s-wave
superconductors: self-consistent calculations

Divya Jyoti  Centro de Física de Materiales CFM/MPC (CSIC-UPV/EHU), 20018
Donostia-San Sebastián, Spain  Donostia International Physics Center (DIPC),
20018 Donostia-San Sebastián, Spain  Deung-Jang Choi  Centro de Física de
Materiales CFM/MPC (CSIC-UPV/EHU), 20018 Donostia-San Sebastián, Spain
Donostia International Physics Center (DIPC), 20018 Donostia-San Sebastián,
Spain  Nicolás Lorente<EMAIL_ADDRESS>Centro de Física de Materiales
CFM/MPC (CSIC-UPV/EHU), 20018 Donostia-San Sebastián, Spain  Donostia
International Physics Center (DIPC), 20018 Donostia-San Sebastián, Spain

###### Abstract

The role of self-consistency in Bogoliubov-de Gennes equations is frequently
underestimated in the investigation of in-gap states created by magnetic
impurities in s-wave superconductors. Our research focuses on the impact of
self-consistency on the in-gap states produced by magnetic stuctures on
superconductors, specifically evaluating the density of states, the in-gap
bands, and their topological attributes. Here, we show results ranging from
single impurity to finite chains, and infinite ferromagnetic spin chains in
wide-band s-wave superconductors. These results show that the order parameter
contains important information regarding quantum phase transitions and their
topological nature, underscoring the importance of self-consistency in such
studies.

##  I Introduction

The Bogoliubov-de Gennes (BdG) equations have played a pivotal role in
simplifying and advancing the theoretical examination of magnetic impurities
within superconductors. The method and examples have been clearly exposed in
Refs. Balatsky  _et al._ , 2006; Zhu, 2016. The BdG framework offers a means
to address the complex many-body problem associated with superconductivity by
employing a mean-field approximation of the pairing interaction. The final
equations amount to a diagonalization of a one-body effective Hamiltonian
Balatsky  _et al._ (2006). It is essential to recognize that the mean-field
nature of the BdG equations necessitates the self-consistent evaluation of the
effective pairing interaction.

In a superconductor, a single magnetic impurity leads to weakening of the
Cooper pair binding energy. This can give rise to the appearance of one-
particle states in the superconducting gap Yu (1965); Shiba (1968); Rusinov
(1969); Balatsky  _et al._ (2006). Approximating the impurity-substrate
interaction by the Kondo Hamiltonian Schrieffer and Wolff (1966), the minimum
energy required for the emergence of quasiparticle excitations decreases when
the Kondo-exchange coupling interaction, $J$ increases. Near a critical value
of the exchange interaction ($J_{c}$), the superconducting condensate becomes
thermodynamically unstable against the formation of quasiparticles Sakurai
(1970); Salkola  _et al._ (1997); Nagai  _et al._ (2013). Henceforth, a
superconductor undergoes a quantum phase transition (QPT) and achieves a spin-
polarized state. Theoretical studies Salkola  _et al._ (1997); Flatté and
Byers (1997, 1997); Flatté (2000); Hoffman  _et al._ (2015); Björnson  _et
al._ (2017); Theiler  _et al._ (2019); Villas  _et al._ (2020) and
experimental realizations of QPT have been reported by several groups Hatter
_et al._ (2015); Farinacci  _et al._ (2018); Huang  _et al._ (2020); Liebhaber
_et al._ (2021); Karan  _et al._ (2022); Liu  _et al._ (2023); Zhou  _et al._
(2022); Uldemolins  _et al._ (2023) where a range of Kondo-exchange couplings
are achieved based on how adatoms attach to specific sites on the substrate
Hatter  _et al._ (2015); Franke  _et al._ (2011); Ji  _et al._ (2010).
Alternatively, by establishing an electrical control in Josephson junction
systems Maurand  _et al._ (2012); Delagrange  _et al._ (2015, 2016); Deacon
_et al._ (2010) or by changing tip-sample distances in scanning tunneling
microscopy (STM) experiments Farinacci  _et al._ (2018), these couplings can
be tuned. The strength and nature of such couplings can have a direct impact
on the superconducting state and therefore incorporating self-consistency may
be necessary to obtain an accurate description of such systems.

Many calculations have demonstrated that the influence of self-consistency is
predominantly quantitative, suggesting that the underlying physics can be
accurately captured using simpler and more computationally efficient non self-
consistent approaches von Oppen and Franke (2021); Salkola  _et al._ (1997);
Flatté and Byers (1997); Schecter  _et al._ (2016); Karan  _et al._ (2022).
However, there exist specific scenarios where the adoption of self-consistent
solutions is essential Žonda  _et al._ (2015); Rozhkov and Arovas (1999). Levy
Yeyati et al. Levy Yeyati  _et al._ (1995) illustrated that concerns regarding
the violation of particle-number conservation in the BdG formalism are
resolved when self-consistency is incorporated in electric current
evaluations. Particularly critical are instances involving QPTs, where the
interaction between impurity and superconductor is so substantial that it
alters the ground state of the entire system, indicating a zero-temperature
QPT and underscoring the significance of self-consistency in these analyses.
Notably, the closure of the superconducting gap, a precursor to a QPT, is an
example of behavior that qualitatively differs when including self-consistency
Salkola  _et al._ (1997); Franke  _et al._ (2011); Hatter  _et al._ (2015);
Heinrich  _et al._ (2018). This gap closure also forms a necessary criterion
for the emergence of topological quantum phase transitions (TQPTs) observed in
spin chains on s-wave superconductors Choy  _et al._ (2011); Pientka  _et al._
(2013); Choi  _et al._ (2019). Consequently, examining topological invariants
and phases through self-consistent methods becomes a critical aspect when
employing BdG equations Björnson  _et al._ (2015); Christensen  _et al._
(2016); Björnson  _et al._ (2017); Theiler  _et al._ (2019).

In this study, we expand on our prior theoretical work on spin impurities in
s-wave superconductors Mier  _et al._ (2021a, 2022), focusing on the self-
consistent evaluation of the order parameters. Our previous research delved
into topological transitions in wide-band s-wave superconductors Mier  _et
al._ (2021a, 2022), relating them to experimentally observed non-zero edge
states Schneider  _et al._ (2022). In our current work, we utilize self-
consistent solutions not only to validate and reinforce the findings of our
earlier studies but also to deepen our comprehension of QPTs. This approach is
inspired by literature Björnson  _et al._ (2015, 2017); Theiler  _et al._
(2019) that underscores the significance of self-consistency in determining
the topological properties of spin chains on superconductors. However,
contrasting views, especially from Ref. Christensen  _et al._ , 2016, suggest
a more indirect role of self-consistency, impacting the phenomena through
renormalization of critical exchange interactions and consequent gap closure,
leading to QPTs.

The structure of our paper is as follows: Initially, we present our numerical
method, based on the techniques developed by Flatté and Byers Flatté and Byers
(1997, 1997) in a discretized form Mier  _et al._ (2021a). We first test the
correctness of the approach by retrieving the BCS temperature dependence of a
bulk wide-band superconductor. Next, we add a magnetic impurity and recover
the properties of the induced in-gap states. Subsequently, we conduct a
systematic investigation into the impact of progressively increasing the
number of impurities on a superconductor, leading to an in-depth analysis of
ferromagnetic infinite chains and their associated topological invariants. We
show that the use of the order parameter as a diagnostic tool proves to be
exceptionally insightful for detecting the emergence of QPTs, particularly in
the realm of TQPTs. A principal finding of our research is that self-
consistency does not alter the topological properties previously identified in
non-self-consistent studies and has a minimal effect on the density of states
in realistic wide-band s-wave superconductors.

##  II Theoretical methodology

###  II.1 Bogoliubov-de Gennes equations using Green’s functions

We solve the BdG using the Nambu formalism and a discretization of the spatial
continuumPientka  _et al._ (2013); Mier  _et al._ (2021b); Cuevas  _et al._
(1996). We start from the Hamiltonian for a pristine BCS superconductor given
by,

$\hat{H}_{BCS}=\xi_{k}\tau_{0}\sigma_{3}+\Delta\tau_{2}\sigma_{2},$ (1)

where $\sigma_{i}$ ($\tau_{i}$) are the Pauli matrices acting on the spin
(particle) subspace, $\xi_{k}$ is the energy from the Fermi level
($\xi_{k}=\epsilon_{k}-E_{F}$) and $\Delta$ is the superconducting gap or
order parameter. The previous Hamiltonian is written in the 4-dimensional
Nambu basis:
$\Psi=(\psi_{\uparrow},\psi_{\downarrow},\psi_{\uparrow}^{\dagger},\psi_{\downarrow}^{\dagger})^{T}$.

To model the experimental system, we add the Hamiltonian describing the
magnetic impurities Flatté (2000); Flatté and Byers (1997, 1997), but assuming
strictly localized interactions Pientka  _et al._ (2013).

$\hat{H}=\hat{H}_{BCS}+\hat{H}_{impurity}=\hat{H}_{BCS}+\sum_{j}^{N}(U_{j}\tau_{3}+J_{j}\vec{S_{j}}\cdot\vec{\alpha})\;$
(2)

with
$\vec{\alpha}=\frac{1+\tau_{3}}{2}\vec{\sigma}+\frac{1-\tau_{3}}{2}\sigma_{2}\vec{\sigma}\sigma_{2}$,
where $\vec{\sigma}$ is the spin operator Shiba (1968).

The interaction contains an exchange coupling, with strength $J_{j}$, and a
non-magnetic potential scattering term, $U_{j}$, per impurity $j$. The
impurity spin is assumed to be a classical vector, $\vec{S_{j}}$, within the
classical-spin approximation, see II.3.

The Hamiltonian is completed by a Rashba term:

$\displaystyle\hat{H}_{Rashba}$ $\displaystyle=$ $\displaystyle
i\frac{\alpha_{R}}{2a}\sum_{i,j,\alpha,\beta}[\hat{c}^{\dagger}_{i+1,j,\alpha}(\sigma_{2})_{\alpha,\beta}\hat{c}_{i,j,\beta}$
(3) $\displaystyle-$
$\displaystyle\hat{c}^{\dagger}_{i,j+1,\alpha}(\sigma_{1})_{\alpha,\beta}\hat{c}_{i,j,\beta}+h.c.]$

where ${\alpha,\beta}$ are the spin indexes. Here, ${i,j}$ represent the
impurity index in a two-dimensional lattice. This interaction couples spins on
next-nearest-neighbour sites. The lattice parameter of the substrate is $a$,
and the factor of $2a$ comes from a finite-difference scheme to obtain the
above discretized version of the Rashba interaction.

We solve the BdG equations using Green’s functions, ${G}^{\nu,\mu}_{i,j}$
evaluated for the sites $i$ and $j$ for the Nambu components $\nu$ and $\mu$
by solving Dyson’s equation:

$\hat{G}=[\hat{G}_{BCS}^{-1}-\hat{H}_{I}]^{-1}$ (4)

where $\hat{G}_{BCS}$ is the retarded Green’s operator for the BCS Hamiltonian
from Eq. (1) and $\hat{H}_{I}=\hat{H}_{impurity}+\hat{H}_{Rashba}$. In the
present study, we will be interested in the effect of the impurities on the
local gap $\Delta_{i}$, for this, we use the local Green’s function. The
Green’s function that is diagonal on spatial indices, $i$ and $j=i$, is given
by the usual local BCS Green’s function Vernier  _et al._ (2011) (following
usual convention $\hbar=1$ unless otherwise specified),

$\displaystyle{G_{i,i}}_{BCS}(\omega)$ $\displaystyle=$
$\displaystyle-\frac{\pi
N_{0}Sgn[Re(\omega)Im(\omega)]}{\sqrt{\Delta^{2}-\omega^{2}}}$ (5)
$\displaystyle\times$ $\displaystyle\begin{pmatrix}\omega&0&0&-\Delta\\\
0&\omega&\Delta&0\\\ 0&\Delta&\omega&0\\\ -\Delta&0&0&\omega\end{pmatrix}.$

Here, $N_{0}$ corresponds to the normal-metal density of electronic states
evaluated at the Fermi energy. The non-local Green’s function can be
approximated by a well-known analytical expression, see for example Refs.
Flatté and Byers, 1997; Meng  _et al._ , 2015; Schecter  _et al._ , 2016; Mier
_et al._ , 2021a, yielding

$\displaystyle{G}_{BCS}(\vec{r},\vec{r^{\prime}},\omega)$ $\displaystyle=$
$\displaystyle-\frac{\pi
N_{0}}{k_{F}r}e^{\frac{-\sqrt{\Delta^{2}-\omega^{2}}r}{\pi\xi\Delta}}\times(\cos(k_{F}r)\tau_{0}\sigma_{3}\;$
(6) $\displaystyle+$
$\displaystyle\frac{\sin(k_{F}r)}{\sqrt{\Delta^{2}-\omega^{2}}}(\omega\tau_{0}\sigma_{0}+\Delta\tau_{2}\sigma_{2})),$

where $r=|\vec{r}-\vec{r^{\prime}}|$, $k_{F}$ is the Fermi wave-vector and
$\xi=\frac{\hbar v_{F}}{\pi\Delta}$ is coherence length.

Within this formalism the condition for self-consistency comes from the mean-
field value of the gap function Zhu (2016):

$\Delta(\vec{r},\vec{r^{\prime}})=V(\vec{r}-\vec{r^{\prime}})\langle\psi_{\downarrow}(\vec{r})\psi_{\uparrow}(\vec{r^{\prime}})\rangle$
(7)

where the unknown pairing potential $V(\vec{r}-\vec{r^{\prime}})$ is
approximated by a local potential. Additionally, a cutoff in energies is
applied in the evaluation of the ground-state average, just taking quasi-
particle states that are within the Debye frequency, $\omega_{D}$. The average
is easy to evaluate using Green’s functions, in particular we calculate:

$\Delta_{i}=\tilde{V}\int^{\hbar\omega_{D}}_{-\hbar\omega_{D}}f(\omega)Im\\{G^{1,4}_{i,i}(\omega)-G^{2,3}_{i.i}(\omega)\\}d\omega.$
(8)

Here, $f(\omega)$ is the Fermi function at temperature $T$, evaluated for
energy $\hbar\omega$, $G^{1,4}$ and $G^{2,3}$ are the anomalous Green’s
function components for Nambu indices 1, 4 and 2, 3, respectively. The
effective interaction with all possible multiplicative constants, $\tilde{V}$,
is assumed to be homogeneous and independent of the impurity interactions. It
can be easily obtained from the value $\Delta_{0}$ for the pristine
superconductor and the BCS Green’s function, Eq. (5). See appendix A.

The new $\Delta$ is computed by using the Green’s functions solving Dyson’s
equations for the full system. The change in the order parameter,
$\delta\Delta(\vec{r})$ is then used as a self-energy to compute the new
inhomogeneous Green’s function following Refs. Flatté and Byers, 1997, 1997.
This procedure is repeated until $\Delta(\vec{r})$ changes less than a certain
tolerance.

Now for the evaluation of infinite chains, we rewrite the Green’s functions in
terms of their k-space values. For this, we take into account that

$G_{i,j}(\omega)=\frac{1}{N_{k}}\sum_{k}G(k,\omega)e^{ik(R_{i}-R_{j})},$

where $N_{k}$ refers to the number of k-points of the calculation and
$R_{i},R_{j}$ describes discrete position coordinates.

Then from Eq. (8), we obtain

$\displaystyle\Delta_{i}=\frac{\tilde{V}}{N_{k}}\int^{\hbar\omega_{D}}_{-\hbar\omega_{D}}f(\omega)\sum_{k}Im\\{G^{1,4}(k,\omega)-$
$\displaystyle G^{2,3}(k,\omega)\\}d\omega,$ (9)

where $\tilde{V}(\vec{r}-\vec{r}^{\prime})$ is assumed homogeneous in an
isotropic s-wave type superconductor, and constant. Since we are representing
an infinite and homogeneous system, the above order parameter $\Delta_{i}$
does not depend on the site $i$ and is constant. We can actually obtain the
k-resolved expression by removing the sum over k points in Eq. (9), this
yields

$\displaystyle\Delta_{k}=\tilde{V}\int^{\hbar\omega_{D}}_{-\hbar\omega_{D}}f(\omega)Im\\{G^{1,4}(k,\omega)-$
$\displaystyle G^{2,3}(k,\omega)\\}d\omega,$ (10)

such that the average order parameter is
$\Delta_{i}=\frac{1}{N_{k}}\sum_{k}\Delta_{k}$. The advantage of the above
expression, Eq. (9) is that the k-space Green’s functions can be easily
evaluated Peng  _et al._ (2015); Schecter  _et al._ (2016).

The self-consistent calculation for the infinite chain is performed by
computing $\tilde{V}$ using the homogeneous superconductor. This means setting
all interactions to zero except the Rashba coupling. Then the value of
$\Delta_{i}$ fixes the coupling $\tilde{V}$. Next, we proceed to define a
zero-iteration $\Delta_{k}$. This is done by evaluating Eq. (10) for the free
Green’s functions, i.e. without interactions and only the Rashba coupling.
Next, the interactions are included using Dyson’s equations and the above
self-consistent loop is performed by defining a new self-energy that contains
the change of $\Delta_{k}$ per k-point and per iteration.

###  II.2 Numerical implementation

We use a lattice approach discretizing the spatial dependence. The size of the
discrete lattice step, $\Delta r$, is quite robust against different values.
The approximations leading to the analytical BCS Green’s function, imply that
$\Delta r\times k_{F}>1$, where $k_{F}$ is the Fermi wave vector.
Additionally, special care must be paid in converging the Rashba interaction,
due to the evaluation of the Rashba gradients on the grid. Yet, the robustness
of the full approach comes from using a $G_{BCS}$ that is defined everywhere
in space. Here, we use 3-D Green’s functions because using 2-D Green’s
functions just bring small numerical changes Brydon  _et al._ (2015) despite
the different in-gap state decay Ménard  _et al._ (2015), that does not affect
the convergence of the order parameter, the main object of our present study.

A consequence of these approximations for the real-space BCS Green’s function
is that $G_{BCS}(k,\omega)$ is defined only for $k<k_{F}$. As a consequence,
the Brillouin-zone $\pi/a$ where $a$ is the distance between impurities,
should be smaller than $k_{F}$ to have the Green’s function defined for all
k-points. A simple compromise is to take $\Delta r\approx a$ for a large
enough $k_{F}$, please see the corresponding discussion in Ref. Mier  _et al._
, 2021a. As in Ref. Pientka  _et al._ , 2013; Schneider  _et al._ , 2021, our
calculations are best suited for spin chains in the dilute-impurity limit.

###  II.3 Classical-spin approximation

The classical-spin approximation was used by Yu, Shiba and Rusinov Yu (1965);
Shiba (1968); Rusinov (1969) to prove and characterize the appearance of in-
gap states in the presence of a magnetic impurity. However, it was soon
recognized that the impurity spin is a quantum operator and spin fluctuations
occur that can eventually produce the Kondo effect Matsuura (1977). In
addition to spin fluctuations, spin-orbit coupling can lift the degeneracy of
the impurity’s spin states, presenting different states that can be accessible
at different energies leading to a different picture of the low-energy physics
žitko  _et al._ (2011); Schmid  _et al._ (2022).

The classical spin approximation gives sensible results in the large-spin
limit given by $S_{i}\rightarrow\infty$ with $J\times S_{i}$ finite. Reference
žitko  _et al._ , 2011 shows that for transition-metal spins, this limit is
not satisfied. And despite the quenching of spin fluctuations by the presence
of the superconducting gap, at $J\times S_{i}$ large enough to yield bound in-
gap states, spin fluctuations become available again.

In the present work, we focus on the properties of in-gap states without
studying the nature of the many-body ground state Moca  _et al._ (2008) or
spin transitions in the impurity von Oppen and Franke (2021). In these
conditions, the classical-spin approximation yields a sufficient description
that is easily accessible with mean-field theories such as the BdG approach
used here.

##  III Temperature effects

The self-consistent BdG approach readily captures the well-established BCS
behavior of the superconducting gap as a function of temperature. In Fig. 1,
we present the outcomes for the Bi2Pd s-wave superconductor, utilizing our
free-electron BdG equations to simulate the BCS superconductor. To enhance
comparability, we include a BCS curve, derived from the subsequent expression,

$\Delta(T)=\Delta_{0}\tanh\Biggl{(}\ 1.74\sqrt{\frac{T_{c}}{T}-1}\Biggl{)}.$
(11)

From Fig. 1, we obtain a critical temperature of 5.6 K, which is in good
agreement with the experimental one $T_{c}=5.4$ K as reported in Ref. Imai
_et al._ , 2012. To obtain this value, the $T=0$ value of the gap,
$\Delta_{0}$, has been fixed to the experimental one, 0.76 meV. The only other
parameter is the free-electron density fixed by the Fermi wave vector $k_{F}$,
that has been taken as 0.183 $a_{0}^{-1}$ following the experimental data of
Ref. Herrera  _et al._ , 2015. The results are largely independent of the
Dynes parameter used as imaginary part of $\omega$ in the Green’s function
expressions Dynes  _et al._ (1978). For the present calculation a Dynes
parameter of $\Gamma=0.05$ meV was used.

Figure 1: Self-consistently evaluated superconducting gap, $\Delta$ computed
as a function of temperature for a BCS model of an s-wave superconductor,
Bi2Pd. A BCS curve following Eq. (11) is also shown for comparison. To obtain
a critical temperature $T_{c}=5.6$ K (in good agreement with the experimental
$T_{c}=5.4$ K Imai  _et al._ (2012)), we only used $\Delta_{0}=0.76$ meV and
the Bi2Pd Fermi wave vector, $k_{F}=0.183~{}a_{0}^{-1}$, as inputs.

##  IV Single magnetic impurity on a superconductor

Let us first study a single magnetic impurity using the above formalism. Here,
we will briefly characterize the in-gap state as with respect to its extension
in the superconductor and to the electron-hole symmetry built-in the
equations.

###  IV.1 Extension of the in-gap states

In simple arguments, there are two different length scales associated to a
single Shiba state Salkola  _et al._ (1997); Flatté and Byers (1997). The
longer scale is due to the actual extension of correlations in the
superconductor. It is related to the coherence length, $\xi$, given by the
decay length $\xi^{\prime}=\xi/\sqrt{1-\frac{\omega^{2}}{\Delta^{2}}}$, where
$\omega$ is the pole of the Shiba state or the quasienergy of the Shiba
excitation.

The second length scale is the immediate decay of the Shiba wavefunction as
soon as we move away from the impurity center. Our model is particularly
unrealistic on the impurity site due to the presence of a $\delta$-function in
the interaction that describes a very localized magnetic interaction (Eq. 5).
Within this approximation, we see that the Shiba wave function has a decay
that is actually given by $(\lambda_{F}/r)^{2}$, where the Fermi wavelength is
given by $\lambda_{F}=2\pi/k_{F}$ and $r$ is the spatial distance to the
impurity’s center.

We can easily retrieve this behavior from our theoretical model. The
wavefunction can be obtained from the Green’s function using the Lippmann-
Schwinger equation:

$\psi(\vec{r})=\psi_{0}(\vec{r})+\int
d\vec{r^{\prime}}G_{BCS}(\vec{r},\vec{r^{\prime}},\omega)H_{impurity}(\vec{r^{\prime}})\psi(\vec{r^{\prime}}).$
(12)

For Shiba states, $\psi_{0}(\vec{r})=0$ because there are no single-
quasiparticle states in the gap of the pristine superconductor. Using the
locality of the interaction with the impurity,
$H_{impurity}(\vec{r})\propto\delta(\vec{r})$, we easily simplify Eq. (12) to

$\psi(\vec{r})\propto G_{BCS}(\vec{r},0,\omega)\psi(0),$ (13)

where $G_{BCS}(\vec{r},0,\omega)$ can be calculated from Eq. (6).

Then, away from the impurity’s center the wavefunction decays as the bare BCS
Green’s function. Using the typical BCS expression for the Green’s function
Flatté and Byers (1997); Meng  _et al._ (2015); Schecter  _et al._ (2016);
Mier  _et al._ (2021a), we see that this decay is of the form
$(\lambda_{F}/r)\,\exp(-r/\xi^{\prime})$. The density scales as the square of
the wave function and hence, as $(\lambda_{F}/r)^{2}\,\exp(-2r/\xi^{\prime})$
in agreement with the discussions found in Refs. Salkola  _et al._ , 1997;
Flatté and Byers, 1997; Ménard  _et al._ , 2015. For a truly 2-D system, the
decay of the Shiba states would rather scale as $1/r$, following the arguments
of Ref. Ménard  _et al._ , 2015.

Then, two length scales are expected for the Shiba states found inside a
superconducting gap. The shorter length scale of the Shiba state is then
related to the inverse of the Fermi wave vector. In BCS superconductors, we
can expect that larger electron densities (larger $k_{F}$) leads to short-
ranged Shiba states. This scale is the one responsible for the hybridization
of Shiba states and the formations of bands needed for the appearance of a
topological phase. Then, impurities on less dense superconductors, can be
placed at larger distances and still create a rich structure of in-gap bands.
Indeed, Ref. Mier  _et al._ , 2021b shows that Majoranas can be found at the
edges of spin chains under these conditions, that perfectly match Cr chains on
the s-wave superconductor Bi2Pd.

###  IV.2 Particle-hole symmetry for in-gap states

Precise experiments carried on with the STM have shown that in-gap states are
not electron-hole symmetric Ruby  _et al._ (2016); Choi  _et al._ (2017).
Indeed, despite having an electron and a hole component, the in-gap states
show very different spatial distributions for the electron and hole components
at the exact same energies (with opposite sign). However, BdG equations are
electron-hole symmetric. This means Zhu (2016) that for every solution at
energy $E$, with wavefunction
${\phi_{i}}=(u_{i\uparrow},u_{i\downarrow},v_{i\uparrow},v_{i\downarrow})^{T}$,
there is a solution at energy $-E$ with wavefunction
${\phi_{i}}^{\prime}=(-v^{*}_{i\uparrow},u^{*}_{i\downarrow},v^{*}_{i\downarrow},-u^{*}_{i\uparrow})^{T}$,where
the electron, $u_{\uparrow}$, and hole, $v_{\downarrow}$, components of the
wave functions are in principle different as experimentally found Ruby  _et
al._ (2016); Choi  _et al._ (2017).

Before the experimental realizations, Ref. Flatté and Byers, 1997 showed ample
evidence of cases where the electron and hole components of Shiba states are
not equal in magnitude. The authors of Ref. Flatté and Byers, 1997, 1997;
Yazdani  _et al._ , 1997 argue that even for $U_{i}=0$ in Eq. (2) the electron
and hole components can have different magnitudes. This statement is contrary
to what is stated in Ref. Salkola  _et al._ , 1997 and contrary to what is
widely accepted that electron-hole symmetry is only attained for $U_{i}=0$.

As argued in Ref. Flatté and Byers, 1997, the reason behind the $U_{i}=0$
electron-hole symmetry is due to the local $\delta$-like approximation for the
impurity interactions. It is easy to see that this is indeed, the case from
Eq. (5) that a local potential will not lift the electron-hole symmetry at the
impurity site. However, as long as the Green’s function is not evaluated for
the same site ($\vec{r}\neq\vec{r}^{\prime}$), the electron-hole symmetry is
lifted. Indeed, the particle component in Eq. (13) depends on

$\cos(k_{F}r)+\frac{\omega}{\sqrt{\Delta_{0}^{2}-\omega^{2}}}\sin(k_{F}r),$

while the hole component depends on

$-\cos(k_{F}r)+\frac{\omega}{\sqrt{\Delta_{0}^{2}-\omega^{2}}}\sin(k_{F}r),$

basically showing a phase shift between electron and hole for strongly bound
in-gap states ($\omega\approx 0$). Let us say a word of caution about using
the above expression very close to the impurity where the Green’s function,
Eq. (6), fails as explained in section II.2 and Ref. Mier  _et al._ , 2021a.
For a more realistic non-local impurity potential, such as the one’s of Ref.
Flatté and Byers, 1997, the particle and hole components can be different even
at the impurity site and with zero non-magnetic scattering, $U_{j}=0$.

###  IV.3 Quantum phase transition

A quantum phase transition (QPT) refers to a fundamental change in the ground
state of a many-body system. To observe and study such transitions, we often
track specific system parameters, such as the order parameter or gap, denoted
as $\Delta$ in Eq. (7) in our current context.

When a single magnetic impurity is present on a superconductor, the cooper
pair binding energy is weakened and the emergence of in-gap quasiparticle
states follows. As the magnetic interaction of the impurity with the
superconducting electrons becomes prominent, the order parameter undergoes an
abrupt change and a QPT is realized. The discontinuous change in the order
parameter characterizes these transitions as first-order. In the present case,
the QPT signals a new magnetic ground state Salkola  _et al._ (1997); Sakurai
(1970).

Figure 2: Order parameter, $\Delta$ as a function of exchange coupling
strength, $J$, for a range of temperature values below $T_{c}$. The
calculations shown here are for a single magnetic impurity on a 2D
superconducting lattice. These calculations have been performed without
potential scattering, $U=0$, or spin-orbit coupling, $\alpha=0$. The used
Fermi wave vector is $k_{F}=0.183\,a_{0}^{-1}$.

Figure 2 shows the order parameter as a function of effective impurity-sample
interactions $J$ given in Eq. (2). At zero temperature, a QPT takes place at a
critical value of the exchange coupling, $J_{c}=3.46$ eV, when the gap at the
impurity site becomes zero. As before, we have particularized to Bi2Pd as
superconductor. We have used the following parameters Herrera  _et al._
(2015), lattice constant of 3.36 Å, $k_{F}$ of 0.183 $a_{0}^{-1}$ and a Dynes
parameter of 0.05 meV, which leads to a substantial smearing of the sharp
transition near to the closing of the gap. We used a Chromium-atom impurity,
with spin S = 5/2. The impurity potential scattering, $U_{i}$, and the
surface’s Rashba coupling have been neglected for Fig. 2. As the temperature
is increased, the gap of the pristine superconductor becomes smaller following
Fig. 1. However, the QPT can still be detected by the change of sign of the
gap, even for temperatures close to the critical temperature.

##  V Several impurities

The above single impurity will produce an in-gap state with particle and hole
components in the BdG formalism Yu (1965); Shiba (1968); Rusinov (1969);
Sakurai (1970). The QPT that we just described for one impurity can be
described in BdG terms as a crossing of the electron and hole characters due
to the closing of the onsite gap Sakurai (1970); Salkola  _et al._ (1997). As
the Kondo exchange interaction, $J$, of Hamiltonian, Eq. (2), increases, the
binding energy of the in-gap state increases. The BdG quasiparticles show that
the particle and hole components approach in energy. When the gap closes, the
energy of both components is the same (zero) and a crossing of the two
components takes place leading to a crossing of the ground state with the
first excited state of the system in a mean-field description Sakurai (1970).

When two impurities are present, a hybridization of the in-gap states can take
place Flatté (2000); Morr and Yoon (2006); Meng  _et al._ (2015). As a
consequence, the in-gap state split in energy and the above picture becomes
more complex. For ferromagnetically aligned states (or in general, a
superconductor with a large spin-orbit coupling) the hybridization between in-
gap states can be sizeable leading to crossings of the two electron and the
two hole components for different $J$ values. This translates into two abrupt
transition of the gap as a function of $J$. However, only one of the
transitions leads to a closing of the gap. The gap, closes for the larger $J$,
when both in-gap states have exchanged their electron-hole character. This is
clearly seen in Fig. 3. There the local density of states (LDOS) at site $i$,
$\rho_{i}(E)$, evaluated at a given quasi-particle energy $E$ is given by,

$\rho_{i}(E)=-\frac{1}{\pi}Im[G_{i,i}(E)],$ (14)

where $G_{i,i}(E)$ is the local Green’s function evaluated under one of the
impurities. The zero of energy corresponds to the Fermi energy. There are four
peaks corresponding to the two in-gap states. As before, we have used a
complex energy, $\omega+i\Gamma$ using a Dynes parameter, $\Gamma$, of 0.01
meV this time. The imaginary value leads to an effective broadening of the
density of states. As $J$ increases, the particle and hole peaks start
reducing their energy difference until they cross. Each of the crossings
correspond to one of the two transitions in Fig. 4 for the dimer curve, in
agreement with previous studies Meng  _et al._ (2015).

Figure 3: (a)-(d) shows the local density of states measured on one of the
impurities of a magnetic dimer, displaying a sequence of in-gap states as the
Kondo coupling strength, $J$, is changed to (a) 2.5 eV, (b) 3.1 eV, (c) 3.4 eV
and (d) 3.6 eV. The plots show zero-energy level crossings leading to the
appearance of two quantum phase transitions for the case of a ferromagnetic
dimer. The $J$-values at which the electron and hole parts of the shiba states
cross each other are identified with critical $J$-points at which abrupt
changes in the order parameter occurs (Fig. 4). In this calculation for
$k_{F}$ = 0.5 $a_{0}^{-1}$, the potential scattering is taken as $U$ = -5.5
eV, the spin is $S$ = 5/2, and we do not include spin-orbit coupling, $\alpha$
= 0.

Figure 4 (a) shows the result of the self-consistent order parameter as a
function of the Kondo exchange coupling, $J$, for 1, 2, 5 and 10 impurities
located at an inter-impurity distance of $3.36$ Å for an s-wave superconductor
with $k_{F}=0.5~{}a_{0}^{-1}$ and $\Delta_{0}=0.76$ meV. We retrieve the
previous results of a single QPT for a single impurity and two QPT for two
impurities ferromagnetically coupled. As the number of impurity increases,
more QPTs take place, leading to several values of $J_{c}$. The larger value
of $J_{c}$ does not seem to increase for more than 5 impurities as we can
corroborate by looking at the infinite-chain order parameter in Fig. 4 (c)
that will be analyzed in the next section.

The antiferromagnetic spin chains, Fig. 4 (b), have a very different behavior.
The transition does not seem to change its critical $J$ value, and it smooths
out for an increasing number of atoms. This behavior can be understood
considering the localized structure of the in-gap states for antiferromagnetic
chains. The hybridization of the local in-gap states to form a band is more
difficult leading to a small dispersion of the in-gap bands and to a more
uniform character of the in-gap evolution.

Figure 4: Order parameter, $\Delta$ measured on the first atom of a finite
spin-chain as a function of magnetic coupling strength, $J$, for impurities
aligned in (a) ferromagnetic and (b) antiferromagnetic ordering. The spin-
chains shown here comprise 1, 2, 5 and 10 atoms. Panels (c) and (d) show order
parameter evaluated for an infinite-chain of ferromagnetic impurities. The
abrupt changes in slope correspond to a first-order QPT. Depending on the
number of atoms in the impurity chain, more than one QPT can be realized. In
long enough chains, these abrupt changes smooth out into two QPT, one of which
is identified as the topological one (where $\Delta$ crosses the zero). All
calculations are performed without potential scattering ($U=0$) and without
spin-orbit coupling ($\alpha$=0) except in (d) where $U$ is 5.5 eV and the
Rashba coupling is $\alpha=3$ eV-Å.

###  V.1 Infinite spin chains

The infinite ferromagnetic spin chain can be calculated using Eq. (9). The
corresponding order parameters are plotted as a function of the exchange
interaction $J$ in Fig. 4 (c). Contrary to the previous finite ferromagnetic
chains, a very low number of transitions can be discerned. This has to do with
the creation of bands of in-gap states that have a very small number of gap
closings Björnson  _et al._ (2015). We have explored the order parameter as a
function of $J$, for three different values of $k_{F}$ to probe low
$k_{F}=0.183\;a_{0}^{-1}$, mid $k_{F}=0.5\;a_{0}^{-1}$ and high
$k_{F}=0.7\;a_{0}^{-1}$ electronic densities. For low and mid densities, we
find that there is a single QPT in the range of studied $J$ values. The lower
the electronic density, the lower $J_{c}$ gets, we understand this by the
increasing screening of the impurity’s magnetic moment at larger electron
densities, leading to larger values of $J_{c}$ to be able to produce an effect
on the superconductor.

At even larger electron densities ($k_{F}=0.7\;a_{0}^{-1}$), there are two
clear transitions appearing, showing that the induced in-gap band structure
matters and can give rise to complex behavior.

Figure 4 (d) shows two cases for topological QPT (TPQT) in agreement with
previous calculations Mier  _et al._ (2021a), that indeed are very similar to
the cases in (c), except for the presence of potential scattering and spin-
orbit coupling. The effect of spin-orbit coupling bridges the order-parameter
behavior between the ferro- and antiferro-magnetically ordered chains. Thus,
the transitions become smoother and closer in values of $J$. At $J=0$ the two
calculations show different $\Delta_{0}$ results because the self-consistent
procedure has been tuned to reproduce the topological results of Ref. Mier
_et al._ , 2021a at the transition value. At low electronic densities, lower
$k_{F}$, the spin-orbit coupling not only washes out the discontinuous jumps
in the order parameter but also delays the occurrence of the QPT to a higher
$J_{c}$ Björnson and Black-Schaffer (2016); Björnson  _et al._ (2017),
however, for high $k_{F}$ regimes, the Rashba coupling is inefficient in
changing the hybridization of the magnetic impurity with the electrons in the
bulk superconductor and thus the results are unaltered.

Figure 5 shows the order parameter resolved in $k$-space as a function of the
exchange interaction, $J$. For lower electron densities we find that there is
not a large variation with $k$, reflecting the very local behavior of the gap.
As a consequence, the gap closes for virtually all k-points at the same
$J_{c}$, leading to a well defined single $J_{c}$. As $J$ increases, the gap
does not seem to close again. This is consistent with the single QPT found in
Fig. 4 (c).

As the electron density increases, the gap closes at different $J_{c}$ as $k$
changes, but the $J_{c}$’s remain relatively close to each other. As a
consequence, there is basically a single effective $J_{c}$ with a small
dispersion leading to a smoother transition as compared to the lower-density
case.

At higher densities, the folding of the bands due to a Fermi wave vector
larger than the Brillouin-zone edge, $k_{BZ}$, leads to two distinct regions,
$k>k_{BZ}$ and $k<k_{BZ}$ and thus, there are two well-defined $J_{c}$ values,
seen in Figs. 4 (c) and (d) for $k_{F}=0.7\;a^{-1}_{0}$ as two abrupt changes
in the slope of the $\Delta$ vs $J$ curves. More detailed information can be
obtained from Fig. 6. This figure permits us to analyze the behavior of
$\Delta_{k}$ at high densities ($k_{F}=0.7\;a_{0}^{-1}$) as a function of $k$,
for a few $J$ values around the critical values, $J_{c}$. $\Delta_{k}$ closes
for all values of $J$ between $J=3.2$ eV, and $J=5.9$ eV, however the average
order parameter $\Delta=1/N_{k}\sum_{k}\Delta_{k}$ is zero for a unique
$J=4.6$ eV. Figures 4, 5 and 6 yield that there are two QPT for an infinite
ferromagnetic spin chain at $k_{F}=0.7\;a_{0}^{-1}$ and the average gap only
closes once.

Figure 5: Plots showing the distribution of the order parameter, $\Delta_{k}$,
in k-space as the Kondo exchange coupling strength, J is varied between 0 and
6 eV at different electronic densities represented by the Fermi wave vectors,
$k_{F}$ taking values of (a) 0.183 $a_{0}^{-1}$, (b) 0.5 $a_{0}^{-1}$ (c) 0.7
$a_{0}^{-1}$, (d) 0.3 $a_{0}^{-1}$ while the Rashba coupling and the
scattering potential are absent. In the following cases, $k_{F}$ is (e) 0.3
$a_{0}^{-1}$, and (f) 0.7 $a_{0}^{-1}$ wherein the influence of Rashba
coupling strength, $\alpha$ of 3 eV-Å and potential scattering, U of 5.5 eV is
considered. The color bar on the right denotes $\Delta_{k}$ measured in units
of meV. Figures (a), (b) and (c) corresponds to the three curves shown in
Fig.4(c) while figures (e) and (f) correspond to the two curves of Fig. 4 (d).

###  V.2 Topological quantum phase transition

The cases in Fig. 4 (d) and Fig. 5 (e) and (f) correspond to systems that
undergo TQPT as shown in Ref. Mier  _et al._ , 2021a where finite potential
scattering and Rashba coupling are operative. As the exchange coupling $J$ is
increased, Ref. Mier  _et al._ , 2021a shows that there are two TQPT creating
a topological zone between two trivial regions at larger and smaller $J$. Let
us choose two different electronic densities in the phase space of Ref. Mier
_et al._ , 2021a, namely $k_{F}=0.3\,a_{0}^{-1}$ and 0.7 $a_{0}^{-1}$. From
Fig. 4 (d) we find two $J_{c}$ delimiting the topological region in $J$ for
$k_{F}=0.7\,a_{0}^{-1}$, however for $k_{F}=0.3\,a_{0}^{-1}$ only one $J_{c}$
is clearly identified.

As before, Fig. 5 contains the information on the evolution of the order
parameter. For the low-density case ($k_{F}=0.3\,a_{0}^{-1}$) we find that
there is a small interval of $J$ values between $-k_{F}$ and $+k_{F}$ where
the gap becomes zero. The dispersion of these critical $J$ values with $k$ is
small, leading to a single transition in Fig. 4, where the dispersion leads to
a small broadening of the transition region, as we show in the preceding
section.

At higher densities ($k_{F}=0.7\,a_{0}^{-1}$), we find two transitions due to
the two different regions in $k$ space as discussed in the preceding section.
To gain insight into the behavior of the order parameter during these TQPT, we
plotted the k-resolved order parameter, $\Delta_{k}$ as a function of k for
different values of J corresponding to trivial and topological regions for
large electron-density case, $k_{F}=0.7\;a_{0}^{-1}$. The results are shown in
Fig. 6. At $J=3.2$ eV, the value of $\Delta_{k}$ is zero only at the edges of
the first Brillouin zone. As J is ramped, the system remains in this quantum
phase until a second phase transition occurs at $J=4.6$ eV. This transition is
identified when $\Delta_{k}=0$ at $k=0$. For larger $J$ values, the system
stays in the topologically trivial region and the gap does not close anymore
as can be seen in Fig. 6 for $J=6$ eV.

In brief, we can conclude that for large electronic densities, the first TQPT
takes place at large $k$ ($J=3.2$ eV) and the second TPQT occurs at a smaller
$k$ value ($J=4.6$ eV).

The topological character of these transitions is revealed by studying the in-
gap bands. Figure 7 shows the band structure and the winding number for the
above example, Ref. Mier  _et al._ , 2021a. The insets of Fig. 7(a)-(d)
clearly identifies closing and reopening of the topological band gap at large
$k$ for small $J$ and small $k$ for large $J$ as we just saw. This result
agrees with the corresponding winding number calculation (bottom set of
figures (a)-(d) and figures (f) - (j)).

The winding number is a topological invariant which is characterized by
counting the number of complete turns in anticlockwise direction that the
winding vector takes around zero. In Ref. Mier  _et al._ , 2021a a complete
account of the theory that we are using here is given. The winding vector,
$\vec{d}$ is ($d_{x},d_{y}$) that corresponds to the real and imaginary part
of the Hamiltonian in the local basis set. The vector draws a closed
trajectory in the first Brillouin zone. Depending on whether the zero is
enclosed in the vector’s trajectory, the topology will be trivial or not. The
reason after this is that taking the zero out or in the vector’s trajectory
implies to make the Hamiltonian matrix elements zero and hence to close the
band gap. Then, we can plot the evolution of ($d_{x},d_{y}$) as a function of
$k$ as done in the bottom set of Figs. 7 (a)-(d) or we can plot $d_{x}$ vs
$d_{y}$ as $k$ is changed, Figs. (f) - (j). Both sets of figures contain the
same information, but it is easy to discern when the $d$-vector really turns
around zero by comparing both sets. Furthermore, the lower-symmetry $Z_{2}$
invariant can be obtained from the values of $d_{x}$ at the Brillouin zone
center and edges, see Ref. Mier  _et al._ , 2021a.

From these figures we obtain that indeed there are two $J_{c}$, the lower one
closing the gap for $k$ values close to the Brillouin-zone edge, while the
larger one closing the gap at the center of the Brillouin zone. And that the
topological $Z_{2}$ invariant changes between the regions delimited by the
$J_{c}$. The $d_{x}$ changes sign between Brillouin zone edge and center for
$J$ values between the two $J_{c}$ showing the topological character of the
transition. The same information is confirmed by the winding number. This
shows that between the two $J_{c}$ there is a topological non-trivial phase,
while outside this region the phase is topologically trivial.

Equivalent results for the low-density case ($k_{F}=0.3\,a_{0}^{-1}$) show
that the gap only closes about $k=0$ in agreement with Fig. 5 (e), explaining
that indeed, there is only one TQPT for this case at $J_{c}=2.6$ eV. The
actual evaluation of the order parameter and its self-consistent evaluation
has permitted us to clearly show the TQPT. In Ref. Mier  _et al._ , 2021a, an
attempt to plot the regions where the gap closed by following the in-gap band
structure Pientka  _et al._ (2013) was made, becoming increasingly difficult
in certain areas of phase space. Here, we show that the self-consistent order
parameter accurately gives this information in good agreement with the
analysis of the topological invariant.

Figure 6: Order parameter as a function of wave-vectors, k for 5 different
values of J such that the system is either in trival or topological region and
at the junction of quantum phase transition. It is observed that first TQPT
corresponds to closing the gap at two k-points. This is the curve
corresponding to $J=3.2$ eV where $\Delta_{k}=0$ at the edges of the first
Brillouin zone. For the second TQPT, the order parameter is zero at a single
point $k=0$, which corresponds to the curve for $J=4.6$ eV. Here, the two
values of the critical exchange couplings, $J_{c}$, agree well with our non
self-consistent results previously calculated in Mier  _et al._ (2021a).The
parameters used for this calculation corresponds to the data presented in Fig.
5(f). Figure 7: Energy band diagrams ((a)-(e) top) and evolution of winding
vector $\vec{d}$ = ($d_{x},d_{y}$) ((a)-(e) bottom) as a function of wave-
vectors, k . The inset of the top set of figures simply shows a zoomed in cut-
out of the bands near zero energy. From left to right, J = 2.0, 3.2, 4.0, 4.6
and 5.0 eV respectively. (f)-(j) illustrates the trajectory of normalized
winding vector, $\vec{d}(k)$ as k is varied in the first Brillouin zone. The
corresponding plots showcase the evolution of topological quantum phase
transition as J is tuned through a range of values. The first TQPT is marked
by a zero band gap at the boundary of the Brillouin zone where the effect of
band folding is observed. This is the onset of the topological region which is
verified by a completion of one full turn of the winding vector, shown in (g).
The second TQPT is observed when topological band closes at k=0. It can be
observed that the winding vector already ceases to be 1 (g) and this leads to
transition to the trivial region where the winding number goes back to 0 (j).
A potential scattering of U = 5.5 eV, Rashba spin-orbit coupling of $\alpha$
=3 eV-Å and fermi wave-vector, $k_{F}=0.7a_{0}^{-1}$ is used for this
calculation.

##  VI Effects of the self-consistent procedure on the local density of states

Closing the local gap should have important effects on the LDOS and hence on
the result of STM studies of in-gap structure. Our calculations show that the
LDOS before and after self-consistency are identical. Hence, the conclusions
of the present study are that non-self-consistent calculations are enough to
describe most experiments based on scanning tunneling microscope studies of
the induced in-gap structure. Indeed, self-consistency has impact on the
evaluation of the order parameter of superconducting gap, but it has no
bearings on the actual crossings of in-gap states and hence on the in-gap
bands. As a consequence the topology of BdG equations is the same with or
without self-consistency.

All the above results have been obtained for wide-band superconductors, where
the bandwidth of the metal phase is orders of magnitude larger than the
superconducting gap as is the case for most s-wave superconductors. As a
consequence, $\lambda_{F}$ is much shorter than the coherence length $\xi$.
Since $\lambda_{F}$ is the scale of distances where the self-consistent gap
changes, and $\xi$ is the usual superconducting length scale, the self-
consistency does not change most of the superconductor, and does not affect
its main properties. Examples are the above in-gap states, leading to DOS and
in-gap band structures that are largerly unaffected during self-consistency.

Previous results in the literature Björnson  _et al._ (2015, 2017); Theiler
_et al._ (2019) were computed in the case when $\lambda_{F}\sim\xi$ that
corresponds to superconducting gaps in the order of magnitude of the normal-
metal bandwidth. We could not reach this very low electron densities, but the
results from these works show that self-consistency qualitatively alter the
in-gap structure and the topological character of the induced bands.

##  VII Discussion and conclusions

We have studied the self-consistent order parameter in the context of magnetic
impurities in realistic bandwidth superconductors. The self consistency is
computed by introducing the variation of the order parameter as a self-energy
in the Dyson equation and calculating the complete Green’s function of the
system. This allows us to evaluate the new order parameter from the Green’s
function and its variation with respect to the previous iteration. This
procedure is iterated until the variations in the order parameter are
negligible. We have reproduced well-known BCS results such as the temperature
behavior of the order parameter and rationalized the induced in-gap structure
by a single impurity, including the spatial extension of the in-gap states and
the temperature behavior of the QPT as a function of the impurity-substrate
coupling.

As we increase the number of impurities on the substrate, the order-parameter
discontinuities give precise information on the appearance of QPT associated
with the crossings of the chemical potential of the different in-gap states as
the impurity-substrate magnetic interaction increases. For ferromagnetic spin
chains of increasing number of impurities the discontinuities in the order
parameter with the exchange interaction seem to yield a maximum number of
QPTs. A careful study in $k$-space show that even if there is a large number
of chemical potential crossings, the behavior of the QPTs is rather given by
the order-parameter as a function of $k$-point. Thus, the order parameter
contains very relevant information on the presence of QPTs.

Antiferromagnetic chains present a single QPT in short chains due to the
localized character of the in-gap states and of the induced bands. When spin-
orbit coupling is introduced, the situation becomes less clear and the
transitions become very smooth, trending to a mixed behavior of the order
parameter for ferromagnetic and antiferromagnetic spin chains.

In the presence of spin-orbit coupling, we can study the topological phases of
the spin chain through the order parameter. In addition to understanding the
phase space where QPTs are produced, studying the topological invariants at
the same time as the self-consistent order parameter gives clear onsets of the
different topological phases. Hence, the self-consistent order parameter
becomes an interesting object to study the phase space for QPTs and
particularly for TQPTs.

For the present systems, we find fast convergence and virtually no impact of
self-consistency on the evaluation of the density of states and on the
topological character of the induced in-gap bands. This is a consequence of
the short-range changes in the order parameter that die within Fermi-length
scale, $\lambda_{F}$ from the impurity. However, the superconducting
properties are set within the coherence length, $\xi$, which is typically
orders of magnitude larger in realistic s-wave superconductors. Thus, non-
self-consistent calculations in large-bandwidth superconductors give accurate
results on the induced in-gap structure contrary to the case of small-bandwith
superconductors Björnson  _et al._ (2015, 2017); Theiler  _et al._ (2019).

In summary, computing the order-parameter is a small computational overhead
that gives valuable information on QPT, temperature dependence of the
superconducting properties and the topological phase space of the in-gap
structure.

###### Acknowledgements.

We are pleased to thank Dr. Cristina Mier for discussions and assistance in
the early stages of this work. Special thanks to Prof. Michael Flatté for
insightful discussions on the topic. This work was supported by the project
PID2021-127917NB-I00 funded by MCIN/AEI/10.13039/501100011033, QUAN-000021-01
funded by the Gipuzkoa Provincial Council and IT-1527-22 funded by the Basque
Government.

##  Appendix A Calculation of the e-h coupling constant

The calculation of Eq. (8) is difficult because despite the natural cutoff of
the Debye frequency, $\omega_{D}$, the superconducting gap, $\Delta$ is still
one or two orders of magnitude smaller than $\omega_{D}$, for typical
realistic-bandwidth BCS superconductors. Here, we assume a second cutoff,
$\omega_{c}$ that is smaller than $\omega_{D}$, but large enough to avoid
problems in the self-consistency. If the cutoff is too small, the self-
consistent gap shows an increase in value as $J$ increases, which disappears
if the cutoff is increased, Fig. 8.

The procedure to attain self-consistency Flatté and Byers (1997, 1997) is to
include $\delta\Delta$ as a self-energy in Dyson’s equation, Eq. (4), where
$\delta\Delta$ is the change in the gap between self-consistency iterations.
The integral of Eq. (8) between both cutoffs, $\omega_{c}$ and $\omega_{D}$ is
assumed to be a constant, because for large-enough energy, we expect to
recover the normal-metal electronic structure. The advantage is that we do not
need to evaluate the integrals between cutoffs in $\delta\Delta$.

Figure 8: The self-consistent gap as a function of the impurity’s Kondo
exchange interaction for four different cutoffs
$\hbar\omega_{c}=1.5\Delta_{0},3\Delta_{0}$, $10\Delta_{0}$ and
$20\Delta_{0}$, where $\Delta_{0}=\Delta(J=0)=0.76$ meV. A low cutoff leads to
a qualitatively wrong behavior of the order parameter, although the critical
Kondo exchange interaction, $J_{c}$, does not change.

In order to evaluate $\tilde{V}$ of Eq. (8), we use that
$\omega_{c}\gg{\Delta}_{0}$ such that

$Im(G^{1,4}_{i,i}(\omega)-G^{2,3}_{i.i}(\omega))=-Re(\frac{2\pi
N_{0}\Delta_{0}}{\sqrt{(\omega+i\gamma)^{2}-\Delta_{0}^{2}}})$ (15)

for temperatures small enough to be in the superconducting phase
($k_{B}T\ll\hbar\omega_{c}$) and for small Dynes parameters
($\gamma\ll\Delta$), we have that the correction to the gap because the cutoff
is actually $\omega_{D}$ is

$\displaystyle\delta\Delta_{i}$ $\displaystyle=$
$\displaystyle\tilde{V}\int^{-\hbar\omega_{c}}_{-\hbar\omega_{D}}f(\omega)Im\\{G^{1,4}_{i,i}(\omega)-G^{2,3}_{i.i}(\omega)\\}d\omega$
(16) $\displaystyle\approx$ $\displaystyle-\tilde{V}\;2\pi
N_{0}{\Delta}_{0}\;ln\left(\frac{\omega_{c}}{\omega_{D}}\right).$

Given the qualitative nature of the present study, we have ignored the
correction and assumed that is partially taken care of by the actual value of
$\tilde{V}$, because we obtain $\tilde{V}$ from Eq. (8), for the clean-
superconductor gap. However, the actual critical values of the QPT will depend
on the values taken for $\omega_{c}$ and $\omega_{D}$. All calculations have
been performed for $\omega_{c}=10\,{\Delta}_{0}$, where the behavior of
${\Delta}_{0}$ is consistent with the results of previous studies Salkola  _et
al._ (1997); Flatté and Byers (1997); Karan  _et al._ (2022) as can be seen in
Fig. 8.

##  Appendix B Role of non-magnetic potential scattering

According to Anderson’s theorem, non-magnetic impurities do not alter s-wave
superconductivity Anderson (1959). However, a non-magnetic impurity does have
an effect although it cannot produce in-gap states Soda  _et al._ (1967).
Within our theory, it is clear that the magnetic part is not the only one to
contribute to the modification of the order parameter. Indeed, the non-
magnetic coupling can have a significant contribution. Figure 9 shows that
there is a reduction of about 48% in the superconducting order parameter as
the coupling strength is tuned from 0 up to 6 eV. Our results are in good
agreement with the observations previously reported by Flatte & Byers in Ref.
Flatté and Byers, 1999.

Figure 9: (a) Plot showing impact of non-magnetic coupling strength, U on the
superconducting order parameter, $\Delta$. (b) 2d map shows the distribution
of $\Delta_{k}$ in k-space as U is ramped from 0 to 6 eV. Color bar
corresponds to $\Delta$ given in meV units. The Rashba spin-orbit coupling is
fixed to 3 eV-Å, $k_{F}$ = 0.3 $a_{0}^{-1}$ and Kondo exchange coupling, $J$,
is absent in this calculation.

## References

  * Balatsky  _et al._ (2006) A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Rev. Mod. Phys. 78, 373 (2006).
  * Zhu (2016) J.-X. Zhu,  _Bogoliubov-de Gennes Method and its Applications_ , Lecture Notes in Physics (Springer International Publishing, 2016).
  * Yu (1965) L. Yu, Acta Physica Sinica 21, 75 (1965).
  * Shiba (1968) H. Shiba, Progress of Theoretical Physics 40, 435 (1968), https://academic.oup.com/ptp/article-pdf/40/3/435/5185550/40-3-435.pdf .
  * Rusinov (1969) A. I. Rusinov, Soviet Journal of Experimental and Theoretical Physics Letters 9, 85 (1969).
  * Schrieffer and Wolff (1966) J. R. Schrieffer and P. A. Wolff, Phys. Rev. 149, 491 (1966).
  * Sakurai (1970) A. Sakurai, Prog of Theo Phys 44, 1472 (1970).
  * Salkola  _et al._ (1997) M. I. Salkola, A. V. Balatsky, and J. R. Schrieffer, Phys. Rev. B 55, 12648 (1997).
  * Nagai  _et al._ (2013) Y. Nagai, Y. Shinohara, Y. Futamura, Y. Ota, and T. Sakurai, arXiv: Superconductivity (2013).
  * Flatté and Byers (1997) M. E. Flatté and J. M. Byers, Phys. Rev. Lett. 78, 3761 (1997), https://link.aps.org/doi/10.1103/PhysRevLett.78.3761 .
  * Flatté and Byers (1997) M. E. Flatté and J. M. Byers, Phys. Rev. B 56, 11213 (1997).
  * Flatté (2000) M. E. Flatté, Phys. Rev. B 61, R14920 (2000).
  * Hoffman  _et al._ (2015) S. Hoffman, J. Klinovaja, T. Meng, and D. Loss, Phys. Rev. B 92, 125422 (2015).
  * Björnson  _et al._ (2017) K. Björnson, A. V. Balatsky, and A. M. Black-Schaffer, Phys. Rev. B 95, 104521 (2017).
  * Theiler  _et al._ (2019) A. Theiler, K. Björnson, and A. M. Black-Schaffer, Phys. Rev. B 100, 214504 (2019).
  * Villas  _et al._ (2020) A. Villas, R. L. Klees, H. Huang, C. R. Ast, G. Rastelli, W. Belzig, and J. C. Cuevas, Phys. Rev. B 101, 235445 (2020).
  * Hatter  _et al._ (2015) N. Hatter, B. W. Heinrich, M. Ruby, J. I. Pascual, and K. J. Franke, Nature Communications 6 (2015).
  * Farinacci  _et al._ (2018) L. Farinacci, G. Ahmadi, G. Reecht, M. Ruby, N. Bogdanoff, O. Peters, B. W. Heinrich, F. von Oppen, and K. J. Franke, Phys. Rev. Lett. 121, 196803 (2018).
  * Huang  _et al._ (2020) H. Huang, R. Drost, J. Senkpiel, C. Padurariu, B. Kubala, A. L. Yeyati, J. C. Cuevas, J. Ankerhold, K. Kern, and C. R. Ast, Communications Physics 3, 199 (2020).
  * Liebhaber  _et al._ (2021) E. Liebhaber, L. M. Rütten, G. Reecht, J. F. Steiner, S. Rohlf, K. Rossnagel, F. von Oppen, and K. J. Franke, Nature Communications 13 (2021).
  * Karan  _et al._ (2022) S. Karan, H. Huang, C. Padurariu, B. Kubala, A. Theiler, A. M. Black-Schaffer, G. Morrás, A. L. Yeyati, J. C. Cuevas, J. Ankerhold, K. Kern, and C. R. Ast, Nature Physics 18, 893 (2022).
  * Liu  _et al._ (2023) Y. Liu, C. Li, F.-H. Xue, W. Su, Y. Wang, H. Huang, H. Yang, J. Chen, D. Guan, Y. Li, H. Zheng, C. Liu, M. Qin, X. Wang, R. Wang, L. Deng-Yuan, P.-N. Liu, S. Wang, and J. Jia, Nano Letters (2023), 10.1021/acs.nanolett.3c02208.
  * Zhou  _et al._ (2022) Y. Zhou, J. Guo, S. Cai, J. Zhao, G. Gu, C. Lin, H. Yan, C. Huang, C. Yang, S. Long, _et al._ , Nature Physics 18, 406 (2022).
  * Uldemolins  _et al._ (2023) M. Uldemolins, A. Mesaros, G. D. Gu, A. Palacio-Morales, M. Aprili, P. Simon, and F. Massee, (2023), arXiv:2310.06030 [cond-mat.supr-con] .
  * Franke  _et al._ (2011) K. J. Franke, G. Schulze, and J. I. Pascual, Science 332, 940 (2011).
  * Ji  _et al._ (2010) S.-H. Ji, Y.-S. Fu, T. Zhang, X. Chen, J.-F. Jia, Q.-K. Xue, and X.-C. Ma, Chinese Physics Letters 27, 087202 (2010).
  * Maurand  _et al._ (2012) R. Maurand, T. Meng, E. Bonet, S. Florens, L. Marty, and W. Wernsdorfer, Phys. Rev. X 2, 011009 (2012).
  * Delagrange  _et al._ (2015) R. Delagrange, D. J. Luitz, R. Weil, A. Kasumov, V. Meden, H. Bouchiat, and R. Deblock, Phys. Rev. B 91, 241401 (2015).
  * Delagrange  _et al._ (2016) R. Delagrange, R. Weil, A. Kasumov, M. Ferrier, H. Bouchiat, and R. Deblock, Phys. Rev. B 93, 195437 (2016).
  * Deacon  _et al._ (2010) R. S. Deacon, Y. Tanaka, A. Oiwa, R. Sakano, K. Yoshida, K. Shibata, K. Hirakawa, and S. Tarucha, Phys. Rev. Lett. 104, 076805 (2010).
  * von Oppen and Franke (2021) F. von Oppen and K. J. Franke, Phys. Rev. B 103, 205424 (2021).
  * Schecter  _et al._ (2016) M. Schecter, K. Flensberg, M. H. Christensen, B. M. Andersen, and J. Paaske, Phys. Rev. B 93, 140503 (2016).
  * Žonda  _et al._ (2015) M. Žonda, V. Pokorný, V. Janiš, and T. Novotný, Scientific Reports 5, 8821 (2015).
  * Rozhkov and Arovas (1999) A. V. Rozhkov and D. P. Arovas, Phys. Rev. Lett. 82, 2788 (1999).
  * Levy Yeyati  _et al._ (1995) A. Levy Yeyati, A. Martín-Rodero, and F. J. García-Vidal, Phys. Rev. B 51, 3743 (1995).
  * Heinrich  _et al._ (2018) B. W. Heinrich, J. I. Pascual, and K. J. Franke, Progress in Surface Science 93, 1 (2018).
  * Choy  _et al._ (2011) T.-P. Choy, J. M. Edge, A. R. Akhmerov, and C. W. J. Beenakker, Physical Review B 84, 195442 (2011).
  * Pientka  _et al._ (2013) F. Pientka, L. I. Glazman, and F. von Oppen, “Topological superconducting phase in helical shiba chains,”  (2013).
  * Choi  _et al._ (2019) D.-J. Choi, N. Lorente, J. Wiebe, K. von Bergmann, A. F. Otte, and A. J. Heinrich, Reviews of Modern Physics 91, 041001 (2019).
  * Björnson  _et al._ (2015) K. Björnson, S. S. Pershoguba, A. V. Balatsky, and A. M. Black-Schaffer, Phys. Rev. B 92, 214501 (2015).
  * Christensen  _et al._ (2016) M. H. Christensen, M. Schecter, K. Flensberg, B. M. Andersen, and J. Paaske, Phys. Rev. B 94, 144509 (2016).
  * Mier  _et al._ (2021a) C. Mier, D.-J. Choi, and N. Lorente, Phys. Rev. B 104, 245415 (2021a).
  * Mier  _et al._ (2022) C. Mier, D.-J. Choi, and N. Lorente, Phys. Rev. Res. 4, L032010 (2022).
  * Schneider  _et al._ (2022) L. Schneider, P. Beck, J. Neuhaus-Steinmetz, L. Rózsa, T. Posske, J. Wiebe, and R. Wiesendanger, Nature Nanotechnology , 1 (2022).
  * Mier  _et al._ (2021b) C. Mier, J. Hwang, J. Kim, Y. Bae, F. Nabeshima, Y. Imai, A. Maeda, N. Lorente, A. Heinrich, and D.-J. Choi, Phys. Rev. B 104, 045406 (2021b).
  * Cuevas  _et al._ (1996) J. C. Cuevas, A. Martín-Rodero, and A. L. Yeyati, Physical Review B 54, 7366 (1996).
  * Vernier  _et al._ (2011) E. Vernier, D. Pekker, M. W. Zwierlein, and E. Demler, Phys. Rev. A 83, 033619 (2011).
  * Meng  _et al._ (2015) T. Meng, J. Klinovaja, S. Hoffman, S. Pascal, and D. Loss, Phys. Rev. B 92, 064503 (2015), https://link.aps.org/doi/10.1103/PhysRevB.92.064503 .
  * Peng  _et al._ (2015) Y. Peng, F. Pientka, L. I. Glazman, and F. von Oppen, Physical Review Letters 114, 106801 (2015).
  * Brydon  _et al._ (2015) P. M. R. Brydon, S. Das Sarma, H.-Y. Hui, and J. D. Sau, Phys. Rev. B 91, 064505 (2015).
  * Ménard  _et al._ (2015) G. C. Ménard, S. Guissart, C. Brun, S. Pons, V. S. Stolyarov, F. Debontridder, M. V. Leclerc, E. Janod, L. Cario, D. Roditchev, P. Simon, and T. Cren, Nature Physics 11, 1013 (2015).
  * Schneider  _et al._ (2021) L. Schneider, P. Beck, T. Posske, D. Crawford, E. Mascot, S. Rachel, R. Wiesendanger, and J. Wiebe, Nature Physics 17, 943 (2021).
  * Matsuura (1977) T. Matsuura, Prog. Theor. Phys. 57, 1823 (1977).
  * žitko  _et al._ (2011) R. žitko, O. Bodensiek, and T. Pruschke, Phys. Rev. B 83, 054512 (2011).
  * Schmid  _et al._ (2022) H. Schmid, J. F. Steiner, K. J. Franke, and F. von Oppen, Phys. Rev. B 105, 235406 (2022).
  * Moca  _et al._ (2008) C. Moca, E. Demler, B. Jankó, and G. Zaránd, Phys. Rev. B 77, 174516 (2008).
  * Imai  _et al._ (2012) Y. Imai, F. Nabeshima, T. Yoshinaka, K. Miyatani, R. Kondo, S. Komiya, I. Tsukada, and A. Maeda, Journal of the Physical Society of Japan 81, 113708 (2012), https://doi.org/10.1143/JPSJ.81.113708 .
  * Herrera  _et al._ (2015) E. Herrera, I. Guillamón, J. A. Galvis, A. Correa, A. Fente, R. F. Luccas, F. J. Mompean, M. Garcia-Hernandez, S. Vieira, J. P. Brison, and H. Suderow, Physical Review B 92, 054507 (2015).
  * Dynes  _et al._ (1978) R. C. Dynes, V. Narayanamurti, and J. P. Garno, Phys. Rev. Lett. 41, 1509 (1978).
  * Ruby  _et al._ (2016) M. Ruby, Y. Peng, F. von Oppen, B. W. Heinrich, and K. J. Franke, Phys. Rev. Lett. 117, 186801 (2016).
  * Choi  _et al._ (2017) D.-J. Choi, C. Rubio-Verdú, J. de Bruijckere, M. M. Ugeda, N. Lorente, and J. I. Pascual, Nature Communications 8, 15175 (2017).
  * Yazdani  _et al._ (1997) A. Yazdani, B. A. Jones, C. P. Lutz, M. F. Crommie, and D. M. Eigler, Science 275, 1767 (1997).
  * Morr and Yoon (2006) D. K. Morr and J. Yoon, Phys. Rev. B 73, 224511 (2006).
  * Björnson and Black-Schaffer (2016) K. Björnson and A. M. Black-Schaffer, Phys. Rev. B 94, 100501 (2016).
  * Anderson (1959) P. Anderson, Journal of Physics and Chemistry of Solids 11, 26 (1959).
  * Soda  _et al._ (1967) T. Soda, T. Matsuura, and Y. Nagaoka, Prog. Theor. Phys. 38, 551 (1967).
  * Flatté and Byers (1999) M. E. Flatté and J. M. Byers (Academic Press, 1999) pp. 137–228.