Bulk-mediated interaction between impurities in 1D atomic chains

TL;DR

This paper investigates how electron-mediated interactions between impurities in 1D atomic chains depend on separation and doping, revealing sign changes and magnitude variations that inform understanding of complex material interactions.

Contribution

It introduces combined numerical and analytical methods to analyze impurity interactions in 1D chains, providing new insights into bulk-mediated effects.

Findings

Interaction energy varies with impurity separation and doping.

Sign of interaction can change depending on parameters.

Provides a simplified model for understanding complex materials.

Abstract

A combination of numerical and analytical methods is employed to study a one-dimensional chain of identical atoms with adsorbates. We show that the electron-mediated interaction energy between two impurities can change sign and magnitude depending on the adatom-adatom separation, as well as the system doping. By focusing on this simple system, we provide insight into the bulk-mediated interaction for more complex materials.

Equations

$${\mathcal{E}}_q=\frac{E-2t\cos \left(qd\right)}{1+2g\cos \left(qd\right)},$$

$$H=$$

$$+$$

$$S$$

$$+\sum _{{\omega }_n}\left(\begin{array}{c}{\bar{\varphi }}_{{\omega }_n},{\bar{\psi }}_{{\omega }_n}\end{array}\right)\left[-1G_0^{-1}\left(i{\omega }_n+\mu \right)\right]\left(\begin{array}{c}{\varphi }_{{\omega }_n}\\ {\psi }_{{\omega }_n}\end{array}\right)$$

$$+\frac{V}{\sqrt{N}}\sum _{q\in \mathrm{BZ}}\sum _{{\omega }_n}{\bar{\gamma }}_{q,{\omega }_n}\left[{\varphi }_{{\omega }_n}+{\psi }_{{\omega }_n}{e}^{-iqdl}\right]$$

$$+\frac{V}{\sqrt{N}}\sum _{q\in \mathrm{BZ}}\sum _{{\omega }_n}\left[{\bar{\varphi }}_{{\omega }_n}+{\bar{\psi }}_{{\omega }_n}{e}^{iqdl}\right]{\gamma }_{q,{\omega }_n}.$$

$$\mathcal{Z}={\mathcal{Z}}_{\mathrm{C}}\int \mathcal{D}\left(\dots \right)e^{-{\sum }_{{\omega }_n}\left(\begin{array}{c}{\bar{\varphi }}_{{\omega }_n},{\bar{\psi }}_{{\omega }_n}\end{array}\right)\left[-{\mathcal{G}}_{{\omega }_n}^{-1}\right]\left(\begin{array}{c}{\varphi }_{{\omega }_n}\\ {\psi }_{{\omega }_n}\end{array}\right)},$$

$${\mathcal{G}}_{{\omega }_n}^{-1}$$

$$-V^2{\Xi }^l\left(i{\omega }_n+\mu \right)\times \left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$$

$${\Xi }^m\left(z\right)$$

$$=-\frac{t+gE}{2{\left(t+gz\right)}^2}\frac{{\left(\Lambda -\sqrt{\Lambda -1}\sqrt{\Lambda +1}\right)}^{|m|}}{\sqrt{\Lambda -1}\sqrt{\Lambda +1}},$$

$$\Lambda$$

$${\mathcal{G}}_{\omega }^{-1}=\left(\omega -ϵ\right)1+i\pi V^2{D}_{\omega }\left(\begin{array}{cc}1& e^{i|q_{\omega }dl|}\\ e^{i|q_{\omega }dl|}& 1\end{array}\right).$$

$${\mathcal{G}}_{\omega }^{-1}$$

$$+i\pi V^2{D}_{\omega }\left(\begin{array}{cc}1& \cos \left|q_{\omega }dl\right|\\ \cos \left|q_{\omega }dl\right|& 1\end{array}\right).$$

$${\mathcal{G}}_{\omega }^{11}=\frac{G_0^{-1}\left(\omega \right)-V^2{\Xi }^0\left(\omega \right)}{{\left[G_0^{-1}\left(\omega \right)-V^2{\Xi }^0\left(\omega \right)\right]}^2-{\left[V^2{\Xi }^l\left(\omega \right)\right]}^2}$$

$${\mathcal{Z}}_{\mathrm{A}}$$

$${\mathcal{Z}}^j$$

$$F_{\mathrm{A}}=-\frac{2}{\beta }\sum _{{\omega }_n}\ln \left[\beta \left(i{\omega }_n+\mu -ϵ-V^2{\Xi }^0\left(i{\omega }_n+\mu \right)\right)\right]$$

$$-\frac{1}{\beta }\sum _{{\omega }_n}\ln \left[1-{\left(\frac{V^2{\Xi }^l\left(i{\omega }_n+\mu \right)}{i{\omega }_n+\mu -ϵ-V^2{\Xi }^0\left(i{\omega }_n+\mu \right)}\right)}^2\right].$$

$$F_I=-\int \frac{d\omega }{2\pi }\ln \left[1-{\left(\frac{V^2{\Xi }^l\left(i\omega +\mu \right)}{i\omega +\mu -ϵ-V^2{\Xi }^0\left(i\omega +\mu \right)}\right)}^2\right].$$

$$F_I=E_b^{H_0{H}_l}-E_b^{H_0}-E_b^{H_l},$$

Full text

I Introduction

Impurities in a material can interact via their host’s electrons when the impurity energy levels hybridize with the electronic states of the bulk. One example is spin-spin interaction, commonly studied using the RKKY formalism, which can lead to (anti-) ferromagnetic ordering. Another is the potential interaction whose nature guides the arrangement of the adatoms by causing them either to cluster or disperse.

With the rise of two-dimensional materials, the interest of the condensed matter community turned to graphene-hosted adatoms. In their pioneering work, the authors of Ref. Shytov et al., 2009 claimed that the nature of the impurity-impurity interaction was determined by the sublattice configuration. The prediction was that atoms on the same sublattice repel while those on the opposite sublattices attract. More recent studies (LeBohec et al., 2014; Agarwal and Mishchenko, 2019), however, showed that while the sublattice arrangement is important, there are also other factors that play a role in the interaction. It was demonstrated, for example, that the interaction can change sign with the separation between the impurities and their on-site energies. First principles calculations have also shown that adsorbed atoms can attract regardless of the host sublattice. (González-Herrero et al., 2016)

Despite the relative simplicity of the graphene lattice, the mathematics involved in calculating the interaction between impurities analytically is substantial and this can obscure the physical picture. This complexity has, in part, led to the conflicting results of earlier publications. In order to gain a better understanding of the physical mechanisms behind the impurity-impurity interaction, it is useful to focus on a system that is more tractable mathematically.

Here, we consider a one-dimensional chain of identical atoms with two adsorbed impurities. By using a combination of analytical and numerical methods, we show that even this simple system exhibits a rich behavior where the interaction can change sign based on the doping of the system, as well as the separation between impurities. The simplicity of the model allows us to employ a non-perturbative approach to construct an intuitive picture of the processes involved in the interaction which can be generalized to other, more complex, systems.

II Model

II.1 Atomic Chain

The system in question is composed of two coupled components: the atomic chain and the adsorbed impurities. To model the chain, we use the tight-binding formalism with nearest neighbor hopping for an equally-spaced arrangement of single-orbital atoms, see Fig. 1. The dispersion for such a system is given by

[TABLE]

where $-t$ is the hopping parameter, $g$ is the orbital overlap between the neighboring atoms, $E$ is the on-site energy, and $d$ is the interatomic distance.

In this work, we compare analytic results with ab initio density functional theory (DFT) calculations. Therefore, we need to ensure that Eq. (1) adequately captures the relevant features of the ab initio band structure of our model system. To this end, we constructed a model 1D chain of C atoms with equal spacing of 1.29 Å.

The DFT band structure of the 1D chain is given in Fig. 1, with the doubly-degenerate band of interest lying between $-10$ and $6.5$ eV. The two bands are composed of $\pi$-bonded $p$ orbitals perpendicular to the direction of the chain. By fitting Eq. (1) to the ab initio data at $qd=0$, $\pi/2$, and $\pi$, we obtain $E$, $g$, and $t$ and plot the analytic formula along with the numerical results.

II.2 Chain-Impurity Coupling

We assume that each adatom has a single $s$ orbital with energy $\epsilon$. The impurity binds to a $p$ orbital of an individual chain atom with the hybridization energy $V$. Positioning one impurity at the 0th atom in the chain and another at the $l$th atom (see Fig. 1) yields the following Hamiltonian:

[TABLE]

where $N$ is the number of states in the one-dimensional Brillouin zone. The operators $c_{q}$ correspond to the chain states, while $a$ and $b$ are the impurity operators.

Using Eq. (2), we can write down the Matsubara-frequency action:

[TABLE]

By comparing $S$ to $H$, one can deduce that $\gamma$ fields correspond to the states in the atomic chain, and $\phi$ and $\psi$ are $a$ and $b$ impurities, respectively. $\Gamma^{-1}_{q}\left(i\omega_{n}+\mu\right)=i\omega_{n}+\mu-\mathcal{E}_{q}$ is the inverse Green’s function for the adsorbate-free atomic chain, while $G_{0}^{-1}\left(i\omega_{n}+\mu\right)=i\omega_{n}+\mu-\epsilon$ is the inverse Green’s function for an isolated adatom.

Action $S$ is used to calculate the partition function $\mathcal{Z}=\int\mathcal{D}\left(\dots\right)e^{-S}$, where $\int\mathcal{D}\left(\dots\right)$ denotes the integration over all the fields. In the following sections, we will employ $\mathcal{Z}$ to explore the impurity-chain coupling and how it leads to the impurity-impurity interaction.

III Green’s and Spectral Functions

When the impurity energy levels hybridize with the bulk, their spectral function becomes broadened. This spectral weight reorganization plays a key role in the impurity interaction. Therefore, we perform a careful analysis of the adatom spectral function starting with the calculation of the impurity Green’s functions.

III.1 Impurity Green’s Function

According to Eqs. (2) and (3), the two adsorbates do not interact directly. Instead, they both couple to the atomic chain which performs a mediating role. The interaction can be made explicit using the Green’s function obtained from the partition function $\mathcal{Z}$. By integrating out the $\gamma$ (atomic chain) fields in the expression of $\mathcal{Z}$, we obtain

[TABLE]

where $\mathcal{Z}_{\mathrm{C}}$ is the partition function for the atomic chain without the adsorbates. The integral on the right hand side is the partition function for the two adatoms $\mathcal{Z}_{\mathrm{A}}$ with the interaction due to the chain included. $\mathcal{G}_{\omega_{n}}$ is the desired $2\times 2$ Green’s function for the coupled impurities whose inverse is

[TABLE]

and

[TABLE]

To appreciate the physical significance of the $\Xi^{m}\left(z\right)$ term, notice that it is the Fourier transform of the chain free-particle propagator $\Gamma_{q}\left(z\right)$. As such, $\Xi^{m}\left(\mathcal{E}+i0\right)$ is the real-space propagator for a particle with energy $\mathcal{E}$ inside the adatom-free chain. For $\mathcal{E}$ outside the band, $|\Lambda|>1$ so that $\Xi^{m}\left(\mathcal{E}+i0\right)$ is real and monotonically decaying with $|m|$ because $\Lambda-\sqrt{\Lambda-1}\sqrt{\Lambda+1}<1$. Physically, this is the consequence of there being no propagating states outside the band.

If $\mathcal{E}$ is inside the band, on the other hand, $\Lambda=\cos\left(q_{\mathcal{E}}d\right)-i0$ so that $\Xi^{m}\left(\mathcal{E}+i0\right)=-i\pi D_{\mathcal{E}}e^{i\left|q_{\mathcal{E}}dm\right|}$, where $q_{\mathcal{E}}$ is the momentum of the chain state with energy $\mathcal{E}$, and $D_{\mathcal{E}}$ is the density of states. Hence, the propagation by distance $dm$ is given by the plane wave phase multiplied by the density of states at $\mathcal{E}$. The amplitude of the propagation does not decay with distance as it does in higher dimensions. Since the propagator functions as the coupling term between the two impurities (see Eq. (5)) the magnitude of the interaction between the adatoms at a given energy does not diminish with increased separation.

In the single-impurity case, the retarded Green’s function is given by the diagonal terms of Eq. (5) with $i\omega_{n}+\mu\rightarrow\omega+i0$ so that $\mathcal{G}_{\mathcal{E}}=\left[G^{-1}_{0}\left(\omega+i0\right)-V^{2}\Xi^{0}\left(\omega+i0\right)\right]^{-1}$. For $\omega$ within the band, this is a standard single-state Green’s function with frequency-dependent broadening due to the imaginary self-energy correction $iV^{2}D_{\omega}$. Outside the band, $\mathcal{G}_{\omega}$ has two poles: one above and one below the band, corresponding to two localized states. Because these two states are non-propagating, they play a minor role in the impurity interaction. Moreover, their spectral weight decreases sharply as $\epsilon$ is moved away from band edges. Therefore, we will neglect them in our qualitative discussion. In the following section, we explore how interaction between two states with simple broadening gives rise to a highly non-trivial spectral function.

III.2 Spectral Function

In the two-impurity case, the inverse of the retarded Green’s function for $\omega$ inside the band, in accordance with Eq. (5), is

[TABLE]

As discussed above, the diagonal term $i\pi V^{2}D_{\omega}$ describes the coupling between a single impurity and the chain. The consequence of this is that the impurity state becomes distributed over a range of composite impurity-chain states of energy $\omega$. The off-diagonal term in Eq. (8) is the coupling between these composite states for the two impurities. This interaction lifts the degeneracy between the mixed impurity-chain states for both adatoms at each $\omega$ and results in two split energy levels. This mixed-level splitting lies at the heart of the impurity interaction.

It is illustrative to separate the real and imaginary parts of Eq. (8):

[TABLE]

Focusing on the real part, we see that the adatom-adatom coupling splits the hybridized impurity-chain levels by $\sim\pi V^{2}D_{\omega}\sin\left|q_{\omega}dl\right|$. The oscillatory nature of the coupling term causes the split levels originating from the hybridized states at different $\omega$’s to periodically bunch up in energy.

To explore this periodic splitting, let us treat the density of states $D_{\omega}$ as approximately constant for some range of $\omega$. With this, for a given value of $l$, the amount of level splitting at each $\omega$ is dictated by the argument of the sine function in Eq. (9). Let $W$ be the energy range of non-split states for which one of the split states ends up at the same energy, forming a spectral peak. Very crudely, $W$ can be approximated from $\left(q_{\omega+\frac{W}{2}}-q_{\omega-\frac{W}{2}}\right)dl\sim\pi$. If $W$ is much smaller than the band width of the chain, we get $W\sim\pi/\left(q^{\prime}_{\omega}dl\right)$. Hence, the spacing between the energy ranges with high density of split states is inversely proportional to $l$ with the amount of splitting not affected by the separation. In addition, flatter bands yield smaller $W$ due to the $q^{\prime}_{\omega}$ term in the denominator.

In order to visualize this periodic concentration of states, we plot the spectral function $\mathcal{A}_{\omega}=-2\,\mathrm{Im}\left[\mathcal{G}^{11}_{\omega}\right]$, where

[TABLE]

is the diagonal element of the Green’s function in Eq. (5) with $i\omega_{n}+\mu\rightarrow\omega+i0$, see Fig. 2. The spectral function remains centered around $\epsilon$ and displays oscillations whose period decreases with the impurity separation, in agreement with our discussion.

IV Adatom Interaction Energy

With the qualitative picture of the two-impurity coupling established, we move to their interaction energy. The first step is to integrate Eq. (4) over the remaining fields to obtain

[TABLE]

where $G^{-1}\left(i\omega_{n}+\mu\right)$ is the diagonal term in Eq. (5). $\mathcal{Z}^{\pm}$ are the partition functions for the interaction-split impurity-chain states. From this, one can extract the free energy for the interacting adatoms by recalling that $F=-\left(\ln\mathcal{Z}\right)/\beta$:

[TABLE]

The first term is the energy of the non-interacting adatoms, coupled to the atomic chain. By subtracting it from $F_{A}$, we are left with the second term, which we identify as the interaction energy between the adatoms. The Matsubara frequency summation in Eq. (12) is easiest to perform at $T=0$, where it can be replaced by an integral:

[TABLE]

Before computing the interaction energy from Eq. (13), let us perform a qualitative analysis using Eq. (9). At $T=0$, the energy of the system is calculated by summing the energies of the filled state. For the states located far from the Fermi level, both split and non-split hybridized states are either filled or empty regardless of the separation $l$. Neglecting the finer point of variable spectral weight, this means that they do not contribute to the interaction energy as $l$ (and, therefore, the amount of splitting) is changed. If the split states are close to the Fermi level, however, changing $l$ can cause them to become filled or empty, modifying the energy of the system. Because $l$ enters inside a periodic trigonometric function, the total energy is periodic in $l$. Finally, since the interaction energy is defined as the difference between the coupled and independent-impurity energies, it exhibits the same $l$-periodicity.

From Eq. (9), one can see that the wavenumber of these energy oscillations is $2q_{\mu}d$, since only the levels close to the Fermi level contribute. The factor of $2$ comes from the fact that the two split levels move in opposite directions as $l$ varies with wavenumber $k_{F}d=q_{\mu}d$. From this, for $k_{F}d<\pi/2$, the period of oscillations is $\pi/\left(k_{F}d\right)$, reminiscent of the Friedel oscillations. At $k_{F}d=\pi/2$, the oscillatory period is 2 so that the maxima and minima of $F_{I}$ alternate with every $l$. If $k_{F}d$ exceeds $\pi/2$, the period becomes less than $2$. In this case, the energy oscillation wavenumber is outside its first Brillouin zone since the period is smaller than lattice spacing. Shifting the wavenumber back to the first BZ gives the period $\pi/\left(\pi-k_{F}d\right)$.

If the period of the $F_{I}$ oscillations is much larger than the interatomic separation in the chain, the interaction energy undergoes a smooth variation with $l$. As the period approaches 2 ($k_{F}d\approx\pi/2$), the granularity of $l$ prevents the smooth oscillations of $F_{I}$. Instead, they acquire the form of beats with the period $\pi/\left|\pi-2k_{F}d\right|$. This can be visualized by plotting a sine curve on a lattice where the period of the function is slightly larger or smaller than twice the lattice constant. With this in mind, we plot the interaction as a function of $l$ in Fig. 3. To illustrate the smooth variation and the beat-like pattern of the interaction energy, we choose two values of $\mu$ for which the smooth period and the beat length is 10 lattice constants. These $\mu$’s are obtained by solving for $k_{F}$ using the expressions above and then calculating the chemical potential from Eq. (1).

One can see that the amplitude of energy oscillations decreases with increasing $l$. This does not contradict the point we made earlier about the magnitude of the propagator not changing with the separation. The amplitude of the energy oscillations is related to the number of states that get filled or emptied with varying $l$. This can be understood as follows. When talking about the peaks of the two-impurity spectral function, we introduced a quantity $W$. This quantity describes the energy range of non-split states that contribute to a spectral peak and can also be thought of as the spacing between the peaks. In other words, larger $W$ results in more states in a peak and, therefore, larger energy oscillation as these are filled and emptied with varying $l$. As was shown above, $W$ decreases as $l$ gets bigger, reducing the amplitude of the energy oscillations.

In addition to the smaller amplitude of energy oscillations, the reduction of peak spacing with increased $l$ also makes the system more sensitive to changes in $\mu$ (see Fig. 4). For large $W$, small variation in $\mu$ keeps the Fermi level within the same peak of the spectral function. As $W$ is reduced, changing $\mu$ causes it to cross multiple spectral peaks, leading to oscillations of the interaction energy.

Finally, we address the role of the adatom energy level $\epsilon$. Because of the broadening, its position relative to $\mu$ determines the number of states that participate in the energy oscillations. Thus, having it closer to $\mu$ leads to a larger amplitude, as seen in Fig. 3.

In order to validate the theoretical model, as well as to gain insight into how the adsorbates interact in a concrete case, we create a model 1D C-chain with 2 H adsorbates, and determine $F_{I}$ ab initio. We note that, because of the non-decaying propagator, an infinite (periodic) chain would result in each H adsorbate interacting with an infinite number of other H adsorbates. We therefore choose to model the system as an isolated 36 C-atom fragment, with one H adsorbate fixed to an edge while the location of the second is varied along the chain. The interaction energy of the impurities can be calculated as

[TABLE]

where $E_{b}^{H_{0}H_{l}}$ is the binding energy of the two-impurity system with H adsorbates on C atoms at positions [math] and $l$, while $E_{b}^{H_{0}}$ and $E_{b}^{H_{l}}$ represent the binding energies of one-impurity systems with the H adsorbate at position [math] and $l$, respectively.

Fig. 5 shows $F_{I}$ as a function of $l$ for four doping levels. Exactly at half filling, the value of $F_{I}$ oscillates with a period equal to $2$, while its amplitude decreases with $l$, in agreement with the discussion above. As $\mu$ is lowered, the period becomes slightly larger than 2, leading to the formation of a beat pattern, consistent with the picture shown in Fig. 3. When $\mu$ is moved close to the bottom of the $\pi$ bands, $k_{F}$ decreases sufficiently to yield the period of energy oscillations much larger than the atomic spacing, leading to a smoother variation with $l$.

We observe a strong qualitative agreement between the DFT and analytical results, as seen from Figs. 3 and  5. From the quantitative standpoint, the results are of the same order of magnitude. The differences in the energy scales can be attributed to two main factors. First, the hopping parameter $V$ used in the analytical calculations was not calculated exactly. Instead, we chose a representative value with the right order of magnitude when compared to the C-C hopping element. In addition, the DFT results were obtained using a finite chain where one of the adsorbates was hosted by an edge atom. As such, finite-size effects play a role in determining $F_{I}$.

V Conclusion

We have performed a detailed study of bulk-mediated interaction between adsorbates in 1D atomic chains. We have shown using analytical and ab initio methods that the interaction energy between these adsorbates exhibits an oscillatory behavior which can be tuned by changing the doping of the system. The qualitative understanding gained from studying this simple system can be applied to more complicated configurations. More specifically, our results suggest an effective approach to further exploring adatom interactions in 2D materials where interaction control by gating could potentially be used for adatom manipulation and assembly.

VI Acknowledgements

The authors acknowledge the National Research Foundation, Prime Minister Office, Singapore, under its Medium Sized Centre Programme. A. R. acknowledges the support by Yale-NUS College (through grant number R-607-265-380-121). K. N. and S.Y. Q. acknowledge support from Grant NRF-NRFF2013-07 from the National Research Foundation, Singapore. Computations were performed on the NUS CA2DM cluster.

Appendix A Computational Methods

DFT calculations were performed in VASP Kresse and Furthmuller (1996) using PBE Perdew et al. (1996) PAW potentials Blochl (1994); Kresse and Joubert (1999) and a kinectic energy cutoff of 400 eV for wavefunctions. In the case of the 1D C chain (Fig. 1), periodic images along the directions orthogonal to the chain axis were separated by 15 Å of vacuum in order to avoid spurious Coulomb interactions. The charge density of the perodic chain was calculated with a uniform Brillouin zone sampling of $36\times 1\times 1$ k-points. The corresponding bandstructure was computed at 100 evenly-spaced k-points between each high symmetry point. In the case of the isolated 36 C-atom fragment, a 15 Å vacuum was also imposed in the direction along the chain. Variations in $\mu$ were achieved by removing electrons from the system while applying a compensating background charge.

Numerical integration and plotting were performed using Julia programming language. Bezanson et al. (2017) The code can be found at https://github.com/rodin-physics/chain-impurity-interaction.