1072 Application or Design Study  (A) Projection View compares the projections
in different ranking schemes. (B) Ranking Comparison View supports analysts to
seek explanations of different ranking schemes at an instance level. (C)
Ranking Tabular View displays the data attributes, lists the current ranking
schemes, and enables analysts to adjust the ranking.

# Inspecting the Process of Bank Credit Rating via Visual Analytics

Qiangqiang Liu†<EMAIL_ADDRESS>Quan Li
<EMAIL_ADDRESS>corresponding author Zhihua Zhu∗ Tangzhi Ye∗
Xiaojuan Ma<EMAIL_ADDRESS>†ShanghaiTech University ∗Corporate Development
Group, Tencent ‡The Hong Kong University of Science and Technology

###### Abstract

Bank credit rating classifies banks into different levels based on publicly
disclosed and internal information, serving as an important input in financial
risk management. However, domain experts have a vague idea of exploring and
comparing different bank credit rating schemes. A loose connection between
subjective and quantitative analysis and difficulties in determining
appropriate indicator weights obscure understanding of bank credit ratings.
Furthermore, existing models fail to consider bank types by just applying a
unified indicator weight set to all banks. We propose RatingVis to assist
experts in exploring and comparing different bank credit rating schemes. It
supports interactively inferring indicator weights for banks by involving
domain knowledge and considers bank types in the analysis loop. We conduct a
case study with real-world bank data to verify the efficacy of RatingVis.
Expert feedback suggests that our approach helps them better understand
different rating schemes.

Human-centered computingVisualization

Introduction

Bank credit rating classifies commercial banks into different levels on the
basis of the analysis of both the publicly disclosed information and part of
the internal bank information according to a certain quantitative standard
[20, 24, 28, 25]. These levels will be succinctly represented by a combination
of numbers and letters, which is convenient for supervisory agencies to
conduct an intuitive judgment on the risks of commercial banks and adopt
corresponding supervisory measures [20, 11]. Bank credit rating methodologies
have gone through several stages such as subjective judgments, rule-based
scorecard analysis, and computational modeling [9, 2]. Most of these
approaches leverage rule-based measures combined with automated methods [11,
26] for credit rating assistance. Although these approaches have demonstrated
promising performances, domain experts still have the following concerns when
directly applying the existing credit rating models to their business
scenarios. (1) Obscure understanding. One purpose of the domain experts is to
generate an internal bank rating that can be produced consistently by
different analysts using the same information with high reliability.
Therefore, the methodologies should be intuitive and explainable to benefit
from comments and suggestions from the experts advancing the bank credit risk
assessment methodologies. (2) Loose connection between subjective and
quantitative analysis. Previous bank rating models are extremely biased
towards a heavy quantitative approach that requires many quantitative
indicators. Therefore, they are not suitable for those imperfect credit
systems and banks with limited financial indicators [20, 11]. Meanwhile,
domain experts have their knowledge in the bank rating decision-making
process. Integrating domain knowledge and experience with quantitative
indicators can improve bank credit rating reliability. (3) Unified indicator
weights. Existing credit rating methods apply a unified indicator weight set
to all bank entities; however, according to the experts, different types of
banks should be assigned to different weights even for the same indicator [5].
To tackle these issues, we identify domain experts’ primary concerns regarding
bank credit rating and propose a visual analytics system, namely, RatingVis to
assist them in exploring and comparing different bank credit rating schemes.
Our primary contributions can be summarized as follows: 1) we elicit design
requirements of bank credit rating from the literature and an observational
study; 2) we propose interactive schemes to infer indicator weights; 3) we
design a visual analytics system to help explore and compare different rating
schemes.

##  1 Related Work

Bank credit rating plays an important role in solving information asymmetry
between depositors and banks [22]. For example, CAMEL rating system [31]
evaluates the capital adequacy, asset quality, management, earnings, and
liquidity of financial institutions, and adopts a five-level scoring system to
rate the operation and management level of commercial banks. Moody’s bank
rating system [3] focuses on asset quality evaluation by understanding the
bank’s internal loan approval policies and procedures, internal risk rating
methods, the means of managing asset quality and risk to obtain bank-related
information and rating data. Standard & Poor’s rating [21] of banks takes into
account eight factors such as economic and industry risks, capital, and
returns. However, there is no universal global standard for rating. Some
methods reply on different data indicators, making it impossible to generalize
to other scenarios. Moreover, experts’ knowledge of the credit risk is missing
in these rating models. Our work integrates domain knowledge into rating and
combines quantitative analysis with qualitative evaluation.

Entities for rating are inherently multi-attribute data items and several
visualization systems that provide interactive multi-attribute data item
ranking were proposed [17, 6, 30, 32, 14, 23]. ValueCharts [4] and LineUp [8]
allow analysts to create customized rankings by clicking and dragging
attributes with adjustable attribute weights. However, they assume that users
are able to quantify the importance of particular attributes. To resolve this
issue, researchers studied what weight sets lead to a specific ranking. Podium
[30] allows users to drag the table rows to rank data based on their
perception of data value. WeightLifter [19] facilitates the exploration of
weight spaces, thereby better understand the sensitivity of a decision to
weight changes. The unified attribute weights derived by the above methods are
not applicable to our case since banks with various types should be assigned
different weights even for the same indicator.

##  2 Background and Observational Study

Ensuring the quality of bank rating has attracted attention from many parties.
The world-famous rating agencies such as Moody’s and Standard & Poor’s
evaluate the financial and operational strength and resilience of banks based
on operation-related indicators and generate the letter-grade credit. To
understand the existing practice of bank rating, we worked with one bank
rating expert (E.1, male, age: $32$), one risk management expert (E.2, male,
age: $35$), one financial data analyst (E.3, female, age: $24$), and one bank
credit expert (E.4, male, age: $28$), who attempted to address the rating
problem by ranking. They selected a set of attributes and employed certain
multi-criteria decision-making schemes such as Analytic Hierarchy Process
(AHP) [1] and Technique for Order Preference by Similarity to Ideal Solution
(TOPSIS) [13] in Statistical Product and Service Solutions (SPSS) for ranking.
Based on the distribution of certain indicators, they divided the ranking into
several segments. The banks in a higher up segment are considered better than
those in a lower segment. Nevertheless, several issues arise. First, the
ranking relies on distance functions to obtain the pairwise similarity.
However, different indicators have different numerical scales, thereby
requiring normalization. Nevertheless, even the experts who are very familiar
with the data have difficulty in estimating the attribute weights after
normalization. “When we classify banks, we can tell the approximate level of
banks by observing certain attributes.” However, the subjective perception of
attribute importance is “difficult to quantify”, making it challenging to take
this intuition to verify the appropriateness of distance functions. Second,
banks are of different types. For example, banks in China are classified into
“large state-owned commercial banks”, “joint-stock commercial banks”, “city
commercial banks”, and “rural commercial banks”. Generally, the credit rating
of the four categories of banks shows a decreasing trend as a whole, but some
low-level banks may perform better than high-level banks. In existing bank
ratings, experts apply the same weight to the same indicator of different
types of banks. However, they realized that “the importance of the same
indicator among different types of banks should be different.” For instance,
banks with a high liquidity ratio of assets indicate a strong financing
ability, e.g., large state-owned commercial banks have many means to obtain
funds while the financing channels of city commercial banks are narrow.
Therefore, experts should pay attention to the liquidity ratio of city
commercial banks. To sum up, we should meet three requirements: R.1 Infer
indicator weights by interaction. Conventionally, the experts rank banks by
assigning weights to quantify attribute contributions, by which they cannot
efficiently determine which and to what extent certain attributes are
important since they only have a holistic understanding of the data.
Therefore, they required that the indicator importance could be interactively
inferred; R.2 Refine indicator weights for different types of banks. Applying
uniform weights to all kinds of banks may result in an irrational credit
rating result. As indicated in [5], ratings can communicate externally as well
the overall riskiness of banks’ assets, and “different risk weights are
assigned to assets with different ratings according to Basel capital rules”.
The experts required an interactive mechanism to refine indicator weights that
consider bank types; R.3 Facilitate comparison among rating schemes. The
experts wished to preserve previous schemes for comparison so they can
understand whether the adjustment leads to a better result. Thus, our approach
should facilitate comparison among different schemes.

##  3 Back-End Engine

In this section, we derive indicator weights and constraints to generate
ranking scores for all items, followed by three ranking schemes.

Modeling Ranking SVM. We adopt Ranking SVM [32, 12] to derive a set of
indicator weights on the basis of the rankings in the data table [30], which
optimizes a SVM hyperplane to the ranking problem with pairwise constraints.
Instead of a full set of data with labels, a limited set of pairs of data
$d_{i}$ and $d_{j}$ and a label are available to derive whether $d_{i}$ is
better or not. The input of the model is difference vectors for data pairs,
e.g., $d_{i}-d_{j}$, and the output of the model is which point is better
[12]. Specifically, Ranking SVM transfers the pair of $(d_{i},d_{j})$ and
their relative ranks to a tuple: $d_{i}-d_{j}=1$ if $d_{i}$ is preferred;
otherwise $d_{i}-d_{j}=-1$. The generated model predicts which one is better
given a pair of points. To avoid an empty result, all constraints are modeled
as soft constraints instead of hard ones so that user interaction can always
generate a set of attribute weights that model constraints as much as possible
[12].

Deriving Constraints. To obtain SVM’s linear separator, the ranking problem is
transferred into a binary classification problem. Particularly, the labeled
data is generated for Ranking SVM using the data items with which the user has
interacted, e.g., dragging to a new position and there are $k$ marked rows.
Without loss of generality, for $k$ points $\\{d_{l_{1}},...,d_{l_{k}}\\}$
with indices $[l_{1},...,l_{k}]$, the set of all combinations of pairwise
difference vectors are created as training instances [12]. In other words, for
$i,j\in\\{1...k\\}$, where $i\neq j$, a training tuple is derived, i.e., each
training instance is the difference between a pair of rows $d_{i}$ and
$d_{j}$, classified as $y=1$ if $d_{i}$ is ranked higher than $d_{j}$, or
$y=-1$ if $d_{i}$ is ranked lower than $d_{j}$. In this study, we set $k=6$ to
ensure a minimum amount of training data for deriving the attribute weight
vector after experimental analysis.

Calculating Ranking Score. After transforming user interaction and learning
the Ranking SVM model, a weight vector w is obtained for us to rank the data
items. We calculate individual dot products of $w$ with each data item to
generate a rank score as: $r(d_{i})=\textbf{w}\cdot
d_{i}=\sum_{j=1}^{m}w_{j}d_{ij}$ with the highest one corresponding to the top
rank.

We introduce three ranking schemes based on Ranking SVM: Scheme 1: Local
Weight. We take original Ranking SVM and only consider the local information
as the ranking scheme; Scheme 2: Global Weight. We generate a global weight
sequence for all banks. When analysts drag the bank $b_{j}$ to the position of
$k$, we sample a part of banks according to the proportion of each bank type
in the original ranking from $1$ to $k-1$ and pair it with $b_{j}$ to form a
positive pair, i.e., the label of $d_{q}-d_{k}$ ($q\in[1,k-1]$) is $1$.
Similarly, we sample banks to form a negative pair with $b_{j}$ from $k+1$ to
$n$ according to the proportion of each bank type, i.e.,the label of
$d_{p}-d_{k}$ ($p\in[k+1,n]$) is $-1$. In this way, we generate a general
weight set suitable for different types of bank; Scheme 3: Type Weight. When
generating pairwise training data, we take bank types into account and
generate a set of weights for each bank type. Supposedly we have
$b=\\{r_{1},...,r_{i},...,r_{m}\\}$ that indicates bank $b$ has $m$ indicators
and $s(b_{j})$ indicates the ranking score of bank $b_{i}$ and we define
$rankSeq=<b_{1},...,b_{j},...,b_{n}>$ as an increasing ranking order of banks,
i.e., the number of banks is $n$, among which
$s(b_{1})>...>s(b_{j})>...s(b_{n})$. We separate the ranking of different
types of banks in $rankSeq$ to form a ranking for each bank type. When
analysts drag a bank $b_{j}$ to the position of $k$, we check the position of
$k$ of $b_{j}$ in the ranking of each bank type and take $6$ pairs near the
position of $k$, e.g., $b_{k-1}$ is better than $b_{k}$ and $b_{k+1}$ is worse
than $b_{k}$. For each type, we generate a set of weights to calculate the
ranking score of the banks of that type.

We adapt an entropy discretization method to transfer ranking to rating [7,
16]. We first sort ranking scores and treat each score as a segmentation
point, and then calculate the entropy of the left and right parts of each
point. We consider the division with the minimum entropy as the first
division. We repeat the above procedure until we have $5$ divisions (we
determine the value of $5$ after discussion with experts). We round down the
score of each bank to multiples of $5$. We use $X$ to indicate the random
variable of the score and sort the bank scores as $(x_{1},x_{2},...,x_{n})$.
$P(x_{i})$ means the probability of score $x_{i}$. The entropy of $X$ is:
$H(X)=E[-logP(x_{i})]=-\sum_{i=1}^{N}P(x_{i})logP(x_{i})$. Assume there are
$k$ different scores among all the ranking scores of banks and $k<n$. We sort
the $k$ scores as $(u_{1},u_{2},...,u_{k})$, which can be considered as
breakpoints of the continuous values of $x_{1},x_{2},...,x_{n}$. Then, we
choose one point that minimizes the entropy from the candidates. We repeat the
process until we have $4$ breakpoints, forming $5$ ratings and complete the
discretization.

##  4 Front-End Visualization

We propose RatingVis to infer indicator weights by interactions and facilitate
exploring and comparing bank credit rating schemes.

Inspired by Lineup [8], we design Ranking Tabular View to help derive
indicator weights from user interaction, evaluate the indicator contribution
to the ranking, and support the comparison of different ranking schemes (R.1 –
R.2). As shown in Inspecting the Process of Bank Credit Rating via Visual
Analytics(C), we leverage a table to present the raw data. It displays the
name of data items, rank, institution types, and the associated attributes.
Analysts can perform a drag-and-drop operation to manually rank data based on
their perception of the relative value of the data and domain knowledge. Green
indicates the ranking is adjusted higher and red indicates the opposite
adjustment; the deeper the color, the more the adjustment. The system then
infers a new indicator weight set using the above-mentioned Ranking SVM
approach. The attained weight vector will be applied to multiply the
normalized indicator values for calculating the indicator contribution. The
sum of the contributions of all the indicators for one data item, i.e., the
rank score, is used to determine the order of the data item. If analysts are
satisfied with the current ranking scheme, they can click on the “Save Weight
Scheme” button and the current ranking scheme will be added on the right side
of Inspecting the Process of Bank Credit Rating via Visual Analytics(C). The
associated weight vector will be applied to other views as well. As shown on
the right side of Inspecting the Process of Bank Credit Rating via Visual
Analytics(C), we present each ranking scheme as a separate column. These
columns adapt bars with different colors to represent different indicators and
the length of the bars indicates the contribution of that indicator. To
compare different ranking schemes, we arrange them horizontally and connect
identical data items across different rankings with blue lines. When selecting
a particular data item, a blue bold line connects all the identical data items
sequentially. The contributions are divided into positive and negative
contributions. Positive contribution (the right side) indicates that the
higher the indicator value, the higher the ranking score, and the negative
contribution indicates that the higher the indicator value (the left side),
the lower the ranking score.

Projection View creates low-dimensional projections and preserves local
similarities to convey neighborhood structure [15] (R.3). When designing the
projection view, we leverage the same weight set and normalized attribute
values as used in the ranking tabular view to obtain 2D projection by e.g.,
t-SNE [10, 27, 29, 18] (Inspecting the Process of Bank Credit Rating via
Visual Analytics(A)). Four subviews represent the projections corresponding to
the attribute weight vectors from default, local, global, and type weight
schemes. The colors indicate different ratings and we encode the bank asset
size as the circle size. Analysts can lasso circles on any projection space
and all identical ones will be connected via curves.

Ranking Comparison View. In previous schemes, different item and the dragged
item form pairs of either positive or negative samples for training the
Ranking SVM model, allowing analysts to understand what the impacts of
different pairs on the ranking result are and how they affect the rating
result (R.3). In Inspecting the Process of Bank Credit Rating via Visual
Analytics(B), the ranking differences among schemes are represented by a
parallel coordinate-like design, in which each axis represents a ranking
scheme and each dot represents a bank. The black dot (dragged bank) is the
item dragged by the analyst. Red dots indicate the negative samples while blue
ones are the positive samples. We arrange the ranking sequentially (the
smaller the better) from left to right and use lines to link all the identical
banks across different schemes. For two adjacent ranking schemes, blue lines
indicate an increasing ranking, and red lines indicate a decreasing ranking.
We encode the bank types as background rectangles with different colors and
discretize the ranking into $5$ ratings. On the right side of Inspecting the
Process of Bank Credit Rating via Visual Analytics(B), we use box plots to
show indicator distribution of all positive and negative pairs. The color of
each box plot represents one indicator. For one indicator, the left box plot
indicates the indicator value distribution of negative pairs and the right box
plot shows the indicator value distribution of positive pairs. The curve shows
the attribute weight.

##  5 Case Study

Adjusting rankings for some banks. The experts first observed that Beijing
Rural Commercial Bank belonging to “Rural Commercial Bank” in Default Scheme
ranks higher with the rating of $2$ (Figure 1(A)). Bank of Communications
belongs to “Large State-owned Commercial Bank” but has a lower rating of $2$.
From Figure 1(B), the experts witnessed that although the asset size of
Beijing Rural Commercial Bank is low ($8,811$), the provision coverage is high
($1068.87\%$), and that is why it is relatively high in the ranking. On the
contrary, although the asset size of Bank of Communications is relatively high
($95,312$), its provision coverage is much lower ($173.13\%$), leading to a
low ranking. The experts commented that according to the bank’s previous
performance and their expertise, the current ranking of Beijing Rural
Commercial Bank is higher than expected, although it has excellent provision
coverage, the asset size is worse than that of the banks ranked near it.
Similarly, the ranking of Bank of Communications is also lower than expected.
“After all, Bank of Communications is one of the six most famous state-owned
banks in China,” said E.1, and “its asset size is strong”. The experts then
manually adjusted the ranking of the two banks, “I will drag the bank to a
ranking that meets my psychological expectation and also consider the bank
type.” E.1 witnessed that the banks ranked $15$ and $16$ are “Rural Commercial
Bank” which belong to the same bank type with Beijing Rural Commercial Bank.
Meanwhile, E.1 thought that Beijing Rural Commercial Bank is better than them
(banks at the position of $15$ and $16$) and is also better than Guangzhou
Bank (at the position of $14$) in terms of asset size. Therefore, E.1 decided
to move Beijing Rural Commercial Bank from its current position of $5$ to $13$
and Bank of Communications from its current position of $8$ to $5$. They then
clicked on the “Save Weight Scheme” button and the system generated a new
weight set, which is applied to generate three new rankings. The adjusted
result is shown in Inspecting the Process of Bank Credit Rating via Visual
Analytics(C).

Figure 1: Beijing Rural Commercial Bank belongs to “Rural Commercial Bank” but
ranks higher with rating of $2$. Bank of Communications belongs to “Large
State-owned Commercial Bank” with a low rating.

Comparing ranking schemes. The experts found that in Type Weight some banks
with a rating of $2$ are close to the banks with a rating of $1$ but they are
not very close in the ranking tabular view (Inspecting the Process of Bank
Credit Rating via Visual Analytics(A)). With curiosity, E.1 selected Beijing
Rural Commercial Bank in Inspecting the Process of Bank Credit Rating via
Visual Analytics(C), and observed the ranking differences in different ranking
schemes and how indicator contribution affects the overall ranking in each
scheme. In Figure 2, the experts further identified that due to different
schemes of sampling positive and negative pairs, the value distribution (the
distance between the upper and lower edges of the box plot) in Local Weight
and the height of the box plot are more concentrated than those in Global
Weight, which makes sense since in Local Weight, both positive and negative
samples are selected near the dragged bank, compared with the random selection
of positive and negative samples in Global Weight. Particularly, in Local
Weight, indicators such as capital adequacy ratio are more important while in
Global Weight, indicators such as asset size, capital adequacy ratio, and
capital profit ratio are emphasized.

Figure 2: Capital adequacy ratio is more important in Local Weight. In Global
Weight, asset size, capital adequacy ratio, and capital profit ratio are
emphasized. In Type Weight, asset size and non-performing loans ratio is
mainly considered in large state-owned commercial banks and rural commercial
banks, respectively.

Identifying critical indicators for ranking different types of banks. From
Figure 3, in Default Scheme, Beijing Rural Commercial Bank belongs to Rural
Commercial Bank and it is located around large state-owned commercial banks,
which is not reasonable. From Figure 2, the experts confirmed that the
indicator weights of large state-owned commercial banks mainly lie in asset
size, while for rural commercial banks, the non-performing loans ratio
matters. E.1 explained “smaller banks tend to lend more loans with risks,
leading to a high non-performing loans ratio, while the risk management of
large banks is good and lending is conservative. So non-performing loans ratio
is smaller”. Thus, the non-performing loans ratio are not significant for
distinguishing large state-owned commercial banks, leading a smaller weight.
The experts further noticed that in the axis of Type Weight, Bank of
Communications has a rating of $1$, ranking the $5^{th}$ and Beijing Rural
Commercial Bank has a rating of $2$ with ranking $21^{st}$; however, in
Default Scheme, Bank of Communications has a rating of $2$, ranking the
$8^{th}$ and Beijing Rural Commercial Bank has a rating of $2$, ranking
$5^{th}$. This observation demonstrates that the performance of Bank of
Communications in Type Weight is improved while the performance of Beijing
Rural Commercial Bank is worse, which is consistent with the adjustment by the
experts, i.e., E.1 adjusted the ranking of Beijing Rural Commercial Bank to
the $13^{rd}$, although in Type Weight, the ranking is recommended as the
$21^{st}$. The experts also observed that the ranking of Taizhou Bank varies
among different schemes. As shown in Inspecting the Process of Bank Credit
Rating via Visual Analytics(C), E.1 found that although the performance of
Taizhou Bank in asset size is on average, it performs well in terms of
provision coverage, asset profit ratio, and capital profit ratio, leading to
overall good performance in Default Scheme (Figure 3). Similarly, the rankings
of Taizhou Bank in Local Weight and Global Weight are all highly ranked.
However, E.1 thought the rankings of Taizhou Bank were seriously overestimated
in the first three schemes and he was satisfied with the ranking in Type
Weight. On the contrary, the experts identified the rankings of Ping’An Bank
in the first three schemes were underrated (Figure 3), but “it actually has a
good credit rating,” said E.1. He commented that the rating in Type Weight
(rating of $2$) is more in line with his expectation, confirming his intuition
that “the weights of different bank types should be treated separately.”

Figure 3: Ranking comparison of four banks among Default Scheme,
Local/Global/Type Weight rankings.

##  6 Conclusion and Future Work

We propose RatingVis to facilitate domain experts exploring and comparing
different bank credit rating schemes. A ranking tabular view deduces the
indicator weights based on user interaction; a projection view visualizes the
ranking data distributions generated by different ranking schemes, and a
ranking comparison view compares the ranking schemes at an instance level. A
case study verifies the efficacy of RatingVis. In the future, we will consider
attributes with various types and explore more ranking schemes for comparison.

###### Acknowledgements.

We thank the anonymous reviewers. This work is partially supported by the
research start-up fund of ShanghaiTech and HKUST-WeBank Joint Laboratory
Project Grant No.: WEB19EG01-d.

## References

  * [1] K. M. A.-S. Al-Harbi.  Application of the ahp in project management.  International journal of project management, 19(1):19–27, 2001\.
  * [2] P. Brusov, T. Filatova, N. Orekhova, and M. Eskindarov.  Rating methodology: new look and new horizons.  In Modern Corporate Finance, Investments, Taxation and Ratings, pp. 475–510. Springer, 2018.
  * [3] R. Cantor.  Moody’s investors service response to the consultative paper issued by the basel committee on bank supervision “a new capital adequacy framework”.  Journal of Banking & Finance, 25(1):171–185, 2001.
  * [4] G. Carenini and J. Loyd.  Valuecharts: analyzing linear models expressing preferences and evaluations.  In Proceedings of the working conference on Advanced visual interfaces, pp. 150–157, 2004.
  * [5] S. Claessens, A. Law, and T. Wang.  How do credit ratings affect bank lending under capital constraints?  BIS Working Papers 747, Bank for International Settlements, Sept. 2018\.
  * [6] D. Cossock and T. Zhang.  Subset ranking using regression.  In International Conference on Computational Learning Theory, pp. 605–619. Springer, 2006.
  * [7] C. R. de Sá, C. Soares, A. Knobbe, P. Azevedo, and A. M. Jorge.  Multi-interval discretization of continuous attributes for label ranking.  In J. Fürnkranz, E. Hüllermeier, and T. Higuchi, eds., Discovery Science, pp. 155–169. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.
  * [8] S. Gratzl, A. Lex, N. Gehlenborg, H. Pfister, and M. Streit.  Lineup: Visual analysis of multi-attribute rankings.  IEEE transactions on visualization and computer graphics, 19(12):2277–2286, 2013.
  * [9] S. Gumparthi, S. Khatri, and V. Manickavasagam.  Design and development of credit rating model for public sector banks in india: Special reference to small and medium enterprises.  Journal of Accounting and Taxation, 3(3):105–124, 2011.
  * [10] N. L. Gustavo and A. Michael.  Multidimensional projection for visual analytics: Linking techniques with distortions, tasks, and layout enrichment.  IEEE Transactions on Visualization & Computer Graphics, PP:1–1, 2018.
  * [11] S. Jadhav, H. He, and K. Jenkins.  Information gain directed genetic algorithm wrapper feature selection for credit rating.  Applied Soft Computing, 69, 2018.
  * [12] T. Joachims.  Optimizing search engines using clickthrough data.  In Proceedings of the eighth ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 133–142, 2002.
  * [13] Y.-J. Lai, T.-Y. Liu, and C.-L. Hwang.  Topsis for modm.  European journal of operational research, 76(3):486–500, 1994.
  * [14] Q. Li, Q. Liu, C. Tang, Z. Li, S. Wei, X. Peng, M. Zheng, T. Chen, and Q. Yang.  WarehouseVis: A Visual Analytics Approach to Facilitating Warehouse Location Selection for Business Districts.  Computer Graphics Forum, 2020. doi: 10 . 1111/cgf . 13996
  * [15] Q. Li, K. S. Njotoprawiro, H. Haleem, Q. Chen, C. Yi, and X. Ma.  Embeddingvis: A visual analytics approach to comparative network embedding inspection.  In 2018 IEEE Conference on Visual Analytics Science and Technology (VAST), pp. 48–59. IEEE, 2018.
  * [16] C. Liu and M. Zhu.  The study and implementation of discretization algorithm based on information entropy.  Journal of Software Engineering and Application, 8(6):358–363, 2019\.
  * [17] T.-Y. Liu.  Learning to rank for information retrieval.  Springer Science & Business Media, 2011.
  * [18] L. v. d. Maaten and G. Hinton.  Visualizing data using t-sne.  Journal of machine learning research, 9(Nov):2579–2605, 2008.
  * [19] S. Pajer, M. Streit, T. Torsney-Weir, F. Spechtenhauser, T. Möller, and H. Piringer.  Weightlifter: Visual weight space exploration for multi-criteria decision making.  IEEE transactions on visualization and computer graphics, 23(1):611–620, 2016.
  * [20] P. Pasricha, D. Selvamuthu, and V. Arunachalam.  Markov regenerative credit rating model.  Journal of Risk Finance, 18(3):311–325, 2017.
  * [21] J. Peng.  Do investors look beyond insured triple–a rating? an analysis of standard & poor’s underlying ratings.  Public Budgeting & Finance, 22(3):115–131, 2002.
  * [22] A. A. Peresetsky and A. M. Karminsky.  Models for Moody’s bank ratings.  MPRA Paper 34864, University Library of Munich, Germany, 2011.
  * [23] A. Puri, B. K. Ku, Y. Wang, and H. Qu.  Rankbooster: Visual analysis of ranking predictions.  arXiv preprint arXiv:2004.06435, 2020.
  * [24] L. Rojas-Suarez.  Rating banks in emerging markets: What credit rating agencies should learn from financial indicators.  In Ratings, rating agencies and the global financial system, pp. 177–201. Springer, 2002.
  * [25] C.-H. Shen, Y.-L. Huang, and I. Hasan.  Asymmetric benchmarking in bank credit rating.  Journal of International Financial Markets, Institutions and Money, 22(1):171–193, 2012.
  * [26] Shi, Baofeng, Chen, Nan, Wang, and Jing.  A credit rating model of microfinance based on fuzzy cluster analysis and fuzzy pattern recognition: Empirical evidence from chinese 2,157 small private businesses.  Journal of intelligent & fuzzy systems: Applications in Engineering and Technology, 2016.
  * [27] C. O. S. Sorzano, J. Vargas, and A. P. Montano.  A survey of dimensionality reduction techniques.  arXiv preprint arXiv:1403.2877, 2014.
  * [28] W. F. Treacy and M. Carey.  Credit risk rating systems at large us banks.  Journal of Banking & Finance, 24(1-2):167–201, 2000.
  * [29] L. Van Der Maaten, E. Postma, and J. Van den Herik.  Dimensionality reduction: a comparative.  J Mach Learn Res, 10(66-71):13, 2009.
  * [30] E. Wall, S. Das, R. Chawla, B. Kalidindi, E. T. Brown, and A. Endert.  Podium: Ranking data using mixed-initiative visual analytics.  IEEE transactions on visualization and computer graphics, 24(1):288–297, 2017.
  * [31] S. Yuksel, H. Dincer, and U. Hacioglu.  Camels-based determinants for the credit rating of turkish deposit banks.  International Journal of Finance & Banking Studies (2147-4486), 4(4):1–17, 2015.
  * [32] S. H. Zanakis, A. Solomon, N. Wishart, and S. Dublish.  Multi-attribute decision making: A simulation comparison of select methods.  European journal of operational research, 107(3):507–529, 1998\.


# Tunable correlation effects of magnetic impurities by the cubic Rashba spin-
orbit couplings

Xiong-Tao Peng  Fang Lin  Department of Physics, Ningbo University, Ningbo
315211, China  Liang Chen  School of Mathematics and Physics, North China
Electric Power University, Beijing, 102206, People’s Republic of China  Lin Li
College of Physics and Electronic Engineering, and Center for Computational
Sciences, Sichuan Normal University, Chengdu 610068, China  Dong-Hui Xu
Department of Physics and Chongqing Key Laboratory for Strongly Coupled
Physics, Chongqing University, Chongqing 400044, China  Center of Quantum
Materials and Devices, Chongqing University, Chongqing 400044, China  Jin-Hua
Sun<EMAIL_ADDRESS>Department of Physics, Ningbo University, Ningbo
315211, China

###### Abstract

We theoretically study the influence of the $k$-cubic Rashba spin-orbit
coupling (SOC) on the correlation effects of magnetic impurities by combining
the variational method and the Hirsch-Fye quantum Monte Carlo (HFQMC)
simulations. Markedly different from the normal $k$-linear Rashba SOC, even a
small cubic Rashba term can greatly alter the band structure and induce a Van
Hove singularity in a wide range of energy, thus the single impurity local
moment becomes largely tunable. The cubic Rashba SOC adopted in this work
breaks the rotational symmetry, but the host material is still invariant under
the operations $\mathcal{R}^{z}(\pi)$, $\mathcal{IR}^{z}(\pi/2)$,
$\mathcal{M}_{xz}$, $\mathcal{M}_{yz}$, where $\mathcal{R}^{z}(\theta)$ is the
rotation of angle $\theta$ about the $z$-axis, $\mathcal{I}$ is the inversion
operator and $\mathcal{M}_{xz}$ ($\mathcal{M}_{yz}$) is the mirror reflection
about the $x$-$z$ ($y$-$z$) principal plane. Saliently, various components of
spin-spin correlation between the single magnetic impurity and the conduction
electrons show three- or six-fold rotational symmetry. This unique feature is
due to the triple winding of the spins with a $2\pi$ rotation of $\mathbf{k}$,
which is a hallmark of the cubic Rashba effect, and can possibly be an
identifier to distinguish the cubic Rashba SOC from the normal $k$-linear
Rashba term in experiments. Although the cubic Rashba term drastically alters
the electronic properties of the host, we find that the spatial decay rate of
the spin-spin correlation function remains essentially unchanged. Moreover,
the carrier-mediated Ruderman–Kittel–Kasuya–Yosida interactions between two
magnetic impurities show twisted features, the ferromagnetic diagonal terms
dominate when two magnetic impurities are very close, but the off-diagonal
terms become important at long distances.

##  I Introduction

Spin-orbit coupling (SOC) is a relativistic effect that locks the spin of a
charge carrier with its angular momentum, and intense efforts have been made
over the past decades to investigate and utilize SOCs in condensed matter
physics. There exist two representative SOCs, namely the Dresselhaus SOC
caused by the bulk inversion asymmetry,Dresselhaus (1955) and the Rashba SOC
due to the spatial inversion asymmetry.Rashba (1960); J. Ohkawa and Uemura
(1974) In low-dimensional systems, the Rashba SOC becomes more important
because it is stronger in the heterointerface,Vas’ ko (1979); Bychkov and
Rashba (1984) and it is often described by the $k$-linear Rashba term, which
can be written as $\propto(k_{-}\sigma_{+}-k_{+}\sigma_{-})$, where
$k_{\pm}=k_{x}\pm ik_{y}$ denote the wave vectors and
$\sigma_{\pm}=\sigma_{x}\pm i\sigma_{y}$ are the spin Pauli matrices.Bihlmayer
_et al._ (2015); Manchon  _et al._ (2015); Bihlmayer  _et al._ (2022)

Besides the normal $k$-linear Rashba SOC, there also exists a higher-order
term, namely the $k$-cubic Rashba SOC which has received continuous attention
these years.Gerchikov and Subashiev (1992); Winkler (2003, 2000); Zhao  _et
al._ (2020) The cubic Rashba SOC, which is often described by the Hamiltonian
$\propto(k_{-}^{3}\sigma_{+}-k_{+}^{3}\sigma_{-})$, can greatly alter the
dispersion relation and the effective field symmetry, and is predicted to
induce larger spin Hall conductivity.Schliemann and Loss (2005); Bleibaum and
Wachsmuth (2006); Ma and Liu (2006) The cubic Rashba SOC has been reported in
a two-dimensional hole gas in inversion asymmetric semiconductors InGaAs and
GaAs heterostructures,Minkov  _et al._ (2005); Winkler  _et al._ (2002) and a
quasi- two-dimensional electron gas formed at a surface of SrTiO3 single
crystal,Nakamura  _et al._ (2012) and in rare-earth ternary materials
TbRh2Si2Usachov  _et al._ (2020a) and EuIr2Si2.Usachov  _et al._ (2020b)

As a prototypical strong correlation problem, Kondo effect in normal metals
has been widely studied and well understood.Anderson (1961); Kondo (1964);
Wilson (1975) The Kondo effect is accompanied by the formation of Kondo cloud,
which is characterized by the antiferromagnetic spin-spin correlation between
the magnetic impurity and the conduction electrons. This spin-spin correlation
function oscillates fast in space, and decays as $\sim 1/r^{D}$ when
$r<\xi_{K}$, while it decays as $\sim 1/r^{D+1}$ if $r>\xi_{K}$,Ishii (1978);
Barzykin and Affleck (1998); Borda (2007) where $\xi_{K}$ is the Kondo length
that extends to $\sim 1\ \mu m$ in typical metals,Moca  _et al._ (2021) and
has been confirmed recently via Fabry–P$\acute{e}$rot oscillations in
conductance.V Borzenets  _et al._ (2020)

The influence of $k$-linear Rashba SOC on the Kondo temperature $T_{k}$ has
been studied previously using various methods, some indicate that $T_{k}$ is
not significantly changed by Rashba SOC,Malecki (2007); Žitko and Bonča
(2011); Isaev  _et al._ (2012) while others claim an exponential enhancement
of $T_{K}$.Zarea  _et al._ (2012); Chen  _et al._ (2017) Later, numerical
renormalization group study found that for a fixed Fermi energy, the Kondo
temperature $T_{K}$ varies weakly with Rashba SOC. If instead, the band
filling is low and held constant, increasing the Rashba SOC can drive the
system into a helical regime where $T_{K}$ is exponentially enhanced.Wong  _et
al._ (2016) Basically, one important reason to change the Kondo temperature is
the divergence of density of states (DOS) which appears close to the band edge
in the presence of $k$-linear Rashba SOC.Chen  _et al._ (2016) On the other
hand, in two-dimensional superconductors, it is found that $T_{K}$ is
determined by the interplay between the Rashba SOC and superconducting energy
gap, that the quantum phase transition between the magnetic doublet and Kondo
singlet ground states is significantly affected by the Rashba SOC.Li  _et al._
(2018)

Moreover, the Kondo screening cloud shows anisotropy in both spatial and spin
spaces in the presence of SOC.Feng  _et al._ (2010); Feng and Zhang (2011) It
has also been proposed to use a magnetic impurity as a way to detect the
Rashba effect through the local magnetization density of states. Chirla  _et
al._ (2013) Taking into account of the indirect exchange couplings between
magnetic impurities, the Ruderman–Kittel–Kasuya–Yosida (RKKY)Kasuya (1956);
Ruderman and Kittel (1954); Yosida (1957) couplings become twisted in the
presence of SOC. The RKKY interaction in two-dimensional systems with SOC can
be written in a general form with three terms: Heisenberg, Ising, and
Dzyaloshinskii-Moriya (DM) interactions, and this general form is valid for
the Rashba SOC, the Dresselhaus SOC, and even when the two types of SOC are
mixed.Imamura  _et al._ (2004); Mross and Johannesson (2009); Zhu  _et al._
(2011)

In this paper, we combine the variational method and the Hirsch-Fye quantum
Monte Carlo (HFQMC)Hirsch and Fye (1986) simulations to study the correlation
effects of the impurities induced by the cubic Rashba SOC. The variational
method has been widely used in the ground states of Anderson impurity problems
in normal metals,Gunnarsson and Schönhammer (1983); Varma and Yafet (1976)
systems with SOCs,Feng  _et al._ (2010); Sun  _et al._ (2015, 2018); Ma  _et
al._ (2018); Wang  _et al._ (2019); Yang  _et al._ (2021) and
superconductors.Simonin and Allub (1995); Simon and Varma (1999); Rozhkov and
Arovas (2000); Huang  _et al._ (2022) The HFQMC technique is a numerically
exact method which has been used to study magnetic impurities in metals,
Hirsch and Fye (1986); Haldane and Anderson (1976); Uchoa  _et al._ (2008);
Fye and Hirsch (1988); Fye  _et al._ (1987); Hirsch and Lin (1987) dilute
magnetic semiconductors,Bulut  _et al._ (2007) graphene based systemsHu  _et
al._ (2011); Sun  _et al._ (2013a, b); Sun and Tang (2018) and in the presence
of SOCs.Sun  _et al._ (2014); Hu  _et al._ (2013) By combining the two
methods, we can obtain not only a heuristic physical picture, but also the
numerically exact results about the correlations. The rest of the paper is
organized as follows. In Sec. II we introduce the model Hamiltonian and
discuss the influence of the cubic Rashba term on the electronic properties of
the host material. In Sec. III, we show the results obtained using the
variational method and the Hirsch-Fye quantum Monte Carlo simulations for
single impurity case. The spin-spin correlation between two magnetic atoms,
which is mediated by the conduction electrons, are given in Sec. IV. Finally,
discussions and conclusions are given in Sec. V.

Figure 1: (Color online) (a) The two-fold degenerate band (black line) splits
into two bands due to the cubic Rashba SOC $\alpha$. The cubic Rashba term has
marginal influence on the upper band (blue line). The red lines correspond to
the lower bands with different $\alpha$ values, and they are drastically
modified by $\alpha$. $\mu$ is the chemical potential, and $\varepsilon_{d}$
is the magnetic impurity energy level. (b) The density of states varies with
$\alpha$, and Van Hove singularity emerges. (c) Spin texture caused by the
cubic Rashba SOC in momentum space. (d) Schematic of the spin-spin interaction
as a function of the distance $R$ between two magnetic impurities. The red
curve means the correlation strength, and the arrows denote the rotation of
spin-spin interaction due to the cubic Rashba SOC.

##  II Model Hamiltonian

We use the Anderson impurity model to study the properties of magnetic
impurities in a system with the cubic Rashba SOC term, the total Hamiltonian
is given by

$\displaystyle H=H_{0}+H_{d}+H_{V}.$ (1)

$H_{0}$ describes the host material with the cubic Rashba SOC, $H_{d}$ is the
magnetic impurity part, and $H_{V}$ denotes the hybridization between the
local impurities and the conduction electrons. The low-energy effective
Hamiltonian of a host system with the cubic Rashba SOC is given by

$\displaystyle
H_{0}=\sum_{\mathbf{k}}c_{\mathbf{k}}^{\dagger}\left[h_{0}(\mathbf{k})-\mu\right]c_{\mathbf{k}},$
(2)

with

$\displaystyle
h_{0}(\mathbf{k})=\frac{\hbar^{2}\mathbf{k}^{2}}{2m}+\frac{i\alpha}{2}(k_{-}^{3}\sigma_{+}-k_{+}^{3}\sigma_{-}).$
(3)

$h_{0}(\mathbf{k})$ is the single particle Hamiltonian incorporating cubic
Rashba SOC,Schliemann and Loss (2005); Moriya  _et al._ (2014); Zarea and
Ulloa (2006) and
$c_{\mathbf{k}}^{\dagger}=(c_{\mathbf{k}\uparrow}^{\dagger},c_{\mathbf{k}\downarrow}^{\dagger})$
is the creation operator in spinor representation. The notations
$k_{\pm}=k_{x}\pm ik_{y}$, $\sigma_{\pm}=\sigma_{x}\pm i\sigma_{y}$ are used
to denote the wave vectors and Pauli spin matrices. $\mu$ is the chemical
potential, and $\alpha$ is the cubic Rashba term which can be adjusted
experimentally. Liu  _et al._ (2018); Xiao  _et al._ (2019); Liu and Lei
(2005); Karwacki  _et al._ (2018) Due to the SOC, the single particle
eigenenergy splits from simple degenerate parabolic bands to two branches,

$\displaystyle\epsilon_{k\pm}=\frac{\hbar^{2}k^{2}}{2m}\pm\alpha k^{3}.$ (4)

The magnetic impurity part is given by

$\displaystyle
H_{d}=\sum_{j,s=\uparrow,\downarrow}(\epsilon_{d}-\mu)d_{js}^{\dagger}d_{js}+\sum_{j}Ud_{j\uparrow}^{\dagger}d_{j\uparrow}d_{j\downarrow}^{\dagger}d_{j\downarrow}.$
(5)

$j$ represents the magnetic impurity index. We study two cases, namely the
single impurity doping and the two-impurity doping. When only one magnetic
impurity is doped in the host, $j=1$. Otherwise if two impurities exist,
$j=1,2$. $d_{s}^{\dagger}$ and $d_{s}$ are the creation and annihilation
operators of the spin-$s$ ($s=\uparrow,\downarrow$) state on the impurity
site. $\epsilon_{d}$ is the impurity energy level which is beneath $\mu$ in
our calculations, and $U$ is the on-site Coulomb repulsion.

Finally, the hybridization term between the localized state and the conduction
electrons reads

$\displaystyle H_{V}=\sum_{\mathbf{k},j,\atop
s=\uparrow,\downarrow}\left(e^{i\mathbf{k}\cdot\mathbf{R_{j}}}V_{k}d_{js}^{\dagger}c_{\mathbf{k}s}+H.c.\right).$
(6)

$V_{k}$ is the hybridization strength, and $\mathbf{R}_{j}$ is the coordinate
of the $j$-th impurity. For two-impurity doping, we assume the two local atoms
and conduction electrons have the same exchange coupling strength $V_{k}$ for
simplification.

In Fig. 1(a), we show the dispersion relation of the single particle energy
bands given in Eq. (4). The cubic Rashba SOC splits the degenerate parabolic
band (black solid line) into two branches. One is $\epsilon_{k+}$ (blue solid
line) and the other one is $\epsilon_{k-}$ (red lines). The cubic Rashba SOC
term has minor effect on $\epsilon_{k+}$, but it can alter $\epsilon_{k-}$
significantly, as we can see from the three red lines, which correspond to
$\epsilon_{k-}$ for slightly different $\alpha$ values. The DOS for different
$\alpha$ values are plotted in Fig. 1(b). The cubic Rashba term largely
modifies $\epsilon_{k-}$, consequently induce Van Hove singularity (VHS) into
the host system, which is expected to greatly influences the local moment
formation of magnetic impurities. Besides, the cubic Rashba SOC also breaks
the rotational symmetry, but the system remains invariant under operations
such as $\mathcal{R}^{z}(\pi)$, $\mathcal{IR}^{z}(\pi/2)$, $\mathcal{M}_{xz}$,
$\mathcal{M}_{yz}$, where $\mathcal{R}^{z}(\theta)$ is the rotation of angle
$\theta$ about the $z$-axis, $\mathcal{I}$ is the inversion operator and
$\mathcal{M}_{xz}$ ($\mathcal{M}_{yz}$) is the mirror reflection about the
$x$-$z$ ($y$-$z$) principal plane. The spin texture given in Fig. 1(c)
reflects all these symmetries, which can be exhibited by the Kondo effect.
Given in Fig. 1(d) is the schematic of our two-impurity case calculation. One
magnetic impurity is fixed at the origin, and the other is located at a
distance $R$ along the $x$-axis. In our calculations, the length unit is
chosen as $k_{0}^{-1}$ which in typical metal is $k_{0}^{-1}\approx 10^{-9}m$.
Correspondingly, the energy unit is $\frac{\hbar^{2}k_{0}^{2}}{2m^{*}}\approx
1.8\times 10^{-2}eV$ and the values of parameters $\alpha$, $U$, $\mu$,
$V_{k}$, $\epsilon_{d}$ are given in units of
$\frac{\hbar^{2}k_{0}^{2}}{2m^{*}}$.Feng and Zhang (2011)

##  III Single impurity correlation effects

###  III.1 The variational method

We can easily diagonalize $H_{0}$ and obtain a quasiparticle operator

$\displaystyle\gamma_{\mathbf{k}\pm}$
$\displaystyle=\dfrac{1}{\sqrt{2}}\left(e^{i\frac{3}{2}\theta_{\mathbf{k}}}c_{\mathbf{k}\uparrow}\pm
ie^{-i\frac{3}{2}\theta_{\mathbf{k}}}c_{\mathbf{k}\downarrow}\right),$ (7)

where $\tan\theta_{\mathbf{k}}=k_{y}/k_{x}$, $\pm$ denotes the upper and lower
energy bands. First let’s discuss the simplest case when $H_{V}=0$ in which
the magnetic impurity state decouples from the host material. Thus the ground
state wave function of $H_{0}$ is given by

$\displaystyle|\Psi_{0}\rangle=\prod_{\\{\mathbf{k}\pm\\}\in\Omega}\
\gamma_{\mathbf{k}\pm}^{\dagger}|0\rangle,$ (8)

where $|0\rangle$ is the vacuum, and the product runs over all the states
within the Fermi sea $\Omega$. As for the impurity part, we assume that the
Coulomb repulsion $U$ is large enough, and the impurity energy level
$\epsilon_{d}$ is below the chemical potential $\mu$, so that the impurity
site is always singly occupied by a local moment. The total energy of the
system under this decoupled case is

$E_{0}=\epsilon_{d}-\mu+\sum_{\\{\mathbf{k}\pm\\}}\left(\epsilon_{\mathbf{k}\pm}-\mu\right).$
(9)

Then we consider the case with hybridization, where the trial wave function of
the ground state is

$\displaystyle|\Psi\rangle=\left(a_{0}+\sum_{\\{\mathbf{k}\pm\\}}a_{\mathbf{k}\pm}d_{\mathbf{k}\pm}^{\dagger}\gamma_{\mathbf{k}\pm}\right)|\Psi_{0}\rangle,$
(10)

where
$d_{\mathbf{k}\pm}=\dfrac{1}{\sqrt{2}}\left(e^{i\frac{3}{2}\theta_{\mathbf{k}}}d_{\uparrow}\pm
ie^{-i\frac{3}{2}\theta_{\mathbf{k}}}d_{\downarrow}\right)$. $a_{0}$ and
$a_{\mathbf{k}\pm}$ are variational parameters to be determined by optimizing
the ground state energy. The energy of the system in the trial state
$|\Psi\rangle$ is given by

$\displaystyle E$
$\displaystyle=\dfrac{\langle\Psi|H|\Psi\rangle}{\langle\Psi|\Psi\rangle}$
(11)
$\displaystyle=\dfrac{\sum_{\\{\mathbf{k}\pm\\}}\left(E_{0}-\epsilon_{\mathbf{k}\pm}+\mu\right)a_{\mathbf{k}\pm}^{2}+2V_{k}a_{0}a_{\mathbf{k}\pm}+\left(\epsilon_{\mathbf{k}\pm}-\mu\right)a_{0}^{2}}{a_{0}^{2}+\sum_{\\{\mathbf{k}\pm\\}}a_{\mathbf{k}\pm}^{2}}.$

The variational method requires $\partial E/\partial a_{0}=\partial E/\partial
a_{\mathbf{k}\pm}=0$, leading to

$\displaystyle\left(E-\sum_{\\{\mathbf{k}\pm\\}}\left(\epsilon_{\mathbf{k}\pm}-\mu\right)\right)a_{0}$
$\displaystyle=\sum_{\\{\mathbf{k}\pm\\}}V_{k}a_{\mathbf{k}\pm},$ (12)
$\displaystyle\left(E-E_{0}+\epsilon_{\mathbf{k}\pm}-\mu\right)a_{\mathbf{k}\pm}$
$\displaystyle=V_{k}a_{0}.$

We can define the binding energy as $\Delta_{b}=E_{0}-E$, then the self-
consistent equation is given by

$\displaystyle\left(\epsilon_{d}-\mu\right)=\Delta_{b}=\sum_{\\{\mathbf{k}\pm\\}}\dfrac{\left|V_{k}\right|^{2}}{\epsilon_{\mathbf{k}\pm}-\mu-\Delta_{b}}.$
(13)

If $\Delta_{b}>0$, the hybridized state is stable against the decoupled state.

Figure 2: (Color online) Self-consistent results of the binding energy
$\Delta_{b}$ for different combinations of parameter values. $\Delta_{b}$
versus $\alpha$ for (a) different $V_{k}$ when $\mu=1.85$ and (b) different
$\mu$ when $V_{k}=0.20$. (c) $\Delta_{b}$ versus $V_{k}$ for various
combinations of $\mu$, $\alpha$. $\Gamma$ is the energy cutoff chosen to be
far away from $\mu$.

In our variational method calculations, the impurity energy level is fixed
slightly below the chemical potential, $\epsilon_{d}=\mu-0.001$, and the
energy cutoff $\Gamma$ is chosen to be far away from $\mu$, that the low-
energy physical properties will not be affected by the choice of $\Gamma$.

We show the self-consistent results of the binding energy $\Delta_{b}$ for
various combinations of $V_{k}$, $\mu$ and $\alpha$ in Fig. 2. In Fig. 2(a) we
show the binding energy with respect to the cubic Rashba term $\alpha$ for
different values of $V_{k}$ when $\mu=1.85$. $\alpha$ greatly alters the band
structure thus the DOS as is shown in Fig. 1. We find that $\Delta_{b}$ shows
a peak around $\alpha=0.1$. This is because for $\alpha=0.1$, the VHS occurs
close to the chemical potential $\mu=1.85$. Besides, $\Delta_{b}$ is larger
for stronger hybridization strength $V_{k}$, implying that the bound state is
more easily formed for strong $V_{k}$. In Fig. 2(b), we fix $V_{k}=0.2$ and
show similar results for various values of $\mu$. The energy corresponds to
the VHS decrease monotonically with $\alpha$, and in a wide range of $\mu$ we
can always observe the peak of $\Delta_{b}$. Fig. 2(c) shows the results of
$\Delta_{b}$ versus $V_{k}$ for different combinations of $\mu$ and $\alpha$.
In general, $\Delta_{b}$ is always positive due to the finite DOS in this
system, and this is consistent with the previous results obtained using the
same method.Feng  _et al._ (2010); Sun  _et al._ (2015, 2018) Larger values of
$\Delta_{b}$ imply that the bound state is more stable. When $\alpha=0.10$,
the VHS lies around $\mu=1.85$, such that the binding energy $\Delta_{b}$ (the
blue straight line) becomes much larger than other cases.

Next, we study the effect of the cubic Rashba SOC on the correlation between
the local spin and the conduction electrons spins. This spin-spin correlation
function measures the spatial Kondo screening cloud. The spin operator of the
magnetic impurity spin is defined as
$\mathbf{S}_{d}=\dfrac{1}{2}\sum_{s,s^{\prime}}d_{s}^{\dagger}(\hat{\bm{\sigma}})_{s,s^{\prime}}d_{s^{\prime}}$
and the conduction electron spin is
$\mathbf{S}_{c}(\mathbf{r})=\dfrac{1}{2}\sum_{s,s^{\prime}}c^{\dagger}_{s}(\mathbf{r})(\hat{\bm{\sigma}})_{s,s^{\prime}}c_{s^{\prime}}(\mathbf{r})$,
where $s,s^{\prime}=\uparrow,\downarrow$. By assuming the magnetic impurity
location as the origin $\mathbf{r}=0$ and the conduction electron position as
$\mathbf{r}$, the spin-spin correlation function is given by

$\displaystyle J_{uv}\left(\mathbf{r}\right)=\langle
S_{c}^{u}(\mathbf{r})S_{d}^{v}\left(0\right)\rangle,$ (14)

where $\langle\dots\rangle$ is the ground state average, and $u,v=x,y,z$ are
the spin indices. $J_{uv}\left(\mathbf{r}\right)$ can be calculated by using
the trial wave-function in Eq. 10.

Figure 3: (Color online) The spatial pattern of the spin-spin correlation
between magnetic impurity and conduction electrons with $(\alpha=0.1)$.
(a)-(f) are for $J_{zz}(\mathbf{r})$, $J_{xx}(\mathbf{r})$,
$J_{yy}(\mathbf{r})$, $J_{xy}(\mathbf{r})$, $J_{xz}(\mathbf{r})$,
$J_{yz}(\mathbf{r})$. The parameters are $\Delta_{b}=0.02$, $\alpha=0.10$,
$V_{k}=1.0$ and $\mu=1.50$. $k_{c}$ is the momentum cutoff chosen with respect
to the energy cutoff $\Gamma$. Figure 4: (Color online) Diagonal components of
the Spin-spin correlation functions $J_{uu}\left(\mathbf{r}\right)$. The
subfigure shows that the spin–spin correlation decays as $\propto
1/\mathbf{r}^{2}$ along the $x$-axis. The parameters are the same as those in
Fig. 3.

The diagonal and the off-diagonal terms take the form

$\displaystyle J_{zz}\left(\mathbf{r}\right)$
$\displaystyle=-\dfrac{1}{8}\left|B(\mathbf{r})\right|^{2}+\dfrac{1}{8}\left|A(\mathbf{r})\right|^{2},$
(15) $\displaystyle J_{xx}\left(\mathbf{r}\right)$
$\displaystyle=-\dfrac{1}{8}\left|B(\mathbf{r})\right|^{2}-\dfrac{1}{8}\operatorname{Re}\left[A(\mathbf{r})\right]^{2},$
$\displaystyle J_{yy}\left(\mathbf{r}\right)$
$\displaystyle=-\dfrac{1}{8}\left|B(\mathbf{r})\right|^{2}+\dfrac{1}{8}\operatorname{Re}\left[A(\mathbf{r})\right]^{2},$
$\displaystyle J_{xy}\left(\mathbf{r}\right)$
$\displaystyle=-\dfrac{1}{8}\operatorname{Im}\left[A(\mathbf{r})\right]^{2},$
$\displaystyle J_{xz}\left(\mathbf{r}\right)$
$\displaystyle=\dfrac{1}{4}\operatorname{Im}\left[B(\mathbf{r})A(\mathbf{r})\right],$
$\displaystyle J_{yz}\left(\mathbf{r}\right)$
$\displaystyle=-\dfrac{1}{4}\operatorname{Re}\left[B(\mathbf{r})A(\mathbf{r})\right],$

where $A(\mathbf{r})=\sum_{\\{\mathbf{k}\pm\\}}\pm
e^{i(\mathbf{k}\cdot\mathbf{r}+3\theta_{\mathbf{k}})}a_{\mathbf{k}\pm}$ and
$B(\mathbf{r})=\sum_{\\{\mathbf{k}\pm\\}}e^{i\mathbf{k}\cdot\mathbf{r}}a_{\mathbf{k}\pm}$.
Due to the phase factors of $\gamma_{\mathbf{k}\pm}$ given in Eq. 7,
$A(\mathbf{r})$ contains the phase factor $3\theta_{\mathbf{k}}$ thus becomes
three-fold rotational symmetric about the $z$-direction while $B(\mathbf{r})$
is isotropic in the $x$-$y$ plane.

In Fig. 3, we plot the spatial patterns of the spin-spin correlation function
$J_{uv}(\mathbf{r})$ ($u,v=x,y,z$), and $k_{c}$ is the momentum cutoff chosen
with respect to the energy cutoff $\Gamma$. $J_{zz}\left(\mathbf{r}\right)$
given in Fig. 3(a) is always isotropic about the origin while
$J_{xx}\left(\mathbf{r}\right)$ and $J_{yy}\left(\mathbf{r}\right)$ given in
Fig. 3(b), (c) are anisotropic because of the SOC in the $x$-$y$ plane. We
find that the diagonal components $J_{xx}\left(\mathbf{r}\right)$ and
$J_{yy}\left(\mathbf{r}\right)$ have three-fold rotational symmetry about the
$z$-direction. Note that the host system given in Eq. 2 is not three-fold
rotational symmetric. However, due to the phase factors of eigenstates given
in Eq. 7, the components of spin-spin correlation function show unique
rotational symmetry. In addition, the host system is invariant under
$\mathcal{R}^{z}(\pi)$, and consequently $J_{xx}\left(\mathbf{r}\right)$ and
$J_{yy}\left(\mathbf{r}\right)$ also satisfy the six-fold rotational symmetry.
All the diagonal terms are negative around $r=0$, indicating the
antiferromagnetic coupling between the magnetic impurity spin and the
conduction electron spins. The off-diagonal terms are merely induced by the
SOC, and we find that
$J_{xy}\left(\mathbf{r}\right)=J_{xy}\left(\mathbf{-r}\right)$ in Fig. 3(d),
which can be analyzed using the $\mathcal{R}^{z}(\pi)$ symmetry of the host
material. The other two off-diagonal components have the property
$J_{xz}\left(x,y\right)=-J_{xz}\left(-x,y\right)$, and
$J_{yz}\left(x,y\right)=-J_{yz}\left(x,-y\right)$. Except for the isotropic
$J_{zz}\left(\mathbf{r}\right)$, all the components of spin-spin correlation
show either three- or six-fold rotational symmetry on the $x$-$y$ plane. The
underlying reason for these unique symmetries is the triple winding of the
spins with a complete $2\pi$ rotation of $\mathbf{k}$,Usachov  _et al._
(2020a); Manchon  _et al._ (2015) which is a hallmark of the cubic Rashba
effect, and can possibly be an identifier to distinguish the cubic Rashba SOC
from the normal $k$-linear Rashba term in experiments.

All the components of the spin-spin correlation function oscillate and decay
in space. To analyze the spatial decay rate of the correlations, in Fig. 4 we
show the diagonal components of the spin-spin correlation function along the
$x$-axis. The parameters are chosen as $\Delta_{b}=0.02$, $\alpha=0.10$,
$V_{k}=1.0$ and $\mu=1.50$.
$J_{zz}\left(\mathbf{r}\right)=J_{xx}\left(\mathbf{r}\right)\neq
J_{yy}\left(\mathbf{r}\right)$ along the $x$-axis. Shown in the subfigure is
the results of $r^{2}J_{uu}(\mathbf{r})$ ($u=x,y$) along the $x$-axis.
According to previous studies, the spin-spin correlation between the magnetic
impurity and the conduction electrons decays as $\approx 1/r^{D+1}$ if
$r>\xi_{K}$,Ishii (1978); Barzykin and Affleck (1998); Borda (2007) where
$\xi_{K}$ is the Kondo length. However, our variational calculations support a
$1/{r}^{2}$ decay for finite $\alpha$ at long distances. Even for the simple
two-dimensional electron gas with $\alpha=0$, the decay rate of the spin-spin
correlation function is still proportional to $1/{r}^{2}$ unless
$\Delta_{b}>0.2$, which is unrealistically larger than the results of
$\Delta_{b}$ obtained in Fig. 2. We presume this is caused by the limitation
of the variational method, and it is necessary to perform the unbiased HFQMC
simulations to get more accurate results.

###  III.2 Quantum Monte Carlo simulations

The Hirsch-Fye algorithm naturally returns the imaginary-time Green’s
functions $g_{dd}^{ss^{\prime}}(\tau)=-\langle
T_{\tau}d_{s}(\tau)d_{s^{\prime}}^{{\dagger}}\rangle$,
$G_{cd}^{ss^{\prime}}(\mathbf{r},\tau)=-\langle
T_{\tau}c_{\mathbf{r}s}(\tau)d_{s^{\prime}}^{{\dagger}}\rangle$,
$G_{dc}^{ss^{\prime}}(\mathbf{r},\tau)=-\langle
T_{\tau}d_{s}(\tau)c_{\mathbf{r}s^{\prime}}^{{\dagger}}\rangle$,
$G_{cc}^{ss^{\prime}}(\mathbf{r},\tau)=-\langle
T_{\tau}c_{\mathbf{r}s}(\tau)c_{\mathbf{r}s^{\prime}}^{{\dagger}}\rangle$,
where $s,s^{\prime}=\uparrow,\downarrow$. In the HFQMC simulations,
$\langle\cdots\rangle$ means taking the average over the discrete auxiliary
field. $\tau$ is the imaginary time ranges from 0 to $\beta$. All the
information about the host material is included in the input non-interacting
Green’s functions $(U=0)$ which can be obtained analytically.

Figure 5: (Color online) The HFQMC results of (a) $n_{d}$, (b)
$n_{d_{\uparrow}}n_{d_{\downarrow}}$ , (c) $m_{d}^{2}$, and (d) $T\chi$ for
various values of the chemical potential $\mu$ and the cubic Rashba term
$\alpha$. We choose $U=0.8$, $V_{k}=1.0$ and the temperature is $K_{B}T=1/32$.

By using the Green’s function returned from the HFQMC simulations, we can
calculate various quantities such as the expectation values of the total
charge:

$\displaystyle n_{d}=\langle n_{d\uparrow}+n_{d\downarrow}\rangle,$

the local moment squared:

$\displaystyle m_{d}^{2}=\langle(n_{d\uparrow}-n_{d\downarrow})^{2}\rangle,$

the double occupancy:

$\displaystyle n_{d\uparrow\downarrow}=\langle
n_{d\uparrow}n_{d\downarrow}\rangle,$

and the spin susceptibility:

$\displaystyle\chi=\int_{0}^{\beta}d\tau\langle(n_{d\uparrow}(\tau)-n_{d\downarrow}(\tau))(n_{d\uparrow}(0)-n_{d\downarrow}(0))\rangle.$

$\beta=1/k_{B}T$ is the inverse temperature. Note that the local moment
squared on the impurity site is given by
$m_{d}^{2}=n_{d}-2n_{d\uparrow\downarrow}$, the closer this value is to one,
the more fully developed is the local moment. In all our QMC simulations, we
fix $\epsilon_{d}-\mu=-U/2$, namely the symmetric case in which the local
moment formation is favored.Anderson (1961)

Figure 6: (Color online) The HFQMC results of local moment squared for (a)
$U=3.0$ with various $V_{k}$ values and (b) $V_{k}=1.0$ for different $U$
values. The chemical potential is $\mu=1.0$ and the temperature is
$K_{B}T=1/32$.

In Fig. 5, we show the thermodynamic quantities with respect to the chemical
potential $\mu$ and the strength of the cubic Rashba term $\alpha$. The
parameters are chosen as $U=0.8$, $V_{k}=1.0$ and the temperature is
$K_{B}T=1/32$. The results for different parameter values shall remain
qualitatively unchanged. As is given in Fig. 1, small values of $\alpha$ can
drastically modify the dispersion relation thus induce VHS. The energy
corresponds to the VHS decreases as $\alpha$ increases, and this will
influence the single magnetic impurity local moment. In Fig. 5(a) we can see
that the occupation number on the impurity site becomes larger in some
regions, which corresponds to the case that the chemical potential is around
the energy where VHS occurs. The double occupancy given in Fig. 5(b) shows
similar behavior, and the local moment is determined by the competition
between the occupation and the double occupancy. We can see that the local
moment plotted in Fig. 5(c) becomes smaller in the same region. It is natural
that if the DOS at $\mu$ is large, the screening of the local magnetic
impurity spin becomes stronger, so the local moment as well as the spin
susceptibility shown in in Fig. 5(d) are suppressed.

In order to check the tunability of local moment by $\alpha$, we show the
results of $m_{d}^{2}$ for different combinations of $U$ and $V_{k}$ in Fig.
6. The chemical potential is fixed at $\mu=1.0$ and the temperature is
$K_{B}T=1/32$. In Fig. 6(a) we choose $U=3.0$, and change the hybridization
strength $V_{k}$. For all the values of $V_{k}$, we can find a dip of
$m_{d}^{2}$ as we switch $\alpha$. The reduction of local moment is caused by
the increase of DOS due to the cubic Rashba term $\alpha$. The change in
$m_{d}^{2}$ is more obvious if $V_{k}$ is larger. In Fig. 6(b), we show the
local moment for different $U$ values while the hybridization is chosen as
$V_{k}=1.0$. We can still see a dip of $m_{d}^{2}$ with as $\alpha$ varies,
and the changes in local moment are more obvious for small $U$ values. In
general, for a magnetic impurity with strong $V_{k}$ and relatively weak $U$,
the local moment is largely tunable by switching the cubic Rashba SOC.

In the following, the spin-spin correlation between the local magnetic
impurity and the conduction electron is studied for different combinations of
$\alpha$, $V_{k}$ and $U$ values. The spin-spin correlation between the
magnetic impurity and conduction electron can be calculated from the Green’s
functions asGubernatis  _et al._ (1987):

$\displaystyle J_{zz}(\mathbf{r})=\langle S_{d}^{z}S_{c}^{z}\rangle$
$\displaystyle=$
$\displaystyle\langle(g_{dd}^{\uparrow\uparrow}-g_{dd}^{\downarrow\downarrow})\times(g_{cc}^{\uparrow\uparrow}-g_{cc}^{\downarrow\downarrow})-g_{dc}^{\uparrow\uparrow}\cdot
g_{cd}^{\uparrow\uparrow}$ $\displaystyle-$ $\displaystyle
g_{dc}^{\downarrow\downarrow}\cdot
g_{cd}^{\downarrow\downarrow}+g_{dc}^{\uparrow\downarrow}\cdot
g_{cd}^{\downarrow\uparrow}+g_{dc}^{\downarrow\uparrow}\cdot
g_{cd}^{\uparrow\downarrow}\rangle,$ $\displaystyle J_{xx}(\mathbf{r})=\langle
S_{d}^{x}S_{c}^{x}\rangle$ $\displaystyle=$
$\displaystyle\langle(g_{dd}^{\uparrow\downarrow}+g_{dd}^{\downarrow\uparrow})\times(g_{cc}^{\uparrow\downarrow}+g_{cc}^{\downarrow\uparrow})-g_{dc}^{\uparrow\downarrow}\cdot
g_{cd}^{\uparrow\downarrow}$ $\displaystyle-$ $\displaystyle
g_{dc}^{\downarrow\uparrow}\cdot
g_{cd}^{\downarrow\uparrow}-g_{dc}^{\uparrow\uparrow}\cdot
g_{cd}^{\downarrow\downarrow}-g_{dc}^{\downarrow\downarrow}\cdot
g_{cd}^{\uparrow\uparrow}\rangle,$ $\displaystyle J_{yy}(\mathbf{r})=\langle
S_{d}^{y}S_{c}^{y}\rangle$ $\displaystyle=$
$\displaystyle\langle-(g_{dd}^{\uparrow\downarrow}-g_{dd}^{\downarrow\uparrow})\times(g_{cc}^{\uparrow\downarrow}-g_{cc}^{\downarrow\uparrow})+g_{dc}^{\uparrow\downarrow}\cdot
g_{cd}^{\uparrow\downarrow}$ $\displaystyle+$ $\displaystyle
g_{dc}^{\downarrow\uparrow}\cdot
g_{cd}^{\downarrow\uparrow}-g_{dc}^{\uparrow\uparrow}\cdot
g_{cd}^{\downarrow\downarrow}-g_{dc}^{\downarrow\downarrow}\cdot
g_{cd}^{\uparrow\uparrow}\rangle,$ $\displaystyle J_{xy}(\mathbf{r})=\langle
S_{d}^{x}S_{c}^{y}\rangle$ $\displaystyle=$ $\displaystyle\langle
i[(g_{dd}^{\uparrow\downarrow}+g_{dd}^{\downarrow\uparrow})\times(g_{cc}^{\uparrow\downarrow}-g_{cc}^{\downarrow\uparrow})-g_{dc}^{\uparrow\downarrow}\cdot
g_{cd}^{\uparrow\downarrow}$ $\displaystyle+$ $\displaystyle
g_{dc}^{\downarrow\uparrow}\cdot
g_{cd}^{\downarrow\uparrow}+g_{dc}^{\uparrow\uparrow}\cdot
g_{cd}^{\downarrow\downarrow}-g_{dc}^{\downarrow\downarrow}\cdot
g_{cd}^{\uparrow\uparrow}]\rangle,$ $\displaystyle J_{xz}(\mathbf{r})=\langle
S_{d}^{x}S_{c}^{z}\rangle$ $\displaystyle=$
$\displaystyle\langle(g_{dd}^{\uparrow\downarrow}+g_{dd}^{\downarrow\uparrow})\times(g_{cc}^{\uparrow\uparrow}-g_{cc}^{\downarrow\downarrow})+g_{dc}^{\uparrow\downarrow}\cdot
g_{cd}^{\downarrow\downarrow}$ $\displaystyle-$ $\displaystyle
g_{dc}^{\downarrow\uparrow}\cdot
g_{cd}^{\uparrow\uparrow}-g_{dc}^{\uparrow\uparrow}\cdot
g_{cd}^{\uparrow\downarrow}+g_{dc}^{\downarrow\downarrow}\cdot
g_{cd}^{\downarrow\uparrow}\rangle,$ $\displaystyle J_{yz}(\mathbf{r})=\langle
S_{d}^{y}S_{c}^{z}\rangle$ $\displaystyle=$ $\displaystyle\langle
i[(g_{dd}^{\uparrow\downarrow}-g_{dd}^{\downarrow\uparrow})\times(g_{cc}^{\uparrow\uparrow}-g_{cc}^{\downarrow\downarrow})+g_{dc}^{\downarrow\uparrow}\cdot
g_{cd}^{\uparrow\uparrow}$ $\displaystyle+$ $\displaystyle
g_{dc}^{\uparrow\downarrow}\cdot
g_{cd}^{\downarrow\downarrow}-g_{dc}^{\uparrow\uparrow}\cdot
g_{cd}^{\uparrow\downarrow}-g_{dc}^{\downarrow\downarrow}\cdot
g_{cd}^{\downarrow\uparrow}]\rangle.$

We assume that the magnetic impurity is located at the origin of the
coordinate, and $\mathbf{r}$ is the position of the conduction electron. In
Fig. 7 we show the results of $J_{zz}(\mathbf{r}=\\{1.0,0\\})$ with various
combinations of $V_{k}$, $\mu$ and $U$. Given in Fig. 7(a) are the results of
$J_{zz}(\mathbf{r})$ versus $\alpha$ for a fixed value of $U=3.0$. For all the
parameters, we see that $J_{zz}(\mathbf{r})$ becomes stronger in a region as
$\alpha$ increases. This region corresponds to the cases when the VHS emerges
around the chemical potential $\mu$. In general, the values of
$J_{zz}(\mathbf{r})$ grows with the hybridization strength $V_{k}$. In Fig.
7(b), we fix the chemical potential as $\mu=1.5$, $V_{k}=1.0$, and present the
results for different $U$. For all the values of $U$, we see similar behavior,
that is the increase of $J_{zz}(\mathbf{r})$ in a certain region of $\alpha$.
This indicates that the increase of spin-spin correlation is common for
intermediate $U$ values. Note that the relative magnitude of the
$J_{zz}(\mathbf{r})$ is not always decrease monotonically with $U$, and it
depends on the choice of $\mathbf{r}$.

Figure 7: (Color online) The HFQMC results of $J_{zz}(\mathbf{r}=\\{1.0,0\\})$
with various combinations of $V_{k}$, $\mu$ and $U$. (a) The results of
$J_{zz}(\mathbf{r})$ versus $\alpha$ for (a) fixed value of $U=3.0$, (b) for
fixed $\mu=1.5$ and $V_{k}=1.0$. The temperature is chosen as $K_{B}T=1/16$.
Figure 8: (Color online) The HFQMC results of the spin-spin correlation
between the magnetic impurity and the conduction electrons in the $x$-$y$
plane. (a)-(f) are for $J_{zz}(\mathbf{r})$, $J_{xx}(\mathbf{r})$,
$J_{yy}(\mathbf{r})$, $J_{xy}(\mathbf{r})$, $J_{xz}(\mathbf{r})$,
$J_{yz}(\mathbf{r})$. The parameters are fixed as $\mu=1.5$, $\alpha=0.1$,
$V_{k}=1.0$, $U=3.0$, and the temperature is $K_{B}T=1/16$.

In Fig. 8 we show the HFQMC results of the spin-spin correlation between the
magnetic impurity and the conduction electrons in the $x$-$y$ plane. The
parameters are fixed as $\mu=1.5$, $\alpha=0.1$, $V_{k}=1.0$ and $U=3.0$, and
$k_{c}$ is the momentum truncation. We can see that the HFQMC results of the
spin-spin correlation exhibit basically the same symmetry with those obtained
using the variational method as in Fig. 3. $J_{zz}(\mathbf{r})$ given in Fig.
8(a) is isotropic in the $x$-$y$ plane. The other two diagonal components
$J_{xx}(\mathbf{r})$ and $J_{yy}(\mathbf{r})$ given in Fig. 8(b) and (c) are
six-fold rotational symmetric, and
$J_{xx}(\mathbf{r})=J_{yy}(\mathcal{R}^{z}(\frac{\pi}{2})\mathbf{r})$ which is
consistent with the symmetry property of the host material, except for minor
statistical errors caused in the QMC simulations. $J_{xy}(\mathbf{r})$ given
in (d) generally has the same symmetry property with that obtained from the
variational method. Although $J_{xz}(\mathbf{r})$ and $J_{yz}(\mathbf{r})$ in
(e) and (f) shows opposite signs in comparison to the counterparts in Fig. 3,
all of them follows the three-fold rotational symmetry.

Figure 9: (Color online) The red solid lines show the results of
$r^{3}J_{zz}(\mathbf{r})$ along the $x$-axis while $\alpha=0.1$ for (a)
$\mu=-1.0$, (b) $\mu=0.5$, (c) $\mu=1.5$ and (d) $\mu=3.0$. We use the
symmetric case $\epsilon_{d}-\mu=-U/2$ in our HFQMC simulations, so
$\varepsilon_{d}$ varies with respect to $\mu$. The blue dashed line in (d) is
the results for $10\times r^{3}J_{zz}(\mathbf{r})$ in a 2DEG with $\alpha=0$
for comparison. The parameters are $U=3.0$, $V_{k}=1.0$ and the temperature is
$K_{B}T=1/16$.

In Fig. 9, the red solid lines show the results of $r^{3}J_{zz}(\mathbf{r})$
along the $x$-axis while $\alpha=0.1$. The parameters are $U=3.0$, $V_{k}=1.0$
and the temperature is $K_{B}T=1/16$. When $\alpha=0.1$ and $\mu=-1.0$, as
shown in Fig. 9(a), only the lower band $\epsilon_{k-}$ involves in the
screening process. We can see the spatial decay rate of the spin-spin
correlation is about $r^{-3}$. As $\mu$ gradually increases, as in (b) and
(c), both bands $\epsilon_{k-}$ and $\epsilon_{k+}$ take part in the Kondo
screening, and the oscillation becomes more complicated. If $\mu=3.0$ as given
in Fig. 9(d), only the upper band $\epsilon_{k+}$ is responsible for the Kondo
screening, and the decay rate of the $J_{zz}(x,0)$ is still proportional to
$r^{-3}$, with different period of oscillation. For comparison, the results of
spin-spin correlation in a 2DEG for $\alpha=0$ is plotted as the blue dashed
line in (d). Note that if $\alpha=0$, the spin-spin correlation is much
smaller than the $\alpha=0.1$ case, so $r^{3}J_{zz}(x,0)$ is multiplied by
$10$ for clarity. Our results support the $1/r^{3}$ decay of the Kondo
screening cloud at long distances, which is consistent with previous
studies.Ishii (1978); Barzykin and Affleck (1998); Borda (2007) We can see
that the spatial decay rate obtained by the HFQMC results is more reliable
than those given by the variational method shown in Fig. 4. Our HFQMC results
show that the decay rate of the spin-spin correlation remains essentially
unchanged in the presence of cubic Rashba term. However, the oscillation
pattern and period are clearly affected by the cubic Rashba SOC.

##  IV Indirect coupling between TWO magnetic IMPURIties

Taking into account of the indirect coupling between two magnetic impurities,
one natural question shall be how the RKKY interaction is influenced by the
cubic Rashba term. For simplicity, we assume that one impurity is located at
the origin, and the other impurity is on the $x$-axis with a distance $R$, as
schematically plotted in Fig. 1(d). HFQMC returns the imaginary time Green’s
functions $G_{jj^{\prime}}^{ss^{\prime}}(\mathbf{R},\tau)$, where
$j,j^{\prime}=1,2$ mark the two magnetic atoms and
$s,s^{\prime}=\uparrow,\downarrow$ are the spin indices. The spin-spin
correlation between two magnetic impurities measures the RKKY interaction
mediated by the conduction electrons. The non-zero components of the spin-spin
correlation function along the $x$-axis areSun  _et al._ (2014)

$\displaystyle\langle S_{1}^{z}S_{2}^{z}\rangle$ $\displaystyle=$
$\displaystyle\langle S_{1}^{x}S_{2}^{x}\rangle$ $\displaystyle=$
$\displaystyle\langle(g_{11}^{\uparrow\uparrow}-g_{11}^{\downarrow\downarrow})\times(g_{22}^{\uparrow\uparrow}-g_{22}^{\downarrow\downarrow})-g_{12}^{\uparrow\uparrow}\cdot
g_{21}^{\uparrow\uparrow}$ $\displaystyle-$ $\displaystyle
g_{12}^{\downarrow\downarrow}\cdot
g_{21}^{\downarrow\downarrow}+g_{12}^{\uparrow\downarrow}\cdot
g_{21}^{\downarrow\uparrow}+g_{12}^{\downarrow\uparrow}\cdot
g_{21}^{\uparrow\downarrow}\rangle,$ $\displaystyle\langle
S_{1}^{y}S_{2}^{y}\rangle$ $\displaystyle=$
$\displaystyle\langle-(g_{11}^{\uparrow\downarrow}-g_{11}^{\downarrow\uparrow})\times(g_{22}^{\uparrow\downarrow}-g_{22}^{\downarrow\uparrow})+g_{12}^{\uparrow\downarrow}\cdot
g_{21}^{\uparrow\downarrow}$ $\displaystyle+$ $\displaystyle
g_{12}^{\downarrow\uparrow}\cdot
g_{21}^{\downarrow\uparrow}-g_{12}^{\uparrow\uparrow}\cdot
g_{21}^{\downarrow\downarrow}-g_{12}^{\downarrow\downarrow}\cdot
g_{21}^{\uparrow\uparrow}\rangle,$ $\displaystyle\langle
S_{1}^{x}S_{2}^{z}\rangle$ $\displaystyle=$ $\displaystyle-\langle
S_{1}^{z}S_{2}^{x}\rangle$ $\displaystyle=$
$\displaystyle\langle(g_{11}^{\uparrow\downarrow}+g_{11}^{\downarrow\uparrow})\times(g_{22}^{\uparrow\uparrow}-g_{22}^{\downarrow\downarrow})+g_{12}^{\uparrow\downarrow}\cdot
g_{21}^{\downarrow\downarrow}$ $\displaystyle-$ $\displaystyle
g_{12}^{\downarrow\uparrow}\cdot
g_{21}^{\uparrow\uparrow}-g_{12}^{\uparrow\uparrow}\cdot
g_{21}^{\uparrow\downarrow}+g_{12}^{\downarrow\downarrow}\cdot
g_{21}^{\downarrow\uparrow}\rangle.$

Figure 10: (Color online) The spin-spin correlation between two magnetic
impurities with respect to the distance between them $R$. (a)-(c) are listed
in the order of increasing $\mu$, and $\epsilon_{d}-\mu=-U/2$. $k_{c}$ is the
momentum truncation, and the parameters are chosen as $U=3.0$ and $V_{k}=1.0$,
$k_{B}T=1/8$. The cubic Rashba SOC term is $\alpha=0.1$, and VHS occurs at
energy $1.85$.

In Fig. 10 we show the spin-spin correlation between the two magnetic
impurities with respect to the distance $R$ between them, and $k_{c}$ is the
momentum truncation. The parameters are chosen as $\alpha=0.1$, $U=3.0$ and
$V_{k}=1.0$, $k_{B}T=1/8$. We consider the symmetric case, with
$\epsilon_{d}-\mu=-U/2$. Along the $x$-axis, we can see that $\langle
S_{1}^{z}S_{2}^{z}\rangle=\langle S_{1}^{x}S_{2}^{x}\rangle\neq\langle
S_{1}^{y}S_{2}^{y}\rangle$. This is due to the cubic Rashba term $\alpha$,
without which all the three components shall be exactly the same. When
$\alpha=0.1$, the VHS emerges at energy value $\mu=1.85$. Figs. 10(a)-(c) are
listed in the order of increasing $\mu$. $\mu=0.5$ given in Fig. 10(a)
corresponds to relatively low DOS, while $\mu=1.50$ and $\mu=1.85$ given in
Figs. 10(b) and (c) are close to the energies where VHS occurs. We can see
that the diagonal terms $\langle S_{1}^{z}S_{2}^{z}\rangle$ and $\langle
S_{1}^{y}S_{2}^{y}\rangle$ are suppressed when $\mu=0.5$. As $\mu$ approaches
the VHS point, the DOS increases, and so the diagonal terms $\langle
S_{1}^{z}S_{2}^{z}\rangle$ and $\langle S_{1}^{y}S_{2}^{y}\rangle$. For all
the cases, the diagonal terms are positive and dominant when the two
impurities are close, indicating that the two magnetic impurities are
ferromagnetically correlated, and the values oscillate and decay in space. The
only non-zero off-diagonal term $\langle S_{1}^{x}S_{2}^{z}\rangle=-\langle
S_{1}^{z}S_{2}^{x}\rangle$ also changes with respect to the values of $\mu$.
The off-diagonal terms correspond to the DM interaction,Imamura  _et al._
(2004); Mross and Johannesson (2009); Zhu  _et al._ (2011) and it is a
manifestation of the SOC in the host material. At distance $k_{c}R\sim 10$,
the off-diagonal correlation $\langle S_{1}^{x}S_{2}^{z}\rangle$ is of the
same order of magnitude as the diagonal terms.

##  V conclusions

In this paper, we apply the variational method and the HFQMC technique to
study the influence of the $k$-cubic Rashba SOC on the correlation effects of
magnetic impurities. The cubic Rashba SOC greatly alters the band structure
and induces a VHS to the host material. The $k$-linear Rashba SOC can also
cause the divergence of DOS, but the divergence occurs at the bottom of the
bands. However, the VHS induced by the cubic Rashba SOC occurs in a very wide
range of energy, and the single impurity local moment becomes largely tunable,
especially for strong $V_{k}$ and relatively weak $U$. Both the variational
method and the HFQMC simulations support the three- or six-fold rotational
symmetry of the various components of spatial spin-spin correlation. This
unique feature is a manifestation of the cubic Rashba SOC, and can possibly be
used in experiments to distinguish the cubic Rashba SOC from the normal
$k$-linear Rashba term. The HFQMC calculations show that the $1/r^{3}$ decay
rate of this spin-spin correlation is essentially unchanged by the cubic
Rashba SOC term $\alpha$. Moreover, the RKKY couplings between two magnetic
impurities displays very complicated form. Besides the normal diagonal
components, we still obtain the finite off-diagonal components, which
corresponds to the DM interaction between two magnetic impurities, and they
become the same order of magnitude as the diagonal terms at distance
$k_{c}R\approx 10$.

##  VI Acknowledgments

J.-H.S. acknowledges financial support from the Zhejiang Provincial Natural
Science Foundation of China (Grant No. LY19A040003) and K.C.Wong Magna Fund in
Ningbo University. D.-H.X. was supported by the NSFC (under Grant Nos.
12074108 and 12147102) and the Natural Science Foundation of Chongqing (Grant
No. CSTB2022NSCQ-MSX0568). L.C. was supported by the NSFC ( under Grant No.
12174101) and the Fundamental Research Funds for the Central Universities
(Grant No. 2022MS051)

## References

  * Dresselhaus (1955) G. Dresselhaus, Phys. Rev. 100, 580 (1955).
  * Rashba (1960) E. Rashba, Sov. Phys. Solid State 2, 1109 (1960).
  * J. Ohkawa and Uemura (1974) F. J. Ohkawa and Y. Uemura, J. Phys. Soc. Jpn. 37, 1325 (1974).
  * Vas’ ko (1979) F. Vas’ ko, P. Zh. Eksp. Teor. Fiz. 30, 574 (1979).
  * Bychkov and Rashba (1984) Y. A. Bychkov and É. I. Rashba, P. Zh. Eksp. Teor. Fiz. 39, 66 (1984).
  * Bihlmayer  _et al._ (2015) G. Bihlmayer, O. Rader, and R. Winkler, New J. Phys. 17, 050202 (2015).
  * Manchon  _et al._ (2015) A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov, and R. A. Duine, Nat. mater. 14, 871 (2015).
  * Bihlmayer  _et al._ (2022) G. Bihlmayer, P. Noël, D. V. Vyalikh, E. V. Chulkov, and A. Manchon, Nat. Rev. Phys 4, 642 (2022).
  * Gerchikov and Subashiev (1992) L. Gerchikov and A. Subashiev, Sov. phys. Semiconduct. 26, 73 (1992).
  * Winkler (2003) R. Winkler,  _Spin-orbit coupling effects in two-dimensional electron and hole systems_ , Vol. 191 (Springer, 2003).
  * Winkler (2000) R. Winkler, Phys. Rev. B 62, 4245 (2000).
  * Zhao  _et al._ (2020) H. J. Zhao, H. Nakamura, R. Arras, C. Paillard, P. Chen, J. Gosteau, X. Li, Y. Yang, and L. Bellaiche, Phys. Rev. Lett. 125, 216405 (2020).
  * Schliemann and Loss (2005) J. Schliemann and D. Loss, Phys. Rev. B 71, 085308 (2005).
  * Bleibaum and Wachsmuth (2006) O. Bleibaum and S. Wachsmuth, Phys. Rev. B 74, 195330 (2006).
  * Ma and Liu (2006) T. Ma and Q. Liu, Phys. Rev. B 73, 245315 (2006).
  * Minkov  _et al._ (2005) G. M. Minkov, A. A. Sherstobitov, A. V. Germanenko, O. E. Rut, V. A. Larionova, and B. N. Zvonkov, Phys. Rev. B 71, 165312 (2005).
  * Winkler  _et al._ (2002) R. Winkler, H. Noh, E. Tutuc, and M. Shayegan, Phys. Rev. B 65, 155303 (2002).
  * Nakamura  _et al._ (2012) H. Nakamura, T. Koga, and T. Kimura, Phys. Rev. Lett. 108, 206601 (2012).
  * Usachov  _et al._ (2020a) D. Y. Usachov, I. A. Nechaev, G. Poelchen, M. Güttler, E. E. Krasovskii, S. Schulz, A. Generalov, K. Kliemt, A. Kraiker, C. Krellner, K. Kummer, S. Danzenbächer, C. Laubschat, A. P. Weber, J. Sánchez-Barriga, E. V. Chulkov, A. F. Santander-Syro, T. Imai, K. Miyamoto, T. Okuda, and D. V. Vyalikh, Phys. Rev. Lett. 124, 237202 (2020a).
  * Usachov  _et al._ (2020b) D. Y. Usachov, M. Güttler, S. Schulz, G. Poelchen, S. Seiro, K. Kliemt, K. Kummer, C. Krellner, C. Laubschat, E. V. Chulkov, and D. V. Vyalikh, Phys. Rev. B 101, 245140 (2020b).
  * Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Kondo (1964) J. Kondo, Prog. Theor. Phys. 32, 37 (1964).
  * Wilson (1975) K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
  * Ishii (1978) H. Ishii, J. Low Temp. Phys. 32, 457 (1978).
  * Barzykin and Affleck (1998) V. Barzykin and I. Affleck, Phys. Rev. B 57, 432 (1998).
  * Borda (2007) L. Borda, Phys. Rev. B 75, 041307 (2007).
  * Moca  _et al._ (2021) C. P. Moca, I. Weymann, M. A. Werner, and G. Zaránd, Phys. Rev. Lett. 127, 186804 (2021).
  * V Borzenets  _et al._ (2020) I. V Borzenets, J. Shim, J. C. Chen, A. Ludwig, A. D. Wieck, S. Tarucha, H.-S. Sim, and M. Yamamoto, Nature 579, 210 (2020).
  * Malecki (2007) J. Malecki, J. Stat. Phys. 129, 741 (2007).
  * Žitko and Bonča (2011) R. Žitko and J. Bonča, Phys. Rev. B 84, 193411 (2011).
  * Isaev  _et al._ (2012) L. Isaev, D. F. Agterberg, and I. Vekhter, Phys. Rev. B 85, 081107 (2012).
  * Zarea  _et al._ (2012) M. Zarea, S. E. Ulloa, and N. Sandler, Phys. Rev. Lett. 108, 046601 (2012).
  * Chen  _et al._ (2017) L. Chen, B. Hu, and R.-S. Han, J. Phys. Condens. Matter. 30, 025601 (2017).
  * Wong  _et al._ (2016) A. Wong, S. E. Ulloa, N. Sandler, and K. Ingersent, Phys. Rev. B 93, 075148 (2016).
  * Chen  _et al._ (2016) L. Chen, J. Sun, H.-K. Tang, and H.-Q. Lin, J. Phys. Condens. Matter. 28, 396005 (2016).
  * Li  _et al._ (2018) L. Li, J.-H. Sun, Z.-H. Wang, D.-H. Xu, H.-G. Luo, and W.-Q. Chen, Phys. Rev. B 98, 075110 (2018).
  * Feng  _et al._ (2010) X.-Y. Feng, W.-Q. Chen, J.-H. Gao, Q.-H. Wang, and F.-C. Zhang, Phys. Rev. B 81, 235411 (2010).
  * Feng and Zhang (2011) X.-Y. Feng and F.-C. Zhang, J. Phys. Condens. Matter 23, 105602 (2011).
  * Chirla  _et al._ (2013) R. Chirla, C. P. Moca, and I. Weymann, Phys. Rev. B 87, 245133 (2013).
  * Kasuya (1956) T. Kasuya, Prog. Theor. Phys. 16, 45 (1956).
  * Ruderman and Kittel (1954) M. A. Ruderman and C. Kittel, Phys. Rev. 96, 99 (1954).
  * Yosida (1957) K. Yosida, Phys. Rev. 106, 893 (1957).
  * Imamura  _et al._ (2004) H. Imamura, P. Bruno, and Y. Utsumi, Phys. Rev. B 69, 121303 (2004).
  * Mross and Johannesson (2009) D. F. Mross and H. Johannesson, Phys. Rev. B 80, 155302 (2009).
  * Zhu  _et al._ (2011) J.-J. Zhu, D.-X. Yao, S.-C. Zhang, and K. Chang, Phys. Rev. Lett. 106, 097201 (2011).
  * Hirsch and Fye (1986) J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986).
  * Gunnarsson and Schönhammer (1983) O. Gunnarsson and K. Schönhammer, Phys. Rev. Lett. 50, 604 (1983).
  * Varma and Yafet (1976) C. M. Varma and Y. Yafet, Phys. Rev. B 13, 2950 (1976).
  * Sun  _et al._ (2015) J.-H. Sun, D.-H. Xu, F.-C. Zhang, and Y. Zhou, Phys. Rev. B 92, 195124 (2015).
  * Sun  _et al._ (2018) J.-H. Sun, L.-J. Wang, X.-T. Hu, L. Li, and D.-H. Xu, Phys. Rev. B 97, 035130 (2018).
  * Ma  _et al._ (2018) D. Ma, H. Chen, H. Liu, and X. C. Xie, Phys. Rev. B 97, 045148 (2018).
  * Wang  _et al._ (2019) L.-J. Wang, X.-T. Hu, L. Li, D.-H. Xu, J.-H. Sun, and W.-Q. Chen, Phys. Rev. B 99, 235108 (2019).
  * Yang  _et al._ (2021) X.-R. Yang, Z.-Z. Huang, W.-S. Wang, and J.-H. Sun, Chin. Phys. B 30, 067103 (2021).
  * Simonin and Allub (1995) J. Simonin and R. Allub, Phys. Rev. Lett. 74, 466 (1995).
  * Simon and Varma (1999) M. E. Simon and C. M. Varma, Phys. Rev. B 60, 9744 (1999).
  * Rozhkov and Arovas (2000) A. V. Rozhkov and D. P. Arovas, Phys. Rev. B 62, 6687 (2000).
  * Huang  _et al._ (2022) Z.-Z. Huang, X.-T. Peng, W.-S. Wang, and J.-H. Sun, Chinese Physics B 31, 107101 (2022).
  * Haldane and Anderson (1976) F. D. M. Haldane and P. W. Anderson, Phys. Rev. B 13, 2553 (1976).
  * Uchoa  _et al._ (2008) B. Uchoa, V. N. Kotov, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 101, 026805 (2008).
  * Fye and Hirsch (1988) R. M. Fye and J. E. Hirsch, Phys. Rev. B 38, 433 (1988).
  * Fye  _et al._ (1987) R. M. Fye, J. E. Hirsch, and D. J. Scalapino, Phys. Rev. B 35, 4901 (1987).
  * Hirsch and Lin (1987) J. E. Hirsch and H. Q. Lin, Phys. Rev. B 35, 4943 (1987).
  * Bulut  _et al._ (2007) N. Bulut, K. Tanikawa, S. Takahashi, and S. Maekawa, Phys. Rev. B 76, 045220 (2007).
  * Hu  _et al._ (2011) F. M. Hu, T. Ma, H.-Q. Lin, and J. E. Gubernatis, Phys. Rev. B 84, 075414 (2011).
  * Sun  _et al._ (2013a) J. Sun, F. Hu, H. Tang, W. Guo, and H. Lin, J. Appl. Phys 113, 17B515 (2013a).
  * Sun  _et al._ (2013b) J. Sun, F. Hu, H. Tang, and H. Lin, Int. J. Mod. Phys. B 27, 1362039 (2013b).
  * Sun and Tang (2018) J.-H. Sun and H.-K. Tang, Chin. Phys. B 27, 077502 (2018).
  * Sun  _et al._ (2014) J. Sun, L. Chen, and H.-Q. Lin, Phys. Rev. B 89, 115101 (2014).
  * Hu  _et al._ (2013) F. M. Hu, T. O. Wehling, J. E. Gubernatis, T. Frauenheim, and R. M. Nieminen, Phys. Rev. B 88, 045106 (2013).
  * Moriya  _et al._ (2014) R. Moriya, K. Sawano, Y. Hoshi, S. Masubuchi, Y. Shiraki, A. Wild, C. Neumann, G. Abstreiter, D. Bougeard, T. Koga, and T. Machida, Phys. Rev. Lett. 113, 086601 (2014).
  * Zarea and Ulloa (2006) M. Zarea and S. E. Ulloa, Phys. Rev. B 73, 165306 (2006).
  * Liu  _et al._ (2018) H. Liu, E. Marcellina, A. R. Hamilton, and D. Culcer, Phys. Rev. Lett. 121, 087701 (2018).
  * Xiao  _et al._ (2019) C. Xiao, Y. Liu, Z. Yuan, S. A. Yang, and Q. Niu, Phys. Rev. B 100, 085425 (2019).
  * Liu and Lei (2005) S. Y. Liu and X. L. Lei, Phys. Rev. B 72, 155314 (2005).
  * Karwacki  _et al._ (2018) L. Karwacki, A. Dyrdał, J. Berakdar, and J. Barnaś, Phys. Rev. B 97, 235302 (2018).
  * Gubernatis  _et al._ (1987) J. E. Gubernatis, J. E. Hirsch, and D. J. Scalapino, Phys. Rev. B 35, 8478 (1987).