We will now work with Raman spectroscopy data from the luminescence of a ruby particle under high pressure: Data/rubis_xx.txt. This measurement is performed in a diamond anvil cell with a laser of wavelength \(\lambda_L\)=568.189 nm, and the pressure is known by calibrating the frequency of a diamond line.

The ruby luminescence contains two lines, and we want to use the position of the R1 line, which is the most intense and has the highest energy, to determine the pressure. We therefore need to determine the parameters \(A\) and \(B\) of the equation of state linking wavelength (in nm) and pressure (in GPa):

\[ P~[GPa] = \frac{A}{B}\left[\left(\frac{\lambda}{\lambda_0}\right)^B-1\right], \] where \(\lambda_0\) = 694.24 nm is the wavelength of the R1 line at zero pressure, and \(\lambda\) is the wavelength of the R1 line at \(P\) pressure. For a given line, the relationship between its measured Raman shift, \(\omega\), and its wavelength is given by \(\lambda = \frac{10^7}{\omega_L - \omega}\), with \(\omega_L= \frac{10^7}{\lambda_L}\) the frequency of the probe laser.

Open one of the files, for example Data/rubis_01.txt to see its structure - they’re all the same. They are made up of two columns, which are the Raman shifts in cm-1 and the intensities. The nomenclature of the files "Data/rubis_xx.txt" is such that xx is the pressure in GPa.

First, consider the steps involved in determining the parameters \(A\) and \(B\) of the equation of state.

  • What data do you need?
  • What steps do you need to take to obtain them?

  • Let’s start by defining the function Prubis(w,A,B), which returns the pressure in GPa as a function of the Raman shift in cm-1, and the parameters \(A\) and \(B\).

Find in the folder “Data” the list of files with the pattern “rubis_” in their name (check out the list.files() function) and store it in flist


What is the length of flist?


Initialize an empty tibble called spec that will store all spectra:


Using a for loop, load each file in flist and store them into spec by making it tidy.

Another version of this that does not use a for loop would be like so (try to understand what’s happening in this code by running each line one after the other and looking at the help on the functions used):


Modify spec to get the pressure in GPa from the file name, and so that the intensity column is normalized to [0,1].


Plot the first spectrum in spec to see the x range


Using ggplot2, plot all spectra normalized to [0,1] (with points) and stacked on top of each other with a vertical shift corresponding to the pressure in GPa


What we want to do now is to fit this data by the sum of two Lorentzians to find the precise position of the R1 line as a function of pressure. To do this, we’ll use the function nls(), which allows us to make a non-linear fit of the data. We also need to define the function lor(x,x0,gamma) which returns a Lorentzian centered at x0 and of width at half-height gamma. Start by defining the function lor(x,x0,gamma).


Now create an empty tibble called results, in which we’ll store all the fitting parameters in a tidy tibble. Then fit each spectrum by the sum of two Lorentzians, and save all fit parameters in results, along with the file name and pressure. For the fit to work, you’ll need to give reasonable initial values for the fit parameters using the start=list(...) argument of nls(). For example, you can use the which.max() function to find the index of the maximum of a vector.

Another way of doing that is to use the map() function from the purrr package, which allows us to apply a function to each element of a list. Here, we want to apply the nls() function to each spectrum in spec, and store the result in a new column of spec. We can do this like so:


Plot the result of the fit on the graph above.


Now all you have to do is plot the position of the peak as a function of pressure, and model this data by the equation of state defined at the beginning. What values of A and B do you obtain? Compare with the literature values: A = 1876 GPa and B = 10.71.