We can solve for linear response to various different perturbations. The most straight-forward case is the response of a finite system to an electric field \(\mathcal{E}_{i,\omega}\) with frequency \(\omega\) in the direction \(i\), where the perturbation operator is \(\delta \hat v = \hat r_i\)~\cite{Andrade_2007}. In this case the polarizability can be calculated as
%
\begin{equation}
\alpha_{ij}\left( \omega \right) = - \sum_n^{\rm occ} \left[ \langle \varphi_n \vert \hat{r}_i \vert \frac{\partial \varphi_n}{\partial \mathcal{E}_{j, \omega}} \rangle + \langle \frac{\partial \varphi_n}{\partial \mathcal{E}_{j, -\omega}} \vert \hat{r}_i \vert \varphi_n \rangle \right] \label{eq:sternheimer_polarizability}\ .
\end{equation}
%
The calculation of the polarizability yields optical response properties (that can be extended to nonlinear response)~\cite{Andrade_2007,Vila_2010} and, for imaginary frequencies, van der Waals coefficients~\cite{Botti_2008}.

It is also possible to use the formalism to compute vibrational properties for finite and periodic systems~\cite{Baroni_2001,Kadantsev_2005}. Currently Octopus implements the calculations of vibrations for finite systems. In this case the perturbation operator is an infinitesimal ionic displacement
$\partial \hat{H}/\partial R_{i \alpha} = \partial \hat{v}_{\alpha}/\partial R_{i \alpha}$, for each direction $i$ and atom $\alpha$.
The quantity to calculate is the dynamical matrix, or Hessian, given by
%
\begin{multline}
  D_{i \alpha, j \beta} = \frac{\partial^2 E}{\partial R_{i \alpha} \partial R_{j \beta}} = D_{i \alpha, j \beta}^{\rm ion-ion} \\
  - \sum_n^{\rm occ} \left[ \left< \varphi_n \left| \frac{\partial \hat{v}_{\alpha}}{\partial R_{i \alpha}} \right| \frac{\partial \varphi_n}{\partial R_{j \beta}} \right> + {\rm cc.}
      + \delta_{\alpha \beta} \left< \varphi_n \left| \frac{\partial^2 \hat{v}_{\alpha}}{\partial R_{i \alpha} \partial R_{\alpha j}} \right| \varphi_n \right> \right] \label{eq:dynmatrix}
\end{multline}
%
The contribution from the ion-ion interaction energy is
\begin{align}
    D_{i \alpha, j \beta}^{\rm ion-ion} =
    \begin{cases}
    Z_\alpha \sum_{\gamma \ne \alpha} Z_\gamma \left[ \frac{\delta_{ij}}{\left|R_\alpha - R_\gamma \right|^3} - 3 \frac{\left( R_{i \alpha} - R_{i \gamma} \right)}{\left|R_\alpha - R_\gamma \right|^4} \right] & \textrm{if}\ \alpha = \beta \\
    -Z_\alpha Z_\beta \left[ \frac{\delta_{ij}}{\left|R_\alpha - R_\beta \right|^3} - 3 \frac{\left( R_{i \alpha} - R_{i \beta} \right)}{\left|R_\alpha - R_\beta \right|^4} \right]              & \textrm{if}\ \alpha \ne \beta
    \end{cases}
\end{align}
where $Z_\alpha$ is the ionic charge of atom $\alpha$.
We have found that an alternative expression for the perturbation operator yields more accurate results when discretized. This is discussed in section~\ref{sec:forces}.

Vibrational frequencies $\omega$ are obtained by solving the eigenvalue equation
%
\begin{equation}
\frac{1}{\sqrt{m_\alpha m_\beta}} D_{i \alpha, j \beta} x_{j \beta} = -\omega^2 x_{j \beta}\ ,
\end{equation}
%
where \(m_\alpha\) is the mass for ion \(\alpha\). For a finite system of $N$ atoms, there should be 3 zero-frequency translational modes and 3
zero-frequency rotational modes. However, they may appear at positive or imaginary frequencies,
due to the finite size of the simulation domain, the discretization of the grid, and finite precision
in solution of the ground state and Sternheimer equation. Improving convergence brings them closer to zero.

The Born effective charges can be computed from the response of the dipole moment to ionic displacement:
\begin{align} \label{eq:Born}
  Z^{*}_{i j \alpha} = -\frac{\partial^2 E}{\partial \mathcal{E}_i \partial R_{j \alpha}}
  = \frac{\partial \mu_i}{\partial R_{j \alpha}} = Z_{\alpha} \delta_{ij} - \sum_n^{\rm occ} \left< \varphi_n \left| \hat{r}_i \right| \frac{\partial \varphi_n}{\partial R_{j \alpha}} \right>\,.
\end{align}
The intensities for each mode for absorption of radiation polarized in direction $i$, which can be used to predict infrared spectra, are calculated by multiplying by the normal-mode eigenvector $x$
\begin{align}
I_i = \sum_{j \alpha} Z^{*}_{ij \alpha} x_{j \alpha}\ .
\end{align}

The Born charges must obey the acoustic sum rule, from translational invariance
\begin{align}
\sum_{\alpha} Z^{*}_{i j \alpha} = Q_{\rm tot} \delta_{ij}\ .
\end{align}
For each $ij$, we enforce this sum rule by distributing the discrepancy equally among the atoms, and thus obtaining corrected Born charges:
\begin{equation}
\tilde{Z}^{*}_{i j \alpha} = Z^{*}_{i j \alpha} + \frac1N\left( Q_{\rm tot} \delta_{ij} - \sum_{\alpha} Z^{*}_{i j \alpha} \right)\ .
\end{equation}
The discrepancy arises from the same causes as the non-zero translational and rotational modes. 

By mixing the response to ionic displacements and electric perturbations it is possible to calculate vibrational Raman response coefficients~\cite{PhysRevB.66.100301}. This feature however, it is still not implemented in Octopus.

The Sternheimer equation can be used in conjunction with $\vec{k} \cdot \vec{p}$ perturbation theory~\cite{Cardona_1966} to obtain
band velocities and effective masses, as well as to apply electric fields via the quantum theory of polarization. In this case the perturbation is a displacement in the \(k\)-point. Using the effective Hamiltonian for the \(k\)-point
\begin{equation}
\hat{H}_{\vec{k}} = \mathrm{e}^{-i \vec{k} \cdot \vec{r}} \hat{H} \mathrm{e}^{i \vec{k} \cdot \vec{r}}\ ,
\end{equation}
the perturbation is represented by the operator
%
\begin{equation}
  \frac{\partial \hat{H}_{\vec{k}}}{\partial \vec{k}} = -i \nabla + \vec{k} + \left[\hat{v}_{\alpha}, \vec{r} \right]\ ,
\end{equation}
%
including the effect on the non-local pseudopotentials.
The first-order term gives the band group velocities in a periodic system,
\begin{equation}
  \mathrm{v}_{n\vec{k}} = \frac{\partial \epsilon_{n\vec{k}}}{\partial \vec{n}} 
  = \left< \varphi_{n\vec{k}} \left| \frac{\partial \hat H_{\vec{k}}}{\partial \vec{k}} \right| \varphi_{n\vec{k}} \right>\ .
\end{equation}
%adding in occupied contributions

Inverse effective mass tensors can be calculated by solving the Sternheimer equation for the \(\vec{k}\cdot\vec{p}\) perturbation. The equation is not solved self-consistently, since the variation of \(k\)-point is not a physical perturbation to the system; a converged \(k\)-grid should give the same density even if displaced slightly. The perturbation $\partial \hat{H}_{\vec{k}}/\partial \vec{k}$ is purely anti-Hermitian. We use instead $-i \partial \hat{H}_{\vec{k}}/\partial \vec{k}$ to obtain a Hermitian perturbation, which allows the response to real wavefunctions to remain real.
The effective mass tensors are calculated as follows:
%
\begin{multline}
  m^{-1}_{ijn\vec{k}} = \frac{\partial^2 \epsilon_{n\vec{k}}}{\partial k_i \partial k_j}
  = \delta_{ij} + \left< \varphi_{n\vec{k}} \left| \frac{\partial \hat{H}_{\vec{k}}}{\partial k_i} \right| \frac{\partial \varphi_{n\vec{k}}}{\partial k_j} \right> + \mathrm{cc.}
  \\
  + \left< \varphi_{n\vec{k}} \left| \left[ \hat{r}_i, \left[ \hat{r}_j, \hat{v}_{\alpha} \right] \right] \right| \varphi_{n\vec{k}} \right>\ .
\end{multline}

The $\vec{k} \cdot \vec{p}$ wavefunctions can be used to compute the response to electric fields in periodic systems.
In finite systems, a homogeneous electric field can be represented simply via the position operator $\vec{r}$.
However, this operator is not well defined in a periodic system and cannot be used. According to the quantum theory of polarization,
the solution is to replace $\vec{r} \varphi$ with $-i \partial \varphi/\partial k$~\cite{Resta_2007}, and then use this as the perturbation on the right hand side in the Sternheimer equation \cite{Gonze1997}.
While this is typically done with a finite difference with respect to $\vec{k}$ \cite{Corso_1996,QE-2009}, we use an analytic derivative from a previous $\vec{k} \cdot \vec{p}$ Sternheimer calculation.
Using the results in eq.~(\ref{eq:sternheimer_polarizability}) gives a formula for the polarization of the crystal:
\begin{align}
  \alpha_{ij} \left( \omega \right) = i \sum_n^{\rm occ} \left[ \left. \left< \frac{\partial \varphi_{n\vec{k}}}{\partial k_i} \right| \frac{\partial \varphi_{n\vec{k}}}{\partial \mathcal{E}_{j, \omega}} \right> + \left. \left< \frac{\partial \varphi_{n\vec{k}}}{\partial \mathcal{E}_{j, -\omega}} \right| \frac{\partial \varphi_{n\vec{k}}}{\partial k_j} \right> \right]\ .
\end{align}
The polarizability is most usefully represented in a periodic system via the dielectric constant
\begin{align}
  \epsilon_{ij} = \delta_{ij} + \frac{4 \pi}{V} \alpha_{ij}\ ,
\end{align}
where $V$ is the volume of the unit cell. This scheme can also be extended to non-linear response.

We can compute the Born charges from the electric-field response in either finite or periodic systems (as a complementary approach to using the vibrational response):
\begin{multline}
  Z^{*}_{i j \alpha} = -\frac{\partial^2 E}{\partial \mathcal{E}_i \partial R_{j \alpha}} = \frac{\partial F_{j \alpha}}{\partial \mathcal{E}_i} \\
  = Z_\alpha \delta_{ij} - \sum_n^{\rm occ} \left[ \left< \varphi_n \left| \frac{\partial \hat{v}_{\alpha}}{\partial R_{i \alpha}} \right| \frac{\partial \varphi_n}{\partial \mathcal{E}_j} \right> + {\rm cc.} \right]
\end{multline}
This expression can be evaluated with the same approach as for the dynamical matrix elements, and is easily generalized to non-zero frequency too.
We can also make the previous expression eq. \ref{eq:Born} for Born charges from the vibrational perturbation usable in a periodic system with the replacement $\vec{r} \varphi \rightarrow -i \partial \varphi/\partial k$.

Unfortunately, the $\vec{k} \cdot \vec{p}$ perturbation is not usable to calculate the polarization \cite{Resta_2007}, and a sum over strings of \(k\)-points on a finer grid is required. We have implemented the special case of a $\Gamma$-point calculation for a large super-cell, where the single-point Berry phase can be used \cite{Yaschenko1998}. For cell sizes $L_i$ in each direction, the dipole moment is derived from the determinant of a matrix whose basis is the occupied KS orbitals:
\begin{align}
\mu_i = - \frac{L_i}{2 \pi} \mathrm{Im}\ {\rm ln}\ {\rm det}\ \left< \varphi_n \left| e^{- 2 \pi i x_i/L_i} \right| \varphi_m \right>\ .
\end{align}

% magnetic is non-self-consistent if wfns are real.