Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘phy487’

Fourier coefficient integral for periodic function

Posted by peeterjoot on November 5, 2013

[Click here for a PDF of this post with nicer formatting]

In phy487 we’ve been using the fact that a periodic function

\begin{aligned}V(\mathbf{r}) = V(\mathbf{r} + \mathbf{r}_n),\end{aligned} \hspace{\stretch{1}}(1.1)

where

\begin{aligned}\mathbf{r}_n = a_1 \mathbf{a}_1 + a_2 \mathbf{a}_2 + a_3 \mathbf{a}_3,\end{aligned} \hspace{\stretch{1}}(1.2)

has a Fourier representation

\begin{aligned}V(\mathbf{r}) = \sum_\mathbf{G} V_\mathbf{G} e^{ i \mathbf{G} \cdot \mathbf{r} }.\end{aligned} \hspace{\stretch{1}}(1.3)

Here \mathbf{G} is a vector in reciprocal space, say

\begin{aligned}\mathbf{G}_{rst} = r \mathbf{g}_1 + s \mathbf{g}_2 + t \mathbf{g}_3,\end{aligned} \hspace{\stretch{1}}(1.4)

where

\begin{aligned}\mathbf{g}_i \cdot \mathbf{a}_j = 2 \pi \delta_{ij}.\end{aligned} \hspace{\stretch{1}}(1.5)

Now let’s express the explicit form for the Fourier coefficient V_\mathbf{G} so that we can compute the Fourier representation for some periodic potentials for some numerical experimentation. In particular, let’s think about what it meant to integrate over a unit cell. Suppose we have a parameterization of the points in the unit cell

\begin{aligned}\mathbf{r} = u \mathbf{a}_1 + v \mathbf{a}_2 + w \mathbf{a}_3,\end{aligned} \hspace{\stretch{1}}(1.6)

as sketched in fig. 1.1. Here u, v, w \in [0, 1]. We can compute the values of u, v, w for any vector \mathbf{r} in the cell by reciprocal projection

Fig 1.1: Unit cell

\begin{aligned}\mathbf{r} = \frac{1}{{2 \pi}} \left(  \left(  \mathbf{r} \cdot \mathbf{g}_1 \right) \mathbf{a}_1 + \left(  \mathbf{r} \cdot \mathbf{g}_2 \right) \mathbf{a}_2 + \left(  \mathbf{r} \cdot \mathbf{g}_3 \right) \mathbf{a}_3 \right),\end{aligned} \hspace{\stretch{1}}(1.7)

or

\begin{aligned}\begin{aligned}u(\mathbf{r}) &= \frac{1}{{2 \pi}} \mathbf{r} \cdot \mathbf{g}_1 \\ v(\mathbf{r}) &= \frac{1}{{2 \pi}} \mathbf{r} \cdot \mathbf{g}_2 \\ w(\mathbf{r}) &= \frac{1}{{2 \pi}} \mathbf{r} \cdot \mathbf{g}_3.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.8)

Let’s suppose that \mathbf{V}(\mathbf{r}) is period in the unit cell spanned by \mathbf{r} = u \mathbf{a}_1 + v \mathbf{a}_2 + w \mathbf{a}_3 with u, v, w \in [0, 1], and integrate over the unit cube for that parameterization to compute V_\mathbf{G}

\begin{aligned}\int_0^1 du\int_0^1 dv\int_0^1 dwV( u \mathbf{a}_1 + v \mathbf{a}_2 + w \mathbf{a}_3 ) e^{-i \mathbf{G}' \cdot \mathbf{r} }=\sum_{r s t}V_{\mathbf{G}_{r s t}}\int_0^1 du\int_0^1 dv\int_0^1 dwe^{-i \mathbf{G}' \cdot \mathbf{r} }e^{i \mathbf{G} \cdot \mathbf{r} }\end{aligned} \hspace{\stretch{1}}(1.9)

Let’s write

\begin{aligned}\begin{aligned}\mathbf{G} &= r \mathbf{g}_1 + s \mathbf{g}_2 + t \mathbf{g}_3 \\ \mathbf{G} &= r' \mathbf{g}_1 + s' \mathbf{g}_2 + t' \mathbf{g}_3,\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.10)

so that

\begin{aligned}e^{-i \mathbf{G}' \cdot \mathbf{r} } e^{i \mathbf{G} \cdot \mathbf{r} }=e^{ 2 \pi i (r - r') u } e^{ 2 \pi i (s - s') u } e^{ 2 \pi i (t - t') u } \end{aligned} \hspace{\stretch{1}}(1.11)

Picking the u integral of this integrand as representative, we have when r = r'

\begin{aligned}\int_0^1 du e^{ 2 \pi i (r - r') u } =\int_0^1 du= 1,\end{aligned} \hspace{\stretch{1}}(1.12)

and when r \ne r'

\begin{aligned}\int_0^1 du e^{ 2 \pi i (r - r') u } ={\left.{{   \frac{    e^{ 2 \pi i (r - r') u }   }   {   2 \pi i (r - r')    }}}\right\vert}_{{u = 0}}^{{1}}=\frac{1}{{2 \pi i (r - r') }} \left(  e^{ 2 \pi i (r - r') } - 1  \right).\end{aligned} \hspace{\stretch{1}}(1.13)

This is just zero since r - r' is an integer, so we have

\begin{aligned}\int_0^1 du e^{ 2 \pi i (r - r') u } = \delta_{r, r'}.\end{aligned} \hspace{\stretch{1}}(1.14)

This gives us

\begin{aligned}\int_0^1 du\int_0^1 dv\int_0^1 dwV( u \mathbf{a}_1 + v \mathbf{a}_2 + w \mathbf{a}_3 ) e^{ -2 \pi i r' u } e^{ -2 \pi i s' v } e^{ -2 \pi i t' w } =\sum_{r s t}V_{\mathbf{G}_{r s t}}\delta_{r s t, r' s' t'}= V_{\mathbf{G}_{r' s' t'}}.\end{aligned} \hspace{\stretch{1}}(1.15)

This is our \textAndIndex{Fourier coefficient}. The \textAndIndex{Fourier series} written out in gory but explicit detail is

\begin{aligned}\boxed{V( u \mathbf{a}_1 + v \mathbf{a}_2 + w \mathbf{a}_3 ) = \sum_{r s t}\left(  \int_0^1 du' \int_0^1 dv' \int_0^1 dw' V( u' \mathbf{a}_1 + v' \mathbf{a}_2 + w' \mathbf{a}_3 ) e^{ -2 \pi i (r u' + s v' + t w') }  \right)e^{ 2 \pi i (r u + s v + t w) }.}\end{aligned} \hspace{\stretch{1}}(1.16)

Also observe the unfortunate detail that we require integrability of the potential in the unit cell for the Fourier integrals to converge. This prohibits the use of the most obvious potential for numerical experimentation, the inverse radial V(\mathbf{r}) = -1/\left\lvert {\mathbf{r}} \right\rvert.

Posted in Math and Physics Learning. | Tagged: , , , , | Leave a Comment »

Huygens diffraction

Posted by peeterjoot on October 23, 2013

[Click here for a PDF of this post with nicer formatting]

We were presented with a diffraction result, that the intensity can be expressed as the Fourier transform of the aperture. Let’s review a derivation of that based on the Huygens principle. Consider the aperture of fig. 2.1.

Fig 1.1: Diffraction aperture

The Huygens principle expresses the amplitude of a wave U(\mathbf{r}) in terms of it’s amplitude U_\circ at \mathbf{r} = 0 as

\begin{aligned}U(\mathbf{r}) \propto \frac{U_\circ e^{i k \left\lvert {\mathbf{r}} \right\rvert}}{\left\lvert {\mathbf{r}} \right\rvert}.\end{aligned} \hspace{\stretch{1}}(2.1)

For multiple point diffraction, the diffracted wave is a superposition of such contributions from all points in the aperture. For example, two exponentials are summed when considering a two slit diffraction apparatus. For a more general aperture as above, we’d form

\begin{aligned}U(\mathbf{R}') \propto U_\circ \int_{\mathrm{A}} d^2\mathbf{r} \frac{e^{i k \left\lvert {\mathbf{r} - \mathbf{R}} \right\rvert}}{\left\lvert {\mathbf{r} - \mathbf{R}} \right\rvert} \frac{e^{i k \left\lvert {\mathbf{R}' - \mathbf{r}} \right\rvert}}{\left\lvert {\mathbf{R}' - \mathbf{r}} \right\rvert}.\end{aligned} \hspace{\stretch{1}}(2.2)

Note that this the Huygens result is an approximation in many ways. Fresnel later fixed up the proportionality factor and considered the angular dependence of the between the source and diffracted rays (the obliquity factor). That corrected result is known as the Huygens-Fresnel principle. Kirchhoff later considered solutions to the scalar wave equations for the electromagnetic field components \mathbf{E} and \mathbf{B}, ignoring the Maxwell coupling of these fields. See section 8.3.1 [1], section A.2 section 10.4 of [3], or section 5.2 of [2] for details. See section 9.8 [4] for a vector diffraction treatment.

Let’s proceed with using eq. 2.2 to obtain our result from class. For simplicity, first position the origin in the aperture itself as in fig. 2.2.

Fig 1.2: Diffraction aperture with origin in aperture

Now we are set to perform a first order expansion of the vector lengths \left\lvert {\mathbf{r} - \mathbf{R}} \right\rvert and \left\lvert {\mathbf{r} - \mathbf{R}'} \right\rvert. It’s sufficient to consider just one of these, expanding to first order

\begin{aligned}\left\lvert {\mathbf{r} - \mathbf{R}} \right\rvert=\sqrt{(\mathbf{r} - \mathbf{R})^2} &= \sqrt{\mathbf{R}^2 + \mathbf{r}^2 - 2 \mathbf{r} \cdot \mathbf{R}} \\ &= R\sqrt{1 + \frac{\mathbf{r}^2}{\mathbf{R}^2} - 2 \mathbf{r} \cdot \frac{\mathbf{R}}{\mathbf{R}^2}} \\ &= \approx R \left( 1 + \frac{1}{{2}} \frac{\mathbf{r}^2}{\mathbf{R}^2} - \mathbf{r} \cdot \frac{\mathbf{R}}{\mathbf{R}^2} \right) \\ &= R + \frac{1}{{2}} \frac{\mathbf{r}^2}{R} - \mathbf{r} \cdot \frac{\mathbf{R}}{R}.\end{aligned} \hspace{\stretch{1}}(2.3)

Assume that both \mathbf{R} and \mathbf{R}' to be far enough from the aperture that we can approximate the \left\lvert {\mathbf{R} - \mathbf{r}} \right\rvert and \left\lvert {\mathbf{R}' - \mathbf{r}} \right\rvert terms downstairs as R = \left\lvert {\mathbf{R}} \right\rvert and R' = \left\lvert {\mathbf{R}'} \right\rvert respectively. Additionally, ignore the second order term above, significant for Fresnel diffraction where the diffraction pattern close to the aperture is examined, but not in the far field. This gives

\begin{aligned}U(\mathbf{R}') &\sim \frac{U_\circ }{R R'}\int_{\mathrm{A}} d^2\mathbf{r} e^{i k \left( R - \mathbf{r} \cdot \hat{\mathbf{R}} \right)} e^{i k \left( R' - \mathbf{r} \cdot \hat{\mathbf{R}}' \right)} \\ &= \frac{U_\circ }{R R'} e^{i k (R + R')}\int_{\mathrm{A}} d^2\mathbf{r} e^{-i k \mathbf{r} \cdot \hat{\mathbf{R}}} e^{-i k \mathbf{r} \cdot \hat{\mathbf{R}}'}.\end{aligned} \hspace{\stretch{1}}(2.4)

Finally write

\begin{aligned}\begin{aligned}\mathbf{k}_s &= - k\hat{\mathbf{R}} \\ \mathbf{k} &= k\hat{\mathbf{R}}',\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.5)

for

\begin{aligned}U(\mathbf{R}') \sim\frac{U_\circ }{R R'}e^{i k (R + R')}\int_{\mathrm{A}} d^2\mathbf{r} e^{-i \mathbf{r} \cdot (\mathbf{k} - \mathbf{k}_s)}.\end{aligned} \hspace{\stretch{1}}(2.6)

Finally, write

\begin{aligned}\mathbf{K} = \mathbf{k} - \mathbf{k}_s,\end{aligned} \hspace{\stretch{1}}(2.7)

and introduce an aperture function \rho(\mathbf{r}) (we called this the scattering density) that is unity for any regions of the aperture that light is allowed through, and zero when light is blocked, and some value in [0,1] for translucent regions of the aperture. We can now expand the integral over the surface containing the aperture

\begin{aligned}\boxed{U(\mathbf{R}') \sim\frac{U_\circ }{R R'} e^{i k (R + R')}\int\rho(\mathbf{r}) e^{-i \mathbf{r} \cdot \mathbf{K}}.}\end{aligned} \hspace{\stretch{1}}(2.8)

References

[1] M. Born and E. Wolf. Principles of optics: electromagnetic theory of propagation, interference and diffraction of light. Cambridge university press, 1980.

[2] G.R. Fowles. Introduction to modern optics. Dover Pubns, 1989.

[3] E. Hecht. Optics. 1998.

[4] JD Jackson. Classical Electrodynamics Wiley. John Wiley and Sons, 2nd edition, 1975.

Posted in Math and Physics Learning. | Tagged: , | Leave a Comment »

Exponential solutions to second order linear system

Posted by peeterjoot on October 20, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

Motivation

We’re discussing specific forms to systems of coupled linear differential equations, such as a loop of “spring” connected masses (i.e. atoms interacting with harmonic oscillator potentials) as sketched in fig. 1.1.

Fig 1.1: Three springs loop

 

Instead of assuming a solution, let’s see how far we can get attacking this problem systematically.

Matrix methods

Suppose that we have a set of N masses constrained to a circle interacting with harmonic potentials. The Lagrangian for such a system (using modulo N indexing) is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} \sum_{k = 0}^{2} m_k \dot{u}_k^2 - \frac{1}{{2}} \sum_{k = 0}^2 \kappa_k \left( u_{k+1} - u_k \right)^2.\end{aligned} \hspace{\stretch{1}}(1.1)

The force equations follow directly from the Euler-Lagrange equations

\begin{aligned}0 = \frac{d{{}}}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{u}_{n, \alpha}}}- \frac{\partial {\mathcal{L}}}{\partial {u_{n, \alpha}}}.\end{aligned} \hspace{\stretch{1}}(1.2)

For the simple three particle system depicted above, this is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_0 \dot{u}_0^2 +\frac{1}{{2}} m_1 \dot{u}_1^2 +\frac{1}{{2}} m_2 \dot{u}_2^2 - \frac{1}{{2}} \kappa_0 \left( u_1 - u_0 \right)^2- \frac{1}{{2}} \kappa_1 \left( u_2 - u_1 \right)^2- \frac{1}{{2}} \kappa_2 \left( u_0 - u_2 \right)^2,\end{aligned} \hspace{\stretch{1}}(1.3)

with equations of motion

\begin{aligned}\begin{aligned}0 &= m_0 \dot{d}{u}_0 + \kappa_0 \left( u_0 - u_1 \right) + \kappa_2 \left( u_0 - u_2 \right) \\ 0 &= m_1 \dot{d}{u}_1 + \kappa_1 \left( u_1 - u_2 \right) + \kappa_0 \left( u_1 - u_0 \right) \\ 0 &= m_2 \dot{d}{u}_2 + \kappa_2 \left( u_2 - u_0 \right) + \kappa_1 \left( u_2 - u_1 \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.4)

Let’s partially non-dimensionalize this. First introduce average mass \bar{{m}} and spring constants \bar{{\kappa}}, and rearrange slightly

\begin{aligned}\begin{aligned}\frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_0 &= -\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} \left( u_0 - u_1 \right) - \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} \left( u_0 - u_2 \right) \\ \frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_1 &= -\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} \left( u_1 - u_2 \right) - \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} \left( u_1 - u_0 \right) \\ \frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_2 &= -\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} \left( u_2 - u_0 \right) - \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} \left( u_2 - u_1 \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.5)

With

\begin{aligned}\tau = \sqrt{\frac{\bar{{k}}}{\bar{{m}}}} t = \Omega t\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

\begin{aligned}\mathbf{u} = \begin{bmatrix}u_0 \\ u_1 \\ u_2\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

\begin{aligned}B = \begin{bmatrix}-\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} - \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} &\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} &\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} \\ \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} &-\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} - \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} &\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} \\ \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} & \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} &-\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} - \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

Our system takes the form

\begin{aligned}\frac{d^2 \mathbf{u}}{d\tau^2} = B \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

We can at least theoretically solve this in a simple fashion if we first convert it to a first order system. We can do that by augmenting our vector of displacements with their first derivatives

\begin{aligned}\mathbf{w} =\begin{bmatrix}\mathbf{u} \\ \frac{d \mathbf{u}}{d\tau} \end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

So that

\begin{aligned}\frac{d \mathbf{w}}{d\tau} =\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} \mathbf{w}= A \mathbf{w}.\end{aligned} \hspace{\stretch{1}}(1.0.9)

Now the solution is conceptually trivial

\begin{aligned}\mathbf{w} = e^{A \tau} \mathbf{w}_\circ.\end{aligned} \hspace{\stretch{1}}(1.0.9)

We are however, faced with the task of exponentiating the matrix A. All the powers of this matrix A will be required, but they turn out to be easy to calculate

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^2=\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}B & 0 \\ 0 & B\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(1.0.11a)

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^3=\begin{bmatrix}B & 0 \\ 0 & B\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}0 & B \\ B^2 & 0\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(1.0.11b)

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^4=\begin{bmatrix}0 & B \\ B^2 & 0\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}B^2 & 0 \\ 0 & B^2\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

allowing us to write out the matrix exponential

\begin{aligned}e^{A \tau} = \sum_{k = 0}^\infty \frac{\tau^{2k}}{(2k)!} \begin{bmatrix}B^k & 0 \\ 0 & B^k\end{bmatrix}+\sum_{k = 0}^\infty \frac{\tau^{2k + 1}}{(2k + 1)!} \begin{bmatrix}0 & B^k \\ B^{k+1} & 0\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

Case I: No zero eigenvalues

Provided that B has no zero eigenvalues, we could factor this as

\begin{aligned}\begin{bmatrix}0 & B^k \\ B^{k+1} & 0\end{bmatrix}=\begin{bmatrix}0 & B^{-1/2} \\ B^{1/2} & 0\end{bmatrix}\begin{bmatrix}B^{k + 1/2} & 0 \\ 0 & B^{k + 1/2}\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

This initially leads us to believe the following, but we’ll find out that the three springs interaction matrix B does have a zero eigenvalue, and we’ll have to be more careful. If there were any such interaction matrices that did not have such a zero we could simply write

\begin{aligned}e^{A \tau} = \sum_{k = 0}^\infty \frac{\tau^{2k}}{(2k)!} \begin{bmatrix}\sqrt{B}^{2k} & 0 \\ 0 & \sqrt{B}^{2k}\end{bmatrix}+\begin{bmatrix}0 & B^{-1/2} \\ B^{1/2} & 0\end{bmatrix}\sum_{k = 0}^\infty \frac{\tau^{2k + 1}}{(2k + 1)!} \begin{bmatrix}\sqrt{B}^{2 k + 1} & 0 \\ 0 & \sqrt{B}^{2 k+1} \end{bmatrix}=\cosh\begin{bmatrix}\sqrt{B} \tau & 0 \\ 0 & \sqrt{B} \tau\end{bmatrix}+ \begin{bmatrix}0 & 1/\sqrt{B} \tau \\ \sqrt{B} \tau & 0\end{bmatrix}\sinh\begin{bmatrix}\sqrt{B} \tau & 0 \\ 0 & \sqrt{B} \tau\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

This is

\begin{aligned}e^{A \tau}=\begin{bmatrix}\cosh \sqrt{B} \tau & (1/\sqrt{B}) \sinh \sqrt{B} \tau \\ \sqrt{B} \sinh \sqrt{B} \tau & \cosh \sqrt{B} \tau \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

The solution, written out is

\begin{aligned}\begin{bmatrix}\mathbf{u} \\ \mathbf{u}'\end{bmatrix}=\begin{bmatrix}\cosh \sqrt{B} \tau & (1/\sqrt{B}) \sinh \sqrt{B} \tau \\ \sqrt{B} \sinh \sqrt{B} \tau & \cosh \sqrt{B} \tau \end{bmatrix}\begin{bmatrix}\mathbf{u}_\circ \\ \mathbf{u}_\circ'\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

so that

\begin{aligned}\boxed{\mathbf{u} = \cosh \sqrt{B} \tau \mathbf{u}_\circ + \frac{1}{{\sqrt{B}}} \sinh \sqrt{B} \tau \mathbf{u}_\circ'.}\end{aligned} \hspace{\stretch{1}}(1.0.11c)

As a check differentiation twice shows that this is in fact the general solution, since we have

\begin{aligned}\mathbf{u}' = \sqrt{B} \sinh \sqrt{B} \tau \mathbf{u}_\circ + \cosh \sqrt{B} \tau \mathbf{u}_\circ',\end{aligned} \hspace{\stretch{1}}(1.0.11c)

and

\begin{aligned}\mathbf{u}'' = B \cosh \sqrt{B} \tau \mathbf{u}_\circ + \sqrt{B} \sinh \sqrt{B} \tau \mathbf{u}_\circ'= B \left( \cosh \sqrt{B} \tau \mathbf{u}_\circ + \frac{1}{{\sqrt{B}}} \sinh \sqrt{B} \tau \mathbf{u}_\circ' \right)= B \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

Observe that this solution is a general solution to second order constant coefficient linear systems of the form we have in eq. 1.5. However, to make it meaningful we do have the additional computational task of performing an eigensystem decomposition of the matrix B. We expect negative eigenvalues that will give us oscillatory solutions (ie: the matrix square roots will have imaginary eigenvalues).

Example: An example diagonalization to try things out

{example:threeSpringLoop:1}{

Let’s do that diagonalization for the simplest of the three springs system as an example, with \kappa_j = \bar{{k}} and m_j = \bar{{m}}, so that we have

\begin{aligned}B = \begin{bmatrix}-2 & 1 & 1 \\ 1 & -2 & 1 \\ 1 & 1 & -2 \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.20)

A orthonormal eigensystem for B is

\begin{aligned}\left\{\mathbf{e}_{-3, 1}, \mathbf{e}_{-3, 2}, \mathbf{e}_{0, 1} \right\}=\left\{\frac{1}{{\sqrt{6}}}\begin{bmatrix} -1 \\ -1 \\ 2 \end{bmatrix},\frac{1}{{\sqrt{2}}}\begin{bmatrix} -1 \\ 1 \\ 0 \end{bmatrix},\frac{1}{{\sqrt{3}}}\begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}\right\}.\end{aligned} \hspace{\stretch{1}}(1.21)

With

\begin{aligned}\begin{aligned}U &= \frac{1}{{\sqrt{6}}}\begin{bmatrix} -\sqrt{3} & -1 & \sqrt{2} \\ 0 & 2 & \sqrt{2} \\ \sqrt{3} & -1 & \sqrt{2} \end{bmatrix} \\ D &= \begin{bmatrix}-3 & 0 & 0 \\ 0 & -3 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix}\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.22)

We have

\begin{aligned}B = U D U^\text{T},\end{aligned} \hspace{\stretch{1}}(1.23)

We also find that B and its root are intimately related in a surprising way

\begin{aligned}\sqrt{B} = \sqrt{3} iU \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix}U^\text{T}=\frac{1}{{\sqrt{3} i}} B.\end{aligned} \hspace{\stretch{1}}(1.24)

We also see, unfortunately that B has a zero eigenvalue, so we can’t compute 1/\sqrt{B}. We’ll have to back and up and start again differently.
}

Case II: allowing for zero eigenvalues

Now that we realize we have to deal with zero eigenvalues, a different approach suggests itself. Instead of reducing our system using a Hamiltonian transformation to a first order system, let’s utilize that diagonalization directly. Our system is

\begin{aligned}\mathbf{u}'' = B \mathbf{u} = U D U^{-1} \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.25)

where D = [ \lambda_i \delta_{ij} ] and

\begin{aligned}\left( U^{-1} \mathbf{u} \right)'' = D \left( U^{-1} \mathbf{u} \right).\end{aligned} \hspace{\stretch{1}}(1.26)

Let

\begin{aligned}\mathbf{z} = U^{-1} \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.27)

so that our system is just

\begin{aligned}\mathbf{z}'' = D \mathbf{z},\end{aligned} \hspace{\stretch{1}}(1.28)

or

\begin{aligned}z_i'' = \lambda_i z_i.\end{aligned} \hspace{\stretch{1}}(1.29)

This is N equations, each decoupled and solvable by inspection. Suppose we group the eigenvalues into sets \{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0 \}. Our solution is then

\begin{aligned}\mathbf{z} = \sum_{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0}\left( a_n \cos \sqrt{-\lambda_n} \tau + b_n \sin \sqrt{-\lambda_n} \tau \right)\mathbf{e}_n+\left( a_p \cosh \sqrt{\lambda_p} \tau + b_p \sinh \sqrt{\lambda_p} \tau \right)\mathbf{e}_p+\left( a_z + b_z \tau \right) \mathbf{e}_z.\end{aligned} \hspace{\stretch{1}}(1.30)

Transforming back to lattice coordinates using \mathbf{u} = U \mathbf{z}, we have

\begin{aligned}\mathbf{u} = \sum_{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0}\left( a_n \cos \sqrt{-\lambda_n} \tau + b_n \sin \sqrt{-\lambda_n} \tau \right)U \mathbf{e}_n+\left( a_p \cosh \sqrt{\lambda_p} \tau + b_p \sinh \sqrt{\lambda_p} \tau \right)U \mathbf{e}_p+\left( a_z + b_z \tau \right) U \mathbf{e}_z.\end{aligned} \hspace{\stretch{1}}(1.31)

We see that the zero eigenvalues integration terms have no contribution to the lattice coordinates, since U \mathbf{e}_z = \lambda_z \mathbf{e}_z = 0, for all \lambda_z = 0.

If U = [ \mathbf{e}_i ] are a set of not necessarily orthonormal eigenvectors for B, then the vectors \mathbf{f}_i, where \mathbf{e}_i \cdot \mathbf{f}_j = \delta_{ij} are the reciprocal frame vectors. These can be extracted from U^{-1} = [ \mathbf{f}_i ]^\text{T} (i.e., the rows of U^{-1}). Taking dot products between \mathbf{f}_i with \mathbf{u}(0) = \mathbf{u}_\circ and \mathbf{u}'(0) = \mathbf{u}'_\circ, provides us with the unknown coefficients a_n, b_n

\begin{aligned}\mathbf{u}(\tau)= \sum_{ \lambda_n < 0, \lambda_p > 0 }\left( (\mathbf{u}_\circ \cdot \mathbf{f}_n) \cos \sqrt{-\lambda_n} \tau + \frac{\mathbf{u}_\circ' \cdot \mathbf{f}_n}{\sqrt{-\lambda_n} } \sin \sqrt{-\lambda_n} \tau \right)\mathbf{e}_n+\left( (\mathbf{u}_\circ \cdot \mathbf{f}_p) \cosh \sqrt{\lambda_p} \tau + \frac{\mathbf{u}_\circ' \cdot \mathbf{f}_p}{\sqrt{-\lambda_p} } \sinh \sqrt{\lambda_p} \tau \right)\mathbf{e}_p.\end{aligned} \hspace{\stretch{1}}(1.32)

Supposing that we constrain ourself to looking at just the oscillatory solutions (i.e. the lattice does not shake itself to pieces), then we have

\begin{aligned}\boxed{\mathbf{u}(\tau)= \sum_{ \lambda_n < 0 }\left( \sum_j \mathbf{e}_{n,j} \mathbf{f}_{n,j}^\text{T} \right)\left( \mathbf{u}_\circ \cos \sqrt{-\lambda_n} \tau + \frac{\mathbf{u}_\circ' }{\sqrt{-\lambda_n} } \sin \sqrt{-\lambda_n} \tau \right).}\end{aligned} \hspace{\stretch{1}}(1.33)

Eigenvectors for eigenvalues that were degenerate have been explicitly enumerated here, something previously implied. Observe that the dot products of the form (\mathbf{a} \cdot \mathbf{f}_i) \mathbf{e}_i have been put into projector operator form to group terms more nicely. The solution can be thought of as a weighted projector operator working as a time evolution operator from the initial state.

Example: Our example interaction revisited

{example:threeSpringLoop:2}{

Recall that we had an orthonormal basis for the \lambda = -3 eigensubspace for the interaction example of eq. 1.20 again, so \mathbf{e}_{-3,i} = \mathbf{f}_{-3, i}. We can sum \mathbf{e}_{-3,1} \mathbf{e}_{-3, 1}^\text{T} + \mathbf{e}_{-3,2} \mathbf{e}_{-3, 2}^\text{T} to find

\begin{aligned}\mathbf{u}(\tau)= \frac{1}{{3}}\begin{bmatrix}2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2\end{bmatrix}\left( \mathbf{u}_\circ \cos \sqrt{3} \tau + \frac{ \mathbf{u}_\circ' }{\sqrt{3} } \sin \sqrt{3} \tau \right).\end{aligned} \hspace{\stretch{1}}(1.34)

The leading matrix is an orthonormal projector of the initial conditions onto the eigen subspace for \lambda_n = -3. Observe that this is proportional to B itself, scaled by the square of the non-zero eigenvalue of \sqrt{B}. From this we can confirm by inspection that this is a solution to \mathbf{u}'' = B \mathbf{u}, as desired.
}

Fourier transform methods

Let’s now try another item from our usual toolbox on these sorts of second order systems, the Fourier transform. For a one variable function of time let’s write the transform pair as

\begin{aligned}x(t) = \int_{-\infty}^{\infty} \tilde{x}(\omega) e^{-i \omega t} d\omega\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

\begin{aligned}\tilde{x}(\omega) = \frac{1}{{2\pi}} \int_{-\infty}^{\infty} x(t) e^{i \omega t} dt\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

One mass harmonic oscillator

The simplest second order system is that of the harmonic oscillator

\begin{aligned}0 = \dot{d}{x}(t) + \omega_\circ^2 x.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Application of the transform gives

\begin{aligned}0 = \left( \frac{d^2}{dt^2} + \omega_\circ^2 \right)\int_{-\infty}^{\infty} \tilde{x}(\omega) e^{-i \omega t} d\omega= \int_{-\infty}^{\infty} \left( -\omega^2 + \omega_\circ^2 \right)\tilde{x}(\omega) e^{-i \omega t} d\omega.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

We clearly have a constraint that is a function of frequency, but one that has to hold for all time. Let’s transform this constraint to the frequency domain to consider that constraint independent of time.

\begin{aligned}0 =\frac{1}{{2 \pi}}\int_{-\infty}^{\infty} dt e^{ i \omega t}\int_{-\infty}^{\infty} \left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') e^{-i \omega' t} d\omega'=\int_{-\infty}^{\infty} d\omega'\left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') \frac{1}{{2 \pi}}\int_{-\infty}^{\infty} dt e^{ i (\omega -\omega') t}=\int_{-\infty}^{\infty} d\omega'\left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') \delta( \omega - \omega' )=\left( -{\omega}^2 + \omega_\circ^2 \right)\tilde{x}(\omega).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

How do we make sense of this?

Since \omega is an integration variable, we can’t just mandate that it equals the constant driving frequency \pm \omega_\circ. It’s clear that we require a constraint on the transform \tilde{x}(\omega) as well. As a trial solution, imagine that

\begin{aligned}\tilde{x}(\omega) = \left\{\begin{array}{l l}\tilde{x}_\circ & \quad \mbox{if latex \left\lvert {\omega – \pm \omega_\circ} \right\rvert < \omega_{\text{cutoff}}$} \\ 0 & \quad \mbox{otherwise}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)$

This gives us

\begin{aligned}0 = \tilde{x}_\circ\int_{\pm \omega_\circ -\omega_{\text{cutoff}}}^{\pm \omega_\circ +\omega_{\text{cutoff}}}\left( \omega^2 - \omega_\circ^2 \right)e^{-i \omega t} d\omega.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Now it is clear that we can satisfy our constraint only if the interval [\pm \omega_\circ -\omega_{\text{cutoff}}, \pm \omega_\circ + \omega_{\text{cutoff}}] is made infinitesimal. Specifically, we require both a \omega^2 = \omega_\circ^2 constraint and that the transform \tilde{x}(\omega) have a delta function nature. That is

\begin{aligned}\tilde{x}(\omega) = A \delta(\omega - \omega_\circ)+ B \delta(\omega + \omega_\circ).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Substitution back into our transform gives

\begin{aligned}x(t) = A e^{-i \omega_\circ t}+ B e^{i \omega_\circ t}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

We can verify quickly that this satisfies our harmonic equation \dot{d}{x} = -\omega_\circ x.

Two mass harmonic oscillator

Having applied the transform technique to the very simplest second order system, we can now consider the next more complex system, that of two harmonically interacting masses (i.e. two masses on a frictionless spring).

F1

Our system is described by

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_1 \dot{u}_1^2+ \frac{1}{{2}} m_2 \dot{u}_2^2-\frac{1}{{2}} \kappa \left( u_2 - u_1 \right)^2,\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

and the pair of Euler-Lagrange equations

\begin{aligned}0 = \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{u}_i}}- \frac{\partial {\mathcal{L}}}{\partial {u_i}}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

The equations of motion are

\begin{aligned}\begin{aligned}0 &= m_1 \dot{d}{u}_1 + \kappa \left( u_1 - u_2 \right) \\ 0 &= m_2 \dot{d}{u}_2 + \kappa \left( u_2 - u_1 \right) \end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Let

\begin{aligned}\begin{aligned}u_1(t) &= \int_{-\infty}^\infty \tilde{u}_1(\omega) e^{-i \omega t} d\omega \\ u_2(t) &= \int_{-\infty}^\infty \tilde{u}_2(\omega) e^{-i \omega t} d\omega.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Insertion of these transform pairs into our equations of motion produces a pair of simultaneous integral equations to solve

\begin{aligned}\begin{aligned}0 &= \int_{-\infty}^\infty \left( \left( -m_1 \omega^2 + \kappa \right) \tilde{u}_1(\omega) - \kappa \tilde{u}_2(\omega) \right)e^{-i \omega t} d\omega \\ 0 &= \int_{-\infty}^\infty \left( \left( -m_2 \omega^2 + \kappa \right) \tilde{u}_2(\omega) - \kappa \tilde{u}_1(\omega) \right)e^{-i \omega t} d\omega.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

As with the single spring case, we can decouple these equations with an inverse transformation operation \int e^{i \omega' t}/2\pi, which gives us (after dropping primes)

\begin{aligned}0 = \begin{bmatrix}\left( -m_1 \omega^2 + \kappa \right) &- \kappa \\ - \kappa &\left( -m_2 \omega^2 + \kappa \right) \end{bmatrix}\begin{bmatrix}\tilde{u}_1(\omega) \\ \tilde{u}_2(\omega)\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Taking determinants gives us the constraint on the frequency

\begin{aligned}0 = \left( -m_1 \omega^2 + \kappa \right) \left( -m_2 \omega^2 + \kappa \right) - \kappa^2=m_1 m_2 \omega^4-(m_1 + m_2) \omega^2 =\omega^2 \left( m_1 m_2 \omega^2 -\kappa (m_1 + m_2) \right).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Introducing a reduced mass

\begin{aligned}\frac{1}{{\mu}} = \frac{1}{m_1}+\frac{1}{m_2},\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

the pair of solutions are

\begin{aligned}\begin{aligned}\omega^2 &= 0 \\ \omega^2 &=\frac{\kappa}{\mu}\equiv \omega_\circ^2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

As with the single mass oscillator, we require the functions \tilde{u}{\omega} to also be expressed as delta functions. The frequency constraint and that delta function requirement together can be expressed, for j \in \{0, 1\} as

\begin{aligned}\tilde{u}_j(\omega) = A_{j+} \delta( \omega - \omega_\circ )+ A_{j0} \delta( \omega )+ A_{j-} \delta( \omega + \omega_\circ ).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

With a transformation back to time domain, we have functions of the form

\begin{aligned}\begin{aligned}u_1(t) &= A_{j+} e^{ -i \omega_\circ t }+ A_{j0} + A_{j-} e^{ i \omega_\circ t } \\ u_2(t) &= B_{j+} e^{ -i \omega_\circ t }+ B_{j0} + B_{j-} e^{ i \omega_\circ t }.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Back insertion of these into the equations of motion, we have

\begin{aligned}\begin{aligned}0 &=-m_1 \omega_\circ^2\left( A_{j+} e^{ -i \omega_\circ t } + A_{j-} e^{ i \omega_\circ t } \right)+ \kappa\left( \left( A_{j+} - B_{j+} \right) e^{ -i \omega_\circ t } + \left( A_{j-} - B_{j-} \right) e^{ i \omega_\circ t } + A_{j0} - B_{j0} \right) \\ 0 &=-m_2 \omega_\circ^2\left( B_{j+} e^{ -i \omega_\circ t } + B_{j-} e^{ i \omega_\circ t } \right)+ \kappa\left( \left( B_{j+} - A_{j+} \right) e^{ -i \omega_\circ t } + \left( B_{j-} - A_{j-} \right) e^{ i \omega_\circ t } + B_{j0} - A_{j0} \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Equality requires identity for all powers of e^{i \omega_\circ t}, or

\begin{aligned}\begin{aligned}0 &= B_{j0} - A_{j0} \\ 0 &= -m_1 \omega_\circ^2 A_{j+} + \kappa \left( A_{j+} - B_{j+} \right) \\ 0 &= -m_1 \omega_\circ^2 A_{j-} + \kappa \left( A_{j-} - B_{j-} \right) \\ 0 &= -m_2 \omega_\circ^2 B_{j+} + \kappa \left( B_{j+} - A_{j+} \right) \\ 0 &= -m_2 \omega_\circ^2 B_{j-} + \kappa \left( B_{j-} - A_{j-} \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

or B_{j0} = A_{j0} and

\begin{aligned}0 =\begin{bmatrix}\kappa -m_1 \omega_\circ^2 & 0 & - \kappa & 0 \\ 0 & \kappa -m_1 \omega_\circ^2 & 0 & - \kappa \\ -\kappa & 0 & \kappa -m_2 \omega_\circ^2 & 0 \\ 0 & -\kappa & 0 & \kappa -m_2 \omega_\circ^2 \end{bmatrix}\begin{bmatrix}A_{j+} \\ A_{j-} \\ B_{j+} \\ B_{j-} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Observe that

\begin{aligned}\kappa - m_1 \omega_\circ^2=\kappa - m_1 \kappa \left( \frac{1}{{m_1}} + \frac{1}{{m_2}} \right)=\kappa \left( 1 - 1 - \frac{m_1}{m_2} \right)= -\kappa \frac{m_1}{m_2},\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

(with a similar alternate result). We can rewrite eq. 1.0.35.35 as

\begin{aligned}0 =-\kappa\begin{bmatrix}m_1/m_2 & 0 & 1 & 0 \\ 0 & m_1/m_2 & 0 & 1 \\ 1 & 0 & m_2/m_1 & 0 \\ 0 & 1 & 0 & m_2/m_1 \end{bmatrix}\begin{bmatrix}A_{j+} \\ A_{j-} \\ B_{j+} \\ B_{j-} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

It’s clear that there’s two pairs of linear dependencies here, so the determinant is zero as expected. We can read off the remaining relations. Our undetermined coefficients are given by

\begin{aligned}\begin{aligned}B_{j0} &= A_{j0} \\ m_1 A_{j+} &= -m_2 B_{j+} \\ m_1 A_{j-} &= -m_2 B_{j-} \end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

\begin{aligned}\boxed{\begin{aligned}u_1(t) &= a+ m_2 b e^{ -i \omega_\circ t }+ m_2 c e^{ i \omega_\circ t } \\ u_2(t) &= a- m_1 b e^{ -i \omega_\circ t }- m_1 c e^{ i \omega_\circ t }.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Observe that the constant term is not really of interest, since it represents a constant displacement of both atoms (just a change of coordinates).

Check:

\begin{aligned}u_1(t) - u_2(t) = + (m_1 + m_2)b e^{ -i \omega_\circ t }+ (m_1 + m_2)c e^{ i \omega_\circ t },\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

\begin{aligned}m_1 \dot{d}{u}_1(t) + \kappa (u_1(t) - u_2(t) )=-m_1 m_2 \omega_\circ^2 \left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)+(m_1 + m_2) \kappa \left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)=\left( -m_1 m_2 \omega_\circ^2 + (m_1 + m_2) \kappa \right)\left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)=\left( -m_1 m_2 \kappa \frac{m_1 + m_2}{m_1 m_2} + (m_1 + m_2) \kappa \right)\left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)= 0.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Reflection

We’ve seen that we can solve any of these constant coefficient systems exactly using matrix methods, however, these will not be practical for large systems unless we have methods to solve for all the non-zero eigenvalues and their corresponding eigenvectors. With the Fourier transform methods we find that our solutions in the frequency domain is of the form

\begin{aligned}\tilde{u}_j(\omega) = \sum a_{kj} \delta( \omega - \omega_{kj} ),\end{aligned} \hspace{\stretch{1}}(1.63)

or in the time domain

\begin{aligned}u(t) = \sum a_{kj} e^{ - i \omega_{kj} t }.\end{aligned} \hspace{\stretch{1}}(1.64)

We assumed exactly this form of solution in class. The trial solution that we used in class factored out a phase shift from a_{kj} of the form of e^{ i q x_n }, but that doesn’t change the underling form of that assumed solution. We have however found a good justification for the trial solution we utilized.

Posted in Math and Physics Learning. | Tagged: , , , , , , , | Leave a Comment »

Discrete Fourier transform

Posted by peeterjoot on October 11, 2013

[Click here for a PDF of this post with nicer formatting]

For decoupling trial solutions of lattice vibrations we have what appears to be a need for the use of a discrete Fourier transform. This is described briefly in [1], but there’s no proof of the inversion formula. Let’s work through that detail.

Let’s start given a periodic signal of the following form

\begin{aligned}a_n = \sum_{k = 0}^{N-1} A_k e^{-2 \pi i k n/N },\end{aligned} \hspace{\stretch{1}}(1.1)

and assume that there’s an inversion formula of the form

\begin{aligned}A_k' \propto \sum_{n = 0}^{N-1} a_n e^{2 \pi i k' n/N }.\end{aligned} \hspace{\stretch{1}}(1.2)

Direct substitution should show if such a transform pair is valid, and determine the proportionality constant. Let’s try this

\begin{aligned}\sum_{n = 0}^{N-1} a_n e^{2 \pi i k' n/N }=\sum_{n = 0}^{N-1} e^{2 \pi i k' n/N }\sum_{k = 0}^{N-1} A_k e^{-2 \pi i k n/N }=\sum_{k = 0}^{N-1} A_k \sum_{n = 0}^{N-1} e^{2 \pi i (k' - k)n/N }.\end{aligned} \hspace{\stretch{1}}(1.3)

Observe that the interior sum is just N when k' = k. Let’s sum this for the values k \ne k', writing

\begin{aligned}S_N(k'-k) = \sum_{n = 0}^{N-1} e^{2 \pi i (k' - k)n/N }.\end{aligned} \hspace{\stretch{1}}(1.4)

With a = \exp( 2 \pi i (k' - k)/N) we have

\begin{aligned}S_N(k'-k) = \sum_{n = 0}^{N-1} a^n= \frac{a^N - 1}{a - 1}= \frac{a^{N/2}}{a^{1/2}} \frac{a^{N/2} - a^{-N/2}}{a^{1/2} - a^{-1/2}}= e^{ \pi i (k' - k)(1 - 1/N)} \frac{ \sin \pi (k'-k) }{\sin \pi (k' - k)/N}.\end{aligned} \hspace{\stretch{1}}(1.5)

Observe that k' - k \in [-N+1, N-1], and necessarily takes on just integer values. We have terms of the form \sin \pi m, for integer m in the numerator, always zero. In the denominator, the sine argument is in the range

\begin{aligned}[ \pi \left( -1 + \frac{1}{{N}} \right), \pi \left( 1 - \frac{1}{{N}} \right)],\end{aligned} \hspace{\stretch{1}}(1.6)

We can visualize that range as all the points on a sine curve with the integer multiples of \pi omitted, as in fig. 1.1.

Fig 1.1: Sine with integer multiples of pi omitted

Clearly the denominator is always non-zero when k \ne k'. This provides us with the desired inverse transformation relationship

\begin{aligned}\sum_{n = 0}^{N-1} a_n e^{2 \pi i k' n/N }= N \sum_{k' = 0}^{N-1} A_k' \delta_{k, k'}= N A_k'.\end{aligned} \hspace{\stretch{1}}(1.7)

Summary

Now that we know the relationship between the discrete set of values a_n and this discrete transformation of those values, let’s write the transform pair in a form that explicitly expresses this relationship.

\begin{aligned}\boxed{\begin{aligned}a_n &= \sum_{k = 0}^{N-1} \tilde{a}_k e^{-2 \pi i k n/N } \\ \tilde{a}_k &= \frac{1}{{N}} \sum_{n = 0}^{N-1} a_n e^{2 \pi i k n/N }.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(1.8)

We have also shown that our discrete sum of exponentials has an delta function operator nature, a fact that will likely be handy in various manipulations.

\begin{aligned}\boxed{\sum_{n = 0}^{N-1} e^{2 \pi i (k' - k)n/N } = N \delta_{k, k'}.}\end{aligned} \hspace{\stretch{1}}(1.9)

References

[1] S.S. Haykin. Communication systems. Wiley, 1994.

Posted in Math and Physics Learning. | Tagged: , , | Leave a Comment »

First post of phy487 (Condensed Matter Physics) notes posted.

Posted by peeterjoot on October 4, 2013

I’ve posted a first (incomplete) draft of my lecture notes for the Winter 2013, University of Toronto Condensed Matter Physics course (PHY487H1F), taught by Prof. Stephen Julian.

Official course description:

“Introduction to the concepts used in the modern treatment of solids. The student is assumed to be familiar with elementary quantum mechanics. Topics include: bonding in solids, crystal structures, lattice vibrations, free electron model of metals, band structure, thermal properties, magnetism and superconductivity (time permitting).”

In this course there will be an attempt to cover

  • Bonding and crystal structure: Why solids form, describing periodic order mathematically, diffraction.
  • Lattice vibrations: .elementary excitation. of a periodic array of atoms (periodicity allows a 10 23 body problem to be solved).
  • Electrons in solids: Explaining this introduces a need for quantum mechanics.  Periodicity -> band structure -> insulators vs metals.
  • Electrical conduction in solids. (metals and semiconductors).

This document contains:

  • Plain old lecture notes. These mirror what was covered in class, possibly augmented with additional details.
  • Personal notes exploring details that were not clear to me from the lectures, or from the texts associated with the lecture material.
  • Assigned problems. Like anything else take these as is. I may or may not have gone back and corrected errors, and did not see the graded versions of the last two problem sets.
  • Some worked problems attempted as course prep, for fun, or for test preparation, or post test reflection.
  • Links to Mathematica workbooks associated with these notes.

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , | 4 Comments »