Magnetic nanoparticles¶
We here show how to analyse (a subset) of the data on $\alpha\text{-Fe}_2\text{O}_3$ nanoparticles presented in J. Chem. Phys. 140, 044709 (2014).
$\alpha\text{-Fe}_2\text{O}_3$ is an antiferromagnet with magnetic moments in the hexagonal $ab$ plane. Magnetic nanoparticles are small enough to be single domain, and all magnetic moments therefore move coherently. Because the nanoparticles are so small, the magnetization can spontanously flip. This is called superparamagnetism, and if the characteristic flipping time is given by $\tau$ then the neutron signal from the magnetic structure is given by a Lorentzian with width (half width at half maximum) of $\Gamma = \hbar/\tau$:
$$ I_\text{SPM} = \frac{A_\text{SPM}}{\pi}\frac{\Gamma}{\Gamma^2+\hbar^2 \omega^2}. $$
The nanoparticles also exhibit quantized spin-wave excitations. We here observe only two ${\bf Q}=0$ excitations, labelled $\hbar \omega_+$ and $\hbar \omega_-$, corresponding to movement within the $ab$ plane and perpendicular to it, respectively. These excitations occur at the magnetic Bragg peaks, which are the (003) ($|Q|=1.37$ Å$^{-1}$) and (101) ($|Q|=1.51$ Å$^{-1}$) reflections. The signal from these excitations can be well described by damped harmonic oscillators: $$ I_\pm = \frac{A_\pm}{\pi} \frac{2\gamma_\pm \hbar \omega_\pm^2}{(\hbar^2 \omega^2 - \omega_{\pm}^2)^2 + 4\gamma_\pm^2 \hbar^2 \omega^2}. $$ Here, $\gamma_\pm$ is the intrinsic width of the excitation
We are here interested in determining the flipping time, $\tau$ and the energy of the two excitations, $\omega_\pm$.
Several experiments were carried out, and we here analyse only a subset of the data, obtained at IN5 at ILL in Grenoble, France. We have the measured intensity as function of energy for four different values of $Q$; at the two magnetic peaks and at two positions where only background is present. The background is complicated and will need to be interpolated based on those measurements. We have data at 1.5 K, where the dynamics are almost completely frozen, to determine the resolution function, and at 150 K where we want to fit the data.
The background is complicated because water gets adsorbed on the surface of the nanoparticles. Some of the water is motionless and gives rise to elastic scattering, modelled by a delta function. Part of the water is mobile and therefore gives rise to quasielastic scattering, which is modelled by a Lorentzian. In total, the scattering from water is given by $$ I_W = A_W \delta(\omega) + \frac{A_\text{W}}{\pi}\frac{\Gamma_W}{\Gamma_W^2+\hbar^2 \omega^2}. $$
Enough background, let's start looking at the data! We first import everything that we need.
import pooch
import scipp as sc
from easydynamics.analysis.analysis import Analysis
from easydynamics.experiment import Experiment
from easydynamics.sample_model import ComponentCollection
from easydynamics.sample_model import DampedHarmonicOscillator
from easydynamics.sample_model import DeltaFunction
from easydynamics.sample_model import Gaussian
from easydynamics.sample_model import Lorentzian
from easydynamics.sample_model import Polynomial
from easydynamics.sample_model.background_model import BackgroundModel
from easydynamics.sample_model.instrument_model import InstrumentModel
from easydynamics.sample_model.resolution_model import ResolutionModel
from easydynamics.sample_model.sample_model import SampleModel
from easydynamics.utils.utils import hbar
# Make the plots interactive
%matplotlib widget
We begin with the resolution. We create an Experiment, download the data using pooch and attach it to the Experiment.
resolution_experiment = Experiment(display_name='Nanoparticles, 1.5 K')
file_path = pooch.retrieve(
url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/tutorial2/docs/docs/tutorials/data/nano_1p5K.h5',
known_hash='32935cfaa5e623137ea5f4ba972d0c921bca8d0e067486666cdcf23abf27823e',
)
resolution_experiment.load_hdf5(filename=file_path)
Downloading data from 'https://github.com/easyscience/dynamics-lib/raw/refs/heads/tutorial2/docs/docs/tutorials/data/nano_1p5K.h5' to file '/home/runner/.cache/pooch/c1cea31a02b0d40173eb269468e1115a-nano_1p5K.h5'.
/home/runner/work/dynamics-lib/dynamics-lib/.pixi/envs/default/lib/python3.12/site-packages/requests/__init__.py:113: RequestsDependencyWarning: urllib3 (2.6.3) or chardet (7.0.1)/charset_normalizer (3.4.5) doesn't match a supported version! warnings.warn(
We now take a look at the data. It is helpful to turn on the log scale since any variations in the data are tiny.
resolution_experiment.plot_data(slicer=True, logy=True)
We will use this data to determine the resolution function. There is a large peak at zero energy transfer, and some small excitations away from the elastic line that we do not wish to deal with. We therefore remove those data points using scipp, and plot the result.
E_min = sc.scalar(-0.2, unit='meV')
E_max = sc.scalar(0.2, unit='meV')
resolution_experiment.data = resolution_experiment.data['energy', E_min:E_max]
resolution_experiment.plot_data(slicer=True)
The resolution function can to a good approximation be modelled as a Gaussian. In truth, it has small Lorentzian tails, but we leave it as an exercise for the interested reader to add this.
We define the SampleModel to be a DeltaFunction, since there is essentially only elastic scattering present. We also add a constant background using a Polynomial.
As in Tutorial 1 we place everything in an Analysis object and plot the start guesses.
delta_function = DeltaFunction(area=100)
res_sample_model = SampleModel(components=delta_function)
res_resolution_model = ResolutionModel()
res_components = ComponentCollection()
res_gauss = Gaussian(area=1, width=0.02)
res_gauss.area.fixed = True
res_components.append_component(res_gauss)
res_resolution_model.components = res_components
background_model = BackgroundModel()
polynomial = Polynomial(coefficients=[1.5])
polynomial.coefficients[0].min = 0.0
background_model.components = polynomial
res_instrument_model = InstrumentModel(
resolution_model=res_resolution_model,
background_model=background_model,
)
res_analysis = Analysis(
experiment=resolution_experiment,
sample_model=res_sample_model,
instrument_model=res_instrument_model,
)
res_analysis.plot_data_and_model()
The start guess looks reasonable, so let us fit everything and check out the result. The fit is not perfect, but it's certainly good enough for our purposes, giving excellent agreement over several orders of magnitude.
res_analysis.fit()
res_analysis.plot_data_and_model()
With our resolution in hand, it's time to look at the data. As before, we download the data file using pooch and make an Experiment.
experiment = Experiment(display_name='Nanoparticles, 150 K')
resolution_experiment = Experiment(display_name='Nanoparticles, 1.5 K')
file_path = pooch.retrieve(
url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/tutorial2/docs/docs/tutorials/data/nano_150K.h5',
known_hash='dc3abdaf6d2165f4dbc5df4e0de1890af8d4cd290f2c79ac7b319ccbe44117b8',
)
experiment.load_hdf5(filename=file_path)
experiment.plot_data(slicer=True, logy=True)
Downloading data from 'https://github.com/easyscience/dynamics-lib/raw/refs/heads/tutorial2/docs/docs/tutorials/data/nano_150K.h5' to file '/home/runner/.cache/pooch/72214c0d025c20e2aeef7074fc27e8ea-nano_150K.h5'.
The magnetic excitations are around 1.1 meV and 0.3 meV, so we cut the data to $\pm1.5$ meV and plot the result. At the two lower $Q$, there is mainly the background signal from water as described above, whereas the signal at higher $Q$ is much more complex because of the magnetic scattering.
E_min = sc.scalar(-1.5, unit='meV')
E_max = sc.scalar(1.5, unit='meV')
experiment.data = experiment.data['energy', E_min:E_max]
experiment.plot_data(slicer=True, logy=True)
By now you know the routine: we create a SampleModel with the desired components, a BackgroundModel and an InstrumentModel using the resolution determined before. We furthermore set the temperature of the sample to be 150 K. This turns on detailed balancing, which has a small but non-zero effect on our fit. As always, we check out the start guesses before we fit.
sample_model = SampleModel()
water_delta_function = DeltaFunction(display_name='Water delta function', area=100)
water_lorentzian = Lorentzian(display_name='Water Lorentzian', area=10, width=0.2)
sample_model.append_component(water_delta_function)
sample_model.append_component(water_lorentzian)
sample_model.temperature = 150
background_model = BackgroundModel()
polynomial = Polynomial(coefficients=[0.15])
polynomial.coefficients[0].min = 0.0
background_model.components = polynomial
instrument_model = InstrumentModel(background_model=background_model)
analysis = Analysis(
experiment=experiment, sample_model=sample_model, instrument_model=instrument_model
)
analysis.instrument_model._resolution_model = res_analysis.instrument_model.resolution_model
analysis.instrument_model.resolution_model.fix_all_parameters()
analysis.plot_data_and_model(logy=True)
We carry out the fit and inspect the result. It's a decent fit at low $Q$, and obviously terrible at the magnetic positions. This is expected; we're only modelling the background and not the full signal yet. The fit at low $Q$ could also be improved, but it is good enough for our purposes.
analysis.fit()
analysis.plot_data_and_model(logy=True)
Now we want to figure out what the background from water looks like at the magnetic positions. To do this, we inspect the parameters describing it and look for patterns. It looks like the area of the delta function, the area of the Lorentzian and the width of the Lorentzian are almost identical for the two non-magnetic positions. That is promising, as we can then fix these values at the magnetic positions.
analysis.plot_parameters(names=['Water delta function area', 'Water Lorentzian area'])
analysis.plot_parameters(names=['Water Lorentzian width'])
We calculate a simple average of the relevant parameters, and fix these in our Analysis object at all Q. We plot the resulting model and see that it indeed fits well at low $Q$. At higher $Q$ it obviosuly does not describe all the signal, since there is magnetic scattering there as well.
The DeltaFunction is the first component (index 0), and the Lorentzian is the second component (index 1). It will soon be possible to refer to components by name as well as index. It will also be made easier to fix parameters at multiple $Q$, but for now we do it manually.
delta_0 = analysis.sample_model.get_component_collection(Q_index=0).components[0]
delta_1 = analysis.sample_model.get_component_collection(Q_index=1).components[0]
delta_area = (delta_0.area + delta_1.area) / 2
lorz_0 = analysis.sample_model.get_component_collection(Q_index=0).components[1]
lorz_1 = analysis.sample_model.get_component_collection(Q_index=1).components[1]
lorz_area = (lorz_0.area + lorz_1.area) / 2
lorz_width = (lorz_0.width + lorz_1.width) / 2
for Q_index in range(analysis.sample_model.Q.size):
delta = analysis.sample_model.get_component_collection(Q_index=Q_index).components[0]
delta.area = delta_area.value
delta.area.fixed = True
lorz = analysis.sample_model.get_component_collection(Q_index=Q_index).components[1]
lorz.area = lorz_area.value
lorz.width = lorz_width.value
lorz.area.fixed = True
lorz.width.fixed = True
analysis.plot_data_and_model(logy=True)
We are now ready to fit the magnetic signal, fixing the background from water as described above. We create a new Analysis object and add the water components and the components describing the magnetic signal as described above.
# Now make a new analysis with this sample model
mag_sample_model = SampleModel()
water_delta_function = DeltaFunction(display_name='Water delta function', area=100)
water_lorentzian = Lorentzian(display_name='Water Lorentzian', area=100, width=0.2)
mag_sample_model.append_component(water_delta_function)
mag_sample_model.append_component(water_lorentzian)
# Add all the magnetic components
DHO1 = DampedHarmonicOscillator(display_name='DHO1', area=5, center=0.35, width=0.2)
DHO2 = DampedHarmonicOscillator(display_name='DHO2', area=1, center=1.1, width=0.1)
mag_lorz = Lorentzian(display_name='Magnetic Lorentzian', area=30, width=0.01)
mag_sample_model.append_component(DHO1)
mag_sample_model.append_component(DHO2)
mag_sample_model.append_component(mag_lorz)
background_model = BackgroundModel()
polynomial = Polynomial(coefficients=[0.15])
background_model.components = polynomial
instrument_model = InstrumentModel(background_model=background_model)
# Create the analysis object
mag_analysis = Analysis(
experiment=experiment, sample_model=mag_sample_model, instrument_model=instrument_model
)
mag_analysis.instrument_model._resolution_model = res_analysis.instrument_model.resolution_model
mag_analysis.instrument_model.resolution_model.fix_all_parameters()
We now fix all the parameters descibing non-magnetic scattering to the values we found in the previous fit.
for Q_index in range(mag_analysis.sample_model.Q.size):
delta = mag_analysis.sample_model.get_component_collection(Q_index=Q_index).components[0]
delta.area = delta_area.value
delta.area.fixed = True
lorz = mag_analysis.sample_model.get_component_collection(Q_index=Q_index).components[1]
lorz.area = lorz_area.value
lorz.width = lorz_width.value
lorz.area.fixed = True
lorz.width.fixed = True
We also fix all the parameters describing the magnetic signal at low Q. We set the areas to 0 and the center and width to fixed values so that they will not appear in the fit. We also take a look at the model before the fit.
for Q_index in [0, 1]:
DHO1 = mag_analysis.sample_model.get_component_collection(Q_index=Q_index).components[2]
DHO2 = mag_analysis.sample_model.get_component_collection(Q_index=Q_index).components[3]
lorz = mag_analysis.sample_model.get_component_collection(Q_index=Q_index).components[4]
DHO1.area = 0.0
DHO1.center = 1.0
DHO2.width = 0.1
DHO1.fix_all_parameters()
DHO2.area = 0.0
DHO2.center = 1.0
DHO2.width = 0.1
DHO2.fix_all_parameters()
lorz.area = 0.0
lorz.width = 0.1
lorz.fix_all_parameters()
mag_analysis.plot_data_and_model(logy=True)
It looks reasonable, so let's fit it and see what happens.
mag_analysis.fit()
mag_analysis.plot_data_and_model(logy=True)
The fit is very good, but we can make one improvement. The high-energy excitation is very weak at $Q=1.37$ Å$^{-1}$. We therefore make the width and center depend on the width and center of the same excitation at $Q=1.51$ Å$^{-1}$. We do this by using the EasyScience make_dependent_on method, where the first argument is the dependency expression (in this case just a variable), and the second is a dictionary telling EasyScience what the parameters are - in this case we say, for example, that the width of DHO2_lowQ should be equal to a, where a is the width of DHO2_highQ.
DHO2_highQ = mag_analysis.sample_model.get_component_collection(Q_index=3).components[3]
DHO2_lowQ = mag_analysis.sample_model.get_component_collection(Q_index=2).components[3]
DHO2_lowQ.width.make_dependent_on('a', {'a': DHO2_highQ.width})
DHO2_lowQ.center.make_dependent_on('a', {'a': DHO2_highQ.center})
mag_analysis.fit(fit_method='simultaneous')
mag_analysis.plot_data_and_model(logy=True)
With the fits done, it's time to inspect the results and answer the questions in the beginning. We start with extracting $\tau$ using an average of the two measurements.
width1 = mag_analysis.sample_model.get_component_collection(Q_index=2).components[4].width
width2 = mag_analysis.sample_model.get_component_collection(Q_index=3).components[4].width
width = (width1 + width2) / 2
print(width1)
print(width2)
print(width)
tau = hbar / width
tau.convert_unit('ns')
print(tau)
<Parameter 'Magnetic Lorentzian width': 0.0013 ± 0.0001 meV, bounds=[1e-10:inf]> <Parameter 'Magnetic Lorentzian width': 0.0014 ± 0.0002 meV, bounds=[1e-10:inf]> <Parameter 'Parameter_194': 0.0013 ± 0.0001 meV, bounds=[1e-10:inf]> <Parameter 'Parameter_195': 0.4947 ± 0.0426 ns, bounds=[0.0:6582119.569509066]>
Finally, we look at the excitation energies. Since we fixed the high energy excitation to be the same for both Q, we in principle only need to look at one of them. We look at both just to double-check that they are indeed the same. We find that the high energy excitaion is at 1.078(8) meV, and the low energy excitation at 0.257(3) meV.
E_minus_003 = mag_analysis.sample_model.get_component_collection(Q_index=2).components[3].center
E_minus_101 = mag_analysis.sample_model.get_component_collection(Q_index=3).components[3].center
E_plus_003 = mag_analysis.sample_model.get_component_collection(Q_index=2).components[2].center
E_plus_101 = mag_analysis.sample_model.get_component_collection(Q_index=3).components[2].center
print(E_minus_003)
print(E_minus_101)
print(E_plus_003)
print(E_plus_101)
print((E_plus_003 + E_plus_101) / 2)
<Parameter 'DHO2 center': 1.0783 ± 0.0085 meV, bounds=[1e-10:inf]> <Parameter 'DHO2 center': 1.0783 ± 0.0085 meV, bounds=[1e-10:inf]> <Parameter 'DHO1 center': 0.2602 ± 0.0030 meV, bounds=[1e-10:inf]> <Parameter 'DHO1 center': 0.2527 ± 0.0046 meV, bounds=[1e-10:inf]> <Parameter 'Parameter_197': 0.2565 ± 0.0028 meV, bounds=[1e-10:inf]>
More analysis and interpretation of the parameters could be done. Some examples are:
- We could look at the ratio of the areas of the excitations and relate it to the direction of the movement compared to the direction of ${\bf Q}$, making use of the fact that neutrons only see the component of the moment that is perpendicular to ${\bf Q}$.
- We could look at the relative area of the excitations and the superparamagnetism and relate that to the anisotropies of the system
- We could also relate the excitation energies to the exchange constants and anisotropies of the material, and we can relate those to the superparamagnetic spin flips described by $\tau$.
However, the main purpose of this tutorial is to show how to use EasyDynamics, not how to do general data analysis. We instead refer the interested reader to the paper cited in the introduction.