Molecular molds for regularizing Kondo states at atom/metal interfaces

Adsorption of magnetic transition metal atoms on a metal surface leads to the formation of Kondo states at the atom/metal interfaces. However, the significant influence of surrounding environment presents challenges for potential applications. In this work, we realize a novel strategy to regularize the Kondo states by moving a CoPc molecular mold on an Au(111) surface to capture the dispersed Co adatoms. The symmetric and ordered structures of the atom-mold complexes, as well as the strong dπ–π bonding between the Co adatoms and conjugated isoindole units, result in highly robust and uniform Kondo states at the Co/Au(111) interfaces. Even more remarkably, the CoPc further enables a fine tuning of Kondo states through the molecular-mold-mediated superexchange interactions between Co adatoms separated by more than 12 Å. Being highly precise, efficient and reproducible, the proposed molecular mold strategy may open a new horizon for the construction and control of nano-sized quantum devices. The surrounding environment’s influence on the Kondo states formed around an adsorbed transition metal atom are challenging potential applications. Here, the authors demonstrate, via STM and DFT+HEOM calculations, that a metal phthalocyanine molecule on an Au(111) surface enables a fine tuning of the Kondo characteristics.

M agnetic atoms are important building blocks of nanodevices for applications in quantum information and quantum computation. For realizing on-demand design and fabrication of nano-devices, it is crucial to achieve precise manipulation and tuning of local quantum states in these magnetic nanostructures [1][2][3][4][5][6][7] . This thus requires the characteristic quantum states stay robust upon variation of surrounding environment, and respond regularly and consistently under external manipulation.
Adsorption of a transition metal (TM) atom on a metal surface could lead to the formation of Kondo states at the atom/metal interface under sufficiently low temperatures, which has been affirmed by the substantially enhanced electric conductance across the interface 8,9 . Physically, the Kondo states originate from the screening of the local spin of d electrons on the TM atom by the spins of itinerant electrons in the metal substrate. Recently, the precise control and tuning of Kondo states have become a focus of experimental efforts [10][11][12][13][14][15][16][17] . However, it is generally observed that the Kondo features vary greatly with the local chemical environment surrounding the TM atom, such as the coordinating ligands, the chemical dopants, the metal substrate, etc. This makes it difficult to produce Kondo states at atom/metal interfaces with regular sizes and shapes.
It is well-known that Kondo states emerge at a Co atom adsorbed on the fcc or hcp domain of Au(111) surface, whereas a Co atom adsorbed on the domain wall does not exhibit Kondo signature 8 . This indicates that the interaction between the Co atom and the gold substrate depends sensitively on the atomistic structure of the surface. Moreover, the Kondo states centered at a Co adatom can be significantly influenced by another Co adatom in its vicinity, because of the direct 9 or indirect 14,18 interatomic interactions. For instance, the substrate-mediated Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction between two Co adatoms can change from antiferromagnetic to ferromagnetic as the interatomic distance varies, leading to enhanced or suppressed Kondo screening 4,19 .
In many experiments, a TM adatom is bound to organic ligands to form an organometallic compound [20][21][22][23] , and its local electronic configuration is significantly influenced by the surrounding. However, the related Kondo states do not necessarily have regular shapes, because the ligand field and the local spin distribution are highly susceptible to variations in the surrounding environment. For instance, for an iron phthalocyanine (FePc) molecule adsorbed on the Au(111) surface, the dI/dV spectra measured at the on-top adsorption site exhibit Kondo signatures that are distinctly different from that at the bridge site 20 . In addition, the presence of a nearby FePc molecule can lead to conspicuous broadening or splitting of the Kondo conductance peak, which is again caused by the substrate-mediated RKKY interaction between the two local spin moments 21 .
The above experimental observations indicate that the local chemical environment has a profound influence on the Kondo states formed at atom/metal interfaces. However, first-principlesbased theoretical study on the Kondo states has remained rather scarce. This poses a serious challenge for potential practical applications-is there a way to regularize the Kondo states, so that they could exhibit uniform features that are insensitive to the variation of environment? To solve this problem, in this paper we propose a novel strategy, with which we successfully realize the regularization of Kondo states centered at the Co adatoms on an Au(111) surface.
Our strategy is motivated by the experience that a regular shape can be attained by using a rigid mold. Following this idea, we propose to use a planar molecule with conjugated rings to capture the dispersed TM adatoms on the metal surface. Each conjugated ring also serves as a mold, which regularizes the local spin distribution as well as the resulting Kondo states centered at the captured TM adatom. Consequently, the influence of environment beyond the range of the planar molecule is substantially weakened.
In this work, we choose to use a metal phthalocyanine (MPc) molecule, which possesses four isoindole units. Each isoindole unit can capture one Co adatom through the strong bonding interaction between the conjugated π-orbitals on the isoindole and the d orbitals on the Co adatom. More interestingly, if there is more than one isoindole unit hosting a Co adatom, long-range superexchange interactions may emerge between these spatially separated Co adatoms, which allows for fine-tuning of Kondo states in a controllable manner.

Results
Construction of K n (CoPc) on Au(111) surface. Figure 1 illustrates the experimental procedure using an MPc molecule to capture the Co atoms dispersed on the Au(111) surface. As shown in Fig. 1a, a CoPc molecule exhibits a clear four-lobe pattern in the scanning tunneling microscopy (STM) image, with each lobe representing an isoindole unit. Pushed by the atomically sharp STM tip, the CoPc molecule can be quite freely moved on the surface. When the CoPc is pushed towards and finally in contact with a Co adatom, one of its four lobes captures the Co adatom and thus forms an atom-molecule complex (referred to as K 1 hereafter). In K 1 the Co atom is located underneath the molecular plane, and hence it interacts directly with both the isoindole unit and the gold surface. The K 1 complex can also be moved around on the surface, while its structure remains stable (see Fig. 1b). This also indicates that the Co atom is bound more tightly to the isoindole unit than to the surface.
Manipulated by the STM tip (see Supplementary Movie 1), a CoPc can capture up to four Co adatoms. Figure 1c demonstrates the sequential formation of the K 2 0 , K 3 , and K 4 complexes from K 1 . Here, K n denotes the complex with n Co adatoms captured by a CoPc. K 2 and K 2 0 are isomers, in which the two Co adatoms locate at the ortho and para positions, respectively. During the process, the CoPc molecule plays the role of Pac-Man, while the Co adatoms are "devoured" like Pac-Dots. The atomistic structures of these complexes are sketched in Fig. 1d, and the presence of the Co adatoms are verified by the topographic images depicted in Fig. 1e. Figure 1f depicts the Kondo maps of these complexes (see Supplementary Figs. 1 and 2). Clearly, a Kondo cloud emerges at the outer region of each occupied isoindole unit, while the cobalt ion at the center is Kondo inactive. It is remarkable that all the Kondo clouds have a similar crescent shape, except in K 2 0 where the two clouds appear more isotropic. This thus highlights the strong regularization effect of the CoPc molecular mold.
Experimental measurement of dI/dV spectra and T K . We now examine more closely the Kondo characteristics by investigating the measured differential conductance (dI/dV versus V) spectra. As shown in Fig. 2a, the measured spectra of a bare Co adatom exhibit a Fano line shape depending sensitively on the adsorption domain. The Fano line shape originates from the quantum interference between two electron conduction channels 24 : the through-space and through-atom channels. In contrast, the measured spectra of a K 1 complex exhibit a single-peak structure independent of the adsorption domain. This is because the CoPc molecule blocks the through-space tunneling channel and thus quenches the Fano interference 24,25 . The Kondo temperature T K is extracted from the width of the measured Kondo conductance peak (Γ), followed by deconvolution of thermal and instrumentinduced broadening (see Methods).
From Fig. 2a, an isolated K 1 has an almost constant T K of about 122 K. By manipulating the STM tip, a K 1 can be moved towards a bare Co adatom. Throughout this process, the spectral line shape of the K 1 remains largely unchanged, while the T K exhibits a weak oscillation around the constant value (see Fig. 2b and Supplementary Fig. 3). Such an oscillation is supposed to originate from the substrate-mediated RKKY interaction between the Co adatom in the K 1 and the nearby bare Co adatom 4,19 , and it is greatly suppressed with the nearby bare Co adatom captured by a CoPc mold. This affirms that the strong interaction between a Co adatom and the associated isoindole unit on the CoPc gives rise to highly regular and robust local spin distribution, so that the Kondo states are hardly affected by the variation of environment beyond the phthalocyanine ring.
For a K n complex hosting two or more Co adatoms, the dI/dV spectra measured at each of the Co adatoms are almost identical (see Supplementary Fig. 2). While the spectral line shape of a K n (n ≥ 2) much resembles that of a K 1 , the Kondo peak of the former is somewhat wider, suggesting a higher T K (see Fig. 2c). Figure 2d gives an overview of T K for all members of the K n family. Two main features are easily recognized: (1) T K increases linearly with n for n = 1, 2, 3, 4. (2) T K of K 2 0 is distinctly higher than that of its isomer K 2 , and lies outside the above linear relationship. These features indicate that the Kondo state centered at a captured Co adatom is affected subtly but regularly by the other Co adatoms confined within the same phthalocyanine ring. This lays the foundation for fine-tuning the Kondo states formed at atom/metal interfaces through precise arrangement of the TM adatoms.
Theoretical calculations of dI/dV spectra and T K . The above experimental manipulations and measurements affirm the efficacy of CoPc molecular mold for regularizing the Kondo states at  Co/Au(111) interfaces. For practical purposes, it is crucial to elucidate the regularization mechanisms, e.g., what is the origin of the regular spectral line shapes, and what leads to the intriguing T K versus n relationship. To provide theoretical insights into these questions, we carry out first-principles-based calculations by a combination of density functional theory (DFT) 26,27 and hierarchical equations of motion (HEOM) 28-30 approach. In Fig. 2d, the T K of K n extracted from the calculated dI/dV spectra are compared against the experimental data (to be elaborated later). Figure 3a depicts the optimized geometries and spin density distribution of a bare Co atom and a K n complex (n ¼ 1; 2 0 ; 3) adsorbed on the Au(111) surface. Apparently, the planar CoPc mold is always parallel to the surface, and the captured Co adatom is right beneath the six-member ring of an isoindole unit. Such a geometry remains unchanged upon moving the K n to different sites on the surface (Supplementary Figs. 4-6). The spinunpaired electrons reside predominantly on the Co adatoms with only a small fraction on the central cobalt ion. Thus, each Co adatom becomes the center of a Kondo cloud, whereas the cobalt ion on the CoPc mold is Kondo inactive.
From Fig. 3a it is also noted that the local spin density of a bare Co adatom exhibits an appreciable change at different adsorption sites. In contrast, the local spin of a captured Co adatom in a K n is hardly affected by the surrounding environment (see also Supplementary Figs. 7 and 8). These results are consistent with the experimental findings as displayed in Fig. 2a, b. Therefore, the CoPc mold indeed regularizes and preserves the local spins of the Co adatoms, and thus facilitates the formation of uniform Kondo states at the K n /Au(111) interfaces.
By analyzing the projected density of states (PDOS) of captured Co adatoms, we obtain a schematic diagram showing the relative energy, broadening and electron occupancy of each d orbital; see the upper right inset of Fig. 3b. The local electron configuration of a captured Co adatom is determined to be d 1:6 xy d 1:6 x 2 Ày 2 d 1:8 z 2 d 1:3 xz d 1:2 yz (see Supplementary Table 1), and the local spin moment originates mainly from the d π (d xz and d yz ) orbitals. Thus, a captured Co adatom is considered to be in a local S ¼ 1 2 state, in clear contrast to the local S = 1 state of a bare Co adatom 31 .
Because of their common symmetry, the d π orbitals on a Co adatom and the conjugated π-orbitals on an isoindole unit of the Kondo characteristics of the K n complexes. a Measured dI/dV spectra of a bare Co adatom and a K 1 adsorbed on different domains of the Au(111) surface. b Variation of T K of a K 1 (or a bare Co adatom) versus its distance to another K 1 (or bare Co adatom) on the Au(111) surface. c Measured dI/dV spectra of the K n complexes (n ¼ 1; 2; 2 0 ; 3; 4). The red lines are fits of the measured data to the Fano function (see Methods for details), and they have been shifted vertically for clarity. The measured data are acquired by placing the STM tip at positions marked by the black crosses in Fig. 1a, b under the background temperature T = 5 K. d T K of the K n complexes extracted from experimentally measured and theoretically calculated dI/dV spectra. The red bars mark the standard deviations of experimental data. The dashed line is a linear fit of the measured T K of K n (n = 1, 2, 3, 4). The circles (triangles) represent the values of T K determined by the heights (widths) of the calculated Kondo conductance peaks.  CoPc mold interact strongly with each other. The resulting d π -π bonds further consolidate the atom-mold complex, and give rise to the dumbbell-shaped spin density distribution (see Fig. 3a). The physical origin of the regularization effect is thus clear. Confined by the CoPc mold, the local electronic structure of a Co adatom is significantly reshaped and regularized by the strong d π -π bonds, leading to a local spin state distinctly different from that of a bare Co adatom.
The strength of Kondo correlation depends critically on how strongly the local spin on a Co adatom is screened by the electronic spins in the surrounding environment. For the K n /Au (111) composites, the Kondo screening is realized via the hybridization of the Co d π orbitals with the s orbitals of the nearby gold atoms. It is such s-d hybridization that leads to the broadening and splitting of the Co d π orbitals as shown in Fig. 3b and Supplementary Fig. 9. The hybridization strength is characterized by Δ s , which is taken as the half-width of the split PDOS peaks that constitute the main contribution to the local spin moment.
With more Co adatoms captured by a same CoPc, the split peaks of PDOS are slightly more broadened. This signifies a perturbation by the other Co adatoms confined within the same CoPc mold, possibly through the substrate-mediated RKKY interactions. The magnitude of such a perturbation is thus roughly proportional to the number of Co adatoms inside the atom-mold complex. Intriguingly, the split peaks of a K 2 0 is somewhat broader than a K 2 , implying the symmetry of the complex has a subtle influence on the strength of RKKY interaction. Accordingly, Δ s increases linearly with n except for n ¼ 2 0 (see the lower right inset of Fig. 3b). The apparent resemblance between the Δ s versus n and T K versus n relationships (cf. Fig. 2d) affirms the key importance of Δ s to the characteristic features of the Kondo states 32 .
To quantitatively describe the Kondo effect, we employ quantum impurity models which explicitly include the electron-electron interactions. We first assume the Kondo clouds centered at different Co adatoms are independent of each other, and thus a tip/K n /Au(111) junction can be represented by a single-orbital Anderson impurity model (AIM) 33 . The calculated dI/dV spectra agree remarkably with the experimental measurements (see Supplementary Fig. 10). The T K of K n are determined based on the height or the width of the Kondo conductance peaks (see Methods and also Supplementary Figs. 11 and 12), and the theoretical values accurately and consistently reproduce the experimental data (see Fig. 2d).
With the single-orbital AIM, the T K of a K 2 0 is predicted to be 157 K (labeled by K Ã 2 0 in Fig. 2d), only slightly higher than that of a K 2 , yet considerably lower than the experimental value of 191 K (Supplementary Fig. 13 and Supplementary Table 3). This indicates that the use of a single-orbital AIM is inadequate, because the two Kondo clouds at the K 2 0 /Au(111) interface are not independent to each other.
Influence of long-range superexchange interaction on T K . To gain deeper insights into the unusually high T K of a K 2 0 , we perform additional experiments by replacing CoPc with CuPc and H 2 Pc molecules (see Supplementary Note 12 for details). Similar atom-mold complexes are constructed by manipulating the STM tip. Unlike the case of CoPc, the measured T K of K 2 0 (CuPc) and K 2 0 (H 2 Pc) are only slightly higher than their isomers K 2 (CuPc) and K 2 (H 2 Pc) by 4-5 K (see Supplementary Figs. 14  and 15). These findings highlight the unique role of the central cobalt ion in the CoPc, which may serve as a hub to assist the two Co adatoms in the para positions to interact with each other.
From the calculated electronic structures of K 2 0 (CoPc)/Au(111) and K 4 (CoPc)/Au(111) composites (see Supplementary Figs. 16  and 17), a delocalized Kohn-Sham orbital, which links the d π orbitals on the two distantly separated Co adatoms via the π-orbitals on the isoindole units and a d π orbital on the central cobalt ion, is recognized (see Fig. 4a, b). This thus provides a channel for the superexchange (SE) interaction 34,35 between the two separated local spins. However, such a channel breaks down in a K 2 0 (H 2 Pc) or a K 2 0 (CuPc) on the Au(111) surface, because the bridging orbital at the central hub is absent or energetically incompatible (see Fig. 4c, d). Consequently, the long-range SE interaction is strong in a K 2 0 (CoPc), but either absent or rather weak in a K 2 0 (H 2 Pc) or K 2 0 (CuPc).
To assess the influence of SE on the Kondo states in a K 2 0 , a two-orbital AIM is employed, which explicitly includes the molecular-mold-mediated SE interaction between the two Co adatoms. Accurate computation of the strength of SE interaction (Δ SE ) by ab initio quantum chemistry methods is highly desirable yet rather challenging. Instead, Δ SE is determined by exploring the variation of T K . To retrieve the experimental value of T K = 191 K (see the inset of Fig. 4), we have Δ SE = 123 meV for a K 2 0 (CoPc), which amounts to an effective ferromagnetic coupling 36 of J eff ≈ 4 meV between the two local spins separated by as far as 12.5 Å (see Supplementary Notes 6 and 15 for details).
As shown in Fig. 4, the T K of a K 2 0 (CoPc) exhibits a nonmonotonic dependence on Δ SE . Such a phenomenon can be understood as follows. As Δ SE increases from zero, the local spin on a Co adatom starts to feel the ferromagnetic interaction from the other Co adatom, and thus it is screened less strongly by the surrounding environment, which leads to the attenuated Kondo correlation. However, when Δ SE reaches a certain value (see Supplementary Fig. 18), the whole K 2 0 (CoPc) complex favors an S = 1 state, resulting in a substantially enlarged spin moment for the total complex. Consequently, the Kondo temperature rises again with further increasing of Δ SE . Among all the other K n , only K 4 (CoPc) possesses a delocalized SE channel. This indicates that the structural symmetry of the whole complex is vital for the existence of SE. Nevertheless, the SE in a K 4 is much weaker than in a K 2 0 (see Supplementary  Fig. 19), because the delocalized orbital involves also the two disconnected Co adatoms at the ortho positions. It is estimated that the SE could lead to a minor change in T K by~4 K for a K 4 .
Besides the strength of Kondo correlation, the SE also affects the spatial distribution of Kondo states. As noted earlier, the Kondo clouds at the K 2 0 /Au(111) interface appear to be more isotropic than in the other composites, which is possibly due to the strong long-range SE.

Discussion
To summarize, unlike all the previous studies in which the Kondo states formed around an adsorbed TM atom always vary sensitively to the surrounding environment, we demonstrate that, the proposed strategy of using a symmetric and conjugated molecular mold gives rise to highly robust and regularized Kondo states. Adjusting the chemical composition of the molecular mold further allows for a fine-tuning of the Kondo characteristics. The resulting atom-mold complexes can move freely on the gold surface and exhibit remarkably uniform Kondo features. Therefore, they may serve as standard building blocks for the design and fabrication of novel quantum devices.

Methods
Scanning tunneling microscopy experiments. The experiments are carried out with a low temperature STM (Omicron) using Au(111)/mica substrates in vacuum under a base pressure of 3 × 10 −11 Torr. The Au(111) surface is cleaned by repeat cycles of Ar ion sputtering at 800 V for 15 min and annealed at 600 K. Submonolayers of Co atoms and CoPc molecules are co-evaporate by e-beam evaporation and by sublimation, respectively, on an Au(111) substrate. The substrates are placed in situ on the cryostat of the microscope with a temperature of about 6.7 K. The samples are then investigated in situ at 5 K. The chemically etched tungsten tip is carefully cleaned by circular Argon ion sputtering at 800-1000 V for about 5-10 min and by field emission lasting for about 30 s under a sample bias of around −100 V with emission current of 2 μA in each cycle 37,38 .
In the consecutive formation of K n complexes, the manipulation conditions are 50 mV and 10 nA, while the imaging conditions are 1 V and 0.2 nA for Fig. 1a-c. For the measurement of topographic images and simultaneously acquired Kondo maps in Fig. 1e-f, the imaging conditions are −15 mV and 2 nA with a sinusoidal modulation of 2 mV and 730 Hz. The dI/dV spectra in Fig. 2 are recorded with a lock-in method. The acquisition conditions are −150 mV and 2 nA with a sinusoidal modulation of 2 mV and 730 Hz. The Kondo maps and dI/dV spectra of the K n complexes are acquired at different domains of the Au(111) reconstruction surface, including hcp and fcc domains and the domain walls.
The Kondo line shapes in the dI/dV spectra are fit to a single Fano function 39 as follows: where ε = (eV − ε K )/Γ, ε K is the center of the Kondo resonance, Γ is the half-width at half-maximum (HWHM) of the Kondo resonance peak, and q is the asymmetry parameter. The Kondo temperature T K is extracted from Γ by deconvolution of the thermal and instrument-induced broadening as follows [40][41][42][43][44] : Here, k B is the Boltzmann's constant, T = 5 K is the environmental temperature, the coefficient λ assumes an empirical value of 2.7 as adopted in previous experiments 42 , and V m = 1 mV is the amplitude of the modulation voltage (corresponding to a peak-to-peak value of 2 mV).
It is worth pointing out that in principle the thermal and instrument-induced broadening of any spectroscopic feature should be examined by a convolution of the expected signal with a proper broadening kernel [40][41][42][43] . In practice it is often more convenient to use the approximate formula of Eq. (2).
Density functional theory method. First-principles calculations are performed by using the spin-polarized density functional theory (DFT) method implemented in the Vienna ab initio simulation package (VASP) 45 . The generalized gradient approximation developed by Perdew, Burke and Ernzerholf (PBE) 46 is chosen for the exchange-correlation functional. The projected augmented wave method is adopted with the energy cutoff of 400 eV. The zero damping DFT-D3 method of Grimme 47 is used to improve the description of van der Waals interactions.
During the structural optimizations, a slab model including three layers of Au atoms is used to represent the Au(111) surface, with each layer containing 56 Au atoms, A vacuum space of 16 Å is used to eliminate the interactions between different slabs. All the atoms except those in the bottom two Au layers are fully relaxed until the residual force on every atom was <0.02 eV Å −1 . The convergence criterion is set to 1 × 10 −5 eV for the total energy. Because of the large size of the supercell, the Γ-point approximation is adopted. Additional numerical tests have been carried out, which have verified that a larger slab model consisting of four atomic layers does not change the main results qualitatively.
Anderson impurity models. The AIMs 33 which explicitly include the electron-electron Coulomb interaction and the Kondo effects are adopted to represent the tip/K n /Au(111) junctions. The total AIM Hamiltonian is comprised of three parts: where H imp and H env are the Hamiltonian of the magnetic impurity and that of its environment, respectively. In particular, the environment is modeled by reservoirs of noninteracting electrons, i.e., where H tip and H sub represent the STM tip (α = t) and the gold substrate (α = s), respectively.ĉ y αkσ ðĉ αkσ Þ is the creation (annihilation) operator for a spin-σ electron on the kth orbital in the αth reservoir. The impurity-environment coupling has the form of H coup ¼ P ασki t αkiĉ y αkσdiσ þ H:c:, where {t αki } are the coupling strengths between the ith impurity orbital and the kth reservoir orbital. The influence of electron reservoirs on the impurity is accounted for through the hybridization functions defined by Λ ijα ðωÞ π P k t αki t Ã αkj δðω À ϵ αk Þ. In a single-orbital AIM, the spin-polarized d π orbitals on each captured Co adatom are represented by a single impurity orbital, which is described by Here,n σ ¼d y σdσ , whered y σ (d σ ) creates (annihilates) an electron of spin-σ (σ = ↑ or ↓) on the impurity orbital of energy ϵ d ; and U is the intra-orbital electronelectron interaction energy. The reservoir hybridization functions assume a Lorentzian form of where Δ α is the effective coupling strength between the impurity orbital and the αth reservoir (α = s or t), and Ω α (W α ) is the band center (width) of the αth reservoir.
The K 2 0 /Au(111) composite is also represented by a two-orbital AIM, in which the d π orbitals on each of the two Co adatoms are modeled by an impurity orbital, and the total impurity is described by The coupling between the ith orbital and αth reservoir is characterized by the hybridization function Λ iiα ðωÞ ¼ The values of the involving energetic parameters, including ϵ d , U, Δ α , Ω α , and W α , are extracted from the DFT calculation results (see Supplementary Table 2) 37,48,49 .
Hierarchical equations of motion method. The HEOM method for fermionic baths 29,30 implemented in the HEOM-QUICK program 50 is employed to solve the stationary states of the AIMs, and to compute the key physical observables to compare directly with the experimental measurements. The detailed form of the HEOM can be found in ref. 29 . Its basic variables are the system reduced density matrix and a hierarchical set of auxiliary density operators (ADOs).
To ensure the numerical accuracy for ultralow temperatures, a recently developed Fano spectrum decomposition scheme 51 is adopted to accurately unravel the reservoir correlation functions. Limited by the computer resources at our disposal, the lowest temperature accessed by the HEOM method is 30 K for the single-orbital AIM, with the truncation tier set to L = 4. The tunneling current across a tip/K n /Au(111) junction under a given bias voltage V is extracted from the first-tier ADOs 29 , and then the dI/dV spectrum is obtained with a finite difference analysis.
We employ two approaches to theoretically determine the Kondo temperature T K of K n /Au(111) composites. The first one makes use of an empirical scaling relation between the zero-bias conductance G dI dV À Á V¼0 of a Kondo impurity and the environmental temperature T as follows 52,53 : Here, G 0 is the conductance at T → 0, g b is the background conductance due to electron transport through the non-Kondo states, and s is a parameter whose value depends on the local spin state of the impurity. Particularly, s = 0.22 is adopted as suggested by previous calculations for a spin-1 2 impurity 52-54 . The Kondo temperature T K can thus be determined by fitting the calculated G versus T to Eq. (6).
The second approach is based on the relation of Eq. (2) from the Fermi liquid theory 55 . This requires the calculation of the full dI/dV spectra for all the tip/K n /Au (111) junctions. Since the environmental temperature T adopted in the calculation is somewhat higher than the experimental counterpart, calculation needs to be done at a series of temperatures. The Kondo temperature T K is then extracted by fitting the resulting Γ versus T to Eq. (2) while taking λ as a tunable parameter (see Supplementary Note 10 for more details).

Data availability
The data that support the findings of this study are available from the authors on reasonable request; see author contributions for specific data sets.

Code availability
The codes that were employed in this study are available from the authors on reasonable request.