A hydrothermal route with parameters optimized for larger nanocrystals was selected. The hydrothermal route (compared to the co-precipitation technique) provides better crystallinity and higher luminescence quantum yield (Runowski and Lis, 2016; Vanetsev et al., 2017). Rare earth (RE) oxides: La₂O₃, Ce₂O₃, Eu₂O₃, and Gd₂O₃ (99.99%, Stanford Materials, United States) were dissolved separately via dropwise addition of HNO₃ to the oxide powder. The resulting solutions were dried at 60°C, transferred to volumetric flasks, and solved in distilled water.

Ammonium fluoride, NH₄F (98+%, Sigma-Aldrich, Poland), was used as a source of fluoride ions. Citric acid trisodium salt dihydrate (Sigma-Aldrich, 97%, Poland) was used as received without further purification. Deionized water was used for the synthesis.

Appropriately calculated (for 0.5 g of the product) amounts of aqueous solutions of metal nitrate were mixed in a beaker. Water was added to give a total volume of 50 mL (Solution 1). In another beaker, 150% of the stoichiometric amount of NH₄F was dissolved by stirring in 50 mL of water at 60°C–70°C (Solution 2). After the complete dissolution of NH₄F, Solution 1 was added dropwise to Solution 2, under continuous stirring, without temperature control. After the addition, the mixture was left under stirring for another 15 min. The mixture was then allowed to precipitate, and the clear portion of the solution was decanted. The precipitate was purified via sedimentation in a centrifuge, decantation, and washing with distilled water: sedimentation for 3 min at 3,000 rpm, decanting, adding 35–40 mL of water, sedimentation for 6–7 min at 4,000–4500 RPM, decanting, addition of 35–40 mL water, and sedimentation for 12 min at 5,500 rpm. The purified precipitate was redispersed in approximately 50 mL of water and subjected to hydrothermal treatment in a Teflon vessel (20 h, 200°C). Then, the product was sedimented, the water was decanted, and the sediment was dried in an oven at 60°C–70°C.

X-ray diffraction (XRD) patterns were recorded using a Bruker AXS D8 Advance diffractometer in the Bragg–Brentano setup, with Cu K_α1 radiation (λ = 1.5406 Å), in the 6°–60° of 2Θ. Transmission electron microscopy (TEM) images were recorded using an FEI Tecnai G2 20 X- TWIN transmission electron microscope, which used an accelerating voltage of 200 kV. The photoluminescence decay curves and excitation spectra were measured under pulsed laser excitation (Opolette 355LD UVDM tunable pulsed laser (Opotek Inc.) operating at 20 pulses per second) QuantaMaster™ 40 spectrophotometer (Photon Technology International) with an R928 photomultiplier as a detector (from Hamamatsu). All spectra were measured at room temperature (293 K) and were appropriately corrected for the instrumental response. The excitation laser was characterized by a pulse length of 7–10 ns.

All of the fitting operations were performed using the custom code we developed. The code was written in the Python programming language, and the SciPy module (version 0.19.0) was used for the actual fitting routine. The code is available in its current form upon request. The code was used to read the input files, normalize the curves, downsample the data if necessary, and plot the fitted functions and the fit residual. Most importantly, the code provides a command-line-based user interface, where a guess for the fit can be typed in manually, copied and pasted, or called from the history of previous fits. In a typical procedure, the user loads a file and specifies the guess for the fit in the form of function name, list of parameters, another function name, list of parameters, and so on. The different functions (if there is more than one) are summed. It is possible to specify an ambiguously long set of any of the functions known to the code and to easily create new functions.

In the interface, fit results appear in the history together with guesses, meaning that the result of any of the fits can easily be used as a guess for a new fit; below, we shall call such an operation re-fitting. Re-fitting was used to test the results for stability. If a given set of fitted parameters is indeed optimal, then using it as a guess must, first, result in nearly identical parameters after fitting, and second, the fitting must converge after a small number of function evaluations (typically three to seven). Least squares fitting is a variational method of minimizing the difference between the experimental data and the fitted curve (i.e., minimizing the residual) and involves searching for a local minimum in an abstract space of fitting parameters. The fact that the result does not change significantly between cycles means that a given set of parameters corresponds to a certain stable minimum and is considered optimal. A small number of function evaluations means that the minimum is “deep” enough and is indeed a minimum. Another way to look at it is that a “correct” result must be repeatable; otherwise, it is not “correct.” The quotation marks are used because no fitted result is exact: all experimental data contains noise and uncertainties that can affect the fits.

Of the available least squares methods, the setup described below proved to be optimal. Other methods seemed more sensitive to the initial guess and more prone to sticking to local minima corresponding to obvious bad fits. In particular, least squares (rather than some more advanced minimization/optimization approach) was selected for historical and reverse-compatibility purposes: it is something people have used, are still using, and will likely use for the foreseeable future. Furthermore, least squares is the most basic approach and is readily available in commonly used software such as Origin, Matlab, Gnuplot, and Python. Finally, for the problem at hand, least squares gets the job done quite quickly.

Prior to the fitting, the decay curves were normalized by dividing by the peak value of the raw curve. The fitting was performed using the least_squares module of SciPy, with the Trust Region Reflective (“trf”) method (Branch et al., 1999; Byrd et al., 1988), 3-point Jacobian, lsmr solver (Fong and Saunders, 2011), and x_scale option. A description of the module can be found in SciPy (2011). In the least_squares module, the input is a vector of guess values of the fitting parameters; for example, x = [ t ₀ , I ₀ , A ₀ , τ ₁ , τ ₂ ]. The module calculates the gradient in each parameter, changing it by 1 × 10⁻⁸ and recalculating the residual function. Because the absolute initial values of the variables may differ by several orders of magnitude (e.g., τ _rise = 10 and A ₃ = 0.0002), they are changed in relatively different steps. However, if we define the x_scale vector equal to x, the problem is redefined in xs = x/x_scale. The guess values become coefficients, and the code manipulates the initial vector of ones, (Dirac, 1927), changing them (initially) by 1 × 10⁻⁸. Therefore, each parameter is changed proportionally to its value and has a similar effect on the cost function (Moré, 1978). The use of x_scale resulted in much greater stability of the fits, that is, faster convergence, faster change toward a “correct” fit, fewer incorrect fits, and smaller chance of getting “strange” (way too small or way too big) values of the fitted parameters.

Most often, the sum of the squared residuals (the residual being the fitted curve minus the experimental curve) was used as the loss function. However, in a few problematic cases, the Cauchy loss function was used.

Because there were 50,000 data points on each curve, and the functions included many exponential components, it would take a long time to fit completely from scratch. That is why downsampling was used: for each group of subsequent n points, the average values of both time and intensity were calculated, resulting in one data point instead of the original n. Noise was also effectively reduced in this procedure. The fits were performed at n = 30 (the original number of points was reduced by a factor of 30), then refits were performed at n = 20, 10, 5, and finally at n = 1, that is, without downsampling. In other words, the hardest part of trial-and-error search for optimal function was performed using smaller data sets, and the resulting fit was then corrected by re-fitting with a gradually increasing number of points.

One of the crucial indicators of fit consistency is the residual plot. Ideally, such a plot cannot contain any signal; only noise must remain (that is, points randomly oscillating up and down from zero intensity). Such residuals are referred to as “flat” because they are essentially (noisy) horizontal lines, with I = 0. In the case of incomplete fits, the residual plot contains some signal in addition to the noise; some pattern (e.g., rise, decay, rise, and decay, waves) is obvious (Lakowicz, 2006). Often, the pattern can be described as sine-like waves (variable-amplitude chirps, also called chirplets): a noisy residual curve going up and down with a gradually changing period.

Summarizing, for each curve considered, the fitting procedure and search for the optimal function was continued until a flat (signal-free, noise-only) residual was obtained. The fitting started with downsampled curves; once a suitable solution was found, refits were performed with gradually smaller downsampling and, finally, for the entire dataset. At each step, subsequent refits were performed several times until it was clear that the solution was stable; that is, the result very similar to the guess is obtained after a small number of function evaluations, namely, three to seven.

The functions used to fit the curve were linear combinations of the pulse functions of the kind I = A (1 − exp (−t/τ ₁)) exp (−t/τ ₂). In order to justify the use of such functions, several simple model systems were proposed and described using differential equations. The sets were solved analytically, and the solutions were pulse functions. We have obtained the solutions using the WolframAlpha mathematical engine (Wolfram Alpha LLC, 2018).

In the examples, we considered processes with monoexponential rise and monoexponential decay. There may be more complicated cases where there are several emitting sites characterized by different decay rates. Considering the same rise lifetimes, the total observed emission intensity can be described by Equation 14, which is basically an exponential rise multiplied by the sum of two exponential decays (a biexponential decay): I = A 0   1   –   exp   –   t   /   τ rise     ⁡ exp –   t   /   τ 1   decay + A 2 ⁡ exp –   t   /   τ 2   decay . (14)

Note the missing A ₁ coefficient before the τ _{1 decay} exponent (which is actually held constant and equal to unity). It can be shown that A ₁ is absorbed by the A ₀ coefficient and is thus redundant from the fitting point of view; freezing A ₂ = 1 results in essentially the same fit result with one less degree of freedom, thus improving the fit stability.

Similarly, we can define a two-rise-one-decay kinetics, which can be interpreted as two sites that differ in population rates but have the same decay rates: I = A 0   1   –   exp   –   t   /   τ 1   rise + A 2   1   –   exp –   t   /   τ 2   rise     exp –   t   /   τ decay . (15)

Again, one of the coefficients is redundant.

Functions like Equations 14, 15 have been successfully used, for example, by Bartkowiak et al. (2019) and Runowski et al. (2018).

Things become more complicated if we consider two (or more) rise components and two (or more) decay components. Consider Equation 16 and compare it with Equations 17, 18. Although Equation 16 is a multi-rise-multi-decay function, the other two are sums of one-rise-one-decay functions. I = A 0   1   –   exp   –   t / τ 1   rise + A 2   1   – ⁡ exp –   t / τ 2 rise       ⁡ exp –   t / τ 3   decay + A 4 ⁡ exp –   t / τ 4   decay . (16) I = A 0   1   –   exp   –   t / τ 1   rise   exp –   t / τ 3   decay + A 1   1   –   exp   –   t / τ 2   rise   exp –   t / τ 4   decay . (17) I = A 0   1   –   exp   –   t / τ 1   rise   exp –   t / τ 4   decay + A 1   1   –   exp   –   t / τ 2   rise   exp –   t / τ 3   decay . (18)

While Equations 17, 18 look essentially similar, they represent two ambiguous solutions. Fitting the same curve using Equations 16, 17 or Equation 18 gives approximately the same τ _1–4, R ², and residual pattern. However, there are two different ways to distribute the rise-and-decay lifetimes between the two one-rise-one-decay functions, namely, (rise-1 times decay-3; rise-2 times decay-4) or (rise-1 times decay-4; rise-2 times decay-3). They may correspond to a site with τ ₁ rise lifetime, τ ₃ decay lifetime, and another site with τ ₂ rise lifetime, τ ₄ decay lifetime. Alternatively, there might be a site with τ ₁ rise lifetime, τ ₄ decay lifetime, and the other with τ ₂ rise lifetime, τ ₃ decay lifetime. Equation 16 is apparently less ambiguous, stating that there are two rise and two decay components in the fitted kinetics.

Opening the brackets shows clearly that Equation 16 is the sum of four pulse functions, where the interdependence of the parameters is obvious: all of the amplitudes affect more than one of the exponent components (Equation 19). Consequently, none of the components can reach 0 and not affect the others at the same time. Such functions, although they result in stable fits, can be very difficult to interpret. Functions of this kind will be referred to as convoluted below. I = A 0   1   –   exp   –   t / τ 1   rise   exp –   t / τ 3   decay + A 0   A 4   1   –   exp   –   t / τ 1   rise   exp –   t / τ 4   decay + A 0   A 2   1   –   exp   –   t / τ 2   rise   exp –   t / τ 3   decay + A 0   A 2   A 4   1   –   exp   –   t / τ 2   rise   exp –   t / τ 4   decay . (19)

Alternatively, consider functions of the following kind: I = A 13   1   –   exp   –   t / τ 1   rise   exp –   t / τ 3   decay + A 14   1   –   exp   –   t / τ 1   rise   exp –   t / τ 4   decay + A 23   1   –   exp   –   t / τ 2   rise   exp –   t / τ 3   decay + A 24   1   –   exp   –   t / τ 2   rise   exp –   t / τ 4   decay . (20)

As in Equation 16, in Equation 20, there are four lifetimes—two rises and two decays. Amplitudes are now independent. Functions of this kind will be referred to as non-convoluted below.

Equation 20 has eight independent parameters instead of seven in Equations 16, 19; however, it is much more flexible and clear. In Equation 20 (taking into account the special cases where some amplitudes are close to 0), it is possible to distinguish between the cases of Equations 17, 18. The presence of one-rise-two-decay or two-rise-one-decay sites can be identified using amplitude similarities. However, fitting the same experimental curve with four independent one-rise-one-decay functions (12 parameters in total) may prove to be overfitting.

Consequently, in this study, we used two kinds of functions: the type of functions from Equation 20 and the type of functions from Equation 16 (non-convoluted and convoluted, respectively). The functions are labeled by the number of rise-and-decay components. For instance, the pulse function is r1d1 (one-rise-one-decay), Equation 14 is r1d2 (one-rise-two-decays), Equation 15 is r2d1 (two-rise-one-decay), and so on. A full list is provided in Supplementary Data Sheet S2, Section SI.1.3.

From a practical point of view, it is convenient to start with the functions of Equation 16 kind (the convoluted ones, fewer parameters), find a stable solution, and express the solution in the form of the corresponding function of Equation 20 (the non-convoluted kind), and refit. As the functions of the convoluted kind have fewer degrees of freedom, in the present study, it was usually more difficult to find a solution using the functions of the non-convoluted kind from scratch—hence the two-step approach.

Systems of rate equations were used to build a function (subroutine) that calculated the dn _i/dt derivatives based on the current values of n _i (at a given time) and the considered process rates. Taking into account the derivative subroutine and rates, the equations were integrated using the scipy.integrate.odeint module from SciPy (SciPy, 2011) with default settings [which is a wrapper of the ODEPACK (Hindmarsh, 1983) library], and call the Livermore Solver for Ordinary Differential Equations (LSODE) Adams/Backward Differentiation Formula method with automatic stiffness detection and switching (Hindmarsh, 1983; Petzold, 1983). As a result of integration, time evolution curves were obtained for all levels of the model. The selected curves could thus be compared with their respective experimental counterparts. One can, therefore, treat the system of equations variationally, manipulating the values of the rates in order to minimize the residual between the simulated and the experimental curves. In such a case, integration is performed from scratch at each fitting step with appropriate parameters. The time axis is defined as an array of t values with a specified step. In the case where rates were found variationally, the time axis was the one from the experimental sample, extrapolated to t = 0 with the same step as in the experimental curve.

According to X-ray diffraction data, the obtained materials corresponded to hexagonal P 6 ₃ /m m c (space group nr. 194) LaF₃. In the material, cell angles α and β are 90°, while γ is 120°. The unit cell comprises two fluoride formula units, with one kind of La position and two kinds of F positions. The La atoms are surrounded by 11 F atoms in a D_3h local symmetry. The main C₃ axis coincides with the cell Cartesian z-axis. Two fluoride atoms occupy the axial positions, three form an equilateral triangle in the equatorial plane, and the other six form two equilateral triangles in the planes parallel to the equatorial one. These triangles are rotated 30° with respect to the equatorial triangle around the C₃ axis. The site has no inversion symmetry. As this geometry is not further discussed in this paper, the visualization is provided in Supplementary Data Sheet S2, Figure S3. Substitution of the dopant ions in the La sites is assumed.

Scherrer analysis estimates the average particle size to be approximately 50 nm. The XRD plots are shown in Supplementary Data Sheet S2, Figure S4. According to TEM, the materials are composed of well-defined nanocrystals. The widths start at 20–30 nm and span up to 50 nm. The lengths also start at 20–30 nm and span up to 80–100 nm. The material is thus a mixture of somewhat spherical particles and variable-length thick rods. Among other particles, rhombic-shaped nanocrystals (characteristic of hexagonal structures) are detectable. The TEM images are shown in Supplementary Data Sheet S2, Figures S5, S6.

Particle sizes allow for the estimation of surface site percentage. Depending on the crystal surface, two layers of metal sites correspond to 5–7 Å layer of atoms. For a spherical nanoparticle of 50 nm diameter, the outer layer of 0.6 nm (considered as the surface) would correspond to approximately 7% of the nanoparticle volume, which is considerable. A similar surface layer size was used by Vanetsev et al. (2017). On the other hand, as the metal nitrate solution was added dropwise to the ammonium fluoride solution, the initial precipitation occurred in the effective excess of the fluoride. It is thus likely that the surface metal cations are fully caped with fluoride anions and not with OH groups. The lack of capping agents and drying of the samples should have reduced the amount of OH groups and surface water as well. It is likely that surface quenching should not be strong in the materials. In addition, due to the different crystallization rates of different lanthanide fluorides (Gulina et al., 2017), the dopants might have a tendency to be located closer to the center of the particle (Gulina et al., 2017; Grzyb et al., 2013).

In the case under study, there were basically two types of curves—decay curves and rise-and-decay curves. The decay curves are the sum of the following components: A _i exp (− t/τ _i). Given the curve that is clearly decaying (that is, no apparent rise is present), the first guess was a monoexponential decay. If the result was unsatisfying (the residual was not “flat”), a second decay was added, then a third, and so on, until the residual contained no apparent signal. However, a curve that looks like a decay is actually the late part of a rise-and-decay. Fitting such a curve with a monoexponential decay will be unsatisfactory. With more components, any “hidden” (not obvious in shape) rising part of the curve would be visible as an exponential component with negative amplitude.

Given the curve with an apparent rise, the first guess was r1d1. The fitting process was then continued until a flat residual was obtained by adding more rise and/or decay components. In such a process, the experimenter’s intuition guides the process. It can be said that if the residual goes up at the beginning, it is probable that another decay component will be required, while a rise component will be required if the opposite is true. However, this was not the rule.

In many cases, it was impossible to correctly guess the entire curve. In such cases, the curve was cropped at some point. For instance, only the first 2 ms, 5 ms, or 10 ms were fitted with simple functions such as r1d1 or r2d1. In a separate fit, the last 10 ms, 20 ms, and 30 ms were fitted with mono- or biexponential decay. This step-by-step approach allowed for the correct detection of rise-and-decay lifetimes, as well as the number of the respective components. Then, taking into account the obtained values, a guess was constructed for the entire curve. It is essential that the resulting whole-curve lifetimes match (at least by the order of magnitude) the lifetimes of the partial fits. If the late part of the curve has monoexponential decay with a lifetime of, for example, 10 ms, then the entire curve fitted with a rise-and-decay function must exhibit a decay lifetime of approximately ∼10 ± 1–2 ms. A decay lifetime of ∼2–3 or 40–1,000 ms should raise suspicions.

It turned out that the decay components of the many-rise-many-decay functions cause the most problems. In other words, the rising part of the curves was relatively easy to guess and fit, while the decay parts were unconvincing. Thus, all of the fittings began with partial fits of the latter parts of the curves, where, in most cases, no rise was present—that is, where the curves appeared as decays and could be fitted with one to three exponential functions with positive amplitudes.

Finally, for some whole-curve fits with rise-and-decay functions, some of the resulting amplitudes were negative. Such a situation indicates an improperly selected type of component (lifetime) with a negative amplitude: the respective rise component (lifetime) must be, in fact, a decay component (lifetime) or vice versa. For example, an r2d3 fit exhibiting a negative amplitude for one of the decay components indicates that r3d2 must be used instead.

It also happened that some of the rise lifetimes approached some of the decay lifetimes. Both effectively cancel out, which in practice means that either the respective rise is redundant and should be removed; (for example, r2d2 should be used instead of r3d2), or both the respective rise and decay are redundant (e.g., r2d1 must be used instead of r3d2). Another indication of this situation is a rapid increase (especially between the consecutive refits) of the selected amplitudes, sometimes by several orders of magnitude, along with a synchronous change in some rise lifetime and some decay lifetime, both having similar and rather erroneous values (e.g., milliseconds or seconds, when microseconds were expected).

Sometimes, finding a good fit meant looking for trends and patterns. Similar systems must result in similar kinetics. Particularly, the Eu³⁺ emission at a given wavelength must exhibit similar rise-and-decay lifetimes (at least of the same order of magnitude), as well as the same (or similar) numbers of rise-and-decay components. In most cases, such correspondence occurred naturally: a series of good fits post factum turned out to exhibit trends and similarities in the values of the variables. However, in the selected cases, there were several options for good fits, while only one had to be selected as “correct”—the one most similar to the other fits of the profiles of emission at the wavelength in question.

One pulsed-excitation measurement was performed for the LaF₃:Gd³⁺ sample with λ_ex. = 272 nm and λ_em. = 312 nm. The curve was a rise-and-decay type. After the excitation pulse, some of the Gd³⁺ ions are in the ⁶I_{17/2, 15/2, 13/2, 11/2, 9/2,7/2} excited states, which for the convenience of readers will be abbreviated and referred to here as the ⁶I manifold/state/level. Radiative decay to the ⁸S_7/2 ground state or a non-radiative process to the ⁶P_7/2,5/2,3/2 levels (⁶P manifold) may occur. Because the energy gap in the latter process is neither large nor small (∼4,000 cm⁻¹), it is difficult to say which process will prevail; however, it is clear that the latter process occurs, resulting in the ⁶P manifold population and the ⁶P → ⁸S_7/2 emission at 312 nm.

Given the populating process: Gd³⁺: ⁶I → ⁶P, the kinetics of the ⁶P manifold must follow a rise-and-decay pattern, as discussed in Section 2.3.2: I = I 0 + A 0   1   – ⁡ exp –   t   –   t 0   /   τ rise exp –   t   –   t 0   /   τ decay . (21)

Here, τ _rise = 1/(W _rad. + W _nrad.), τ_decay = 1/W _rad., W_rad. is the rate of the ⁶P → ⁶S_7/2 radiative relaxation, and W _nrad. is the rate of the ⁶I → ⁶P non-radiative process. The vertical (intensity) offset I ₀ and horizontal (time) offset t ₀ are also featured, giving the equation actually used in the fitting.

The r1d1 function (Equation 21, solution S1; Table 1) can be fit to the experimental curve, obtaining a rise time of 363 μs and a decay time of 12,839 μs. However, the result (Figures 2a, c) indicates that a more complex function should be used. Note the wavy residual in Figure 2c. Such a residual clearly contains some signal (significant within 800–5,000 μs after the pulse), indicating an incomplete (bad) fit. A good fit can be unambiguously obtained using the following equation (r2d1 function), with the rise lifetimes of 98.7 µs and 700 µs (the corresponding amplitude ratio was ∼1:0.235, contributions of 81% and 19%, respectively) and the 12,700 μs decay (solution S2, Table 1): I = I 0 + A 0 1   – ⁡ exp –   t   –   t 0   /   τ rise 1 + A 2     1   – ⁡ exp –   t   –   t 0   /   τ rise 2 exp –   t   –   t 0   /   τ decay . (22)

The corresponding residual in Figure 2d contains no noticeable signal—at least with amplitudes significantly larger than the noise. It can also be seen in Figure 2b that Equation 22 sits on the experimental points much better. Thus, there are two exponential components in the rise in Figure 2. One of them must be the Gd^{3+ 6}I → ⁶P relaxation, while the nature of the second is discussed below. All of the fits discussed below resulted in the overall picture shown in Figures 2a, c or Figures 2b, d and were considered as bad or good fits, respectively.

There is another, more complex way to fit a given kinetics. Namely, two independent r1d1 equations (Equation 21) can be used. In this case, the I ₀ and A ₀ variables are replaced by I ₁, A ₁, and I₂, A ₂ in the first and second functions, respectively. The S3 solution (Table 1) was obtained from the S1 solution by adding another r1d1 component. The associated solution S4 was obtained from S3, using the τ_decay2 < τ_decay1 in the initial guess. The fit of the S4 is only as good as the S3. Finally, the solution S5 guess was obtained from solution S3 by swapping the decay lifetimes in the two components and setting both amplitudes to 0.5. Once again, it was a good fit of the Figures 2b, d kind. The resulting decay lifetimes are very close and consistent with the S2 solution. The decay lifetimes of solutions S3 and S5 differ from the decay lifetimes of solutions S1 and S2.

In solution S4, τ _decay1 is very close to τ _decay1 in S1, which means that τ _decay2 of solution S4 is redundant. Thus, this solution is incorrect. This was shown as an example of an error. Another bad solution is S5: its similar values of τ _decay indicate that only one decay component had to be used in the summary function (i.e., r2d1 had to be used), which is consistent with solution S2.

Solution S3 provides a different interpretation than S2. In S2, the two species of Gd³⁺ are considered to interact with each other, as described in the next section. Two independent equations are used in S3, both corresponding to an isolated Gd³⁺ site. They can be considered a surface site with faster ⁶I → ⁶P relaxation due to the presence of OH groups and a shorter radiative lifetime originating from the less symmetric coordination geometry (τ_rise = 99.3 μs, τ_decay = 12,466 μs). The second site, with a longer radiative lifetime and weaker non-radiative quenching (τ_rise = 703 μs, τ_decay = 13,748 μs), must be the bulk site. However, this interpretation implies that the ratio of surface to bulk sites is 0.914:0.214 ≈ 4.27:1. In other words, 81% of emissions come from the surface sites, which is possible albeit unexpected. However, considering the S4 and S5 solutions, the S3 solution is likely a redundant solution and not a true indication of the dominant surface sites.

It is worth noting that the S2 and S3 fits are equally good. The S2 solution contains six fit parameters, while S3 contains seven fit parameters. From the point of view of numerical complexity, S2 is preferred. However, the two solutions assume different physics of the process. In S2, there is only one decay lifetime, while in S3, there are two. As mentioned in Section 3.2, of two equally good solutions, the one with the smaller number of decay components should be selected.

This section illustrates the care that should be taken when fitting multiexponential functions. The functions are guess-sensitive and (to some extent) flexible, providing opportunities for mistakes and misinterpretations (including intentional ones). Care should be taken during fitting routines, during interpretation, and when describing the result. In particular, as we illustrate here, it is a good and recommended practice to show the residuals (Debasu et al., 2011) and describe the fitting peculiarities.

The sample, excited in either the 272 nm f-f band of Gd³⁺ or the 250 nm f-d band of Ce³⁺, shows a sharp emission peak at 312 nm. Much weaker, broad emission bands are observed in the range of 290–400 nm. The bands coincide with the Ce³⁺ d-f emission in LaF₃:Ce³⁺,Eu³⁺ (Figure 4). Consequently, in the case of 250 nm excitation, not all the energy from Ce³⁺ is transferred to Gd³⁺. Some of the energy is emitted as light, which is not surprising given the allowed nature of the d-f transition of Ce³⁺ and its very short lifetime. In the case of 272 nm emission, the broadband emission indicates that some portion of Ce³⁺ ions become excited, either by direct absorption of the excitation light or by energy transfer.

In contrast to LaF₃:Gd³⁺, in LaF₃:Ce³⁺,Gd³⁺ (excited both in the 272 nm f-f band of Gd³⁺ or in the 250 nm f-d band of Ce³⁺, λ_em. = 312), there are several decay components, while the rise is completely missing. Both decay profiles can be fitted with three or four exponential components; see Table 2.

The fitting results are quite sensitive to the fitting procedure. In particular, there are several options for treating the time axis offset. Recorded curves do not start at zero time; the first point has a time value of 831 μs. One of the options is to not use time offset whatsoever, simply because each piece of the exponential decay profile is the same decay profile, except for amplitude. However, this logic applies under ideal conditions, where noise is exactly 0. Because experimental noise is present, a short-lived component with a lifetime of, say, approximately 50 μs would be practically unnoticeable at 831 μs and would require a ridiculously large amplitude to have any effect. Another option is a variable time offset. However, this option introduces an undesirable degree of freedom, which also renders lifetimes and amplitudes dependent on its value, and results in a form of internal dependence between parameters. Therefore, we used a fixed offset of 831 μs; that is, we assumed that the experimental decay profiles start at the time point of 831 μs. Such fits are shown in bold in Table 2 and are considered “correct.” Note the differences with the other options. This approach also eliminates ambiguity in the number of components. With variable offset, as well as without offset, it is possible to achieve an acceptable fit with only three components. With a fixed offset of 831 μs, four components are clearly visible. We emphasize that all the fits were good, that is, stable, reproducible (although guess-dependent), and resulting in the flat residual.

In the LaF₃: 1%Eu³⁺, 1%Ce³⁺ sample obtained in the same series of samples as the discussed LaF₃: 1%Gd³⁺, 1%Ce³⁺ sample, a broad band of the Ce³⁺ f-d emission is observed in the 280–320 nm range, peaking at 303 nm. Consequently, a strong overlap of the Ce³⁺ f-d excitation band with the 272 nm ⁶I ↔ ⁸S_7/2 band of Gd³⁺ is expected. The Ce³⁺ f-d emission overlaps with the 312 nm ⁶P ↔ ⁸S_7/2 band of Gd³⁺. Thus, the following energy transfer is possible: Gd^{3+ 6}I → Ce³⁺ f-d → Gd^{3+ 6}P. This mechanism is supported by the fact that the 272/312 nm decay profiles of Gd³⁺ in LaF₃:Gd³⁺,Ce³⁺ do not show any rise, indicating the involvement of Ce³⁺ in the Gd^{3+ 6}P dynamics. Energy transfer processes with Ce³⁺ as one of the parts must be very fast due to the allowed f-d transitions of the latter.

The Ce³⁺ f-d radiative decay is not forbidden and is characterized by a very short lifetime [we used the value of 29.2 ns (Kroon et al., 2014)], while the Gd³⁺ radiative lifetime is much longer, approximately 12.5–12.6 ms. In other words, the decay rate of the Ce³⁺ f-d transition is approximately half a million times greater than the radiative decay rate of Gd³⁺. Roughly the same ratio is expected for the respective absorption probabilities. Consequently, when LaF₃:Ce³⁺,Gd³⁺ is excited at 272 nm, mostly Ce³⁺ ions are excited. At 250 nm excitation, only Ce³⁺ ions are excited. Thus, both cases must be characterized by very similar overall kinetics despite the fact that particular excited atoms and groups of atoms may depend on the excitation wavelength. It is worth noting that the excited spot of the sample also varied slightly, depending on the wavelength. Thus, “correct” decay fits satisfy both quality and consistency requirements; that is, the two excitation cases exhibit similar lifetimes and similar amplitudes.

In the LaF₃:Eu³⁺; LaF₃:Ce³⁺,Eu³⁺; LaF₃:Gd³⁺,Eu³⁺; and LaF₃:Ce³⁺,Gd³⁺,Eu³⁺ samples, the Eu³⁺ dopant exhibits similar, albeit excitation-dependent, emission spectra. In particular, the 613 nm band (hypersensitive), the 618 nm band, and, to some extent, the 583 nm band exhibit different relative intensities when excited at either 250 nm or 272 nm excitation (Figure 6). In the case of 395 nm excitation, the variation is much smaller. In other words, emission is also dependent on the sensitization path. Noteworthy, the bands in the 500–550 nm range (⁵D₁→⁷F₀ at 534–536 nm, ⁵D₀→⁷F₁ at 525.7 nm, and ⁵D₂→7F₃ at 509–510 nm) are much less distinct if a Ce co-dopant is present (Figure 6c).

This property is shown in Figure 6d, which shows the integrated areas of the emission spectra ranges 586–606 nm and 606–640 nm, in % of the total integrated spectrum. The very presence of such a phenomenon indicates some inhomogeneities in the structure. At least two spectroscopically different (although similar) Eu³⁺ sites are present. The sites differ not only in the first coordination sphere (judging by the shapes of the spectra) but also in their Ln³⁺ nearest neighbors (judging by the dependence on the sensitization path).

It is obvious that the excitation dependence is most prominent in samples containing Ce³⁺. Because the ionic radius of Ce³⁺ (1.196 Å, c. n. 9) is similar to that of La³⁺ (1.216 Å, c. n. 9), doping of LaF₃ with Ce³⁺ should not result in any significant defect formation. On the other hand, the core-shell structure of TbF₃ and CeF₃ (instead of the mixed-lanthanide system) may form spontaneously (Grzyb et al., 2013). It is thus plausible that the introduction of Ce³⁺ somehow changes the character of the formed dopant ion clusters, resulting in increased asymmetry of the coordination geometry of some Eu³⁺ ions.

This paper presents four kinds of Eu-containing samples and six wavelengths at which the emission temporal evolution profiles were measured. There are thus 24 rise-and-decay profiles, and each of them can tell a different story. While some fits were unambiguous and simple, others required a lot of work to complete. Only the most noteworthy and illustrative cases of fitting peculiarities are discussed below, with plots where necessary. For clarity, in many cases, only lifetimes are shown below. The respective full data for the mentioned cases are provided in the supplementary tables (Supplementary Data Sheet S1).

The profiles were fitted with the function of the following kind: I = I 0 + A 11   1   –   exp   –   t / τ 1 r exp –   t / τ 1 d + A 12   1   –   exp   –   t / τ 1 r exp –   t / τ 2 d + A 21   1   –   exp   –   t / τ 2 r exp –   t / τ 1 d + A 22   1   –   exp   –   t / τ 2 r exp –   t / τ 2 d + A 31   1   –   exp   –   t / τ 3 r exp –   t / τ 1 d + A 32   1   –   exp   –   t / τ 3 r exp –   t / τ 2 d + A 13   1   –   exp   –   t / τ 1 r exp –   t / τ 3 d + A 23   1   –   exp   –   t / τ 2 r exp –   t / τ 3 d + A 33   1   –   exp   –   t / τ 3 r exp –   t / τ 3 d . (39)

This function was used for all samples containing the Eu³⁺ dopant. It thus contains a total of nine rise-and-decay components (r3d3), as that was the highest complexity required. For many samples, only some of the components were used. The general form of the function remained the same to facilitate comparison. See Supplementary Data Sheet S2, Section SI.1.3 for the full list of functions. In several cases, a function with one rise and four decay components was used (r1d4, Equations 40, A ₅ = 0): I = I 0 + A 0   1   –   exp   –   t / τ r   ⁡ exp –   t / τ 1 d + A 2 ⁡ exp –   t / τ 2 d + A 3 ⁡ exp –   t / τ 3 d + + A 4 ⁡ exp –   t / τ 4 d + A 5 ⁡ exp –   t / τ 5 d   ) . (40)

The functions are shown in a slightly simplified form. The actual functions used also had time offset as a variational parameter, that is, (t−t ₀), instead of only t in Equations 39, 40. Equation 40 has five decay components and is thus the r1d5 function. With A ₂-A ₅ set to 0, the simplest pulse function (r1d1) is obtained. With A ₃-A ₅ set to zero, the function is r1d2. With A₃–A₅ set to 0, the function is r1d2. With A ₄ and A ₅ set to zero, the function is r1d3.

Exponential decay has a general form: I = I 0 + Σ i     A i ⁡ exp –   t / τ i   , i = 1 , 2 , 3 , . (41)

The fitting results for the temporal evolution profiles of the LaF₃:Eu³⁺ emission are shown in Tables 12, 13. To improve the readability of the tables, the coefficients have been rounded to three decimals and the lifetimes to at least three digits. Tables with more significant digits are shown in the supplementary tables (Supplementary Data Sheet S1).

Transition manifolds were identified using the Carnall tables. The 618 nm curve has noticeably different lifetimes than other transitions from the ⁵D₀ manifold. According to the tables, for the ⁵D₀→⁷F₂, transition, the wavelength is 614.7 nm, while the wavelength of the ⁵D₁→⁷F₄ transition is 617.1 nm. The corresponding A ₃₂ coefficient of the 618 nm fit (Table 12) is much lower than for the 591 nm, 613 nm, and 694 nm traces, which indicates a small participation of the rise-3-decay-2 component (5.5 ms rise time, 7.5 ms decay time). Therefore, it was concluded that there is an admixture of ⁵D₁→⁷F₄ emission in the dynamics of 618 nm photoluminescence. The resulting decay-2 component falls off the trend with its value of approximately 7.5 ms, while the remaining three are approximately 1 ms.

The 552 nm emission unambiguously corresponds to ⁵D₁→⁷F₂ of Eu³⁺. The 583 nm emission can be assigned to ⁵D₀→⁷F₀. However, Carnall’s tables show that the latter must be located at 578.3 nm. The energy of the ⁵D₁→⁷F₃ transition corresponds to 582.7 nm, while the dynamics of the 552 and 583 nm emission are quite similar, except for the decay-3 component of the latter. The coefficients involving decay-3 are rather small. Thus, the 583 nm emission is mainly ⁵D₁→⁷F₃, with some addition of ⁵D₀ radiative decay. This decay may be due to the optical overlap of ⁵D₀→⁷F₀ emission from an independent site but may also originate from energy transfer interactions of a pair of neighboring Eu³⁺ ions.

One can speculate that the average lifetime of ⁵D₀ is approximately 13 ms, while the average lifetime of ⁵D₁ is approximately 1–5 ms. However, due to the complexity of the observed rise-and-decay kinetics and the fact that dopant-dopant interactions were observed in the LaF₃:Gd³⁺ and LaF₃:Ce³⁺,Gd³⁺ samples, there must be several interacting Eu³⁺ species in the studied system. The observed lifetimes may, therefore, result from the system of the interacting levels as a whole and not correspond to the actual radiative rates (W _rad. ) or total decay rates (W _rad. + W _nrad. ) at the levels.

The Ce and Gd samples mentioned above could be treated using relatively simple model systems and empirical/fittable parameters. In contrast, the Eu³⁺ dopant comprised dozens of manifolds below the used excitation energies (corresponding to 272 nm for Gd-doped samples and 250 nm for Ce-doped samples) (Peijzel et al., 2005). Even a system of a single Eu³⁺ ion is a challenge, while a model system with two Eu³⁺ ions would contain many more degrees of freedom (approximately half a hundred) than can be reliably inferred from the experimental kinetics in question (in particular, six amplitudes and five lifetimes). To properly describe such a system, it is necessary to find transition dipoles for the transitions in question, calculate energy transfer rates, and estimate non-radiative transfer rates, as, for example, in Shyichuk et al. (2016b). In this paper, we thus leave the Eu³⁺ data without deeper analysis and keep the focus on temporal profile fitting.