**Superconducting- and Graphene-Based Devices**

Editor **Filippo Giubileo**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Filippo Giubileo CNR-SPIN Salerno via Giovanni Paolo II 132 Italy

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Nanomaterials* (ISSN 2079-4991) (available at: https://www.mdpi.com/journal/nanomaterials/ special issues/graphene nano).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4761-9 (Hbk) ISBN 978-3-0365-4762-6 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editor**

#### **Filippo Giubileo**

Filippo Giubileo is a Senior Researcher working at the Institute for Superconductors, Innovative Materials and Devices (SPIN) of the Italian National Research Council (CNR) since 2004. He has a National Scientific Qualification as Full Professor in Experimental Condensed Matter Physics. He has been a member of the scientific evaluation panel for the Flag-Era Joint Transnational call on graphene basic research. His research activity deals with electrical properties of 1D and 2D materials and their application as field effect transistors, photodetectors, chemical sensors, and field emitters. His expertise includes advanced surface analysis on the nanometer scale via scanning probe microscopy/spectroscopy in high magnetic fields and variable temperatures, fabrication and characterization of superconducting heterostructures, and pairing symmetry and electronic properties in superconductors. FG is listed in the World Top 2% Scientist classification released by PlosBiology in 2021. He has published more than 100 papers in peer-reviewed journals, and he was Chair-Organizer of 6 international conferences (ISMN08, SM-2010, S4E-2014, GM-2016, TTN-2018, NANO M&D-2019). He also serves as a member of editorial board for several international journals (Nanomaterials, Nano Express, Journal of Nanomaterials, Chemosensors).

### *Editorial* **Superconducting- and Graphene-Based Devices**

**Filippo Giubileo**

CNR-SPIN Salerno, 84084 Fisciano, Italy; filippo.giubileo@spin.cnr.it

This Special Issue has been organized to collect new or improved ideas regarding the exploitation of superconducting materials, as well as graphene, aiming to develop innovative devices. For instance, several graphene applications can be enhanced by modifying their surface to introduce a non-zero bandgap, tune adhesion and/or hydrophobicity/hydrophilicity, etc. Similarly, contact resistance in graphene-based field effect transistors can be improved by irradiation [1], leading to an improved device performance. Wei Qin et al. [2] conducted a detailed theoretical investigation using first-principles calculations of "Lithium Diffusion in Silicon Encapsulated with Graphene". They considered monolayer graphene on silicon substrate to simulate the Si microparticles that were encapsulated in a graphene cage, which can be exploited as anodes in lithium-ion batteries. They demonstrated that defective graphene strongly reduces the energy barriers for Li diffusion in Gr or Gr/Si. Abid et al.'s [3] report, entitled "Interface Kinetics Assisted Barrier Removal in Large Area 2D-WS2 Growth to Facilitate Mass Scale Device Production", employed chemical vapor deposition technique to synthesize mono- and few-layer WS2 with areas up to cm2 on graphene-oxide-coated Si/SiO2 substrates. They show that as-developed WS2 layers are polycrystalline (mono- and few-layer), with single-crystal domains that are triangular and hexagonal in shape. Alejandro Toral-Lopez et al.'s [4] report, on "GFET Asymmetric Transfer Response Analysis through Access Region Resistances", aimed to exploit graphene-based devices to increase the functionality of Si-technology in the field of radio-frequency electronics. They conducted an in-depth investigation of the role of access regions on the performance of graphene-based field effect transistors (GFETs). They demonstrated that the access region conductivity can be tuned by the back-gate bias, improving the RF performance. Graphene represents a prototype of 2D materials and is still widely investigated. Many layered materials, such as the transition metal dichalcogenides, are largely studied for their use as a conducting channel in nanometric field effect transistors, including MoS2 [5–7], ReSe2 [8], WSe2 [9], etc. Regarding superconducting-based devices, Jose C. Verde et al.'s [10] report is entitled "Calculations of Some Doping Nanostructurations and Patterns Improving the Functionality of High-Temperature Superconductors for Bolometer Device Applications". They propose that high-temperature superconductors (HTS) can be nanostructured (and patterned) to obtain an increased functionality as sensing materials for resistive transition-edge bolometer devices (TES). Calculations have been performed to consider the spatial variations in carrier doping into the CuO2 planes of the YBaCuO perovskite superconductor, demonstrating an improvement in the bolometric parameters with respect to conventional, nonstructured HTS materials. Paola Romano et al. [11] report, "Transport and Point Contact Measurements on Pr1−xCexPt4Ge12 Superconducting Polycrystals" demonstrated that the material has a collective pinning regime with a quasi-2D character for a Ce-doping of about x = 0.07. Moreover, while investigating the properties of metal/superconductor nano-junctions, they showed that the observed conductance features are explained in terms of a superconducting-order parameter with nodal directions, as well as a sign change in the momentum space. Indeed, numerical simulations reported in the framework of Blonder–Tinkham–Klapwijk model show that s-wave pairing and anisotropic s-wave are unsuitable for the reproduction of experimental data obtained at a low temperature. Carlo Barone et al.'s [12] report is entitled "Current-Resistance Effects

**Citation:** Giubileo, F.

Superconducting- and Graphene-Based Devices. *Nanomaterials* **2022**, *12*, 2055. https://doi.org/10.3390/ nano12122055

Received: 7 June 2022 Accepted: 10 June 2022 Published: 15 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Inducing Nonlinear Fluctuation Mechanisms in Granular Aluminum Oxide Nanowires". They measured electric transport and voltage fluctuations in the normal state and in the temperature range of 8–300 K, observing both nonlinear resistivity and two-level tunneling fluctuators. This study helps to improve the fabrication process, therefore reducing the possible sources of decoherence in the superconducting state. This is crucial in quantum technology applications. Sergio Pagano et al.'s [13] report is entitled "Iron-Based Superconducting Nanowires: Electric Transport and Voltage-Noise Properties". In this work, they fabricated ultra-thin Co-doped BaFe2As2 nanowires and characterized their transport and intrinsic noise properties. They also investigated the ageing effect on device degradation by means of noise spectroscopy. Interestingly, iron-based superconducting nanowire detectors have several advantages, due to their high operating temperature, when used as innovative single-photon detectors working in the visible and infrared spectral region.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Lithium Diffusion in Silicon Encapsulated with Graphene**

**Wei Qin 1,2,\*, Wen-Cai Lu 1,2,3,\*, Xu-Yan Xue 1,2, Kai-Ming Ho <sup>4</sup> and Cai-Zhuang Wang 4,\***


**Abstract:** The model of a graphene (Gr) sheet putting on a silicon (Si) substrate is used to simulate the structures of Si microparticles wrapped up in a graphene cage, which may be the anode of lithiumion batteries (LIBS) to improve the high-volume expansion of Si anode materials. The common low-energy defective graphene (*d*–Gr) structures of DV5–8–5, DV555–777 and SV are studied and compared with perfect graphene (*p*–Gr). First-principles calculations are performed to confirm the stable structures before and after Li penetrating through the Gr sheet or graphene/Si-substrate (Gr/Si) slab. The climbing image nudged elastic band (CI-NEB) method is performed to evaluate the diffusion barrier and seek the saddle point. The calculation results reveal that the *d*–Gr greatly reduces the energy barriers for Li diffusion in Gr or Gr/Si. The energy stability, structural configuration, bond length between the atoms and layer distances of these structures are also discussed in detail.

**Keywords:** perfect graphene (*p*–Gr); defective graphene (*d*–Gr); Gr/Si slab; diffusion barrier; CI-NEB calculation

#### **1. Introduction**

Facing the ever-growing demands for electrical vehicles and portable electronic devices [1–3], rechargeable lithium-ion batteries (LIBS) with high energy density, long cycle life and fast charge rate have become the focus of intense research. Silicon is an attractive and promising anode material for LIBS due to its high theoretical capacity (~4200 mAh/g) [4], over ten times higher than conventional graphite, and abundance on Earth. However, the large-volume expansion of silicon as an anode material upon lithiation (~300%) [5,6] in practical applications leads to fracture and loss of inter-particle electrical contact, consequently causing early capacity fading, thus blocking the further improvement of LIBS. Attempting to avoid the mechanical fracture via decreasing the material size, the use of nanostructure silicon (nano–Si) has been shown to be fruitful, such as in silicon crystalline amorphous core-shell nanowires [7], interconnected Si hollow nanosphere [8], Si nanotubes [9], and so on. In particular, Si–C nanostructure materials combining Si with carbonaceous materials [10–24] can retain capacity well, further enhancing the cycling stability and charging rate [10,25] due to the buffering effect and high electrical conductivity of carbon. On the other hand, a large number of theoretical calculations have been focused on the interactions between Li and Si crystals or between Li and graphite separately [26–34], or on the influence of Si–C composites on lithiation [35]. Despite the impressive improvement of the stability of Si anodes after adopting nano-Si, the Coulombic efficiency is still far lower than conventional graphite due to the large surface area. In addition, the complex synthesis process makes nano-Si costly. These are the major obstacles to the further development and mass production of LIBS technology. Subsequently, Si microparticles were explored to replace nano-Si as an anode material due to the lower cost in manufacture, but the cycling performance is far poorer than that of nano-Si. After the first few cycles, Si microparticles

**Citation:** Qin, W.; Lu, W.-C.; Xue, X.-Y.; Ho, K.-M.; Wang, C.-Z. Lithium Diffusion in Silicon Encapsulated with Graphene. *Nanomaterials* **2021**, *11*, 3397. https://doi.org/10.3390/ nano11123397

Academic Editors: Henrich Frielinghaus and Francisco Javier García Ruiz

Received: 12 November 2021 Accepted: 10 December 2021 Published: 15 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

are pulverized into smaller sizes, turn into amorphous particles from Si crystalline, and then lose electrical contact [36–39].

In recent years, a strategy of Si microparticles surrounded by a porous carbon matrix was used to counteract these inadequacies. The combination maintains both Li absorption capacity and electrical conductivity well for many battery-charging cycles. For example, a conformal growth of the conductive graphene cage, which can wrap up Si microparticles, was reported as a promising encapsulation material [40]. The similar surface chemistry to graphite allows graphene cage to form a stable solid electrolyte interphase [41,42]. Being mechanically strong and flexible, graphene cages remain undamaged with Si microparticles fracturing and confine all the fractured Si pieces within the cages, thus maintaining the essential electrical contact between broken Si particles within [40]. Graphene-encapsulated Si microparticles thus provide a new architecture of materials for the development of LIBS technology. A theoretical understanding of Li diffusion in such materials could therefore provide information and guidance for the design and development of LIBS, but the relevant calculated studies are scarce. Odbadrakh et al. reported that cavities along the reconstructed Si surface provide diffusion paths for Li [43]. Chou and Hwang investigated the role of the interface in lithiation of silicon/graphene (Si/Gr) [44] composites and demonstrated the charge transfer from Li to both silicon and graphene. They also showed that Li cations exhibit substantially higher mobility along the Si/Gr interface than that of bulk Si [44]. For Li diffusion, most studies focused on diffusion paths of Li on the graphene surface [45–47]. However, the difficulty of Li intercalating through the graphene sheet to bind to Si particles, i.e., the diffusion barrier during lithiation, is deeply influenced by the surface characteries of graphene for such a combination. It is well known that many defects will inevitably occur in the process of graphene producing. Can these defects be used to facilitate the diffusion of Li? If so, which common defects are easier for Li to intercalate through to reduce the diffusion barrier of Li? Do Si particles affect the diffusion? Does the diffusion barrier increase or decrease in the presence of Si particles? Studying these issues could assist in exploring more suitable encapsulation materials.

Usually, single–vacancy (SV), Stone–Wales (SW) and divacancy (DV) defects are typical point defects in graphene and the subjects of intense research [48–54]. SV and SW defects are created by removing a C atom and rotating a C–C bond by 90◦ in the plane, respectively. The DV defect is formed by removing a C–C dimer; when the dimer is removed from graphene, a series of reconstructions around the DV defect take place upon annealing. In this work, the Gr/Si (a monolayer graphene is put on the silicon substrate) is used to simulate the Si microparticles encapsulated with the graphene cage. The Si substrate and defective graphene (*d*–Gr) are adopted to simulate Si particles and porous graphene, respectively. The perfect graphene (*p*–Gr) is also studied as a comparison. The energy stability of these models (both Gr and Gr/Si, Gr = *p*–Gr, DV5–8–5, DV555–777 and SV) are investigated by first-principles calculations. Their energy barriers of Li diffusion through them are also determined.

#### **2. Computational Methods**

Considering both the calculation cost and accuracy, a 6 × 6 hexagonal supercell of *d*–Gr is employed to model the porous graphene cage. Three common defects [55], including the SV defect, DV5–8–5 defect (containing two pentagons and one octagon) and DV555–777 defect (three heptagons in the center surrounded by three pentagons) are considered. A 4 × 4 supercell of the Si (111) surface containing two double atomic layers are employed as the Si substrate. The silicon dangling bonds on the bottom are passivated by hydrogen atoms, and these hydrogen atoms and silicon atoms on the bottom are fixed. The Si–Si bond lengths in Si substrate are adjusted by about 4% to fit the 6 × 6 Gr sheet to minimize the effect of the lattice mismatch on graphene when the in-plane periodic boundary conditions are used. To eliminate the interaction between Gr (or Gr/Si) and their periodic replicas, the vacuum layers of 20 Å and 18 Å are adopted in the perpendicular direction to the Gr surface for monolayers Gr and Gr/Si, respectively. The structures of Gr and Gr/Si are displayed

in Figure 1. All calculations are performed using the projector-augmented wave (PAW) with Perdew, Burke, and Ernzerh of (PBE) GGA functional for the exchange-correlation energy [56,57], as implemented in the Vienna ab initio simulation package (VASP) [58]. We also consider the Van der Waal's force followed by the well-established Grimme's DFT-D2 formula [59,60] in the VASP calculation for all the Gr/Si. The Brillouin zone (BZ) integration is approximated by using the Gamma centered k-grid. For structural relaxation and energy evaluation, the 8 × 8 × 1 k-grid and 4 × 4 × 1 k-grid are adopted for the monolayer Gr and Gr/Si slab, respectively. We also calculate the energy barriers of Li diffusion in Gr and Gr/Si, and the climbing image nudged elastic band (CI-NEB) method implemented in the VASP code is used to seek the saddle points and the minimum energy path. As a reference, we study the adsorption sites about Li on the surfaces of perfect graphene (*p*–Gr) and *p*–Gr/Si. Three typical adsorption sites are considered: the top (T) site (Li is above a C atom), the bridge (B) site (Li is above the midpoint of a C–C bond) and the hollow (H) site (Li is above the center of a hexagon ring). The results show that the H site is the most stable adsorption site. Similarly, in *d*–Gr and *d*–Gr/Si (Gr = DV5–8–5, DV555–777 and SV), the structures are more stable when Li is above the center of the polygon rings. These lowest-energy structures are shown in Figures 2–4 as the initial configurations of the CI-NEB calculations. We call them Li/Gr and Li/Gr/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), respectively. When Li penetrates through the Gr sheet, Li is basically on the symmetric sites of the other surface of the graphene due to the planar structures of the monolayer Gr sheet. Their configurations are similar to the initial points from the top view, while for Gr/Si, Li is between the Gr sheet and Si substrate after passing through the Gr layer. We consider various relative sites between Li and Si atoms, such as Li above a Si atom, Li above a Si–Si bond, Li above a folded Si hexagon ring, and so on. All of the candidates are relaxed to find the lowest-energy structures. These lowest-energy structures after Li intercalates into Gr or Gr/Si are adopted as the end points and called Gr/Li and Gr/Li/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), respectively. In the CI-NEB calculation, four images are employed in addition to the initial and end points. The structures of the saddle points are found and called Li@Gr and Li@Gr/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), respectively. They are also displayed in Figures 2–4.

In order to evaluate the energetic stability of different *d*–Gr and *d*–Gr/Si (Gr = DV5–8–5, DV555–777 and SV), we firstly calculate their formation energies (*E*f). For monolayer Gr sheets, their *E*<sup>f</sup> is defined as

$$
\dot{E}\_{\rm f} = \dot{E}\_{\rm d} - N \,\mu\_{\rm C} \tag{1}
$$

and for *d*–Gr/Si slab, the *E*<sup>f</sup> is defined as

$$E\_\mathbf{f} = E\_\mathbf{d} - E\_{\text{subb}} - N\,\mu\_\mathbb{C} \tag{2}$$

where *Ed* and *E*sub are the total energy of *d*–Gr (or *d*–Gr/Si) and Si substrate in the given supercells, respectively; *N* is the number of C atoms in the supercell, and *μ*<sup>C</sup> is the chemical potential of carbon, which is leveled from perfect graphene with the same unit cell size. Our calculation results showed that the formation energies (*E*f) of DV5–8–5, DV555–777 and SV are 7.645, 6.386 and 7.579 eV, respectively, which is in agreement with previous studies [51,59]. The formation energies (*E*f) of DV5–8–5/Si, DV555–777/Si and SV/Si are 2.621, 1.254 and 1.906 eV, respectively. The order of *E*<sup>f</sup> does not alter with the existence of the Si substrate.

To examine the stability of the structures upon Li adsorption, we also calculate the adsorption energies (*E*ads) of Li/Gr and Li/Gr/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV) which are defined as

$$E\_{\text{ads}} = E\_{\text{t}} \text{ (Li-S)} - E\_{\text{t}} \text{ (S)} - E\_{\text{Li}} \tag{3}$$

where *E*t (Li–S) and *E*t (S) are the total energies of Gr (or Gr/Si) with and without a Li atom adsorbing on their surfaces, respectively. *E*Li is the energy of putting one Li atom into the same unit cell size. The results of the Li adsorption energies (*E*ads) are listed in the second column of Table 1. The lower the *E*ads, the more stable the adsorption structure (Li/Gr and Li/Gr/Si). From Table 1, we can see that for monolayer Gr (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), Li prefers to adsorb on the *d*–Gr (Gr = DV5–8–5, DV555–777 and SV), compared with *p*–Gr. The *E*ads of *d*–Gr is more than 1 eV lower than that of *p*–Gr. In particular, Li/SV is the most stable adsorption structure among the four sheets, whose *E*ads is nearly 1.8 eV lower than that of *p*–Gr. For Li/Gr/Si, the order of their *E*ads is consistent with that of Li/Gr, except Li/SV/Si. Comparing with *p*–Gr/Si, Li prefers to adsorb on DV5–8–5/Si and DV555–777/Si. The differences of their *E*ads are not as large as that between monolayer *p*–Gr and *d*–Gr, about 0.5 eV. It indicates that the differences of *E*ads are reduced between *p*–Gr/Si and *d*–Gr/Si (*d*–Gr = DV5–8–5 and DV555–777) due to the Si substrate. However, it is difficult to adsorb on SV/Si for Li. The *E*ads of Li/SV/Si is even higher, about 1 eV, than that of Li/*p*–Gr/Si. It may be related to the structural deformation of SV/Si when Li adsorbs on it (detailed in Results and Discussion).

**Figure 1.** Structures of (**a**) *p*–Gr, (**b**) DV5–8–5, (**c**) DV555–777, (**d**) SV, (**e**) *p*–Gr/Si, (**f**) DV5–8–5/Si, (**g**) DV555–777/Si and (**h**) SV/Si; the left and right configurations (**e**–**h**) are the top and side views of Gr/Si, respectively.

The formation energies (*E*f) of Gr/Li/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV) are also calculated when Li intercalates the Gr sheet and is located between Gr and the Si substrate, which is defined as

$$E\_{\rm f} = E\_{\rm t} \left( \text{Gr/Li/Si} \right) - E\_{\rm t} \left( \text{Gr} \right) - E\_{\rm t} \left( \text{Si} \right) - E\_{\rm Li} \tag{4}$$

where *E*t (Gr/Li/Si) is the total energies of the Gr/Li/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV) intercalated structures; *E*<sup>t</sup> (Gr) and *E*<sup>t</sup> (Si) are the total energies of the Gr sheet and Si substrate, respectively. In order to measure the energy stability of the structures before and after Li intercalation, we also calculate their energy differences (*E*diff), which are defined as

$$E\_{\rm diff} = E\_{\rm f} \left( \rm Li/Gr/Si \right) - E\_{\rm f} \left( \rm Gr/Li/Si \right) \tag{5}$$

where *E*<sup>f</sup> (Li/Gr/Si) and *E*<sup>f</sup> (Gr/Li/Si) are the formation energies of Li/Gr/Si and Gr/Li/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV) calculated in the same method (Formula (4)). The values of *E*<sup>f</sup> of Gr/Li/Si and *E*diff are shown in the third and fourth columns of Table 1, respectively. The results show that the structures after Li passes through the Gr sheet are more stable than those of the adsorption on the Gr surface. It may be because Li combines with both C atoms in the Gr and Si atoms on the surface of the Si substrate when Li passes through the Gr sheet, which make the structure more stable. The difference value of *p*– Gr/Si is the largest. It indicates that the stability of the intercalated structure (*p*–Gr/Li/Si) after Li passes through *p*–Gr is much higher than that of the adsorption on the Gr surface (Li/*p*–Gr/Si). By comparison, the *E*diff value for *d*–Gr/Si is relatively low.

**Figure 2.** Structures of (**a**) Li/*p*–Gr (**b**) Li@*p*–Gr (**c**) Li/DV5–8–5 (**d**) Li@DV5–8–5 (**e**) Li/DV555–777 and (**f**) Li@DV555–777; their main bond lengths (abbreviated "b" in the list, Å) are shown in the list at the bottom of the figure. The red numbers in these figures correspond to the subscripts in the list.

**Figure 3.** Structures of (**a**) Li/*p*–Gr/Si (**b**) Li@*p*–Gr/Si (**c**) *p*–Gr/Li/Si (**d**) Li/DV5–8–5/Si (**e**) Li@DV5–8–5/Si (**f**) DV5–8– 5/Li/Si (**g**) Li/DV555–777/Si (**h**) Li@DV555–777/Si and (**i**) DV555–777/Li/Si. The left and right configurations (**c**,**f**,**i**) are the top and side views of the end structures in the CI-NEB calculation, respectively. The main bond lengths (abbreviated as "b" in the list, Å) are shown in the list at the bottom of the figure. The red numbers in these figures correspond to the subscripts.

**Figure 4.** Structures of (**a**) Li/SV (**b**) Li@SV (**c**) SV/Li (**d**) Li/SV/Si (**e**) Li@SV/Si and (**f**) SV/Li/Si. The left and right configurations are the top and side views of the six structures in the CI-NEB calculation, respectively. The main bond lengths (abbreviated as "b" in the list, Å) are shown in the list at the bottom of the figure. The red numbers in these figures correspond to the subscripts.



#### **3. Results and Discussion**

#### *3.1. Energetic Stability*

Figures 2–4 show the most stable structures of Li adsorbed on different Gr and Gr/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV) surfaces. From the figures, we can see that Li prefers to absorb above the center of the polygon ring with the most sides so that it can have more carbon neighbors. For instance, the adsorption energies (*E*ads) are lower when Li absorbs above the hexagon ring center on the *p*–Gr and *p*–Gr/Si surfaces. Similarly, the lowest-energy adsorption sites are above the centers of the heptagon and octagon rings for DV555–777 (or DV555–777/Si) and DV5–8–5 (or DV5–8–5/Si) structures, respectively. For SV (or SV/Si), the most stable adsorption site is just above the C atom, removed from the SV defect. When Li is diffused into the Gr (or Gr/Si), Li is in the plane of Gr for the saddlepoint structures. After penetrating through the Gr sheet, Li is basically on the symmetric sites of the other surface of the Gr sheet, and in Gr/Si, the positions of Li are slightly shifted from the top view. The Gr (Gr = *p*–Gr, DV5–8–5 and DV555–777) layer maintains a planar structure throughout, whether there is the Si substrate or not (see Figures 2 and 3). The SV defect is a little special. When Li is adsorbed on the surface of SV, the C atoms in the SV are no longer in the same plane. One of the C atoms is obviously deviated from the SV plane and approaches the Li atom, whether in Li/SV or Li/SV/Si. After Li penetrates through SV, the SV layer is also non-planar, whether for SV/Li or SV/Li/Si. In SV/Li, one C atom is close to the Li atom; more C atoms deviate from the SV plane and approach the Li atom and Si surface in SV/Li/Si. For the saddle-point structures, the Li atom is basically in the initial plane of SV, but the C atoms around Li are almost deviated from the original plane. In Li@SV, these C atoms are distributed on both sides of the plane, while in Li@SV/Si, they are all close to the side of the Si surface (see Figure 4).

In the bottom of Figures 2–4, the relevant bond lengths, including b(C–C), b(Li–C) and b(Li–Si) are also listed. The C–C bond length (b1(C–C) in Figure 2a) of the hexagon ring absorbed by Li in Li/*p*–Gr is 1.428 Å. The corresponding value is 1.429 Å (b1(C–C) in Figure 3a) when there is the Si substrate, which is basically unchanged. Similar behavior is seen in DV5–8–5 and DV555–777 (see Figures 2 and 3). We also find the similar regular pattern in their saddle-point structures Li@Gr/Si (Gr = *p*–Gr, DV5–8–5 and DV555–777). The maximum variation of the C–C bond length is about 0.03 Å when there is the Si substrate. It indicates that the Si substrate has negligible effects on the C–C bond lengths. Furthermore, the C–C bond lengths in the saddle-point structures are larger than the ones of the corresponding initial structures Li/Gr or Li/Gr/Si (Gr = *p*–Gr, DV5–8–5 and DV555– 777). When Li penetration is complete, the C–C bond lengths are basically restored to the values of the initial structures. It may be due to Li being in the plane of Gr in the saddlepoint structures, while Li leaves the plane of Gr in Gr/Li or Gr/Li/Si (Gr = *p*–Gr, DV5–8–5 and DV555–777). Similarly, the Li–C bond lengths also vary very little when there is the Si substrate. Among the three Gr sheets (Gr = *p*–Gr, DV5–8–5 and DV555–777), the Li–C bond lengths' value of DV555–777 changes most. Despite all this, its average Li–C bond length changes less than 0.03 Å when there is the Si substrate, compared with the values b5–8(Li–C) in Li/DV555–777 (from 2.233 Å to 2.438 Å in Figure 2e) and Li/DV555–777/Si (from 2.249 Å to 2.425 Å in Figure 3g). Comparing with the three Gr (or Gr/Si) above, there are some different results due to the C atoms around Li deviating from the SV plane when Li diffuses in the SV (or SV/Si). The average bond lengths of C–C and Li–C in the saddlepoint structure Li@SV do not increase although the Li atom is in the SV plane but decrease. After Li penetrates the SV plane, both C–C and Li–C bond lengths completely return to the initial values. For SV/Si, the results are slightly complicated due to the existence of Si substrate. The average bond length of C–C in the saddle-point structure Li@SV/Si increases slightly (~0.006 Å) compared to the value of the initial structure (Li/SV/Si), while the average bond length of Li–C decreases over 0.3 Å compared to the initial structure. After Li penetrates the SV plane, the positions of the C atoms around Li change greatly, and the average bond length of C–C does not return to the initial value, but is increased by about 0.03 Å. Similar to the other three Gr/Si, both the C and Si atoms bond with the Li atom at this time. In addition, the shortest bond length of Li–C in Li/SV (b8(Li–C) in Figure 4a) is 2.049 Å, which is lower than the values in other Li/Gr (b2(Li–C) in Figure 2a, b4(Li–C) in Figure 2c and b6(Li–C) in Figure 2e). The closer combination between Li and C in Li/SV can reduce *E*ads. The order of *E*ads does not change when putting these Gr sheets on the Si substrate, except SV/Si. It is the most difficult to adsorb on the surface of SV/Si for the Li atom in these four slabs, which may be due to the structural deformation of SV/Si when Li adsorbs on it.

As a reference, we also calculate the layer distances of each structure, displayed in Table 2. We first calculate the average coordinates of all the C atoms in the Gr (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV) sheet and Si atoms on the surface of the Si substrate as the positions of the Gr sheet and Si surface, respectively. Then the distances between the Li atoms, Gr sheet and Si surface are calculated and labeled as *d*C–Si, *d*Li–C and *d*Li–Si in Table 2. The distances between the Gr sheet and the surface of Si substrate (*d*C–Si) in most Gr/Si slabs change little upon Li adsorption. Especially for DV555–777/Si, its *d*C–Si value is almost unchanged before and after Li adsorption. However, the *d*C–Si value of SV/Si decrease by about 0.1 Å upon Li adsorption, which is about 10 times that of other structural variation. It indicates that the adsorption of Li has the greatest effect on the bonding between SV and the surface of the Si substrate among the four Gr above, while the effect is the smallest in DV555–777/Si. In Gr/Li/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), Li is located between the Gr and the surface of the Si substrate as Li diffuses into the Gr. Although there is one more Li atom between them, compared with Gr/Si, most *d*C–Si values of Gr/Li/Si do not increase but decrease, except Li/DV555–777/Si. In sum, the adsorption of Li has a slight effect on the *d*C–Si values in *p*–Gr/Si and DV5–8–5/Si, and the effect disappears after Li diffuses into *p*–Gr/Si (or DV5–8–5/Si). For DV555–777/Si, when Li is adsorbed on the DV555–777 surface, the *d*C–Si value is almost unchanged. However, when Li diffuses into DV555–777, the distance between DV555–777 and the surface of Si substrate is widened. Nevertheless, its *d*C–Si value is still the smallest of the corresponding four slabs, whether in Gr/Si, Li/Gr/Si, or Gr/Li/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV). In these Gr/Si slabs, the adsorption of Li has the greatest influence on the *d*C–Si value of SV/Si. The *d*C–Si value decreases when Li is adsorbed on the SV/Si surface; its value further decreases when Li passes through the SV. The total reduction is more than 0.2 Å. It may be due to the larger change of the C positions in SV during Li diffusion to SV/Si. The existence of the Si substrate also slightly affects the distances (*d*Li–C) from the Li to Gr sheet. With the Si substrate, the values of *d*Li–C increase for Li/*p*–Gr/Si but decrease for Li/*d*–Gr/Si (Gr = DV5–8–5, DV555–777 and SV). The value of *d*Li–C in Li/SV/Si reduces most, about 0.1 Å. For monolayer Gr (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), the *d*Li–C values are

basically unchanged due to Li being in the symmetrical position of Li/Gr after Li passes through the Gr sheet. Nevertheless, the *d*Li–C values in Gr/Li/Si increase, except DV555– 777/Li/Si, when there is the Si substrate. The *d*Li–C value in DV5–8–5/Li/Si increase the most, over 0.4 Å. Considering the distance between Gr and the surface of the Si substrate (*d*C–Si), Li is generally closer to the Si surface, except DV555–777/Li/Si. In proportion, the distance from Li to Gr is SV > *p*–Gr > DV5–8–5 > DV555–777. Comparing with the other three Gr/Si slabs, Li in DV555–777/Li/Si is closer to DV555–777. From this point of view, Li seems to prefer to bond to DV555–777, or it may be less easily bonded to the Si substrate in DV555–777/Si. For *d*Li–Si in Gr/Li/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), in proportion, the order is the opposite of *d*Li–C mentioned above, i.e., SV < *p*–Gr < DV5–8–5 < DV555–777. The layer distances indirectly show the bonding between Li and Gr (or Si substrate) in the four Gr/Si slabs above, which can help us better understand the energy stability of these structures.

**Table 2.** The layer distances (Å) of Gr and Si surface of Si substrate (*d*C–Si) and the distances (Å) between Li atom and Gr (*d*Li–C) or Si surface (*d*Li–Si) during Li diffusion in Gr and Gr/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV).


From the bond lengths and layer distances analysis above, we can see that the SV defect is special. Although the structural stability of Li/SV is high, the adsorption and diffusion of Li have a great impact on the planar structure of SV when the Si substrate exists. In the process of Li diffusion, the position of the C atom around the defect, structural configuration, bond lengths and layer distances in SV/Si changes greatly. In other words, the structure of SV/Si changes greatly, and its structural stability is low during Li diffusion in it. For the other three Gr/Si slabs, the structural configurations, bond lengths and layer distances are relatively stable during Li diffusion in them, that is, the influence of the Si substrate on them can be basically ignored.

#### *3.2. Li Diffusion in Gr and Gr/Si (Gr = p–Gr, DV5–8–5, DV555–777 and SV)*

In order to better illustrate the difficulty of Li diffusion in Gr and Gr/Si (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), we calculate their diffusion barrier as shown in Figure 5. As a reference, we firstly study Li diffusion in *p*–Gr and *p*–Gr/Si. It is the most stable structure when Li is adsorbed on the H site of *p*–Gr and the *d*Li–C is about 1.755 Å as shown in Table 2, which is chosen as the initial point (Li/*p*–Gr) in the CI-NEB calculations. The end point (*p*–Gr/Li) is the symmetric adsorption position on the other side of the *p*–Gr plane. The calculation results show that the energy barrier of the Li diffusion in *p*–Gr is as high as 7.559 eV, shown in Figure 5a, which indicates that it is quite difficult for Li to penetrate

through the *p*–Gr sheet. For *p*–Gr/Si, the most stable structure when Li is adsorbed on it (Li/*p*–Gr/Si shown in Figure 3a) is still chosen as the initial point, and *p*–Gr/Li/Si in Figure 3c is the end point in the CI-NEB calculation after Li intercalates *p*–Gr. From Figure 5, we can see that the energy barrier for Li diffusion in *p*–Gr/Si is higher (0.176 eV) than in the case of the *p*–Gr sheet. That is, the presence of the Si substrate increases the energy barrier for Li penetration through *p*–Gr. These calculation results suggest that in the model (Si microparticles encapsulated with graphene cage), it is very difficult for Li to diffuse into the *p*–Gr layer when Si particles are encapsulated with *p*–Gr. Furthermore, it is even more difficult for Li to diffuse out once it is intercalated between *p*–Gr and the Si substrate since the energy barrier of diffusion is higher, shown in Figure 5b. It can be seen that it does not seem to be a good choice to encapsulate Si particles by *p*–Gr. In that way, what about the common defective graphene (*d*–Gr)?

Defects will inevitably occur in the process of graphene producing, and DV5–8–5, DV555–777 and SV are the common low-energy defect structures on graphene. Therefore, we also study the diffusion behavior of Li diffusion in *d*–Gr and *d*–Gr/Si (Gr = DV5–8–5, DV555–777 and SV), which is shown in Figure 5c–h. In the CI-NEB calculations, the lowest-energy structures of Li adsorbed on *d*–Gr (or *d*–Gr/Si) are used as the initial point structures. Similar to *p*–Gr, for *d*–Gr, their end points are the symmetric adsorption position on the other side of the plane, and for *d*–Gr/Si, they are the *d*–Gr/Li/Si structures shown in Figures 3 and 4. The calculation results show that the defective graphenes greatly reduce the energy barriers when Li diffusion in graphene. Comparing with the energy barrier of Li diffusion in *p*–Gr, the corresponding values of DV555–777 and SV are significantly lower by ~4.723 eV and ~4.92 eV, respectively. The value of DV5–8–5 decreases the most, and its diffusion barrier is only 1.406 eV, which is lower (~6.153 eV) than the value of *p*–Gr. How about encapsulating Si particles using these defective graphenes since defects can reduce the energy barrier of Li diffusion into graphene? *d*–Gr/Si is used to simulate the encapsulating structures. When Li diffuses in DV555-777/Si and DV5–8–5/Si, the energy barriers are slightly increased compared with the case of monolayer DV555-777 and DV5–8–5. Nevertheless, their diffusion barriers decrease by 4.705 eV and 6.276 eV than the value in *p*–Gr/Si, respectively. For SV/Si, the energy barrier of Li diffusion to it is unexpectedly decreased to 0.587 eV. It is probably because the position of the C atom around the SV defect changes greatly during Li diffusion, and SV no longer maintains planar configuration. In addition, the structures of Gr/Li/Si are more stable after Li penetrating in the Gr sheet than adsorption on the surface of Gr/Si, whether for *p*–Gr/Si or *d*–Gr/Si. It indicated that it is more difficult for Li to diffuse out from Gr/Si. However, it is much easier for Li to diffuse out from *d*–Gr/Si than from *p*–Gr/Si, shown in Figure 5. From the perspective of the diffusion barrier, the SV defective graphene seems to be a good material for encapsulating Si particles. However, our calculation results suggest that the C atoms and surface Si atoms around SV defects in SV/Si greatly shift during Li diffusion. The times of charge and discharge are reduced and the battery life is shortened if Si particles are encapsulated with SV. On the other hand, the energy barrier of Li diffusion from SV/Si is not low; it is even slightly higher than the value of DV5–8–5/Si. In brief, the probability of Li penetration through graphene would be greatly enhanced if Si particles are encapsulated with defective graphene instead of perfect graphene. Considering the energy barriers of Li diffusion into and out from Gr/Si and the stability of the structures of Gr/Si during Li diffusion, DV5–8–5 defective graphene may be a good choice as the material for encapsulated Si particles.

**Figure 5.** Energy barrier of Li diffusion.

In order to better understanding the electronic structures of DV5–8–5, we also calculate the Bader charge, band structures and PDOS of these structures during Li diffusion in DV5–8–5 and DV5–8–5/Si. As a reference, the electronic structures of corresponding *p*–Gr are also calculated, which are shown in the supplementary material. The DV5–8–5 defect exhibits a Dirac-like cone band structure, but a gap opens between the Dirac-like cones as shown in Figure S2a. There is a nearly flat band above the EF level and within the gap between the Dirac-like cones. When the Li atom adsorbs on it (i.e., Li/DV5–8–5), the band in the gap becomes more dispersive and partially occupied as shown in Figure S2b. Meanwhile, the gap between the Dirac-like cones reduces. With the diffusion of the Li atom (Li@DV5–8–5), the gap between the Dirac-like cones becomes wider as shown in Figure S2c. After Li penetration is completed, the band structure is recovered as shown in Figure S2b. During Li diffusion, the Li atom maintains the Li<sup>+</sup> ion. When the Si substrate is included, the DV5–8–5 still exhibits a Dirac-like cone structure. The gap between the cones becomes wider. There are three flat bands in the gap due to the influence of the Si substrate, and these bands are all above EF as shown in Figure S2d. During Li diffusion in DV5–8–5/Si, the system remains as p–type doping; the Li atom maintains the Li+ ion; the bands in the gap become more complex; and some of the bands become more dispersive and partially occupied. Similar to the case of DV5–8–5, the gap between the Dirac-like cones experiences a process of first reducing (Li/DV5–8–5/Si), then increasing (Li@DV5–8–5/Si), and finally basically recovering (DV5–8–5/Li/Si).

#### **4. Conclusions**

During Li diffusion in the Gr sheet and Gr/Si slab (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), the initial points, saddle points and end structures are studied by first-principles calculation. Firstly, the initial and end points of different structures are screened, and their energetic stability is also discussed. Then the CI-NEB calculations are employed to evaluate the diffusion barrier by seeking the saddle point and the minimum energy path. For monolayer Gr, Li prefers to adsorb on *d*–Gr (Gr = DV5–8–5, DV555–777 and SV), comparing with *p*–Gr. The *E*ads of Li/DV5–8–5 and Li/DV555–777 is more than 1 eV lower than that of Li/*p*–Gr; the *E*ads of the most stable adsorption structure Li/SV is nearly 1.8 eV lower than that of Li/*p*–Gr. The *E*ads order of Li/Gr/Si is consistent with that of Li/Gr, except Li/SV/Si. It is difficult to adsorb on SV/Si for Li, as its *E*ads is even higher about 1 eV than that of Li/*p*–Gr/Si. It may be related to the structural deformation of SV/Si. During Li diffusion in Gr, the Gr sheets basically maintain the planar whether for Gr or Gr/Si, except for SV. When Li diffuses in SV, one of the C atoms is obviously deviated from the SV plane and approaches the Li atom whether in the initial point Li/SV or the end point SV/Li. More C atoms deviate from the SV plane and approach the Li atom when Li diffuses in SV/Si. Even in the saddle point Li@SV and Li@SV/Si, the C atoms around Li are almost deviated from the original plane, i.e., during Li diffusion in SV and SV/Si, the SV plane structure is seriously deformed; moreover, several Si atoms on the surface of Si substrate deviate from their original positions in SV/Li/Si. It may be the reason why many properties of Li diffusing in SV or SV/Si are different from those of the other Gr, including the order of *E*ads, bond lengths between atoms, layer distance and even the diffusion barrier. For monolayer Gr (Gr = *p*–Gr, DV5–8–5 and DV555–777), the average bond lengths of C–C and Li–C in the saddle-point structures (Li@Gr) are larger than the ones of the corresponding initial structures (Li/Gr) due to Li being in the plane of Gr. When Li penetration is complete, the bond lengths are basically restored to the values of the initial structures. In Gr/Si, the Si substrate has negligible effects on these bond lengths. Comparing with the three defects above, the average C–C and Li–C bond lengths of the saddle-point structure Li@SV are not increased but decreased, although the Li atom is also in the SV plane. In Li@SV/Si, the average C–C bond length increases slightly, while its Li–C average bond length decreases over 0.3 Å compared to the initial structure. After Li penetrates the SV plane, the C–C average bond length does not return to the initial value but is increased. The distances between the Gr sheet and the surface of

the Si substrate (*d*C–Si) in Gr/Si (Gr = *p*–Gr, DV5–8–5 and DV555–777) change little during Li diffusion in them. However, the *d*C–Si value of SV/Si decreases about 0.1 Å upon Li adsorption, which is about 10 times that of other structural variation. When Li passes through the SV, its *d*C–Si value further decreases, and the total reduction is more than 0.2 Å. The existence of the Si substrate also slightly affects the distances (*d*Li–C) from Li to the Gr sheet. Comparing with Li/Gr (Gr = *p*–Gr, DV5–8–5, DV555–777 and SV), the *d*Li–C in Li/*p*–Gr/Si increases, while decreasing in Li/*d*–Gr/Si. The *d*Li–C in Li/SV/Si decreases most, about 0.1 Å. After passing through Gr, Li is generally closer to the Si surface, except for DV555–777/Li/Si. Considering the distance between Gr and the Si surface (*d*C–Si), the distance from Li to Gr is SV > *p*–Gr > DV5–8–5 > DV555–777 in proportion. Accordingly, the order of *d*Li–Si is opposite to that above. It indirectly provides evidence that Li binds more closely to the Si surface in SV/Li/Si while bonding more tightly to the DV555–777 layer in DV555–777/Li/Si. The CI-NEB calculation results show that *d*–Gr greatly reduces the energy barriers whether Li diffuses into Gr or Gr/Si. The energy barrier of Li diffusion in DV5–8–5 is just 1.406 eV, which is lower by ~6.153 eV than the value of *p*–Gr. When Li diffuses into Gr/Si, the energy barriers are slightly increased comparing with the case of monolayer Gr. Nevertheless, the diffusion barriers of *d*–Gr/Si are still much lower than the value in *p*–Gr/Si. The energy barrier of Li diffusion to SV/Si is unexpectedly decreased to 0.587 eV. In addition, the structures of Gr/Li/Si are more stable than Li/Gr/Si whether for *p*–Gr or *d*–Gr, which indicates that it is more difficult for Li to diffuse out from Gr/Si. However, it is much easier for Li to diffuse out from *d*–Gr/Si than from *p*–Gr/Si. Our results suggest that the Si particles encapsulated by defective graphene may be a candidate for anode materials of LIBS. From the perspective of the diffusion barrier, both SV and DV5–8–5 defective graphenes seem to be good materials for encapsulated Si particles, and SV is better. However, the SV/Si structure is deformed greatly during Li diffusion, which reduces the times of charge and discharge and shortens the battery life. Moreover, the energy barrier of Li diffusion out from SV/Si is not low and is even slightly higher than the value of DV5–8–5/Si. In contrast, DV5–8–5/Si remains a relatively stable structure during Li diffusion within it. Therefore, DV5–8–5 defective graphene may be a good choice as the material of encapsulated Si particles.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/nano11123397/s1: Table S1: The bader charges of Li atom, Gr (Gr = *p*–Gr and DV5–8–5) sheet and Si substrate of the initial, saddle and end structures during Li diffusion in Gr and Gr/Si. Figure S1: Band structures of (a) *p*–Gr (b) Li/*p*–Gr (c) Li@*p*–Gr (d) Si Substrate (e) *p*–Gr/Si (f) Li/*p*– Gr/Si (g) Li@*p*–Gr/Si (h) *p*–Gr/Li/Si. Figure S2: Band structures of (a) DV5–8–5 (b) Li/DV5–8–5 (c) Li@DV5–8–5 (d) DV5–8–5/Si (e) Li/DV5–8–5/Si (f) Li@DV5–8–5/Si and (g) Li/DV5–8–5/Si. Figure S3: The line configurations of the initial, saddle and end structures during Li diffusion in Gr and Gr/Si (Gr = *p*–Gr and DV5–8–5); and the typical C atoms discussed in Figures S4–S11 are labeled. The green ball represents Li atom; and the brown and blue lines represent C atoms and Si atoms and their bonding, respectively. Figure S4: PDOS of the typical C atoms (labeled in Figure S3) of p–Gr and all Si atoms of Si substrate. Figure S5: PDOS of Li atom and the typical C atoms (labeled in Figure S3) in (a) Li/p–Gr and (b) Li@p–Gr. Figure S6: PDOS of the typical C atoms (labeled in Figure S3) and all Si atoms in p–Gr/Si. Figure S7: PDOS of Li atom, all Si atoms and the typical C atoms (labeled in Figure S3) in (a) Li/p–Gr/Si (b) Li@p–Gr/Si and (c) p–Gr/Li/Si. Figure S8: PDOS of the typical C atoms in DV5–8–5 labeled in Figure S3. Figure S9: PDOS of Li atom and the typical C atoms (labeled in Figure S3) in (a) Li/DV5–8–5 and (b) Li@DV5–8–5. Figure S10: PDOS of all Si atoms of the typical C atoms (labeled in Figure S3) in DV5–8–5/Si. Figure S11: PDOS of Li atom, all Si atoms and the typical C atoms (labeled in Figure S3) in (a) Li/DV5–8–5/Si (b) Li@DV5–8–5/Si and (c) DV5–8–5/Li/Si.

**Author Contributions:** Data curation, investigation, writing—original draft preparation, reviewing and editing, W.Q.; supervision, W.-C.L.; investigation, X.-Y.X.; reviewing and editing, K.-M.H.; reviewing and editing, supervision, C.-Z.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Nos. 21603114 and 21773132) and US DOE-BES: DE-AC02-07CH11358.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Work at Ames Laboratory was supported by Department of Energy, Office of Science, Basic Energy Sciences, Division of Materials Science and Engineering including a grant of computer time at the National Energy Research Scientific Computing Center (NERSC) in Berkeley. Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. DE-AC02-07CH11358.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Interface Kinetics Assisted Barrier Removal in Large Area 2D-WS2 Growth to Facilitate Mass Scale Device Production**

**Abid 1, Poonam Sehrawat 1, Christian M. Julien 2,\* and Saikh S. Islam 1,\***


**Abstract:** Growth of monolayer WS2 of domain size beyond few microns is a challenge even today; and it is still restricted to traditional exfoliation techniques, with no control over the dimension. Here, we present the synthesis of mono- to few layer WS2 film of centimeter2 size on graphene-oxide (GO) coated Si/SiO2 substrate using the chemical vapor deposition CVD technique. Although the individual size of WS2 crystallites is found smaller, the joining of grain boundaries due to *sp*2 bonded carbon nanostructures (~3–6 nm) in GO to reduced graphene-oxide (RGO) transformed film, facilitates the expansion of domain size in continuous fashion resulting in full coverage of the substrate. Another factor, equally important for expanding the domain boundary, is surface roughness of RGO film. This is confirmed by conducting WS2 growth on Si wafer marked with few scratches on polished surface. Interestingly, WS2 growth was observed in and around the rough surface irrespective of whether polished or unpolished. More the roughness is, better the yield in crystalline WS2 flakes. Raman mapping ascertains the uniform mono-to-few layer growth over the entire substrate, and it is reaffirmed by photoluminescence, AFM and HRTEM. This study may open up a new approach for growth of large area WS2 film for device application. We have also demonstrated the potential of the developed film for photodetector application, where the cycling response of the detector is highly repetitive with negligible drift.

**Keywords:** chemical; vapor deposition; graphene oxide; transition-metal dichalcogenides; WS2

#### **1. Introduction**

Lamellar two-dimensional (2D) WS2 has immense potential for electronic applications because of its remarkable properties such as tunable direct bandgap [1–3], high photoluminescence intensity [4,5], high emission quantum field [6], attractive spin-orbit coupling [6], substantial exciton/trion binding energies [7–9], etc., which give WS2 superior leverage over other transition-metal dichalcogenides (TMDCs). Unfortunately, the research on WS2 is not matured yet, especially large-scale production of single-crystal monolayers [10–12]. Challenges lie in reliable synthesis of atomically thick 2D layers and controlled manipulation of electronic properties [13]. Though mechanical cleaving still remains a relatively facile technique to prepare high quality single to few layered WS2, problems such as absence of control over thickness, relatively small lateral size and low yield limit its applications and commercialization [14,15]. The most popular techniques to produce few-layered WS2 include intercalation driven exfoliation [16], liquid phase exfoliation [17], physical vapor deposition [18], hydrothermal reaction [19], and heat treatment of W and S containing precursors [20]. However, the domain size of developed WS2 via aforementioned methods is generally restricted to a few microns (μm), and the production of WS2 monolayer films with desired domain size remains a distant dream. To achieve this, several strategies (e.g., atomic layer deposition [21], pulsed-laser deposition [22,23], metal/metal oxide thin films [24], suitable metal substrates [25], etc.), have been suggested

**Citation:** Abid; Sehrawat, P.; Julien, C.M.; Islam, S.S. Interface Kinetics Assisted Barrier Removal in Large Area 2D-WS2 Growth to Facilitate Mass Scale Device Production. *Nanomaterials* **2021**, *11*, 220. https://doi.org/10.3390/ nano11010220

Received: 30 December 2020 Accepted: 13 January 2021 Published: 16 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to achieve uniform dispersion of precursors over substrates. Presently, large scale continuous growth of TMDCs is highly complex and expensive, and results in poor quality films [9].

Chemical vapor deposition (CVD) remains a challenge for a long time for growing two-dimensional TMDC materials with high crystallinity, desired thickness, and sufficient domain size [10,12,26–28]. For WS2 growth, typical synthesis approach involves sulfidation of WO3 powders at sufficiently high temperature. The growth conditions such as precursor type, process step, and their manipulations during growth, are found to greatly affect the product quality. Most of the reports on CVD grown WS2 solely bank on vapor phase reaction between S and W precursors at suitably high temperatures, where a single furnace is used for both precursors. In such systems, there is no control over the temperature of individual precursor, thereby limiting the parameter space of growth reaction [28]; leading to poor uniformity and repeatability. High melting point of WO3 powders and sulfurization rate are found to influence the growth process. Achieving large triangle shaped growth is challenging, because of lower vapor pressure of WO3 due to its high melting temperature (1473 ◦C) [29]. This is a serious problem and severely lowers the partial pressure of WO3 vapors [29]. Further, the low vapor pressure of WO3 also hampers the availability of W atoms on substrate surface. Thus, enhancing the partial pressure of WO3 vapors is a natural approach to obtain large sized WS2 flakes [12]. This can be accomplished by reducing the pressure of the furnace. However, this increases the S vapor transfer speed, thereby increasing the sulfurization rate [10]. Fast sulfurization is not desirable for the transport and diffusion of precursor vapors on a substrate and may result in aggregation of precursors at certain locations on the substrate, thereby limiting the growth of large sized WS2 crystals [10,30,31]. An efficient way involves a trade-off between the increase in partial pressure of WO3 vapors and lower the transfer speed of S vapors.

Recently, many researchers have focused on developing monolayer WS2 via CVD. Zhang et al. employed low-pressure chemical vapor deposition (LPCVD) technique to synthesize atomically thin triangle shaped WS2 crystals on sapphire substrate having single domain size ~50 μm [32]. An improvement was introduced in the CVD process by Cong et al. to effectively increase the concentration of precursors and furnace pressure [33]. They achieved single-domain growths of ~178 μm [33]. Li et al. introduced further changes in the CVD reaction by employing alkali metal halides as growth promoters [34]. The groups of Zhang [32] and Li [34] reported that due to strong reducing nature, the mixing of suitable amount of H2 enhances the sublimation and sulfurization of WO3 precursor. Fu et al. [35] investigated the influence of CVD reaction temperature and gas flow rate (comprising a mixture of 97% Ar and 3% H2) on the growth size and morphology of WS2 films. The optimized growth conditions yielded a domain size of ~52 μm. Rong's group [28] achieved precise control over S introduction time by utilizing a two-zone furnace and obtained very large area WS2 films of 370 μm domain size. Recently, Liu et al. [36] reported sequential synthesis of several one-, and two-dimensional nanomaterials by intentionally creating initial low sulfur conditions, to divide the growth process in two stages: the first stage is a partial reduction and the second one is sulfurization. By this, they could exploit a WO3-*<sup>x</sup>* intermediate route, and were able to grow large size WS2 (150 μm).

TMDCs growth is very susceptible to substrate treatment prior to formation [37]. To assist the nucleation by seeding the substrate with various aromatic molecules, such as -3,4,9,10-tetracarboxylic acid tetrapotassium salt (PTAS), perylene-3,4,9,10-tetracarboxylic dianhydride (PTCDA), and reduced graphene oxide (RGO), perylene encourages the lateral growth of TMDCs. However, this seeding technique helps to grow the large area mono-to few layer TMDCs but continuous film from the abovementioned technique is still limited. Although the CVD method has many advantages, the coordination among various growth parameters is highly complex and requires further elucidation.

In this work, we present a novel approach for large area synthesis of WS2 on graphene oxide (GO) surface via CVD under high temperature conditions of 1070 ◦C. Upon transformation to RGO during growth process, GO becomes a mixed *sp*3-*sp*<sup>2</sup> bonded network, where nano-sized *sp*2-patches of size ~3–6 nm, are randomly distributed with more than 60–70% coverage [38–40]. The *sp*2-patches are six-fold aromatic carbon atoms with C=C bonds, whereas major percentage of *sp*<sup>3</sup> bonds constitutes epoxides (C–O) and hydroxyl (−OH) groups in the basal plane and somewhat lesser amount of carboxyl (−COOH) and carbonyl (−C=O) functionalities in the edges [38–40]. These *sp*2-chemical bonds act as a seed to promote the growth of large area WS2 film. The as-developed WS2 layers are polycrystalline comprising mainly mono- to few layers. The single-crystal domains are triangular and hexagonal in shape, with sizes up to 60 μm on Si/SiO2 and ~15 μm on the GO coated Si/SiO2, respectively. The domain size of WS2 film on GO coated Si/SiO2 is much smaller than that on non-GO substrate, but the development of large continuous surface makes a great footprint in WS2 growth, attributed to C=C assisted coalescence of high-density polycrystalline grain boundaries into large domain area. Suitable explanation is given how *sp*2-patches comprised of C=C bonds are pivotal in the expansion of domain size of WS2 film. Surface roughness of the top layer of substrate is a well-recognized technique for nucleation of nanostructure growth and is cross-checked with intentionally created rough surface on silicon wafer. The proposed device fabrication technique is well suited for mass production of identical devices, and therefore bids for commercial scale development. Moreover, this technique can be extended to the fabrication of other TMDCs. We have further demonstrated the use of the developed film for photodetector applications.

#### **2. Materials and Methods**

WO3 powders (≥99.995% trace metals basis), sulfur (≥99.998% trace metals basis), natural graphite powders (flakes size: 45 μm, ≥99.99% trace metals basis), and sodium nitrate (NaNO3, ≥99.995% trace metals basis) were commercially secured from Sigma Aldrich, India. Sulfuric acid (H2SO4) and potassium permagnate (KMnO4) have been supplied by Thermo Fisher Scientific, India and DI water was obtained from Merck (Darmstadt, Germany). For experimental purposes, only analytical grade materials/chemicals were procured and used as received.

The morphology, size and structure of as-synthesized WS2 were examined via a transmission electron microscope (TEM) (JEOL JEM F-200, Tokyo, Japan) working at 200 kV. The samples for TEM investigations were kept on the conventional holey carbon Cu grid. Highresolution transmission electron microscopy (HRTEM, JEOL, Tokyo, Japan) studies reveal the uniformity and high quality of WS2 samples grown on both Si/SiO2 and GO coated Si/SiO2 substrates. Basic studies of WS2 growth and the domain size were performed on an optical microscope (Leica DM4 P, Wetzlar, Germany). Raman spectrophotometers (WITec Alpha 300RA, Ulm, Germany and LabRAM HR800 HORIBA JY, Kyoto, Japan) fitted with a Peltier cooled CCD detector and confocal microscope were operated for verifying the crystallinity, defects and number of layers. A diode laser of excitation wavelength 532 nm was used to excite the samples. A 100X objective lens with spot size of 1 μm was used to focus the excitation onto the samples. The thickness of grown sample was analyzed by atomic force microscopy (AFM) fitted with the Raman system. Surface morphology of the CVD grown WS2 on Si/SiO2 and GO coated Si/SiO2 was studied on a scanning electron microscope (SEM) (Nova Nano SEM 450, FEI, Oregon, USA). The photoluminescence (PL) measurements were performed by fluorescence spectrometer (Agilent Technologies, Carry Eclipse fluorescence spectrometer, California, USA). For PL investigations, the prepared samples were dispersed in ethanol by simply ultrasonicating the growth substrate for few minutes and isolating the upper one-third solution into a quartz cuvette.

#### **3. Results and Discussion**

Two approaches have been adopted to synthesize monolayer WS2 of large domain size, where CVD is employed to grow WS2 on two substrates, i.e., Si/SiO2 and GO-coated Si/SiO2 (Si/SiO2/GO), with the aim to extend the domain size of WS2 and uniform coverage of the entire substrate area (see Sections 3.1 and 3.2).

#### *3.1. First Approach: Mono- to Few-Layered WS2 Growth on Si/SiO2*

WS2 growth was carried out inside a two-zone tube furnace having separate highand low-temperature zones, denoted as HT and LT, respectively, as shown in Figure 1a. Ultrahigh purity argon at flow rate of 100 sccm under atmospheric pressure was used as carrier gas to transport sulfur vapors to the HT zone. In a typical experiment, WO3 powders (200 mg of pure grade 99.9%) was placed in an alumina crucible at the center of HT-zone of the furnace. Si/SiO2 substrate (10 mm × 10 mm) was placed on the alumina crucible in inverted position so as to allow WO3 vapor to hits the substrate and deposit on it. The sulfur powder (400 mg of pure grade 99.5%) was loaded upstream, and heated in LT-zone of the furnace. The HT- and LT-zone temperatures were controlled in such a way that the temperature in both the zones simultaneously reach at 1070 ◦C and 220 ◦C, respectively [28]. Monolayer WS2 grains were grown on Si/SiO2 substrate having small domain size and the growth remained scattered. The reaction occurring at the surface of the substrate can be summarized as:

$$2\text{WO}\_3 + 7\text{S} \to 2\text{WS}\_2 + 3\text{SO}\_2.\tag{1}$$

**Figure 1.** (**a**) Experimental setup. (**b**–**e**) Optical images of synthesized WS2 film showing the presence of mono-to-few layered WS2 on Si/SiO2 substrate.

Thus, the vapor phase reaction between tungsten oxide and sulfur, produces tungsten disulfide and sulfur dioxide as end products. The optimized time duration for WS2 growth is found to be 45 min, once the temperature of HT zone reaches 1070 ◦C. After completing the growth process, the CVD system was allowed to cool naturally to the room temperature.

Figure 1b–e shows the optical images of the WS2 sample grown on Si/SiO2 substrate. Two types of crystallites are visible: (a) equilateral triangles, and (b) hexagons; both are characteristic of WS2 single crystal formation. Most of the crystals include monolayers. However, bi- and tri-layered regions are also visible at the center of large monolayer crystals, which is supposed to be the nucleation site, even for further growth of multi-layers [41].

Formation of hexagons in case of multilayer grains originates from the random orientation of individual WS2 layers and differences in crystal orientation, where each orientation acts as nucleation site [41–43]. The side length of the crystals in case of equilateral triangles is found to be ~60 μm, which is reasonably large. The CVD growth of TMDCs is generally believed to be a self-limiting process [44]. The lateral size of WS2 increases with the partial pressure of chalcogen precursors [45], while its thickness is determined by the partial pressure of transition metal precursors [46]. We have kept the weight of sulfur powders to be the twice of the weight WO3 (chalcogen precursor). Thus, the lateral growth dominates in our case and could explain the growth of large size single crystals [46]. Furthermore, we could find WS2 single crystals in variable sizes, ranging from 100 nm to 60 μm on the same Si/SiO2 substrate. The growth duration is found to directly influence the area

of grown crystals, before the crystal growth reaches its limit. Therefore, large sized WS2 crystals can be formed if the growth time is kept sufficiently long.

• Electron Microscope Analysis

Figure 2a–c shows the FESEM images of WS2 samples deposited on Si/SiO2 substrate, corroborating the findings of optical microscopy, indicating both equilateral triangles and hexagons. In case of WS2, the growth of equilateral triangles suggests single crystallinity; whereas the hexagonal or star-like structures result from several rotationally symmetric grains [4,9,26,46,47]. The size distribution histogram of hexagonal and triangle flakes of WS2 is shown in Figure 2d–e. Where the average size of the hexagonal and triangular flake is found to be 71.5 μm and 44.5 μm, respectively. The ratio of the hexagon to the triangular flakes over the substrate is 2:1. Notably, the growth products, though consisting of monolayers, are scattered over the substrate and do not cover the entire substrate surface.

**Figure 2.** (**a**–**c**) FESEM micrographs of the WS2 layers grown on Si/SiO2 substrate at various magnification levels and (**d**,**e**) Size distribution histogram of hexagonal and triangle WS2 flakes.

Figure 3a–c shows a low magnification bright field TEM imaging of WS2 crystals suspended on a holey carbon Cu grid. For TEM imaging, first the WS2 grown onto the Si/SiO2 substrate was ultrasonicated in ethanol for 20 min to get the WS2 flakes. After this, the suspension was centrifuged at a speed of 3000 rpm for 15 min to remove unwanted contamination or to allow settling down of heavy particles. Finally, top onethird supernatant was used for TEM analysis. We used HRTEM and selected area electron diffraction (SAED) to characterize the WS2 domains. Figure 3a shows the hexagonal flakes with lateral dimension of about ~2.5 μm. In Figure 3b, the distorted cum truncated WS2 triangle is clearly observable. At higher resolutions (Figure 3c,d), the crystalline nature of the flakes is clearly displayed. In Figure 3e, the selected area electron diffraction (SAED) pattern is indicated, which further confirms the highly crystalline morphology of as-grown WS2.

**Figure 3.** The TEM images of WS2 flakes shown in (**a**,**b**) both at low magnification, (**c**) at high magnification. (**d**) zoomed view of (**c**), and (**e**) SAED pattern of the zoomed area.

• Raman Analysis

The *E-k* diagram and the phonon dispersion in Brillouin zone for WS2 are shown in Figure 4a–c. Raman spectroscopy was employed to detect both the in-plane (*E*2g1) and out-of-plane (*A*1g) phonon modes. Figure 4d shows Raman spectrum of the WS2 flakes deposited over the Si/SiO2 substrate, indicating a strong peak at 354 cm−<sup>1</sup> (2LA(*M*)), signifying a longitudinal acoustic mode at *M* point. This signals the monolayer growth of WS2 [9,26]. The peak at 421 cm−<sup>1</sup> corresponding to the *A*1g mode, although weak, is also observed. As is widely reported, the separation between these two peaks increases from 67 to 72 cm−<sup>1</sup> with increasing the number of layers from single- to tri-layer due to a redshifting of 2LA(*M*) mode [26,48–51]. However, the observed peak separation of ~67 cm−<sup>1</sup> is the fingerprint of WS2 monolayer growth. Its topography, revealing a triangular shape probed via Atomic Force Microscopy (AFM), is shown in Figure 4e,f. The whole crystal is almost atomically flat and clean. The height profile in Figure 4f indicates the as-grown WS2 on Si/SiO2 having a step size of 0.89 nm; reaffirming the growth of a WS2 monolayer. To verify the uniformity in the flakes, Raman mapping studies were performed from one edge to the other (Figure 5a,b). The synthesized WS2 flakes were further analyzed by recording multiple Raman spectra across the area of post WS2 grown substrate. Typically, seven measurement sites were selected on the top surface, and it can be concluded that we have achieved the monolayered WS2 flakes. As seen in Figure 5a,b, the spots 1 and 7 only show the Raman-active LO mode of Si, while spots 2, 3, 4, 5, and 6 express the Raman features of WS2 monolayers, and interestingly, the Raman result is also found supplemented by AFM studies.

**Figure 4.** (**a**) Real space lattice structure and Brillouin zone in reciprocal space, (**b**) electronic band structure (*E-K* diagram), (**c**) phonon dispersion, (**d**) Raman spectrum, (**e**) AFM image of WS2 monolayer crystal grown on Si/SiO2 substrate and (**f**) AFM cross-section height profile for the deposited sample revealing a thickness of the grown film to be ~0.89 nm.

**Figure 5.** (**a**) Optical microscope image of WS2 deposited over Si/SiO2, indicating the selected measurement spots, and (**b**) Raman spectra of the spots identified in (**a**) where peak at 520 cm−<sup>1</sup> in all the spectra is the LO mode of the silicon wafer.

• Photoluminescence (PL) Analysis

The PL properties of 2D-TMDCs are mainly governed by their excitonic features [1,7– 9,52]. In the limiting case of 2D materials, as shown in Figure 6a, the lines of electric field between the hole and electron start extending beyond the material, thereby reducing the dielectric screening and enhancing coulomb interactions [36,38,39]. A strong coulombic force is the main reason behind strong photoemission efficiency. This enhanced coulomb interaction suggests formation of bound e-h pairs or "excitons", having greater binding energies making them stable even at room temperature [7,52–54]. In case of monolayer 2D material, significant deviation is expected from the hydrogenic model, prevalent to describe Wannier-excitons in inorganic semiconducting materials [7,51]. Therefore, the calculation of binding energy of excitons in case of single-layered WS2 needs to consider modified potential [54]. High intense peak is observed at 630 nm (1.968 eV) in the as-synthesized

monolayer WS2, which mainly originates from neutral-exciton (*X*A) emission due to direct transition between the lowest conduction band (CB) and the highest valence band (VB) at the same *K*-point in the Brillouin zone [7,49,55–58], and schematically represented in Figure 4c. The fundamental bandgap of a single layer WS2 is ~2.1 eV, but the large excitonic binding energy of 0.032 eV renders the optical bandgap measured at 1.968 eV [8,50–52] as shown schematically in Figure 6a,c shows the PL spectra of WS2 flakes on the Si/SiO2 substrate. Figure 6b represents the Neutral A-excitonic transition (XA) in the energy level diagram. The FWHM values of the PL spectra lie between 25 to 35 meV for both the approaches which is in accordance with the literature and further establish the fact of monolayer formation [59,60]. The presence of enhanced coulombic interaction at higher distant points suggests more likelihood of formation of higher-order excitons, i.e., trions and bi-excitons [51].

**Figure 6.** (**a**) Schematic diagram showing the reduced dielectric screening on the e-h pair bonding in a typical 2D material [38–40]. (**b**) Neutral A-excitonic transition (*X*A) in the energy level diagram, and (**c**) PL spectra of WS2 grown over Si/SiO2 substrate.

#### *3.2. Second Approach: Uniform Growth of WS2 Flakes on Spin-Coated GO*

In this approach, we have spin coated pre-synthesized graphene oxide (GO) solution over the Si/SiO2 substrate prior to CVD growth. Hummers' method [38–40] has emerged as a well-established technique to prepare graphite oxide which can then be ultrasonicated to yield graphene oxide (GO). Briefly, natural graphite powders (5.0 g), H2SO4 (120 mL) and NaNO3 (2.5 g) were mixed via continuous stirring. During stirring, the temperature was kept below 10 ◦C in order to check overheating of the mixture. After obtaining a uniform mixing, KMnO4 (15.0 g) was added to this mixture by stirring for 24 h at room temperature. A condensed light brown colored slurry was obtained, which was diluted with deionized water (150 mL) and stirring for 2 h at 98 ◦C. The resulting mixture was vigorously washed with HCl and deionized water, to yield graphite oxide. Subsequently, graphite oxide was ultrasonicated (37 kHz, 500 W) for 3 h. The resulting suspension was centrifuged at 5000 rpm and GO was obtained by collecting the top one-third supernatant from the centrifugation product. A suspension of GO was prepared by mixing 5 mg GO in deionized water via ultrasonication. The substrate for WS2 growth was prepared by spin coating (2500 rpm) 100 μL GO suspension on a Si/SiO2 wafer. Figure 7a includes the SEM micrograph of the GO film coated over Si/SiO2 substrate. The GO flakes appear uniformly dispersed over the surface of the substrate, and are expected to act as a seed layer for the growth of WS2 film. In Figure 7b, the Raman spectrum of the GO film is displayed in 500–3500 cm−<sup>1</sup> range, in which the D-, G- and 2D-band are clearly defined at 1354, 1592, and 2712 cm−1, respectively. Further, XRD analysis was performed to confirm the quality of prepared GO. The XRD pattern in Figure 7c shows the peak at 2θ ≈ 11◦, which validates the high quality of the GO flakes.

**Figure 7.** (**a**) FESEM micrograph, (**b**) Raman spectrum, and (**c**) XRD pattern of the GO film deposited on Si/SiO2 substrate before the growth of the WS2 layer.

The WS2 growth was performed similar to that described in Section 3.1 (i.e., WS2 grown at Si/SiO2 substrate). The only difference is the coating of Si/SiO2 substrate with GO, which is expected to get converted into RGO due to in-situ heating encountered in the HT-zone during the growth process. The synthesis process is schematically illustrated in Figure 8a–c. Figure 8d–g shows the Optical images of the WS2 grown on the GO coated Si/SiO2, where in Figure 8d clear interface is shown between WS2 film and the Si substrate. Figure 8h–k shows the FESEM micrographs at different magnification. In this approach, sufficient amount of sulfur vapor is taken into consideration for the reaction, similar to the first approach. Here, it is noteworthy to say that the *sp*2-carbon nanostructures in the RGO matrix fulfill three-pronged objectives: (i) they improve the surface adhesion, thereby enhancing local adsorption of WO3 molecules; (ii) they reduce WO3 to yield WO3-*x*; and finally, (iii) they also act as active sites for heterogeneous nucleation, as demonstrated in Figure 9.

**Figure 8.** (**a**) Spin coating GO solution over Si/SiO2 substrate, (**b**) GO coated Si/SiO2 substrate, (**c**) CVD experimental setup used for WS2 growth, and (**d**–**g**) optical microscope images of WS2 grown on GO coated Si/SiO2 substrate, and (**h**–**k**) FESEM micrographs of WS2 grown over GO coated Si/SiO2 substrate.

**Figure 9.** Schematic diagram to represent WS2 growth mechanism on GO coated Si/SiO2 substrate.

Upon heating, WO3 vapors reach the surface of the substrate. A portion of these vapors diffuses and reacts with *sp*2-bonded carbon patches of RGO, resulting in the formation of an intermediate phase, WO3-*x*, along with the release of CO2. Thus, an increase in the concentration of WO3 vapors, results in simultaneous increase in the CO2 and a corresponding decrease in the size of *sp*2-carbon patches. Further, at elevated temperatures, WO3-*<sup>x</sup>* molecules gain the tendency to migrate and aggregate around the *sp*2-bonded carbon patches. These carbon patches or 'carbon nanostructures' are expected to act as heterogeneous nucleation sites; where nucleation takes place via sulfurization of adjoining WO3-*<sup>x</sup>* molecules. The synthesis kinetics involve heterogeneous and homogeneous nucleation reactions, occurring simultaneously and competing with each other, as the main growth processes [61,62]. Initially, due to a lower surface free-energy barrier, heterogeneous nucleation outpaces the homogeneous nucleation [59]. In addition, the local high concentration of WO3-*<sup>x</sup>* facilitates the formation of bilayer nuclei. In case of perfectly overlapping bilayer formation, the nucleation shifts from the bottom layer to the second layer almost instantly. If growth does not stop at this point, few-or multi-layer WS2 flakes are formed. Nonetheless, a complete understanding of growth mechanics warrants further investigations. The growth process can be summarized in the following chemical reactions [61,62]:

$$\rm{WO\_3 + (x/2)C \to WO\_{3-x} + (x/2)CO\_2} \tag{2}$$

$$\rm{WO\_3 + (x/2)S \to WO\_{3-x} + (x/2)SO\_{2-}}\tag{3}$$

$$\rm{WO\_{3-x} + (7-x)/2S \to WS\_2 + (3-x)/2SO\_2} \tag{4}$$

The use of GO produces substantially different growth characteristics than the isolated islands observed during the first approach. The WS2 islands now coalesce to form nearly continuous films on centimeter scale as is evident in the optical images and SEM micrograph (Figure 8). Apart from the monolayers, multilayered growth is also visible. The multilayered growth is more favored at either the nucleation sites or grain boundaries (the locations shown by pointing arrows). Continuous monoto-few layer growths on centimeter scale are representative of the samples grown via the second approach, where the Si/SiO2 substrate is coated with GO. We have found the size of the individual WS2 crystals less than 15 μm, which is much larger than those reported for h-BN substrates, where a domain size of less than 1 μm is observed [63–65]. Although it is smaller than the WS2 grown on pristine Si/SiO2 substrate, its continuous expansion on the entire substrate has strong merits and fulfils the objectives of this work, as well as the requisite need for device applications.

Figure 10a shows the Raman spectrum of a WS2 film grown on GO coated substrate. The peaks present at 351.2 and 421.4 cm−<sup>1</sup> is relatable with the *E*<sup>1</sup> <sup>2</sup>*<sup>g</sup>* and *A*1g Raman modes. Besides, the *E*<sup>1</sup> <sup>2</sup>*<sup>g</sup>* and *A*1g modes are found to be slightly red- and blue-shifted, respectively. Our concern is to study the shift in *E*<sup>1</sup> <sup>2</sup>*<sup>g</sup>* and *A*1g peaks due to thickness variation in a mixed multi-layered WS2 flake. The increase in number of layers strongly enhances the out-ofplane vibrations, while coulomb interactions tend to decrease the frequency of in-plane vibrations, thereby causing monotonous increase in frequency separation between *E*<sup>1</sup> <sup>2</sup>*<sup>g</sup>* and *A*1g peaks [58]. While there is little variation in peak positions, small modifications in the

relative intensities between *E*<sup>1</sup> <sup>2</sup>*<sup>g</sup>* and *A*1g modes are clearly observed. Samples synthesized via first approach have a slightly higher *E*2g<sup>1</sup> intensity, with the ratio *E*<sup>1</sup> <sup>2</sup>*g*/*A*1g of ~3.2. Conversely, the *A*1g intensity is higher for the second approach, where *E*<sup>1</sup> <sup>2</sup>*g*/*A*1g ratio is of ~1.25. Moderate changes in the Raman intensity ratio are believed to indicate changes in the number of layers for exfoliated WS2 [46,66]. We have performed Raman mapping analysis at multiple spots of the same sample to confirm the uniformity of WS2 across the substrate area as shown in Figure 10b,c. It is evident that there is slight variation in the separation as well as intensity ratio of in-plane and out-plane vibrational peaks as a function of sample spot location; and a summary of the peak position, difference in Raman peak frequency and intensity ratio of the peaks is tabulated in Table 1. The shift in the frequency of 2LA (M) and A1g modes has been employed as a measure of the number of layers. However, this approach to determine the number of layers is erroneous due to the close proximity of 2LA (M) and *E*<sup>1</sup> <sup>2</sup>*<sup>g</sup>* phonon modes. It has been found that with increasing the number of layers, the intensity ratio (*I*2LA(M)/*I*A1g) decreases from 2.2 (for single layer) to 0.47 (for bulk) [64]. At the first instance, it looks the possible growth of mono- to few-layered WS2. However, in case of monolayer MoS2, such small intensity variations in *E*2g<sup>1</sup> and *A*1g peaks were correlated to the differences in electronic doping levels or strain as WS2 and RGO effectively combine to form heterostructures [67]. Generally, TMDCs are highly prone to significant number of chalcogenide vacancies [59,64,65]. Sulphur vacancies exist in WS2, thus introducing localized donor states deep inside the bandgap, which can be subsequently filled by environmental impurities to make it doped. As of now, various competitive approaches exist to find out the origin of anomaly in the interpretation of peak intensity ratio and shift in peak position of Raman modes. Nonetheless, the contribution to the intensity changes due to sulfur content in WS2 films as a dopant, and local strain in between WS2 and RGO, cannot be ignored [66].

**Figure 10.** (**a**) Raman spectrum and (**b**) Raman spectra as a function of spot on different positions (**c**) Intensity ratio and difference of the peak as function of spot (**d**) Photoluminescence (PL) of WS2 grown on GO coated Si/SiO2 substrate, and (**e**) Evolution of PL spectra due to the transformation of sp3- to *sp*2-bonded region during the reduction of GO to RGO. Reproduced with permission from [67]. Copyright 2012 Wiley.


**Table 1.** Summary of Raman studies of WS2 nanoflakes deposited on different spot of GO coated substrate.

The PL spectra, shown in Figure 10d, comprise with neutral A-exciton (*X*A) at 1.968 eV as also observed in case of WS2 growth on Si/SiO2 substrate. A striking difference observed in WS2 growth on GO-coated Si/SiO2 substrate, is the appearance of a low intensity peak at 2.03 eV (611 nm), attributable to the weak transition of carriers from disorder induced states within π-π\* band structure of RGO [67] (schematically shown in Figure 10e). Our laser excitation visible wavelength does not cover the intense emission peak around 440 nm. A comparative study in terms of PL spectral positions, FWHM and shift in peak positions are put in Table 2. Surface roughness, in particular, the sharp edges, is reported to help in catalyst-free CNT growth [68]. To ascertain the role of surface roughness on the continuous growth of WS2, we conducted CVD growth in identical conditions: first, on polished surface, marked with scratches on polished side of silicon wafer; and second, on backside, the unpolished surface of the same wafer. In the polished wafer side, WS2 debris having no shape and sizes of flakes, are observed on and around the scratch marks. SEM images (Figure 11a–d) show few vertically aligned and flower petal shaped flakes with some horizontal flakes on the scratch marked region having negligible presence of triangular as well as hexagonal WS2 crystallites.

**Table 2.** Summary of Photoluminescence (PL) studies of WS2 nanoflakes deposited on different substrates.


In contrast, we have achieved horizontally grown WS2 crystallites, majority in hexagonal shape on the unpolished surface (Figure 11e–g). In both cases, thickness of WS2 flake(s) is multilayered. It may be inferred that surface roughness does also play an important role in the WS2 growth. The optical profilometer image (surface microscopy) of the GO coated sample surface before WS2 deposition is shown in Figure 11h. It is to be noted that surface of RGO film is highly rough in nature; and this may be one of the sources of nucleation sites too, in addition to the major contribution from high density *sp*2-bonded patches, to enhance the joining of grain boundaries, resulting in continuous WS2 film growth over GO coated Si/SiO2 substrate. Therefore, it may be understood that continuous growth

is possible in GO coated substrate, in which *sp*2-patches of carbon purely act as a seed or catalyst that enhance growth origins.

**Figure 11.** FESEM micrographs of WS2 grown on scratch marks made on polished surface of Si wafer (**a**–**d**) and WS2 grown on unpolished Si wafer (**e**–**g**). In both cases, other CVD growth parameters were constant. (**h**) The optical profilometer image (surface microscopy) of the GO coated sample surface before WS2 deposition. The scale bar indicates the surface roughness of the sample close to 3.72 μm.

The present work has shown the successful growth of mono to few layer WS2 film of square-centimeter size and triangular shape on graphene-oxide-coated Si/SiO2 substrate using CVD technique. The real issue addressed in this work is the large area coverage of the WS2 continuous film; however, many of the researchers have also tried to resolve this issue by using various approaches such as atomic layer deposition, pulse laser deposition, and metal/metal oxide thin film or noble metal substrates, etc. Although these strategies yield large area continuous TMDCs, these techniques also suffer from high costs, increased complexity, and sometimes poor quality. Further, to be useful, WS2 grown over inert metal needs etching via acidic or basic (e.g., concentrated HF, KOH and NaOH) solutions, before transferring to Si or flexible substrates. Usually, the transfer processes not only end up in damaging the grown material (WS2), but produce substrate residues and negative environment impact, as well. We successfully grew the continuous WS2 film at very low cost by simply coating the Si/SiO2 substrate with graphene oxide (GO) and this technique is highly reproducible. Table 3 enlists the optical properties of WS2 nanoflakes deposited on different substrates. However, HRTEM and Raman spectroscopy confirm the high crystallinity and quality of the grown WS2 flakes. Note that the Raman peak frequencies difference (66–72 cm−1) and the PL peak position (2.00–1.92 eV) as well as intensity fluctuates over the WS2 triangles as enlisted in Table 2. The change of frequency and/or width of Raman modes of 2H-type transition-metal dichalcogenides can be used to evaluate the sample quality because they are strongly sensitive to external electrostatic doping (*A1g* mode) and lattice strain (*E*2g<sup>1</sup> mode) [69,70]. Wang et al. found that *E*2g<sup>1</sup> mode exhibits obvious red-shift when increasing strain [70]. In contrast, Chakraborty demonstrated that a softening and broadening of the *A*1*<sup>g</sup>* mode occur with electron doping, whereas the *E*2g<sup>1</sup> mode remains essentially inert [69].


**Table 3.** Summary of optical properties of WS2 nanoflakes deposited on different substrates.

• Photoconductive response of sensor prepared with GO coated Si/SiO2

To check the continuity of the grown film, I-V characteristics were monitored in the voltage range of −12 V to +12 V using Keithley SCS 4200 for variable channel length (from 0.5 to 5.0 mm). The inset of Figure 12a shows the digital photographs of the samples, where the electrical probes are set at different channel lengths. The observed I-V curves indicate that flakes grown over the GO modified substrate are interconnecting and form a continuous network. Figure 12b demonstrates the possible use of developed film as photodetector where the change in photocurrent under illumination of laser wavelength of 635 nm, keeping the laser power constant at 40 mW, is shown. The cyclic response is highly reproducible with negligible drift. The photoresponsivity and specific detectivity are found to be 0.9 <sup>μ</sup>A W−<sup>1</sup> and 0.15 × <sup>10</sup><sup>5</sup> jones.

**Figure 12.** (**a**) Current vs. voltage characteristics at variable channel length and (**b**) photoconductive response under laser wavelength 635 nm.

#### **4. Conclusions**

Efforts have been made to understand the growth mechanics of large domain size of mono-to-few layer growth of WS2. Nature of substrate is utilized to understand the intricate mechanism behind continuous film growth. Large area (domain size of around 60 μm), mono-to-few layer WS2 crystals of triangular shape, were grown on pristine Si/SiO2 and GO-coated Si/SiO2 substrates, using CVD technique with a two-zone furnace. For GO coated Si/SiO2 substrate, uniform mono-to-few layer growth of WS2 crystals is achieved, although the size of WS2 crystals is smaller than those grown on pristine Si/SiO2 substrate. In case of growth on Si/SiO2 substrate, monolayer WS2 grain boundaries are in isolation, and highly scattered on the substrate surface. Further tests were conducted to examine the role of surface roughness, as reported in the past for catalyst-free CNT growth on pristine silicon surface. We conducted CVD growth in identical conditions, one on polished surface, marked with scratches; and the other one, on the unpolished surface of the same wafer. Strikingly, WS2 growth exclusively on rough surface in both the cases, confirmed that surface roughness does play an important role. In addition, high density

*sp*2-bonded patches enhance the joining of grain boundaries, to enable continuous WS2 film growth over GO coated Si/SiO2 substrate. It is concluded that the combination of high-temperature CVD with reduced graphene oxide surface is supposed to be a favorable condition for TMDCs growth with full coverage of the substrate. The idea of GO coated Si/SiO2 substrate may be employed for mass production of identical devices, and therefore bids for commercial scale development and can be extended to other TMDCs.

**Author Contributions:** Conceptualization and Supervision, S.S.I.; formal analysis, A. and P.S.; investigation, A. and P.S.; writing—original draft preparation, S.S.I.; writing—review and editing, C.M.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors acknowledge the Inter-University Accelerator Centre (IUAC), New Delhi, India for providing the TEM facility. We are thankful to Debdulal Kabiraj and Ambuj Mishra for the TEM measurements. The author would also like to acknowledge Central Instrumentation Facility (CIF), Jamia Millia Islamia, New Delhi, India for Photoluminescence spectroscopic studies.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Photoresponsivity (Rλ)*

The ratio of photocurrent to the incident intensity of the laser wavelength λ is defined as photoresponsivity of the photodetector and it can be calculated as [71]:

$$R\_{\lambda} = \frac{\Delta I\_{\text{ph}}}{P \cdot A}.\tag{A1}$$

*Photosensitivity (S)*

It is defined as the ratio of change in current under illumination and dark current [72]:

$$\text{S (\%)} = \frac{I\_{\text{photo}} - I\_{\text{Dark}}}{I\_{\text{Dark}}} \times 100. \tag{A2}$$

*Specific Detectivity*

It represents the ability of the detector to capture the weakest signal, which can be mathematically expressed as [71]:

$$D = \frac{R\_{\lambda}\sqrt{A}}{\sqrt{2eI\_{Dark}}},\tag{A3}$$

where *A* represents the effective area of photodetector, and *e* is the electronic charge.

#### **References**


## *Article* **GFET Asymmetric Transfer Response Analysis through Access Region Resistances**

#### **Alejandro Toral-Lopez 1,\*, Enrique G. Marin 1,2, Francisco Pasadas 3, Jose Maria Gonzalez-Medina 1, Francisco G. Ruiz 1,4, David Jiménez <sup>3</sup> and Andres Godoy 1,4,\***


Received: 23 June 2019; Accepted:15 July 2019; Published: 18 July 2019

**Abstract:** Graphene-based devices are planned to augment the functionality of Si and III-V based technology in radio-frequency (RF) electronics. The expectations in designing graphene field-effect transistors (GFETs) with enhanced RF performance have attracted significant experimental efforts, mainly concentrated on achieving high mobility samples. However, little attention has been paid, so far, to the role of the access regions in these devices. Here, we analyse in detail, via numerical simulations, how the GFET transfer response is severely impacted by these regions, showing that they play a significant role in the asymmetric saturated behaviour commonly observed in GFETs. We also investigate how the modulation of the access region conductivity (i.e., by the influence of a back gate) and the presence of imperfections in the graphene layer (e.g., charge puddles) affects the transfer response. The analysis is extended to assess the application of GFETs for RF applications, by evaluating their cut-off frequency.

**Keywords:** GFET; RF; access region

#### **1. Introduction**

Two-dimensional materials (2DMs) have awakened the great interest of the nanotechnology community during the last decade [1]. Their striking physical properties, intrinsically different from their 3D counterparts, open a vast field of opportunities only partially exploited so far. Among these alternatives, 2DMs find a natural spot in electronics, where their monoatomic thickness makes them especially attractive to overcome the hurdles related to the transistor scaling-down [2].

Graphene is not only the pioneer, but also the most singular member of the 2DM family [3]. It is characterized by a gapless Dirac-cone bandstructure, where electrons and holes have symmetric dispersion relationships. The literature is abundant in Graphene Field-Effect Transistors (GFETs) [4–6], where this particular band structure is manifested in an ambipolar behaviour and a poor *I*ON/*I*OFF ratio (direct consequence of the easiness to switch the carrier transport from electrons to holes and vice versa). This issue jeopardizes the use of GFETs in digital electronics, although a successful demonstration has been achieved in [7]. In radio-frequency (RF), however, graphene has revealed itself as an interesting candidate [8], and devices with cut-off frequencies of hundreds of GHz have already been demonstrated [9,10], even reaching wafer scale integration [11], or being applied for flexible electronics [12,13]. The main strategies to boost GFETs performance have consisted of the scaling-down of the gate oxide thickness [4,14], the encapsulation in hexagonal boron nitride [15] or the improvement in the quality of the graphene-insulator stack [7,16]. In particular, clean self-aligned

fabrication, based in pre-deposited gold, has been proposed in [17]; while the self-aligned transfer of the gate stack (processed in a sacrificial substrate) has been detailed in [18].

The transfer characteristic of experimental GFETs is V-shaped, but very often shows an asymmetry with respect to the Dirac voltage [19], usually associated with different electron and hole mobilities. These mobility dissimilarities are the common path to handle the device response asymmetry, leaving out of the spot the relevance of the gate underlapped areas [15,20,21]. These access regions (intended to minimize the capacitance coupling between the gate and the source and drain) impact, however, strongly on the GFET electrical behaviour, as they constitute a noticeable resistance pathway for carrier transport. Partial attempts on the modelling of this issue have been discussed from an analytical resistance-based perspective in [20,22], but a comprehensive study of their impact in the GFET performance is still lacking [18]. In this work, we direct our attention to this asymmetric response of GFETs and, by means of detailed numerical simulations, we explain such effect studying the impact of the access regions in the transfer characteristic as well as in the RF performance of such devices.

The rest of the document is organized as follows. Section 2.1 presents the numerical model employed for this study. To check and validate it we compare, in Section 2.2, the simulated transfer response of two GFETs against the corresponding experimental measurements. Section 2.3 contains a thorough analysis of the access resistances and a discussion of its influence on the cut-off frequency, *f*T. Finally, the main conclusions are drawn in Section 3.

#### **2. Results**

#### *2.1. Device Simulation*

A schematic depiction of the physical structure of the simulated GFET is shown in Figure 1. The graphene flake is sandwiched in between a top insulator layer, with thickness *t*TOX and dielectric permittivity *ε*TOX, and an insulating substrate, with thickness *t*BOX and dielectric permittivity *ε*BOX. Both oxides are assumed thick enough as to neglect any tunnelling current through them. A four-terminal device is considered, with a front gate extending over a length *L*Chn (the device channel length), giving rise to two under-lapped regions of length *L*Acc (the access region length) that connect it with the source and drain terminals. The back gate, when considered, extends all along the structure including the channel as well as the access regions. *V*FG, *V*BG, and *V*<sup>D</sup> stand for the front gate, back gate, and drain terminal biases respectively, while the source terminal, *V*S, is assumed to be grounded. The total resistance of this structure, *R*T, can be schematically split into the series combination of three resistances corresponding to the source access region (*R*S,Acc), the channel region (*R*Chn) and the drain access region (*R*D,Acc).

To determine the *I* − *V* response of GFET devices, we have self-consistently solved the coupled Poisson, Drift-Diffusion and continuity equations [23,24]. For the device modelling, we have considered a longitudinal *x* − *y* section of the GFET, assuming invariance along the device width (*z*). The resulting 2D Poisson equation is given by:

$$\nabla \left( \varepsilon \left( \mathbf{x}, \mathcal{Y} \right) \nabla V \left( \mathbf{x}, \mathcal{Y} \right) \right) = -\rho \left( \mathbf{x}, \mathcal{Y} \right) \tag{1}$$

where *V* is the electrostatic potential; *ρ* is the net charge density in the structure, that comprises the mobile (electrons and holes) and fixed (dopants) charges; and *ε* is the dielectric permittivity.

The Drift-Diffusion transport equation is formulated in terms of the pseudo-Fermi level (*E*F) as proposed in [25]:

$$J(\mathbf{x}) = q \left[ \mu\_{\rm n} n\_{\rm 1D}(\mathbf{x}) + \mu\_{\rm P} p\_{\rm 1D}(\mathbf{x}) \right] \frac{dV\_{\rm E\_F}}{d\mathbf{x}} \tag{2}$$

where *VEF* is the potential associated with this level and *n*1D (*p*1D) is the graphene electron (hole) 1D density profile. Here, *μ*<sup>n</sup> (*μ*p) stands for the electron (hole) mobility. Due to the extreme confinement, the carriers are supposed to move only along the transport direction (*x*). *J* must comply with the continuity equation that, under steady-state conditions, is formulated as: ∇ · **J** = 0. Ohmic contacts

are assumed at the source and drain terminals, with the Fermi level at the source grounded, *E*F,S = 0, and at the drain given by *E*F,D = −*qV*DS. The equation system is then iteratively solved for each set of terminal biases, until a convergence threshold is achieved for the potential and charge concentrations.

In addition to the mobile charge and dopants in the graphene layer, we account for the existence of puddles [26,27]. Their associated charge density, *N*p, is assumed constant and added to both electron and hole charge densities [28]. In this way, puddles impact on the overall graphene layer conductivity while conserving a neutral net charge character.

**Figure 1.** Schematic of the simulated GFET and the characteristic resistances of the device. The dashed and dotted rectangles indicate the regions used for the different simulations. While the dotted rectangle only encompasses the channel region, the dashed one includes the access regions.

#### *2.2. Validation*

To assess the capability of the numerical simulator to reproduce and explain the experimental results, we have first validated it against the devices fabricated in [29,30]. Both are GFETs based on monolayer graphene embedded between a SiO2 layer, which acts as a substrate, and a Y2O3 layer, which acts as a front gate dielectric. In both cases, this Y2O3 layer is 5 nm thick while the substrate is 300 nm thick in [29], and 286 nm thick in [30]. For the device presented in [29], the distance between the source and drain contacts is 1.5 μm and the front gate length is 600 nm, while in [30] the device is 8.2 μm long and its front gate is 7 μm long. In other words, in both experimental devices the gate contact does not cover the whole region between source and drain contacts, thus creating two symmetrical under-lapped regions at both channel edges; namely, the device access regions. To reproduce the data reported in [29], the same mobility is assumed for both types of carriers, electrons and holes (*<sup>μ</sup>* = *<sup>μ</sup>*<sup>n</sup> = *<sup>μ</sup>*p) with a value of 90 cm2/Vs, and a puddle charge density of <sup>7</sup>·1011 cm−<sup>2</sup> is considered. N-type chemical doping of 10<sup>12</sup> cm−<sup>2</sup> is defined for the graphene layer. To account for the graphene-metal contact resistances, which are in series with the total resistance of the structure, *<sup>R</sup>*T, we include two additional 100 nm long N-type doped regions (5·1010 cm−2) in both source and drain ends [31]. The back gate is grounded and *V*DS is set to 0.1V. To fit the data presented in [30], the values used are *<sup>μ</sup>* =1091 cm2/Vs, *<sup>N</sup>*<sup>p</sup> = <sup>8</sup> · <sup>10</sup><sup>11</sup> cm−<sup>2</sup> and the graphene layer chemical doping is set to 10<sup>11</sup> cm−2. The back gate is also grounded and *V*DS is set to 0.05 V. The experimental and simulated transfer characteristics are shown in Figure 2a [29] and Figure 2b [30]. The simulated I-V characteristics match very accurately with the experimental results in the whole range of biases and are able to catch the transfer response of the electron and hole branches, especially in Figure 2b.

**Figure 2.** Comparison between the simulation results and the experimental data extracted from [29] (**a**) and [30] (**b**).

#### *2.3. Access Region Analysis*

As mentioned in Section 1, the existence of access regions and puddles is a very common scenario in the experimental realization of GFETs due to the difficulties to precisely control the fabrication process in this early stage of the technology. They modify the behaviour of the transistors, in many cases determining their performance, and therefore deserving a particular attention that is usually obliterated. Hence, once the numerical simulator has been validated, we now proceed to analyse the effect of the access regions.

#### 2.3.1. Including the Access Regions

To begin with, we have considered a test structure where the front gate covers the whole device length (i.e., suppressing the access regions) and compared the results with those obtained later when access regions are included. These scenarios are illustrated in Figure 1 by the dotted and dashed frames respectively. The material stack comprises a monolayer graphene sandwiched between a 3 nm thick HfO2 layer (front gate insulator) and a 27 nm thick SiO2 layer (back gate insulator). The front gate, which determines the channel length (*L*Chn), is 100 nm long and both access regions are 35 nm long (*L*Acc). Electron and hole mobilities are equal (*μ* = 1500 cm2/Vs) and no chemical doping or puddle charge density is considered in the graphene layer.

The transfer characteristic of the device without access regions is depicted in Figure 3a for different values of *V*DS. As can be observed, the device exhibits the ambipolar V-shaped *I* − *V* response of an ideal GFET. The minimum of the *I* − *V* curve defines the Dirac voltage (*V*Dirac) that is shifted to larger *V*FG when *V*DS increases. The behaviour is perfectly symmetric with respect to *V*Dirac, reflecting the symmetry between electron and hole properties.

Next, the GFET including the access regions is investigated. The resulting transfer characteristic is shown in Figure 3b. Comparing Figure 3b and Figure 3a, a marked variation of the GFET response is observed. First, there is a notable decrease in the values of *I*DS, around a factor ×100. Second, the transfer characteristic shows a saturation trend for high |*V*FG| which resembles much better the experimental response. Third, and more important, the *I* − *V* characteristic is no longer symmetric with respect to *V*Dirac, though the mobility is identical for both kinds of carriers.

**Figure 3.** *I*DS − *V*FG curves of the device without (**a**) and with (**b**) access regions.

To provide insights into these changes, the resistance of the different regions of the device are calculated. Figure 4 shows their values for *V*DS = −0.1 V and *V*DS = −0.2 V. Mirror symmetric behaviour is observed for positive *V*DS. The access region resistances, *R*S,Acc and *R*D,Acc, show values comparable with the channel resistance, *R*Chn. At the Dirac voltage, where the channel resistivity is the highest, *R*Chn commands the series association, but still the access regions have a noticeable contribution. For |*V*FG − *V*Dirac| > 0.1 V the total resistance is mainly determined by *R*S,Acc and *R*D,Acc. Consequently, the total resistance (*R*T) is not controlled just by the channel conductivity and, therefore, by the gate terminal. The weak dependence of *R*S,Acc and *R*D,Acc on *V*FG is reflected in the *I*DS trend to saturation. As the values of *R*S,Acc and *R*D,Acc are higher than the channel resistance, a larger fraction of *V*DS drops in the access regions. This fact reduces the potential at the channel edges with respect to the no-access-regions scenario, reducing the output current. In addition, the *R*Acc − *V*FG dependence is not symmetric, so neither are the access region potential drops, resulting into a non-symmetric reduction of the output current, that is, an asymmetric *I*DS − *V*FG curve shown in Figure 3b. This lack of equivalence between the source and drain access regions is explored in detail in the following section.

**Figure 4.** Resistance of the three device regions (channel, source and drain access regions) compared with the total resistance as a function of the gate potential, for two *V*DS values: −0.1 V (solid) and −0.2 V (dashed).

#### 2.3.2. Gate Misalignment

In the previous section, we assumed that the gate is perfectly aligned in the middle of the channel leading to identical source and drain access regions (*L*<sup>S</sup> = *L*<sup>D</sup> = *L*Acc) at both ends. A more realistic scenario should consider the impact of having non-equal *L*<sup>S</sup> and *L*D, enabling us to test the non-equivalent role of *R*S,Acc and *R*D,Acc on the GFET response. For this purpose, we have analysed GFETs where the top gate contact is not placed in the centre of the structure, resulting in access regions of different length. In particular, we have kept *L*<sup>S</sup> (or *L*D) equal to 35 nm while *L*<sup>D</sup> (or *L*S) is modified. Specifically, we considered four scenarios: (i) short source, (ii) short drain, (iii) long source and (iv) long drain. The length of the short and long regions is set to 17.5 nm and 70 nm, respectively. The *I*DS − *V*FG curves, along with the resistances *R*S,Acc, *R*D,Acc and *R*Chn obtained in each case, are depicted in Figure 5.

**Figure 5.** Transfer response (**a**,**b**) and structure resistances (**c**,**d**) as a function of the gate bias. These results are obtained reducing the length of either the source (**a**,**c**, solid lines) or drain access region (**b**,**d**, dashed lines) down to 17.5 nm, and increasing the length of either the source (**a**,**c**, solid lines) or the drain access region (**b**,**d**, dashed lines) up to 70 nm.

As expected, there are significant differences between devices. Shortening either the source or the drain access regions results in a higher output current (Figure 5a) and reduces both its saturation and its asymmetry with respect to the elongated scenario (Figure 5b). When comparing the shorter regions (Figure 5a) it is clearly observable that the *L*<sup>S</sup> = 17.5 nm device (solid lines) has a more symmetric response than the *L*<sup>D</sup> = 17.5 nm (dashed lines). This is more evident for *V*DS = 0.1 V and emphasizes the role of the source access region with respect to the drain access region. An equivalent conclusion can be achieved from the elongated devices (Figure 5b). The longer *L*<sup>S</sup> results in an increased asymmetry between both branches. These results can be explained by analysing the resistances of the structure. Figure 5c,d show *R*S,Acc, *R*D,Acc and *R*Chn as a function of *V*FG for *V*DS = 0.1 V. When any access region is shortened (Figure 5c), its resistance is similar or lower than the channel resistance regardless *V*FG. The longer region resistance controls the total current (except for *V*FG close to zero). When one of the regions is enlarged this effect is emphasized. The transfer responses in Figure 5b are clearly saturated due to the dominant role in the total conductivity of the longer access region.

#### 2.3.3. Impact of Electrostatic Doping and Puddles

To reduce the impact of the access regions in the overall device performance, it is possible to increase their conductivity by means of an electrostatic doping using the back-gate terminal. In the following we analyse how the back gate influences the GFET behaviour. Figure 6 shows the transfer characteristic for three different values of *V*BG: −1 V, 0 V and 1 V (solid lines). For *V*BG = 0 V the results are quite similar to the scenario without back gate. In the other two cases, depending on the polarity of *V*BG, electrons or holes are accumulated in the graphene layer. As a result, the P-type (N-type) branch is enhanced for *V*BG = −1V(*V*BG = 1 V), regardless the value of *V*DS. As in the previous scenario, the origin of this behaviour can be traced back to the resistance associated with the access regions.

**Figure 6.** *I*DS − *V*FG characteristics of the complete structure when three different back gate potentials are used (−1V(**a**), 0 V (**b**) and 1 V (**c**)). Solid lines correspond to the device without puddles and dashed lines to the device with *N*<sup>p</sup> = 1012 cm<sup>−</sup>2.

Figure 7 depicts the device resistances for different *V*BG and *V*DS = −0.1 V (without puddles, solid lines). For |*V*BG| = 1 V the total resistance near the Dirac voltage is dominated by *R*Chn. When *V*FG is increased above *V*Dirac, the symmetry of *R*Chn is kept since it is mostly controlled by the front gate, while the asymmetry of *R*S,Acc and *R*D,Acc is exacerbated due to the electrostatic doping, giving rise to the large asymmetry observed in the transfer response, in Figure 6. In particular, the asymmetric step-like dependence of the access resistances on *V*FG (for *V*BG = 0 V) is the result of the electrostatic competition between the front and back gates to control the access regions closer to the channel. When *V*FG and *V*BG have the same polarity, they add their electric forces to increase the carrier density in the aforementioned zones, increasing the conductivity and therefore lowering the whole access resistance. However, if *V*FG is opposite to *V*BG, both gates compete to accumulate different types of charges, resulting in a depleted region close to the channel edges that decreases the conductivity and increases the overall access region resistances. An equivalent conclusion was achieved in [26] where a strong modulation of the total resistance by two additional gates is observed, as in Figure 7.

An additional aspect that cannot be overlooked is the effect of the presence of puddles in the graphene layer [27,32]. To shed light on this issue Figure 6 includes the *I*DS − *V*FG response when a puddle charge density of *N*<sup>p</sup> = 1012 cm−<sup>2</sup> is considered (dashed lines). Two major changes are observed after including the puddles: (i) the total current is increased, and (ii) the asymmetry is clearly reduced. These changes derive from the equal contribution of puddles to the conductivity of both electrons and holes, and explain why the *I* − *V* curves of some experimental devices are reasonably symmetric close to the Dirac voltage, where the conductivity of puddles is dominant. In this situation, the conductivity of the whole graphene layer is increased for electrons and holes, in contrast with the electrostatic doping generated by the back gate. This non-selective improvement of the conductivity is translated into the resistances of the device: Figure <sup>7</sup> includes the *<sup>R</sup>* − *<sup>V</sup>*FG relation for *<sup>N</sup>*<sup>p</sup> = 1012 cm−<sup>2</sup> (dashed lines). The step-like behaviour of *R*S,Acc and *R*D,Acc is softened when the puddles are included, resembling the *V*BG = 0 V case.

**Figure 7.** Total (**a**), channel (**b**), source (**c**) and drain (**d**) resistances for different back gate biases and *V*DS = −0.1 V. Solid lines (referred to the left axis) correspond to the no puddles scenario while dashed lines (referred to the right axis) depict the values obtained when a puddle concentration of *N*<sup>p</sup> = 10<sup>12</sup> cm−<sup>2</sup> is considered.

#### 2.3.4. RF Performance

Finally, to determine the impact of the access regions in the RF performance, we evaluate the cut-off frequency, *f*T, as a RF figure of merit (FoM). The value of *f*<sup>T</sup> is calculated as in [33,34]:

$$f\_{\rm T} = \frac{1}{2\pi} \frac{\rm g\_{\rm m}}{\rm C\_{\rm fg}} \tag{3}$$

where *g*<sup>m</sup> is the transconductance and *C*fg the front gate capacitance.

Figure 8 shows *f*<sup>T</sup> as a function of *V*FG under two scenarios: no puddles (solid lines) and *N*<sup>p</sup> = 1012cm−<sup>2</sup> (dash-dotted lines). To assess the impact of the access regions, the performance of the intrinsic device (structure indicated by the dotted rectangle in Figure 1) is depicted too (dashed lines). In addition, to evaluate the magnitude of the calculated values, the experimental measurements of *f*<sup>T</sup> reported in [35] and [36] are indicated by the arrows on the right side axis of Figure 8. Despite the device structure and the bias conditions are different, the channel lengths of these experimental devices are similar to the ones simulated here (144 nm [35] and 140 nm [36]), and therefore constitute a good reference. Importantly, a de-embedding procedure was carried out for the RF measurements of these experimental devices by using specific "short" and "open" structures with identical layouts in order to remove the effects of the parasitics associated with the pads and connections, but not the contact and access region resistances.

**Figure 8.** *f*<sup>T</sup> of the back-gated device with access regions under two scenarios: no puddles (solid lines) and *N*<sup>p</sup> = 10<sup>12</sup> cm−<sup>2</sup> (dash-dotted lines). The values obtained for the intrinsic device are depicted by the purple dashed line. The arrows labelled by marks on the right side axis indicate the values of *f*<sup>T</sup> extracted from [35] (circle) and [36] (square and triangle). The yellow line indicates the physical limit for graphene *<sup>v</sup>*F/2*πL*, determined by the transit time *<sup>L</sup>*/*v*F, with the Fermi velocity *<sup>v</sup>*<sup>F</sup> <sup>≈</sup> <sup>10</sup><sup>8</sup> cm/s and *L* =100 nm (squares).

Including the access regions results in a quite different response compared with the intrinsic device, as the associated parasitic resistances provoke a bias dependent decay of *f*T. Considering the scenario without puddles, when the back gate is properly biased, *f*<sup>T</sup> is considerably improved. If we analyse Figure 8 in combination with Figure 7, those combinations of *V*FG, *V*BG for which the *R*<sup>S</sup> − *V*FG (*R*<sup>D</sup> − *V*FG) curve shows its minimum values, are those for which *f*<sup>T</sup> shows a greater improvement. When *R*<sup>S</sup> (*R*D) is higher, *f*<sup>T</sup> is spoiled with respect to the *V*BG = 0 V case. This relation between the access region conductivity and the improvement of the RF performance was experimentally observed in [21] where a higher *f*<sup>T</sup> was demonstrated when a GFET with two additional electrodes was properly biased to control such conductivity. When puddles are included, the channel conductivity increases, what reduces the control of the back-gate bias, and simultaneously results in a more symmetric *f*<sup>T</sup> − *V*FG dependence.

#### **3. Conclusions**

GFETs have been thoroughly studied in order to assess the impact of the access regions in the device performance. The validation of our approach against two experimental devices spotlights the importance of these regions as well as the presence of puddles to reproduce the state-of-the-art technology. When the access regions are considered, the transfer response reveals a lower, saturated and asymmetric *I*DS − *V*FG characteristic that is not observed in their absence. To explore the impact of a variable conductivity of these regions we have included a back gate in the structure able to introduce an electrostatic doping. The back gate increases the output current as well as the asymmetry of the transfer characteristic. The latter effect is explained in terms of the competition of the back and front gates that results in a depletion of the amount of carriers close to the channel edges when both biases have an opposite polarity. The influence of puddles is also theoretically investigated, observing that they reduce the asymmetry of *I*DS − *V*FG.

The analysis of the impact of the access regions and puddles have been extended to the prediction of the cut-off frequency to assess the properties of GFETs for potential RF applications. Our results reveal an important degradation of the *f*<sup>T</sup> − *V*FG relation due to access regions. The application of an appropriate back gate bias can tune the access region conductivity generating a remarkable improvement in the RF performance. The presence of puddles also mitigates this degradation, but neglects the possibility of tuning the access regions conductivity.

**Author Contributions:** A.T.-L., E.G.M., F.G.R. and A.G. conceived the work. A.T.-L., F.P., and J.M.G.-M. performed the numerical simulations under the supervision of E.G.M., F.G.R., D.J. and A.G. All authors analysed the results, contributed to the discussion and wrote the manuscript.

**Funding:** This research was founded by Spanish government grant numbers TEC2017-89955-P (MINECO/AEI/FEDER, UE), TEC2015-67462-C2-1-R (MINECO), IJCI-2017-32297 (MINECO/AEI), FPU16/04043 and FPU14/02579, and the European Union's Horizon 2020 Research and Innovation Program under Grant GrapheneCore2 785219.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

2DM Two-dimensional material MOSFET Metal-Oxide-Semiconductor Field-Effect Transistor GFET Graphene Field-Effect Transistors

RF Radio-Frequency

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Calculations of Some Doping Nanostructurations and Patterns Improving the Functionality of High-Temperature Superconductors for Bolometer Device Applications**

#### **Jose C. Verde 1, Alberto S. Viz 1,2, Martín M. Botana 1,2, Carlos Montero-Orille 2,3 and Manuel V. Ramallo 1,2,\***


Received: 26 November 2019; Accepted: 30 December 2019; Published: 3 January 2020

**Abstract:** We calculate the effects of doping nanostructuration and the patterning of thin films of high-temperature superconductors (HTS) with the aim of optimizing their functionality as sensing materials for resistive transition-edge bolometer devices (TES). We focus, in particular, on spatial variations of the carrier doping into the CuO2 layers due to oxygen off-stoichiometry, (that induce, in turn, critical temperature variations) and explore following two major cases of such structurations: First, the random nanoscale disorder intrinsically associated to doping levels that do not maximize the superconducting critical temperature; our studies suggest that this first simple structuration already improves some of the bolometric operational parameters with respect to the conventional, nonstructured HTS materials used until now. Secondly, we consider the imposition of regular arrangements of zones with different nominal doping levels (patterning); we find that such regular patterns may improve the bolometer performance even further. We find one design that improves, with respect to nonstructured HTS materials, both the saturation power and the operating temperature width by more than one order of magnitude. It also almost doubles the response of the sensor to radiation.

**Keywords:** superconducting devices; photodetectors; nanostructured materials, nanostructured and microstructured superconductors; high temperature superconductors; bolometers

#### **1. Introduction**

Bolometers are radiation sensors that detect incident energy via the increase of the temperature *T* caused by the absorption of incoming photons [1–16]. Bolometers are often used, e.g., for thermal infrared cameras (see, e.g., [1–10]), mm-wave sensing [1–6,8–13], space-based [12,13], laboratory far-infrared spectroscopy [13], etc. [1–16]. Superconductors are among the best candidate materials for bolometers, due to their extreme sensitivity to *T* near the superconducting transition, measurable for instance through the sharp variations of the electrical resistance *R* (resistive transition-edge bolometer—TES). For resistive TES bolometers, a key figure for performance is the so-called "temperature coefficient of resistance" (TCR), given by [1–19]:

$$\text{TCR} = \left| \frac{1}{R} \frac{\text{d}R}{\text{d}T} \right|. \tag{1}$$

High bolometric sensitivity requires a large value of TCR. For instance, structures of vanadium oxides <sup>V</sup>*x*O*y*, commonly used in semiconductor-based bolometers, present TCR <sup>∼</sup> 0.025 <sup>K</sup>−<sup>1</sup> [20]. Much larger TCR may be achieved with superconductor materials kept at base temperatures coincident with their normal–superconducting transition, *Tc*. This is the case mainly when using conventional low-temperature superconductors with *Tc* ≤ 1 K (the so-called low-*Tc* TES bolometers), that achieve TCR ∼ 1000 K−<sup>1</sup> or even more [17–19], making them a technology of choice for detecting the most faint radiations, as the cosmic infrared background or in quantum entanglement and cryptography applications [17–19]. Note that for these measurements the very low temperature required to operate the low-*Tc* TES is often not seen as a major problem, because cryogenizing the sensor below a few Kelvin is required anyway in order to minimize the thermal noise coming from the bolometer itself. However, the requirement of a highly-stabilized liquid-helium-based cryogenics is a serious difficulty for adoption of low-*Tc* TES in other applications.

After the discovery of high-*Tc* cuprate superconductors (HTS), various authors have explored their use for resistive bolometers with simpler liquid-nitrogen-based cryogenics (the so-called resistive HTS TES bolometers [1–16]). The compound YBa2Cu3O*<sup>δ</sup>* (YBCO) is the HTS material usually considered for this application, usually with maximum-*Tc* doping, i.e., stoichiometry *δ* 6.93. Such YBCO thin films provide TCR <sup>∼</sup> 1.5 <sup>K</sup><sup>−</sup>1, low noise at an operational temperature *Tc* <sup>∼</sup> 90 K, and also favorable values for the rest of parameters contributing to good bolometric performance (thermal conductivity, infrared absorbance, response time, etc.) [1–16].

Besides the cryogenics, the other difference with respect to low-*Tc* TES is that in actual implementations [1–16] the resistive HTS TES operate under current bias and usually in the ohmic regime (instead of non-ohmic resistance and the voltage bias employed in low-*Tc* TES to avoid thermal runaways [17–19]).

However, the HTS TES until now proposed still share some of the significant shortcomings of low-*Tc* TES: First, thermal stability of the cryogenic bath is still challenging (liquid-nitrogen systems are simpler but tend to thermally oscillate more than those based on liquid helium). Secondly, both types of TES have useful TCR only at the superconducting transition, corresponding to operational temperature intervals Δ*T* of just ∼0.1 K or less for low-*Tc* TES, and ∼1 K for the resistive HTS TES proposed until now [1–19].

The HTS TES systems proposed until today are homogeneous in nominal composition and critical temperature [1–16]. However, in the recent years different novel techniques have been developed to impose regular patterns on HTS thin films, creating custom designs, down to the micro- and the nano-scales [21–30]. This allows custom-engineering regular variations of the critical temperature over the film surface. Realization of these regular and controlled patterning has been experimentally achieved using, e.g., local ferroelectric field-effect [21], nanodeposition [22], focused ion beam [23], etc. In fact, the nanostructuring of HTS has become the specific subject of recent conferences [24–27] and networks [28,29] funded by the European Union.

However, the use of nanostructured films for optimizing HTS TES has been considered only very marginally up to now, the only precedent to our knowledge being Reference [6] by Oktem et al., who considered films with random distributions of nonsuperconducting incrustations producing limited increases of Δ*T* up to only ∼2 K (and also small, and not always favorable, TCR variations).

Our aim in the present work is to propose that certain custom nanostructurating and patterning of HTS materials may improve their functionality for resistive HTS TES sensors. In particular, we calculate the case of nanostructuring and patterning of the local carrier doping level *p* (the number of carriers per CuO2 unit cell) in the prototypical HTS compound YBCO via local variations of oxygen stoichiometry (as realizable, e.g., via local desoxigenation, ion bombardement with different masks, etc.). Our main objective will be to obtain an increase of the operational temperature interval, Δ*T*, in which (*i*) the TCR is large and (*ii*) *R* is linear with *T* (i.e., d*R*/d*T* constant with *T*, that is another desirable feature that simplifies both the electronic control of the bolometer and the required stability of the cryogenic setup). Accompanying this Δ*T* increase we will also obtain improvements of other bolometric characteristics, such as the saturation power and in some cases the TCR itself.

We organize our studies of the structured HTS materials in two parts: First we study the simplest case of carrier doping nanostructuring, namely the random nanoscale structuration that appears by just using oxygen stoichiometries that do not maximize *Tc*. We present our methods for those randomly structured HTS in Section 2; these consists of finite-element computations (and also, to confirm their validity, analytical estimates using effective-medium approximations [31]) that we apply to calculate the performance of the material in various example bolometer device implementations. The results following these methods in random nanoscale structurations are presented in Section 4. These results indicate that this first simple structuration may already improve some of the bolometric parameters with respect to conventional, nonstructured HTS materials

The rest of the paper considers structurations that include not only the unavoidable random disorder but also the additional imposition of custom regular arrangements of zones with different nominal doping levels (patterning), studying different examples aimed to progressively improve the bolometric performance. The additional methods needed to calculate this added patterning are presented in Section 3, and the results are discussed in Sections 5 to 7 for various custom pattern designs, each of them improving the previous one. The most optimized pattern design (Section 7) is a four-step discretized exponential-like dependence of nominal doping with the longitudinal position. This arrangement should be also the easier one to fabricate. With respect to conventional nonstructured HTS TES materials, it improves by more than one order of magnitude the Δ*T* and the saturation power, and it also doubles the TCR sensitivity.

#### **2. Methods for Structured Nonpatterned Resistive HTS TES**

Our methodology consists of computing the electrical resistance *R* versus temperature *T* of each of the structured HTS materials considered by us, and then using such *R*(*T*) to calculate the corresponding performance for bolometric operation.

For completeness we will, in fact, consider various example bolometer-device designs, including simple square-shaped sensors such as those in [7,16] (that may be a micrometric size as appropriate for building megapixel cameras; we shall also consider two different substrates for completeness) and also larger sensors for millimeter-wavelength sensing using a meander geometry (as those built, always with nonstructured HTS, in [1–4,6,8–10,12,13]). While naturally we could not calculate in this work the whole range of possible device designs for a bolometer, our results in these example implementations show that our proposed nano optimizations of the materials should lead to improvements in at least some popular types of resistive HTS TES device designs.

Also, for our *R*(*T*) calculations we shall use two alternative calculations, so as to be confident about the validity of the results: Finite-element computations first, and then effective-medium formulae (both paths have been successful in other studies of structured HTS [30–35], and we shall also include some confirming example comparisons with real data). Let us provide the details of all such procedures in the following sections.

#### *2.1. Main Operational Parameters for Resistive HTS TES Devices*

We consider in this work, a HTS TES of resistive type, i.e., the temperature increase is sensed through the measurement of the electrical resistance, as in [1–16]. Contrarily to the most common case of low-*Tc* TES, the measurement is in current-bias (*I*-bias) configuration in all experimentally implemented resistive HTS TES published to our knowledge (the experimental difficulties for voltage-bias (*V*-bias) sensing in HTS TES were explained, e.g., by Khrebtov et al. [15]). Also the *I* value usually employed [1,7,9,10] is sufficiently small as to correspond to the ohmic regime (*R* = *V*/*I* constant with *I*) in all the operational *T*-range (see also below; this is also in contrast to low-*Tc* TES). As already mentioned in the Introduction, for ohmic resistive TES a main parameter of merit is the TCR, that may be also expressed as

$$\text{TCR} = \frac{R(T^+) - R(T^-)}{\Delta T \, R(T^-)},\tag{2}$$

where *T*<sup>−</sup> is the base operation temperature (the one in absence of radiation), *T*<sup>+</sup> is the maximum temperature up to which the ohmic *R*(*T*) maintains the strong and constant slope with temperature, and

$$
\Delta T = T^+ - T^- \tag{3}
$$

will be henceforth called the operational temperature interval.

The other important parameter is *P*max, the maximum power measurable without saturation. In a *I*-bias resistive TES, it is possible to obtain *P*max at good aproximation [18,19] by just applying the heat flow equilibrium condition at saturation:

$$P^{\text{max}} + I^2 R(T^+) = G \Delta T,\tag{4}$$

where *I*2*R*(*T*+) is the heat rate due to the Joule effect, *G*Δ*T* is the power dissipated towards the cryobath, and *G* is the thermal conductance between the film and the bath.

We shall consider in this work three example resistive HTS TES device designs, to probe the effects of our proposed material optimizations in them. In particular, we consider two cases of microsensor bolometers, plus one case adapted to millimeter-wavelength sensing that we specify in the following section.

#### 2.1.1. Microsensor Device Design

The first resistive HTS TES device design consists of depositing a thin layer of YBCO HTS material over a substrate, the area of the HTS and the substrate being micrometric. In particular, we consider the convenient area (6 μm)2, that makes possible building a 1 megapixel array of sensors in ∼1 cm2. Each substrate is in direct contact with a cryogenic liquid-nitrogen bath and we shall consider two possible substrate compositions: SrTiO3 (STO, most appropriate to grow HTS films) and a CMOS-type substrate of interest for technological integration (note that HTS TES over CMOS substrates, in particular silicon/Yttria-stabilized zirconia (YSZ)/zirconia, have been already fabricated [10], using nonstructured YBCO). The experimental value of *G* for both types of substrate can be obtained from [10,36]. We consider YBCO thickness 100 nm and substrate thickness 1 mm. Also, we consider a bias current of *I* = 6 μA, that corresponds to a current density *j* = 10<sup>3</sup> A/cm2, a value used in [13,15] and that corresponds to the ohmic range in all the *T*-range of operation [13,15,37].

#### 2.1.2. Millimeter-Wave Sensor Device Design

The second device design we shall consider corresponds to the one most employed by experimentalists having produced resistive HTS TES [1–6,8–10,12,13]. It corresponds to a larger design using, as substrate, a suspended membrane of millimetric surface and CMOS-type composition; on top of such thin (micrometric thickness) membrane substrate a single meander of YBCO material is deposited. The larger area precludes building small megapixel sensors, but this is not important, e.g., for sensing millimeter wavelengths (that could not be constrained in smaller pixel sizes anyway, and that are among the main applications of bolometers). The meander geometry allows instead the ability to optimize the so-called "static voltage responsivity", *S*V, an important parameter defined by:

$$S\_{\rm V} = \frac{\varepsilon L\_0}{I(1 - L\_0)'} \tag{5}$$

where *ε* is the absorbance of the sample, *I* is the bias current, and

$$L\_0 = \frac{\text{TCR } I^2 \text{ R}}{\text{G}} \tag{6}$$

is the so-called loop gain coefficient, which is a relative measure of the positive electrothermal feedback of the device [15] and *G* is again the thermal conductance towards the bath (whose experimental value may be found, e.g., in [13]); we consider (3 mm)<sup>2</sup> membranes for that evaluation. For a stable operation, the loop gain coefficient *L*<sup>0</sup> should be smaller than 1, the value *L*<sup>0</sup> = 0.3 being usually taken as optimal. Therefore, the maximum static responsivity is obtained by tuning the geometry of the meander so to tune the *I*2*R* contribution to *L*<sup>0</sup> in Equation (6) (Let us also note here that for *I*-bias the *L*<sup>0</sup> may remain constant, and *L*<sup>0</sup> = 0.3, only if *R* is linear with *T*, so that maximizing Δ*T* is also interesting in this respect).

In our calculations, we will use for the meander section the same size 6 μm × 100 nm as previously for μm-sensors. We then choose the meander length for each sensing material so that always *L*<sup>0</sup> ∼ 0.3. We also consider the same bias current and current density, I = 6 μm and *j* = 103 A/cm2, than for μm-sensors (again corresponding to the ohmic range [13,15,37] and comparable to values used in experimental meander resistive HTS TES [13,15]). These choices not only are realistic but also they will allow us to use the same computer calculations of *R*(*T*) for both μm- and mm-device designs (as only a geometric correction prefactor is needed to change *R* from one design to the other) and the same pair of bias current and current density, *I* = 6 μA and *j* = 10<sup>3</sup> A/cm<sup>2</sup> .

#### *2.2. R*(*T*) *in the Normal State of Nonstructured HTS*

In HTS materials, the ohmic resistance versus temperature, *R*(*T*), markedly varies with the doping level *p* (number of carriers per CuO2 unit cell that for instance in YBCO may be changed through the oxygen content). This is true both for the value of the critical temperature *Tc*(*p*) below which, the superconductivity transition occurs, and for the *R*(*T*) magnitude and *T*-dependence in the normal state *T* > *Tc*(*p*). The *R*(*T*)-versus-*p* phase diagram has been extensively studied in many works such as the review [38] (see also, e.g., [39,40]). Here, let us recall that the superconducting critical temperature is maximum at *p* ∼ 0.155 (separating the so-called underdoped *p* < 0.155 and overdoped *p* > 0.155 compositions). Above *Tc*(*p*), the material presents a normal-state background electrical resistivity, *ρ*b(*T*, *p*), that is linear on *T* above a certain so-called pseudogap temperature *T*∗, and is pseudoparabolic semiconducting-like [38] for *Tc* < *T* < *T*∗. In YBCO, it is *T*∗(*K*) ≈ 270 − 3000(*p* − 0.1) [38], so that for *p* - 0.16, it is *T*<sup>∗</sup> < *Tc* and the semiconducting-like region disappears. Instead of these rapid crude approximations, we will use in our analysis, all through the present work, the detailed quantitative results for *Tc*(*p*), *T*∗(*p*), and *ρ*b(*T*, *p*) given in reference [38] for YBCO.

Near *Tc*(*p*), obviously *R*(*T*) undergoes the superconducting transition towards *R*(*T*) = 0. This transition is not fully sharp, instead, a sizable rounding of *R*(*T*) occurs in the vicinity of *Tc*(*p*). This rounding is known to have two contributions: critical fluctuations and doping inhomogeneities, that we describe in the following subsections.

#### *2.3. Rounding of R*(*T*) *Near the Superconducting Transition Due to Critical Phenomena*

The critical fluctuations around the transition play an important role in HTS and have been studied in detail, e.g., in [32,41–48]. The effects of critical fluctuations in the resistance curves are commonly summarized via the so-called paraconductivity, Δ*σ*, defined as the additional contribution to the electrical conductivity due to fluctuations: In particular, the total conductivity *σ*(*T*) near the transition becomes

$$
\sigma = \Lambda \sigma + 1/\rho\_{\flat\flat} \tag{7}
$$

where *ρ<sup>b</sup>* is the normal-state background resistivity (see previous subsection). Because Δ*σ* follows critical-divergence laws near the transition, its effect far from *Tc* (for *T* -1.7*Tc*) is totally negligible. Closer to *Tc*, however, Δ*σ* becomes progressively important and two *T* ranges may be distinguished. For 1.01*Tc T* 1.7*Tc*, i.e., the so-called Gaussian fluctuations region, Δ*σ* is well described by the Lawrence–Doniach paraconductivity equation for layered superconductors: [41–45,49]

$$
\Delta \sigma = \frac{\varepsilon^2}{16\hbar d} \left\{ \frac{1}{t} \left( 1 + \frac{B\_{\rm LD}}{t} \right)^{-1/2} - \frac{2}{c} + \frac{t + B\_{\rm LD}/2}{c^2} \right\}, \tag{8}
$$

where *e* is the electron charge, *h*¯ is the reduced Planck constant, *t* = ln(*T*/*Tc*) is the reduced temperature, *B*LD = (2*ξc*(0)/*d*)<sup>2</sup> is the Lawrence–Doniach [49] layering parameter, *ξc*(0) is the Ginzburg–Landau coherence length amplitude in the out-of-plane direction (1.1 Å in YBCO [44,45]), *d* is the superconducting layer periodicity length (5.85 Å in YBCO [44,45]) and *c* is a high-temperature cutoff constant 0.7 [43,44,50,51].

Closer to *Tc*, for *TBKT T* 1.01*Tc* we find the strong phase fluctuation regime, dominated by the Berezinskii–Kosterkitz–Thouless (BKT) transition temperature *TBKT* (∼ *Tc* − 2*K* in YBCO) [42,52]. In this regime, the paraconductivity can be obtained using the equation: [33–35,47,53]

$$
\Delta \sigma = A\_{\rm BKT} \exp \left[ 4 \sqrt{\frac{T\_{\rm c} - T\_{BKT}}{T - T\_{BKT}}} \right], \tag{9}
$$

where *A*BKT is a constant, obtainable by requiring continuity of Equations (8) and (9) at the intersection of the Gaussian and BKT regimes, i.e., at *T* = 1.01*Tc*.

#### *2.4. R*(*T*) *Transition Rounding Due to Intrinsic Structuration of the Carrier Doping Level; Nominal vs. Local Doping*

As it has been explicitly demonstrated in various relatively recent experimental and theoretical works [32,46,47,54,55], an additional (and crucial for some doping levels) ingredient to understand the phenomenology of the resistive transition in HTS is to take into account the random *Tc*-inhomogeneities associated with the intrinsic disorder of the doping level. This intrinsic structuration is due to the fact that HTS compounds are non-stoichiometric, and therefore each dopant ion has at its disposal various lattice positions to occupy. For concreteness, we focus our present article in the case of the YBa2Cu3O*<sup>δ</sup>* superconductor with oxygen as a dopant ion. Experimental measurement indicates that a typical size of each local inhomogeneity is about (30 nm)2 for HTS [32,54,55]. This produces, therefore, a certain randomness in the doping at the local scale, unavoidably present in even the more carefully grown HTS samples. A relatively easy geometrical calculation [32,54,55] reveals that this intrinsic structuration shall produce a Gaussian distribution of local dopant levels, as

$$
\omega(p,\overline{p}) = \frac{2\sqrt{\ln 2}}{\sqrt{\pi}\Delta p} \exp\left[-\left(\frac{p-\overline{p}}{\Delta p/(2\sqrt{\ln 2})}\right)^2\right],\tag{10}
$$

where *ω*(*p*, *p*) is the fraction distribution of local doping levels, *p*, for a film with average doping level *p* (henceforth called nominal doping), and Δ*p* is the FWHM of the Gaussian distribution. This Δ*p* may be obtained, in turn, on the grounds of coarse-graining averages (see, e.g., Equation (6) of References [32,54]) and for YBCO it is Δ*p* ∼ 0.006 (with a small dependence on *p* that may be considered in excellent approximation linear Δ*p* = 0.0032 + 0.0189*p*) [32,54].

Due to the *Tc*(*p*) dependence in HTS, the above distribution of local *p* values leads, in turn, to a corresponding distribution of local critical temperatures around the nominal value *Tc* = *Tc*(*p*). The corresponding full width at half maximum (FWHM) for such intrinsic *Tc* structuration has been also considered, e.g., in [32,54]. Not surprisingly, it becomes quite negligible (∼1 K in YBCO) for the nominal dopings *p* 0.155 that maximize *Tc* (and that has been used up to now for HTS TES; *p* 0.155 corresponds to YBa2Cu3O6.93 stoichiometry at which *Tc*(*p*) is maximum and *dTc*/*dp* ∼ 0). However,

for other dopings the situation may become very different and the *Tc* distribution can reach FWHM values as large as, e.g., ∼5 K for *p* = 0.13, significantly influencing the *R*(*T*) roundings [32,54,55].

#### *2.5. Obtainment of the R*(*T*) *Curve of Nonpatterned HTS TES Using Finite-Element Computations*

To calculate the resistance transition curves, *R*(*T*), of the *Tc*-structured HTS material, we have used software (TOSERIS, available by request to authors) that numerically solves the electrical mesh-current matrix equations of a film modeled as a 200 × 200 square lattice of monodomains, where each domain *i* may have its own doping *pi*, and thus its own *Tci* and resistivity curve. We have used Equations (7)–(9) for the *ρi*(*T*, *p*) functionality of each monodomain *i*. The model also includes an I-bias power source and a voltmeter connected with zero-resistance contacts to opposite edges of the sample (see, e.g., scheme in Figure 1) and the *R*(*T*) of the film will be obtained as the external *V*/*I*. Calculating the circuit requires to numerically invert, for each temperature, the sparse matrix with dimensions 40, 001 × 40, 001 that defines the mesh-current equations. This is a parallelizable computation for which we employed a 31 Tflops supercomputer (LBTS-*ε*psilon, that comprises about 12,000 floating-point units and is described in [56]). It was 100% allocated to run our software during several weeks.

**Figure 1.** Electrical resistivity *ρ* versus temperature *T* obtained for YBa2Cu3O*<sup>δ</sup>* (YBCO) films with a single, uniform value for the nominal doping level *p*, including the case with negligible *Tc* nanostructuration (or maximum-*Tc* doping, *p* = 0.155, in which *Tc* saturates near its maximum value and the *Tc* disorder is negligible) and various cases in which the *p* value corresponds to significant *Tc* nanostructuration (see Section 4 for details). The data points correspond to the finite-element computations, and the continuous lines to the analytical effective-medium (EM) equations.The shaded gray region signals the operational temperature range Δ*T* (in which *ρ* is strongly dependent and linear in *T*). In the upper drawing we illustrate the simulated sample and setup, also including a zoom at smaller length scales, imaging the spatial variation of the local doping level *p*(*x*, *y*) (each *p*(*x*, *y*) monodomain has typical size (30 nm)<sup>2</sup> and the distribution is Gaussian around the average *p*(x), see main text for details).

We have performed our calculations with numerical values representative of the HTS compound YBCO and therefore for the area of a finite-element monodomain *i* we used (30 nm)2, that is expected to correspond to the size of a doping *Tc* inhomogeneity in YBCO. [32,54,55] Therefore, the surface of the simulated HTS film is going to be (6 μm)2, in agreement with the microsensor HTS TES device implementation of Section 2.1.1.

In the case of the nonpatterned HTS considered in this section, the only spatial variation of doping and *Tc* is the unavoidable intrinsic dopant ion structuration and, thus, we assign the local *pi* and *Tci* value to each of our 200 × 200 monodomains *i* as follows: We first build a set of 40,000 values of dopings following the Gaussian distribution given by Equation (10). We then assign each of those *p*-values to each node *i* randomly. Finally, those *pi* are transformed to *Tci* values (and corresponding *ρi*(*T*) functions) following the quantitative results of [38]. A scheme of an example of the resulting spatial distribution is provided in Figure 1 (note the random nanostructuration in the zoomed area).

#### *2.6. Analytical Estimates Using an Effective-Medium Approximation*

To additionally probe the consistency of our computations, we will use, as a useful test, semi-analytical results that we calculate using the so-called effective-medium equations (EM approximation). The EM approximation was first introduced by Bruggeman [57] for general random inhomogeneous materials, and then adapted, e.g., by Maza and coworkers [31] for HTS with Gaussian random *Tc* distributions. As shown in those early works, the EM approximation is a coarse-averaging model that may be considered accurate for temperatures not too close to the *R* → 0 point (at which percolation effects may be expected to invalidate the approximation). In the case of our 2D media, the EM equations may be summarized as the following implicit condition for the conductivity *σ* of each region with random doping inhomogeneities [31,33–35]:

$$\int\_{0}^{\infty} \frac{\sigma\_{p} - \sigma}{\sigma\_{p} + \sigma} \, \omega(p, \overline{p}) \mathrm{d}p = 0. \tag{11}$$

Here, *p*, *p* and *ω*(*p*, *p*) retain the same meaning as in Equation (10), and *σ<sup>p</sup>* is the electrical conductivity corresponding to doping level *p*. The above equation has to be numerically solved to obtain *σ*; however, the computational weight is much lower than the finite-element computation method (seconds versus hours or even days in our parallel computer [56]).

#### **3. Additional Methods for Structured and Patterned Resistive HTS TES**

We now describe the additional methods needed to obtain the *R*(*T*) curve of HTS films in which, additional to the random nanostructuration considered in previous section, also a regular pattern of nominal doping levels is imposed, with the aim to obtain designs that optimize the bolometric functionality. In these films, a regular spatial variation of the nominal doping level *p* is created by the samples' grower by using any of the different methods for micro- and nanostructuration developed in the recent years by experimentalists in HTS films (see, e.g., [21–30]; for instance, for YBCO, this is possible by local deoxygenation using cover masks, ion bombardment, etc.). In particular, all the specific example patterns considered in this work will be expressible as functions *p*(*x*), where *x* is the coordinate in the direction parallel to the external bias current (see scheme in Figure 1) and thus it will be useful to introduce the corresponding function *λ*(*p*), or relative length weight of each nominal *p* value in the film, defined as

$$
\lambda(\overline{p}) = \frac{1}{L} \frac{d\mathbf{x}}{d\overline{p}} \tag{12}
$$

where *L* is the total length of the film in the *x*-direction. Crucial for our studies, one has still to add to these nominal *p*(*x*) variations the unavoidable nanometric-scale randomness of the doping level (considered in the previous sections), i.e.,

$$p(\mathbf{x}, y) = \overline{p}(\mathbf{x}) + p\_{\text{random}}(\mathbf{x}, y), \tag{13}$$

with *p*random(*x*, *y*) consistent with Equation (10) evaluated using the local *p*(*x*).

#### *3.1. Obtainment of the R*(*T*) *Curve of Patterned HTS Using Finite-Element Computations*

To obtain the *R*(*T*) curves of patterned resistive HTS TES, we use the finite-element software TOSERIS already described in Section 2.5. We again use a 200 × 200 simulation mesh and now we assign to each of those finite elements a local doping as follows: First, we associate to each element *i* a nominal doping *pi* corresponding to the pattern to be simulated. Then we randomly calculate the local doping *pi* following the Gaussian distribution given by Equation (10), evaluated with the nominal doping *pi* of each node. Finally, to each node we assign the *Tci* and *ρi*(*T*) corresponding to their local *pi* as per the quantitative results of [38] for the HTS material YBCO (see Sections 2.2 and 2.3).

We also tested that the sets of nodes sharing the same *pi* value follow the Gaussian distribution, and each *R*(*T*) simulation was repeated for several so-generated samples to verify their reproducibility. These checks indicate that our choice of a 200 × 200 node mesh provides enough statistical size. If we attribute to each node the size (30 nm)<sup>2</sup> corresponding to each *Tc*-monodomain in YBCO (see Section 2.4 and [32,54]), the whole 200 × 200 sample corresponds to (<sup>6</sup> <sup>μ</sup>m)2, that is realistic for a microbolometric pixel.

Unless stated otherwise, we will again use in our calculations the numerical values in Sections 2.2 to 2.5 for the common material characteristics, such as, e.g., a film thickness of 100 nm or values for the critical-fluctuation parameters as per Section 2.3.

#### *3.2. Analytical Estimates Using an Extended-EM Approximation*

Besides performing finite-element computations, we will test our results against estimates based on the EM approach. For that purpose, we must suitably extend this approximation to account for the 1D gradient of nominal dopings corresponding to each example pattern to be considered in this work. For that, we consider the film as an association in a series of domains, each one with its own resistance and nominal doping *p*, so that:

$$R(T) = \int\_{x=0}^{x=L} \frac{d\overline{p}}{\sigma(\overline{p}(x), T)S},\tag{14}$$

where *L* is again the total length of the superconductor, *S* is its transversal surface, and *σ*(*p*(*x*), *T*) is the ohmic conductivity obtained using the monodomain-EM approach, i.e., using Equation (11) for each doping *p*(x). Equation (14) can be also written as the following integration over nominal doping, with the help of the *λ*(*p*) function defined in Equation (12):

$$R(T) = L \int\_{\overline{p}\_0}^{\overline{p}\_L} \lambda(\overline{p}) \frac{\mathrm{d}x}{\sigma(\overline{p}(\mathbf{x}), T)S},\tag{15}$$

For a discrete distribution (stepwise function *p*(*x*)), the above equation becomes instead a summation:

$$R(T) = \sum\_{i=1}^{N} \frac{L\_i}{\sigma(\overline{p}\_i, T)S},\tag{16}$$

where *N* is the number of discrete domains, each with its own nominal doping *pi* and length *Li*.

Note that Equations (14) to (16) do not explicitly take into account the transverse currents when associating the different domains *pi* and *Li*. Non-longitudinal transport inside each domain is built-in by using the EM approach for each *σ*(*pi* , *T*). However, this approximation may be expected to fail when it has to describe percolations (because both the sum in a series of domains and the EM approximation do not take them into account). Therefore it could be expected to overestimate the value of *R*(*T*) in the very close proximity to the fully superconducting *<sup>R</sup>*(*T*) <sup>→</sup> <sup>0</sup><sup>+</sup> state.

#### **4. Results for Structured Nonpatterned HTS Materials**

In the reminder of this article, we describe the results of applying our methods to different structured HTS materials. We consider first the case of nonpatterned HTS materials, i.e., those with a uniform nominal doping *p*. As already mentioned, the doping level *p* 0.155 (corresponding to the stoichiometry YBa2Cu3O6.93) is the one that has been used up to now to experimentally produce HTS TES, and, due to the saturation of *Tc* near such doping level, it corresponds to a HTS film without a *Tc*-nanostructure. However, the cases with uniform *p* < 0.155 correspond to films with random *Tc*-nanostructuration.

#### *4.1. Results for the R*(*T*) *Profile and Operational Parameters for Resistive TES Case*

In Figure 1 we show the *R*(*T*) resulting from our finite-element computations for the case of YBCO films with nominal dopings *p* = 0.135, 0.140, 0.147 and 0.155, represented as circles, squares, diamonds, and triangles, respectively. As previously mentioned, to these nominal dopings, a random intrinsic structuration has to be added (following a Gaussian distribution as per Equation (10)) in order to obtain the local *p*(*x*, *y*) values; this is illustrated by the zoom in the pictured drawn in Figure 1.

As it is clearly shown in that figure, the reduction of *p* induces not only a shift of the transitions towards somewhat lower temperatures (as expected from the *Tc*-vs-*p* phase diagram of HTS), but also a widening of the *T*-width of the transition region. The gray areas in Figure 1 are the *T*-range, Δ*T*, in which the *R*(*T*) transition occurs and *R*-vs-*T* is linear, i.e., the operational interval Δ*T* of Equation (3). This Δ*T* increases as *p* decreases, as may be noticed both in Figure 1 and also in Table 1. In this table we summarize the bolometric operational parameters obtained by using such *R*(*T*) results and the methods in Section 2 for three example resistive HTS TES devices designs (described in Sections 2.1.1 and 2.1.2).

It is clear that the improvement that the *Tc*-nanostructure may provide for such bolometric characteristics: Notably, the increase of Δ*T* (about five times higher for the more nanostructured case than in the nonstructured one) is translated in an enhancement of *P*max of the bolometer device. This means the ability to receive a higher amount of radiation without saturating the sensing material. Note that a larger Δ*T* also implies a less demanding cryogenic setup in terms of required stability.

Therefore, our studies suggest that this first simple *Tc* nanostructuration already improves some of the bolometric operational parameters with respect to the conventional, nonstructured HTS materials proposed until now.

**Table 1.** Summary of the main operation parameters for resistive high-*Tc* cuprate superconductors (HTS) transition-edge bolometer devices (TES) devices using, as sensing materials, YBCO superconductors with the various doping nanostructurations and patterns explored in this work. These include the usual nonstructured YBCO (*p* = 0.155) and the novel options considered by us: structured nonpatterned YBCO (*p* ≤ 0.155, Section 4) and structured patterned YBCO (Sections 5 to 7). The Δ*T* and TCR follow from the *ρ*(*T*) curves presented in Figures 1 and 3–5. The *P*max values were calculated for the three example device implementations described in Sections 2.1.1 and 2.1.2 (microbolometers with direct cooling and SrTiO3(STO) or silicon/Yttria-stabilized zirconia (YSZ)/Zirconia (CMOS)-type substrates, plus millimeter-wave sensors using a meander geometry).


#### *4.2. Verification Using the Analytical EM Approximation and Against Existing Measurements*

Figure 1 also shows (as continuous lines) the results obtained by applying the EM approach, i.e., Equation (11), to the same parameter values and doping levels as used in the previous subsection. It can be seen in that figure that the coincidence between the finite-element computation and the EM approximation is excellent, even in the *<sup>R</sup>* <sup>→</sup> <sup>0</sup><sup>+</sup> tails (a log–log zoom of such tails evidences moderate deviations in relative values, negligible in the absolute scale of Figure 1, as coherent with the expectation that the EM approximation is less accurate when percolative current paths appear [31,57]; see also Figure <sup>2</sup> for better evidence of these *<sup>R</sup>* <sup>∼</sup> <sup>0</sup><sup>+</sup> deviations. This comparison gives, then, a first argument supporting the validity, at least concerning the main features, of our calculation methods for the doping structuring effects.

A second argument supporting such validity is the comparison with actual measurements for *R*(*T*). In the case of structured nonpatterned HTS films considered in this section, measurements valid for such a comparison do exist. In particular, we have plotted in Figure 2 the data measured in reference [48] in high-quality YBCO (YBa2Cu3O*δ*) films comprised by a single zone of nominal oxygen stoichiometries *O*6.78 and *O*6.85 (i.e., *p* 0.140 and *p* 0.156 respectively; to get the relations between oxygen ratio and doping we have interpolated the experimental data described in reference [58]). We can see in that figure the good accuracy of the theoretical (EM approach) and computational (finite-element) methods used in this article to reproduce the experimental resistance curves of nonpatterned YBCO films in the studied doping range.

**Figure 2.** Comparison between The electrical resistance vs. temperature curves resulting from our methodology (open symbols for the finite-element computations and solid line for the EM analytical approximation) and the measurements of reference [48] (solid symbols), for YBCO films with uniform nominal oxygen stoichiometry corresponding to: in panel a) to *p* = 0.156, i.e., very near the maximum-*Tc* doping level; in panel b) to *p* = 0.140, i.e., a *Tc*-nanostructured nonpatterned HTS. Note that the normalization of the axes significantly zooms the transition region with respect to Figure 1.

#### **5. Results for Structured HTS Materials Patterned with a Linear** *p*(*x*) **Variation**

Let us now study different instances of doping-patterned HTS, seeking to progressively identify pattern designs optimizing the performance as bolometric sensor materials. We start by considering in

this section, a simple linear variation of *p* along the longitudinal direction (the direction of the overall externally-applied electrical current, see the scheme in Figure 3):

$$\overline{p}(\mathbf{x}) = \overline{p}\_0 + \left(\frac{\overline{p}\_L - \overline{p}\_0}{L}\right) \mathbf{x}\_\prime \tag{17}$$

where *p*<sup>0</sup> and *pL* are the *p*-values at the opposite ends of the film *x* = 0 and *x* = *L*. In terms of the *λ*(*p*) function of Equation (12), this linear-in-*x p*-pattern simply becomes the constant value

$$
\lambda\left(\overline{p}\right) = \frac{1}{\overline{p}\_L - \overline{p}\_0}.\tag{18}
$$

For our computations we chose the rather typical values *p*<sup>0</sup> = 0.135 and *pL* = 0.161.

#### *5.1. Results for the R(T) Profile and Operational Parameters for Resistive TES Use*

The results of our numerical finite-element evaluation for this *p*-pattern are displayed in Figure 3 (see also Table 1 for a comparative summary). As evidenced in the Figure 3, this type of structuring of the HTS film significantly broadens the *R*(*T*) transition (compare, e.g., with Figure 1 that represents, in the same *T*-scale, the results for HTS with comparable, but uniform, *p*-values). However, this structuring does not lead to a linear dependence of *R* vs. *T* in that transition region. This may pose a difficulty in resistive TES applications, that ideally require a *R*(*T*) variation both large (i.e., large TCR) and linear (i.e., constant d*R*/d*T*). The range Δ*T* in which both conditions are met is merely about 1.4 K for this type of *p*-pattern, already suggesting that further structuring optimizations would be desirable (see next section). In Table 1, we summarize the operational parameters obtained for the linear *p*(*x*) resistive bolometer. We can conclude that they are of the same order or worse than the parameters obtained for the typical already existing resistive HTS TES (the case *p* 0.155). This even includes the TCR, that is lower due to the increase of the operational temperature *T*− and then of *R*(*T*−). The maximum energy and power may be somewhat higher due to the small increase of the width of the linear regime.

**Figure 3.** In the upper row, we illustrate a YBCO film patterned following the linear variation of the nominal doping *p*(*x*) studied by us in Section 5 (see Equations (17) and (18)). We also illustrate doping across the film as a 2D color map, taking into account that the local doping level *p*(*x*, *y*) (zoom in the picture) results from accumulating the lineal *p*(*x*) and the random Gaussian disorder at smaller length scales of about (30 nm)<sup>2</sup> (see main text for details). In the lower row, we plot the resistivity vs. temperature *ρ*(*T*) that we obtain for such film (data points for finite-element computations, continuous line for the analytical EM approximation, see Equation (19)). Note that the transition widens considerably with respect to Figure 1, but the operational Δ*T* range (shaded gray region) is small due to the nonlinearity of *ρ* with *T* in most of the transition.

We can conclude that this first patterning does not effectively optimize the operational parameters of the resistive HTS TES, mainly because it broadens the *R*(*T*) transition but does not achieve *R*(*T*) linearity in it.

#### *5.2. Verification Using the Extended-EM Analytical Approximation*

To check the validity of our results for the linear *p*(*x*) variation, we also have used the formulae described in Section 3.2. For that, our first step has been to combine Equation (15) with the *λ*(*p*) formula for this type of pattern (Equation (18)). This leads us to the new equation:

$$R(T) = \frac{L}{S(\overline{p}\_L - \overline{p}\_0)} \int\_{\overline{p}\_0}^{\overline{p}\_L} \frac{d\overline{p}}{\sigma(\overline{p}, T)} \quad \text{(linear } \overline{p}(\mathbf{x}) \text{ pattern)},\tag{19}$$

where *σ*(*p*, *T*) results from Equation (11). The result of this analytical estimate is displayed in Figure 3 as a continuous line. As in the case described for constant nominal doping, this estimate achieves good agreement with the finite-element computation, confirming the basic accuracy of our results.

#### **6. Results for Structured HTS Materials Patterned with a Continuous Exponential-Like Doping Variation**

Seeking to find a *p*(*x*) profile producing a *R*(*T*) transition that improves the bolometric operational characteristics, we have explored numerous *p*(*x*) options beyond the simple linear function discussed above. In the present section we present the results that we obtained with the continuous *p*(*x*) functionality that led us to better bolometric performance (and a step-like, noncontinuous variation will be later discussed, in Section 7). This continuous *p*(*x*) profile is more intuitively described by means of the length weight function *λ*(*p*). In particular, we consider *p*-profiles leading to the following exponential *λ*(*p*) function:

$$\lambda\left(\overline{p}\right) = A \exp\left(\frac{\overline{p}\_0 - \overline{p}}{\delta \overline{p}}\right),\tag{20}$$

where *δp* and A are constants, the latter being easy to obtain by normalization considerations as

$$A = \frac{1}{\delta \overline{\mathcal{p}}} \frac{1}{1 - \exp\left(\frac{\overline{\mathcal{p}}\_0 - \overline{\mathcal{p}}\_L}{\delta \overline{\mathcal{p}}}\right)}.\tag{21}$$

In these equations, *p*<sup>0</sup> and *pL* are, as in the previous sections, the nominal doping at *x* = 0 and *x* = *L* respectively, being *L* the size of the film. Note that, by applying Equation (12), this corresponds to:

$$p(\mathbf{x}) = \overline{p}\_0 - \delta \overline{p} \ln \left\{ 1 - \frac{\mathbf{x}}{L} \left[ 1 - \exp \left( \frac{\overline{p}\_0 - \overline{p}\_L}{\delta \overline{p}} \right) \right] \right\}. \tag{22}$$

For the case of YBCO films considered in this article, and for *p*<sup>0</sup> = 0.135 and *pL* = 0.161 as in the previous section, we found that the *δp* value that best optimizes the bolometric characteristics (most notably Δ*T*) is *δp* = 0.007. We also employed in our evaluations the same common parameter values as described in Sections 2.3 to 2.5.

In the upper row of Figure 4, the corresponding doping profile is pictured, both as a *p*(*x*) representation and as a 2D color density plot. It may be noticed that at the qualitative level the *p*(*x*) function itself is not too dissimilar to an exponential (however a purely exponential dependence of *p* with *x* would produce less optimized bolometric performance).

#### *6.1. Results for the R(T) Profile and Operational Parameters for Resistive TES Use*

The results of our numerical finite-element evaluation for this *p*-pattern are displayed in the second raw of Figure 4 (see also Table 1). As evidenced in that Figure 4, not only the transition is significantly broadened with respect to nonpatterned HTS but also (unlike what happened in the case of a linear *p*(*x*) variation) Δ*T* is highly increased. In particular, the Δ*T* region is increased up to 8.3 K, almost 10 times more than for nonstructured HTS.

As may be seen in Table 1, the improvements also occur in the *P*max parameter, that increase about one order of magnitude with respect to nonstructured HTS. However, note that the TCR value is one order of magnitude worse than in the case of such conventional, nonstructured HTS. This shortcoming and other improvements will be addressed in Section 7 with an evolved *p*-pattern design.

#### *6.2. Verification Using the Extended-EM Analytical Approximation*

We have checked our results also using an analytical estimate, by adapting to this pattern the effective-medium approach described in Section 3.2. For that, we have combined Equation (15) with the Equation (22) defining this *p*-pattern, to obtain the new formula:

$$R(T) = \frac{AL}{S} \int\_{\overline{\overline{p}}\_0}^{\overline{p}\_L} \exp\left(\frac{\overline{p}\_0 - \overline{p}}{\delta \overline{\overline{p}}}\right) \frac{d\overline{p}}{\sigma(\overline{p}, T)} \quad \text{(exp-like pattern)},\tag{23}$$

where *σ*(*p*, *T*) results from the Equation (11). The result of this analytical estimate is displayed in Figure 4 as a continuous line. Again it slightly overestimates *R*(*T*) in the tail of the transition, but basically confirms the finite-element results. As already mentioned, the overestimation is expected to be linked to precursor percolation effects.

**Figure 4.** In the upper row, we illustrate a YBCO film patterned following the exponential-like variation of nominal doping *p*(*x*) studied by us in Section 6 (see Equations (20) to (22)). We also illustrate doping across the film as a 2D color map, taking into account that the local doping level *p*(*x*, *y*) (zoom in the picture) results from accumulating the exponential-like *p*(*x*) and the random Gaussian disorder at smaller length scales of about (30 nm)<sup>2</sup> (see main text for details). In the lower row, we plot the resistivity vs. temperature *ρ*(*T*) that we obtain for such film (data points for finite-element computations, continuous line for the analytical EM approximation, adapted in this work to this *p*(*x*) case, see Equation (23)).

#### **7. Results for Structured HTS Materials Patterned with a Four-Step Exponential-Like Doping Variation**

While the *p*(*x*)-pattern design considered in the previous section produced notable improvements of the bolometric features, at least two concerns may be expressed about it: First, any current structuration experimental technique [21–30] may have difficulties producing such a smooth and exponential-like variation of *p* with *x* (Equations (20) to (22)). Instead, it would be preferable a simpler and, mainly, *discrete*-pattern, i.e., one comprised of a few zones, each with a single *p*. This would ease fabrication, e.g., by means of several stages of deoxygenation of YBCO films using different cover masks in each stage. Secondly, the continuous-pattern of the previous section presents somewhat worsened TCR value with respect to some of the nonpatterned resistive HTS.

To address both issues, we consider now the discrete *p*(*x*) pattern described in the upper row of Figure 5. This pattern defines four zones, each with a single uniform *p* chosen to optimize the linear region of the transition. These values of *p* follow a discretized version of the exponential pattern:

$$L\_i = B \exp\left(\frac{\overline{p}\_0 - \overline{p}\_i}{\delta \overline{p}}\right),\tag{24}$$

where *Li* is the length of the zone of nominal doping *p*<sup>i</sup> , and *B* is a constant so that ∑*<sup>i</sup> Li* = *L*.

We tested the bolometric performance for various doping levels *pi* of the four zones. We obtained the best results with the set *pi* = {0.136, 0.141, 0.145, 0.160}, that we describe next.

#### *7.1. Results for the R(T) Profile and Operational Parameters for Resistive TES Use*

The results for our numerical finite-element evaluation for this *p*-pattern are displayed in Figure 5 (see also Table 1). As evidenced there, the *R*(*T*) transition becomes significantly broad and linear, with such linear region conveniently starting at *T*− = 76.6 K (so that the HTS TES could be operated with the simplest liquid-nitrogen bath, at 77 K). The corresponding Δ*T* is now almost 13 K, the largest obtained in this paper. Also the TCR value > 5 K−<sup>1</sup> is the largest obtained in this work, being almost double than for conventional nonstructured (i.e., maximum *Tc*) YBCO. The *P*max values (see Table 1) also reflect these improvements, being again the best among the HTS options considered in this work and more than one order of magnitude larger than for the nonstructured HTS.

To sum up, these finite-element computations reveal that this relatively simple-to-fabricate *p*-pattern produces order-of-magnitude improvements over nonstructured HTS materials in Δ*T* and *P*max, and also a 66% improvement in TCR.

#### *7.2. Verification Using the Extended-EM Analytical Approximation*

We have checked our results for this four-step *p*(x) pattern also using an analytical estimate, by adapting to this pattern the approach described in Section 3.2. In this case we used their discretized version , given by Equations (16). By combining it with the *p*(*x*)-pattern given by Equation (24) we now obtain:

$$R = \frac{B}{S} \sum\_{i=1}^{N} \frac{1}{\sigma \left(\overline{p}\_{i\prime} \, T\right)} \exp\left(\frac{\overline{p}\_0 - \overline{p}\_i}{\delta \overline{p}}\right),\tag{25}$$

where the *pi* are the nominal dopings of each of the *i* zones, and *σ*(*pi* , *T*) results from Equation (11). This analytical estimate is displayed in Figure 5 as a continuous line. It fully confirms the main features obtained by the finite-element computations. Similarly to the case of the other *p*-patterns considered in our work, the estimate is expected to be less reliable in the lower part of the transition.

**Figure 5.** In the upper row, we illustrate a YBCO film patterned following the discretized four-step exponential-like variation of the nominal doping *p*(*x*) studied by us in Section 7. We also illustrate doping across the film as a 2D color map, taking into account that the local doping level *p*(*x*, *y*) (zoom in the picture) results from accumulating the four-step exponential *p*(*x*) and the random Gaussian disorder at smaller length scales of about (30 nm)<sup>2</sup> (see main text for details). In the lower row, we plot the resistivity vs. temperature *ρ*(*T*) that we obtain for such film (data points for finite-element computations, continuous line for analytical EM approximation, see Equation (25)). This is the most optimized design found in this work for the use in resistive HTS TES devices (see also Table 1).

#### **8. Conclusions**

To sum up, we considered the advantages of structuring and patterning of the doping level (and hence of the critical temperature) in high-temperature superconductors with respect to their operational characteristics as resistive bolometric sensors (resistive HTS TES) of electromagnetic radiation. In particular we studied some chosen examples of spatial variations of the carrier doping into the CuO2 superconducting layers due to oxygen off-stoichiometry. Our main results are (see also Table 1 for a quantitative account):

*(i)* Non-patterned structured HTS materials (i.e., those with a nominal doping level uniform in space but that does not maximize the critical temperature, thus having random *Tc*-nanostructuring) may already provide some benefit for bolometric use with respect to the nonstructured HTS materials used up to now for those devices. In particular, they present a widened transition leading to a larger operational temperature interval Δ*T* and also larger *P*max (corresponding to the larger maximum detectable radiation power before sensor saturation). However, these improvements come at the expense of a certain reduction of the sensibility of the sensor as measured by the TCR value.

*(ii)* The bolometric performance may be significantly more optimized with the use of HTS materials including an additional regular dependence on the position of the nominal doping level, *p*(*x*) (doping-level patterning). In that case, ad-hoc pattern designs may be found by progressively seeking widened and linear *R*(*T*) transitions. Our more optimized design is shown in Figure 5 and consists of just four zones of different sizes and doping levels (related by the exponential-like Equation (24) evaluated at *pi* = 0.136, 0.141, 0.145, and 0.160). With this design the operational temperature is conveniently located at *T*− = 76.6 K, the operational temperature interval Δ*T* is almost 13 K (more than one order of magnitude larger than for the conventional nonstructured YBa2Cu3O6.93), the TCR value is > 5 K−<sup>1</sup> (almost double than for the nonstructured case), and the *P*max values are also optimized about one order of magnitude.

**Author Contributions:** J.C.V. and M.V.R. derived the EM equations. J.C.V., A.S.V., M.M.B., and M.V.R. wrote the manuscript. J.C.V., A.S.V., C.M.-O., and M.V.R. contributed to the design of the studied patterns. All authors contributed to the computations and analyses of the results.

**Acknowledgments:** This work was supported by projects FIS2016-79109-P (AEI/FEDER, UE) and AYA2016-78773-C2-2-P(AEI/FEDER,UE), by the Xunta de Galicia under grants ED431D 2017/06 and ED431C 2018/11, the Consellería de Educación Program for Development of a Strategic Grouping in Materials AeMAT under Grant No. ED431 2018/08, Xunta de Galicia, and by the CA16218 nanocohybri COST Action. JCV thanks the Spanish Ministry of Education for grant FPU14/00838.

**Conflicts of Interest:** "The authors declare no conflict of interest."

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Transport and Point Contact Measurements on Pr1**−**xCexPt4Ge12 Superconducting Polycrystals**

#### **Paola Romano 1,2, Francesco Avitabile 1, Angela Nigro 3,4, Gaia Grimaldi 2, Antonio Leo 3,4, Lei Shu 5,6, Jian Zhang 5, Antonio Di Bartolomeo 2,3 and Filippo Giubileo 2,\***


Received: 30 July 2020; Accepted: 7 September 2020; Published: 10 September 2020

**Abstract:** We performed a detailed investigation of the superconducting properties of polycrystalline Pr1−xCexPt4Ge12 pellets. We report the effect of Ce substitution, for x = 0.07, on magnetic field phase diagram H-T. We demonstrate that the upper critical field is well described by the Ginzburg–Landau model and that the irreversibility field line has a scaling behaviour similar to cuprates. We also show that for magnetic fields lower than 0.4 T, the activation energy follows a power law of the type *H*−1/2, suggesting a collective pinning regime with a quasi-2D character for the Ce-doped compound with x = 0.07. Furthermore, by means of a point contact Andreev reflection spectroscopy setup, we formed metal/superconductor nano-junctions as small as tens of nanometers on the PrPt4Ge12 parent compound (x = 0). Experimental results showed a wide variety of conductance features appearing in the dI/dV vs. V spectra, all explained in terms of a modified Blonder–Tinkham–Klapwijk model considering a superconducting order parameter with nodal directions as well as sign change in the momentum space for the sample with x = 0. The numerical simulations of the conductance spectra also demonstrate that s-wave pairing and anisotropic s-waves are unsuitable for reproducing experimental data obtained at low temperature on the un-doped compound. Interestingly, we show that the polycrystalline nature of the superconducting PrPt4Ge12 sample can favour the formation of an inter-grain Josephson junction in series with the point contact junction in this kind of experiments.

**Keywords:** superconductivity; transport properties; energy gap; superconducting order parameter; proximity effect; nano-junction; Andreev reflection

#### **1. Introduction**

Filled skutterudite materials have attracted a great deal of attention for a large number of properties such as metal–insulator transitions, spin fluctuations, and heavy fermion behaviour [1–4]. Several compounds in the family of filled skutterudites also show the phenomenon of superconductivity [5–8]. They have the chemical formula MT4X12, where M is an electropositive metal (Sr, Ba, La, Pr, Th), T is a transition metal (Fe, Os, or Ru), and X usually represents a pnictogen (Sb, As, or P). The first Pr-based superconductor to be discovered was the heavy-fermion PrOs4Sb12, with a critical temperature Tc = 1.85 K, showing intriguing properties such as a giant electronic specific heat

coefficient [1]. Moreover, experiments of thermal transport [9] on single crystals evidenced the possible existence of a superconducting phase at high magnetic fields in which the energy gap has at least four point nodes, and a second phase at low magnetic fields in which the energy gap is characterized by only two point nodes. Recently, a new Pt-based family of skutterudite, with chemical formula MPt4Ge12, was synthetized, showing superconducting properties at relatively high temperatures. In particular, the compound with praseodymium (Pr) as metal shows a transition temperature Tc = 7.9K, while for the compound with Lanthanum (M = La), Tc = 8.3 K has been reported [10], as confirmed by electrical resistivity, magnetic susceptibility and specific heat measurements. Nuclear magnetic resonance experiments have given indications for conventional superconductivity in LaPt4Ge12 [11].

Lower superconducting critical temperature was previously reported for other MPt4Ge12 compounds such as SrPt4Ge12 (Tc = 5.10 K) and BaPt4Ge12 (Tc = 5.35 K) [12]. The higher critical temperatures for Pr and La compounds with respect to Sr and Ba compounds have been explained as the existence of a larger density of states at the Fermi level, as resulting from 73Ge nuclear quadrupole resonance experiments at zero field [13]. It has also been suggested that PrPt4Ge12 and LaPt4Ge12 can be characterized by two superconducting gaps. Indeed, according to heat capacity measurements as a function of temperature and magnetic field, the superconducting state cannot be explained by considering a single isotropic or anisotropic energy gap [14]. The presence of two distinct linear regions in the magnetic field dependence of the Sommerfeld coefficient of electronic heat capacity was interpreted as a possible indication for two-gap superconductivity in these compounds. The critical current density and pinning force of superconducting PrPt4Ge12 was measured in magnetization experiments [15], revealing that dependence of both quantities with respect magnetic field can be explained using a double exponential model already developed to explain the properties of the two-band superconductor MgB2 [16–19].

μSR experiments on PrPt4Ge12 showed a time-reversal symmetry breaking below Tc [20–23], in contrast to the results on LaPt4Ge12, for which the time-reversal symmetry breaking is absent and a conventional superconductivity with a fully gapped density of states is supposed [21]. The superconducting order parameter of LaPt4Ge12 has also been studied using specific heat and thermal conductivity measurements [24], showing that the sharp transition in the specific heat and its zero-field temperature dependence are well described in a conventional BCS (Bardeen–Cooper–Schrieffer) scenario characterized by a single energy gap and s-wave symmetry [25]. On the other hand, de Haas–van Alphen measurements have been reported with state-of-the-art band-structure calculations showing that LaPt4Ge12 and PrPt4Ge12 have almost identical electronic structures, Fermi surfaces and effective masses [26]. So far, few investigations have been reported that probe the superconducting energy gap in the PrPt4Ge12 compound, with results supporting both nodal and nodeless energy gaps. Recently, electrical resistivity, magnetic susceptibility, specific heat, and thermoelectric power experiments were performed on Pr(1-x)CexPt4Ge12 to investigate the influence of the magnetic state of the Ce ions on the superconducting properties of the compound [27]. Interestingly, the results indicate a crossover from a nodal to a nodeless superconducting energy gap and that PrPt4Ge12 could be a two-band superconductor in which the electron scattering due to Ce substitution can suppress the superconductivity within one of the bands.

In this paper, we perform a detailed study of the superconducting transport properties of the Pr1−xCexPt4Ge12 compound for x = 0.07. We measured the resistive transitions for the sample with x = 0.07 in external applied magnetic fields. We deduced the H-T phase diagram of the Ce-doped sample Pr0.93Ce0.07Pt4Ge12, that is the upper critical field Hc2, as well as the irreversibility line. We also analysed the resistive transition data in the framework of the thermally assisted motion of vortices. We also performed direct measurements of the superconducting energy gap in the parent compound PrPt4Ge12 (x = 0) by means of point contact spectroscopy experiments, which made it possible to realize metal/superconductor nano-junctions with dimensions of few nanometres. We measured the conductance spectra of the point contact junction at low temperature (4.2 K) and we demonstrated that the conductance feature can be reproduced in a theoretical model that takes into account the

symmetry of the superconducting order parameter with nodal directions and change of sign in the momentum space. We also estimated the superconducting energy gap for the sample PrPt4Ge12 in the range 0.55–0.95 meV.

#### **2. Materials and Methods**

Polycrystalline Pr1−xCexPt4Ge12 pellets (with x = 0, 0.07, 0.1) were synthesized in argon atmosphere by arc melting, using Pr ingots, Ce rods, Pt sponge, and Ge pieces as preparation materials, weighed in stochiometric ratios. Arc melting and turning over were repeated five times to obtain high chemical homogeneity. The samples were then annealed in a sealed quartz tube in 200 Torr argon atmosphere at 800 ◦C for 14 days. X-Ray diffractometry (reported elsewhere [27]) confirmed the sample quality evidencing the expected cubic skutterudite crystal structure.

Measurements of sample resistance as a function of the temperature, *R*(*T*), were carried out by standard four-probe offset-compensated technique inside a Cryogenic Ltd. cryogen-free magnet system, equipped with a variable temperature insert (vertically inserted inside a superconducting solenoid for fields up to 9T) in which the temperature can vary in the range 1.6 K–325 K. A Keithley 2430 DC current source was employed to bias the sample, while measuring the voltage by means of a Keithley 2182 nanovoltmeter.

Point Contact Spectroscopy experiments were performed by means of a home-built mechanical inset equipped with a screw-driven chariot to move a metallic tip towards the surface of the superconducting sample, in order to form a metal/superconductor nanometric constriction, the so-called point contact junction. The junction is then cooled down by immersing the inset in a liquid helium cryostat to perform current-voltage (I–V) measurements at low temperature by standard four-probe technique. Differential conductance (dI/dV vs. V) spectra are then obtained performing numerical derivative of the I–V curves.

#### **3. Results and Discussion**

#### *3.1. Transport Properties*

The dynamic behaviour of Abrikosov vortices in type II superconductors determines the transport properties of superconducting samples. In particular, vortices are set in motion when the Lorentz force, due to the bias current, exceeds the strength of vortex pinning forces, causing a nonzero Ohmic resistance. The presence of pinning, caused by any spatial inhomogeneity of the material like impurities, point defects, grain boundaries, etc., allows superconducting material to sustain current without flux motion and dissipation, giving rise to a nonzero critical current density. Different phases can be recognized in the magnetic field phase diagram *H-T* depending on the relative strengths of the pinning potential, the Lorentz driving energy, and the elastic energy of the vortex lattice, thermal energy and dimensionality. The interplay between these interactions results in several phase separation lines in the *H-T* phase diagram [28].

In the following, we characterize the *H-T* phase diagram of the sample Pr0.93Ce0.07Pt4Ge12, the upper critical field Hc2, and the irreversibility line, above which the critical current density becomes zero. Then, the resistive transitions in external applied magnetic fields are analysed in the framework of the thermally assisted motion of vortices.

#### 3.1.1. Magnetic Field Temperature Phase Diagram

In Figure 1a we report the resistance vs. temperature curves, *R*(*T*), for Pr1−xCexPt4Ge12 samples having x = 0, 0.07 and 0.10. The data were normalized to the normal state resistance, *R*N, evaluated just before the onset of the superconducting transition. We measured *R*<sup>N</sup> = 4.3 mΩ for the undoped sample, *R*<sup>N</sup> = 34 mΩ for sample with x = 0.07 and *R*<sup>N</sup> = 63 mΩ for sample with x = 0.10. The superconducting critical temperature *T*c was estimated at 50% of the onset transition resistance, obtaining *T*c = 7.9 K for the undoped sample, *T*<sup>c</sup> = 4.7 K for samples with x = 0.07, and *T*<sup>c</sup> = 3.6 K for samples with *x* = 0.1. The evolution of the critical temperature *T*c as a function of the Ce doping *x* is summarized in Figure 1b. The effect on Tc of the partial substitution of Pr by Nd has been also reported in samples with Nd content xNd up to 0.1 [29]. The critical temperature is weakly dependent by Nd content, being reduced by only 10% at xNd = 0.1. On the other hand, the effect of the Ce substitution is much more important; the Tc is reduced by 59% at x = 0.07 and by 45% at x = 0.1. The effect of externally applied magnetic field up to 1 T on the R(T) curve for sample with x = 0.07 is shown in Figure 1c.

**Figure 1.** (**a**) The resistance as a function of the temperature normalized to the normal state resistance value for the three different Ce doping. (**b**) Evaluated critical temperature values as a function of the doping. The inset shows a schematic diagram of the Pr1−xCexPt4Ge12 crystal structure, showing the Pr,Ce atoms residing in icosahedral cages formed by tilted PtGe6 octahedral. (**c**) Magnetic field dependence of the resistance versus temperature curve measured for the sample with x = 0.07 doping.

In Figure 2a, the resulting magnetic field–temperature (*H-T*) phase diagram is shown, together with the upper critical field line *H*c2(T), which separates the normal and the superconducting state, and the irreversibility line *Hirr*(*T*). In the *H-T* phase diagram, the *Hirr*(*T*) line separates the Abrikosov vortex pinned regime and the vortex liquid regime, which is at a given temperature; the critical current density goes to zero at the irreversibility field. To evaluate the upper critical field, μ0*H*c2, we employed to the 90% of normal state resistance *R*<sup>N</sup> criterion, while the irreversibility line, μ0*H*irr, was obtained using the 10% of *R*<sup>N</sup> criterion. To analyse the temperature dependence of the upper critical field, the μ0*H*c2(*T*) data extracted by the *H-T* phase diagram are shown in Figure 2b.

**Figure 2.** (**a**). The *H-T* phase diagram as obtained by R(T) measurements at different applied magnetic fields for the sample with the doping x = 0.07. The irreversibility line (0.1 RN curve) and the upper critical field behaviour (0.9 RN curve) are indicated. (**b**) The upper critical field *Hc*2(*T*) data are shown as black spheres. The red solid line is a fit of the data by the Ginzburg–Landau Equation (Equation (1) in the text). The blue dashed line is obtained by the single band WHH model (Equation (2) in the text). (**c**) Scaling plots of the normalized upper critical field Hc2(T)/Hc2(0) plotted as function of the reduced temperature T/Tc for the parent compound PrPt4Ge12 [14,29,30], Nd-doped samples Pr1−xNdxPt4Ge12 [29] and for our Pr1−xCexPt4Ge12 sample with x = 0.07. Hc2(T) curves are obtained by resistivity, R, specific heat, C, and magnetization measurements, M. The inset shows the Hc2 data as a function of the temperature for the same samples in the main panel. (**d**) The irreversibility field *Hirr*(*T*) data are shown as open circles. The red solid line is a fit of the data by Equation (5) in the test.

The red solid line in Figure 2b represents the best fit of the data by the Ginzburg–Landau (G-L) formula:

$$H\_{c2}(T) = H\_{c2}(0)\frac{\left(1 - t^2\right)}{\left(1 + t^2\right)}\tag{1}$$

where *Hc*2(0) is the upper critical field at zero temperature and *t* = *T*/*TC* is the reduced temperature. The data are very well described by the G-L Equation with μ0*Hc*2(0) = 1.23 T, as already reported for (Pr,La)Pt4Ge12 and Pr1−xNdxPt4Ge12 compounds [14,24,29,31–34]. For comparison, in Figure 2b, we also show the temperature dependence of the upper critical field derived within the Werthamer–Helfand–Hohenberg (WHH) model (blue dashed line), which includes orbital and Zeeman pair breaking [35,36]. In particular, for a single band superconductor in a dirty limit, the model yields:

$$H\_{c2} = \frac{2\Phi\_0 \, k\_{\rm B} \, T\_c}{hD\_0} \, ht \tag{2}$$

where Φ<sup>0</sup> is the magnetic flux quantum, *kB* the Boltzmann constant, *Tc* the critical temperature in zero applied magnetic field, *h* is the reduced Planck constant, *D*<sup>0</sup> = *h*/2*me* with *me* the electron mass, and *h* is a parameter that runs from 0 to ∞ as T varies from *Tc* to 0, while the reduced temperature *t* = *T*/*Tc* is given by:

$$\ln(t) = \Psi\left(\frac{1}{2}\right) - \mathrm{Re}\Psi\left(\frac{1}{2} + h(i+d)\right) \tag{3}$$

where Ψ is the digamma function, *d* = *D*/*D*0, and *D* the diffusivity of the band related to the slope of the μ0*Hc*2(*T*) experimental data near *Tc*, given by *D* = 4*kB*/π*e*μ<sup>0</sup> *dHc*2(*T*)/*dT* . For our sample, <sup>μ</sup>0*dHc*2(*T*)/*dT* = <sup>−</sup>0.23 TK<sup>−</sup>1. Within this model, the zero-temperature upper critical field *Hc*2(*0*) can be obtained by the well-known WHH formula:

$$H\_{c2}(0) = 0.693 \times T\_c \left| \frac{dH\_{c2}(T)}{dT} \right| \tag{4}$$

We note that the WHH model underestimates the *Hc*<sup>2</sup> field at low temperatures, with the WHH formula giving μ0*Hc*2(0) = 0.99 T, corresponding to the value experimentally measured at T = 2.1 K (μ0*Hc*<sup>2</sup> = 1 T).

The Ce substitution for Pr in PrPt4Ge12 does not seem to modify the behaviour of the temperature dependence of the upper critical field, showing a characteristic positive curvature near *Tc*, as for the parent compounds (Pr,La)Pt4Ge12 and the doped Pr1−xNdxPt4Ge12 material. To analyse the effects on the upper critical field Hc2(T) curve of the partial substitution of Pr in PrPtGe in Figure 2c, we compare our findings for the sample with Ce doping x = 0.07 with the results reported in the literature for the un-doped parent compound PrPt4Ge12 [14,29,30] and for Nd-doped samples Pr1−xNdxPt4Ge12 [29]. A scaling behaviour described by the G-L formula, Equation (1), is obtained when the normalized upper critical field Hc2(T)/Hc2(0) is plotted as a function of the reduced temperature t = T/Tc. We also point out that the figure includes Hc2(T) curves obtained by resistivity, specific heat, and magnetization measurements. Furthermore, the scaling behaviour is the same for the two different doping, Ce and Nd, despite the different strength of critical temperature lowering induced by the two type of doping.

Different mechanisms have been proposed to explain this behaviour, such as strong coupling effects, multi-band electronic structure, and disorder [37,38].

Moreover, the coherence length at zero temperature was evaluated by μ0*Hc*2(0) = Φ0/4πξ(0) 2 , with Φ<sup>0</sup> being the magnetic flux quantum. A value of ξ(0) ≈ 12 nm was obtained, which is close to the value reported for the un-doped compound, ξ*x*=0(0) ≈ 14 nm [39].

In Figure 2d, the temperature dependence of the irreversibility field is shown. The scaling relation,

$$H\_{irr}(T) = H\_{irr}(0)\left(1 - \left(\frac{T}{T\_{CO}}\right)^2\right)^n \tag{5}$$

was adapted to the experimental data with *n* = 1.5 (red solid line in Figure 2c) and μ0*Hirr*(0) = 0.94*T*. *TC*<sup>0</sup> is the zero field transition temperature, and the exponent *n* is determined by the flux pinning mechanism [40]. At temperatures close to*TC*0, the scaling Equation was reduced to *Hirr* ∼ (1 − *T*/*TC*0) 1.5 as observed by Yeshumn et al. for single crystals of YBa2Cu3O7 high-temperature superconductor and interpreted within the thermally activated flux-creep theory [41]. The exponent *n* = 1.5 is also consistent with the values found for the iron-based 122 and 1111-families and the TlSr2Ca2Cu3Oy compound [42–44].

In the field range between *Hirr*(*T*) and *Hc*2(*T*), the thermal fluctuations become important, and the superconducting state loses its zero-resistance behaviour. In type II superconductors, high-field and -current applications are limited by the irreversibility line *Hirr*(*T*); for example, values of *Hirr* at low temperatures from 50% to ∼80% of the *Hc*<sup>2</sup> have been observed in MgB2 and up to 85% in (Y0.77,Gd0.23)Ba2Cu3Oy films [45,46]. In our sample, the irreversibility field was 70% of the upper critical field at low temperatures, and dropped to less than 20% of *Hc*<sup>2</sup> close to the critical temperature *TC*0.

#### 3.1.2. Temperature Dependence of the Vortex Activation Energy

To further analyse the pinning properties of Pr1−xCexPt4Ge12, we analysed the field dependence of the pinning activation energy *U*. This physical parameter was evaluated by a linear fit of the transition region in the *R*(*T*) curves represented in an Arrhenius plot, which are shown in Figure 3a. Indeed, the Arrhenius plots show that the resistivity is thermally activated over about two orders of magnitude at low fields; in this regime, the resistance's dependence on temperature and field can be written in the form: *R*(*T*, *H*) = *R*0*e*−*U*(*T*,*H*)/*kBT* [47].

**Figure 3.** (**a**) The Arrhenius plot of the R(T) curves for the sample with the doping x = 0.07. (**b**) The pinning activation energy as a function of the applied magnetic field for the same sample. The dotted lines are obtained by a linear fit on the data in the log–log plot.

In Figure 3b, the resulting *U*(*H*) curve is shown. In particular, small activation energies, up to 600 K, are reported, close to the values <sup>∼</sup> <sup>10</sup><sup>3</sup> *<sup>K</sup>* observed in Bi2212 superconductor [48]. According to the literature [34,36], the dependence of the pinning activation energy on the applied magnetic field follows a power law of the type *H*−α, with the exponent α, which assumes different values depending on the vortex pinning regime. In our case it is evident a crossover between two power law behaviours with different exponents. This crossover is at about 0.4 T, and the exponents are α ≈ 0.5 for μ0*H* < 0.4 T and α ≈ 5 for μ0*H* > 0.4 T. Both values can be associated with a collective pinning regime, with an exponent value between 0.5 and 1, which is usually found in cuprate superconductors and could be related to the quasi-2D character of these materials [41,49,50]. The existence of a crossover suggests the presence of two different pinning centres with different dimensions within a collective pinning regime, as was observed, for example, in undoped and Nd-doped PrPt4Ge12 samples [29], as well as YBa2Cu3O7 compounds [50] and in Nd2−xCexCuO4−<sup>δ</sup> thin films [51].

#### *3.2. Point Contact Spectroscopy*

Point contact Andreev Reflection spectroscopy (PCAR) is a very powerful technique, widely applied to investigate the fundamental properties of superconductors, such as the superconducting energy gap amplitude, the density of states (DOS) at the Fermi level, and the symmetry of the superconducting order parameter (OP) [52]. PCAR experiments have been reported to study conventional BCS superconductors [53], high Tc cuprates (both hole doped and electron doped) [54–59], multiband superconductors [60–64], ruthenocuprates [65], iron-pnicniteds [66], heavy fermion superconductors [67], non-centrosymmetric superconductors [68], and topological superconductors [69]. This technique has also been successfully applied for precise measurements of the thickness and of the polarization in thin ferromagnetic/superconductor multilayers [70–75]. The PCAR technique consists of realizing a nano-contact between a tip-shaped normal-metal (N) electrode and a superconductor (S), thus forming a N/S nano-junction. By tuning the transparency of the N/S interface (i.e., by changing the tip pressure on the sample surface) one can realize different tunnelling regimes, going from the Andreev

reflection [76,77] regime (in the case of a low potential barrier at the interface, corresponding to high interface transparency) to quasiparticle tunnelling regime for low interface transparency (high potential barrier). In a typical PCAR experiment, an intermediate regime can be achieved, with both Andreev reflection processes and quasiparticle tunnelling contributing to current transport through the N/S interface. If an electron travels from the normal side of the junction, with an energy lower than the superconducting energy gap, towards the N/S interface, it can enter into the superconducting side only as a Cooper pair, i.e., forming a pair with another electron, while originating a reflected hole in N with the opposite momentum. Consequently, a single Andreev reflection event causes a charge transfer to the S side of 2*e*, with *e* the electron charge. From a theoretical point of view, the transport through a point contact junction between a normal metal and a conventional BCS superconductor (with isotropic s-wave symmetry of the superconducting OP) was described in the BTK theory [78], in which the interface barrier is modelled by a dimensionless parameter Z. The case Z = 0 corresponds to an N/S junction, with a completely transparent barrier, in which the Andreev process is the dominant mechanism responsible for the transport current. On the other hand, Z > 1 represents a junction with a low transparent barrier, corresponding to a dominant tunnelling current flowing through the junction.

The BTK theory was subsequently extended to the case of superconductors with amplitude variation of the OP in the k-space (as for an anisotropic s-wave) and to the case of unconventional superconductors, in which the sign of the OP may also change (as for d-wave symmetry) [79]. It has been demonstrated that if an incident quasiparticle at the N/S interface experiences a different sign of the OP, Andreev bound states are formed at the Fermi energy [80]. From an experimental point of view, the formation of Andreev bound states are seen in the differential conductance spectra as a peak at zero bias [81–85]. In the case of d-wave symmetry, conductance depends on both the incident angle ϕ of the quasiparticle at the N/S interface and on the orientation angle α between the *a-*axis of the superconducting order parameter and the crystallographic axis. The conductance expression can be reduced to the model for anisotropic s-wave by simply assuming costant phase and assuming that only the OP amplitude changes in the *k*-space. We mention here that the BTK model and its extended version applied here do not consider the possible effects due to the case of energy-dependent DOS (the effects being mostly expected on the thermoelectric effects in NS junctions [86].

#### PCAR Experiment

The PCAR experiment was performed on the Pr1−xCexPt4Ge12 sample having x = 0 (Tc = 7.9 K). We used a gold tip as a normal metal electrode, which was gently pushed onto the sample surface to realize the N/S nanoconstriction. The setup was then immersed in a helium liquid bath for low-temperature (T = 4.2 K) characterization. We measured the current–voltage characteristics I–V using a standard four-probe configuration, using a dc current supply to bias the junction and measuring the voltage by a nano-voltmeter. The conductance curves, d*I*/d*V–V,* were obtained by numerical derivation of the I–V curves. The PCAR setup also makes it possible to vary the tip pressure on the sample, obtaining a tuning of the barrier transparency and, consequently, different junction resistances. In our experiment, we obtained junction resistances RN in the range 0.1 Ω−50 Ω. We noticed here that we did not have a direct control of the geometrical dimensions of the N/S junction formed in the PCAR experiment. However, we estimated the junction size through the Sharvin formula *RN* = 4ρ/ - 3π*d*<sup>2</sup> , in which the normal resistance of the junction is related to the resistivity ρ = 3.5 μΩcm [87] and the mean free path = 103 nm [87] in the superconducting material, as well as to contact dimension *d*.

In Figure 4, we report normalized d*I*/d*V–V* curves, measured at low temperature (T = 4.2 K) for different nano-junctions, realized by varying the tip pressure and/or position on the sample surface. The conductance curves in Figure 4a,b are characterized by a zero-bias conductance peak (ZBCP) and have normal resistance RN of 20 Ω and 0.2 Ω, respectively. Based on the Sharvin formula, we found the junction size to be *d* = 7 nm and *d* = 65 nm, respectively. In both cases, this confirms that the point contact is in the ballistic regime [52], in which the size of the junction is smaller than the mean free path in the superconductor (*d* << ). This corresponds to the physical conditions for which an electron

can accelerate freely through the point contact, with no heat generated in the contact region, allowing energy-resolved spectroscopy.

**Figure 4.** Normalized conductance spectra, dI/dV–V, measured at low temperature (T = 4.2 K) in different sample locations on the Pr1−xCexPt4Ge12 (with x = 0) sample. Experimental data (empty symbols) in (**a**,**b**) are compared to numerically calculated curves for the three different symmetries. The I–V curve is shown in the inset. Experimental data (empty symbols) in (**c**,**d**) are compared to numerically calculated curves for d-wave symmetry only. Inset in (**d**) represent the schematic of model in which an inter-grain Josephson junction is formed in series with the point contact junction.

The ZBCPs reported in Figure 4a,b have similar height (∼1.2) and energy width (∼ ±3 meV), while presenting a quite different shape. Experimental data are compared to theoretical fittings obtained for the three mentioned symmetries of the OP. We notice that s-symmetries are not able to completely reproduce the observed features. This discrepancy is more evident in Figure 4a, where the blue solid line represents the simulation obtained assuming a d-wave symmetry, with Δ = 0.95 meV, Z = 0.87 and α = 0.32 as fitting parameters. The Z value gives an indication of an intermediate regime, in which both Andreev reflection and quasiparticle tunnelling contribute to the conduction mechanism. The α value is an indication that the current direction is in between the nodal direction (α = π/4) and α = 0 (corresponding to maximum energy gap amplitude). In Figure 4b, very similar fitting curves are obtained for the different symmetries, although in this case, too, the d-wave symmetry seems to better reproduce the experimental behaviour, with fitting parameters Δ = 0.85 meV, Z = 1.3 and α = 0.38.

On the other side, the spectra reported in Figure 4c,d appear very different. Indeed, in Figure 5c the ZBCP is very narrow and its amplitude is above 10, a value that cannot be obtained in s- or s-anisotropic fitting models, the maximum height being limited to 2. Moreover, at the side of the ZBCP, conductance minima are present, at voltages below ±1 mV. Such a feature is usually expected only when the superconducting OP is characterized by a sign change, as in the d-wave symmetry. Accordingly, we succeeded to simulate the experimental data by using the extended BTK model for a d-wave superconductor, assuming Δ = 0.85 meV, Z = 2.5 and α = 0.39. It was not possible to obtain

similar features by applying the s-wave or s-anisotropic symmetry of the OP. We notice here that all experimental data have been simulated without introducing the so-called Γ-Dynes smearing factor [88] typically used to take into account possible pair-breaking effects of various origins (impurities, inelastic scattering, magnetic field, etc.).

**Figure 5.** (**a**) Conductance spectra measured in a different location on the same superconducting sample: the lower (green) spectrum was measured soon after the tip approach on the surface. The upper spectrum was measured after increasing the tip pressure on the surface. The upper (black) spectrum was vertically shifted (+0.2) for clarity. Solid lines represent the numerical fits. (**b**) Evolution of conductance spectra (solid lines) calculated numerically for Δ = 0.55 meV and α = 0.46, and for 0 < Z < 1. The scattered (green) points refer to experimental data of Figure 5a. (**c**) Evolution of conductance spectra (solid lines) calculated numerically for Δ = 0.55 meV, Z = 0.39, α = 0.29, and 0 < *RJJ* < 0.42 Ω. The scattered (black) points refer to experimental data of Figure 5a.

Another completely different shape is observed in the spectrum of Figure 4d, where a wide ZBCP is followed by several features that cannot be reproduced by simply applying the model discussed above. However, we need to take into account that the superconducting sample is a pellet formed by pressed powders. Consequently, the tip pressure on the surface can cause the formation of an inter-grain Josephson junction in series with the point contact junction as depicted in the inset of Figure 4d. The effect of such inter-grain effects in point contact measurements has already been observed in experiments on MgB2 [61], MgCNi3 [89], and Pr1−xLaCexCuO4-y [57]. In this extended model, we need to take into account that the metallic tip (N electrode) forms a point contact junction on a superconducting grain, and this in turn forms a Josephson junction with another superconducting grain. Consequently, the total voltage drop *V* (experimentally measured) is given by the sum of the point contact *VPC* and the Josephson junction *VJJ* contributions *V* = *VPC* + *VJJ*. If the flowing current *I* is lower than *IJJ* there is no voltage drop at the inter-grain junction (*VJJ* = 0). Otherwise, *VJJ* can be calculated according to Lee formula [90] as *VJJ* <sup>=</sup> *RJJIJJ* - *I*/*IJJ*<sup>2</sup> − 1, where *RJJ* and *IJJ* are the resistance and the critical current of the Josephson junction, respectively. If the flowing current *I* is lower than *IJJ* there is no voltage drop at the inter-grain junction (*VJJ* = 0). Then, the total conductance *G* can be calculated from the condition 1/*G* = *dVPC dI* + *dVJJ dI* . We succeeded in simulating the conductance spectrum reported in Figure 4d assuming Δ = 0.55 meV, Z = 0.54 and α = 0.20. The fitting parameters related to the Josephson junctions were *RJJ* = 0.1 Ω and *IJJ* = 3.2 mA. We notice that these parameters are not completely free, being necessarily *RN* = *RPC* + *RJJ* and *RJJIJJ* < Δ. We remark here that neglecting the existence of a Josephson junction in series with the point contact would result in an over-estimation of the superconducting energy gap, because the measured voltage at which the conductance features are observed is larger than the real voltage *VPC* applied to the point contact junction (*V* = *VPC* + *VJJ*). We observe that the superconducting energy gap values obtained from the numerical fittings of most of the experimental spectra are in the range 0.85–0.95 meV that correspond to a ratio 2Δ/kBTC in the range 2.5–2.8, smaller than the BCS value (3.52). In the case of the conductance curve of Figure 4d, we estimate a superconducting energy even smaller (Δ = 0.55 meV, i.e., 2Δ/kBTC = 1.6). This may be an

indication of suppressed superconductivity on the probed surface, with the point contact experiments being sensitive to a thin surface layer of tens of nanometres. Consequently, the correct procedure for estimating the ratio 2Δ/kBTC would be to use the local critical temperature, which should be estimated based on the temperature evolution of the conductance spectra (not available in this experiment); this local Tc could be lower than the bulk sample Tc, giving an increased 2Δ/kBTC ratio. We note that in the case of d-wave symmetry, electrons injected along different directions may experience different pairing amplitudes. Consequently, the shape of the conductance curves does not depend only on the height Z of the potential barrier at the interface, but also on the direction of the current injection. In Figure 5a, we show conductance measurements performed in a different location of the sample. The two spectra are the result of two successive measurements, in which the second spectrum was measured after increasing the tip pressure on the surface to increase the barrier transparency. The lower spectrum shows a ZBCP with limited amplitude (about 1.2) and is reproduced by the extended BTK model by assuming Δ = 0.55 meV, Z = 0.71 and α = 0.46. However, the second spectrum (shifted for clarity) has a much higher and more narrow ZBCP with two relative maxima appearing at the side of the peak. To understand the evolution of the second conductance curve, we show in Figure 5b the expected behaviour of the spectra obtained by keeping fixed the parameters Δ = 0.55 meV and α = 0.46, varying the barrier strength Z in the range 0 < Z < 1. We notice that in this scenario the main effect of the Z parameter is on the ZBCP height. The appearance of further conductance features is explained considering that the increased pressure of the tip on the surface has a double effect: it helps to obtain a more transparent barrier (lower Z), while it favours the formation of an inter-grain junction. Indeed, the numerical simulation of the experimental spectrum is obtained with good agreement by assuming Δ = 0.55 meV, Z = 0.39, α = 0.29 and *RJJ* = 0.22 Ω, *IJJ* = 0.24 mA. In Figure 5c, we show the evolution of conductance spectra numerically calculated by fixing the parameters Δ = 0.55 meV, Z = 0.39, α = 0.29 and varying only *RJJ* in the range 0–0.42 Ω.

#### **4. Conclusions**

We investigated the superconducting properties of polycrystalline Pr1−xCexPt4Ge12 pellets, reporting the magnetic field phase diagram H-T for the Ce-doped compound with x = 0.07, as well as the point contact spectroscopy characterization of the parent compound PrPt4Ge12 (x = 0).

Interestingly, the irreversibility field line, found for Pr0.93Ce0.07Pt4Ge12, shows a scaling behaviour similar to high-temperature superconducting cuprates. The vortex activation energy was also evaluated at different applied magnetic fields. At magnetic fields lower than 0.4 T, the activation energy follows a power law of the type *<sup>H</sup>*−α, with the exponents <sup>α</sup> <sup>≈</sup> 0.5, which could indicate a collective pinning regime with a quasi-2D character.

For the compound with x=0 (PrPt4Ge12), we realized normal metal/superconductor nano-junctions, with lateral dimensions of few nanometres, by pushing a gold tip onto the surface of polycrystalline sample. Several conductance spectra were measured at low temperatures, showing zero bias conductance peak with variable amplitude, height and width. All experimental data for the PrPt4Ge12 sample were consistently interpreted in the framework of extended BTK theory. A small energy gap was observed in the range 0.55 meV–0.95 meV, indicating the possible formation of inter-grain Josephson junctions in series with the point contact.

**Author Contributions:** Conceptualization, F.G. and P.R.; Data curation, F.A., A.L. and A.D.B.; Investigation, P.R., A.N., G.G., A.L. and F.G.; Resources, L.S. and J.Z.; Supervision, F.G.; Writing—original draft, F.G., P.R., A.L. and A.N.; Writing—review & editing, F.G., P.R., A.D.B. and A.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Research and Development Program of China, No. 2017YFA0303104, Shanghai Municipal Science and Technology Major Project (Grant No.~2019SHZDZX01).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Current-Resistance E**ff**ects Inducing Nonlinear Fluctuation Mechanisms in Granular Aluminum Oxide Nanowires**

#### **Carlo Barone 1,2,3,\*, Hannes Rotzinger 4,5, Jan Nicolas Voss 4, Costantino Mauro 1, Yannick Schön 4, Alexey V. Ustinov 4,6,7 and Sergio Pagano 1,2,3**


Received: 17 February 2020; Accepted: 11 March 2020; Published: 14 March 2020

**Abstract:** The unusual superconducting properties of granular aluminum oxide have been recently investigated for application in quantum circuits. However, the intrinsic irregular structure of this material requires a good understanding of the transport mechanisms and, in particular, the effect of disorder, especially when patterned at the nanoscale level. In view of these aspects, electric transport and voltage fluctuations have been investigated on thin-film based granular aluminum oxide nanowires, in the normal state and at temperatures between 8 and 300 K. The nonlinear resistivity and two-level tunneling fluctuators have been observed. Regarding the nature of the noise processes, the experimental findings give a clear indication in favor of a dynamic random resistor network model, rather than the possible existence of a local ordering of magnetic origin. The identification of the charge carrier fluctuations in the normal state of granular aluminum oxide nanowires is very useful for improving the fabrication process and, therefore, reducing the possible sources of decoherence in the superconducting state, where quantum technologies that are based on these nanostructures should work.

**Keywords:** quantum electronics; noise spectroscopy; granular aluminum oxide; superconducting nanowires; current-resistance effects

#### **1. Introduction**

Granular superconductors have attracted great interest for their rich phase diagram and for practical advantages, such as increased critical field [1], critical temperature [2,3], and kinetic inductance [4]. In particular, granular aluminum oxide (AlOx) has recently proven to play a prominent role in quantum electronics and in quantum bit (qubit) design due to its low electric loss and high kinetic inductance [5]. As a matter of fact, the kinetic inductance of an AlOx wire is proportional to its normal state resistance and, therefore, can be orders of magnitude higher than its geometric inductance [6]. Such a super-inductance could be very useful for the realization of a novel type of superconducting quantum circuit with an impedance above the vacuum impedance, and it could easily be constructed from AlOx wires [5,7]. The robustness of the material and the simplicity of fabrication make them a promising candidate to replace conventional Josephson tunnel junction arrays. The main differences with Josephson junctions are observed when AlOx is patterned in the form of a nanowire. The nanoscale structure can be advantageous in terms of electric loss, and in the limit of a short constriction can even be used as a non-linear element for quantum circuit design [8].

Although granular aluminum seems to be a very promising alternative material for qubit architecture, it has a disordered nature, due to a large amount of amorphous aluminum oxide. Previous published structural and morphological studies suggest that the granularity on the size of grains, connectedness, and the intergrain spacing depend on the preparation method [9], which results in films with a disorder distributed on a scale that ranged from few units to few tens of nanometers [10,11]. These types of granular films have different electric transport mechanisms, whose understanding is important, especially for systems realized at the nanoscale level. In this respect, noise spectroscopy has already been used for the sensitive investigation of low-dimensional superconducting films [12] and two-dimensional (2D) oxide interfaces [13,14], allowing for a better understanding of the charge carriers kinetic processes. Moreover, from a technological point of view, the integration of nanowire-based elements with standard circuit processing can be seriously limited by the increase of the 1/f electric noise. Therefore, such a noise component has to be carefully studied in order to have a deep comprehension of the nature of fluctuation mechanisms.

In view of all these considerations, DC electric and magneto-transport measurements, as well as voltage-noise investigations, have been performed on AlOx nanowires. The observed phenomenology, although related to the normal state, deserves to be understood in details, since it could also have an impact on the superconducting state, where qubit-based technology operates.

#### **2. Materials and Methods**

All of the investigated wires, as well as the surrounding structures, were made from a single layer of oxidized aluminum. As a first step, a 20 nm thick layer of disordered oxidized (granular) aluminum was grown by reactive sputter deposition on a *c*-plane sapphire substrate in a controlled oxygen atmosphere, thus allowing for a wide range of sheet resistances. The mean sheet resistance of the AlOx film was determined to be about 2.2 kΩ. However, the sheet resistance can vary up to 50% for different film areas across a 20 mm by 20 mm chip due to inhomogeneities that are related to the sputtering process.

An electron beam lithography process was used to define the wire structures, while using a bilayer of hydrogen silsesquioxane (HSQ) and polymethylmethacrylate (PMMA) resist. Here, the high resolution of the HSQ and the ability to use an organic solvent to lift off the PMMA are beneficial for a clean interface and also to contact the AlOx structures. The PMMA serves as a sacrificial layer, it, however, also influences (limits) the achievable nanowire resolution. Two anisotropic dry etch (RIE) steps were performed to transfer the structures into the AlOx layer after the resist definition: oxygen RIE for the PMMA, and argon/chlorine RIE for the AlOx.

The left panel of Figure 1 shows a typical scanning electron microscope image of another nanowire sample that was fabricated with similar process conditions. Here, the wires have a thickness of 20 nm and a width of few tens of nanometers, in the order of 60 nm. It has been found that small variations in the process parameters can lead to a significant spread in the wire width and resolution of the above described bilayer process. We suspect that the limited temperature stability of the PMMA resist might be responsible for this and, however, further optimizations may be necessary for a better reproducibility of the nanowire shape. Several devices were measured, each of them made by two 150 μm × 150 μm contact pads connected by a narrow stripline, as shown in the right panel of Figure 1. Nanowires of different lengths are fabricated in the middle of each stripline. Devices, named "a", "b", "g", and "h" are the nanowires of length (*WL*) 250, 500, 750, and 1000 nm, respectively. Devices named "i", instead, are microstrips 5 μm-long and, nominally, 100 nm-wide.

**Figure 1.** (Left panel) Colorized optical picture of a device investigated. The colorized scanning electron microscope image of a typical nanowire is shown in the enlargement. The wire is 20 nm thick and it has a typical width up to 60 nm. (Right panel) Layout of the chip. The measured nanowire samples are evidenced by the circles. The contacts are made by 50 μm thick with aluminum wire by ultrasonic bonding.

A simple two-probe measurement method was employed, since all of the measured resistances were much higher than that of the contact wires to investigate the samples transport and noise properties. Additionally, the wires and contacts contribution to the overall electric noise was verified to be negligible [15]. The value of resistance measured involves the whole stripline, not only the submicron part. This implies that, for devices "a", "b", "g", and "h", there are additionally 10 squares contributions to consider for the total resistance. Electric transport measurements were performed in a temperature stabilized closed-cycle refrigerator, between 300 and 8 K. A dipole electromagnet, type 3470 from GMW Associates, generated an external magnetic field (maximum value of 1500 Gauss). The field is applied in the plane of the sample and parallel to the nanowire length which coincides with the bias current direction. Low-noise DC and AC electronic bias and readout have been used [15,16], thus giving very low statistical errors on the experimental data points (in the following plots, the error bars are usually smaller than the points dimension). The range of explored bias currents depends on the sample resistance, being limited by the overall voltage drop that is available at the current source (100 V).

#### **3. Results**

#### *3.1. DC Electrical Transport Measurements*

Figure 2 shows the temperature dependence of the resistance of all the investigated devices. Differently from the previously analyzed case of AlOx thin films [17], the nanowires are characterized by a monotonic resistivity decrease by increasing the temperature between 8 and 300 K. This decreasing behavior, which is not unusual in high resistive samples, seems to be more pronounced at low temperatures, although the resistance value does not scale with the wire length, probably due to differences of the average wire width and resistivity across the nanowire devices. By using the microstrip (nominally 100 nm-wide) as a reference, it is possible to estimate the room-temperature resistivity of the nanowires <sup>ρ</sup>*RT* <sup>≈</sup> 5.2 <sup>×</sup> 103 <sup>μ</sup><sup>Ω</sup> cm. For these values of <sup>ρ</sup>*RT* (higher than <sup>4</sup> <sup>×</sup> 103 <sup>μ</sup><sup>Ω</sup> cm), the granular aluminum specimens are characterized by a resistance that rises monotonically as the temperature is reduced, confirming the experimental behavior shown in Figure 2, as reported in literature [18,19]. Conversely, for <sup>ρ</sup>*RT* values between 1 <sup>×</sup> <sup>10</sup><sup>3</sup> and 4 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>μ</sup><sup>Ω</sup> cm there is an intermediate regime in which the resistance changes very little with temperature and it shows a logarithmic increase at low temperatures (Kondo regime). This is the case of the previously analyzed AlOx thin films [17], which, therefore, have a resistivity temperature dependence different from the nanowires studied here. It is also important to underline that the strictly metallic regime is usually characterized by values of ρ*RT* that are lower than 103 μΩ cm [20].

**Figure 2.** Resistance versus temperature plots. The data refer to five investigated devices. The solid lines are obtained by interpolating the experimental points and they are guides for the eyes.

All of the devices exhibit almost linear current voltage characteristics, with a small curvature that can be observed for bias current above few μA. By analyzing in detail the bias current dependence of the static resistance *R*, defined as the ratio *V*/*I*, the wire resistance shows a tendency to decrease for increasing bias. The repeatability of this effect has been verified in all of the nanostructures analyzed, by performing several current cycles from 0 to 80 μA. Consequently, the reversible resistance change is independently present from the temperature, but is more pronounced at low temperatures, as shown in Figure 3. This has been also observed on manganite compounds and it is known as the negative current-resistance (CR) effect [21–24]. Several mechanisms could be considered to explain the experimental evidence of reversible CR effects in AlOx nanowires. In analogy with the case of magnetic multilayered structures [25,26], the CR effect could be attributed to a local ordering that is produced by a current-induced local magnetic field. Alternatively, the CR effect could be ascribed to a percolative process between metallic AlOx clusters that are separated by an insulating matrix [27,28]. In this framework, an important role is played by the tunneling barrier between conducting grains and by the shape of such a barrier, which the applied current can change.

**Figure 3.** Static resistance versus bias current at different temperatures. The data refer to the nanowires 1000 nm-long at *T* = 300 K (red dots) and at *T* = 9 K (blue triangles). The small jumps visible in the traces are reproducible and are, most probably, due to different current paths in the sample.

Overall, DC measurements alone are not able to distinguish the possible mechanism that is responsible for the CR effect. Therefore, more sensitive investigations, such as electric noise spectroscopy, have been performed, in order to acquire additional information.

#### *3.2. Voltage-Noise Spectral Density Measurements*

The first indication on the nature of the fluctuation processes is obtained by analyzing the frequency composition of the voltage-noise spectral density *SV*. In Figure 4, typical noise spectral traces are reported for the four nanowires under test at different bias current values and temperatures. The peaks are due to external sources and they should not be considered. The relevant information is in the background curves. *SV* is characterized by a frequency dependence that varies on temperature. In particular, the presence of a 1/*f* noise component followed by a constant amplitude spectrum is observed at temperatures down to 150 K (Figure 4a,b). Conversely, below 150 K, Lorentzian components, with a 1/*f* <sup>2</sup> dependence of *SV*, are clearly visible (Figure 4c,d). This behavior is typically associated with physical fluctuation mechanisms that are ascribed to two-level tunneling fluctuators (TLTFs) [29,30]. According to this interpretation, groups of atoms can occupy two configurations; therefore, their energy can be represented in the form of two potential wells that are separated by a barrier. In general, the wells are asymmetric and the atoms can tunnel from one well to the other with a transition rate that is strongly temperature-dependent. By indicating with τ the minimum relaxation time of such processes, for frequencies *f* 1/2πτ the main contribution to *SV* is given by the 1/*f* noise, while, for frequencies *f* 1/2πτ, the spectral density is dominated by Lorentzian noise. By lowering the temperature, τ increases and the cut-off frequency, above which the 1/*f* <sup>2</sup> noise component dominates, decreases, being visible inside the frequency bandwidth of acquisition (see, for details, all of the low-temperature spectra in Figure 4).

**Figure 4.** Voltage-noise spectra. The frequency dependence of *SV* is shown for all of the investigated nanodevices for two bias currents and four different temperatures: 300 K (**a**), 150 K (**b**), 100 K (**c**), and 9 K (**d**). The black trace, which is the lowest spectral trace for each temperature, shows the unbiased noise, which is composed by the instrumental, the external, and the sample Johnson noise.

Spontaneous transitions between the levels of the TLTFs can lead to fluctuations of macroscopic quantities, such as resistance [29]. In the ohmic region, where resistance does not depend on the bias current, resistance fluctuations are characterized by a *SV* versus *I* <sup>2</sup> dependence. This means that the overall noise level, being defined as:

$$\text{Noise Level} = \frac{1}{V^2} \int\_{f\_{\text{min}}}^{f\_{\text{max}}} S\_V(f) df\_\prime \tag{1}$$

with [*fmin* ; *fmax*] as the frequency bandwidth, should be constant as a function of the bias current. A confirmation of this behavior is found in the plots of Figure 5, where the current dependence of the noise level is shown for three devices at different temperatures, for low levels of the applied bias. In Figure 5, it is also evident a noise level increase when the temperature is reduced below 150 K. This is due to the activation of additional low-temperature noise sources, being identified in the framework of the TLTFs model and evidenced in the spectral traces as Lorentzian components (see Figure 4c,d, for details).

**Figure 5.** Current dependence of the nanodevices noise level in the low-bias region. The noise level dependence on the bias current, as evaluated from Equation (1), is shown for currents up to 200 nA and for devices of different lengths: 250 nm (**a**), 500 nm (**b**), and 750 nm (**c**). The investigated temperature range is from 8 to 300 K.

However, for bias currents where CR effects are present, the noise level also changes. In particular, the data in Figure 6 show a noise level reduction at large levels of the applied bias. This behavior, which is evident both at room temperature and at 9 K (black squares and magenta diamonds in Figure 6, respectively), indicates, therefore, that the local ordering that is attributed to CR effects gives origin to a decrease of noise fluctuators. The decrease of the noise level, both for decreasing temperature and for increasing DC bias current, seems to be incompatible with a Joule heating effect. In fact, Joule heating would imply an increase of the nanowire temperature with the bias current, and this should also result in an increased noise level, since it is enhanced by temperature. Instead, Figure 6 shows clearly that the noise level decreases at larger bias currents.

**Figure 6.** Current dependence of the nanodevices noise level in the full-bias region. The overall noise level, as evaluated with Equation (1), is shown as a function of the bias current for the nanowire 1000 nm-long. The temperatures of 300 and 9 K are reported.

The nature of the noise decrease and, consequently, the nature of reversible-type CR phenomena is at this stage unclear, but a better insight could be obtained by looking at its dependence on external parameters, such as external magnetic field or electric field.

#### **4. Discussion**

As said before, in other materials, it has been found that possible mechanisms producing the current-resistance effects and, at the same time, the noise level reduction could be of magnetic origin or could be due to the percolative process. Whether one of these models can be associated to granular AlOx nanowires has to be investigated.

Magnetic resistance fluctuations are usually characterized by a noise level increase in the low-temperature region. This type of fluctuation process, which is already found for granular aluminum oxide thin films, is ascribed to Kondo-like conductivity when the noise level, especially at low temperatures, decreases with an external magnetic field of the order of one thousand Gauss [17,31]. In the case of AlOx nanowires, no magnetic dependence is observed in both DC and AC measurements, even if a noise enhancement is evident by lowering the temperature (for details, see Figure 7). By applying an external magnetic field *H* = 1000 Gauss, the magnitude of the CR effect is unaltered (Figure 7a). Moreover, the normalized 1/*f* noise power (*f*·*SV*) does not change with *H* for all of the bias currents used (Figure 7b). This experimental finding indicates that, contrarily to what observed for the films, in AlOx nanowires, a local ordering of magnetic origin does not produce the noise reduction and, therefore, does not seem to be responsible of the measured CR effect. The different behavior that was observed between granular aluminum film and nanowires could be related to

differences in morphology and grain structure in the two cases, which also show different resistivity, as discussed above.

**Figure 7.** Magnetic field behavior of DC and AC data at 9 K. The current dependence of (**a**) the resistance and of (**b**) the normalized 1/f noise component of device 750 nm-long is plotted in absence and in presence of a magnetic field of 1000 Gauss (black and red curves, respectively).

Alternative to the magnetic effect, a nonlinear transport regime can also result from a random resistor network (RRN) model and, consequently, percolation processes [32–34]. These conduction mechanisms have been already considered in granular films that are similar to the ones investigated here [9], being realized by reactive sputter deposition and with a grain size in the order of 3–4 nm [5]. In particular, a reasonable description of nonlinear conduction is given by the dynamic random resistor network (DRRN) model, which considers a network that consists of conducting and insulating bonds [35]. For sufficiently low applied electric field across the network, the current flows only through the backbone of the percolating system, so that the conduction is linear. As the applied bias is increased, some of the insulating bonds of the network will experience an electric field exceeding a critical field (*Ec*), above which they become conducting [35]. The macroscopic conduction becomes nonlinear when insulating bonds start to be conducting. This is similar to what is shown in Figure 8a, where the nanowire resistance is plotted versus the local electric field *E*, being computed as the voltage difference divided by the nanowire length. In this framework, the conducting clusters forming the network backbone are separated by insulating regions of very small widths. Through such regions, it is expected that tunneling or hopping conduction will take place, thus providing extra paths for electrons. The appearance of extra parallel conduction channels with the increase in bias, which add to an existing network, does lead to a decrease in the total noise of the network, which is consistent with that observed in Figure 8b.

The noise reduction above the onset of nonlinearity has been already observed in conductor-insulator mixtures, such as the carbon-wax system, and, in terms of the DRRN model, has been explained as due to an increase in the total number of fluctuators or, equivalently, in the system size from the fluctuation point of view [36]. A better display of the nonlinear transport in AlOx nanowires can be made by computing the negative reversible CR effect, as:

$$CR\,effect\left(\%\right) = \frac{R(E) - R\_0}{R\_0},\tag{2}$$

where *R*(*E*) is the resistance at a fixed electric field and *R0* is the resistance value in the linear region. The quantity that is defined with Equation (2) is plotted in Figure 8c as a function of *E* for three nanodevices of different length at the temperature of 9 K. Here, it is possible to identify a nonlinearity

threshold at a critical electric field *Ec* <sup>≈</sup> 7.2 <sup>×</sup> 10<sup>4</sup> V m<sup>−</sup>1. Similarly, the noise reduction in percentage can be computed as:

$$\text{NR } (\%) = \frac{\text{NL}(E) - \text{NL}\_0}{\text{NL}\_0},\tag{3}$$

where *NL*(*E*) is the noise level value at a fixed electric field and *NL0* is the noise level value in the linear region. *NR* is shown in Figure 8d for the same nanodevices of Figure 8c and it confirms the existence of a threshold field. This effect has been well documented at 9 K, and is present also at higher temperatures, see Figure S1 (Supplementary Materials). However, further analysis of the temperature dependence of the high current behavior of the nanowires has not been performed, being beyond the scope of this work, as the main applications of AlOx nanowires are at cryogenic temperatures. The noise reduction effect is much more evident than the CR effect, revealing a greater sensitivity of the noise spectroscopy with respect to the resistivity measurements. Moreover, the evidence of a saturated noise level at high bias values finds strict analogy with the saturated transport properties that were observed in carbon-wax composites and described in terms of the DRRN theoretical approach [36]. This gives a further indication on the applicability of this model to explain the reversible nonlinear transport regime in AlOx nanowires.

**Figure 8.** Electric field behavior of DC and AC data at 9 K. For the nanowires 500 nm-long (red circles), 750 nm-long (blue triangles), and 1000 nm-long (yellow diamonds), it is shown the electric field dependence of the: (**a**) static resistance, (**b**) noise level evaluated with Equation (1), (**c**) current-resistance effect evaluated with Equation (2), and (**d**) noise reduction evaluated with Equation (3).

Although a nonlinear conducting behavior is expected for metal-dielectric composite nanostructures falling around the percolation threshold [37,38], the strong effect on the magnitude of voltage-noise is a peculiar feature of AlOx nanowires, not being observed in common nonlinear composites dielectrics. This property of the charge carriers fluctuation mechanisms has a direct impact on the normal state conductance and, as a consequence, on the kinetic inductance of granular aluminum nanodevices, a key requirement for their use in superconducting quantum circuits.

#### **5. Conclusions**

The electric transport and voltage-noise properties of granular aluminum oxide nanowires have been investigated for nanodevices of different width and length and in a temperature range from 8 to 300 K. Differently from thin films, which show a linear ohmic regime, nanowires are characterized by a resistivity reduction with increasing bias. In cycled devices, the negative resistance change is reversible. In this regime, above a critical electric field, the nonlinear current-voltage characteristics are also associated with a strong noise level reduction, whose amplitude exceeds that observed on DC resistance of more than one order of magnitude. This nonlinear electric transport behavior of the nanowires does not seem to have a magnetic origin. Moreover, the magnetic field does not affect the noise amplitude, contrarily to what found in thin films with lower resistivity and at low temperatures. The nonlinear resistance could be explained, instead, in terms of a dynamic random resistor network model, usually considered as a variant of the standard random percolation mechanism. If a DRRN model describes the physics of the AlOx nanowires, then it is expected that noise sources can be activated whenever there is an electric field in the material higher than a threshold critical value, and large enough fields can suppress it. In the superconducting state, one could expect that *E* field is zero inside the nanowire and no noise sources are activated. However, in the presence of RF fields, this is no longer true and the noise sources could play an important role in the dynamics of the circuit attached to these devices. Therefore, the identification and reduction of the charge carrier fluctuations in the normal state, by e.g. improving the fabrication process, could be important for reducing the possible sources of decoherence in the superconducting state.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2079-4991/10/3/524/s1, Figure S1: Electric field behavior of resistance and noise data at different temperatures.

**Author Contributions:** C.B., H.R., A.V.U. and S.P. conceived and designed the experiments; C.B., C.M. and S.P. performed the experiments; C.B., H.R. and J.N.V. analyzed the data; J.N.V., C.M., Y.S. and A.V.U. contributed materials/analysis tools; C.B., H.R. and S.P. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank S. Abate of CNR-SPIN Salerno (Italy) for his technical support. A.V.U. acknowledges partial support from the Ministry of Education and Science of the Russian Federation in the framework of Increase Competitiveness Program of the National University of Science and Technology MISIS (Contract No. K2-2017-081). C.B. and S.P. acknowledge partial support from University of Salerno through grants FARB17PAGAN and FARB18CAVAL, and from INFN through experiment FEEL.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Iron-Based Superconducting Nanowires: Electric Transport and Voltage-Noise Properties**

**Sergio Pagano 1,2,3,\*, Nadia Martucciello 2,3, Emanuele Enrico 4, Eugenio Monticone 4, Kazumasa Iida <sup>5</sup> and Carlo Barone 1,2,3**


Received: 26 March 2020; Accepted: 20 April 2020; Published: 30 April 2020

**Abstract:** The discovery of iron-based superconductors paved the way for advanced possible applications, mostly in high magnetic fields, but also in electronics. Among superconductive devices, nanowire detectors have raised a large interest in recent years, due to their ability to detect a single photon in the visible and infrared (IR) spectral region. Although not yet optimal for single-photon detection, iron-based superconducting nanowire detectors would bring clear advantages due to their high operating temperature, also possibly profiting of other peculiar material properties. However, there are several challenges yet to be overcome, regarding mainly: fabrication of ultra-thin films, appropriate passivation techniques, optimization of nano-patterning, and high-quality electrical contacts. Test nanowire structures, made by ultra-thin films of Co-doped BaFe2As2, have been fabricated and characterized in their transport and intrinsic noise properties. The results on the realized nanostructures show good properties in terms of material resistivity and critical current. Details on the fabrication and low temperature characterization of the realized nanodevices are presented, together with a study of possible degradation phenomena induced by ageing effects.

**Keywords:** iron-based superconductors; nanowires; single-photon detectors

#### **1. Introduction**

Most widespread applications of superconducting nanowires regard their use as single-photon detectors, due to their ability in detecting single photons in the visible and IR spectral region [1,2]. Moreover, their peculiar physics has brought to the development of novel electronic cryodevices, such as pulse discriminators [3,4], logic gates [5], and also memory elements [6].

The interest in superconducting nanowires has also recently increased, due to their possible application in quantum technologies, including quantum sensing and computing. This has been highlighted in the case of YBa2Cu3O7-x (YBCO) nanowires with phase-slip dynamics, where evidence of energy-level quantization in the nanowires has been reported [7]. Moreover, it has also been shown that the absorption of a single photon changes the quantum state of the nanowire, an important result for the development of single-photon detectors with high operating temperature and superior temporal resolution [7].

Although traditional superconducting nanowire single-photon detectors (SNSPDs) have the advantage of offering single-photon sensitivity, combined with low dark count rates [8], low jitter [9], short recovery times, and free-running operation [2], one drawback is their low operating temperature. This is essentially due to the fact that current SNSPDs are fabricated using mostly conventional low-temperature superconductors, such as NbN and WSi [10,11], in order to achieve high sensitivity and to simplify nanofabrication processes. Several efforts have been made to realize nanowires with high critical temperature (*Tc*) materials, including MgB2 and cuprate superconductors [12–15]. The recently discovered iron-based superconductors have attracted great interest to explore their potentialities in the field of large-scale current transport [16,17] and in microelectronics or nanoelectronics applications [18,19]. These compounds could pave a new way to the fabrication of superconducting nanowires, also profiting of intrinsic material properties to improve detection performances (in particular, speed and efficiency). In this respect some preliminary results are reported in [20]. However, there are several issues to be taken into account that are mainly related to the occurrence of non-hysteretic current-voltage characteristics (no switching), to the difficulty in fabricating ultra-thin films (dead layer problem of high-*Tc* compounds), and to the easy surface degradation.

In order to overcome these difficulties and address these challenges, we have realized and tested a number of nanowires made of iron-based superconductors. A detailed investigation of their electric transport properties is reported in this work. In addition, a study of the charge carrier fluctuation effects, in a wide temperature range, has been carried out, in order to get useful information on the physics of the nanowires for the optimization of nanopatterning processes and for the realization of high-quality electrical contacts. To our knowledge, there are no other noise measurements on iron-based nanowires reported in the scientific literature. By means of noise spectroscopy, moreover, a deeper understanding of physical mechanisms and of ageing-induced degradation effects can be obtained, as already demonstrated for graphene [21,22], for iron-chalcogenide [23], and for the same Co-doped BaFe2As2 superconductors [24] used here. The reported experimental results, together with theoretical interpretations, may be useful in view of the fabrication of more performant and usable nanodevices.

#### **2. Materials and Methods**

A necessary starting point to realize a nanowire detector is a very thin superconducting film, ideally with a thickness of the order of the material coherence length, in order to maintain the superconductive state and maximize the photon energy sensitivity. However, the growth mechanism of high-*Tc* superconductors is complex and often results in a "dead" layer at the interface with the substrate. Therefore, a compromise between thickness and film quality has to be achieved. In this respect, Co-doped BaFe2As2 superconducting thin films, with a thickness of about 20 nm, were grown, by using a pulsed laser deposition (PLD) technique, on 0.5-mm-thick CaF2 (001) substrates. The choice of the Co doping level is done in order to optimize the superconducting transition temperature of the material [25]. This resulted in a nominal composition of the PLD target of Ba:Fe:Co:As = 1:1.84:0.16:2. The fabrication parameters, such as laser repetition frequency (7 Hz) and growth temperature (700 ◦C) [26], together with the choice of CaF2 substrate, were selected in order to optimize the *Tc* of the films, as well as other material properties [27]. The high phase purity of the target material and of the samples was observed with X-ray diffraction (Malvern Panalytical Ltd., Malvern, UK) and by using transmission electron microscopy (JEOL Ltd., Tokyo, Japan) the absence of appreciable defects was also verified [26]. Moreover, to reduce possible degradation effects due to ageing, an in-situ passivation process was employed, by depositing a thin layer of magnesium aluminate (MgAl2O4) spinel (MAS). Its high resistance against chemical attack, high thermal shock resistance and compatibility with a large variety of metals are important properties that make MAS very advantageous as a protective cap layer [28].

Subsequently, the superconducting films were patterned in the shape of several nanowires using a combination of optical and e-beam lithography and ion beam etching. The design used for the geometry definition, shown in Figure 1a, consists of a series of nanowires with a nominal width and length of 500 nm and 5 μm, respectively, each connected to 1 mm × 0.5 mm bonding pads. A 100 nm thick resist layer, type ma-N 2401, was used for reproducing the design on the samples (Figure 1b).

**Figure 1.** Microphotographs and images representing the different steps of the nanowires fabrication process: from the used mask design (**a**) to the finished 1 cm2 chip, after electronic and optical lithography (**b**–**f**). The devices, whose experimental results are here reported and discussed, are evidenced with colorized circles.

A subsequent ion milling process (Figure 1c) defined the nanowire geometry. After etching, due to nanolithography process imperfections, the final devices realized showed different lengths, ranging from 1.8 (Figure 1d) to 3.4 μm (Figure 1e), while keeping the same width. Finally, to reduce the contact resistance, a thin layer of Ti/Au (5/50 nm) was deposited after Ar ion cleaning, on the contact pads area, obtaining the finished chip consisting of 27 separated devices. Among these, the ones here investigated are evidenced with circles of different colors and refer to pristine (black circle) and to two years-aged (red and green circles) nanostructures (see Figure 1f, for details).

All the electrical characterizations were performed in a closed cryocooler system, characterized by an operation temperature range from 8 to 325 K, with a temperature stabilization better of 0.2 K. Low-noise DC bias, and DC and AC readout electronics were used to record the sample electrical response, as shown Figure 2 [29]. In particular, the measurements were performed by using a two-probe technique (green box of Figure 2). A low-pass passive filter was inserted at the input of the bias current source, with a user variable series resistance *RS*, much larger than the sample resistance *RM*, and a cutoff frequency of few Hz (blue box of Figure 2). The output AC voltage was amplified by a low-noise instrumentation amplifier, model AD8429 having 1 nV/ √ Hz noise level, (red box of Figure 2), and its spectral analysis was done by a dynamic signal analyzer model HP35670A. The DC voltage was amplified by a low noise instrumentation amplifier, model AD8221, and recorded by a digital voltmeter. The absence of unwanted contact noise contribution was verified by resorting to a specific procedure [30], whose validity has been already tested in the case of superconductors [31], innovative carbon nanotube and photovoltaic devices [32,33], and magnetic compounds [34,35]. The lowest noise level measured, corresponding to the instrumental background, was 1.71 <sup>×</sup> <sup>10</sup>−<sup>18</sup> <sup>V</sup>2/Hz.

**Figure 2.** Schematic of the bias and readout electronics. The colored dashed boxes enclose the bias circuit (blue), the device under test (green), and the instrumentation amplifier (red). See text for details.

#### **3. Results**

#### *3.1. DC Electrical Transport Measurements*

The temperature dependence of the resistance of Co-doped BaFe2As2 nanowires was very similar to what already found for thin films [27]. In particular, the data in Figure 3a show a metallic behavior for temperatures higher than 100 K, both for the thin film (blue circles) and for a nanowire (black squares). A minimum of the resistivity was observed around 100 K, with an upturn down to the onset of superconductivity (*Tonset <sup>c</sup>* ), which was about 30 K for the thin film and 24 K for the nanowire. The difference in critical temperature was due to different film thickness, much larger in the case of the thin film. The resistance increase, observed by decreasing the temperature from 100 K to *Tonset <sup>c</sup>* , can be attributed to the occurrence of localization effects [36]. In this respect, a detailed analysis, reported in [37], limits to variable-range hopping (VRH) or to weak-localization (WL) effects the possible explanation of the resistivity behavior. The DC measurements alone were not able to distinguish between the two mechanisms and, therefore, additional experimental investigations, such as noise spectroscopy, are necessary.

The results obtained by studying the fluctuation processes are shown in the following and will be useful to clarify this issue. The most evident difference between thin films and nanowires concerns the value of *Tc*. More in details, for the single nanowire the superconducting transition, defined at 50% of the normal-state resistance *RN*, occurred at around 16 K and was significantly lower than that of the thin film. This fact, as already discussed above, was due to the choice of fabricating a very thin film for the nanowire sample, in order to boost the device sensitivity. The observed low *Tc* was, however, larger than that of standard NbN nanowires (about 11 K). Moreover, as shown in Figure 3b, where the normalized resistance was reported in the superconducting transition region, after two years of ageing the nanodevices maintained comparable critical temperatures. This finding gives the important indication that iron-based nanowires, realized with the passivation process described above, do not show degradation of their superconducting properties after up to two years of storage at room temperature and normal atmosphere. This result might be very useful in technological applications requiring long-time operation and stable performances, as also strengthened by an evident homogeneous distribution of the critical temperatures all over the chip area. Such a feature is confirmed by the experimental data of Figure 3b, which refer to differently positioned nanowires.

**Figure 3.** Resistance versus temperature plots. (**a**) The data refer to the resistance of thin films (blue circles) and of nanowires (black squares), normalized to their room temperature values. (**b**) An enlargement of the superconducting transition region is shown, where the resistances, normalized to their normal state values at *T* = 26 K, are reported for pristine and two years-aged nanostructures.

All the realized nanowires have non-hysteretic current–voltage (*I–V*) curves down to a temperature of 8 K, as shown in Figure 4a. Although the absence of degradation is a positive feature, the absence of hysteresis in *I–V* curves in the superconducting state is, instead, a negative point, when aiming at developing SNSPDs. It cannot be excluded that, by further lowering the temperature, a hysteretic *I–V* curve could be observed. However, as our aim was to realize high operating temperature nanowires, this is not of interest for this work. A critical current *Ic* of the order of 10 μA is measured at 8 K, as shown in the inset of Figure 4a. A non-hysteretic behavior was also reported in [20] for iron-based superconductors and in [38,39] in the case of YBCO. Some hysteresis was observed for a YBCO microbridge (not a nanowire) at *T* = 4K[40]. Compared to the reported results, the nanowires presented here have, however, the advantage of showing higher resistance values, similar to that of the low-*Tc* nanowires. Figure 4b shows the DC current dependence of the nanowire differential resistance *RD* for various temperatures. A gradual change of *RD* can be observed in the temperature range across the superconducting transition up to the normal resistance value, of about 1200 Ω. In the inset are clearly visible abrupt changes in *RD* occurring at bias current of about 0.6–0.7 mA, both for positive and negative values. These are most probably due to the transition to the normal state of the large area structures (see Figure 1b). Overall, it is evident that iron-based nanowires do not yet show characteristics optimal for light detection. In order to get more information on the transport processes, noise measurement have been employed.

#### *3.2. Voltage–Noise Spectral Density Measurements*

The basic properties of random data in the time domain, known as noise or fluctuations, are described by the autocorrelation function. In the frequency domain, information on the noise mechanisms in action are obtained by the spectral density function *S*x, which is the Fourier transform of the autocorrelation function, as derived from the Wiener–Khintchine theorem [41]. In condensed matter physics, the noise quantity usually measured is the spectral density of voltage fluctuations, whose frequency composition *SV(f)* gives indications on the different electric noise contributions. In particular, the most common types of low-frequency noises are: (I) Johnson or thermal noise, that is generated by the thermal agitation of the charge carriers inside an electrical conductor at equilibrium and is defined as *SV* = *4kBTR* (with *kB* the Boltzmann constant, *T* the temperature, and *R* the electrical resistance); (II) the shot noise, that is originated from the discrete nature of the electric charge and is defined as *SV* = *2eIR2* (with *e* the electron charge, and *I* the bias current); and (III) the 1/*f* or flicker noise, that is related to the ensemble average of thermally activated two levels random resistance fluctuations,

due to intrinsic material aspects, such as grain boundaries, defects, etc. [41]. The noise of type (I) and (II) has frequency-independent spectral density amplitude ("white noise"), while noise type (III) has a clear frequency dependence ("colored noise").

**Figure 4.** Current–voltage (*I–V*) characteristics. (**a**) The different curves refer to temperatures from 8 to 290 K. An enlargement of the *I–V* curve at the temperature of 8 K is shown in the inset. (**b**) The DC current dependence of the differential resistance *RD* is reported at different temperatures in the log-scale. In the inset, the full current range, spanning from negative to positive values, is shown using a linear-scale.

Both white and 1/*f* noise components were evident in the frequency dependence of *SV* for the iron-based nanowires investigated here, as shown in Figure 5 for two devices (nanowire #1 and nanowire #2 indicated in Figure 1f) and at the temperature of 290 K. More in details, the overall experimental noise spectral density can be modeled as:

$$\mathcal{S}\_V(I, T) = \frac{\mathcal{K}(I, T)}{f} + \mathcal{S}\_0 + \mathcal{S}\_B \tag{1}$$

where *K(I,T)* is the 1/*f* noise amplitude, dependent on bias current and temperature; *S*<sup>0</sup> is the "white noise" amplitude, which also may be bias and temperature dependent; *SB* is the instrumental background noise amplitude, in our case estimated to be a constant value of 1.71 <sup>×</sup> 10−<sup>18</sup> V2/Hz (see Section 2). All the physical information on the transport processes resides in*K* and *S*0, whose analysis in temperature and in bias current can be very useful for understanding the dynamic processes of the charge carriers. It is worth noting that a similar approach has been used to study nanowires made of granular aluminum oxide [42].

The experimental data can be interpreted in terms of the model described in Equation (1) to obtain the temperature and bias current dependence of *K* and *S*0. In Figure 6a,b are shown the DC current dependencies of the 1/*f* and "white noise" amplitudes, respectively, at various temperatures near *Tc*. Figure 6a shows a noise peak occurring at both low temperatures and DC currents, with a characteristic slope of two, and a noise growth occurring at high temperatures and currents, with a characteristic slope of 1.4. Figure 6b shows a large noise peak at low currents and temperatures, in correspondence of the superconducting transition of the nanowire, and a noise increase at larger temperatures, almost current independent.

**Figure 5.** Voltage-spectral density traces. The frequency dependence of *SV* is shown for two different devices of the same chip, device #1 (**a**) and device #2 (**b**), at room temperature (290 K) and for bias currents ranging from 0 to 0.8 mA.

**Figure 6.** Electric noise amplitudes behavior. The dependencies on bias current, and at several temperatures, are shown for the 1/*f* (**a**) and for the "white" (**b**) noise components, respectively.

To better appreciate the dependence of the noise on temperature and bias current, the experimental data can be shown in three-dimensional (3D) plots. In Figure 7a the 1/*f* noise amplitude *K* is plotted, with a linear temperature and bias current scale, and three noise peaks, evidenced by arrows, are clearly visible. The peaks indicated by the green and red arrows are expected and are due to percolation fluctuations occurring in the transition region of the nanowire (red arrow) and of the large area region (green arrow), respectively, where the superconductor is a mixture of normal metallic and superconducting phases [41,43,44]. On the contrary, the strong increase of the 1/*f* noise at low temperatures and high currents (blue arrow in Figure 7a) cannot be explained by reported theoretical interpretations and deserves more investigation in the next future.

**Figure 7.** Amplitude of the 1/*f* noise component. The temperature and bias current dependencies of the 1/*f* noise are shown by using linear-axes scales (**a**) and logarithmic-axes scales (**b**) for the noise amplitude and for the bias current. The 1/*f* noise peaks at different temperatures are indicated with colorized arrows in panel (**a**). The bias current slopes are explicitly reported in panel (**b**).

The same behavior of the 1/*f* noise amplitude is shown in Figure 7b by using a linear temperature x-scale and a logarithmic bias current and noise scale. This choice allows evidencing the power law dependence of the noise amplitude on bias current. The figure shows a (standard) quadratic dependence at low temperatures and below *Ic*. This is expected according to a simple percolation model describing the full superconducting regime [41,43,44]. Unexpectedly, a current slope of 1.4 is observed for temperatures up to *Tonset <sup>c</sup>* and for higher currents than *Ic*. This unusual finding may be connected with the unusual 1/*f* noise peak and could reveal the activation of additional noise sources that occurs in a region where the nanowire is already in the normal state but the thin-film part is still superconducting. This "intermediate" regime is characterized by a strong non-uniform conductivity, which could be responsible for both the high noise level and its non-standard bias current dependence. Additionally, the graph of Figure 7b shows that, at temperatures above *Tonset <sup>c</sup>* , where the nanowire normal resistivity is characterized by a metal-insulator behavior, the current slope of the 1/*f* noise is again quadratic. This finding excludes the presence of WL effects, which appear as a linear current dependence of the 1/*f* noise component [45,46], while indicates the occurrence of VRH conductivity in the insulating regime and of resistance fluctuation processes in the metallic region [37,41].

The temperature and bias current dependence of the "white noise" *S*<sup>0</sup> amplitude, shown in Figure 8, has a peculiar behavior. *S*<sup>0</sup> has a maximum at low currents (centered near *Ic*) and low temperatures (when the nanowire is superconducting), and seems to peak when the bias current reaches the critical one of the nanowire. Moreover, *S*<sup>0</sup> has a smaller peak near *Tonset <sup>c</sup>* , where the large area regions start their superconducting transition. At larger current and higher temperature, the "white noise" amplitude tends to become current independent and to approach the standard Johnson noise value. These experimental findings are new and, to our knowledge, not been reported in literature.

**Figure 8.** Amplitude of the "white noise" component. The temperature and bias current dependencies of the "white noise" are shown by using linear-axes scales for the noise amplitude and for the temperature, and logarithmic-axis scale for the bias current. The instrumental background noise is also shown as a gray floor.

#### **4. Conclusions**

Nanowire structures were realized using ultra-thin films of Co-doped BaFe2As2. The sample fabrication and nanolithography techniques were optimized to preserve the superconducting properties of the material. The resulting critical temperature, about 16 K, was reduced of only 1 K after two years of room temperature and ambient atmosphere storage, showing that degradation effects due to ageing were not relevant. Due to their high superconducting transition temperature, iron-based nanowires represent an interesting alternative to NbN-based nanowires, whose critical temperature was about 10 K in ultra-thin film devices, for application as single-photon detectors. However, additional improvements in the fabrication technology are needed in order to achieve hysteretic *I–V* characteristic, a necessary ingredient for SSPD device. In this respect, a deep investigation of the transport mechanism in these devices is important, and it has been shown that noise spectroscopy can provide useful information. The study of the noise sources and their effect on dark counts, combined with the photo-response analysis, is currently in progress with the aim of developing useful single photon detectors operating at temperature above 10 K.

**Author Contributions:** S.P., N.M. and C.B. conceived and designed the experiments; S.P. and C.B. performed the experiments; S.P., N.M., K.I. and C.B. analyzed the data; E.E., E.M. and K.I. contributed materials/analysis tools; S.P. and C.B. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank S. Abate of CNR-SPIN Salerno (Italy) for his technical support. The nanowires patterning has been performed at NanoFacility Piemonte, INRiM, a laboratory supported by Compagnia di San Paolo. S.P., N.M. and C.B. acknowledge partial support from University of Salerno through grants FARB17PAGAN and FARB18CAVAL, and from INFN through experiment FEEL.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Nanomaterials* Editorial Office E-mail: nanomaterials@mdpi.com www.mdpi.com/journal/nanomaterials

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-4762-6