Exploring dynamics, disorder, and anharmonicity in THz spectra of crystalline materials

Terahertz spectral measurements of crystalline materials are known to be particularly sensitive to both the material’s long-range order and lattice dynamics. Here we explore a range of materials including common crystalline organic materials such as α-lactose monohydrate and l-cysteine through to engineered materials including topological insulators (Bi2(Te(1−x)Sex)3) where structure and disorder can be more finely controlled. By comparing variable temperature THz spectral measurements to ab-initio simulations using a range of methods and structures we can begin to unpick the origins of these spectra, and how they are influenced by dynamics and disorder across the lattice.


INTRODUCTION
][5] This has now gone beyond simply understanding the nature of the harmonic phonon vibrations and associated oscillator strengths for a specific material to incorporate other effects including scattering, the influence of sample morphology and orientation 6,7 through to the effects of anharmonicity, 8-10 temperature, 11,12 disorder, and molecular dynamics. 5,13 n this article, we present the work on a range of crystalline materials where a number of computational methods are used to begin to unpick the origins of the measured spectral properties of each.

α-LACTOSE MONOHYDRATE
α-Lactose monohydrate is a crystalline disaccharide that has become commonly used as a standard sample to measure at THz frequencies. 14It is particularly useful because it has a very low-frequency phonon mode (0.52 THz) which is insensitive to temperature, while other phonons above 1 THz change significantly both central position and width as a function of temperature. 14,15 s this material has become an important THz standard, it is crucial that we fully understand the effects of both temperature and sample preparation on the resultant spectra.THz spectra were measured on a THz-TDS instrument that has been previously described in detail elsewhere. 16The sample was a pellet of 10% α-lactose monohydrate by weight, mixed with polytetrafluoroethylene (PTFE) powder (with a particle size 1 µm) and pressed into a solid pellet, with a copper ring support, with a total thickness of ∼400 µm using a force of ∼8 t.The sample was then placed in an Oxford Instruments Microstat-HE cryostat and measurements taken from 4 K to 300 K. Resultant absorption spectra for two temperatures can be seen in Figure 1.
Figure 1: The measured and calculated absorption coefficient (cm −1 ) for a 10% α-lactose monohydrate powder mixed with PTFE at different temperatures.In (a) the spectrum is recorded at 4 K and the calculation is a harmonic phonon calculation and includes no temperature effects.(b) shows the experimental spectrum of the same sample recorded at 300 K where the calculated spectra has used the Quasi-Harmonic Approximation (QHA) to calculate the harmonic phonon spectrum at a fixed unit cell volume of 771.86 Å3 corresponding to QHA determined temperature of 300 K As a starting point for Density functional theory (DFT) calculations a crystal structure of α-lactose monohydrate recorded at 150 K was taken from the CCDC. 17Calculations were performed using VASP 18 with the Perdew-Burke-Ernzerhof (PBE) functional 19 and the Projector Augmented Wave (PAW) pseudopotentials distributed with VASP 5.4.4. 20Dispersion corrections were included using the Grimme DFT-D3 method. 21A Monkhorst-Pack 22 k-point grid of 7 1 5 with a total of 35 k-points and a cut-off energy of 900 eV was used based on convergence testing.Geometry optimisation of the structure resulted in a unit cell with a volume 13% lower than the experimental volume.This is partly owing to the lack of any temperature effects in the DFT calculation which will generally lead to a structure closer to an experiment at 0 K.It should be noted that if no dispersion correction is included this results in a cell expansion of more than 40% during optimisation and demonstrates the importance of the inclusion of dispersion corrections in these calculations.Once optimised, Phonopy 23 was used with VASP as the computational engine to calculate the phonon frequencies at the Γ-point using a finite difference approach while Born Charges were calculated using VASP.These calculations were post-processed using PDielec 6,7 where we assumed that the sample was a powder sample of spherical particles with 10% by weight of α-lactose monohydrate in PTFE using a Maxwell-Garnett (MG) effective medium approximation (EMA) 24,25 with peak widths chosen to best correlate with the experimental measurement.Mie-scattering 6,26,27 of air voids within the sample was also included by changing parameters until the background of the calculation matched with the rising background of the spectra.In this case, this correlated to a 7% volume of air with a void size of 10 µm which seems reasonable compared to previous measurements. 28Figure 1a shows the comparison between the measurement at 4 K and the calculated spectra.
In general the correlation between experiment and calculation in Figure 1a is very good.The majority of experimental peaks can be attributed to the nearest calculated feature with the only real discrepancy in the region between 2.5 THz and 3 THz.Here a peak in the experiment at 2.86 THz does not correspond to a calculated phonon mode and the experimental feature at 2.66 THz is 0.2 THz lower than the calculation.It should also be noted that while the relative intensities for both plots are similar the calculated spectra systematically under estimates the intensity of each mode.
To explore the influence of temperature of the THz spectra of α-lactose monohydrate we performed Quasi-Harmonic Approximation (QHA) 29,30 calculations using VASP and Phonopy.These calculations perform calculations at a number of fixed volumes and use a fit to an equation-of-state (EOS) to determine a volume-temperature plot which can be used to determine the temperature-dependent phonon spectra. 23,31 mportantly, this takes into account the influence of thermal expansion on the phonon spectra but still assumes that each phonon is harmonic in nature.To do this, we first performed a geometry optimisation using the same VASP settings as above for five fixed cell volumes on either side of the optimised volume.For each cell, a finite difference calculation was performed using Phonopy to calculate the phonon frequencies at the Γ-point for each volume, while the Born charges for each optimised cell were calculated using VASP.The tools provided with Phonopy were then used to fit the energy-volume information generated to the Rose-Vinet EOS which allows for a number of parameters to be determined including the temperature-dependent phonon spectra.For comparison with the optimisation above, the experimental unit cell volume determined at 150 K is 768.97Å3 which is a difference of only 0.016% compared to the experiment.
From the fit to the EOS, a final calculation was performed at a cell volume of 771.86 Å3 corresponding to a temperature of 300 K where the calculation is compared to the experimental spectra recorded at the same temperature in Figure 1b.As with the spectra at 4 K this assumes an MG-EMA with 10% α-lactose monohydrate in PTFE.Peak widths were chosen to best correlate with the experiment, with mie-scattering of air voids using the same size and volume used in Figure 1a.As with the low-temperature comparison, the correlation between the experiment and calculation in Figure 1b is very good.Interestingly, the correlation between experimental and calculated intensity is much closer than at low temperature and generally speaking the shift of peak position with temperature is well reproduced with a couple of exceptions.The first is the mode at 0.52 THz here at low temperature we underestimate the measured peak position while at 300 K the calculation overestimates the peak centre.The second is the peak at 2.55 THz which has a poor correlation with calculation both in terms of peak centre and intensity.This is the same peak that correlated poorly with calculation at low temperature.Further work is needed to include additional anharmonic effects to explore the origins of the discrepancy with calculation and experiment in both of these spectral features.

TOPOLOGICAL INSULATORS
Topological insulators (TIs) are a class of materials where the bulk band gap is bridged by conducting surface states. 32Here we explore the ternary alloy Bi 2 (Te (1−x) Se x ) 3 which is a well-known TI material family. 33These materials have a very strong degenerate phonon which can be easily tuned by varying the value of x by varying the ratio of tellurium and selenium fluxes during the molecular-beam epitaxy (MBE) growth. 34Of particular interest in this work is that a) The prominent phonon, somewhat counter-intuitively, red-shifts upon cooling while other Raman active phonons do not 34,35 and b) we observe a reduction in crystalline anharmonicity caused by the preferential ordering with the x = 1 3 showing the lowest anharmonicity. 34The fine control of alloy fraction during growth and the interesting temperature response, therefore, make this an excellent material to explore how anharmonicity (and defects) influence the spectral properties of this ternary alloy.TI samples were grown on [0001] oriented Al 2 O 3 substrates by co-deposition of evaporated bismuth, selenium, and tellurium which has been described in detail previously. 34The THz-TDS 16 measurements were undertaken between 3 K and 300 K in an Oxford Instruments Microstat-HE cryostat.Additional measurements were taken of a blank Al 2 O 3 substrate without TI film as a reference where the THz conductivity of the film can then be determined. 36Measured Real and Imaginary conductivity of the Bi 2 Te 3 film can be seen in Figure 2a and Figure 2c respectively along with Drude-Lorenz fits in red.
Phonon calculations were performed using both VASP 18 and Phonopy 23 for Bi 2 Te 3 and Bi 2 Se 3 using the conventional cell.The PBE functional and PAW pseudopotentials, D3-BJ 21, 37 dispersion corrections.A Monkhorst-Pack 22 k-point grid of 9 9 1 and a cut-off energy of 520 eV was chosen based on convergence testing.PDielec 6,38 was used to calculate the absorbance (Figure 2b) and refractive index (Figure 2d).PDielec uses a generalised transfer matrix method 39 to calculate the response of a thin crystalline film of Bi 2 Te 3 on a Al 2 O 3 with the THz beam incident on the (003) for the conventional cell which is equivalent to the (111) face of the primitive cell.This shows a good overall correlation with the experiment although the central frequency of the strong phonon is overestimated by a little more than 0.5 THz.As discussed above the temperature-dependent phonon shift of these materials shows an interesting red-shift.Figure 3a shows the response of this phonon mode in Bi 2 Se 3 with the central frequency change, (∆ν), as a function of temperature shown in Figure 3b as the black data points.Initially, we began to use the QHA method described in Section 2 with multiple unit cells of different volumes used to fit the Rose-Vinet EOS.While this method worked well for α-lactose monohydrate the temperature dependence determined for Bi 2 Se 3 (Figure 3b, blue data points) shows the incorrect trend.To explore this further we then considered the anharmonicity of the motion associated with the phonon mode in these materials.Here we use the method described by Skelton et al. 8 and begin by mapping out the 1D-potential-energy surface of the phonon mode eigenvector and calculate the energy as a function of displacement along the normal-mode coordinate.By fitting this potential energy surface to a 6-order polynomial and then solving a 1D Schrödinger equation for the polynomial we determine the eigenvalues and anharmonic vibrational partition function.This, in turn, allows an effective renormalised harmonic frequency for the phonon mode at a specific temperature to be determined.While this does not include effects such as direct coupling between modes that more advanced methods like VCI and VSCF can include, 40,41 this is a method that is commonly used to correct for some anharmonic effects of low-frequency modes in many materials. 4,8,9 Te green data points in Figure 3b show the temperature dependence of this renormalised harmonic frequency for the phonon mode and while the magnitude of the change is smaller than seen experimentally, both show appropriate red-shifting of the mode.We have previously used this method to explain that change in crystalline anharmonicity caused by the preferential ordering as x is changed but it is worth noting that the anharmonicity is likely to be much more complicated as this model does not account for phonon-phonon, electron-phonon or spin-orbit coupling, or thermal expansion effects (which QHA does but shows the opposite trend), which play an important role in the temperature-dependent phonon dynamics of these materials.We continue to explore the complexity of anharmonicity in this ternary alloy.

L-CYSTEINE
In the TI samples discussed in Section 3 the substitution of Se with Te is used to influence both the central frequency and anharmonicity 34 of the key phonon mode, which could be viewed as the influence of defects in the crystal structure on the measured spectral parameters.When looking at molecular systems as discussed in Section 2 we tend to assume that we consider a perfect static crystal with no defects.3][44] Previous attempts to calculate the THz spectra of l-cysteine have shown poor correlation with the experiment. 45Here we explore a range of methods to calculate the THz spectra of l-cysteine in order to improve the correlation between measurement and calculation while extending our understanding of how anharmonicity, dynamics and disorder influence the spectral properties.
We have used a combination of static and dynamic calculations to explore the disorder and stability of all four forms of l-cysteine and detailed information about the calculations of all forms can be found in Kendrick et al. 13 For Form I, initially, static calculations of unit cells containing both the SH...S (the most stable form) and SH...O configurations were performed using VASP 18 with similar settings to those described in Section 2 apart from the inclusion of the D3-BJ dispersion corrections instead. 21,37 honon frequencies were also calculated using VASP via Density Functional Perturbation Theory (DFPT). 18A cutoff energy of 600 eV and a reciprocal space k-point resolution of 0.2 Å3 was used based on convergence testing.The resultant spectra are labelled as Static -SH...S and Static -SH...O in Figure 4.
Additionally, Ab Initio molecular dynamics calculations of larger supercells were performed using CP2K. 46As was used in VASP, all calculations employed the PBE functional, 19 with GD3/BJ dispersion correction. 21,37 he Molopt double-ζ valence plus polarization basis sets 47 were used with GTH PBE pseudopotentials. 48Calculations were performed on two different supercells.The first, labelled SC8, is a supercell of Form I containing eight different l-cysteine molecules in the SH...S molecular configuration.The second is a larger supercell, labelled DC32, that contains 32 different l-cysteine molecules where the SH...S and SH...O bonding motifs have been randomly assigned with a distribution of 50% of each motif, similar to what is expected in a room temperature crystal. 43Calculations were run for both supercells at 88 K and 300 K. Initial calculations were performed at constant pressure and temperature (NPT) for at least 3 ps to obtain average cell dimensions, after an initial equilibration at this fixed volume for at least 2 ps, production runs at constant volume and temperature (NVT) were then run for a further 30 ps.The IR spectra for each production run for each trajectory were then calculated from the dipole-dipole autocorrelation function of the cell dipole over the course of the trajectory. 49gure 4 shows a summary of these calculated spectra compared to IR and THz spectra of Form I of l-cysteine.Static spectra of both bonding motifs were calculated by post-processing the VASP calculations using PDielec 6, 7 in a similar way to Section 2. The sample was assumed to be a powder sample of spherical particles with 10% by weight of l-cysteine using a MG-EMA.A constant peak width was assumed for all peaks and chosen to best correlate with the experimental measurement, as such Figure 4a uses a peak width of 10 cm −1 while all other plots including static calculations assume a peak width of 5 cm −1 .Importantly, in the dynamic calculations, all motion is inherently included and therefore the peak widths shown are largely a result of the coherence time of the motions within the calculation. 13gure 4a shows the static calculation of the SH..S bonding motif compared to the dynamic calculations of both supercells at 300 K along with the experimental spectra taken at 300 K digitised from Mink et al. 50It can be seen that all three calculated spectra are largely similar in both peak position and intensity and correlate well with the experiment.The only real difference is emphasised in the inset where the inclusion of anharmonicity and calculated linewidth from the dynamic calculation shifts all modes to lower frequencies and slightly improves correlation with the experiment, particularly below 20 THz. Figure 4b compares the calculated static IR spectra of both bonding motifs and uses a narrower linewidth of 5 cm −1 to highlight differences between the spectra.The two spectra are very similar and it would be difficult to suggest either was better correlated with the experiment apart from a peak in the SH...O motif spectra at 30 THz which is much weaker in both the SH...S calculation and experiment.In short, the IR spectra are largely insensitive to either molecular dynamics or disorder within the crystal.This is clearly not the case for the THz spectra which are shown in Figure 4c and Figure 4d which show a comparison to experimental spectra digitsed from Ren et al. 51 taken at 88 K and 300 K respectively.The static spectra from either bonding motif show little correlation to experimental spectra at either temperature, this is particularly surprising for the SH...S bonding motif which is expected to dominate the crystal structure at low temperature.Correlation between experiment and calculation is improved when compared to the dynamic calculation of the SC8 spectra at 88 K suggesting the dynamics and anharmonicity play a significant role in the origins of the spectral response in l-cysteine, but the correlation is reduced significantly at the higher temperatures.The spectra calculated from the dynamics of the DC32 cell correlates best with the experimental spectra at both temperatures, highlighting the role that crystalline disorder also plays in the THz spectra of this material.

CONCLUSIONS
Here we have explored a number of calculation methods to aid the interpretation of the measured THz spectral parameters of a number of very different crystalline materials.For α-lactose monohydrate we showed that static harmonic phonon calculations correlate well with low-temperature experimental spectra although generally underestimate the spectral intensity of the modes and the changes of spectra with temperature can largely be explained by the quasi-harmonic approximation.For the ternary alloy, Bi 2 (Te (1−x) Se x ) 3 , we show that the main phonon frequency counter-intuitively red-shifts with decreasing temperature.Static DFT calculations correlate well with the experiment but tend to over-predict the frequency of the key mode.Here, QHA fails to aid the interpretation of the temperature dependence of the spectra which is better explained by the shape of the potential-energy surface of the mode which can be influenced by doping in the material.Finally, we showed that for complex molecular crystals like l-cysteine where disorder is present, then both anharmonicity, molecular dynamics and disorder in the crystal must all be taken into account to aid the interpretation of the experiment.What needs to be highlighted here is that interpreting THz spectra can be complex, and a one-size fits all approach to materials is not available, however, great progress is being made where both anharmonic and dynamic effects can be taken into account in calculations of larger and larger systems.

Figure 2 :
Figure 2: Shows the measured and calculated spectral parameters for a thin film of Bi 2 Te 3 on a Al 2 O 3 substrate.a) shows the measured real conductivity in black with Drude-Lorenz fit in red.b) shows the calculated Absorbance for the same film.c) shows the measured imaginary conductivity in black with Drude-Lorenz fit in red.d) shows the calculated real refractive index for the same film.

Figure 3 :
Figure 3: a) Shows the temperature-dependent real conductivity of a thin film of Bi 2 Se 3 with an inset schematic of the relevant phonon mode.The positions C 1 and C 2 labelled are sites for Se and Te.In Bi 2 Se 3 , all three sites are occupied by Se but as Te is added this preferentially fills site C 2 before C 1 .(b) Shows the change in central frequency (∆ν) of the phonon mode for the experiment (black) shown in (a), calculations using QHA (blue) and the fit of the anharmonic potential (green) discussed in the text.

Figure 4 :
Figure 4: a) Shows the IR experimental spectra of l-cysteine 50 compared to the static calculation with the SH...S bonding motif and dynamic spectra calculated from supercells SC8 and DC32 at 300 K. b) shows the static IR spectra calculated from the SH...S and SH...O bonding motifs.c) shows the experimental THz spectra of l-cysteine recorded at 88 K 51 compared to the static calculations of both bonding motifs and dynamic calculations for both supercells at 88 K and d) shows the same as (c) but with experimental 51 and dynamic calculations at 300 K