Linear Response Theory & Electric Susceptibility

Background: Kubo’s Formula

Since we are going to apply time-dependent perturbations in the linear response regime, we need to remind ourselves, how the values of an observable vary when the perturbation comes into play. That is to say we need to remind ourselves of:

Derivation of Kubo formula

{\displaystyle \langle {\hat {A}}(t)\rangle =\left\langle {\hat {A}}\right\rangle _{0}-{\frac {i}{\hbar }}\int _{t_{0}}^{t}dt'\left\langle \left[{\hat {A}}(t),{\hat {V}}{\mathord {\left(t'\right)}}\right]\right\rangle _{0}}

And thus in linear response theory, the change in expectation value of an operator \hat A is \delta\langle \hat A(t)\rangle = \frac{i}{\hbar}\int_{-\infty}^t dt' \, \langle [\hat H_I(t'), \hat A(t)] \rangle_0 , \tag{1}

with \hat H_I(t') the interaction Hamiltonian.

In our case we have a perturbation due to the interaction between electric polarization and the time dependent electric field.

\hat H(t) = H_0 - \int d\mathbf r'\,\hat{\mathbf P}(\mathbf r') \cdot \mathbf E(\mathbf r',t), \tag{2}

where \hat{\mathbf P}(\mathbf r) is the polarization operator.

For \hat A(t) = \hat P_i(\mathbf r,t),\ Equation 1 becomes: \delta \langle \hat P_i(\mathbf r,t)\rangle = \frac{i}{\hbar}\int_{-\infty}^t dt' \int d\mathbf r' \, \Big\langle \big[ -\hat P_j(\mathbf r',t') E_j(\mathbf r',t'), \, \hat P_i(\mathbf r,t) \big] \Big\rangle_0 .

Since E_j(\mathbf r',t') is a classical field (a c-number, not an operator), it commutes with operators and can be pulled out: = -\frac{i}{\hbar}\int_{-\infty}^t dt' \int d\mathbf r' \, E_j(\mathbf r',t') \, \langle [\hat P_j(\mathbf r',t'), \hat P_i(\mathbf r,t)] \rangle_0 .

Using [X,Y] = -[Y,X], this becomes \delta \langle \hat P_i(\mathbf r,t)\rangle = \frac{i}{\hbar}\int_{-\infty}^t dt' \int d\mathbf r' \, \langle [\hat P_i(\mathbf r,t), \hat P_j(\mathbf r',t')] \rangle_0 \, E_j(\mathbf r',t') .

We can now introduce the retarded susceptibility (response function): \chi_{ij}^R(\mathbf r,t;\mathbf r',t') = \frac{i}{\hbar} \, \Theta(t-t') \, \langle [\hat P_i(\mathbf r,t), \hat P_j(\mathbf r',t')] \rangle_0 .

Then the response takes the compact Green’s-function form: \delta \langle \hat P_i(\mathbf r,t)\rangle = \int d\mathbf r' \int_{-\infty}^\infty dt' \, \chi_{ij}^R(\mathbf r,t;\mathbf r',t') \, E_j(\mathbf r',t').

From Kubo’s Formula to Susceptibility

We start from the general constitutive relation between the polarization P and the electric field E in a linear medium:

P_i(\mathbf{r},t) = \varepsilon_0 \int d^3r' \int_{-\infty}^{\infty} dt'\; \chi(\mathbf{r}, t; \mathbf{r}',t') \, E(\mathbf{r}',t'). \tag{3}

This Green’s function relation emerges out of the Kubo formula. Below we make two changes to it: (1) We use the component notation with the summation convention in place, and (2) We impose the constraint that time dependence is stationary (i.e. is a function of t-t') and the space dependence exhibits translational invariance or homongeneity (is a function of \mathbf{r}-\mathbf{r'}):

P_i(\mathbf{r},t) = \varepsilon_0 \int d^3r' \int_{-\infty}^{\infty} dt'\; \chi_{ij}(\mathbf{r}-\mathbf{r}',\,t-t') \, E_j(\mathbf{r}',t').

So we have gone from the Kubo formula, to a space-time convolution because the susceptibility depends only on the differences \mathbf{r}-\mathbf{r}' and t-t'. The Fourier transform is going to turn it into a multiplication:

\boxed{\tilde{P}_i(\mathbf{k},\omega)=\varepsilon_0\,\tilde\chi_{ij}(\mathbf{k},\omega)\,\tilde{E}_j(\mathbf{k},\omega)}.

Fourier Transform Conventions

We define the Fourier transform as

f(\mathbf{r},t) = \int \frac{d^3k}{(2\pi)^3} \int \frac{d\omega}{2\pi} \; \tilde{f}(\mathbf{k},\omega) e^{i(\mathbf{k}\cdot\mathbf{r}-\omega t)},

with the inverse transform

\tilde{f}(\mathbf{k},\omega) = \int d^3r \int dt\; f(\mathbf{r},t) e^{-i(\mathbf{k}\cdot\mathbf{r}-\omega t)}.

The Fourier transform P_i(\mathbf{r},t):

\tilde{P}_i(\mathbf{k},\omega) = \int d^3r \int dt\; e^{-i(\mathbf{k}\cdot\mathbf{r}-\omega t)} P_i(\mathbf{r},t).

Inserting the constitutive relation gives

\tilde{P}_i(\mathbf{k},\omega) = \varepsilon_0 \int d^3r \int dt\; e^{-i(\mathbf{k}\cdot\mathbf{r}-\omega t)} \int d^3r' dt' \, \chi_{ij}(\mathbf{r}-\mathbf{r}',t-t') E_j(\mathbf{r}',t').

Let \mathbf{R}=\mathbf{r}-\mathbf{r}' and \tau=t-t', so \mathbf{r}=\mathbf{R}+\mathbf{r}' and t=\tau+t'. Then

\tilde{P}_i(\mathbf{k},\omega) = \varepsilon_0 \int d^3R d\tau \int d^3r' dt'\; e^{-i[\mathbf{k}\cdot(\mathbf{R}+\mathbf{r}') -\omega (\tau+t')]} \chi_{ij}(\mathbf{R},\tau) E_j(\mathbf{r}',t').

The exponential factor splits into two parts:

e^{-i(\mathbf{k}\cdot\mathbf{R}-\omega \tau)}\, e^{-i(\mathbf{k}\cdot\mathbf{r}'-\omega t')}.

The integrals factor, giving

\tilde{P}_i(\mathbf{k},\omega) = \varepsilon_0 \underbrace{\int d^3R d\tau\;\chi_{ij}(\mathbf{R},\tau) e^{-i(\mathbf{k}\cdot\mathbf{R}-\omega\tau)}}_{\tilde\chi_{ij}(\mathbf{k},\omega)} \underbrace{\int d^3r' dt'\; E_j(\mathbf{r}',t') e^{-i(\mathbf{k}\cdot\mathbf{r}'-\omega t')}}_{\tilde{E}_j(\mathbf{k},\omega)}.

So we obtain the diagonal relation:

\boxed{\tilde{P}_i(\mathbf{k},\omega)=\varepsilon_0\,\tilde\chi_{ij}(\mathbf{k},\omega)\,\tilde{E}_j(\mathbf{k},\omega)}.

In doing so we have taken the following path: \chi(\mathbf{r}, t; \mathbf{r'}, t') \rightarrow \chi(\mathbf{r} - \mathbf{r'}, t-t')\rightarrow \chi(s, \tau)\rightarrow \chi(\mathbf{k}, \omega)

There is another way of seeing how \chi goes over to the Fourier domain. We transform ALL the variables independently:

\chi(\mathbf{r}, t; \mathbf{r'}, t') \rightarrow \chi(\mathbf{k}, \mathbf{k'}, \mathbf{\omega}, \mathbf{\omega'})

But the imposition of the stationarity and homogeneity conditions introduces delta functions:

\chi(\mathbf{r} - \mathbf{r'}, t-t') \rightarrow \chi(\mathbf{k}, \mathbf{k'}, \mathbf{\omega}, \mathbf{\omega'}) \delta(\mathbf{k}- \mathbf{k'})\delta(\mathbf{\omega}- \mathbf{\omega'})

To see this we start with the most general case (before assuming translational invariance). \chi depends on two points and two times:

\chi_{ij}(\mathbf{r},\mathbf{r}',t,t').

We could Fourier transform it using two independent sets of variables:

\chi_{ij}(\mathbf{r},\mathbf{r}',t,t') = \int \frac{d^3k\,d^3k'}{(2\pi)^6} \frac{d\omega\,d\omega'}{(2\pi)^2}\; \tilde\chi_{ij}(\mathbf{k},\mathbf{k}',\omega,\omega') e^{i(\mathbf{k}\cdot\mathbf{r}-\omega t)}e^{-i(\mathbf{k}'\cdot\mathbf{r}'-\omega' t')}.

But if the medium is homogeneous in space and stationary in time, then

\chi_{ij}(\mathbf{r},\mathbf{r}',t,t') = \chi_{ij}(\mathbf{r}-\mathbf{r}',t-t').

Letting \mathbf{R}=\mathbf{r}-\mathbf{r}' and \tau=t-t', we have

\tilde\chi_{ij}(\mathbf{k},\mathbf{k}',\omega,\omega') =\int d^3r\, d^3r'\, dt\, dt' \; \chi_{ij}(\mathbf{r}-\mathbf{r}',t-t') e^{-i(\mathbf{k}\cdot\mathbf{r}-\omega t)} e^{i(\mathbf{k}'\cdot\mathbf{r}'-\omega' t')}.

Changing variables to (\mathbf{R},\mathbf{r}',\tau,t') gives

\int d^3R\, d^3r'\, d\tau\, dt' \; \chi_{ij}(\mathbf{R},\tau) e^{-i(\mathbf{k}\cdot\mathbf{R}-\omega\tau)} e^{-i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}'} e^{i(\omega-\omega')t'}.

Now do the \mathbf{r}' and t' integrals:

\int d^3r'\; e^{-i(\mathbf{k}-\mathbf{k}')\cdot\mathbf{r}'} =(2\pi)^3\delta(\mathbf{k}-\mathbf{k}')

\int dt'\; e^{i(\omega-\omega')t'} =2\pi\delta(\omega-\omega').

So we find

\boxed{\tilde\chi_{ij}(\mathbf{k},\mathbf{k}',\omega,\omega') =(2\pi)^4 \delta(\mathbf{k}-\mathbf{k}')\delta(\omega-\omega')\, \chi_{ij}(\mathbf{k},\omega)}.

Physically

The delta functions are a direct mathematical consequence of space and time translational invariance.
They enforce momentum and energy conservation: a plane wave at (\mathbf{k},\omega) can only drive a polarization at the same (\mathbf{k},\omega).
This is why the constitutive relation becomes diagonal in Fourier space:

So we know that Equation 3 would become the product of Fourier transformed versions of \chi and \mathbf{E}. We can then substitute them from above to find once again: P_i(\mathbf{k},\omega)=\varepsilon_0\,\chi_{ij}(\mathbf{k},\omega)\,E_j(\mathbf{k},\omega).

Now to find susceptibility \chi(\mathbf{k}, \omega)

We now want to find an expression for components of \chi(\mathbf{k}, \omega). If we follow Kubo linear-response formula step by step we will see that resonance terms (with denominators \hbar\omega_{\ell\ell'}\mp\hbar\omega\mp i\gamma) appear.

So now the task is to:

Find \chi(s, \tau)

Step 1: Kubo formula for the polarization response

Suppose the system is perturbed by a classical electric field that couples to the polarization by

\color{green} H'(t) = -\int d^3r'\; P_j(\mathbf r')\,E_j(\mathbf r',t).

To first order in the field, the induced expectation value of the polarization is (Kubo linear response)

\delta\langle P_i(\mathbf r,t)\rangle = \frac{i}{\hbar}\int_{-\infty}^{t} dt'\, \big\langle \big[\color{green}H'(t')\color{black},\,P_i(\mathbf r,t)\big]\big\rangle_0,

where \langle\cdots\rangle_0 is the expectation value in the unperturbed equilibrium state. Inserting H'(t') and using linearity gives the response kernel

\delta\langle P_i(\mathbf r,t)\rangle = \int d^3r'\int_{-\infty}^t dt'\; \chi_{ij}(\mathbf r-\mathbf r',t-t')\,E_j(\mathbf r',t'),

with the susceptibility (Kubo form):

\boxed{\; \chi_{ij}(\mathbf s,\tau) = \frac{i}{\hbar}\,\Theta(\tau)\,\langle\,[\,P_i(\mathbf r+\mathbf s,\;t+\tau),\;P_j(\mathbf r,t)\,]\,\rangle_0 \;}

(here \mathbf s=\mathbf r-\mathbf r', \tau=t-t'; because of homogeneity the kernel depends only on differences).

\newcommand{\green}[1]{\color{green}{#1}} \newcommand{\red}[1]{\color{red}{#1}} \newcommand{\gray}[1]{\color{gray}{#1}}

2) Move to the energy eigenbasis and Heisenberg time dependence

Let \{|\ell\rangle\} be eigenstates of the unperturbed Hamiltonian H with energies E_\ell. In the Heisenberg picture

P_i(\mathbf r,t) = e^{\frac{i}{\hbar}Ht}\,P_i(\mathbf r)\,e^{-\frac{i}{\hbar}Ht}.

Insert a complete set of eigenstates for the commutator and evaluate the expectation value as a thermal (or general diagonal) ensemble with probabilities w_\ell (e.g. w_\ell=e^{-\beta E_\ell}/Z for thermal equilibrium). Then

\chi_{ij}(\mathbf s,\tau) = \frac{i}{\hbar}\Theta(\tau)\sum_{\ell} w_\ell \red{\sum_{m}}\Big(\langle\ell|P_i(\mathbf r+\mathbf s,t+\tau)\red{|m\rangle\langle m|}P_j(\mathbf r,t)|\ell\rangle \\ - \langle\ell|P_j(\mathbf r,t)\red{|m\rangle\langle m|}P_i(\mathbf r+\mathbf s,t+\tau)|\ell\rangle\Big).

\begin{align*} \chi_{ij}(\mathbf s,\tau) = \frac{i}{\hbar}\Theta(\tau)\sum_{\ell} w_\ell \sum_{m}\Big( \color{#F4442E}\langle\ell|P_i(\mathbf r+\mathbf s,t+\tau)|m\rangle\color{black} \color{#188FA7}\langle m|P_j(\mathbf r,t)|\ell\rangle \color{black}\\ - \langle\ell|P_j(\mathbf r,t)|m\rangle\langle m|P_i(\mathbf r+\mathbf s,t+\tau)|\ell\rangle\Big). \end{align*} \tag{4}

Where the colored bits can be reduced to their Schrodinger counterparts by using the Heisenberg time dependence:

\begin{align*} \color{#F4442E}\langle\ell|P_i(\mathbf r+\mathbf s,t+\tau)|m\rangle\color{black} &= e^{i\color{green}\tfrac{(E_\ell-E_m)}{\hbar}\color{black} (t+\tau)}\langle\ell|P_i(\mathbf r+\mathbf s)|m\rangle \\ &=e^{i\color{green}\omega_{\ell m}\color{black}(t+\tau)}\langle\ell|P_i(\mathbf r+\mathbf s)|m\rangle, \end{align*}

\begin{align*} \color{#188FA7}\langle m|P_j(\mathbf r,t)|\ell\rangle\color{black} &= e^{i\color{green}\tfrac{(E_\ell-E_m)}{\hbar}\color{black}t}\langle m|P_j(\mathbf r)|\ell\rangle\\ &= e^{i\color{green}\omega_{\ell m}\color{black}(t)}\langle m|P_j(\mathbf r)|\ell\rangle \end{align*}

Where we have introduced the Bohr frequencies

\omega_{m\ell}\equiv\frac{E_m-E_\ell}{\hbar} = -\omega_{\ell m}.

Note how the overall t-dependence cancels (as it must for a time-translation invariant equilibrium average), leaving \tau-dependence through e^{\tfrac{i}{\hbar}(E_\ell-E_m)\tau}.

So the commutator contribution becomes (for \tau>0)

\begin{align*} \chi_{ij}(\mathbf s,\tau) = \frac{i}{\hbar}\Theta(\tau)\sum_{\ell,m} w_\ell\Big( &e^{i\omega_{\ell m}\tau} \langle\ell|P_i(\mathbf r+\mathbf s)|m\rangle\langle m|P_j(\mathbf r)|\ell\rangle\\ - &e^{i\omega_{m\ell}\tau}\langle\ell|P_j(\mathbf r)|m\rangle\langle m|P_i(\mathbf r+\mathbf s)|\ell\rangle \Big). \end{align*}

Convert it to \chi(\mathbf{k}, \omega) using the recipes above.

Now form the space–time Fourier transform used in the text (including the adiabatic factor e^{-\gamma\tau}, \gamma\to0^+, to ensure convergence):

\chi_{ij}(\mathbf k,\omega) = \lim_{\gamma\to0^+}\int d^3s\int_0^\infty d\tau\; \chi_{ij}(\mathbf s,\tau)\,e^{-i\mathbf k\cdot\mathbf s + i\omega \tau -\gamma\tau}.

Insert the expression for \chi_{ij}(\mathbf s,\tau). Define the spatial Fourier components of the polarization operator

P_{i,\mathbf k} \equiv V^{-1/2}\int d^3r\;P_i(\mathbf r)\,e^{-i\mathbf k\cdot\mathbf r},

so that matrix elements of the spatially shifted operator pick up factors e^{\pm i\mathbf k\cdot\mathbf s}. After doing the \mathbf s-integral one obtains matrix elements of the Fourier components P_{i,\mathbf k} and P_{j,-\mathbf k} (this is the reason the final expression involves P_{i,\mathbf k} and P_{j,-\mathbf k}).

Carrying out the \tau-integral for a generic exponential factor:

\int_0^\infty d\tau\; e^{i\omega_{\ell m}\tau}\,e^{i\omega\tau}\,e^{-\gamma\tau} = \int_0^\infty d\tau\; e^{i(\omega_{\ell m}+\omega)\tau-\gamma\tau} = \frac{i}{\omega_{\ell m}+\omega + i\gamma}.

Be careful with signs: if the integrand was e^{i\omega_{\ell m}\tau}\,e^{ - i\omega \tau} then the denominator would be \omega_{\ell m}-\omega + i\gamma. Using the combination that arises from the commutator gives the two resonance denominators.

After performing the spatial Fourier transform (producing matrix elements (P_{i,\mathbf k})_{\ell m} and (P_{j,-\mathbf k})_{m\ell}) and the \tau-integrals, and reorganizing the two terms from the commutator, you obtain

\boxed{\; \chi_{ij}(\mathbf k,\omega) = \sum_{\ell,m} w_\ell\Big[ \frac{(P_{i,\mathbf k})_{\ell m}\,(P_{j,-\mathbf k})_{m\ell}}{\hbar\omega_{m\ell}-\hbar\omega - i\hbar\gamma} \;+\; \frac{(P_{j,-\mathbf k})_{\ell m}\,(P_{i,\mathbf k})_{m\ell}}{\hbar\omega_{m\ell}+\hbar\omega + i\hbar\gamma} \Big], \;}

where \omega_{m\ell}=(E_m-E_\ell)/\hbar and the infinitesimal \gamma\to0^+ implements the adiabatic switching (and gives the usual causal i0^+ prescription).

This is exactly the structure of equation (6.1.31): two contributions appear corresponding to absorption-like (denominator \hbar\omega_{m\ell}-\hbar\omega - i0^+) and emission-like (denominator \hbar\omega_{m\ell}+\hbar\omega + i0^+) processes, each weighted by the equilibrium population w_\ell and the relevant transition matrix elements of the polarization operator in \mathbf k-space.

Remarks / interpretation

The first fraction resonates when \hbar\omega\approx\hbar\omega_{m\ell}=E_m-E_\ell (absorption from \ell to m).
The second fraction resonates when \hbar\omega\approx -\hbar\omega_{m\ell} (the opposite process).
The small positive \gamma enforces causality (the Fourier transform of \Theta(\tau) gives poles in the lower/upper half-plane consistent with retarded response).
If w_\ell is thermal, you can combine terms using detailed-balance relations to express \chi in alternative forms.