On the nature of as-grown and irradiation-induced Ga vacancy defects in β-Ga₂O₃

β-Ga₂O₃ is an attractive next-generation semiconductor for high-power and high-temperature electronics due to its excellent physical properties, such as an ultra-wide bandgap of (4.6–4.9 eV) and a high empirical estimated critical electric field (6–8 MV/cm).^1–4 It has the advantage of being available as inexpensive bulk substrates that can be grown from the melt in a large size via Czochralski,⁵ edge-defined film-fed growth,⁶ and floating zone methods.⁷ Relatively high electron mobility as well as excellent thermal and chemical stability have been achieved.⁸ Controlled n-type doping is routinely performed by adding group IV elements (Si, Ge, and Sn).^2,9,10

The presence of point defects in the lattice, such as vacancies, interstitials, and impurity atoms, which can appear during the material synthesis and/or post-growth processing, strongly affects the electrical and optical properties of semiconductors. Theoretical calculations predict that Ga vacancy-related defects exhibit very low formation energies thanks to their multiply negative charge states in the upper half of the bandgap of β-Ga₂O₃.⁹ The monoclinic β-Ga₂O₃ lattice hosts two distinct Ga sites: a fourfold coordinated Ga(1) site and a sixfold coordinated Ga(2) site. In addition, the lowest-energy configurations of the two types of vacancies are the so-called split configurations $V Ga ia$⁠, $V Ga ib$⁠, and $V Ga ic$⁠, where a fourfold coordinated Ga(1) atom has relaxed inward into the interstitial space creating a structure with two “half-vacancies” with an interstitial in between.⁹ 

The split configurations appear as the dominant form of existence of Ga vacancies in β-Ga₂O₃, as found in a wide variety of experiments: FTIR spectroscopy, EPR spectroscopy, positron annihilation spectroscopy, and transmission electron microscopy.^11–24 We note that these split Ga vacancies are found to form complexes with O vacancies and hydrogen impurities.^9,12,20,25 Limited experimental evidence has been presented on unrelaxed Ga vacancies, although they have been suggested as the defect created in particle irradiation.^11,26 Electrical and optical characterization performed in β-Ga₂O₃ after irradiation with protons, electrons, or neutrons finds V_Ga-related defects.^27–33 The V_Ga defects introduced by the irradiation have been interpreted to be present as complexes with O vacancies and/or H impurities, but revealing structural details with these experiments is difficult.

In this work, we present results obtained with 4D (3D orientation + temperature dependence) Doppler broadening spectroscopy and temperature-dependent positron lifetime spectroscopy in Fe, Mg, and Sn doped β-Ga₂O₃ crystals before and after irradiation with 6 MeV protons at room temperature. The unusually strong overall anisotropy of the Doppler broadening signals combined with the positron lifetime values across the 30–600 K temperature range confirms the presence of split Ga vacancies in various configurations and complexes in the as-grown crystals. Interestingly, this unusual anisotropy is reduced after the 6 MeV proton irradiation, accompanied by an increase in the positron lifetime. These are interpreted as clear evidence of unrelaxed Ga vacancies being introduced as radiation damage at room temperature.

This article is organized as follows. Section II describes the experimental approach in detail. In Sec. III, we present the key points of our experimental results obtained with 4D Doppler broadening spectroscopy and positron lifetime spectroscopy. In Sec. IV, we discuss the experimental results in terms of split and unrelaxed Ga vacancies by comparing them to the positron annihilation signals predicted by first-principles calculations.

The β-Ga₂O₃ single crystals were grown by the edge-defined film-fed growth (EFG) and Czochralski (CZ) methods.^5,6 The EFG crystals were doped during growth with either Fe for semi-insulating behavior (resistivity > 10¹⁰ Ω cm) or Sn for n-type behavior (free electron concentration 2 × 10¹⁸ cm⁻³). The CZ crystals were doped with Mg for semi-insulating behavior (resistivity > 10⁵ Ω cm). For further details on the material and the initial positron annihilation results on as-grown samples, see Ref. 13. The samples were subjected to 6 MeV proton irradiation at room temperature to a fluence of 1 × 10¹⁶ cm⁻². Based on SRIM simulations,³⁴ the projected range of the implanted hydrogen is 130 μm in β-Ga₂O₃ and the average concentration of irradiation-induced V_Ga on the ion track is expected to be 5 × 10¹⁸ cm⁻³. The exponential stopping profile of positrons injected directly from a ²²Na source into β-Ga₂O₃ leads to less than 0.1% of the positrons reaching this depth; hence, the positron experiments probe the ion track. Irradiating with 5–10 MeV protons has been shown to be a powerful way to create a rather homogeneous profile of mono-vacancy defects for positron studies in metals and semiconductors.^35–39 

We characterize the vacancy defects in the as-grown and irradiated β-Ga₂O₃ single crystals with positron annihilation spectroscopy.^40–42 We employ the so-called fast positron approach, where the positron source in the form ²²NaCl wrapped in 1.5 μm thick Al foil is sandwiched between two identical sample pieces. We measured the positron lifetime with a conventional temperature-controlled digital spectrometer with two collinear plastic scintillation detectors, with a timing resolution of 260 ps (FWHM). The measurements were performed in the temperature range of 30–600 K for as-grown samples and 30–300 K for irradiation samples to avoid thermal recovery of the irradiation-induced defects. We recorded 1 × 10⁶ counts in each spectrum. Annihilations in the 1 MBq source (210 ps—1.7%; 400 ps—1.5%; and 1500 ps—0.1%) were subtracted from the data before analysis where we employed the standard approach of fitting one to two exponential decay components.⁴¹ 

The unusually strong anisotropy, representing the dependence of data obtained for different measurement directions relative to the crystal orientation, of the Doppler broadening of positron–electron annihilation radiation in β-Ga₂O₃ is discussed in detail in the seminal work by Karjalainen et al.¹² In most other semiconductors, this anisotropy is small enough to be ignored in the analyses and interpretations, but in β-Ga₂O₃, it is larger than the differences between the positron annihilation signals in the delocalized state in the lattice and the localized states at various vacancy defects.^12,13,26,43 This leads to the necessity of determining the main crystallographic orientations of the samples ([100], [001], and [010]) and carefully aligning the samples concerning the high-purity Ge detector in the experiments. Karjalainen et al. developed the so-called 3D-Doppler approach where the sample–source–sample sandwich is continuously rotated with respect to the detector to map the changes in the annihilation signals. Here, we expand this approach to “4D Doppler,” adding the measurement temperature as another variable. In the interest of experiment time, measurements were performed only in the three main crystallographic orientations.

The energy resolution of the HPGe detector was 1.15 keV (FWHM) at 511 keV, and 1 × 10⁶ counts were accumulated in each spectrum. The Doppler broadening of the annihilation radiation is characterized by the conventional S and W parameters that are used to describe the shape of the spectrum. The S parameter was defined as the fraction of the counts in the central region of the annihilation line (corresponding to electron momenta 0–0.45 a.u.), and the W parameter is the fraction of counts in the wing region of the annihilation line (corresponding to electron momenta 1.5–4.1. a.u.).⁴¹ In addition, as proposed in Refs. 12 and 13, we also monitor a narrower version of the W parameter, W2, defined here as the fraction of counts in the 2.0.–4.1. a.u. momentum range. The sample experimental conditions were applied in the 4D Doppler experiments as in the positron lifetime experiments.

The average positron lifetime as a function of temperature in as-grown and 6 MeV proton irradiated EFG:Fe, CZ:Mg, and EFG:Sn samples is shown in Fig. 1. The data for the as-grown samples are discussed in some detail already in Ref. 13. The main features to be noted are the decreasing trend of the average positron lifetime above 300 K in the CZ:Mg and EFG:Fe samples and the increasing trend in the same temperature range in the EFG:Sn sample. In addition, it is instructive to observe that at room temperature, the difference between the three samples is at its smallest, highlighting the necessity to perform temperature-dependent experiments. This is of particular importance as the lifetime spectra could be analyzed with only one decay component at any of the temperatures. These observations together with the room-temperature 3D-Doppler results led to the conclusion in Ref. 13 that the as-grown samples contain high concentrations (>1 × 10¹⁸ cm⁻³) of various kinds of split Ga vacancies, leading to saturation trapping of positrons. It was also found that the EFG:Fe appears to have the lowest concentration of these defects without full saturation at elevated temperatures.

After the 6 MeV proton irradiation, the average positron lifetime clearly increases in the EFG:Fe sample by 5–6 ps below 200 K and by more than 10 ps at 300 K. The average positron lifetime in the EFG:Sn after irradiation is essentially identical to the lifetime in the EFG:Fe after irradiation, with the lifetime increasing above 200 K. However, there is no change at the low temperatures compared to the as-grown state of EFG:Sn. The increase in the average positron lifetime is a clear indication of additional vacancy defects created in the 6 MeV proton irradiation in the EFG:Fe and EFG:Sn samples. In the CZ:Mg sample, the average positron lifetime does not change in the irradiation at temperatures below 200 K, and only a very slight increase can be seen at 200–300 K. However, as the irradiation conditions are the same for all the samples, the same concentration and nature of irradiation-introduced defects should be expected. We conclude that the concentration of pre-existing vacancy defects is significantly higher in the CZ:Mg sample than in the EFG:Fe and EFG:Sn samples and also significantly higher than the concentration of irradiation-induced defects. Hence, in the CZ:Mg sample, positrons annihilate predominantly as trapped at the in-grown vacancy defects even after irradiation. Only a single lifetime component can be resolved in the spectra also after irradiation in all of the samples, except the 300 K point measured in the 6 MeV proton irradiated EFG:Fe sample. There, a second lifetime component of τ₂ = 244 ± 20 ps with an intensity of I₂ = 5 ± 1% could be extracted.

Figure 2 shows the temperature dependences of the S, W, and W2 parameters measured along the three crystallographic orientations [100], [001], and [010] before and after 6 MeV irradiation in the EFG:Fe, CZ:Mg, and EFG:Sn samples. As with the average positron lifetime, the change in the temperature for each individual sample is relatively small in all of the parameters. In all the samples, the S parameter increases in all the samples in all measurement directions after irradiation. This is conventionally interpreted as an increase in vacancy concentrations, as the S parameter is more sensitive to the open volume and the W parameter to the identities of the atoms surrounding the annihilation site. Correspondingly, the W and W2 parameters decrease after irradiation in the EFG:Fe and EFG:Sn samples. Interestingly, the W and W2 parameters increase after irradiation in the CZ:Mg samples. We note also that the difference between as-grown and irradiated data is the smallest in the CZ:Mg samples, as in the case of positron lifetime.

To better analyze the changes in the S, W, and W2 parameters in the three crystallographic orientations as a function of temperature in the as-grown samples, we present the data from three distinct temperature ranges (low 30–200 K, mid 200–400 K, and high 400–600 K) in Fig. 3. The presented parameters are obtained by averaging in the respective temperature ranges, giving a single value for each temperature range. The smallest overall anisotropy is observed in the EFG:Fe samples across the entire temperature range, while the CZ:Mg samples exhibit the largest overall anisotropy with the [010] always located at the furthest lower-right corner. With increasing temperature, the overall anisotropy is reduced in all the samples. We define overall anisotropy here as the relative difference between the largest and smallest parameter values spanned by the data obtained in each sample. As typically the vacancy values exhibit a larger S and smaller W parameter, the reference point is taken as the one with the smallest S and largest W parameter, resulting in positive relative differences in S parameters and negative relative differences in W parameters. Table I shows the numerical values of the overall anisotropies for as-grown and irradiated samples.

The red arrows in Figs. 3(a)–3(d) highlight that the order of the [010] points for all samples remains the same at all temperatures and for both the (S, W) and (S, W2) parameters. It can be seen that in some cases, the anisotropy is more “triangular” in shape [e.g., EFG:Sn in Fig. 3(c)] and in some cases more “linear” [e.g., EFG:Fe in Fig. 3(e)]. The differences are subtle but in some cases systematic. For EFG:Fe, the anisotropy is more linear in (S, W2) and more triangular in (S, W) at all temperatures. A similar observation can be made for EFG:Sn at low and high temperatures, although the effect appears even more subtle. In addition, EFG:Sn exhibits a more triangular anisotropy in the mid-temperature range compared to the other samples, observed in both (S, W) and (S, W2). We also take note of the drop of the [010] point of the EFG:Sn sample in the (S, W2) parameter [Fig. 3(f)]. The shape of the anisotropy in CZ:Mg behaves oppositely to the two other samples: it is more triangular anisotropy at all temperatures in (S, W2) than in (S, W).

Figure 4 shows the comparison between as-grown and irradiated 4D Doppler data for the EFG:Fe, CZ:Mg, and EFG:Sn samples. The averaging over the temperature ranges is similar to that employed in Fig. 3, with only two ranges (30–200 and 200–300 K) for the irradiated samples. The changes in the (S, W) and (S, W2) parameters caused by the defects introduced in the irradiation are clear. The overall anisotropy shrinks in all samples, but most significantly in EFG:Fe at both low and mid-temperature ranges. A clear drop in the (S, W2) plot is seen for EFG:Fe, and in addition, the data for [100] and [001] directions have essentially merged. In the EFG:Sn samples, all irradiated data are shifted downward and significantly to the right. A characteristic drop in the (S, W2) plot also appears in EFG:Sn after the irradiation, accompanied by a strong shift of the [010] point toward the lower right. The effect of the irradiation is stronger in the mid-temperature range than at low temperatures. In the case of the CZ:Mg samples, the data shift somewhat upward in both (S, W) and (S, W2) plots. The shape of the anisotropy also changes after irradiation in the CZ:Mg samples and becomes more linear in both (S, W) and (S, W2), but overall, the changes caused by the irradiation are smaller than in EFG:Fe and EFG:Sn.

First-principles theoretical calculations have been employed to predict the (S, W) and (S, W2) parameters for the elementary vacancy defects and a handful of their complexes with H atoms and oxygen vacancies, published in Refs. 12, 13, and 26. We have compiled them in Fig. 5 to allow comparison with the present experimental data and more detailed interpretation. In the figure, the (S, W, W2) parameters are shown as normalized by the parameters (S, W, W2) of β-Ga₂O₃ in the [001] lattice direction and the direction of the smallest calculated S parameter in the β-Ga₂O₃ lattice. It can be clearly seen that the greatest overall anisotropy belongs to split vacancy configurations $V Ga ib$ and $V Ga ic$ as well as to their complexes with hydrogen, $V Ga ib−1 H$⁠, and $V Ga ic−1 H$⁠. Both unrelaxed Ga vacancies $V Ga 1$ and $V Ga 2$⁠, as well as $V Ga ia$ and the larger complexes with oxygen vacancies $V Ga ib−2 V O$ and $V Ga ia−3 V O$⁠, show the smallest overall anisotropies, comparable to that of the lattice. The anisotropy of complexes with two hydrogen atoms $V Ga ib−2 H$ and $V Ga ic−2 H$ is smaller than their “pure” split configurations $V Ga ib$ and $V Ga ic$⁠. The shape of the anisotropy of the split configuration $V Ga ia$ is significantly different from $V Ga ib$ and $V Ga ic$⁠. It should also be noted that there is clearly more overlap in the parameters of the different vacancy defects and the lattice (S, W) parameters than in (S, W2) parameters. Similarly, as in other metal–oxide semiconductors,⁴⁴ the distinction between a metal vacancy and a cation–anion vacancy pair is difficult as these produce almost identical positron annihilation data. Only when complexes with two or more O vacancies do the (S, W) and (S, W2) parameters of the split Ga vacancies change significantly, shifting to the right. It is also worth noting that there is a “triangular” vs “linear” shape of the anisotropy, and also, the difference between (S, W) and (S, W2) in this respect is the difference for the different defects. For example, the lattice data are mostly linear in both (S, W) and (S, W2). Interestingly, $V Ga ib$ becomes more triangular when going from (S, W) to (S, W2), while $V Ga ic$ becomes more linear. The anisotropy of the unrelaxed V_Ga does not change shape significantly.

In addition, positron lifetimes have been predicted using theoretical calculations.^12,13 The main findings are that the split Ga vacancies exhibit positron lifetimes 25–35 ps longer than that in the β-Ga₂O₃ lattice. Complexes with O vacancies or H atoms do not change the predicted lifetime significantly: a single additional O vacancy increases the lifetime by 3 ps while adding an H atom to the vacancy reduces the lifetime by the same amount. Unrelaxed Ga vacancies on the other hand exhibit clearly longer positron lifetimes, 55 ps longer than in the lattice. We note that the latter is close to typical values for metal vacancies in oxide semiconductors. The split Ga vacancies have a reduced apparent open volume.

Comparison of the experimental observations and theoretical predictions allows us to make several interpretations concerning the nature of the Ga vacancies in the as-grown and irradiated β-Ga₂O₃ crystals. First, we note that even the smallest overall anisotropy found in experiments in the as-grown crystals, measured in the EFG:Fe crystal at high temperatures, is significantly larger than that predicted for the β-Ga₂O₃ lattice: 2.7% vs 2.0% for the S parameter and −9% vs −3% for the W parameter. This confirms the earlier interpretation¹³ that the positron data obtained in any of the as-grown crystals are dominated by split Ga vacancies, as these are the only defects with sufficiently large anisotropies predicted in calculations. Furthermore, the slight reduction of the anisotropy with increasing temperature in the EFG:Fe samples is consistent with the behavior of the positron lifetime as a function of temperature. Hence, we confirm the interpretation from Ref. 12 that at high temperatures, the positron trapping is no longer at saturation in EFG:Fe samples, with part of the annihilations originating from the β-Ga₂O₃ lattice. This conclusion places an upper limit of the concentration of the negatively charged vacancy defects in the EFG:Fe samples at 5 × 10¹⁸ cm^–3 as higher concentrations would lead to full saturation also above 300 K.⁴¹ The lower limit for saturation trapping at low temperatures is roughly 1 × 10¹⁸ cm^–3. Judging by the fact that the anisotropy in the EFG:Fe samples is more linear than triangular in the (S, W2) plot when compared to the (S, W) plot strongly suggests that the split Ga vacancies are of the $V Ga ic$ type in these samples. Whether they are complexed with a (single) H atom or an O vacancy cannot be resolved from the data. We note that the $V Ga ic$ split Ga vacancy is also predicted to be the lowest-energy configuration of the Ga vacancy in β-Ga₂O₃.⁴⁵ 

At low and high temperatures, the Doppler broadening data appear similar to those in EFG:Sn, but the positron lifetime at the same temperatures is longer. This suggests that the split Ga vacancies are of a similar nature, and the longer lifetime suggests complexes with O vacancies. Also, the split Ga vacancy—Sn complex, where Sn substitutes for the Ga atom in the interstitial position, is predicted to produce a slightly longer lifetime than the “clean” split Ga vacancy.²⁶ Hence, some of the split Ga vacancies may be associated with Sn. We note that the comparison between scanning transmission electron microscopy (STEM) and deep-level optical spectroscopy (DLOS) data revealed an increase in the concentration of split V_Ga-Sn defect vacancies with increasing Sn doping in the EFG:Sn substrate.²⁴ In the mid-temperature range, the positron lifetime drops, and the anisotropy is both reduced and more triangular compared to low and high temperatures. This could be interpreted as a contribution from $V Ga ia$ split Ga vacancies based on comparison with theoretical predictions, but this suggestion requires further experiments to be verified. We note that the temperature-dependent positron data do not indicate that any of these defects would be in the negative charge state. This is in line with the samples being n-type and not compensated. The varying efficiency of positron trapping at these defects at different temperatures is likely to be caused by thermal escape from the split Ga vacancies that have relatively small open volumes. Thermal escape from “deep” localized states at vacancy defects has been observed in, e.g., Si and the experiments on the V-As₃ complex.⁴⁶ Further investigations will be necessary to elucidate this issue.

The CZ:Mg samples exhibit the largest overall anisotropy that is more triangular in the (S, W2) than in (S, W). We found that the positron lifetime is slightly longer in CZ:Mg than in EFG:Fe, and the low-to-high-temperature behavior clearly indicates that two distinct defects are contributing to the data. The anisotropy is consistent with that predicted for $V Ga ib$ split Ga vacancies, probably complexes with O vacancies and/or H impurities. Secondary ion mass spectrometry (SIMS) experiments have shown that hydrogen is much more efficiently incorporated into Mg-doped β-Ga₂O₃.⁴⁷ Therefore, it is possible that hydrogen indeed plays a role in causing the highest anisotropy in our Doppler measurements in the CZ:Mg samples. It is clear from the positron lifetime data that positron trapping is at saturation at temperatures up to 500 K. This indicates that the vacancy concentrations are at least—two to three times higher than in EFG:Fe. The temperature behavior strongly suggests that the split Ga vacancies are in the negatively charged state, similarly as in EFG:Fe.

In the EFG:Sn samples, the effects of the irradiation on the positron data are similar to EFG:Fe in the sense that the average positron lifetime increases, the overall anisotropy is reduced, and the (S, W2) data are shifted downward. This is consistent with unrelaxed V_Ga being introduced in the irradiation. Additionally, the shift of the (S, W) and (S, W2) data to the right suggests that V_Ga-nV_O complexes are involved, with n > 1. As the irradiation is predicted to create more V_Ga than V_O, the EFG:Sn crystals may contain pre-existing V_O that are invisible to positrons. This is also consistent with the as-grown samples containing split Ga vacancy complexes with V_O, as discussed in Sec. IV A.

Hardly any change can be seen in the average positron lifetime in the CZ:Mg samples after irradiation, but the slight increase at room temperature suggests that additional V_Ga defects with longer lifetimes are also found. As the concentration of the in-grown split Ga vacancies is significantly higher in CZ:Mg than in the EFG:Fe and EFG:Sn, this is as expected. The slight reduction in the overall anisotropy in (S, W) and (S, W2) is consistent with unrelaxed V_Ga being introduced in the irradiation, but the upward shift cannot be explained. Further investigations are necessary to elucidate this behavior.

The vast majority of the experimental literature focusing on the detection and identification of Ga vacancies in β-Ga₂O₃ agrees that in-grown Ga vacancies are found in their split configurations, also in complexes with H impurities and O vacancies. This is particularly so for experiments based on analyzing local vibrational modes with Fourier transform infrared reflectivity (FTIR)^13,20–23 and the details of positron annihilation.^12–16,26 Our present experiments strongly suggest that any of the three split configurations can be present in as-grown β-Ga₂O₃, although not necessarily all in the same samples. Theoretical calculations^9,45 predict that V^ic_Ga has the lowest formation energy, but the differences are of the order of 0.2–0.3 eV, small enough so that at the relatively high growth temperatures all are present at high concentrations. The details of the migration and stabilization by other defects and impurities during cooling down are likely to play a decisive role in determining which of the split Ga vacancies are present in the final product. We note that the fact that electron paramagnetic resonance (EPR) experiments have only been reported for irradiation or high-temperature annealing-induced defects in β-Ga₂O₃^{11,14,17–19} supports the interpretation that the in-grown Ga vacancies are mainly present as complexes rather than as isolated defects.

There is an ongoing debate on whether the EPR signals in β-Ga₂O₃ originate from unrelaxed or split Ga vacancies. Comparing experiments with state-of-the-art theoretical calculations suggests that the irradiation-induced Ga vacancies are in the split configuration.¹⁷ Our positron annihilation experiments indicate the opposite: both the positron lifetime and Doppler broadening show that unrelaxed Ga vacancies are introduced in the irradiation. We note that while the theoretical prediction for the formation energies of the unrelaxed Ga vacancies is 0.5–1.0 eV higher than for the split Ga vacancies,⁴⁵ irradiation should by definition create unrelaxed Ga vacancies. The barriers predicted for the relaxation into split Ga vacancy configurations are in the range of 0.3–1.4 eV depending on the initial and final configurations,⁴⁵ suggesting that the unrelaxed V_Ga may be stable at room temperature. In fact, in Ref. 26, it is suggested based on positron annihilation experiments that irradiation could introduce unrelaxed Ga vacancies that relax into split configurations in post-irradiation thermal treatments above 300 °C.

We have applied positron annihilation spectroscopy to study the vacancy defects in β-Ga₂O₃ crystals doped with Fe, Mg, and Sn before and after 6 MeV proton irradiation. By measuring the positron lifetime and the details of the anisotropy of the Doppler broadening signals as a function temperature and comparing them with theoretical calculations, we confirm that in-grown Ga vacancies are in the split Ga vacancy configuration irrespective of growth method or doping. We suggest that any of the three different split Ga vacancy configurations may be present, depending on the doping of the crystals and the growth method (EFG of CZ). We also report on the observation of unrelaxed Ga vacancies being introduced as the main defects in irradiation.

This material is based upon work supported by the Air Force Office of Scientific Research under Award No. FA8655-23-1-7057. This work has been partially supported by the Finnish Academy of Science and Letters through a personal grant (I.Z.) and the Finnish Cultural Foundation through the Additional Million-euro Funding to the Science program.

The authors have no conflicts to disclose.

Iuliia Zhelezova: Formal analysis (lead); Investigation (lead); Project administration (lead); Software (equal); Writing – original draft (equal). Ilja Makkonen: Software (supporting); Writing – review & editing (supporting). Filip Tuomisto: Conceptualization (lead); Data curation (lead); Formal analysis (equal); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.