PyMossFit: A Google Colab Option for Mössbauer Spectra Fitting

1. Introduction

Mössbauer effect spectroscopy [1,2] is a highly specialized technique that investigates the resonant absorption of $\gamma$ rays by atomic nuclei. The Mössbauer effect, observed primarily in isotopes such as iron-57 (⁵⁷$Fe$), provides information on hyperfine interactions, including isomer shifts, quadrupole splittings, and magnetic hyperfine fields. These parameters offer detailed information on the electronic, magnetic and structural environment of the sample, making the technique invaluable in materials science, condensed matter physics, chemistry, biology [3], planetary science [4], and also in cosmology [5]. The gamma ray of resonant absorption in ⁵⁷$Fe$ nuclei is 14.4 keV and hyperfine interactions promote the break of degeneracy at ground and excited nuclear levels. The above cited parameters are described in the following subsections:

1.1. Isomer Shift, $IS$ or $\delta$

This phenomenon corresponds to a certain probability that electrons have to be inside the nucleus [6], an exclusive quantum event. The next equation describes this parameter.

$$IS=C·{\displaystyle \frac{\delta R}{R}}·\left(\right|{\psi }_A{\left(0\right)|}^2-|{\psi }_S\left(0\right){|}^2)$$

(1)

where ${\psi }_A\left(0\right)$, ${\psi }_B\left(0\right)$ are the s-wavefunctions at the origin of the absorbent (A) and the $\gamma$-ray source (S) nucleus, while ${\displaystyle \frac{\delta R}{R}}$ is the relative nuclear radius R change during recoilless $\gamma$ resonant absorption.

1.2. Quadrupole Splitting, $QS$ or $\mathsf{\Delta }$

This quantum effect gives account of the interaction of nuclear charge with an electric field gradient (EFG) originated in an asymmetric electronic orbital distribution. The analytical expression is

$$QS={\displaystyle \frac{1}{2}}e·Q·V_{ZZ}\sqrt{(1+{\displaystyle \frac{1}{3}}·\eta )}$$

(2)

being Q the quadrupolar moment, $V_{ZZ}$ the principal component of the EFG, and $\eta$ the asymmetric factor.

1.3. Magnetic Hyperfine Field, $B_{HF}$

Nucleons have an intrinsic quantum property, the spin. When the spin is coupled with the quantum orbital moment, a total moment can interact with the local magnetic field induced by the electronic orbitals, giving rise to a more complex splitting (Zeeman effect) and, according to the nuclear selection rules with some forbidden transitions, and consequently, certain decays by $\gamma$-ray emission are not physically possible.

The Figure 1 shows the energy levels schemes of all of these hyperfine interactions for the ⁵⁷$Fe$ isotope. As can be observed, IS shifts the excited and ground energy levels between the source nucleus and the sample one; meanwhile, QS promotes the splitting of the excited level and $B_{HF}$ the same in both the ground and excited levels, with the $I_e=+3/2\to I_g=-1/2$ and $I_e=-3/2\to I_g=+1/2$ forbidden transitions by the selection rules.

In a typical Mössbauer experiment, a radioactive source emits gamma rays, which are absorbed by the nuclei in the sample (see Figure 2). Resonant absorption occurs, in a small fraction, when $\gamma$ rays hit the probe with recoilless. The detector measures the intensity of the transmitted radiation as a function of the velocity of the gamma-ray source (with an additional energy added by Doppler effect). This results in a Mössbauer spectrum, typically characterized by sharp peaks or dips at resonance frequencies corresponding to different hyperfine interactions within the sample. The most common experimental setup corresponds to Transmission Geometry (in this case, the observed lines come from a reduction of 14.4 keV gamma counting in the detector, as compared to the background signal). Another possible option is the Backscattering Mössbauer Spectroscopy (BMS) which corresponds to detection of back-scattered electrons after gamma absorption. In recent years, the use of synchrotron radiation was extended, helping to collect Mössbauer spectra with deeper details and in a short time [7].

1.4. Statement of Need

Recently, in a review by Grandjean et al. [8], the authors made a series of suggestions about how good measurements should be taken, which could be a good practice for Mössbauer data treatment and its corresponding fitting presentation. Usually, scientists working with Mössbauer spectroscopy manage their data and fit them with non-open source software that runs on a Windows OS with poor upgrade services [9]. The needs of an open source option with less OS or package dependency has motivated this work.

In this sense, Google Colab is a useful tool for a collaborative team job in data analysis with the advantage of no additional packages locally installed, also compatible with the fact that the user can run their codes from multiple devices.

2. Mössbauer Spectra and Curve Shapes

The shape of a Mössbauer spectrum varies depending on the nature of the hyperfine interactions in the sample. For example, a simple paramagnetic material might produce a single absorption peak (Lorentzian or Gaussian) due to the isomer shift. Materials with quadrupole splitting generate a doublet (two peaks), while materials experiencing magnetic hyperfine interactions often exhibit more complex sextet patterns. These spectral features often overlap, making it difficult to isolate and quantify each individual component.

2.1. Curve Fitting and Least Squares Method

Curve fitting plays a crucial role in Mössbauer spectral analysis, as it allows for the decomposition of these complex spectra into individual contributions, each associated with specific hyperfine parameters. The challenge lies in accurately reproducing the spectral shapes using mathematical models and adjusting the parameters until a satisfactory fit is achieved.

To accurately extract physical parameters from Mössbauer spectra, fitting procedures are applied to the experimental data. One of the most common approaches is the least squares method, which minimizes the difference between the experimental data points and the theoretical model curve. The objective is to adjust the parameters of the theoretical model (such as isomer shift, quadrupole splitting, and line broadening) so that the calculated spectrum fits the experimental data as closely as possible. Then, a minimization of the ${\chi }^2$ is defined as:

$${\chi }^2=\sum _{i=1}^N{\displaystyle \frac{{(y_i^{exp}-y_i^{model}\left(\overrightarrow{p}\right))}^2}{{ϵ}_i^2}}$$

(3)

beeing $\overrightarrow{p}$ a set of selected parameters for fitting, which minimizes the distance between the experimental data, $y_i^{exp}$ and the modeled data points $y_i^{model}$.

In Python, this process can be implemented using libraries like Lmfit [10], SciPy [11], or NumPy [12], which provide robust tools for processing data, graphics, and least-squares curve fit. The general approach involves defining a model function that represents the expected shape of the Mössbauer spectrum, which could be a sum of multiple Lorentzian or Gaussian functions, depending on the number and type of spectral components. Lorentzian functions are preferred with crystalline samples, whereas PseudoVoigts (sum of Lorentzian and Gaussian functions) are appropriate for Fe sites in disordered materials.

The least-squares fitting algorithm iteratively adjusts the model parameters until the sum of squared residuals between the experimental and calculated spectra is minimized. An example of a Lorentzian function definition for Python, can be seen next,

⠀

from scipy.optimize import curve_fit

import numpy as np

⠀

# Define the model function (e.g., a sum of Lorentzian components)

def lorentzian(x, amplitude, center, width):

    return amplitude * (width**2 / ((x - center)**2 + width**2))

⠀

# Generate synthetic data for fitting

x_data = np.linspace(-10, 10, 500)  # Example velocity range

y_data = lorentzian(x_data, 1, 0, 1) + np.random.normal(0, 0.02,

        size=len(x_data))

⠀

# Perform least-squares curve fitting

initial_guess = [1, 0, 1]  # Initial guess for parameters

params, covariance = curve_fit(lorentzian, x_data, y_data,

        p0=initial_guess)

⠀

# params contains the optimized parameters: amplitude, center,

        and width

⠀

This example demonstrates fitting a single Lorentzian function to synthetic data, though more complex models involving sums of Lorentzians or other spectral shapes can be implemented to represent real Mössbauer spectra.

3. Description of the Python code for Google Colab (PyMossFit)

The Python code is available in several Jupyter notebooks as typical examples found in practice [13]. In their names, some text is added that indicates the number and type of interactions (i. e. 1Q means one quadrupolar interaction, 2X for two magnetic sites, and so on). Some selected parts of it are described throughout this page.

The code is structured in three cells. The first includes the installation in the Google Colab environment of Lmfit, the core of data fitting. Likewise, it imports some packages of Scipy, Pandas, Matplotlib and Numpy, among others. The next step includes a Drive connection that asks the user for permission.

The second cell reads the datafile (format should be inspected previously to define "delimiter", "columns" and "skiprows" parameters). The required inputs are date (in a YYYYMMDD format) and maximum velocity associated to the extreme channels, in order to define the velocity scale.

In this cell, spectrum folding is performed with a Discrete Fourier Transforming (DFT) routine, the Numpy fft. The theory of this procedure corresponds to the Nyquist-Shannon Sampling Theorem that helps to determine a folding channel from the symmetry of discrete Fourier spectra [14].

Data can also be smoothed using a Savitsky-Golay package (I use savgol from Scipy). After folding, a new datafile is saved for fit with the use of the next cell.

The next figures illustrate the partial plots of data management with this procedure. This spectrum corresponds to zinc-ferrite nanoparticles obtained by a wet chemical method [15], which present a spinel structure with a tetrahedral site for $F{e}^{3+}$ and a octahedral site for both, $F{e}^{2+}$ and $F{e}^{3+}$, cations.

When raw data are plotted, they look as shown in Figure 3 the following because of the collection of absorption resonances converted to signals in a single-channel analyzer output, while the source moves back and forth. Figure 4 and Figure 5 show the DFT plot and the final folded spectrum, respectively. In the final folded spectrum, the horizontal sequence of channels was transformed to a $mm/s$ scale based on a calibration with a metallic ⁵⁷$Fe$ thin foil sample.

The set of physical parameters such as linewidth ($\Gamma$), isomer shifts (IS), quadrupole splitting (QS) and magnetic hyperfine field ($B_{HF}$) are calculated after adjusting the amplitud (a), full width at half maximum (b), centroid (m), line shift (d) and line separation (q). The initial set of these parameters can be fitted or fixed by selecting the "True" or "False" options, respectively.

A typical fit report of the output of this cell looks as follows:

Figure 6. An example of a fit report output of a Jupyter notebook in Google Colab 

Figure 6. An example of a fit report output of a Jupyter notebook in Google Colab 

Finally, the code saves fitted parameters, experimental and modeled sub-spectra in CSV formatted output files.

3.1. Some Examples

The following figures illustrate the fit results obtained from the Mössbauer spectra of different samples. Each exhibits several characteristics corresponding to the short-range-order neighborhood of ⁵⁷$Fe$ probes. Figure 7 corresponds to the fitted spectrum of an aged cathodic material $LiFeP{O}_4$ (LFP) where iron-rich secondary phases can clearly be identified. These undesired phases, which are responsible for the lowering of storage capacity, precipitated after extended thermal and environmental degradation. The fitting procedure took less than 40 secs after an appropriate selection of range parameters.

In contrast, Figure 8 shows a spectrum of a pure LFP sample which was synthesized by a solid-state reaction from analytical grade precursors.

The Figure 9 corresponds to the fitted spectrum with PyMossFit of a soil sample taken by the Mössbauer spectrometer on board the Mars Opportunity Rover. The data is available from [16] and this set is named the "Bounce Rock Case" (B0RR067). The iron-rich detected phase was identified with the code and corresponds to the pyroxene phase ($(Ca,Mg,Fe)Si{O}_3$), with two different octahedral Fe sites, M1 and M2 respectively [17,18].

The final example of the use of the PyMossFit notebook is the observed in Figure 10 and is the spectrum of an Iguazú fall soil sample. It corresponds to a complex spectrum where the use of a combination of sextets, doublets, and a singlet was required. All these subspectra correspond to different characteristics of iron oxides and hydroxides present in this kind of soils, originating from the erosion of basalt rocks.

3.2. Match of Fitted Parameters with a Local Database

A final cell allows us to identify the phases from an user’s own database or a public ones [19]. This procedure is based on the minimization of Eulerian distances of the set of parameters with the bibliographic collection, which can be updated or extended by the user. Some lines that interpret ranges of hyperfine parameters from the database are also included. The phases cited in the examples above were identified or verified using this code section.

The machine learning algorithm used to guess the present phases is based on a K-Nearest Neighbors (sklearn.neighbors) [20]. The output is a recommendation with a ranking of suggested phases tha match the best.

Figure 11. A match report obtained with a KNN algorithm for phase identification 

Figure 11. A match report obtained with a KNN algorithm for phase identification 

4. Conclusions

PyMossFit is a tool for Mössbauer spectroscopy analysis that explodes cloud computing. It was coded in Python and run in Jupyter notebook. The code requires the installation of the Lmfit package as the core of model fitting. PyMossFit was designed to work in Python environment at cloud computing services, such as Google Colab, being a good choice for user update, free OS and device dependencies. Likewise, it simplifies the collaborative work during spectra analysis with the addition of a code section for phase identification based on a KNN algorithm.