Interpolating Sodium Chloride phonons

Sodium Chloride forms face-centred cubic crystals in Space group \(F m \bar{3} m\) (Hall symbol \(-F 4 2 3\)) with lattice parameter \(a \approx 5.69\) Å.

Euphonic is a project which can take force constant information at gridded positions in \(\mathbf{Q}\) and use Fourier interpolation to approximate the dynamical matrix at arbitrary \(\mathbf{Q}\). It can then solve the eigenvalue problem to determine the eigenvalues, \(\omega_i^2(\mathbf{Q})\), and eigenvectors, \(\hat{\epsilon}_{ij}(\mathbf{Q})\), which are the squared energy and vectors related to the atom-resolved displacements for the \({i=1,…,3 N_\text{atom}}\) phonon branches, where \({j=1,\ldots,N_\text{atom}}\).

In the specific case of NaCl with its basis described in the face-centred lattice, the dynamical matrix is a \(24\times24\) matrix and the eigenvalue problem can be solved relatively quickly on a modern computer. Still, it may be advantageous to minimise the work performed by Euphonic and let brille do the heavy lifting.

The interface between brille and Euphonic is implemented as a class brilleu.BrillEu in the brilleu module.

The first step in creating an appropriate interpolation grid object is to load the force constant data from, e.g., a CASTEP binary file. An appropriate file for NaCl can be used in conjunction with this script to reproduce the content of this tutorial.

nacl = getBrillEuObj('NaCl', max_volume=1e-5, parallel=True, sort=True, scattering_lengths=scattering_lengths)

Constructing the BrillEu object creates a hybrid grid in the irreducible Brillouin zone and then uses the ForceConstants object to store \(\omega_i(\mathbf{q}_\text{ir})\) and \(\hat{\epsilon}_{ij}(\mathbf{q}_\text{ir})\) at each grid vertex. The stored data are then used to determine a sorting permutation for each vertex which attempts to place equivalent modes at equal indexes on all neighbouring vertices. As we will see, this process is not yet perfected.

Dispersion plots

Interpolation can be used to determine the dispersion relation along any path in reciprocal space. One such path starts at the origin and goes to \(\mathbf{Q} = (1\,2\,3)\) which is on the Brillouin zone boundary for \(\mathbf{G} = (1\,1\,3)\) and passes through the \(\mathbf{G} = (1\,1\,1)\) and \(\mathbf{G} = (0\,2\,2)\) first Brillouin zones.

The \((0\,0\,0)\), \((1\,1\,1)\), \((0\,2\,2)\), and \((1\,1\,3)\) first Brillouin zones and the path from \((0\,0\,0)\) to \((1\,2\,3)\)

We can construct points along the path and find \(\omega_i(\mathbf{Q})\) using

xi = np.linspace(0, 1, 100)
qhkl = np.vstack((xi, 2 * xi, 3 * xi)).T
w_q = nacl.w_q(qhkl)

Since this path passes through four first Brillouin zones, an equivalent path must pass through the irreducible Brillouin zone multiple times. We can find the irreducible path using

irqhkl, G, rot, invrot = nacl.grid.BrillouinZone.ir_moveinto(qhkl)

The irreducible path is shown in the irreducible Brillouin zone on the left, note that it passes through the irreducible polyhedron multiple times. The dispersion relation on the right is obtained by linear interpolation within the irreducible polyhedron.

Inelastic Neutron Scattering

Inelastic neutron scattering experiments performed with surveying instruments produce observations of scattering intensity at points in \((\mathbf{Q},E)\) which are not on a regular grid.

The intensity measured in a neutron scattering experiment from one-phonon scattering

\[\begin{aligned} S(\mathbf{Q}, E) & \propto \sum_i S_i(\mathbf{Q}) \delta(E-\hbar\omega_i(\mathbf{Q})) \\ S_i(\mathbf{Q}) & = \frac{1}{\omega_i(\mathbf{Q})} \left|\sum_j \frac{b_j}{\sqrt{M_j}} \mathbf{Q}\cdot\hat{\epsilon}_{ij}(\mathbf{Q}) e^{i \mathbf{Q}\cdot\mathbf{r}^0_j} \right|^2 \end{aligned}\]

with \(S_i(\mathbf{Q})\) and \(\omega_i(\mathbf{Q})\) directly calculable from theory via, e.g., StructureFactor.

If the observations were on a regular \((\mathbf{Q},E)\) grid then it would be straightforward to calculate \(S_i(\mathbf{Q})\) and \(\omega_i(\mathbf{Q})\) for each unique \(\mathbf{Q}\) with the extension to \(S(\mathbf{Q},E)\) involving an approximation to the convolution of the instrumental resolution and a delta function.

Since, the detectors observe along (nonlinear) kinematically constrained paths through \((\mathbf{Q},E)\), each observation is likely to have a unique \(\mathbf{Q}\) and, as a result, the calculation of \(S_i(\mathbf{Q})\) and \(\omega_i(\mathbf{Q})\) is necessary for every observation.

If the \((\mathbf{Q},E)\) observation positions are known for all detectors then brille can be used to find \(S(\mathbf{Q},E)\). Correct handling of the non-uniform observations is non-trivial and special software packages exist for this purpose, e.g., Horace. Here we will avoid these complications by calculating \(S(\mathbf{Q},E)\) for all points in a \((\mathbf{Q},E)\) grid.

mu, nu = np.mgrid[0:1:50j, 0.5:30:60j]
ξ = mu.reshape(mu.size, 1)
En = nu.reshape(nu.size, 1)
Qhkl = np.concatenate((ξ, 2 * ξ, 3 * ξ), axis=1)

brSQE = nacl(Qhkl, En, resfun='sho', param=1., T=1.).reshape(mu.shape)

Inleastic neutron scattering intensity for the same path through recprocal space as above.