Anisotropic ferroelectric distortion effects on the RKKY interaction in topological crystalline insulators

Manipulation of electronic and magnetic properties of topological materials is a topic of much interest in spintronic and valleytronic applications. Perturbation tuning of multiple Dirac cones on the (001) surface of topological crystalline insulators (TCIs) is also a related topic of growing interest. Here we show the numerical evidence for the ferroelectric structural distortion effects on the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between two magnetic impurity moments on the SnTe (001) and related alloys. The mirror symmetry breaking between Dirac cones induced by the ferroelectric distortion could be divided into various possible configurations including the isotropically gapped, coexistence of gapless and gapped, and anisotropically gapped phases. Based on the retarded perturbed Green’s functions of the generalized gapped Dirac model, we numerically find the RKKY response for each phase. The distortion-induced symmetry breaking constitutes complex and interesting magnetic responses between magnetic moments compared to the pristine TCIs. In the specific case of coexisted gapless and gapped phases, a nontrivial behavior of the RKKY interaction is observed, which has not been seen in other Dirac materials up until now. For two impurities resided on the same sublattices, depending on the distortion strength, magnetic orders above of a critical impurity separation exhibit irregular ferromagnetic ⇔ antiferromagnetic phase transitions. However, independent of the impurity separation and distortion strength, no phase transition emerges for two impurities resided on different sublattices. This essential study sheds light on magnetic properties of Dirac materials with anisotropic mass terms and also makes TCIs applications relatively easy to understand.


Scientific Reports
| (2021) 11:5273 | https://doi.org/10.1038/s41598-021-84398-0 www.nature.com/scientificreports/ the surface states of SnTe in the presence of anisotropic ferroelectric distortion in the present work can be readily detected in ARPES and tunneling spectroscopy experiments. The rest of the paper is organized as follows. In "Theoretical background" section, we first introduce the theoretical background of pristine TCIs. Then we present the electronic features of TCIs in the presence of ferroelectric distortion. In "Low-energy RKKY coupling in TCIs in the presence of ferroelectric distortion" section, we calculate the RKKY interaction in gapped TCIs and show the numerical results when two magnetic impurities are on the same and different sublattices. Finally, the paper ends in "Summary" section with a summary of the remarkable findings.

Theoretical background
Electronic features in unperturbed TCIs: from four-band to two-band model. We begin by representing the reported angle-resolved-photoemission-spectroscopy experiments 8,11,35 of the surface states in (001) plane of SnTe and related alloys such as Pb 1−x Sn x Te and Pb 1−x Sn x Se. The "unperturbed" effective Hamiltonians of four Dirac cones at { x and ′ x } near the X 1 point and { y and ′ y } near the X 2 point of SBZ forming surface states in TCIs are given by [we set = 1 for simplicity throughout the present work] where the three terms describe respectively two copies of Dirac cones, the term which shifts the energy of these two Dirac cones to form two new high-energy Dirac cones with opposite energies and the term which shifts two copies of Dirac cones to have opposite chiralities. The Pauli matrices � σ = (σ 0 , σ x , σ y , σ z ) and � τ = (τ 0 , τ x , τ y , τ z ) respectively act in the spin and space. Researches were first dealt only with the first term, while the atomically sharp interface between the system and the vacuum results in two additional term, so-called intervalley scattering terms 39 . Although extra momentum k = (k x , k y )-dependent terms have been added to these models later, they do not affect the experimentally observed physics of TCIs significantly 41,42 . The used parameters in this paper η 1 = 3.53 eV.Å , η 2 = 1.91 eV.Å , n = 0.055 eV, and δ = 0.04 eV are taken from Refs. 38,40 .
TCI (001) plane of SnTe or related alloys possesses C 4 rotation symmetries along the x and y directions, implying that one may focus on X 1 only and use the relation X 2 = C 4 X 1 C −1 4 to capture the physics around the X 2 point. Of other symmetries, C 2 is important as well between Dirac cones themselves so that ′ x = C 2 x C −1 2 and ′ y = C 2 y C −1 2 . Of course, we also allow to conclude y = C 4 x C −1 4 and ′ y = C 4 ′ x C −1 4 [38][39][40] . Figure 1 captures the configurations of four Dirac cones in the SBZ of (001) plane in TCIs. The symmetry breaking perturbations do not enter Eq.(1) and will be introduced at the level of low-energy two-band model.
Focusing the above points, the energy spectrum of surface states in the vicinity of X 1 point from Eq. (1a) is given by where N 2 = n 2 + δ 2 , µ = +(−) refers to the conduction (valence) band and +(−) inside the square root stands for the upper (lower) conduction and lower (upper) valence band. We print the band structure along the lines Ŵ − X 1 and X 1 − M in the SBZ [see Fig. 1] in Fig.2a. The system is gapless at Dirac point � ′ x = (−N/η 1 , 0) , while the gap opens at saddle point S 1 = (0, +n/η 2 ) . We stress that using the rotational symmetries C 2 and C 4 , one simply obtains other Dirac and saddle points.
A physical quantity particularly useful for such dense systems is the density of states (DOS) in which the electronic correlations between charge carriers is translated to the Green's functions using the identity where o + = 5 meV is a very small real number. This Green's function is related to the DOS via (1a) H X 1 (k) = η 1 k x σ y − η 2 k y σ x ⊗ τ 0 + nτ x + δσ y τ y , (1b) H X 2 (k) = η 2 k x σ y − η 1 k y σ x ⊗ τ 0 + nτ x + δσ x τ y , www.nature.com/scientificreports/ The specific nature of Green's function elements will be the subject of the RKKY interaction of perturbed TCIs in "Low-energy RKKY coupling in TCIs in the presence of ferroelectric distortion" section. Having the diagonal elements of the Green's function matrix, one obtains DOS, as shown in Fig. 2b. From Eq. (3), DOS diverges at singularities of the numerator, occurring at energies ±δ corresponding to the saddle points, leading to emerged van Hove singularities. These singularities remind a Lifshitz transition [a deformation of the Fermi surface] 38,40 . This procedure will be altered when the system is perturbed.
Here we adopt one simplification with a view for applying and/or calculating the ferroelectric distortioninduced RKKY coupling in TCIs. We immediately turn to the two-band model by linearizing the Hamiltonian near the Dirac cones by starting with the x point, since the gapless phase of surface states originates from the Dirac cones. Thus, we leave the saddle points here since it is gapped naturally. We achieve the following lowenergy Hamiltonian 43 To this end, p x = k x − x measuring from the x point, p y = k y , η 1 = (2N/n)η 1 and η 2 = (2δ/n)η 2 are required to be produced. Eventually, the energy spectrum is given by describing the linear dispersion energy of massless fermions along the line Ŵ − X 1 in the (001) plane of TCIs [see left side of Fig. 2b]. We adopt one more simplification here with the aid of definitions v F = η 1η2 , K x = η 1 /η 2 p x and K y = η 2 /η 1 p y . Hence, we have We will follow Eq. (7a) hereafter as the effective low-energy Hamiltonian, describing the basic physics of gapless surface states. Using the transformations {σ i � → −σ i and K i → −K i } [i ∈ {x, y} ] for two-fold C 2 and {σ y � → −σ x , σ x � → σ y , K y � → −K x , K x � → K y } for four-fold C 4 rotational symmetries for the pristine TCIs, one may list the Hamiltonian of other Dirac cones in the SBZ as H 0 Let us apply the ferroelectric structural distortion which acts as a perturbation in the system. Electronic features in TCIs in the presence of ferroelectric distortion. In this part, we first express the two-band Hamiltonian for TCI (001) surface states in the presence of ferroelectric distortion. Because the most essential physics of susceptibility is linked to the DOS, we apply our analysis to establish the ferroelectric distortion-induced electronic features in TCIs.
Two-band Hamiltonian of ferroelectric distortion effects is simply given by a mass term originating from the atomic displacements. This displacement may be produced by strain or electric field, resulting in crystal  Adding these terms to the corresponding pristine Hamiltonians, one may find the gapped phase for surface states. Again, we focus on the x here inspiring the presence of two-and four-fold rotational symmetries in pristine TCI. By this, we obtain Using the aforementioned C 2 and C 4 rotational symmetries, the Hamiltonian of other gapped Dirac cones in the SBZ can be listed as 36,40 By evaluating Eq. (9) with the mass term xF , one obtains the gapped dispersing bands at Dirac cones, as represented in Fig. 3 [the bands appear similarly with different slopes if one plots E ( K) as a function of K y ]. One, in the left panel of Fig. 3, observes that the bands get apart from each other to open a gap of 2 xF keeping the Fermi energy at zero energy. It should be pointed out that for the two-band model, van Hove singularities are absent in DOS due to the linear dispersion energy of Dirac fermions as well as due to throwing away saddle points as the origin of the DOS divergences. As shown in the right panel of Fig. 3, the electronic DOS presents a different trend for surface states. It should be mentioned that we have plotted the  www.nature.com/scientificreports/ total DOS here by summing over the Green's functions from both X 1 and X 2 points, i.e. Eq. (4) is rewritten as D (E ) = − 2 π k∈SBZ Im Tr G 0 X 1 (k, E ) + Tr G 0 X 2 (k, E ) , where Eq. (3) can be employed to find the Green's functions around X 1 and X 2 points. As can be seen, the final band of the system is gapless, giving rise to the relevance of RKKY compared to the van Vleck mechanism. Out of the gap, despite of a tiny difference close to the band edges, the value of DOS is almost the same for the pristine TCI compared to the special case, coexistence of gapped and gapless Dirac cones. Similarly, one finds the band structure other Dirac cones, however; the gap sizes may be different applying different xF and yF simultaneously. Simultaneous displacement means that the displacement is not applied purely along the x-or y-direction. This, in turn, leads to valley-dependent electronic features on the (001) surface of TCIs. Since the distortion does not add/remove particles to/from the system, the number of states should be preserved, meaning that the area under DOS is conserved in the absence and presence of ferroelectric distortion. However, to apply this constraint on DOS, we need to work with full band spectrum of SnTe crystal. In this work, Hamiltonian is derived based on an effective theory at low energies, and without loss of any behavior of susceptibility, we leave this constraint in Fig. 3.
Nevertheless, in the case of gapped TCI, the relevance of RKKY response is questionable by definition since there is no any states inside the gap around the Fermi energy. Let us comment on this point before delving into the RKKY results. Of the other well-known interactions for magnetic coupling is given by the van Vleck mechanism. It is demonstrated that for chemical potentials inside the gap, the van Vleck interaction between magnetic impurities located on the surface of a thin film of topological insulator is very large and relevant which is stemming from the strong spin-orbit coupling observed in Bi 2 Te 3 family 22 . This large interaction is responsible for strong ferromagnetism observed in TI thin films 22 . The origin of such ferromagnetic ordering comes back to the large coupling matrix elements of S z between the conduction and valence band states which is sizable in the band inverted phase. It is also demonstrated that the van Vleck interaction is controlled by the structural inversion asymmetry 44 . However, two-band effective Hamiltonian represented in Eq. (7a) is written based on the p-orbitals of Sn and Te atoms 40 without any mixing. Moreover, the gap opened here is originated from the breaking of mirror symmetry induced by the ferroelectric distortion when spin-orbit coupling is ruled out at low energies. As a result, the van Vleck interaction in gapped TCIs is quite small. In such situation, although RKKY coupling is quite small too, RKKY interaction is a relevant quantity. Indeed, as shown in the integration of Eq. (13), valence band states are involved to the response of magnetic impurities even if the Fermi level lies inside the gap. So the value of susceptibility would be very smaller than the case in which the Fermi level falls inside the conduction band. Based on a self-consistent calculation on the surface spectrum of TIs, the structure of RKKY interaction is preserved in the presence of the gap induced by exchange field of magnetic impurities, nevertheless, its strength is considerably suppressed 45 . It should be mentioned that the RKKY interaction is also calculated inside the gap energies of TI induced by proximity of a superconductor on top of TI surface 46 . Furthermore, RKKY interaction is relevant in the presence of local gaps induced by doped magnetic impurities on the surface of TI 47 . We will explain the fingerprints of the gapped Dirac ones in our system in RKKY coupling in the next section.

Low-energy RKKY coupling in TCIs in the presence of ferroelectric distortion
To obtain the coupling J between two magnetic moments with the separation of R = R 2 − R 1 considering different lattice sites {α, β} = 1, 2 , one allows to use the second-order perturbation theory resulting in where J RKKY αβ (R) is called the RKKY exchange interaction strength given by in which one of the impurities is set to the zero position and another one at R . In this equation, χ αβ (0, R) is the spin susceptibility which is directly proportional to the RKKY exchange interaction J RKKY αβ (R). The spin susceptibility for our spin-degenerate system described by a two-band model can be extracted using the retarded Green's functions 48-51 as where i, j ∈ {x, y, z} , E F is the Fermi energy and wherein G 0 αβ (E , R, 0) is the 2 × 2 matrix of lattice non-interacting Green's functions in the spin space. Regarding the Fermi energy, Fig. 3 states that it lies in the zero energy, in both gapless and gapped [middle of band gap] phases. To study the RKKY responses when the system is not doped by donors or acceptors, it is our task to set E F to zero in the following. Taking the Fourier transform of the reciprocal-space Green's functions into account in the vicinity of two in-equivalent X 1 and X 2 points in the SBZ of (001) plane in TCIs, one obtains where q is the momenta around the Dirac cones. These momenta should be small enough so that each Dirac cone involves in the interaction, i.e. |� q| ≪ | � i | and | � ′ i | ( i ∈ {x, y} ). Also, the integration is performed over the entire with the defined R x = R cos(ϕ R ) and R y = R sin(ϕ R ) . All G functions are the real-space Green's function, however, the reciprocal-space Green's function is easily calculated through w h e r e E + io + = iω a n d It is necessary to point out that we are allowed to extend the momentum cutoff q c to ∞ since the RKKY response at long (short) distance between magnetic moments arises mainly from the small (large) momenta. Thus, to cover the most contribution of momenta over the entire SBZ, we use q c → ∞ simply.
After pretty simple calculations, we find where A = ω 2 + � 2 xF and K 0/1 (A R/v F ) is the modified Bessel function of the zero/first kind. Note that these only provide the lattice Green's functions for x point. It is straightforward to deduce that where B = ω 2 + � 2 yF . By inspiration of the two-fold C 2 and four-fold C 4 rotational symmetries, one simply finds lattice Green's functions for other three Dirac points:  (14) gives rise to C αβ ii (ω, R, 0) = G 0 αβ (ω, R, 0)G 0 αβ (ω, 0, R) and C αβ ij (ω, R, 0) = 0 which i and j could be each of x, y and z due to the spin-degeneracy. Eventually, using the feature dE = idω from the definition of E + io + = iω , the spin susceptibility reads Having Green's functions as well as susceptibility we divide the following discussion into two parts: (i) impurities located on the same sublattices and (ii) impurities located on different sublattices. We comment that it is of course possible to locate impurities on the center of square plaquettes or on the bonds of the SnTe (001) surface. For impurities far-enough from each other n compared to lattice parameter, the Hamiltonian can be simply given by where e and e ′ are the itinerant electrons living around nearest neighbor magnetic impurities S 1 and S 2 . By this, we respectively find χ center ] when the impurities are located at the center of square plaquettes and bonds between nearest neighbors. It is clear that having the responses on the same and different sublattices, these two extra cases can easily be extracted. The stronger response of impurities on different sublattices [it will be confirmed later] is dominant to calculate the entire response of these two extra configurations. So the magnetic susceptibility will show an antiferromagnetic ordering in the system. However, we would neglect these configurations in the present paper from a pragmatic point of view. Also, the bulk states do not significantly affect the RKKY coupling because it has been found that although distortion has a negligible effect on the band structure in the bulk, it can dramatically affect the Dirac surface states 7 . We note that Dirac materials that made a revolution in condensed matter physics have a strong impact on magnetic couplings.The decaying rate of RKKY coupling is faster in Dirac materials compared to 2D ordinary metals [24][25][26][27][28][30][31][32][33][34] . This faster decaying is also seen for gapped Dirac cones. In other words, gapped Dirac cones have different results than other gapped semiconductors. Impurities on the same sublattices. When two magnetic moments reside on the same sublattices of the surface (001) of SnTe (the same Sn or Te sublattices), the susceptibility can be calculated as After tedious calculations and defining γ = , C 2 = cos(2� y R y ) , C 3 = 2 cos γ cos ξ , and C 4 = 2 cos γ cos ζ . We comment that the susceptibility is function of both distance R and the angel ϕ R due to the presence of {R x , R y } . Also, we stress that the RKKY interaction between two magnetic moments on the same Sn and/or Te sublattices is the same, i.e. χ 11 Simply, it is straightforward to conclude that the coefficients C i [i ∈ 1, 2, 3, 4 ] are the origins of the sign changing in the susceptibility in the presence of the distortion.  (25), the spin susceptibility when the magnetic moments reside on the same sublattices supports a wide range of possibilities on the spin flipping and eventually on the magnetic phase transitions, which strongly depend on the gap strengths induced by xF and yF . To reveal the ferroelectric distortion effect on the RKKY interaction J RKKY we consider different possible configurations. Thus, to deduce the clear physical insights from such a complex response, we divide the following analysis into four parts, as represented in Fig. 4  Particularly, this classification is necessary due to various possible configurations for the gaps xF and yF . The specific nature of each part will be focused on the magnetic ordering and the role of short-and long-range responses. By this explicit elucidation, one can easily distinguish the novelty of electronic and/or magnetic features in TCIs compared to other gapless and gapped Dirac materials.
For the magnetic ordering, the sign of susceptibility helps, however, in all cases, it is not hard to seek for the short-and long-range behaviors of the RKKY response. It can be simply formulated considering the conditions � x/yF R/v F ≪ 1 and � x/yF R/v F ≫ 1 in Eq. (25). By this, using the simplified modified Bessel functions with γ ′ = 0.5772 being the Euler-Mascheroni constant, one finds corresponding short-and long-range responses.
Before turning to the analysis of the results, we would focus on the periodic feature of the RKKY coupling in the pristine and distorted square lattice of the SnTe (001) surface. By the lattice, one easily requires the property J RKKY αα (R, ϕ R ) = J RKKY αα (R, ϕ R + 2π) at any distance R and any gap, as confirmed by the numerical result represented in Fig. 5, however, different interaction strengths between two magnetic impurities at different positionings are obvious. This anisotropy of the susceptibility could be expected from the cosinusoidal functions (included in C i coefficients) behind the integrals in Eq. (25). A general feature in the RKKY response is the gap-dependent beating pattern with distinct periodicities, meaning that various Fermi wavelengths contribute to form the entire periodic pattern. In our system, multiple surface Dirac cones, two along the x direction and two along the y direction in the SBZ of SnTe (001) and related alloys are responsible for the beating effect in the RKKY oscillations. Now it is time to study each aforementioned part.
Pristine TCI. For the pristine case, i.e. xF = yF = 0 [see the first row of Fig. 4], following Eq. (25) we achieve the expression www.nature.com/scientificreports/ in great agreement with our previous work 43 . This clearly shows that the term F (R , ϕ R ) behind the integrals causes the positive sign of the susceptibility forever, meaning that there is no magnetic phase transition at all and the pristine system possesses a ferromagnetic ordering. Moreover, the RKKY response oscillates with R because of the cosine functions in the interference term F (R , ϕ R ) . On the other hand, since the integral is solvable analytically and no limitation is needed, one observes the decaying rate of the susceptibility serving as R −3 in both short-and long-range impurity separations. The decaying rate treatment is similar to that of undoped graphene [48][49][50][51] . All these are understandable in all black curves in Fig. 6a-c, for which the response shows a ferromagnetic spin arrangement due to the positive sign of RKKY interaction. We avoid the repetition of the results of this part since they are well-established in Ref. 43 .

Isotropically gapped TCI.
To determine the susceptibility of isotropically gap-induced TCIs due to the isotropic ferroelectric distortion, we consider | xF | = | yF | = 0 [see the second row of Fig. 4]. We would stress that for the case of proximity coupling to a ferromagnet, a mass term with the same signs at x and ′ x , is induced by the exchange field on the SnTe (001) surface to align the spin direction. Or it may be regarded as the Zeeman term arising from a perpendicular external magnetic field 40 . However, the theory of the gapped Dirac model presented in this work, refers to a sign change of mass terms for the closed Dirac cones. In this case, combining the coefficients behind the integrals in Eq. (25) leads to the expression where A ′ = ω 2 + � 2 0 and Regarding the magnetic ordering of this configuration, we comment that both ferromagnetic and antiferromagnetic phases are possible to emerge because of the mixture of sine and cosine functions in the interference terms behind the integrals. However, to elucidate the role of short-and long-range impurity separations, we stick to Eq. (26) as well as we use the mathematical identity 52 resulting in the final expressions respectively for the short-range � 0 R/v F ≪ 1 and for the long-range � 0 R/v F ≫ 1 with the help of Laplace method www.nature.com/scientificreports/ Thus the short-range responses demonstrate the ferromagnetic ordering (stemming from the interference term F (R , ϕ R ) in which the direction ϕ R is not important in changing the sign of the susceptibility) with the decaying rate of R −3 like the pristine TCIs, while the long-range response decays as R −3/2 � 3/2 0 e −2 � 0 R/v F and both ferromagnetic and antiferromagnetic orderings are expected to appear due to the interference term G (R , ϕ R ) in which the direction ϕ R is a matter in sign switching of the susceptibility. It is worthwhile mentioning that the RKKY interaction can not be influenced at all with the gap at short-range distances. The latter, in turn, means that the decaying function is a function associated with an extremum at which the critical distance R or gap 0 can be characterized. We will come to this last point later. The phase of the RKKY oscillations in metals becomes random in the presence of static non-magnetic impurities. In this case, averaging over various impurity configurations leads to an exponential factor appearing at distances larger than the mean free path 53,54 . An exponential decay on averaged RKKY couplings was also reported in disordered graphene 55 . In gapped graphene, RKKY coupling also experiences an exponential decay at large distances giving rise to Heisenberg interaction 48,56,57 .
In this configuration which the distortions are applied with the same strengths along different axes, referring as orange curves in our numerical results shown in Fig. 6a-c, the short-range treatment is independent of the gap, as expected from Eq. (32), while long-range one is associated with interesting spatial-dependent magnetic www.nature.com/scientificreports/ phase transitions expected from Eq. (33). Note that, physically point of view, the atomic packing factor for the same sublattices on the SnTe (001) surface along the x or y direction, Fig. 6a, is smaller than the atomic packing factor along the different directions, Fig. 6b and c. Thus, the shorter accessible distances are expected for ϕ R = 0 . Although the short-range response is gap-independent, to describe long-range gap-dependent behaviors, we numerically plot Fig. 7 based on Eq. (25) at large-enough R = 200 Å as a function of the gap potential 0 . Interestingly, the feature concluded in Eq. (33) comes up, implying that the RKKY response is an exponentially decaying function along all directions ϕ R with a minimal gap at which the response starts to switch its trend [see the inset panel of Fig. 7].

Coexistence of gapless and gapped TCI.
For the individual case of | xF | = 0 and | yF | = 0 [see the third and fourth rows of Fig. 4] we have The same expression can be obtained for the case of | yF | = 0 and | xF | = 0 subsituating C 2 → C 1 and C 1 → C 2 in the third first terms. Following Eqs. (26) and (31), we obtain respectively the effective short-range � 0 R/v F ≪ 1 and long-range � 0 R/v F ≫ 1 responses wherein the decaying rates are similar to the previous case, whereas the interference term for the long-range case is quite different. While depending on the C i coefficients, one expects almost the same behaviors for the RKKY interaction of gapped TCI along the x or y direction at short-range distances, the ϕ R -dependent long-range interaction would be different [see blue and red curves in Fig. 6a-c].
This part provides the first interesting novel remark of the current paper for which the SnTe (001) surface confronts two phases at the same time, a gapless Dirac cones along the x/y direction and a gapped one along the y/x direction. As mentioned in the introduction, multiple surface Dirac conses in different topological materials may show such an interesting feature as well [12][13][14][15] . However, the effects of such coexisted gapless and gapped phases on RKKY response can be rarely found in other Dirac materials. Thereby, one expects nontrivial behavior of the RKKY response.
In Eq. (35), we found the underlying analytical expressions behind the short-and long-range responses. For the sake of completeness, it is necessary to discuss the intermediate-range RKKY response as well. Following the above points and considering structural symmetry of the SnTe (001) surface, one can expect some spatial (34)  (R, ϕ R )/J 2 C t with ± xF and ± yF . In addition to the above-mentioned points, we intend to discuss the origin of the decreasing and increasing trend of the susceptibility with the gap, which is that of nontrivial point reported before. In Fig. 8, the susceptibility decreases with |� x/yF | and after a critical gap potential, which is obviously R-and ϕ R -dependent, increases. The reason can be traced back to the metallic phase of the system in one direction, while the gapped phase along another one. This implies that the reason for the observed nontrivial trend can be understood from the fact that the states around the band edges for the gapped Dirac cones also matter. These states which belong to the valence band, contribute to the RKKY interaction. The susceptibility is intensified if one can enhance the density of the band edge states in the valence band.
Indeed, the susceptibility is proportional to the correlation function defining as D 0 (ω, R) = − 1 π ImG 0 αα (ω, R, 0) . To this end, we stick D 0 (ω, R) at the energies equal to the gap size, e.g. xF [the analysis is precisely the same for yF ] in Fig. 9 setting ϕ R = 0 [as the representative test case] and yF = 0 . We intend to find the relation between the correlation function D (� xF ) = D 0 (ω = � xF ) and { xF , R} to understand the physical reason behind the nontrivial behavior of the susceptibility in Fig. 8. Two qualitatively different regions as a function of the distance R may be seen: For short-range R < R c [R c = 8 Å] a high D (� xF ) is observed, which on increasing the distance decreases gradually for all R > R c [see below for the origin of the critical R c ]. A closer look at this evolution of D (� xF ) is provided in the inset panel of Fig. 9, in which we label the critical distance R c . Clearly, D (� xF ) increases slightly with the gap at certain R < R c , while it increases up to a critical gap at certain R > R c and starts to decrease. The critical gaps at R > R c are labeled by the red solid dots. From the contour plot, it is also noteworthy that the size of this critical gap is inversely proportional to the distance R, meaning that it decreases with R. The rapid change in D (� xF ) in terms of the gap edges for R > R c is the main reason of the suddenly changes in the RKKY response at R-and ϕ R -dependent critical gap potentials. The same argument is valid for yF .
It is worthwhile to understand the reason behind the critical R c from the susceptibility in Eq. (25). For simplicity, we restrict ourselves to the band edge ω = � xF and set yF = 0 . The integrands of the terms show a maximum value of R c ≃ 8 Å, for the considered range of distortion potential 0 < � xF < 0.1 , meaning that the slope of this function approaches zero at this critical R c . The same arguments can be set for the case of ω = � yF and xF = 0 . It should be restressed that the critical R c may be different for other directions ϕ R .
To summarize non-trivial behavior of the RKKY interaction in the short and intermediate range of impurity separations, high correlation function at the band edges leads to an enhanced host states in the valence band giving rise to stronger magnetic response. Thereby if magnetic impurities are resided along the gapped Dirac cones, www.nature.com/scientificreports/ the RKKY interaction would be based on the mentioned mechanism which is increasing with the gap size at short-range distances. It should be reminded that at long-range separations, the RKKY interaction exponentially decays with the gap size. Simultaneously, massless Dirac fermions along the perpendicular direction indirectly affects RKKY response especially when Fermi energy lies in the gap. On the other hand, if impurities are aligned along the gapless Dirac cones, they would couple to each other by means of the massless host electrons while massive Dirac fermions also affects RKKY interaction indirectly via the valence band states [please see Fig. 10d]. For a physical interpretation, we would state that the system possesses a new quasifermion in this situations in the presence of simultaneous massless and massive fermions. In fact, one is allowed to define ψ quasifermion (k) = i P i ψ massless fermions i (k) + P ′ i ψ massive fermions i (k) in terms of the orthogonal massless and massive wave functions with corresponding probability distributions {P i , P ′ i } for the i-th orbital, which are responsible for the nontrivial RKKY treatment compared to the individual massless and massive fermions. We just mention this point here to justify the new created fermions when coexisting gapless and gapped phases. More details of the momentum-dependent Hamiltonian and the band structure of these introduced quasifermions can be tracked from the Eq. (3) and corresponding Fourier transforms, which are out of the scope of the present paper.
Anisotropically gapped TCIAnisotropically gapped TCI. In the case of anisotropically gapped TCI, the sign of xF yF is important and one would report the general formulation of the susceptibility in Eq. (25) for this part [see the last row of Fig. 4]. However, the short-and long-range responses can be read as Figure 9. The behavior of the correlation function at E = xF as a function of ferroelectric distortion potential xF and distance R for ϕ R = 0 at yF = 0 . The critical gap for which D (� xF ) turns to the increasing trend at a certain R > R c decreases with the distance R labeled by red solid dots.
It is not easy to deduce a general treatment for this case due to the complexity of the interference term at longrange impurity separations. But, the decaying rates are similar to two previous cases.
To understand the conditions for which the magnetic phase transition occurs by symmetry breaking and band gap opening in such a situation, in Fig. 10 we show the RKKY interaction strength over the full allowed range of gaps for different values of distance R > R c between two magnetic impurities with ϕ R = 0 . The breakdown of the antiferromagnetic and ferromagnetic phases is governed by the positive and negative signs in the color bars. For a quantitative analysis of J RKKY αα (R, ϕ R )/J 2 C t , we label zero response by dotted lines in the color bars. We noted that there is no asymmetric behavior for RKKY response concerning the band inversion effect originating from the gap signs. It should be highlighted that at large separations of Fig. 10d, stronger RKKY coupling is obvious when magnetic moments located along ϕ R = 0 are coupled to each other via the host massless Dirac fermions along the x-axis, xF = 0.
In particular, independent of the distance R, the amplitude of the RKKY interaction attains a maximum in the absence of any mass term induced by the ferroelectric distortion, while decreases, as discussed before, with the gap potentials. Following the long-range RKKY response in Eq. (36), Fig. 10 shows that for the cases of xF yF = 0 , the RKKY response decreases exponentially with the gap such that at large-enough distances, the magnetic phase transitions take place. For R = 10 Å [see Fig. 10a] it is no surprise that the J RKKY αα (R, ϕ R )/J 2 C t gets positive values because this is not the effective distance for the sign switching of the RKKY response. For R = 30 Å i.e. Fig. 10b, again, no transition takes place, while for distance R = 50 Å in Fig. 10c, transitions are appearing for |� x/yF | 0.06 . However, if we look at the magnetic moments on the same sublattices with longer separation R = 70 Å in Fig. 10d, the transition occurs in a wider range of gaps involving the negative RKKY interaction, namely |� x/yF | 0.03 . We comment that for ϕ R = π/2 , one would find the same feature except that in turn, stronger RKKY coupling would be observed if massless Dirac fermions along the y-axis play the role of host carriers, yF = 0.
For other ϕ R angels, we discuss RKKY response to the magnetic moments located along the direction ϕ R = π/6 for R = 25 Å and 50 Å in Fig. 11a and b, respectively, as well as with ϕ R = π/4 for R = 30 Å and 50 Å in Fig. 11c and d. For the same range of gaps we examined, the symmetry property of the RKKY interaction with respect to the gap sign is broken, i.e. J RKKY , −� x/yF ) . However, regarding the Eq. 25 and numerical results in Fig. 11, the RKKY coupling depends on the sign of yF yF [see the first and third as well as the second and fourth quarters in all panels of Fig. 11]. The above interpretation is valid for both considered angles.
Regarding the magnetic phase transition at ϕ R = π/6 for R = 25 Å for the gaps with the property � xF � yF < 0 , the spin flipping occurs, as shown in Fig. 11a. If we seek for such a phase transition at longer distances, however, we find it in the majority of gaps with the property � xF � yF > 0 , as shown in Fig. 11b. It is interesting that based on the amplitude variation of the responses, the transition in the short range of ϕ R = π/6 is stronger than that of long-range one. In the case of angle ϕ R = π/4 , the situation is different, and much fewer amplitudes corresponding to the magnetic phase transitions are obtained. While for the short-range RKKY interaction, Fig. 11c, there is no phase transition, a weak transition emerges for the long-range case when the gaps only satisfy the feature � xF � yF > 0.
We notice that the critical susceptibility at which the magnetic phase transition appears depends strongly on the distance R and angle ϕ R . Thereby, the behavior of spin susceptibility for two magnetic impurities placed www.nature.com/scientificreports/ on different atomic sites (which possesses quite different interaction conceptually) may be quite different and notable. In the next section, we will elucidate this.
Impurities on different sublattices. When two magnetic moments reside on different sublattices of the surface (001) of SnTe, the susceptibility can be calculated as where D 1 = 2 cos 2 (� x R x ), D 2 = 2 cos 2 (� y R y ) and D 3 = 4 cos γ cos(� x R x ) cos(� y R y ).
The emergent conclusions extracting from the above susceptibility, i.e. Eq. (38) could be listed as the following. Firstly, we conclude that the RKKY response is negative anyway, providing an antiferromagnetic spin configuration for impurities sitting on the different sublattices. Secondly, the RKKY interaction is an even function of the gap terms, meaning that the gap term sign is not a matter subject from the point of RKKY response. In other words, it is not possible to detect the magnetic phase transitions when the impurities are placed on different sublattices through the magnetic response, in contrary to the previous the same sublattices case.
Similar to the same sublattices case, the RKKY response is periodic in ϕ R , as expected from the structure symmetry of the SnTe lattice. This fact is numerically verified in Fig. 12 in the absence and presence of the ferroelectric distortion to highlight the distortion effect on the RKKY intensity. From this plot, one immediately would conclude that the presence of both xF and yF simultaneously with different possibilities for the sign can not influence the RKKY response [see the orange and green curves on top of each other]. This is a direct consequence of power 2 of the gaps in Eq. (38). Furthermore, it exhibits a beating type pattern with gap-dependent amplitudes originating from the Fermi wave vectors involving in the response.
The RKKY interaction is again strongly anisotropic concerning the direction of applied displacements. However, for a general discussion of the distortion configuration effects on the magnetic response of impurities resided on different sublattices, we again divide the following analyses into four parts considering (1) pristine TCI, (2) isotropically gapped TCI, (3) coexistence of gapless and gapped TCI, and (4) anisotropically gapped TCI.
Pristine TCI. Considering xF = yF = 0 [see the first row of Fig. 4] in Eq. (38), one achieves www.nature.com/scientificreports/ Interestingly, similar interference term as the case of the same sublattices can be obtained, in an excellent agreement with Ref. 43 , while the magnetic ordering is different due to the negative sign of susceptibility referring to the antiferromagnetic phase. It is worth mentioning that the pristine RKKY amplitude is different in this case compared to the same sublattices due to the difference between the modified Bessel functions of the zero and first kinds. The corresponding results are shown as black curves in Fig. 13a-c for which the short-and long-range responses decay as R −3 with the impurity separation, whereas RKKY interaction oscillates with intermediate R due to oscillating interference term L (R , ϕ R ) . However, no magnetic phase transition is observed at all.
Isotropically gapped TCI. Turning now to the isotropically gapped TCIs with the feature | xF | = | yF | = 0 [see the second row of Fig. 4] leads to the expression With the help of Eqs. (26) and (31), one finds for the short-range � 0 R/v F ≪ 1 and for the long-range � 0 R/v F ≫ 1 using the Laplace method The reason can be understood from the same interference terms behind the integrals in the equation above, while different interference terms came up in Eqs. (32) and (33). These behaviors are greatly confirmed numerically in orange curves (under the green ones) of the Fig. 13a-c describing the decaying rate above-mentioned for the short-and long-range responses as well as oscillation for the intermediate distances.
The exponential decaying of the RKKY interaction with the gap potential −J RKKY αβ (R, ϕ R )/J 2 C t is investigated at long-range regime in Fig. 14 along different directions ϕ R . The process of exponentially decay of the RKKY interaction with the gap potential is the recovery of Eq. (43).

Coexistence of gapless and gapped TCI.
For the individual case of | xF | = 0 and | yF | = 0 [see the third and fourth rows of Fig. 4], we have Substituting D 1 → D 2 and D 2 → D 1 in the third first terms give rise to another case of | yF | = 0 and | xF | = 0 . Following Eqs. (26) and (31), we obtain respectively the effective short-range � 0 R/v F ≪ 1 and long- showing the same decaying rates as the previous case, while different interference terms emerge. Also, a simple comparison between blue and red curves in Fig. 13 and Eq. (45) alarms that the short-range RKKY responses are somewhat gap-independent, whereas the long-range response comes up with a difference in oscillations for different configurations. It is necessary to mention that the decaying rates here also behave similarly to graphene with oscillations depending on the various gap configurations [49][50][51] .
As expected, the difference between the case of {| xF | = 0 , | yF | = 0} and {| yF | = 0 , | xF | = 0} belong to the long-range responses except the case of ϕ R = π/4 [ Fig. 13c] for which both behave similarly independent of the impurity separation R because in this case, we have D 1 = D 2 and there is no priority for these configurations at all.
The effect of gap formation on the RKKY interaction at an intermediate impurity separation of R = 25 Å is studied in Fig. 15 for the complementary angles ϕ R when the distortion is applied along the X 1 − Ŵ − X 1 line or along its perpendicular orientation X 2 − Ŵ − X 2 . The response is monotonically antiferromagnetic for all gap www.nature.com/scientificreports/ terms, as expected. Moreover, a correspondence between the distortion direction and the substitutional orientation of impurities is observed in Fig. 15, i.e. the RKKY response for the distortion applied along the y-axis with sitting impurities along ϕ R = 0 corresponds to the situation in which the distortion is applied along the x-axis with ϕ R = π/2 . Again we observe an increased treatment of the RKKY interaction after a critical gap, which is originated from the high density of the band edges in short-range distances or the role of metallic states along the other direction. The physic behind this enhancement is similar to what we presented for the same sublattices. We comment that this is also a nontrivial phenomenon as the same sublattices case when both gapless and gapped phases coexist in the whole band spectrum. An increase behavior of the RKKY coupling with the gap size at intermediate-range of impurity separations is traced back to those accumulated valence band states at the band edges affecting RKKY response. Interestingly, such nontrivial increasing trend occurs if magnetic impurities are located along arbitrary orientations. Let us turn to the last distortion configuration.
Anisotropically gapped TCI. Referring again to Eq. (38), the situation of interest for the presence of both xF and yF is when they compete to change the RKKY amplitude [see the last row of Fig. 4], in which the effective short-range � 0 R/v F ≪ 1 and long-range � 0 R/v F ≫ 1 responses are given by where We again intend to systematically study the impact of the competition between xF and yF on the RKKY response. To do so, we plot Fig. 16 at an intermediate distance R = 30 Å for different directions, namely (a) ϕ R = 0 , (b) ϕ R = π/6 and (c) ϕ R = π/4 . Independent of the direction ϕ R , the RKKY response is symmetric concerning � xF � yF < 0 or � xF � yF > 0 , in which the RKKY amplitude is weaker than the other areas. Moreover, we once more confirm that there is no magnetic phase transition at all when the magnetic impurities are resided on different sublattices. A quick comparison between different directions shows that the case of ϕ R = π/4 leads to the maximum RKKY amplitude due to the equal contribution of the Fermi wave vectors in D i coefficients.
Regarding different magnetic ordering of two configurations (when two impurities residing on the same or different sublattices), let us briefly mention that atomic orbitals and topology of the structure play the main role in magnetic ordering. Hamiltonian in Eq. (5) is written in the basis set of Te and Sn p-orbitals. When the same sublattices are considered, the packing factor is lower compared to the different sublattice case. It means that atomic separation in the same sublattice case is larger than the different sublattice case. Therefore, it is reasonable if the overlapping of the p-orbitals between neighbors atoms influences RKKY coupling. The effect of the atomic structure on RKKY interaction is seen also in graphene such that magnetic ordering along direction AA-zigzag is ferromagnetic while it is antiferromagnetic for AB-zigzag direction [48][49][50][51] . www.nature.com/scientificreports/ Discussion. The current distortion-induced RKKY response for the magnetic impurities on the same and different sublattices are for the case of undoped Dirac cones. However, we should mention that the presence of electron and/or hole-doping changes the results and the competition between the distortion and doping concentration is an important study because doping may induce new magnetic phases to the system due to the breakdown of the lattice symmetry 43,48,51 . Additionally, surface roughness as a perturbation preserving time-reversal symmetry is always present in real materials. This perturbation breaks the mirror symmetry locally on the surface, however, the mirror symmetry may still be preserved on the average macroscopically point of view, i.e. when the variation of the atomic places is slow enough 7 . Thereby, we argue that, in this case, the gapless results work out well at long-range impurity separations, while one needs to apply a continuous gapped Dirac model for short-range distances on average. Thus, in the smooth variation of roughness for short distances, the gap is opened in the band spectrum, and so our results for RKKY interaction at short ranges can be considered, while for long-range distances the gapless results are applicable. In the presence of sharp variation of the atomic positions, one should similarly look at the mirror symmetry whether it is preserved or not.
Furthermore, we would state that the current results are based on the zero temperature, however; the temperature dependency of short-and long-range magnetic couplings in the absence and presence of symmetry breaking and doping can be studied as well within the finite-temperature self-consistent field approximation 58,59 . Temperature highlights the screening response of the conducting electrons, which strongly affects the indirect exchange interaction between magnetic impurities, resulting in temperature-dependent magnetic phase transitions.
Finally, we emphasize that the magnetic scattering potential may also happen in the presence of magnetic impurities. In such a situation, the RKKY interaction should be approached beyond the linear response theory and the local magnetization should also be considered 60,61 . This revises the theory of linear off-resonant RKKY interaction between magnetic impurities and considers the impact of impurity resonances. We leave the results of all these matters to come up in our future works.

Summary
In summary, we have addressed the ferroelectric distortion effects on the RKKY interaction between two magnetic impurities on the (001) surface of SnTe and related alloys as available TCIs. We have shown that the lowenergy Hamiltonian of TCIs receives a mirror symmetry breaking between multiple Dirac cones in the presence of the distortion, ensuring the gap at Dirac points. The coexistence of the gapless and gapped phases as well as the isotropically and anisotropically gapped phases grant a fascinating interference correlation between magnetic moments, which manifest themselves in different magnetic ordering. While other Dirac materials can be gapped or gapless, the key point of the present paper, especially is the the coexistence of the gapless and gapped phases and their nontrivial and novel fingerprints in the gap-induced RKKY response. The results are significantly different compared to the pristine TCI and other Dirac materials.
Our numerical results present different magnetic phase diagrams depending on the position of impurity magnetic moments on the (001) surface of TCIs. We have proposed a critical impurity separation of R c for which the RKKY coupling shows different behaviors with the distortion strength for separations beneath and above R c . Although the distortion does not lead to a magnetic phase transition for the magnetic moments resided on different sublattices, an irregular ferromagnetic ↔ antiferromagnetic ordering emerges as soon as the impurities reside on the same sublattices with the separation above R c . We further demonstrated that it is feasible to manipulate this magnetic ordering by tuning the distortion strength. Our work paves the way to design protocols in tuning the multiple surface states in TCIs using the external magnetic impurities and ferroelectric distortion for the spintronic and valleytronic applications.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.