Growth of GaP_{1 − x − y}As_yN_x on Si substrates by chemical beam epitaxy

Dilute nitrides are a type of highly mismatched alloys, in which a small N concentration replaces group V elements in III-V semiconductors. The small, highly electronegative N atoms substituting more metallic group V elements introduce highly localized states close to the conduction band edge of the host material. When the N mole fraction, x, increases above the impurity limit (typically x > 0.005), a dilute nitride alloy is formed, whose conduction band structure is drastically modified. The electronic band structure of such alloys has been described by the band anticrossing (BAC) model,¹ which considers an interaction between the localized N states and the extended states of the semiconductor matrix, whose coupling energy given by V_NM= C_NM x^1/2 increases with N mole fraction. In the frame of the BAC model, two new conduction bands result from that interaction, whose energies are given by^1,2

where E_N is the localized nitrogen level energy and E_M is the host semiconductor conduction band extended state energy. The above interaction leads not only to a splitting of the conduction band but also to a large reduction in the fundamental bandgap energy even for low N contents due to the large values observed for C_NM, usually in the range of some electron volts. In the above equation, E_M = E_M (k) can be used when calculating the new bands dispersion relation, or also E_M = E_M (T) for obtaining their temperature dependence.

Band structure modification in highly mismatched alloys led to a large interest in semiconducting alloys such as GaAs_1 − xN_x, Ga_1 − yIn_yAs_1 − xN_x, as well as GaN_xSb_yAs_{1 − x − y} or Ga_1 − yIn_yN_xSb_zAs_{1 − x − z} grown on GaAs substrates, proposed for optoelectronic and photovoltaic applications.^3–5 

While the isolated N energy level in GaAs and Ga_1 − yIn_yAs alloys lies above the conduction band Γ minimum, in the case of GaP and GaP_1 − yAs_y (y < 0.65), this level is located below the conduction band minima. The behavior of N as an isoelectronic impurity in GaP and GaP_1 − yAs_y in the P-rich side of the alloy has been largely exploited to manufacture light emitting diodes.⁶ 

GaP_1 − xN_x and GaP_{1 − x − y}As_yN_x are interesting alloys due to the proximity of GaP and Si lattice parameters. Because GaP_0.98N_0.02 is lattice matched to Si, the integration of this III-V semiconductor on Si is promising for electronic and photovoltaic applications.^7,8 The addition of As for the growth of GaP_{1 − x − y}As_yN_x layers adds further control for both the lattice parameter and the bandgap of the resulting alloy. Since GaP_1 − yAs_y is derived from GaP, it is an indirect bandgap semiconductor for low As mole fraction values (y < 0.5) where the isolated N impurity level lies below the conduction band, but N incorporation converts this material into a direct gap semiconductor with a narrow lowest conduction band derived from the N states.⁹ 

GaP based alloys like GaP_{1 − x − y}As_yN_x are difficult to grow in a controlled way because of the presence of three different anions competing for the incorporation to the group V sublattice sites during growth. It is known that As tends to displace both P and N during growth,^10,11 while P may also replace N. Therefore, large fluxes of P and N—or precursor fluxes—are needed during epitaxial growth.

Several groups have reported previous research on the growth and characterization of GaP_{1 − x − y}As_yN_x layers. While some published works focus on GaP_{1 − x − y}As_yN_x layers grown on GaP substrates by gas source MBE,¹² other papers report results on GaP_{1 − x − y}As_yN_x quantum wells grown by MBE on Si substrates, where only a very thin layer of this material forms the active layer of a quantum well.^13,14 Thick GaP_{1 − x − y}As_yN_x layers grown on Si substrates by metalorganic chemical vapor deposition have also been reported.¹¹ 

As far as we can know, previous attempts to grow Si lattice-matched GaP_{1 − x − y}As_yN_x layers by chemical beam epitaxy (CBE) have not been reported. We have used the CBE technique in this work to demonstrate its ability for the growth of high quality GaP_{1 − x − y}As_yN_x layers on Si(001) substrates with nearly lattice-matched composition; to address the issue of antiphase domain generation associated with the growth of polar semiconductors on Si, we have used misoriented Si(001) substrates with carefully designed buffer layer growth.¹⁵ 

Samples were grown on Si(001) substrates misoriented 4° toward the [111] direction. This misorientation has demonstrated to induce the formation of double steps on the Si substrate,¹⁶ thus reducing the nucleation of antiphase domains and boundaries during epitaxial growth of polar semiconductors.¹⁵ Si substrates were first etched in a diluted HF solution to remove any remaining native oxide layer before they were oxidized again under a controlled UV-ozone exposure; a second etching in the same solution, followed by deionized water rinse and blowing in N₂, was performed immediately before their introduction in the growth system.

Growth was performed by CBE in a Riber CBE32 system using triethylgallium (TEGa), tertiarybutylarsine (TBAs), tertiarybutylphosphine (TBP), and 1,1-dimethylhydrazine (DMHy) as precursors for Ga, As, P, and N, respectively, by means of a pressure-based flux control system. TBAs and TBP were injected through the same high temperature cracking cell (T = 820 °C); this setup allows for As/P flux homogenization inside the cell, but it suffers from a crosstalk pressure around 7%. TEGa and DMHy fluxes were injected through separate low temperature cells (T = 120 °C), avoiding cracking of these precursors inside their cells; Ga and N precursors are then allowed to decompose at the sample surface during growth. Sample temperature was measured using a dual band optical pyrometer.

After Si substrate cleaning at T_s = 750 °C, a (2 × 2) like surface reconstruction was initially observed by reflection high-energy electron diffraction (RHEED), due to the superposition of (1 × 2) and (2 × 1) domains on adjacent Si(001) terraces. A further annealing at the same temperature was performed, until an apparent (2 × 1) surface reconstruction is observed, indicating that terraces separated by double steps are built to accommodate the substrate misorientation. GaP buffer layer growth was then started at T_s = 530 °C, alternating TEGa and TBP fluxes, starting with TEGa exposure to avoid Si-P dimer formation.¹⁷ This migration enhanced epitaxy (MEE) procedure has been proved to grow good quality GaP buffer layers on Si substrates by MBE even at much lower temperatures than used in this work.¹⁵ However, it should be noted that growth temperatures below 500 °C are not useful for CBE growth due to the reduced TEGa cracking efficiency.

The initial GaP buffer layer growth regime turns to be 3D after 4–5 MEE cycles, as indicated by diffraction spots developing on the RHEED screen. As the buffer layer growth proceeds, RHEED pattern is progressively improved, in such a way that after the growth of 60 monolayers a faint (2 × 4) surface reconstruction is observed, whose quality is improved by a short annealing at 600 °C. The above observation suggests a monodomain GaP layer with reduced or even suppressed antiphase domains and antiphase boundaries.¹⁸ A graded buffer layer is afterwards grown on top of the GaP buffer layer, before the growth of the final epitaxial layer (either GaP_1 − xN_x, GaP_1 − yAs_y or GaP_{1 − x − y}As_yN_x) during 1 h at T_s = 550 °C.

Figure 1 is the RHEED pattern image obtained along one of the two azimuths parallel to substrate step edges after the growth of both GaP buffer layer and GaP_{1 − x − y}As_yN_x graded layer. A clear fourfold pattern is observed at this stage, with diffraction streaks that are slightly tilted respect to the surface normal—the bottom shadow represents the substrate surface projection—due to the layer and substrate misorientation angle.

Higher intensity regions in the RHEED pattern correspond to diffraction spots superimposed to the central and first order streaks, indicating that the surface planarity is not optimal.

The above pattern is usually preserved during the growth of layers nearly lattice-matched to the substrate, with the intensity of the spots being progressively reduced as the layer thickness is increased, accompanied by the appearance of a weak specular spot, indicating the good quality of the grown layers.

Large V/III flux ratios were used during growth, keeping always the same TBP flux. Different TBAs flux values were injected only when growing GaP_1 − yAs_y or GaP_{1 − x − y}As_yN_x layers, in order to control the incorporated As mole fraction. Sample series were grown at two different growth rates, either 0.1 or 0.3 μm/h, by changing TEGa pressure in the gas line. A short and moderate postgrowth annealing at T_s = 650 °C for 5 min under V element precursor fluxes was performed in the case of GaP_{1 − x − y}As_yN_x samples.

Samples were analyzed by high-resolution X-ray diffraction (HRXRD) in a Bede D3 diffractometer, using a 4-bounce double Si crystal monochromator, and a 0.5 mm slit as an analyzer, in order to check their crystalline quality.

Sample composition was obtained from Rutherford backscattering spectrometry (RBS) and nuclear reaction analysis (NRA) spectra obtained in a tandem accelerator using a collimated beam of 2 MeV ⁴He⁺ ions for RBS and 3.72 MeV ⁴He⁺ ions for NRA to induce the resonant ¹⁴N(α,p)O¹⁷ nuclear reaction. During NRA acquisition, one detector was covered by an 18 μm thick mylar foil, in order to separate the proton signal from the large background of backscattered α particles. A second uncovered detector was used for simultaneous RBS acquisition at 3.72 MeV in order to obtain a good indicator of the particle dose. As, P, and N mole fractions in the epitaxial layers were obtained after simultaneous fitting of three spectra (RBS at 2 and 3.72 MeV and NRA at 3.72 MeV).

Photoluminescence (PL) spectra were acquired with the sample introduced in a temperature controlled closed-cycle He cryostat (between 12 K and room temperature), using a focused multiline Ar⁺ laser with a 200 mW power output (estimated power density of 1 kW/cm²), a 40 cm focal length monochromator, and a dual Si-GaInAs photodetector operating with a lock-in amplifier and a chopper.

Bandgap calculations for GaP_1 − xN_x, GaP_1 − yAs_y or GaP_{1 − x − y}As_yN_x have been done using the wide set of III-V semiconductor parameters commonly accepted,^19,20 the model-solid theory²¹ and the band anticrossing model for N-containing alloys.¹² 

GaP_{1 − x − y}As_yN_x may be grown lattice matched to Si substrates in a wide range of N and As mole fractions. Figure 2 represents the lattice parameter of this quaternary alloy as a function of the N, As, and P mole fractions, as calculated by the Vegard law from the lattice parameters of the three involved end binaries,

The dashed line indicating the value of the Si substrate lattice parameter [a(Si) = 5.431 Å] is also drawn over the ternary diagram, showing the allowed range of GaP_{1 − x − y}As_yN_x compositions for Si lattice matching. Both ternary end values are also displayed, being the only ternary compounds allowed for Si lattice matching in this system.

More explicitly, the GaP_{1 − x − y}As_yN_x interval of compositions allowing lattice matching to the Si substrate may be obtained from Eq. (2) as

meaning the As mole fraction needed to compensate N induced stress is nearly 5 times (or 4.69) the N excess over the lattice matching value for GaP_1 − xN_x (given by x = 0.02).

Plotted in Fig. 3 is the calculated evolution of the three conduction band minima (CB) associated with Γ, X, and L symmetry points of the Brillouin zone for the ternary alloy GaP_1 − yAs_y as a function of As mole fraction (bottom axis) at T = 300 K. The vertical axis energy origin in this plot is not meaningful, as it was arbitrarily assigned to the energy of the InSb valence band maximum.¹⁹ The averaged energy of the three valence bands or valence band offset (VBO), the valence band maxima for heavy and light holes [VB (hh-lh)], and spin–orbit split-off band [VB (SO)] are also depicted. The position of the N energy level (E_N) is also included in this figure; note that this level moves into the bandgap for y < 0.65, near the crossover bandgap transition from indirect to direct (y = 0.5).

The evolution of the E₋ conduction band minimum for Si lattice-matched GaP_{1 − x − y}As_yN_x, resulting from the coupling of N localized states with Γ minimum of the host GaP_1 − yAs_y conduction band, as predicted from the BAC model is also included in Fig. 3, where the upper axis represents the N mole fraction, x.

It is important to point out that for BAC calculations, the host matrix of the quaternary alloy GaP_{1 − x − y}As_yN_x has been considered as the GaP_{1 − y′}As_y′ ternary alloy whose As and P mole fractions keep the same ratio than in the quaternary alloy,

Because the coupling of N states near the Brillouin zone center has shown to be stronger than the coupling near Brillouin zone boundaries,²² the E₋ bands associated to N coupling with X and L extrema are expected to appear at energies higher than the E₋ values plotted in Fig. 3.

The calculated energy decrease of the E₋ minimum (and bandgap energy) is very pronounced as As and N are added to GaP, to form Si lattice-matched GaP_{1 − x − y}As_yN_x. When compared to GaP_1 − xN_x for the same N mole fraction (calculations not shown here), the E₋ curve plotted in Fig. 3 decreases faster, due to the extra bandgap reduction induced by As incorporation to the GaP_1 − yAs_y host matrix. In this way, the bandgap reduction in Si lattice-matched GaP_{1 − x − y}As_yN_x increases with incorporation of both N and As.

The As mole fraction of the epitaxial layers, obtained by RBS, as a function of the TBAs pressure used during growth in the corresponding gas feeding line is plotted in Fig. 4. Three different data sets are drawn in this figure. Red squares and green circles denote obtained y values for GaP_1 − yAs_y layers grown at rates of 0.1 μm/h and 0.3 μm/h, respectively, while blue diamonds denote y values obtained for GaP_{1 − x − y}As_yN_x layers grown at 0.3 μm/h.

In order to compensate the possible P displacement by As, a large TBP/TBAs flux ratio was used (compare values in the inset and in horizontal axis of Fig. 4), assuring the availability of a much larger P excess on the sample surface during growth. The linear dependence of the As mole fraction on TBAs pressure during growth is clearly seen in this plot, independently of the growth rate.

The addition of DMHy during growth to incorporate N does not seem to affect the As mole fraction incorporated in GaP_{1 − x − y}As_yN_x samples (blue diamonds in Fig. 4). The straight line in this figure shows the fitted linear dependence, intersecting the horizontal axis at 0.13 Torr, value corresponding to the 7% crosstalk pressure at the cell input produced by the TBP gas line working at 1.7 Torr. TBAs pressure in its gas line should be larger than the above crosstalk value in order to inject a useful TBAs flux in the common cell. When TBAs pressure is increased over the crosstalk pressure, As incorporation is controlled by the TBAs gas line pressure excess over the crosstalk threshold value, giving the linear dependence observed in Fig. 4.

The dependence of As mole fraction on TBAs pressure—or flux—should leave the linear behavior for large TBAs values, as y should asymptotically approach the unity value. This behavior is not seen in Fig. 4, because used TBAs fluxes are not large enough to obtain As mole fraction values larger than y = 0.35, which are probably too low to observe a noticeable change on the slope of the data set.

The N mole fraction of the epitaxial layers, x, obtained from the simultaneous fitting of RBS and NRA spectra, as a function of the DMHy pressure used during growth in the corresponding gas feeding line is plotted in Fig. 5. Two different data sets are drawn in this figure. Red squares denote x values for GaP_1 − xN_x layers grown at 0.1 μm/h, while blue diamonds denote x values obtained for GaP_{1 − x − y}As_yN_x layers grown at 0.3 μm/h.

The data in Fig. 5 indicate a generally linear dependence of N mole fraction on the DMHy gas line pressure. N content measured by RBS-NRA combined techniques shows a much larger data scatter than observed for As content in Fig. 4. Two samples even clearly deviate from the general trend. This behavior could be related to some uncontrolled growth parameter not yet identified, such as a substrate temperature too low for a full DMHy decomposition at the sample surface, or sample-to-sample temperature variations due to pyrometer window metallization. No noticeable difference is observed in Fig. 5 on N mole fraction values obtained for ternary or quaternary alloys, neither the growth rate seems to influence N incorporation.

As and N mole fractions measured for the three sample sets studied in this work—both ternaries GaP_1 − xN_x and GaP_1 − yAs_y, and the quaternary GaP_{1 − x − y}As_yN_x—are represented in Fig. 6. Experimental points obtained for both ternaries stay on both axes, as they correspond to y = 0 and x = 0 values, respectively.

Also drawn in this plot is the Si lattice matching line for the quaternary alloy, given by Eq. (3), and the lines for lattice mismatch ɛ = ±0.005 between the lattice parameters of the quaternary epitaxial layer, a_e, and of the Si substrate, a_s, defined as

Lattice mismatch values in the range |ɛ| ≤ 0.005 usually allow the growth of epitaxial layers up to or even thicker than 200 nm without strain relaxation. As can be seen in Fig. 6, the composition of quaternary samples studied in this work is very close to the lattice matching values, always inside the |ɛ| ≤ 0.005 stripe. The data spread observed for N mole fraction in Fig. 5 is also observed for GaP_{1 − x − y}As_yN_x in Fig. 6; whereas TBAs pressure was increased proportionally to DMHy pressure in the sample set in order to obtain near Si lattice-matched layers, As mole fraction increased proportionally (see Fig. 4), but N mole fraction was not tightly controlled during growth.

Figure 7 represents the reciprocal space map (RSM) obtained by HRXRD near the Si [004] diffraction peak of an illustrative epitaxial sample including a GaP_1 − xN_x layer on top of a GaP buffer layer grown on a Si substrate. Incident X-ray beam projected azimuth was parallel to step edges, and sample was tilted by the misorientation angle to bring the Si [004] peak to the diffraction plane.

Three peaks are clearly seen in this figure, which are slightly elongated due to the use of a 0.5 mm slit as an analyzer. The three peaks observed in the RSM are identified—from top to bottom of the diagram—as due to diffraction by (004) planes of the Si substrate, GaP_1 − xN_x, and GaP layers, respectively.

The topmost and most intense peak is associated with the Si [004] diffraction peak, due to its smaller lattice parameter. The lower and less intense peak is due to the thin GaP buffer layer, having the larger lattice parameter of the three layers. The three observed peaks are aligned along the vertical axis over the h = k = 0 position, as any possible epilayer small tilt induced by substrate misorientation cannot be observed under this configuration.

The positions of the two peaks related to Si substrate and GaP buffer layer are insensitive to N mole fraction, while the position of the central peak of Fig. 7, associated with the GaP_1 − xN_x layer, moves its position when x is changed from sample to sample (or also y for GaP_1 − yAs_y and GaP_{1 − x − y}As_yN_x layers). Starting at the GaP position for x = 0, the GaP_1 − xN_x related peak moves upwards when increasing the N content until it overlaps the Si peak at l = 4 for x = 0.02; for larger x values, the GaP_1 − xN_x would be located in the topmost position of the RSM. The relative position of GaP_1 − xN_x peak with respect to GaP and Si peaks allows then to calculate N content in this layer if the relaxation degree of the layer is known.

GaP buffer layer in this map is a good test for structure relaxation information. The vertical peak position in the RSM for a fully relaxed GaP layer may be calculated (in Si [hkl] units) as

giving the value l = 3.986, while for a fully compressively strained GaP layer, it may be calculated using the basic elastic theory as

where C₁₁ and C₁₂ are the elastic constants of the epitaxial layer, giving the value l = 3.973. As the GaP peak in Fig. 7 is very close to this last value, it can be concluded that GaP layer and perhaps the whole epitaxial film are fully strained. In such a case, the position of the GaP_1 − xN_x peak can be used in Eq. (7) to calculate the lattice parameter a_e, obtaining the N content by Vegard's law. For this calculation, the GaP values of the epitaxial layer elastic constants may be used in a first approach, due to the expected low x value; a further refinement by iteration using Vegard's law for C_ij may also be done. Following that procedure, the value x = 0.015 was obtained for the sample shown in Fig. 7, which is only slightly smaller than the value x = 0.019 obtained by combined RBS-NRA measurements. This discrepancy may be attributed to the smaller accuracy of ion-beam based techniques when compared to HRXRD, or to the fact that RBS-NRA is sensitive to the total amount of N present in the layer, while HRXRD is mainly sensitive to substitutional N.

The well defined shape and width of the GaP_1 − xN_x peak for this very thin (100 nm) film confirms the good crystalline quality of these samples, as suggested from observed RHEED patterns during growth. Similar results are observed on GaP_{1 − x − y}As_yN_x samples, but due to the contribution of both x and y mole fractions to the layer diffraction peak position, no straightforward information about sample composition may be obtained, unless another supplementary technique be used; moreover, for GaP_{1 − x − y}As_yN_x layers nearly lattice matched to the Si substrate, the layer [004] peak is almost hidden under the stronger Si [004] peak.

A series of PL spectra measured for a GaP_0.88As_0.09N_0.03 sample at temperatures between 12 K and 250 K is shown in Fig. 8. These spectra are representative for GaP_1 − xN_x and GaP_{1 − x − y}As_yN_x sample sets grown in this work. Two main emission peaks or bands are observed in every spectrum, the main one being observed in the high-energy side of the spectrum. This peak is probably due to near bandgap recombination, close to the E₋ energy of the BAC model. The full width at half maximum (FWHM) of the main peak is estimated to be in the range 120–130 meV for all the samples studied in this work. The peak energies of both features move toward lower energy values when increasing temperature because of alloy bandgap reduction, as it is typically expected and observed in semiconductors. The integrated intensity of both peaks decreases as the temperature is increased from T = 12 K up to room temperature, being the low energy peak the predominant feature of the PL spectra at high temperatures.

The nature of both peaks, separated about 150 meV, remains unresolved and needs further investigation, as one of them could arise from residual carbon doping, always present in the samples grown in our CBE system under no intentional doping growth conditions.²³ Another possible origin for these peaks may be related to localized states in the conduction band, due to N concentration fluctuations. This kind of localized states sometimes exhibits a typical “S” shape in the low temperature shift of the peak^24,25 that has not been observed in our samples.

The detection of PL signal in both series of samples, GaP_1 − xN_x and GaP_{1 − x − y}As_yN_x, is a clear evidence of the direct nature of the fundamental bandgap, as a consequence of the strong coupling between N localized states and the direct Γ minimum of the conduction band. In contrast, no detectable PL signal has been observed in the set of GaP_1 − yAs_y samples (y ≤ 0.35) studied in this work (see Fig. 6). The large value obtained for the FWHM of the measured PL peaks could arise from different factors. It could be attributed to local N mole fraction fluctuations in each sample, inducing bandgap variations broadening the PL peaks. Unresolved valence band splitting due to residual biaxial strain on the epitaxial layers could also contribute to peak broadening.

Figure 9 represents the position of the higher energy PL peak for GaP_1 − xN_x and GaP_{1 − x − y}As_yN_x sample series at T = 12 K as a function of the N mole fraction obtained by combined RBS-NRA. The bandgap energy due to the E₋ band build-up has been calculated at the same temperature using the BAC model, also represented by continuous lines. Bandgap modification induced by biaxial strain in GaP_1 − xN_x due to the lattice mismatch with the Si substrate has not been taken into account as this effect, including valence band splitting, has been estimated under 30 meV in the composition range shown in Fig. 9.

Measured PL peak energies are close to the predicted values, but again some data spread is clearly observed. The limitations on the accuracy for N mole fraction determination should be the most significant factor responsible for data spread. The lack of a hard postgrowth annealing further homogenizing N content in the samples, and probably reducing PL peak width, could also be a second source for data spread.

Data and calculations plotted in Fig. 9 confirm that the observed higher energy PL transition is related to a near bandgap recombination. Results for the GaP_1 − xN_x samples are generally closer to the calculated values than for GaP_{1 − x − y}As_yN_x sample set. This larger data spread in the case of quaternary samples can be explained from the results shown in Fig. 6, because As and N mole fractions do not fit to the values for exact Si lattice matching required by Eq. (3), due to the lack of a precise control on x.

Apart from the above considerations, it should be noted that the difference between calculated bandgap and measured PL peak energy is generally smaller than the large value of the FWHM of measured PL peaks (in the range from 120 to 130 meV).

Good quality epitaxial layers of GaP_1 − xN_x, GaP_1 − yAs_y, and GaP_{1 − x − y}As_yN_x have been successfully grown by chemical beam epitaxy at 550 °C on Si(001) substrates with a 4° misorientation toward the [111] direction. As and N mole fractions of the grown layers, obtained by combined RBS and NRA measurements, show a linear dependence on the flux of their precursor, TBAs or DMHy, respectively. In the case of N mole fraction, some data spread has been observed, probably due to the lack of precise control over some growth parameter or to the accuracy of the measurement technique. The GaP_{1 − x − y}As_yN_x layers grown on Si substrates show a lattice mismatch not larger than |ɛ| = ±0.005 for N contents in the range 0.02 < x < 0.05.

RHEED and HRXRD reciprocal space mapping confirms the good crystal quality of the epitaxial layers. Strong PL signal has been detected in GaP_1 − xN_x and GaP_{1 − x − y}As_yN_x layers even at room temperature, confirming the predicted direct bandgap for these alloys. PL spectra show two wide peaks, the one at the higher energy being most likely related to near bandgap recombination. The position of this peak has been compared to the bandgap energy predicted by BAC as a function of x in the layers; there is good agreement, taking into account some data spread and the large peak width.

This work was supported by the Spanish Ministry of Science, Innovation and Universities, under Project No. TEC2016-78433-R. One of the authors (B. J. García) acknowledges a scholarship from the Spanish Ministry of Education under agreement with the Spanish Fulbright Commission. Technical assistance for ion beam measurements at Centro de Microanálisis de Materiales (CMAM, Madrid, Spain) and Lawrence Berkeley National Laboratories (LBNL, Berkeley, CA, USA) from M. D. Ynsa and Y. Lee, respectively, is appreciated. The work performed at LBNL was supported by the Director, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The work of N. López was partially supported by the Grant ERC-Stg-2017-758885 from the European Commission and by the Ramon y Cajal Program from the Spanish Ministry of Education.