# Power and Skew Reduction Using Resonant Energy Recycling in 14-nm FinFET
Clocks
††thanks: This work was supported in part by the Rezonent Inc. under Grant
CORP-0061 and the UMBC Startup grant.

Dhandeep Challagundla, Mehedi Galib  Computer Science and Electrical
Engineering
University of Maryland Baltimore County
Baltimore, Maryland, USA
{vd58139<EMAIL_ADDRESS>Ignatius Bezzam  IC Design Group
Rezonent Inc.
Milpitas, USA
<EMAIL_ADDRESS>Riadul Islam  Computer Science and Electrical Engineering
University of Maryland Baltimore County
Baltimore, Maryland, USA
<EMAIL_ADDRESS>
###### Abstract

As the demand for high-performance microprocessors increases, the circuit
complexity and the rate of data transfer increases resulting in higher power
consumption. We propose a clocking architecture that uses a series LC
resonance and inductor matching technique to address this bottleneck. By
employing pulsed resonance, the switching power dissipated is recycled back.
The inductor matching technique aids in reducing the skew, increasing the
robustness of the clock network. This new resonant architecture saves over 43%
power and 91% skew clocking a range of 1–5 $GHz$, compared to a conventional
primary-secondary flip-flop-based CMOS architecture.

###### Index Terms:

Clock skew, LC resonance, clock tree architecture, pulsed flip-flops, Power
consumption.

##  I Introduction

Power consumption is one of the major problems faced in the high-performance
microprocessor industry [1, 2, 3]. The need for an increase in performance has
steered the operating frequencies higher, resulting in an increased complexity
among the microprocessor designs [2, 4, 5]. This higher power led designers to
constantly come up with innovative techniques to reduce the power while trying
to meet all the design constraints that impact the performance [6, 7, 8, 9]. A
significant portion of dynamic power consumed in a high-frequency design is
due to the switching activity in the clock network [10]. To address this,
several low power techniques such as dynamic voltage and frequency scaling
(DVFS) [11], clock gating [12] and LC resonant clocking [13, 14, 15, 16, 17],
current-mode clocking [18, 19, 20] are commonly used. Among them, inductor-
based LC resonant clocking techniques have great potential to save switching
power due to their constant phase and magnitude.

Figure 1: LC resonance topologies to reduce dynamic power consumption (a)
series resonance topology can address wide frequency band (b) parallel
resonance topology can address very narrow frequency band.

There are several LC resonant techniques, such as parallel resonance [13]
(please see Fig. 1(a)), intermittent resonance [15] and series resonance [14,
21]. However, most of the LC resonant techniques suffer from a limitation of
narrow frequency band and high skew. Moreover, most of the industry-standard
electronic design automation (EDA) tools do not explicitly support integrating
LC resonance in the clock tree architecture. Additionally, designing a
resonant clock architecture requires the designer to have multiple domain
expertise due to the non-linear behavior of inductors.

To overcome the narrow frequency band, researchers utilize series resonance
techniques [14], as shown in Fig. 1(b). This approach uses an inductor placed
in the discharge path to store the dissipated energy in the form of a magnetic
field. This energy is recycled in the next rising clock edge.

To enable a resonant clock architecture, we need resonant flip-flops (FFs) for
synchronous circuits. However, they occupy a substantial chip area consuming
high power. Researchers proposed many low-power flip-flops, however, not
suitable for resonant operations [22, 23, 24]. In this research, we propose
several conventional register-based pulsed FFs suitable for series resonance
and reduce the overall power consumption in the clock network.

Besides power, skew plays a critical role in enabling high-frequency
operation. To reduce the skew generated by clock trees, we introduce the first
inductor tuning technique to match the series resonance inductor with the load
capacitance of the clock tree, which generates pulses with equal resonant
frequencies. As the resonant frequency depends on the series inductor and the
load capacitance, we generate a constant pulse-width for a wideband (WB) of
input clock frequencies. Therefore, calibrating the inductors once for a clock
tree architecture enable WB frequency operation.

This paper proposes a clocking architecture to recycle the power dissipated by
the synchronous elements using series resonance technology and improve clock
skew by adapting the inductor tuning technique. In particular, the main
contributions of this work are:

  * •

An architecture to recycle the power dissipated by synchronous elements using
series resonance.

  * •

A novel pulse generator with dual-rail booster using series resonance.

  * •

A set of pulsed register-based FFs that exploits the behavior of pulsed series
resonance.

  * •

An inductor tuning technique to compensate for the skew in the clock tree
architecture.

##  II Background

In a traditional clocking method, half of the switching power is utilized to
charge a capacitive node when the clock transitions from 0-to-1. The other
half of switching power is dissipated in the discharge cycle when the clock
transitions from 1-to-0. The pulsed series resonance (PSR) technique recycles
this dissipated energy by placing an inductor in the discharging path. LC
resonant is most widely used among several energy recycling techniques as it
precisely replicates conventional CMOS clocking. However, it suffers from a
higher slew rate while demonstrating great savings in dynamic power
consumption [14] [25].

Due to resonance, the free energy swing obtained as a result of recycling
energy, is the difference between resonant high output ($V_{OH}$) and resonant
low outupt ($V_{OL}$) [14] can be expressed as

$V_{OH}-V_{OL}=\frac{V_{DD}}{2}(1+e^{-\pi/Q})-\frac{V_{DD}}{2}(1-e^{-\pi/2Q})$
(1)

where $Q$ is the quality factor of the inductor, which is given by
$Q=\sqrt{{L}/{(CR^{2})}}$. We utilize an external power source to pull the
output from $V_{OH}$ to $V_{DD}$. The resonant frequency at which the inductor
resonates with the load capacitance can be expressed as
$f_{RES}=\dfrac{1}{2\pi\sqrt{LC}}$.

In this work, we proposed several pulsed-type resonant FFs with resonant clock
trees for PSR operation. Several prior works have focused on low-power FF
designs. In [23] an 18T FF was designed to achieve 40% improvement in
energy/cycle compared to primary-secondary FF (PSFF). However, to mitigate
voltage degradation caused due to the non-complementary topology, the authors
used a poly bias technique [22], which requires extra design effort. In [22]
an 18T single-phase clocked FF was designed for low power operation. It showed
68% lower power consumption at 0.6V supply but had functionality issues when
the voltage was scaled, as reported in [26, 27]. For reliability and robust
operation, we implement widely used traditional pulsed register-based and true
single-phase clock (TSPC)-based FFs [28].

##  III Proposed Clock Architecture

The proposed architecture comprises a pulse generator, clock drivers utilizing
on-chip inductors, and pulsed registers, as shown in Fig. 2.

Figure 2: The proposed wideband resonant clock tree architecture consists of a
system clock source as the root, followed by a pulse generator, multiple PSR
drivers with on-chip inductors, clock gaters, clock buffers, and finally,
various sets of resonant pulsed FFs in the leaf nodes.

###  III-A Pulse Generator

The pulse generator is depicted in Fig. 2 takes the input from the clock
source and generates a pulse with boosted amplitude. The series inductor
$L_{1}$ and the matching capacitance $C_{1}$ generate a delay of
$T_{d}=\pi\sqrt{L_{1}C_{1}}$. The clock and the delayed signal are fed into an
XNOR gate to generate a pulse at both clock edges with a pulse width $T_{d}$.
Now, a voltage doubler circuit is employed to invert the generated dual
triggered pulse resulting in a boosted signal $V_{SR}$. The voltage doubler
circuit uses the pulsed series resonance technique to generate a boosted
signal. When the $\overline{V_{sr}}$ is low, the PMOS transistors $M_{1}$ and
$M_{3}$ are “ON,” and the inductor resonates with the load capacitance $C_{2}$
and the additional PSR capacitance. For large load capacitances, the value of
the series inductors is quite small. The inductor in the voltage doubler
circuit can be adjusted according to the load of the pulse generator to
produce a boosted signal $V_{SR}$. We use a dual-rail booster circuit to
reduce the power consumed by the voltage doubler by decreasing the resistance
of the pull-up network.

###  III-B Pulsed Series Resonance Driver

The boosted signal $V_{SR}$ from the pulse generator stage is provided as the
input to multiple PSR drivers to generate a pulse signal $R_{CLK}$. The
inductors on the PSR drivers resonate with the capacitance of the tree to
generate a pulse signal that is traversed through many levels of transmission-
gate clock gaters and clock buffers. Since we provide a boosted $V_{SR}$
signal as the input, we obtain a rail-to-rail swing at the output of the PSR
driver, and it improves the robustness of the design. Then, the output signal
of the PSR driver $R_{CLK}$ is inverted and supplied to pulsed registers as
the clock input signal $Pclk$.

Figure 3: The proposed series resonant pulsed FFs use (a) A PSR to generate
pulse signal to drive the register stage, (b) 13T register, (c) pulsed
register, (d) TSPC register to implement three pulsed FFs.

###  III-C Proposed Resonant Pulsed Flip-Flops

To utilize the generated $Pclk$ pulse (from the PSR, as shown in Fig. 3(a)),
we propose a 13T pulsed FF (13TPFF), as shown in Fig. 3(b). It takes the input
data and inverts it to provide it to the transistors M2 and M3, respectively.
The M2 and M3 transistors drain are connected to the storage cells where the
data is stored as logic “1” or a logic “0”. If a “1” is stored in the
register, the value at $S=1$ and $S_{B}=0$. If a “0” is stored, the voltages
will be reversed.

When $Pclk$ is “0,” the transistor M1 is turned off, wherein the FF is in
hold/retain state, and the values of $S$ and $S_{B}$ are unaltered. Consider
the case when $Data=1$ and $Pclk=1$, the transistors M2 and M1 turn on
connecting the node $S_{b}$ to ground, which then discharges the node and
makes it 0, making $Q=1$ writing a “1” into the register. When $Data=0$ and
$Pclk=1$, the transistors M3 and M1 turn “ON” and write a “0” at node $Q$. The
FF has an active-low asynchronous reset. The M4 and M5 transistors are turned
“ON” and “OFF,” respectively, when the $Reset$ signal goes to low resulting in
a logic “1” at node $S_{b}$ that writes a “0” the output $Q$.

Along with a 13TPFF, we also design pulsed energy recovery FFs. The pulsed
resonant FF (PRFF) is based on the traditional latch, and the resonant TSPCFF
is based on the TSPC register [28]. Except, they use pulsed series resonance
to take advantage of the input pulse signal $V_{SR}$ to recycle energy. The
energy recovery FFs are positive edge-triggered with asynchronous active low
$Reset$ signal. The use of conventional registers makes an easy integration of
resonance clock trees into existing clock tree architectures.

###  III-D Skew Reduction Methodology

Skew is defined as the spatial variation in the arrival time of clock
transition at two different locations. There are several reasons for this
skew, and one such cause is different loads on clock drivers [28].

In Fig. 5, the eight different branches of the clock tree are having eight
different capacitances $C_{SR1}$, $C_{SR2}$, and up to $C_{SR8}$ due to on-
chip variation (OCV). This capacitance mismatch between different branches of
a clock tree will result in different clock arrival times. Each branch of the
clock tree represents a separate $LC$ resonant tank. Using the resonant
frequency $f_{RES}=1/2\pi\sqrt{LC}$, we match the inductors $L_{SR1}$ to
$L_{SR8}$ with the load capacitances $C_{SR1}$ to $C_{SR8}$, respectively, to
have equal frequencies. This inductor matching would result in equal frequency
signals in all the clock branches, thus, reducing the skew. The resonant
frequency independent of the input clock frequency will not be affected by
wide frequency band operation. The primary reason is the $delay=T_{d}$ of the
pulse generator circuit is independent of clock pulse width, and it works on
the clock edges. For all the clock frequencies less than the resonant
frequency $f_{RES}$, we can have the same inductor value that results in
reduced skew. Our results in Section IV-B2 also supported this claim.

##  IV  Results and Discussion

###  IV-A Experimental Setup

Figure 4: Illustration of Monte-Carlo simulation results for various FFs by
considering 5000 samples with $\pm 10\%$ length variation. (a) PSFF, (b) PRFF,
(c) TSPCFF, and (d) 13TPFF.

The proposed resonant clock tree architecture, as shown in Fig. 5 is
implemented using a standard 14nm FinFET technology. Conventional clock tree
architecture is used as a reference model to compare with the proposed
architecture. Each tree has eight clock drivers, 4K clock gaters, 8K clock
buffers, and 32K FFs. The traditional clock tree makes use of transmission
gate PSFFs, whereas the resonant clock tree uses pulsed FFs, as shown in Fig.
5(a) and Fig. 5(b), respectively. All the FFs layouts are compatible with a
standard cell height of 24 horizontal M2 tracks. All the simulations are
performed for frequencies ranging from 1 ${GHz}$ to 5 ${GHz}$.

Figure 5: Clock tree architectures used for functional simulations, (a)
conventional clock tree architecture with eight branches and different loads
totaling 32k FFs, (b) resonant clock tree architecture replicating the same
number of branches and loads as the conventional one.

Figure 6: Simulation waveforms show a CLK input of 0.5 ${GHz}$ is provided to
generate a 1 ${GHz}$ $Pclk$ clock to assert as a clock input to the pulsed
FFs.

Figure 7: 13TPFF layout is implemented using 14 nm FinFET technology following
standard cell height, which improves setup-time constraints compared to PSFF.
TABLE I: The proposed 13TPFF exhibits better set-up time than the PSFF and
better hold-time than TSPCFF and PRFF while consuming more dynamic power and
area; however, it consume lower static power than PSFF and enables power
saving in overall clock architecture.

Types of FF |  Normalised area | Delay ($ps$) | Static power ($pW$) | Dynamic power ($\mu W$)
---|---|---|---|---
C-Q | ts | th | D=0 | D=1 | 1GHz | 2GHz | 3GHz | 4GHz | 5GHz
PSFF | 1 | 32.5 | 14 | 2 | 1550 | 593 | 8.3 | 14.1 | 21 | 28 | 35.1
| PRFF
---
0.59 | 35.1 | -95 | 96 | 278 | 272 | 7.16 | 13.8 | 20.4 | 27.1 | 33.8
| TSPCFF
---
0.84 | 41.9 | -92 | 93 | 283 | 664 | 12.3 | 20.2 | 28 | 35.9 | 43.7
| 13TPFF
---
1.75 | 37.3 | -25 | 60 | 501 | 538 | 16.2 | 31.1 | 46 | 61 | 76

TABLE II: Our proposed PRFF outperforms all the FFs and consumes 43% less
power than the conventional PSFF with ${11\times}$ improvement in skew, while
TSPCFF and 13TPFF consume 26% and 20.2% less power than the PSFF,
respectively.

Types of flip-flop | Skew ($ps$) | Total tree power ($mW$)
---|---|---
1GHz | 2GHz | 3GHz | 4GHz | 5GHz
PSFF | 51.1 | 30.8 | 60.6 |  89.2 |  116 | 138
PRFF | 4.61 | 17.4 | 34.1 | 50.3 | 65.7 | 78.7
TSPCFF | 2.05 | 22.2 | 43.6 | 64.6 | 84.4 | 102
13TPFF | 3.92 | 23.8 | 46.8 | 69.4 | 90.7 | 110

###  IV-B Implementation Results

####  IV-B1 Power and Performance Comparison of Registers

For measuring the performance and functionality of the proposed 13TPFF under
process variations, we consider 5000 samples of CLK-to-Q ($t_{c-q}$) delay
using Monte-Carlo simulation. $\pm 10$% variation in the length of all devices
is considered while performing the simulations. The $t_{c-q}$ delay
distributions of PSFF, PRFF, TSPCFF, and 13TPFFs are shown in Fig. 4(a), Fig.
4(b), Fig. 4(c), and Fig. 4(d), respectively. Among all the resonant FFs, the
PRFF has lowest mean $t_{c-q}$ of 35 ps with standard deviation of 0.066 $ps$.

The normalized layout area, $t_{c-q}$, setup times $(t_{s})$, hold times
$(t_{h})$, and power for the FFs are listed in Table I. Among all the
competing FFs, the 13TPFF consumed the highest layout area of 9.62 ${um^{2}}$,
which is 1.75$\times$ the area of PSFF whose area is 5.151 ${um^{2}}$, and
2.9$\times$ the area of a PRFF whose area is 3.091 ${um^{2}}$. The proposed
13TPFF has a $t_{s}$ of -25 $ps$ and a $t_{h}$ of 60 $ps$ with a clock-to-q
delay of 37.3 $ps$. Empirically, pulsed register-based FFs exhibit negative
$t_{s}$, which tremendously impacts resolving $t_{s}$ related timing issues.
The PSFF has a $(t_{s})$ of 14 $ps$, $(t_{h})$ of 2 $ps$ and $t_{c-q}$ of 32.5
$ps$. The resonant PRFF has a better $(t_{s})$ of -95 $ps$ but has a high
$(t_{h})$ of 96 $ps$ which is similar to the TSPCFF with -92 $ps$ $(t_{s})$
and 93 $ps$ $(t_{s})$. However, the power consumed by the proposed 13TPFF is
2$\times$ higher than the PSFF. Among all the competing FFs, the PRFF consumes
the lowest dynamic static and dynamic powers.

####  IV-B2 Clock Tree

The functional simulation for the resonant clock is shown in Fig. 6. We
provide 0.5 ${GHz}$ clock as input clock source shown in Fig. 2. The output of
the pulse generator is a boosted signal $V_{SR}$ of 1 ${GHz}$ frequency. This
$V_{SR}$ signal is then provided to a PSR driver whose output is $R_{CLK}$.
This clock signals $IN_{1}$ along with the Data signal generates the output
$Q$, as shown Fig. 3(b). We compare the power and skew of the proposed clock
architecture, as shown in Table II for frequencies ranging from 1 ${GHz}$ to 5
${GHz}$. The power consumed by the proposed architecture while using 13TPFFs
at 1 ${GHz}$ frequency is 23.8 $mW$, compared to a conventional clock tree
architecture with PSFF that consumes 30.8 $mW$. The skew generated by the
conventional clock tree is 51.1 $ps$. As a result of inductor tuning, the skew
generated by the proposed resonant clock architecture is 3.92 $ps$. The
proposed architecture saves 22.7% power using the 13TPFFs compared to
conventional tree. The resonant clock tree using TSPCFFs has a skew of 2.1
$ps$ and saves 27.9% power, whereas the clock with PRFFs has a skew of 4.6
$ps$ and saves 43% power compared to conventional clock.

##  V  Conclusion

This paper proposed a resonant clock architecture to balance the skew and
recycle the power consumed. The proposed architecture with 13TPFF saves 22.7%
power with 92% lower clock skew than the conventional clock tree architecture.
Furthermore, it saves 43% power with a 91% skew reduction while using the PRFF
compared to conventional PSFF-based CMOS clock tree architecture in 14 nm
FinFET technology.

## References

  * [1] T. Fischer, S. Arekapudi, E. Busta, C. Dietz, M. Golden, S. Hilker, A. Horiuchi, K. A. Hurd, D. Johnson, H. McIntyre, S. Naffziger, J. Vinh, J. White, and K. Wilcox, “Design solutions for the bulldozer 32nm soi 2-core processor module in an 8-core cpu,” in _2011 IEEE International Solid-State Circuits Conference_ , 2011, pp. 78–80.
  * [2] A. Cunningham, “Don’t buy a desktop pc with one of intel’s newest processors—here’s why,” https://www.nytimes.com/wirecutter/blog/desktopc-pc-intel-newest-processors/, April 30, 2021.
  * [3] A. Kumar Mishra, D. Vaithiyanathan, and U. Chopra, “Design and analysis of ultra-low power 18t adaptive data track flip-flop for high-speed application,” _International Journal of Circuit Theory and Applications_.
  * [4] N. Kumar and D. P. Vidyarthi, “A novel energy-efficient scheduling model for multi-core systems,” _Cluster Computing_ , vol. 24, no. 2, pp. 643–666, 2021\.
  * [5] A. A. Khan, A. Ali, M. Zakarya, R. Khan, M. Khan, I. U. Rahman, and M. A. Abd Rahman, “A migration aware scheduling technique for real-time aperiodic tasks over multiprocessor systems,” _IEEE Access_ , vol. 7, pp. 27 859–27 873, 2019.
  * [6] H. Jeong, T. W. Oh, S. C. Song, and S.-O. Jung, “Sense-amplifier-based flip-flop with transition completion detection for low-voltage operation,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 26, no. 4, pp. 609–620, 2018.
  * [7] L. Touil, A. Hamdi, I. Gassoumi, and A. Mtibaa, “Design of low-power structural fir filter using data-driven clock gating and multibit flip-flops,” _Journal of Electrical and Computer Engineering_ , vol. 2020, 2020.
  * [8] P. CET, “Review of low power design techniques for flip-flops,” _International Journal of Pure and Applied Mathematics_ , vol. 120, no. 6, pp. 1729–1749, 2018.
  * [9] L. Cherif, M. Chentouf, J. Benallal, M. Darmi, R. Elgouri, and N. Hmina, “Usage and impact of multi-bit flip-flops low power methodology on physical implementation,” in _2018 4th International Conference on Optimization and Applications (ICOA)_ , 2018, pp. 1–5.
  * [10] J. M. Rabaey, _Low power design essentials_ , 2nd ed. Springer, Jan. 2009.
  * [11] K. Nowka, G. Carpenter, E. MacDonald, H. Ngo, B. Brock, K. Ishii, T. Nguyen, and J. Burns, “A 32-bit powerpc system-on-a-chip with support for dynamic voltage scaling and dynamic frequency scaling,” _IEEE Journal of Solid-State Circuits_ , vol. 37, no. 11, pp. 1441–1447, 2002.
  * [12] V. Tirumalashetty and H. Mahmoodi, “Clock gating and negative edge triggering for energy recovery clock,” in _2007 IEEE International Symposium on Circuits and Systems_ , 2007, pp. 1141–1144.
  * [13] F. u. Rahman and V. Sathe, “Quasi-resonant clocking: Continuous voltage-frequency scalable resonant clocking system for dynamic voltage-frequency scaling systems,” _IEEE Journal of Solid-State Circuits_ , vol. 53, no. 3, pp. 924–935, 2018.
  * [14] I. Bezzam, C. Mathiazhagan, T. Raja, and S. Krishnan, “An energy-recovering reconfigurable series resonant clocking scheme for wide frequency operation,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 62, no. 7, pp. 1766–1775, 2015.
  * [15] H. Fuketa, M. Nomura, M. Takamiya, and T. Sakurai, “Intermittent resonant clocking enabling power reduction at any clock frequency for near/sub-threshold logic circuits,” _IEEE Journal of Solid-State Circuits_ , vol. 49, no. 2, pp. 536–544, 2014.
  * [16] R. Islam, “Low-power resonant clocking using soft error robust energy recovery flip-flops,” _Electron Test_ , vol. 34, p. 471–485, 2018.
  * [17] P.-Y. Lin, H. A. Fahmy, R. Islam, and M. R. Guthaus, “Lc resonant clock resource minimization using compensation capacitance,” in _IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2015, pp. 1406–1409.
  * [18] R. Islam and M. R. Guthaus, “Cmcs: Current-mode clock synthesis,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 25, no. 3, pp. 1054–1062, 2017.
  * [19] M. Guthaus and R. Islam, “Current-mode clock distribution,” January 2014.
  * [20] R. Islam, H. Fahmy, P.-Y. Lin, and M. R. Guthaus, “Differential current-mode clock distribution,” in _International Midwest Symposium on Circuits and Systems (MWSCAS)_ , 2015, pp. 1–4.
  * [21] I. Bezzam and N. Rawat, “Digital circuits for radically reduced power and improved timing performance on advanced semiconductor manufacturing processes,” July 2021.
  * [22] Y. Cai, A. Savanth, P. Prabhat, J. Myers, A. S. Weddell, and T. J. Kazmierski, “Ultra-low power 18-transistor fully static contention-free single-phase clocked flip-flop in 65-nm cmos,” _IEEE Journal of Solid-State Circuits_ , vol. 54, no. 2, pp. 550–559, 2019.
  * [23] F. Stas and D. Bol, “A 0.4-v 0.66-fj/cycle retentive true-single-phase-clock 18t flip-flop in 28-nm fully-depleted soi cmos,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 65, no. 3, pp. 935–945, 2018.
  * [24] M.-Y. Tsai, P.-Y. Kuo, J.-F. Lin, and M.-H. Sheu, “An ultra-low-power true single-phase clocking flip-flop with improved hold time variation using logic structure reduction scheme,” in _2018 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2018, pp. 1–4.
  * [25] R. Islam, B. Saha, and I. Bezzam, “Resonant energy recycling sram architecture,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 68, no. 4, pp. 1383–1387, 2021.
  * [26] G. Shin, E. Lee, J. Lee, Y. Lee, and Y. Lee, “A static contention-free differential flip-flop in 28nm for low-voltage, low-power applications,” in _2020 IEEE Custom Integrated Circuits Conference (CICC)_ , 2020, pp. 1–4.
  * [27] H. You, J. Yuan, Z. Yu, and S. Qiao, “Low-power retentive true single-phase-clocked flip-flop with redundant-precharge-free operation,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 29, no. 5, pp. 1022–1032, 2021.
  * [28] J. M. Rabaey, A. Chandrakasan, and B. Nikolic, _Digital integrated circuits- A design perspective_ , 2nd ed. Prentice Hall, 2004.


# Substitutional tin acceptor states in black phosphorus

Mark Wentink  Department of Electronic and Electrical Engineering, University
College London, WC1E 7JE, London, UK  Julian Gaberle  Department of Physics
and Astronomy, University College London, WC1E 6BT, London, UK  Martik
Aghajanian  Departments of Materials and Physics and the Thomas Young Centre
for Theory and Simulation of Materials, Imperial College London, SW7 2AZ,
London, UK  Arash A. Mostofi  Departments of Materials and Physics and the
Thomas Young Centre for Theory and Simulation of Materials, Imperial College
London, SW7 2AZ, London, UK  Neil J. Curson  Department of Electronic and
Electrical Engineering, University College London, WC1E 7JE, London, UK
Johannes Lischner  Departments of Materials and Physics and the Thomas Young
Centre for Theory and Simulation of Materials, Imperial College London, SW7
2AZ, London, UK  Steven R. Schofield  Department of Physics and Astronomy,
University College London, WC1E 6BT, London, UK  Alexander L. Shluger
Department of Physics and Astronomy, University College London, WC1E 6BT,
London, UK  Anthony J. Kenyon  Department of Electronic and Electrical
Engineering, University College London, WC1E 7JE, London, UK
<EMAIL_ADDRESS>
###### Abstract

Nominally-pure black phosphorus (BP) is commonly found to be a p-type
semiconductor, suggesting the ubiquitious presence of impurity species or
intrinsic, charged defects. Moreover, scanning tunnelling microscopy (STM)
images of black phosphorus reveal the presence of long-ranged double-lobed
defect features superimposed onto the surface atomic lattice. We show that
both the p-type doping of BP and the defect features observed in STM images
can be attributed to substitutional tin impurities. We show that black
phosphorus samples produced through two common synthesis pathways contain tin
impurities, and we demonstrate that the ground state of substitutional tin
impurities is negatively charged for a wide range of Fermi level positions
within the BP bandgap. The localised negative charge of the tin impurities
induces hydrogenic states in the bandgap and it is the 2p level that sits at
the valence band edge that gives rise to the double-lobed features observed in
STM images.

London Centre for Nanotechnology, University College London, WC1H 0AH, London,
UK London Centre for Nanotechnology, University College London, WC1H 0AH,
London, UK London Centre for Nanotechnology, University College London, WC1H
0AH, London, UK

##  1 Introduction

Black phosphorus (BP) is a two-dimensional (2D) van der Waals layered material
and is considered a promising candidate for future (opto)electronic and
quantum-electronic devices 1, 2. BP holds a unique position within the family
of 2D van der Waals layered materials, being the only mono-elemental 2D
material besides graphite 3. Electrical and optical characterisation of
nominally pure samples of BP uniformly find the material to be electronically
p-type, dating back to the earliest measurements in the 1960s 4, 5, 6, 7, 8.
This suggests the ubiquitous presence of charged impurities or structural
defects in BP; however, the identity of these impurities/defects, and
therefore the source of the p-type doping in BP, remains unexplained. In
addition to causing shifts in the Fermi level, atomic point defects can
greatly affect the fundamental properties of crystalline materials; for
example they can act as nucleation centers for exotic phenomena such as charge
density waves or superconducting phases 9, 10. Thus, determining the atomic
structure of the intrinsic defects/impurities in BP is of critical importance
for both the fundamental understanding of BP, and for nanoelectronic
applications.

Black phosphorus is comprised of hexagonally bonded sheets of phosphorus
atoms, similar to the hexagonally bonded sheets of carbon atoms that form
graphene. However, because phosphorus has one extra valence electron than
carbon, the phosphorus atoms of BP are $sp^{3}$ hybridised with three bonds to
nearest neighbours and an electron lone pair, rather than $sp^{2}$ hybridised
with a single electron (per carbon atom) in a $\pi$ orbital as is the case for
graphene sheets. This results in the individual phosphorene sheets of BP
forming a buckled configuration, as indicated schematically in Fig. 1a. The
directions parallel and perpendicular to these rows are conventionally
referred to as the zig-zag and armchair directions, respectively. In addition
to this strong structural anisotropy, BP also exhibits both in-plane and out-
of-plane electronic anisotropy, with the in-plane anisotropy characterised by
strong changes in electron effective mass and mobility between the zig-zag and
armchair directions 11, 12.

The presence of point defects in nominally intrinsic BP has been confirmed by
atomic-resolution scanning tunnelling microscopy (STM) measuremts made by
multiple groups 7, 8, 13. The defects appear as long-ranged anisotropic
protrusions superimposed on the atomic lattice, and in the literature they
have largely been attributed to charged single-phosphorus vacancies 7, 8, 14.
However, this assignment is in disagreement with density functional theory
(DFT) calculations that demonstrate very small diffusion barriers for
monovacancies in BP 15, 16, 17, 18. These low barriers mean that monovacancies
can be expected to diffuse much faster than STM scanning times at room
temperature, and remain mobile even at 77 K. However, no evidence of diffusion
of the defects observed by STM have been reported 7, 8, 14, consistent also
with our observations here. DFT calculations also show that monovacancies are
unstable with respect to the formation of divacancies. 17 Since divacancies
are charge neutral, they are not expected to produce the long-ranged features
observed in STM images, nor can they explain the p-type doping observed in BP.
Hence vacancy defects do not provide an explanation for the experimentally
measured p-type electronic behaviour of BP, or the appearance of the defects
observed in STM experiments.

Here, we use a combination of atomic-resolution STM, scanning tunnelling
spectroscopy (STS), X-ray photoelectron spectroscopy (XPS), DFT, and tight-
binding simulations, to investigate the origin of p-type doping in BP and the
long-range defect features observed in STM images of BP. We find that both the
p-type doping and the defects seen in STM images can be attributed to
substitutional tin impurities in the lattice, which are present due to the use
of tin-based catalysts in BP growth, and that produce negatively-charged
acceptor states in BP. The characteristic double-lobed appearance of the
defects in STM images results mainly from imaging a 2p orbital that exists
just inside the bandgap at the valence band edge, with a lesser contribution
also from a 1s orbital that lies deep within the bandgap.

##  2 Methods

###  2.1 Experimental methods

Scanning Tunnelling microscopy and spectroscopy (STM/STS) measurements were
performed using a commercial Scienta Omicron system operating at 77 K with a
base pressure below 10-10 mbar. Black phosphorus flakes were purchased from
SmartElements and 2DSemiconductors. Samples with large single crystal areas
were selected for STM measurement and mounted onto stainless steel sample
holder plates using Epotek H21D conductive epoxy. Samples were outgassed in
vacuum at 150∘C overnight, before being cleaved in-situ to expose an
atomically-clean surface. STM tips were chemically etched from tungsten wire
and had their surface oxide removed by electron bombardment in the STM vacuum
chamber. STS was performed by measuring the current as a function of the
sample voltage with a constant tip-sample separation. X-ray photoemission
spectroscopy (XPS) measurements were performed on a Thermo K-alpha system
operating at a base pressure of 10-9 mbar. Samples were cleaved ex-situ before
rapid transfer into the loadlock and pumping down to 10-6 mbar. Cleaved
samples were exposed to ambient conditions for less than one minute.

###  2.2 Density functional theory calculations

Density functional theory (DFT) was used to calculate the electronic structure
of substitutional tin defects in monolayer and multilayer BP. The calculations
were performed using the CP2K code 19, which employs a mixed Gaussian and
plane wave basis-set (GPW). The DZ_MOLOPT_GTH basis set was used together with
Goedecker-Teter-Hutter (GTH) pseudopotentials 20. The plane wave cutoff was
converged at 400 Ry, SCF convergence was set to $10^{-6}$ a.u. and residual
forces on atoms were smaller than 0.01 eV/Å. Since GGA functionals predict a
metallic behavior for BP, the PBE0-TC-LRC hybrid functional was used with a
cutoff radius of 2 Å and 10% Hartree-Fock (HF) exchange 21. In order to reduce
the computational cost of the hybrid functional calculations, the auxiliary
density matrix method (ADMM) was used, which uses a reduced basis set for the
HF exchange calculation and thus enables the study of supercells up to 1500
atoms 22

The optimized unit cell vectors were found to be $a=3.4$ Å, $b=4.5$ Å, and
$c=10.9$ Å which is in good agreement with experiment and other DFT studies.
The bandgap for bulk BP was calculated to be 0.5 eV, which increased to 1.35
eV for a monolayer. The surface properties of a thick sample used in this work
were calculated using a 1296 atom cell ($9\times 9\times 4$ supercell) and the
monolayer consisted of 324 atoms ($9\times 9$) up to 1296 atoms ($18\times
18$) cells with 20 Å vacuum gap. In order to estimate the band offset with
electrodes, the ionization potential of the slab was calculated with respect
to the vacuum level. For a monolayer, the VBM is at $-5.5$ eV, but increases
to -5.2 eV for a bilayer and $-5.1$ eV for a four layer slab.

To calculate the defect formation energies and charge transition levels, the
chemical potentials of P and Sn atoms in the gas phase were used. All
calculations included the potential alignment; however, we did not include
charge corrections when calculating the charged states of Sn impurity but
rather relied on the fact that these states are strongly delocalized and Sn
atoms are positioned far away from each other due to large supercells. For a
monolayer, we could achieve the formation energy convergence with the cell
size.

###  2.3 Tight binding method

Tight-binding calculations were carried out to study the electronic structure
of BP with a charged Sn substitutional defect. The BP surface is modelled as a
monolayer that sits on a substrate with the dielectric constant of BP. For the
BP monolayer, tight-binding parameters of Rudenko and Katsnelson23 were used.
As the band gap of the monolayer is significantly larger than what is measured
on the BP surface, we first diagonalize the Hamiltonian for the defect-free
system and then apply a ”scissor” shift to correct the band gap. Next, the
defect potential is expressed in the basis of scissor-corrected eigenstates
and the Hamiltonian is diagonalized again. To model the charged defect, the
screened Coulomb potential induced by a point charge located at the position
of the Sn nucleus (in the relaxed DFT geometry) is added as an on-site term to
the Hamiltonian. For the screened potential, we employ an anisotropic Keldysh
model of the form $V_{\text{K}}(x,y)=\frac{Z\pi
e^{2}}{2r_{0}}[H_{0}(\epsilon_{\text{env}}\frac{r(x,y)}{r_{0}})-Y_{0}(\frac{\epsilon_{\text{env}}r(x,y)}{r_{0}})]$,
where Z=-1 denotes the defect charge, $r_{0}=22.25$ Å12 and
$\epsilon_{\text{env}}=10.5$ denotes the environment dielectric constant
(determined by averaging the dielectric constants of bulk BP and of vacuum).
Moreover, $r(x,y)=\sqrt{(\lambda_{x}x)^{2}+(\lambda_{y}y)^{2}+d^{2}}$
describes the anisotropy of the screening with $\lambda_{x}=0.42$,
$\lambda_{y}=1.2$ and $d=0.35$ Å. The values of $d$ (which denotes the height
of the Sn nucleus above the BP sheet), $\lambda_{x}$ and $\lambda_{y}$ were
fitted to obtain good agreement with the measured tunneling spectrum of the
defect. To obtain converged results for a single charged defect, we employ a
$60\times 60$ supercell and use a $3\times 3$ k-point sampling of the
corresponding first Brillouin zone.

##  3 Results

###  3.1 STM imaging and chemical composition of BP samples

Bulk BP crystals were cleaved under ultrahigh vacuum and imaged with STM at 77
K (see Methods). An atomic-resolution STM image showing the zig-zag rows of
the topmost atoms of a BP crystal is shown in Fig. 1b. Fourier analysis of
these images (not shown) yields lattice constants of $3.5\pm 0.1$ Å and
$4.2\pm 0.1$ Å in the zigzag and armchair directions, respectively, in
agreement with our density-functional theory calculations and previous studies
7, 8, 17. A larger area image where 10 defects can be observed is shown in
Fig. 1c. These defects are typical of those in the literature 7, 8, 13 and
present double-lobed protrusions that can extend as far as 10 nm in their
long-axis direction (aligned to the armchair direction of the crystal). A
higher resolution image showing a surface region where two defects can be seen
is shown in Fig. 1d.

Figure 1: (a) Structural model of a single phosphorene sheet shown in top and
side views. The phosphorus atoms at the top and bottom of the sheet are
colored red and black, respectively; typical STM images resolve only the top
(red colored) atoms. Also shown is a single unit cell of BP, which shows two
adjacent phosphorene layers in the AB stacking configuration. (b) Atomic-
resolution STM image of a clean and defect free area of a BP surface (image
parameters: $-0.4$ V, 40 pA). (c) Large area STM image showing multiple ($\sim
10$) defects that appear as long range elongated protrusions ($0.4$ V, 40 pA).
(d) Atomic-resolution image showing a pair of defects ($-0.4$ V, 40 pA). The
zig-zag rows are aligned with the vertical axis of the image in panels (b) and
(d), but run diagonally in panel (c). Figure 2: XPS analysis of vapour-grown
(VG) and high-pressure (HP) grown black phosphorus. Left: survey spectra
resolving P and C peaks. The broad peak at 1230eV is the C KLL Auger peak.
Right: high-resolution spectra of present impurities. O 1s peaks are observed
at 532.5 and 530.2eV, the Sn 3d doublet at 487.0 and 495.2eV, and the I 3d
doublet at 619.2 and 630.1eV

We use commercial BP samples obtained from two different suppliers: samples
obtained from Smart-elements GmbH were fabricated using a vapour transport
method 24, 25, while samples obtained from 2D Semiconductors Inc. were
fabricated via a high pressure conversion of red phosphorus 26, 27. We refer
to these samples as “vapour grown” and “high pressure” BP, respectively. To
test for the presence of chemical impurities, we have analysed freshly cleaved
samples of both sample types using XPS. Figures 2a,b show XPS survey scans for
both sample types. All the major peaks in this scan belong to phosphorus with
the only other feature being trace carbon contamination, likely due to low-
level hydrocarbon background in the high-vacuum chamber used for the
measurement. Additional contamination peaks (O, Sn, and I) were found using
higher resolution and sensitivity XPS measurements, and we present these
results in Figs. 2c–h. The presence of an oxygen peak at 532.5 eV for both
samples (Figs. 2c,d) can be attributed to physisorbed water, since the samples
were cleaved in ambient conditions before being immediately introduced to the
vacuum chamber. The vapor grown sample (Fig. 2c) exhibits a smaller component
at 530.2 eV, which may indicate some trace oxidation of the sample 28; this
difference between the two samples can be attributed to the significantly
smaller crystal sizes of vapour grown BP, producing a higher density of
crystal edges and rendering the samples more vulnerable to oxidisation 29. Tin
and iodine peaks were also found in both sample types (Figs. 2e–h). The
observation of these elements in our samples is in agreement with a previous
report 30 and can be attributed to the use of tin and iodine as catalysts
during synthesis. Similar high resolution scans in other regions (e.g., Fe,
Cu, Ti, Zn) did not produce any measurable peaks. It is notable that for both
of our sample types the tin and iodine peaks are close to the measurement
limit of our instrument.

###  3.2 DFT calculations of tin and iodine impurities in BP

The observation of tin and iodine in our samples suggests the possibility that
one or both of these impurity species is responsible for the p-type doping of
BP and/or the observed long-range defect features seen in STM images. Recent
DFT calculations examined a wide range of impurities in BP in adsorbed,
intercalated, and substitutional lattice positions 17. Iodine was found to
produce the least stable structures of all the elements considered. The most
stable iodine structure was an adsorbed neutral iodine atom, which did not
produce any electronic states in the BP bandgap 17. The absence of a localised
charge or bandgap state rules out iodine being responsible for the features
observed in STM images. In contrast, the substitutional tin impurity was found
to be the most stable of all impurity defects considered 17. To better
understand the defect behaviour as it is moved from the surface into the bulk,
Sn impurities were calculated in various depths in a four layer BP slab. The
slab consisted of 1296 atoms ($9\times 9$ surface unit cell) and a single
substitutional Sn was calculated in the upper row of the top layer, in the
bottom row of the top layer and similarly in the second layer. As the Sn atom
is moved into the slab the formation energy decreases slightly. The positively
charged Sn in the top row of the top layer has the formation energy of -3.50
eV, while that in the second layer is slightly reduced to -3.56 eV.

To investigate whether substitutional tin impurities might be the cause of the
p-type doping in BP and the features observed in STM images, it is necessary
to determine the relative stability of the three possible charge states of the
defect (negative, neutral, and positive) in the presence of an electron
reservoir with a chemical potential equal to the Fermi level of the crystal.
To achieve this, we calculate the formation energies of the positive and
negative charge states of the substitutional tin impurity at the surface of
the four layer slab by adding or subtracting a single electron from the
computational cell 31. The resulting defect states are drawn schematically in
Fig. 3a. Each phosphorus atom in BP is three-fold coordinated with an electron
lone pair. Tin has one less valence electron than phosphorus, such that the
neutral state of tin in BP features a singly-occupied dangling bond orbital.
The negative charge state of the substitutional defect has one additional
electron, forming a lone pair, while the positive charge state has one less
electron than the neutral case, leaving an unoccupied dangling bond orbital on
the tin atom. The relative stability of the three charge states as a function
of the Fermi level, referenced to the valence band maximum, is ploted in Fig.
3(b). From this charge transition level plot, we see that the negative charge
state of the substitutional tin impurity is energetically favourable for Fermi
level values across most of the bandgap; the positively charged state is
favourable only when the Fermi level lies within 0.02 eV of the valence band
edge, and the neutral state is never favoured. This demonstrates that an
isolated substitutional tin impurity in BP will naturally adopt a negatively
charged configuration by capturing an electron from the BP conduction band. As
the number of substitutional tin impurities increases, this will naturally
lead to a shift of the Fermi level toward the valence band in an analogous
manner to acceptor states in cubic semiconductors. The above results suggest
that the p-type doping commonly observed in BP samples can be attributed to
substitutional tin impurities.

Figure 3: DFT modelling of substitutional Sn in black phosphorus. (a):
Substitutional Sn adopts the same configuration as P, replacing the lone pair
with a dangling bond that can hold 0, 1 or 2 electrons. (b): Formation energy
of charged Sn defects for varying Fermi level. 0.0 eV corresponds to the
valence band maximum.

###  3.3 Tunneling spectroscopy of surface defects

Next, we consider whether substitutional tin impurities can explain the long-
ranged defect features commonly observed in STM images of BP. In cubic
semiconductors, such as silicon and gallium arsenide, long-ranged features
that appear superimposed on the surface atomic lattice have been attributed to
presence of bulk substitutional dopant impurities located within several
nanometres of the surface 32, 33, 34. Early reports found isotropic circular
protrusions and depressions that were attributed to the effect of a localised
charge-induced band bending, i.e., imaging the screened Coulomb potential of
the dopant atom 32. Later reports discovered that in some cases it is also
possible to find highly anisotropic protrusions, and these features have been
attributed to directly imaging the probability distributions of the hydrogenic
states of substitutional dopants 33, 34. The highly anisotropic, double-lobed
protrusions that are seen in BP (Figs. 1b-d) are thus suggestive of imaging
anisotropic hydrogenic states induced by a substitutional dopant.

To determine whether this is the case, we have performed tunnelling
spectroscopy measurements on the defects. Figure 4a shows tunnelling
spectroscopy measurements recorded over a defect-free area of the surface
(black trace) and over a double-lobed feature (red trace). From the
measurement over the defect free region we are able to extract a bandgap of
$0.26\pm 0.2$ eV, in agreement with previous reports 35, and we confirm that
our sample is p-type 5, 6, 7, 8. Examining the STS measurement over the
defect, we notice a bandgap state peak at $\sim 0.02$ eV, just above the
valence band maximum. We also notice that the defect measurement shows a
shoulder at $\sim 0.12$ eV, suggestive of a second peak deep within the
bandgap. These features are consistent with previously published spectra of
defects in BP 7.

Figure 4: Comparison of experimental spectroscopy and tight-binding
calculations. (a) and (b): experimental and simulated dI/dV curves on the bare
surface (black) and in the vicinity of a double-lobed defect (red). (c) and
(d): Experimental and simulated spatially resolved STS along a line through a
double-lobed defect, where the colour scale corresponds to conductance.

Figure 4c shows tunnelling spectrum as a function of position (horizontal
axis) along the centre line of a single defect. The left and right regions
(near 0 nm and 10 nm) show the spectrum acquired over the clean surface and
are in agreement with the point spectroscopy measurements in Fig. 4a (black
trace). The central region of the measurement, i.e., between $\sim 2$ and 8
nm, shows the spatial variation of the tunnelling spectroscopy along the
length of the defect. In particular, the bandgap state that was seen at $\sim
0.02$ eV in the point spectrum (Fig. 4a) is resolved into a double-lobed
feature with a node in the centre extending $\sim 8$ nm across the surface;
this is the same spatial variation as the features observed in the topographic
images, demonstrating that this bandgap state is the origin of the double-
lobed features observed in STM images. The broad shoulder at $\sim 0.12$ eV
appears as a diffuse region of enhanced conductivity, notably appearing as a
continuous region with a maximum intensity near the centre rather than a node
as in the case for the 0.02 eV feature. We do also see some upward band
bending of the band edges on either side of the defect, consistent with the
defect carrying a negative point charge.

###  3.4 Calculations of hydrogenic states of tin at the BP surface

In order to explain the physical characteristics of the defect state and
determine whether it can be attributed to substitutional tin, it is necessary
to perform simulations with unit cells large enough to model hydrogenic wave
functions, and to capture the effects of charge-induced band bending. Such
unit cells consist of many thousands of atoms, and are not currently
accessible using DFT. Therefore, we have modelled the hydrogenic states
produced by a substitutional tin impurity using tight-binding calculations. In
these calculations, the defect is modelled as a negative point charge at the
position of the tin nucleus (see Methods), which represents the negative
charge that is localised on the lone pair orbital of the substutional tin
impurity as predicted by our DFT calculations. We plot the simulated local
density of states (LDOS) resulting from the tight-binding calculation in Fig.
4b; the red and black traces show the calculation with and without the
negative point charge, respectively, and these can be compared directly to the
corresponding experimentally measured tunnelling point spectra in Fig. 4a.

The simulated LDOS of the defect-free surface (black trace, Fig. 4b)
reproduces the bandgap seen in the experimental data (black trace, Fig. 4a).
The simulated LDOS of the defect (red trace, Fig. 4b) shows a bandgap state
near the valence band maximum and a second peak near 0.18 eV, in good
agreement with the experimental data (red trace, Fig. 4a). The spatial
variation of the LDOS across a defect is shown in Fig. 4d and can be compared
to the experimental measurement in Fig. 4c. Here we see that the double-lobed
state just inside the bandgap at the valence band edge is reproduced in our
tight bonding model, as is the upward band bending of the conduction band
edge. The tight binding simulation also shows a nodeless state close to mid-
gap, that corresponds well to the nodeless state observed in the experimental
data.

The gap states in the tight-binding simulations are hydrogenic states induced
by the localised negative point charge placed in the lattice to represent the
negatively-charged substitutional tin impurity. This negative charge creates a
partially screened Coulomb potential and induces hydrogenic defect states.
Analysis of the simulation data shows that the two bandgap states in Fig. 4b
can be attributed to the 1s hydrogenic ground state, and the 2px hydrogenic
excited state (where $x$ corresponds to the heavy hole, i.e., armchair,
direction), while the 2s and 2py states are shifted into the BP valence band.
It is the anisotropic electronic structure of BP that results in the
reordering of the state energies compared with the normal spherical hydrogenic
model.

##  4 Conclusion

We have investigated the origin of p-type doping in nominally-intrinsic BP,
and whether this doping can be attributed to the long-ranged defects commonly
observed in STM images of BP. Using the combination of XPS, atomic-resolution
STM/STS, and DFT, we demonstrated that tin impurities are present in BP
samples, and that substitutional tin impurities are both highly stable and act
as acceptor states. Using tight-binding calculations, we have demonstrated
that substitutional tin impurities induce hydrogenic defect states in BP, and
that the 2p state of these defects is the main source of the long-range
features observed in STM images of BP. These observations provide compelling
evidence that substitutional tin impurities are responsible for the p-type
doping in BP, and are the origin of the long-ranged defect features observed
in STM images of BP.

## Acknowledgements

MW and MA were supported through studentships in the EPSRC Centre for Doctoral
Training in Advanced Characterisation of Materials (EP/L015277/1) and EPSRC
Centre for Doctoral Training in Theory and Simulation of Materials
(EP/L015579/1), respectively. AAM, JL, ALS, MA, JG, and AJK acknowledge the
support of the Thomas Young Centre through grant TYC-101. ALS and AJK
acknowledge funding provided by EPSRC (EP/P013503/1), and Leverhulme Trust
(RPG-2016-135). SRS acknowledges funding provided by EPSRC (EP/L002140/1). JG
acknowledges the use of the ARCHER high-performance computing facilities via
the membership in the U.K.s HPC Materials Chemistry Consortium, funded by
EPSRC (EP/L000202) and the computer resources of the UK Materials and
Molecular Modelling Hub, which is partially funded by EPSRC (EP/P020194).

## References

  * Kim et al. 2019 Kim, D. K.; Hong, S. B.; Jeong, K.; Lee, C.; Kim, H.; Cho, M. H. P-N Junction Diode Using Plasma Boron-Doped Black Phosphorus for High-Performance Photovoltaic Devices. _ACS Nano_ 2019, _13_ , 1683–1693
  * Jia et al. 2019 Jia, J.; Jeon, S.; Jeon, J.; Park, J. H.; Lee, S. Versatile Doping Control of Black Phosphorus and Functional Junction Structures. _J. Phys. Chem. C_ 2019, _123_ , 10682–10688
  * Novoselov et al. 2004 Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. _Science_ 2004, _306_ , 666–669
  * Warschauer 1963 Warschauer, D. Electrical and optical properties of crystalline black phosphorus. _J. Appl. Phys._ 1963, _34_ , 1853–1860
  * Akahama et al. 1983 Akahama, Y.; Endo, S.; Narita, S.-i. Electrical Properties of Black Phosphorus Single Crystals. _Journal of the Physical Society of Japan_ 1983, _52_ , 2148–2155
  * Castellanos-Gomez 2015 Castellanos-Gomez, A. Black Phosphorus: Narrow Gap, Wide Applications. _Journal of Physical Chemistry Letters_ 2015, _6_ , 4280–4291
  * Kiraly et al. 2017 Kiraly, B.; Hauptmann, N.; Rudenko, A. N.; Katsnelson, M. I.; Khajetoorians, A. A. Probing Single Vacancies in Black Phosphorus at the Atomic Level. _Nano Letters_ 2017, _17_ , 3607–3612
  * Qiu et al. 2017 Qiu, Z.; Fang, H.; Carvalho, A.; Rodin, A. S.; Liu, Y.; Tan, S. J.; Telychko, M.; Lv, P.; Su, J.; Wang, Y. et al. Resolving the Spatial Structures of Bound Hole States in Black Phosphorus. _Nano Letters_ 2017, _17_ , 6935–6940
  * Li et al. 2012 Li, L.; Lu, W.; Zhu, X.; Sun, Y. Effect of Fe-doping on charge density wave in 1T-TaS2: Instability and induced superconductivity. _Journal of Physics: Conference Series_ 2012, _400_ , 022061
  * McChesney et al. 2010 McChesney, J. L.; Bostwick, A.; Ohta, T.; Seyller, T.; Horn, K.; González, J.; Rotenberg, E. Extended van Hove Singularity and Superconducting Instability in Doped Graphene. _Phys. Rev. Lett._ 2010, _104_ , 136803
  * Xia et al. 2014 Xia, F.; Wang, H.; Jia, Y. Rediscovering black phosphorus as an anisotropic layered material for optoelectronics and electronics. _Nat. Commun._ 2014, _5_ , 4458
  * Rodin et al. 2014 Rodin, A. S.; Carvalho, A.; Castro Neto, A. H. Strain-Induced Gap Modification in Black Phosphorus. _Phys. Rev. Lett._ 2014, _112_ , 176801
  * Zhang et al. 2009 Zhang, C. D.; Lian, J. C.; Yi, W.; Jiang, Y. H.; Liu, L. W.; Hu, H.; Xiao, W. D.; Du, S. X.; Sun, L. L.; Gao, H. J. Surface Structures of Black Phosphorus Investigated with Scanning Tunneling Microscopy. _The Journal of Physical Chemistry C_ 2009, _113_ , 18823–18826
  * Riffle et al. 2018 Riffle, J. V.; Flynn, C.; St Laurent, B.; Ayotte, C. A.; Caputo, C. A.; Hollen, S. M. Impact of vacancies on electronic properties of black phosphorus probed by STM. _Journal of Applied Physics_ 2018, _123_
  * Guo and Robertson 2015 Guo, Y.; Robertson, J. Vacancy and Doping States in Monolayer and bulk Black Phosphorus. _Sci.Rep_ 2015, 1–10
  * Cai et al. 2016 Cai, Y.; Ke, Q.; Zhang, G.; Yakobson, B. I.; Zhang, Y. W. Highly Itinerant Atomic Vacancies in Phosphorene. _JACS_ 2016, _138_ , 10199–10206
  * Gaberle and Shluger 2018 Gaberle, J.; Shluger, A. L. Structure and Properties of Intrinsic and Extrinsic Defects in Black Phosphorus. _Nanoscale_ 2018, _10_ , 19536–19546
  * Hu and Dong 2015 Hu, T.; Dong, J. Geometric and electronic structures of mono- and di-vacancies in phosphorene. _Nanotechnology_ 2015, _26_ , 0–9
  * VandeVondele et al. 2005 VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. _Comp. Phys. Commun._ 2005, _167_ , 103–128
  * VandeVondele and Hutter 2007 VandeVondele, J.; Hutter, J. Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases. _J. Chem. Phys._ 2007, _127_ , 114105–1–9
  * Guidon et al. 2009 Guidon, M.; Hutter, J.; VandeVondele, J. Robust Periodic Hartree Fock Exchange for Large-Scale Simulations Using Gaussian Basis Sets. _J. Chem. Theory Comput._ 2009, _5_ , 3010–3021
  * Guidon et al. 2010 Guidon, M.; Hutter, J.; VandeVondele, J. Auxiliary Density Matrix Methods for Hartree Fock Exchange Calculations. _Journal of Chemical Theory and Computation_ 2010, _6_ , 2348–2364
  * Rudenko and Katsnelson 2014 Rudenko, A. N.; Katsnelson, M. I. Quasiparticle band structure and tight-binding model for single- and bilayer black phosphorus. _Phys. Rev. B_ 2014, _89_ , 1–6
  * Lange et al. 2007 Lange, S.; Schmidt, P.; Nilges, T. Au3SnP7 Black Phosphorus: An Easy Access to Black Phosphorus. _Inorganic Chemistry_ 2007, _46_ , 4028–4035
  * Köpf et al. 2014 Köpf, M.; Eckstein, N.; Pfister, D.; Grotz, C.; Krüger, I.; Greiwe, M.; Hansen, T.; Kohlmann, H.; Nilges, T. Access and in situ growth of phosphorene-precursor black phosphorus. _Journal of Crystal Growth_ 2014, _405_ , 6–10
  * Bridgman 1914 Bridgman, P. W. Two New Modifications of Phosphorus. _JACS_ 1914, _36_ , 1344–1363
  * Shirotani 1982 Shirotani, I. Growth of Large Single Crystals of Black Phosphorus at High Pressures and Temperatures, and its Electrical Properties. _Molecular Crystals and Liquid Crystals_ 1982, _86_ , 203–211
  * Luo et al. 2016 Luo, W.; Zemlyanov, D. Y.; Milligan, C. A.; Du, Y.; Yang, L.; Wu, Y.; Ye, P. D. Surface chemistry of black phosphorus under a controlled oxidative environment. _Nanotechnology_ 2016, _27_ , 434002
  * Kuntz et al. 2017 Kuntz, K. L.; Wells, R. A.; Hu, J.; Yang, T.; Dong, B.; Guo, H.; Woomer, A. H.; Druffel, D. L.; Alabanza, A.; Tomanek, D. et al. Control of Surface and Edge Oxidation on Phosphorene. _ACS Applied Materials & Interfaces_ 2017, _9_ , 9126–9135
  * Mayorga-Martinez et al. 2018 Mayorga-Martinez, C. C.; Sofer, Z.; Sedmidubský, D.; Luxa, J.; Kherzi, B.; Pumera, M. Metallic impurities in black phosphorus nanoflakes prepared by different synthetic routes. _Nanoscale_ 2018, _10_ , 1540–1546
  * Schofield et al. 2013 Schofield, S. R.; Studer, P.; Hirjibehedin, C. F.; Curson, N. J.; Aeppli, G.; Bowler, D. R. Quantum engineering at the silicon surface using dangling bonds. _Nat. Commun._ 2013, _4_ , 1649
  * Liu et al. 2001 Liu, L.; Yu, J.; Lyding, J. W. Atom-resolved three-dimensional mapping of boron dopants in Si(100) by scanning tunneling microscopy. _Appl. Phys. Lett._ 2001, _78_ , 386
  * Yakunin et al. 2004 Yakunin, a.; Silov, A.; Koenraad, P. M.; Wolter, J.; Van Roy, W.; De Boeck, J.; Tang, J.-M.; Flatté, M. Spatial Structure of an Individual Mn Acceptor in GaAs. _Phys. Rev. Lett._ 2004, _92_ , 216806
  * Sinthiptharakoon et al. 2014 Sinthiptharakoon, K.; Schofield, S. R.; Studer, P.; Brázdová, V.; Hirjibehedin, C. F.; Bowler, D. R.; Curson, N. J. Investigating individual arsenic dopant atoms in silicon using low-temperature scanning tunnelling microscopy. _J. Phys. Condens. Matter_ 2014, _26_ , 012001
  * Prytz and Flage-Larsen 2010 Prytz, O.; Flage-Larsen, E. The Influence of Exact Exchange Corrections in van der Waals Layered Narrow Bandgap Black Phosphorus. _J. Phys: Condens. Matter_ 2010, _22_ , 015502