Peapods

Download data files

---

1. All plots should be performed using ggplot2.
2. Load the libraries tidyverse and broom and set the global ggplot2 theme to theme_bw().

---

A bit of context

Important

Don’t worry, you don’t need to know anything related to it to do the exercise!

Figure 1: A TEM micrograph of carbon nanopeapods.

In this exercise, we’re going to determine the activation energy of the rotation of C₆₀ fullerene molecules in carbon nanopeapods (Figure 1), i.e. C₆₀@SWNT. For this, we’ve performed Inelastic Neutron Scattering spectroscopy at the Institute Laue Langevin in Grenoble. The signal we’ll be working on is a quasi-elastic signal arising from the rotational motion of the fullerenes, and you’ll see that the signal is evolving with temperature on a certain range of temperatures, while it’s constant below a critical temperature.

This is due to the fact that the rotation of the molecules is an activated process: below a certain temperature the thermal fluctuations are not enough to overcome the energy barrier, and the molecules rotations are locked. Above the critical temperature, the rotations are allowed, and the molecules show a diffusive rotation with a statistics following an Arrhenius law.

Note

The width of the quasi-elastic signal is directly related to the amplitude of the rotations. In the following, we want to determine:
(1) the evolution of the signal’s width with temperature, and then
(2) determine the activation energy of the process by fitting an Arrhenius law to this evolution.

Data wrangling

1. Using list.files(), find all xxK.csv files in the folder Data/, where xx is the temperature. You’ll store the result in a vector called flist.

2. These files all contain 2 columns, w and int. Using read_csv(), load all these files into a tidy tibble called df, with the columns T, w and int. Using a combination of as.numeric(), basename() and str_replace(), make sure the column T contains the temperatures as numerical values.

Note

In case you didn’t manage to get there, you can continue by running:

df <- read_csv("Data/alldata.csv")

Plotting and fitting the data

3. Now we have loaded our data, let’s take a look at it. Try to make a plot looking like this:

colors <- colorRampPalette(c("black","royalblue","seagreen","orange","red"))
df |> 
    ggplot(aes(x = w, y = int + as.numeric(factor(T))*.1, color=T)) +
       geom_point(alpha=0.5, size=.5)+
       scale_color_gradientn(colours = colors(50))+
       labs(x="Energy Transfer [meV]",
            y="Intensity [arb. units]",
            color="T [K]")

4. Now we see that all signals are Gaussian peaks centered around 0 meV with a varying width. A Gaussian is defined using the function dnorm(w, mean=w0, sd=FW). Modify the following code to fit all the peaks using the nls() function.

5. Now let’s check that our fit is good by adding a red line to the previous plot, the red line being the result of the fit. For this, we will use the .fitted column of the unnested augmented column.

6. Let’s simplify our df_fit tibble so that it only contains 3 columns, T, FW and dFW, FW and dFW being the result of the fit and the corresponding standard error.

Fitting the Arrhenius law

Note

In case you didn’t manage to get there, you can continue by running:

df_fit_sorted <- read_csv("Data/FW_fit_results.csv")

7. Now let’s plot our results to see the evolution of the width with temperature. Add vertical errorbars corresponding to the standard error on the fit.

df_fit_sorted |> 
    ggplot(aes(x = T, y = FW, ymin=FW-dFW, ymax=FW+dFW)) +
       geom_point(size=1)+
       geom_errorbar(width=0)+
       labs(x = "Temperature [K]",
            y = "Width [meV]")

8. You see that below a certain temperature, there is no evolution of the width, as discussed at the beginning. Create a subset of the df_fit_sorted tibble only containing data for which the width evolves with temperature.

9. Now let’s perform the fit of the data using an Arrhenius law like below, where FW is the width of the signal. You can either consider it a nonlinear model (and fit it with nls()) or a linear one (and fit it with lm(). In that case, you’ll need to add a column containing the 1/T values to perform your fit).

\[ FW = FW_0 \times e^{-E_a/k_BT} \]

10. Plot the result of the fit to make sure it’s good, like so:

11. What is the value of the activation energy, in Kelvins (with standard error)? Does it make sense to you when seeing the width evolution with temperature?