Introduction

Quantum phase transitions (QPTs) corresponding to superconductor-insulator transition (SIT) and superconductor-metal transition (SMT) in two-dimensional (2D) superconductors have been a subject of great interest in last few decades1,2,3. QPTs can be initiated by tuning non-thermal external control parameters such as disorder (through film thickness variations)4, magnetic field5, carrier concentrations6, pressure among others7. There exists a quantum critical point (QCP) associated with a continuous QPT which separates the superconducting ground state to the normal insulating/metallic state of the system. In the immediate vicinity of a continuous QPT, a d-dimensional quantum system can be represented by a (d + 1) dimensional classical system with imaginary time as the extra dimension. Upon approaching the QCP, both correlation length and correlation time diverge. The divergence can be expressed using characteristic length scales in both spatial (ξ) and temporal directions (\({\xi }_{\tau }\)) as, \(\xi \sim {|\delta |}^{-\nu },{\xi }_{\tau } \sim {\xi }^{z} \sim {|\delta |}^{-z\nu }\), where \(\delta\) measures the distance from the critical point. The exponents \(\nu\) and \(z\) correspond to the correlation length critical exponent and the dynamical critical exponent, respectively. Thus, a continuous QPT leads to a power-law singularities in the length scales at the QCP indicating a power-law singularities in observables too8 and the corresponding scaling of physical quantities represent universal behavior in the critical exponents \({z}\nu\) 7.

However, in contrast to the QPT associated with a single QCP in the clean limit, the emergence of multiple critical points (MCPs) in magnetic-field driven QPT has been observed in various systems9,10 where the disorder plays a major role in controlling the phase transition. The presence of spatial quenched disorders in the form of impurities, defects, dislocations, grain boundaries and others are inevitable in any real sample and assessing their influence on QPT in d-dimensional system is generally done by the Harris criterion, dν > 2 with \(\nu\) being the correlation length exponent11. On large length scale, when the effective disorder strength and its randomness average out under course graining, the Harris criterion is fulfilled and the system undergoes a pure fixed critical point transition with clean limit critical exponents. In contrast, the quenched disorder becomes relevant for the systems that violate the Harris criterion12. When the disorder strength reaches a finite value upon course graining, the randomness becomes competitive at all scale and the transition is controlled by a random fixed critical point with critical exponents (\({z}^{{\prime} }{\nu }^{{\prime} }\)) differing from those of the clean limit (zν). Here, the modified correlation length exponent (\({\nu }^{{\prime} }\)) fulfils the Harris criterion and the transition resembles the conventional clean critical behavior and the dynamical critical exponent (\({z}^{{\prime} }\)) saturates at a disorder dependent finite value12,13. In this case, the physical properties of the system are dominated by the rare regions that are spatially localized by the influence of the quenched disorder. Further, an extreme case appears when the disorder strength diverges at the transition under course graining, and the Harris criterion is violated. The transition here is governed by a fixed critical point which is of infinitely strong randomness with diverging dynamical critical exponent. Here, the characteristic length scales in the temporal (\({\xi }_{\tau }\)) and spatial direction (\(\xi\)), in contrast to the clean limit power-law scaling (\({\xi }_{\tau } \sim {\xi }^{z}\)), follow an activated scaling as \({\xi }_{\tau } \sim {e}^{{\xi }^{\psi }}\) with \(\psi\) being the tunneling correlation exponent. This extreme case of quenched disorder mediated criticality of infinite randomness with diverging dynamical critical exponent at the transition is known as quantum Griffiths singularity (QGS)14. Therefore, the systematic ways to establish the existence of QGS include the appearance of a continuum of MCPs in field dependent resistance Rs(B) at low temperature, divergence of dynamical exponent at zero-temperature limit while approaching the QPT and the applicability of an activated scaling to describe the physical quantities near the transition.

In 2D superconductors, the first evidence of QGS came from molecular-beam-epitaxy (MBE) grown Ga trilayer15 which was then followed by other low dimensional superconducting systems16,17,18,19,20,21,22,23,24,25,26. The observation of MCPs near the field induced QPTs in these systems led to the emergence of QGS at very low temperature above the mean field upper critical field \({B}_{{c}_{2}}^{{MF}}\)15,25,27. The QGS can be originated by quenched disorder mediated formation of rare ordered superconducting puddles, known as Griffith’s regions that are coupled through long range Josephson coupling and form a vortex glass like phase where the critical field of the rare ordered Griffith’s region becomes larger than the upper critical field corresponding to the clean limit QPT. Though the emergence of QGS has mainly been correlated with systems exhibiting SMT due to stronger Josephson coupling with the metallic normal state15,20,28, of late, its association with SIT based systems has been experimentally demonstrated for a few materials18,20. In order to explore the QGS in both type of systems, a wide range of variations in normal state resistance (RN) is needed among the samples. For example, about 20 fold variation in RN was reported for TiO thin films20, whereas, 3-4 times variation in resistance was considered for the NbN samples18. Here, we demonstrate the observation of QGS in 2D disordered TiN thinfilms with a variation of more than two orders of magnitude in the RN.

Disordered TiN is one of the initial materials which has been able to successfully demonstrate the experimental observation of quantum criticality2, field- and disorder-induced SIT2,29,30, multiple crossing points10 among the other interesting quantum phenomena2,31,32,33. But, till now, no report has shown the existence of QGS in TiN thin films. Here, we report the observation of multiple crossing points in the magnetoresistance isotherms along with a diverging dynamical critical exponent at low temperature while approaching the QPT in disordered TiN thin films of thickness ranging 2–4 nm. A direct activated dynamical scaling analysis on the field and temperature dependent sheet resistance Rs(B,T) for a wide range of temperature further confirms the infinite-randomness fixed critical point to be governing the QPT. Based on the detailed analysis, the emergence of the QGS in TiN thin films is established and the corresponding field-temperature phase diagram features a broad QGS region extending up to a temperature which is very close to the Tc of the sample.

Results

The temperature dependent sheet resistance Rs(T) measurements for five representative TiN samples of thickness 2–4 nm are presented in Fig. 1. Three samples measured down to 2 K are shown in Fig. 1a and the rest two measured down to 300 mK take place in Fig. 1b. With reduction in thickness, the normal state resistance (RN) increases and the transition temperature (Tc) decreases. Here, Tc corresponds to the temperature at which the derivative dRs/dT becomes maximum. The Tc values closely match with their respective mean field transition temperature Tc0 as obtained from the quantum corrections to the conductivity (QCC) theories31,34. Further, while cooling down from room temperature, an upturn with a resistance maximum Rmax appears for all the samples just before the superconducting onset indicating a weakly localized metallic state34. The appearance of the resistance peak Rmax at the temperature Tmax is shown by the arrows in Fig. 1a for the sample TN10A and in Supplementary Fig. 1 for samples TN9 & TN8. It is evident that with reduction in thickness, samples become more disordered and more localized as indicated by the increased RN and the steeper slope for the negative dRs/dT region of the resistance upturn, respectively.

Fig. 1: Temperature dependent sheet resistance Rs(T) measured on TiN thin films with variation in film thickness.
figure 1

Based on the lowest achievable measurement temperature, two sets of samples are presented in (a) and in the main panel of (b). a Rs(T) for the temperature range from room temperature down to 2 K for three representative samples from the first set. Inset: the same set of Rs(T) but for a selected range of temperature to highlight the metal-superconductor transition region. b Rs(T) near the transition region for other two samples from the second set measured down to 300 mK. Here, the samples TN10A & TN10B were prepared in one single run and the samples TN9 & TN9A were grown in the same batch. Inset: The resistance Rmax, measured at Tmax as shown by the arrow in Fig. 1a just above the superconducting onset, is plotted against the resistance measured at room temperature (R300 K) for all the samples. The dashed horizontal line, separating the blue and purple shaded regions, refers the quantum resistance RQ. The resistance values for the samples are compared with respect to the quantum resistance.

Full size image

In the inset of Fig. 1b, the RN values for all the samples presented in this article and in Supplementary Fig. 2 (Supplementary Note 1) are summarized with respect to the quantum resistance RQ which is denoted by the horizontal dashed line separating the two shaded regions in cyan and purple. As for the reference of RN, the resistance Rmax is plotted against the room temperature resistance R300 K. Here, the selected samples offer more than two orders of change in the RN from ~140 Ω to ~18 kΩ. The wide variations in RN are obtained by tuning the thickness and the annealing temperature (Ta) during the growth process. Two different temperatures for Ta, 780 °C and 820 °C, are considered here, and the resistance values for the corresponding samples are shown in olive-green and blue rectangles, respectively. Comparing the resistance with respect to RQ, the samples with resistance less than 500 Ω are considered here as weakly disordered (RS≪RQ), the samples with resistance more than 1.5 kΩ are considered as moderately disordered (RS < RQ) and the samples with resistance much higher than RQ are considered as strongly disordered (RS≫RQ).

Emergence of multiple crossing points in magnetoresistance isotherms

The magnetic-field and temperature dependent sheet resistance Rs(B,T) measurements for the samples measured down to 2 K are shown in Fig. 2 where the left panel represents the Rs(T) measured under different field and the right panel displays the magnetoresistance isotherms Rs(B). The sample thickness varies in descending order from the top to the bottom. For comparison among the samples, the variations of Rs(T) with field are placed in equivalent scaling as evident by their normalized resistance Rs/Rmax. With increasing field, a resistance plateau indicating a field induced SMT (or weakly localized metallic transition)18 is observed for the 4 nm (TN8) and 3 nm (TN9) thick samples, whereas, a relatively strong upturn indicating a SIT or superconductor to weak insulator transition20 appears for the thinnest sample TN10A (2 nm). Further, a selective region of the measured Rs(B) isotherms from each of the samples TN8, TN9 & TN10A is highlighted in Fig. 2b, d, f, respectively. In contrast to the conventional single point crossing, the Rs(B) isotherms, measured between 2 and 5 K with 0.1 K interval, display a series of crossing points (shown in black filled diamonds) that form a continuous line of SMT/SIT critical points. The appearance of a continuum of MCPs hints towards the existence of QGS at very low temperature (\(T\to 0\)) where the critical exponents product \(z\nu\) diverges15,16,21. Furthermore, a quantitative analysis on the relative span ΔBc/ΔT covered by the MCPs is done with respect to the sample thickness. Here, ΔBc is the extent in magnetic field for a specific range of temperature ΔT. The values of ΔBc/ΔT for samples TN8 (4 nm), TN9 (3 nm) and TN10A (2 nm) are found as 0.306 T/K, 0.446 T/K and 0.896 T/K, respectively. Thus, thickness reduction leads to an enhancement in the relative span of the MCPs and from the trend, one can expect a single crossing point in relatively thicker samples. A similar trend is observed in Supplementary Fig. 3 (Supplementary Note 2) for the samples TN4 (4 nm) and TN5A (3 nm) that are annealed at 820 °C.

Fig. 2: Magnetic field dependent sheet resistance for the first set of samples (TN8, TN9A & TN10A) measured down to 2 K.
figure 2

Field dependent Rs(T) for (a) TN8, (c) TN9A & (e) TN10A and isothermal magnetoresistance [Rs(B)] demonstrating multiple crossing points in the temperature range from 2 K to 5 K for (b) TN8, (d) TN9A & (f) TN10A. Here, the black solid diamonds mark the crossing points (Bc) for three consecutive Rs(B) isotherms measured with an interval of 0.1 K in the temperature. Double sided arrows show the span ΔBc in the field covering the multiple crossing region for the temperature window of 2–5 K. Here, it should be noted that a resistance change of about ~ 3 Ω between Rs(T) and Rs(B) for the sample TN8 occurs due to thermal cycling during the measurements.

Full size image

Figure 2 indicates that the field extent covered by MCPs increases with increasing disorder and hence the randomness in the critical behavior gets stronger with increasing disorder. In Fig. 3, we present Rs(B,T) measurements for a highly disordered sample (TN6) with RN≫RQ. The zero-field Rs(T) measured from room temperature down to 2 K as shown in Fig. 3a clearly shows a non-metallic behavior in the normal state above the Tmax with negative dRS/dT in the whole temperature range from 300 K down to Tmax. Here, the upward slope gets much steeper just above the Tmax while approaching the Tmax from higher temperature side and a resistance peak Rmax appears at the Tmax. This leads to a much higher Rmax than the resistance R300 K measured at 300 K. Further, for T < Tmax, a resistance drop of about 25% of the Rmax is observed down to 2 K. The evolution of the transition region of the Rs(T) under magnetic field as shown in the inset of Fig. 3a presents a transition from superconducting to insulating regime upon the field variation from 2 T to 4 T. Further, Fig. 3b presents Rs(B) isotherms for the temperature window of Tmax to 2 K where a continuum of MCPs (shown by the black filled diamonds) is clearly evident. The relative span (ΔBc/ΔT) of the MCPs appears as ~2.37 T/K which is indeed much higher than the samples shown in Fig. 2 with relatively low RN. Therefore, MCPs are observed in ultrathin TiN (≤4 nm) samples spanning in wide range of disorder which indicates the existence of MCPs in both types of field induced QPTs SMT and SIT.

Fig. 3: Field dependent Rs(T) and Rs(B) isotherms for a highly resistive TiN thin film sample (TN6) with Rs≫RQ in its normal state.
figure 3

The magnetoresistance Rs(B) data is extracted from the field dependent Rs(T) by the inverse matrix method. The thickness of the sample is (2 ± 0.5) nm and the annealing temperature (Ta) during its growth was ~ 820 °C. Here, the increased resistance is caused by the formation of elemental Si at Ta ~ 820 °C in addition to TiN. a Zero-field Rs(T) measured from room temperature down to 2 K. Inset: Rs(T) measured under external magnetic field applied perpendicular to the sample plane. b Isothermal Rs(B) obtained for the temperature window of 2–3.2 K. The black filled diamonds represent the crossing points between the neighboring isotherms in small temperature intervals.

Full size image

Emergence of quantum Griffiths singularity

The continuum of MCPs observed in Rs(B,T) measured down to 2 K (shown in Fig. 2, 3) indicates towards the existence of QGS in these samples at zero temperature limit19 where the observation of diverging critical exponent \(z\nu\) provides the solid evidence of QGS28. Consequently, two samples TN9A (~3 nm) and TN10B (~2 nm) exhibiting field-induced SMT & SIT, respectively, are measured down to 300 mK. The Rs(B,T) measurements for TN9A are presented in Fig. 4. With increasing magnetic field from 3.5 T to 3.75 T in the field dependent Rs(T) as shown in Fig. 4a, transition from superconducting to weakly localized metallic state is observed at the lowest measurement temperature. Further, Fig. 4b presents the Rs(B) isotherms measured in the temperature interval 0.3 K ≤ T ≤ 3 K where the neighboring Rs(B) isotherms cross at different critical fields (shown in black diamonds) that form a continuous line of MCPs indicating the existence of QGS19,28. This continuum of MCPs spans over a wide range of temperature and magnetic field. The temperature dependence of the crossing field Bc(T) is shown in the inset of Fig. 4b. Here, the average temperature is considered at each crossing point and the error bars represent the range.

Fig. 4: Probing quantum Griffiths singularity by magnetoresistance measurements carried out at lower temperature down to 300 mK for the sample TN9A ( ~3 nm).
figure 4

a A selected region of field dependent Rs(T) measured under perpendicular magnetic fields from 0 to 4.5 T in 0.25 T steps. b Isothermal Rs(B) measured in the temperature window of 0.3 K ≤ T ≤ 3 K. The magnetoresistance isotherms are measured with a temperature interval of 0.05 K between two neighboring isotherms from 0.3 K to 0.85 K, with 0.1 K interval from 0.85 K to 1.85 K and with 0.25 K interval from 2 K to 3 K. The crossing points are shown by the black diamonds. Inset: The crossing field (Bc)—temperature (T) phase boundary near the quantum superconductor-to metal transition. The red curve is the fit using Eq. (5) with \({B}_{c}^{* }\), \({u}\) and \(p\) as adjustable parameters and the characteristic temperature T0 fixed at 4.41 K (see the text for details). c Divergent behavior of critical exponent \(z\nu\) as a function of magnetic field. The \(z\nu\) values are obtained through finite size scaling (FSS) analysis for a set of adjacent magnetoresistance isotherms. The solid red curve is the fit based on activated scaling law using Eq. (2) and two black dashed lines represent \(z\nu =1\) (horizontal) & \({B}_{c}^{* }=3.78{T}\) (vertical). d Sheet resistance as a function of the scaling parameter \(\delta {\left({{{{\mathrm{ln}}}}}\left({T}_{0}/T\right)\right)}^{1/\nu \psi }\) related to the activated dynamical scaling as described in Eq. (4) for the temperature span 0.3 K ≤ T ≤ 1.75 K. Error bars denote the temperature range for a set of magnetoresistance isotherms.

Full size image

The existence of QGS can be established by the divergence of critical exponent \(z\nu\) upon approaching the QPT at T → 0 and B → \({B}_{c}^{* }\), where \({B}_{c}^{* }\) is the characteristic critical field15,16,21. Here, the critical exponent \(z\nu\) is obtained by adopting the finite size power law scaling (FSS) analysis as is generally used for a continuous QPT associated with a fixed single QCP7. Under the power law scaling, the temperature and magnetic field dependent sheet resistance Rs(B,T) near a field induced QPT can be described as7,

$${R}_{s}(B,T)={R}_{s}^{c} \, f(\delta {T}^{-1/z\nu })$$
(1)

where, \(\delta =|B-{B}_{c}|\) is the distance from the critical field \({B}_{c}\), \({R}_{{{{{{\rm{s}}}}}}}^{{{{{{\rm{c}}}}}}}\) is the sheet resistance at the critical point and \({f}(x)\) is the scaling factor with \(f\left(0\right)=1\). \(z\) and \(\nu\) are the dynamical and corelation length critical exponents, respectively. To perform the FSS analysis on the Rs(B,T) data shown in Fig. 4b, a set of three adjacent Rs(B) isotherms crossing at a particular critical field Bc in each small temperature interval is considered for extracting \(z\nu\) values by using the power law scaling in Eq. (1). The detailed FSS analysis is illustrated in Supplementary Figs. (4)–(6) in the Supplementary Note 3. The extracted \(z\nu\) values are plotted against the critical fields in Fig. 4c. Initially, zν increases slowly with magnetic field from 3 K down to ~ 1.55 K, with further lowering temperature, \(z\nu\) increases rapidly with field and shows a diverging behavior while approaching the characteristic critical field \({B}_{c}^{* }\). This reflects an activated scaling behavior of \(z\nu\) with respect to the field near the characteristic critical field \({B}_{c}^{* }\) and is expressed as35,

$${z\nu \, \approx \, C(B-{B}_{c}^{* })}^{-\upsilon \psi },$$
(2)

where C is a constant and \(\upsilon =1.2\) & ψ = 0.5 are the corelation length exponent and the tunneling exponent associated with the 2D infinite-randomness fixed critical point (IRFCP), respectively35,36,37. The activated scaling law is in good agreement with the experimental data as shown by the red solid curve in Fig. 4c which yields \({B}_{c}^{* }\) = 3.78 T. This indicates towards the existence of IRFCP justifying the presence of QGS in the sample.

Further, a QPT governed by an IRFCP can be described by a direct activated dynamical scaling as a whole where the infinite-randomness critical exponent \(\upsilon \psi\) can replace the critical exponent \(z\nu\) associated with the conventional power law scaling. Here, at \(T\to 0\) the critical exponent products from power law scaling (zν) and activated dynamical scaling (\(\upsilon \psi\)) are related as28,

$$\left(\frac{1}{z\nu }\right)=\left(\frac{1}{\upsilon \psi }\right)\cdot \frac{1}{ln({T}_{0}/T)}$$
(3)

where the exponent \(\upsilon \psi\) (\(\upsilon\) is the corerelation length critical exponent and \(\psi\) is the tunneling exponent) represents the universality class of the QPT associated with IRFCP. The fit using Eq. (3) to the temperature dependent \(z\nu\) (shown in Supplementary Fig. 7 in the Supplementary Note 4) leads to the characteristic temperature T0 and the critical exponent \(\upsilon \psi\) which are used as free parameters. Further, Eq. (3) indiates a diverging \(z\nu\) while approaching to \(T\to 0\), confirming the QCP is of IRFCP type where physical quantities can be described by the direct activated dynamical scaling. Consequently, the sheet resistance Rs can be expressed as28,

$${R}_{s}\left(\delta ,{{{{\mathrm{ln}}}}}\frac{{T}_{0}}{T}\right)=\Phi \left[\delta {\left({{{{\mathrm{ln}}}}}\frac{{T}_{0}}{T}\right)}^{1/\upsilon \psi }\right]$$
(4)

where, \(\delta =\left|{B}_{c}-{B}_{c}^{* }\right|/{B}_{c}^{* }\) is the distance from the critical field \({B}_{c}^{* }\) associated to the IRFCP and \(\Phi\) corresponds to the scaling function. Here, T0 and \(\upsilon \psi\) are the same as have been used in Eq. (3). Further, the crossing field Bc(T) near the transition is also expected to follow the activated dynamical scaling with irrelevant scaling correction as28,

$${B}_{c}\left(T\right)={B}_{c}^{* }\times \left[1-u{\left({{{{\mathrm{ln}}}}}\frac{{T}_{0}}{T}\right)}^{-p}\right]$$
(5)

with the exponent \(p=\frac{1}{\upsilon \psi }+\frac{\omega }{\psi }\) and u is the leading irrelavant scaling variable responsible for the correction to the scaling. \(\omega\) is the associated irrelevant exponent which is always positive22. The red solid curve in the inset of Fig. 4b represents the fit to the Bc(T) phase boundary using Eq. (5) where \({B}_{c}^{* }\), \(u\) and \(p\) are considered as the fitting parameters. T0 = 4.41 K and \(\upsilon \psi =0.63\) are obtained from the fit using Eq. (3) [Supplementary Fig. 7a in the Supplementary Note 4]. Here, \(\upsilon \psi =0.63\) closely matches with the theoretically predicted value of ~ 0.6 with \(\upsilon \sim 1.2\) & ψ ∼ 0.535,36,37. The best fit values for \({B}_{c}^{* }\), \(u\) and \(p\) are obtained as (3.77 ± 0.02) T, (0.115 ± 0.007) and (1.92 ± 0.15), respectively.

Further, in Fig. 4d, the direct activated dynamical scaling using Eq. (4) has been applied on a set of Rs(B) isotherms and the best collapse is obtained for \({B}_{c}^{* }=3.77{{{{{\rm{T}}}}}}\), T0 = 4.56 ± 0.23 K and \(\upsilon \psi =0.63\) in the temperature interval 0.3 K ≤ T ≤ 1.75 K. Here, \({B}_{c}^{* }\) is obtained by minimizing the variance in magnetoresistance isotherms when plotted against the scaling parameter18,28. The same value for \({B}_{c}^{* }\) is obtained from the phase boundary fitting of Bc(T) using Eq. (5) [inset of Fig. 4b]. Furthermore, \({B}_{c}^{* }=3.77{{{{{\rm{T}}}}}}\), agrees closely to 3.78 T which is obtained from the fit to the field dependence of \(z\nu\) using Eq. (2) [shown in Fig. 4b]. The characteristic temperature T0, used as an adjusting parameter to achieve the best collapse, is obtained as 4.56 ± 0.23 K which is close to the temperature 4.41 K as found from the fit of \(z\nu\) versus temperature using Eq. (3) (Supplementary Fig. 7 in the Supplementary Note 4). It is observed in Fig. 4d that for large values of the scaling parameter away from \({B}_{c}^{* }\), the scaling collapse starts to break down at ~ 4.29 T and ~1.85 T, for the upper and lower branches, respectively. These breakdowns can be attributted to the negligible quantum fluctuations occuring in the normal state at high magnetic field (upper branch) and in the superconducting state at very low magnetic field and low temperature (lower branch)18,28. It should be noted that the temperature window of 0.3 K ≤ T ≤ 1.75 K remains the same for (i) the direct activated scaling of the MR isotherms using Eq. (4) [in Fig. 4d], (ii) the Bc(T) phase boundary fitting using Eq. (5) [in the inset of Fig. 4(b)] and (iii) the \(z\nu\) vs. temperature fit using Eq. (3) [in Supplementary Fig. 7 in the Supplementary Note 4]. However, at high temperature, the quantum fluctuation becomes insignificant compared to the thermal fluctuations and the data set starts to deviate from the activated scaling.

Next in Fig. 5, we present the Rs(B,T) measurements carried out at temperature down to 300 mK for a more resistive (RN ~ 2 kΩ) sample TN10B (~2 nm) [shown in Fig. 1b] with weak insulating background where the formation of rare ordered superconducting puddles (Griffiths regions) via long-range Josephson coupling becomes more challenging compared to a metallic background. A selected region from the set of Rs(T), measured under fixed perpendicular fields from 0 to 5 T in 0.25 T steps, is shown in Fig. 5a. The trace of superconductivity with positive downward slope is observed for field up to 3.75 T. But, a clear upward transition with negative slope (dRs/dT < 0) occurs for an increased field of 4 T. Here, as the RN is lower than the RQ, the logarithmic derivative of the sheet conductance \(\left({\sigma }_{s}=1/{R}_{s}\right)\) is considered for examining the state for B ≥ 4 T. Though dRs/dT is negative but often this does not ascertain whether a sample is metallic or insulating18,38. It is shown that the logarithmic derivative of sheet conductance \(w={{{{{\rm{dln}}}}}}{\sigma }_{s}/{{{{{\rm{dln}}}}}}T\) at zero-temperature limit (\(T\to 0\)) is more sensitive than the derivative dRs/dT in determining whether a sample is metallic or insulating. For \({w} \, > \, 0\), a positive value of dw/dT indicates a metallic state whereas, a negative value of dw/dT indicates the sample to be very likely insulating39. In the inset of Fig. 5a, \(w={{{{{\rm{dln}}}}}}{\sigma }_{s}/{{{{{\rm{dln}}}}}}T\) corresponding to the Rs(T) measured under 4 T is plotted against T1/2. Here, \(w\) is positive and towards the lowest measurement temperature \(w\) increases with decreasing temperature leading to a negative value of\(\,{{{{{\rm{d}}}}}}w/{{{{{\rm{d}}}}}}T\). This indicates a field induced superconducting to insulating quantum phase transition (SIT) at zero-temperature limit for B ≥ 4 T. The insulating phase near the transition has been investigated further by Rs(B) isotherms measured at low temperature down to 300 mK and for a magnetic field up to 10 T. The corresponding magnetoresistance (MR) data is shown in Supplementary Fig. 10b in the Supplementary Note 6 which unveils a broad resistance peak for B ≥ 4 T. The appearance of resistance peak near the transition indicates the field-induced localization of Cooper pairs2,29,40,41. However, as the RN is much lower than the quantum resistance RQ of Cooper pairs, a parallel contribution from fermionic channel of quasiparticles is most likely to be present in addition to the bosonic contribution from localized Cooper pairs. Further, the logarithmic temperature dependence of the conductance [Supplementary Fig. 10c, d in the Supplementary Note 6] at higher field confirms the fermionic scenario5,34,42,43.

Fig. 5: Temperature and magnetic field dependent resistance Rs(B,T) measurements carried out at low temperature down to 300 mK for the sample TN10B (thickness 2 ± 0.5 nm).
figure 5

a A selected region of Rs(T) measured under perpendicular magnetic fields from 0 to 5 T in 0.25 T steps. Inset: Logarithmic derivative of sheet conductance (dlnσs/dlnT) vs. T1/2 obtained from the Rs(T) measured under 4 T field. b A detailed view of isothermal Rs(B) measured in the temperature window of 0.3 K ≤ T ≤ 3 K. The black diamonds represent the crossing points. Inset: The crossing field (Bc) – temperature (T) phase boundary near the quantum superconductor-to metal (SMT) transition. The red curve is the fit using Eq. (5) with \({{{{{{\rm{B}}}}}}}_{{{{{{\rm{c}}}}}}}^{* }\), \({{{{{\rm{u}}}}}}\) and \({{{{{\rm{p}}}}}}\) as adjustable parameters while the characteristic temperature T0 was fixed at 3.11 K (see the text for details). c Field-dependence of the critical exponent product zν which shows diverging behavior with magnetic field while approaching the characteristic critical field \({{{{{{\rm{B}}}}}}}_{{{{{{\rm{c}}}}}}}^{* }\) at the lowest measurement temperature. The zν values are obtained through finite size scaling (FSS) analysis for a set of adjacent magnetoresistance isotherms in each small interval of temperature. The solid red curve is the fit based on the activated scaling law using Eq. (2). The black dashed lines represent \({{{{{\rm{z}}}}}}{{{{{\rm{\nu }}}}}}=1\) (horizontal) & \({{{{{{\rm{B}}}}}}}_{{{{{{\rm{c}}}}}}}^{* }=4.03{{{{{\rm{T}}}}}}\) (vertical). The horizontal dashed-dotted broken line indicates a plateau region where \({{{{{\rm{z}}}}}}{{{{{\rm{\nu }}}}}}\) remains almost unchanged with changing magnetic field in the high temperature regime. d Sheet resistance as a function of the scaling parameter \({{{{{\rm{\delta }}}}}}{\left({{{{\mathrm{ln}}}}}\left({{{{{{\rm{T}}}}}}}_{0}/{{{{{\rm{T}}}}}}\right)\right)}^{1/{{{{{\rm{\nu }}}}}}{{{{{\rm{\psi }}}}}}}\) related to the activated dynamical scaling as described in Eq. (4) for the temperature span of 0.3 K ≤ T ≤ 1.75 K. Error bars denote the temperature range for a set of magnetoresistance isotherms.

Full size image

Further, a detailed view near the field induced transition to the insulating state is shown in Fig. 5b by the Rs(B) isotherms measured in the temperature interval 0.3 K ≤ T ≤ 3.0 K. A continuum of MCPs (shown by the black diamonds) is evident. The phase boundary as represented by the temperature dependent crossing field Bc(T) is displyed in the inset of Fig. 5b. The red solid curve represents the activated dynamical fit using Eq. (5) for the temperature range 0.3 K ≤ T ≤ 1.75 K with T0 = 3.11 K and \(\upsilon \psi =0.62\) as obtained from the fit using Eq. (3) to the temperature dependent \(z\nu\) (Supplementary Fig. 7b in the Supplementary Note 4). The critical field \({B}_{c}^{* }\) is obtained as 4.03 ± 0.02 T. The best fit values for the parameters \(u\), \(p\) and the critical field \({B}_{c}^{* }\) are obtained as 0.031 ± 0.006, 2.15 ± 0.26 and 4.03 ± 0.02 T, respectively.

The field dependent \(z\nu\) is plotted in Fig. 5(c) which indicates a diverging \(z\nu\) while approaching the characteristic critical field \({{B}}_{c}^{* }\) at the low temperature limit. The experimental data is in good agreement with the activated scaling law as described by the red solid curve using Eq. (2) indicating towards the existence of IRFCP transition which justifies the presence of QGS in the systems. The best fit value for the characteristic field \({B}_{c}^{* }\) is obtained as 4.03 T. Further in Fig. 5d, a direct activated dynamical scaling using Eq. (4) is performed on the Rs(B,T). Here, the best collapse is observed for \({{B}}_{c}^{* }=4.0{{{{{\rm{T}}}}}}\) in the temperature range 0.3 K ≤ T ≤ 1.75 K with \(\upsilon \psi\) fixed at the value of 0.62. The best collapse, obtained by using the characteristic temperature T0 as an adjusting parameter, leads to T0 = 3.27 ± 0.25 K. The critical field \({{B}}_{c}^{* }=4.0{{{{{\rm{T}}}}}}\) is very close to the value of 4.03 T which is obtained from the phase boundary fitting of Bc(T) [inset of Fig. 5b] as well as from the activated scaling law as presented in Fig. 5c. Further, the characteristic temperature T0 is in good agreement too with the value obtained from the temperature dependence of \(z\nu\) shown in Supplementary Fig. 7b in the Supplementary Note 4. Here, the scaling collapse breaks down at high temperature and at fields far from the critical field \({{B}}_{c}^{* }\) due to the diminishing quantum fluctuations28.

Here, it is noteworthy to mention that Eq. (4) associated with the direct activated dynamical scaling of the MR isotherms predicts a single critical point \({B}_{c}^{* }\) which does not account for the temperature dependence of the crossing fields as observed by the MCPs. The scaling behaviors, shown in Fig. 4d and in Fig. 5d, confirm the existence of a single critical point in the magnetic field for a wide range of temperature from 300 mK to about 1.75 K for both the samples TN9A and TN10B. Further, Eq. (3) indicates a diverging behavior of the critical exponent \(z\nu\) with temperature at \(T\to 0\) and the corresponding fit shown in Supplementary Fig. 7 establishes the diverging behavior of \({z}\nu\) at low temperature. These observations confirm the QPT to be governed by an IRFCP. However, there is a clear difference in the temperature dependent critical field Bc(T) for these two samples TN9A (SMT category) and TN10B (SIT category). For the former, the critical field moves towards upward direction while approaching the zero-temperature limit [inset of Fig. 4b], wheras, for the latter, Bc increases initially with decreasing temperature in the high temperature regime (\(T \, \gtrsim \, 1\,{{{{{\rm{K}}}}}}\)) but finally it tends to saturate around \({B}_{c}^{* }\,(4.03\,{{{{{\rm{T}}}}}})\) as the temperature approaches to 0 K [inset of Fig. 5b]. The upward trend of \({B}_{c}(T)\) as observed for TN9A is similar to the other reported SMT systems exihibiting QGS15,18,20,24,25,26, whereas, a saturating Bc at zero-temperature limit is one of the features observed in systems showing SIT20. The differences in the behavior of Bc at \(T\to 0\) might be related to the different Josephson interaction strength among the superconducting puddles with insulating & metallic normal states, for the SIT and SMT systems, respectively.

Discussion

Based on the low temperature magnetotransport measurements, a comprehensive B-T Phase diagram has been constructed for the 2-nm thick sample TN10B for which the normal state is weakly insulating and is considered in the category of SIT system. The phase diagram is shown in Fig. 6 where TcOnset(B) (the half-filled green diamonds) and the crossing field Bc(T) (the pink squares) follow the same track. Here, the TcOnset(B) separates the weakly localized insulating regime with negative dRs/dT from the superconducting fluctuation regime with positive dRs/dT which is mainly driven by the thermal fluctuation at the high temperature limit for T > Tm = 1.75 K. The separation of the weakly localized insulating normal state from the thermal fluctuation regime is shown by the blue broken curve, connecting the TcOnset(B) and the Bc(T) above 1.75 K. The mean field critical field \({B}_{c}^{{MF}}(T)\) (shown by the blue spheres) is obtained from the field dependent Rs(T) by using the Ullah-Dorsay (UD) scaling method21,44,45,46 by which in actuality, the mean field critical temperature \({T}_{c}^{{MF}}(B)\) is extracted for a particular field. The UD scaling is done for the field range 1.5 T ≤ B ≤ 3.25 T and outside of this window the scaling deviates (The details are shown in Supplementary Fig. 9 in the Supplementary Note 5). It is found that the obtained \({T}_{\!\!{c}}^{{MF}}\) closely matches with the Tc. Accordingly, for field lower than 1.5 T, \({T}_{\!\!{c}}^{{MF}}(B)\) is represented by the Tc(B). A fit using the empirical formula16,18 \({B}_{c}^{{MF}}\left(T\right)={B}_{c}^{{MF}}\left(0\right)[1-{\left(T/{T}_{c}\right)}^{2}]\) appears to be in good agreement with the mean field critical field \({B}_{c}^{{MF}}(T)\) as represented by the solid back curve. From the fit, the upper critical field \({B}_{c}^{{MF}}(0)\) and the Tc(0) are obtained as 3.65 T and 1.9 K, respectively. At temperature lower than the mean field critical temperature \({T}_{c}^{{MF}}(B)\), phase fluctuation is the dominant mechanism which relates to the thermally activated flux flow (TAFF) at high temperature and to quantum creep at low temperature21. These two regions are separated by the characteristic temperature TTAFF as represented by the purple triangles. TTAFF is extracted from the Arrhenius plot analysis of the logarithmic of Rs(T) measured under a particular field when plotted against inverse temperature (T−1). The TTAFF represents the crossover temperature between the TAFF regime (pink region) and the quantum phase fluctuation regime (blue region). The quantum phase fluctuation regime below the TTAFF (B) is marked as the condensed superconducting (SC) phase.

Fig. 6: B-T Phase diagram constructed from the Rs(B,T) measurements for the sample TN10B ( ~ 2 nm).
figure 6

The half-filled green diamonds represent the temperature TcOnset corresponding to the temperature at which the slope dRs/dT changes its sign from negative to positive in the field dependent Rs(T) while approaching the superconducting transition. The pink half-filled squares are the crossing critical fields Bc(T) obtained from the crossing points of the neighboring MR isotherms in each small temperature interval. The blue spheres represent the mean-field critical field \({B}_{c}^{{MF}}\left(T\right)\) obtained from the field dependent Rs(T) by using the Ullah-Dorsay scaling method for the field range 1.5 T ≤ B ≤ 3.25 T. For fields lower than 1.5 T, the Tc values are considered at the place of TcMF. The purple triangles represent the characteristic temperature TTAFF which separates thermally activated flux flow (TAFF) regime from the quantum phase fluctuation regime (denoted here as SC regime) and is obtained from the Arrhenius plot analysis of the field dependent Rs(T). The black solid curve is a fit for the \({B}_{c}^{{MF}}\left(T\right)\) using the empirical formula \({B}_{c}^{{MF}}\left(T\right)={B}_{c}^{{MF}}\left(0\right)\left[1-{\left(T/{T}_{c}\right)}^{2}\right]\). The green solid curve is the fit to the Bc-T phase boundary at the high field regime using activated dynamical scaling, as expressed in Eq. (5). Here the SC region, as separated from the thermal creep region by TTAFF (purple triangles), is extended up to the characteristic critical field \({B}_{c}^{* }\) in such a way that the extended region above the \({B}_{c}^{{MF}}\) corresponds to zν ≥ 1. The blue broken curve following the TcOnset (Green diamonds) is the guide to the eye which separates the weakly localized insulating state from the thermal fluctuation regime at high temperature (T > 1.75 K). At T < 1.75 K, the green solid curve serves as the boundary between the weakly localized insulating region and the QGS region.

Full size image

Finally, the QGS region is marked based on the activated dynamical scaling analysis associated with the IRFCP as presented in Fig. 5. Here, the magnetoresistance data can be adequately described by the activated dynamical scaling [Eq. (4)] for the temperature span of 0.3 K ≤ T ≤ 1.75 K. In the same temperature window, the Bc(T) can be portrayed by the IRFCP governed activated dynamical scaling with the irrelevant scaling corrections [Eq. (5)] as shown by the green solid curve in Fig. 6. Thus, the activated dynamical scaling analysis on the measured data suggests for a wide range of temperature up to 1.75 K above the \({B}_{c}^{{MF}}(0)\) as the region of QGS. This is consistent with the field dependence of \(z\nu\) presented in Fig. 5c which displays a plateau region where \(z\nu\) remains almost constant with magnitude < 1 for B\(\,\lesssim \, 3.8\,{{{{{\rm{T}}}}}}\) and T\(\,\gtrsim \, 1.75\,{{{{{\rm{K}}}}}}\). For temperature below 1.75 K, \(z\nu\) starts to rise sharply with the field and finally it diverges while approaching the critical field \({B}_{c}^{* }\). The similar trend is observed also in the literature where, initially at high temperature, \(z\nu\) shows a plateau region (\(z\nu\)  ~ 0.5 < 1) with increasing magnetic field, but with decreasing temperature, \(z\nu\) changes very rapidly with the field and diverges while approaching the characteristic critical field \({B}_{c}^{* }\) associated with the IRFCP transition, thus, establishing the presence of QGS15,17,19,20,21,24,25. Accordingly, the region extending in temperature from absolute zero to 1.75 K above the mean field upper critical field \({B}_{c}^{{MF}}(0)\) can be denoted as the region of QGS which is highlighted by the light orange shade in Fig. 6. In this region, field variation is very weak (as it tends to saturate towards zero temperature) compared to a relatively wide variation in temperature. At the higher temperature side within the QGS region, though the exponent \(z\nu\) increases with the field but it remains less than unity (\(z\nu \, < \, 1\)). Whereas, close to the critical field \({B}_{c}^{* }\) for T ≤ 0.6 K in the QGS region, a diverging behavior of \(z\nu\) with field is observed and the condition \(z\nu \, > \, 1\) is achieved. Accordingly, the quantum phase fluctuation regime below the TTAFF (B) is extended up to the critical field \({B}_{c}^{* }\) by the purple dashed curve which separates the QGS region with \(z\nu \, > \, 1\) from that with \(z\nu \, < \, 1\).

In the QGS regime, the dominant quenched disorders lead to the formation of superconducting puddles in which the superconducting order parameter locally gets concentrated and transition to a stable macroscopic phase-coherent state can be achieved when the Josephson coupling between these rare ordered puddles becomes strong enough47. The coherent coupling between the puddles depends on the characteristic correlation length through which the phase coherence is preserved9. The correlation length varies inversely with the temperature for long range coupling7. Near QPT, the correlation length is smaller than the distance between the puddles but comparable to the size of the individual puddle and the transport is dominated by the quantum fluctuation of the superconducting order parameter and can be considered as the region of quantum fluctuation47. The existence of QGS at relatively high temperature is not commonly observed due to the smearing out of the effect of the quenched disorder by the thermal fluctuations and the rare ordered regions of superconducting puddles may face survival issues against the thermal fluctuations15. However, in the present study as determined by the direct activated scaling, a reasonably large portion in the phase diagram is designated as the QGS regime where the region with \(z\nu \, > \, 1\) represents the absolute dominance of QGS with diverging critical exponents and the portion with \(z\nu \, < \, 1\) extending up to 1.75 K indicates the onset of QGS with a significant influence from the thermal and quantum fluctuations of the order parameter. The wider span in temperature has been reported recently for nonepitaxial NbN19. Here, the appearance of QGS for a wide temperature range and with a large variation in the RN strongly indicate towards the robustness of QGS in these nonepitaxial polycrystalline TiN thin film samples that are inherently disordered by the presence of other non-superconducting elements and phases48.

Theoretically, the QGS associated with SMT/SIT has been studied by the similar superconducting droplet model where dilute but large droplets/puddles of phase-coherent rare ordered regions are separated by the disordered metallic/insulating background acting as the dissipative medium which helps to slow down the dynamics of the rare ordered regions. Here, the model assumes an isolated droplet picture where the distance between the neighboring droplets is more than the superconducting correlation length so that no percolating path is available. In order for QGS to emerge, the susceptibility of an isolated droplet should diverge while approaching the QCP. Here, the probability \({P}_{{droplet}}\) for the formation of large superconducting droplets falls exponentially with its size as, \({P}_{{droplet}} \sim \exp \left(-b{L}^{d}\right)\) whereas, the susceptibility of the individual droplet \({\chi }_{{droplet}}\) grows exponentially with the size as, \({\chi }_{{droplet}} \sim \exp \left(a{L}^{d}\right)\). L is the lateral dimension of the droplet and d refers to the space dimensionality. b relates with the disorder strength and c represents the distance from the criticality47. Therefore, the existence of QGS depends strongly on the balance between \({P}_{{droplet}}\) and \({\chi }_{{droplet}}\)49. Further, the dissipation, caused by the gapless electrons that can penetrate the superconducting droplets entirely, can change the afore-mentioned balance dramatically. The gapless electrons in the dissipative bath can penetrate these large droplets only up to a maximum distance of the coherence length. As it is argued47,50, for the droplets that are large compared to the coherence length, the coupling to the dissipative bath grows at most in proportion to the surface length in 2D (surface area in 3D) and hence, the susceptibility varies exponentially with the perimeter instead of area in 2D (the surface area rather than the volume in 3D). Thus, the critical dimensionality of the rare region is below the lower critical dimension of the problem, which leads to the conventional QPT51.

Here, the main issue relates with two length scales that are size and the coherence length of the puddles. It has been shown that the energy scale related to the rare-ordered region falls exponentially with the size of the droplet which leads to the singular density of states for anomalously large droplets52,53. Hence, the dimensional crossover for anomalously large droplets is expected to occur at extremely low temperature28,49 where the characteristic correlation length becomes larger than the separation between the superconducting puddles and the transition can still be represented by a particular universality class as described by appropriate critical exponents in the case of conventional QPT47. Whereas, for the case of QGS, there exists an optimal size which is of the order of coherence length with an optimal range for the radius of the droplets. Above this optimal range, the probability of finding a rare-ordered droplet becomes extraordinarily rare and below this range, the susceptibility of the droplets becomes small47. Hence, QGS can be expected for the droplets having the dimensions falling within this optimal range which extends up to the limit where the droplets still follow the quantum dynamics and the QGS phase arises from the quantum fluctuations of the superconducting order parameter in the droplets49,52.

Now, according to the BCS theory, the coherence length \({\xi }_{{BCS}}=\frac{\hslash {\nu }_{F}}{\pi \Delta (0)}\), becomes anomalously large while approaching the QCP where superconducting energy gap Δ(0) → 0. Here, \({\nu }_{F}\) is the Fermi velocity. Further, from the fluctuation spectroscopy studied in 2D superconductor under perpendicular magnetic field, the existence of quantum liquid phase with long coherence length for the field range \({B}_{c2} \, < \, B \, < \, {B}_{c}^{* }\) is expected due to the quantum fluctuations among the coherently rotating fluctuating Cooper pairs (FCPs) that form the rare-ordered superconducting droplets54,55. In this field region, the quantum fluctuation driven spatial coherence length \({\xi }_{{QF}}\) becomes much longer than \({\xi }_{{BCS}}\) \(\left({\xi }_{{QF}}{ \, \gg \, \xi }_{{BCS}}\right)\) and \({\xi }_{{QF}}\) increases with increasing field above the \({B}_{c2}\) till it reaches to the QCP. Here, the coherence length \({\xi }_{{QF}}\) measures the characteristic size of the superconducting droplets55. Therefore, the assumption that the gapless electron bath causing the dissipation can penetrate into the droplet is self-consistent in this region of field near the QCP where the coherence length is much larger due to the quantum fluctuations.

Conclusions

In conclusion, a systematic investigation through low temperature magnetotransport measurements in disordered TiN thinfilms of thickness in the range of 2–4 nm demonstrates field induced superconductor to metal (SMT) and/or superconductor to (weak) insulator (SIT) type QPTs that are associated with multiple crossing points (MCPs) in magnetoresistance isotherms for a number of samples varying more than two orders of magnitude in their RN values. Measurements down to 300 mK manifest a diverging behavior of the dynamical critical exponent \(z\nu\) upon approaching the QPT indicating the emergence of QGS during the transition. Further, the applicability of the activated dynamical scaling to describe the magnetoresistance isotherms for a wide range of temperature suggests the universality class of the transition which is governed by the infinite-randomness fixed critical point. The association of QGS with SMT and SIT can be distinguished by the appearance of an upward trend of Bc(T) and a saturating Bc(T), respectively. Finally, a comprehensive phase diagram is established for the SIT system which displays a wide QGS region. The robustness of the QGS in the presented nonepitaxial TiN thin films is confirmed by its existence over a wide temperature range and also over a wide variation in the disorder strength as evident in their RN values.

Methods

We have employed an intrinsic Si (100) substrate covered with Si3N4 dielectric spacer layer of 80 nm thickness for the thin films samples to be prepared on. The Si3N4 topping layer was grown by using low pressure chemical vapor deposition (LPCVD) technique and in this study, Si3N4 is the only source of nitrogen for the nitridation of Ti to produce TiN thin films56. After adopting a standard cleaning process for the Si3N4/Si (100) substrates, Ti films were deposited on the substrate by dc magnetron sputtering using a Ti target of 99.995% purity in the presence of high purity Ar (99.9999%) gas. Sputtering of Ti was performed with a base pressure less than 1.5 × 10−7 Torr. Subsequently, the sputtered Ti films were transferred in situ to an attached UHV chamber for annealing. Different batches were prepared by varying the annealing temperature (Ta) and the thickness of the Ti films. The annealing temperatures used to prepare the samples presented in this study were 780 °C and 820 °C with a variation of ± 10 °C. The annealing was done for 2 hours at a pressure less than 5 × 10−8 Torr. During the annealing process, Si3N4 substrate decomposed into Si (s) and N (g) atoms and due to high reactive nature of Ti towards both Si and N, the annealing process led to the formation of superconducting TiN as a majority phase along with metallic TiSi2 as a minority phase. The detailed chemical kinetics involved in the substrate mediated nitridation technique is reported elsewhere48. For low temperature transport measurements, thin films in multi-terminal device geometry were grown by using shadow mask made of stainless-steel. Further, electrical contact leads of Au (80-100 nm)/Ti (5 nm) were deposited by dc magnetron sputtering by using complimentary shadow masks.

Low temperature resistivity measurements from room temperature down to 2 K were carried out in a Physical Property Measurement System (PPMS) equipped with a 16 T magnet from Quantum Design at UGC-DAE CSR Indore. The same measurement system was used to measure down to 300 mK with the help of a dilution refrigerator insert/port/unit. The electrical resistivity were measured with an excitation of 100 nA.