\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{a4margins,hyperref,amssymb,amsmath}

\newcommand{\comment}[1]{}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\sq}[1]{\left(#1\right)^2}

\newcommand{\dx}{\,\mathrm d x}
\newcommand{\dy}{\,\mathrm d y}
\newcommand{\dz}{\,\mathrm d z}
\newcommand{\dV}{\,\mathrm d V}

\newcommand{\mtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\detmtwo}[4]{\left|\begin{array}{cc}#1&#2\\#3&#4\end{array}\right|}

\newcommand{\ang}[1]{\langle#1\rangle}

\begin{document}

\title{Wave Equations with the Exponential of a Quadratic as a Solution}
\author{Stephen Brooks}
\maketitle

\section{Motivation}\label{-motivation}
Numerical simulations are typically done with points (macroparticles or grid meshes), or infinite plane waves (Fourier space) as a basis.  However, for some applications involving wave equations, something in between is needed, which is infinite in neither frequency space nor real space.  This motivates the use of wavelets or wavepackets.  In this note, a simple wavepacket that can be displaced in both frequency and position space (while being localised in both) is described.  By letting the parameters of the wavepacket vary in time, it is already an exact solution of many non-interacting wave equations.

\section{Wavepackets in Wave Equations}\label{-definition}
\newcommand{\re}{\mathrm{Re}\,}
\newcommand{\im}{\mathrm{Im}\,}
In one dimension ($x$), consider functions of the form $e^{ax^2+bx+c}$, where $a,b,c$ may be complex.  Typically $\re a\ge 0$ so the function does not tend to infinity.  The centre of the wavepacket is at $x=-\re b/(2\re a)$.  The spatial frequency $k=(2\im a)x+\im b$, where $a$ can be chosen to be real if no chirp is required.  The overall constant multiplier of $e^c$ may be used for normalisation purposes.

If the parameters vary as a function of time, functions of the following form are obtained:
\[ f(x,t) = e^{a(t)x^2+b(t)x+c(t)}. \]
Using prime for $\partial_x$ and dot for $\partial_t$, derivatives of $f$ have a fairly simple form:
\begin{eqnarray*}
\dot f &=& (\dot a x^2+\dot b x+\dot c)f \\
f' &=& (2ax+b)f \\
f'' &=& (2a+(2ax+b)^2)f = (4a^2x^2+4abx+(2a+b^2))f.\\
\end{eqnarray*}
Consider differential equations of the general form below, which is first-order in time:
\[ \dot f = Af''+(B_1x+B_0)f'+(C_2x^2+C_1x+C_0)f. \]  
Substituting the above expressions for the derivatives and dividing throughout by $f$ gives:
\[ \dot a x^2+\dot b x+\dot c = A(4a^2x^2+4abx+(2a+b^2)) + (B_1x+B_0)(2ax+b) + C_2x^2+C_1x+C_0, \]
which must be true for all $x$, so coefficients of powers of $x$ may be equated to give
\begin{eqnarray*}
\dot a &=& 4Aa^2+2B_1a + C_2 \\
\dot b &=& 4Aab + B_1b + 2B_0a + C_1 \\
\dot c &=& A(2a+b^2) + B_0b + C_0.
\end{eqnarray*}
The equations above determine the time evolution of the parameters $a,b,c$ for $f(x,t)$ to be a solution of the differential equation; they are mildly nonlinear but are easily integrated numerically on a computer.  Simulating the evolution of three numbers takes much less computational effort than a whole mesh of numbers, so if the state is well approximated by the sum of a few wavepackets (the original equation was linear in $f$), then this is a good basis to use.

The above derivation also explains why certain length polynomials (such as $B_1x+B_0$) were allowed as coefficients in the differential equation: if the highest order term of the derivative of $f$ in question is below $x^2$, additional factors of $x$ are allowed.  The zero-order derivative allows a quadratic function to be multiplied by $f$, such as a quadratic potential used in quantum theory.

\subsection{Example: Time-Dependent Schr\"odinger Equation in 1D}
The wavefunction $\psi(x,t)$ of a single, non-relativistic quantum (scalar) particle in 1D satisfies the equation below:
\[ i\hbar\dot\psi = \left(\frac{-\hbar^2}{2m}\partial_x^2+V(x,t)\right)\psi, \]
which can be rearranged into the form
\[ \dot\psi = \frac{i\hbar}{2m}\psi''+\frac{-i}{\hbar}V\psi. \]
This is an example of the general form in the previous section if 
\[ A=\frac{i\hbar}{2m}, \qquad B_1=B_0=0, \qquad C_2x^2+C_1x+C_0=\frac{-i}{\hbar}V. \]
Note that this is only an exact solution for potentials $V$ that are quadratic in $x$.  Behaviour in other potentials could be approximated by using small wavepackets and approximating $V$ by its local second-order Taylor series.  However, time-varying potentials are allowed, since having the $C_n$ vary with time does not invalidate the derivation of $\dot a, \dot b, \dot c$.

The equation for the evolution of $a$ is
\[ \dot a = \frac{2i\hbar}{m}a^2 + \frac{-i}{\hbar}V_2, \]
where $V_2$ is the coefficient of $x^2$ in the potential.

\subsection{Higher-order Generalisation}
The wavepacket definition may be extended so that the exponent is a polynomial of order $N$, in which case the method for calculating the time derivatives of the coefficients is analogous.  The expression for $f^{(n)}$ contains an $n$th order polynomial, so its coefficient in the differential equation can be an $(N-n)$th order polynomial.

\subsection{Generalisation to Multiple Dimensions}\label{-multid}
\newcommand{\x}{\mathbf x}
\newcommand{\av}{\mathbf a}
\newcommand{\bv}{\mathbf b}
\newcommand{\cv}{\mathbf c}
\newcommand{\dv}{\mathbf d}
\newcommand{\sv}{\mathbf s}
\newcommand{\pv}{\mathbf p}
\newcommand{\qv}{\mathbf q}
\newcommand{\rv}{\mathbf r}
\newcommand{\tr}{\,\mathrm{tr}}
\newcommand{\sym}{\,\mathrm{sym}}
In multiple dimensions, where $\x$ is a vector, wavepackets can be defined in the following way:
\[ \exp\left(\sum_{i,j} a_{ij}x_ix_j + \sum_i b_i x_i + c\right) = e^{\x^TA\x + \bv\cdot\x + c}. \]
Here $A=(a_{ij})$ is a matrix, which can be chosen to be symmetric, while $\bv=(b_i)$ is a vector.  Defining
\[ f(\x,t) = e^{\x^TA(t)\x + \bv(t)\cdot\x + c(t)} \]
gives analogous expressions for the partial derivatives of $f$:
\begin{eqnarray*}
\dot f &=& (\x^T\dot A\x + \dot\bv\cdot\x + \dot c)f \\
\partial_i f &=& \left(2\sum_{j} a_{ij}x_j + b_i\right)f = (2\av_i\cdot\x + b_i)f \\
\partial_i \partial_j f &=& \partial_i ((2\av_j\cdot\x + b_j)f) \\
&=& (2a_{ij}+(2\av_i\cdot\x + b_i)(2\av_j\cdot\x + b_j))f \\
&=& (\x^T(4\av_i\av_j^T)\x + (2b_i\av_j+2b_j\av_i)\cdot\x + (2a_{ij}+b_ib_j))f.
\end{eqnarray*}

The general differential equation in 3D to the same order as the one used in 1D is:
\begin{eqnarray*}
\dot f &=& \sum_{i,j}p_{0,ij}\partial_i\partial_j f
+ \sum_i \left( \sum_jq_{1j,i}x_j
+ q_{0,i} \right) \partial_i f
+ \left( \sum_{i,j}r_{2ij}x_ix_j
+ \sum_{i}r_{1i}x_i
+ r_0 \right) f \\
&=& \nabla^TP_0\nabla f + \x^TQ_1\nabla f + \qv_0\cdot\nabla f + \left( \x^TR_2\x + \rv_1\cdot\x + r_0 \right) f.
\end{eqnarray*}
Equating terms of the same order in $\x$ (after dividing by $f$) gives:
\begin{eqnarray*}
\dot A &=& 4AP_0A + 2\sym(Q_1A) + R_2 \\
\dot\bv &=& 4AP_0\bv + Q_1\bv + 2A\qv_0 + \rv_1 \\
\dot c &=& 2\tr(AP_0) + \bv^TP_0\bv + \qv_0\cdot\bv + r_0,
\end{eqnarray*}
where $\tr(A) = \sum_i a_{ii}$ and $\sym(A) = \frac12(A+A^T)$.  The symmetrisation is used to keep the derivative of $A$ symmetrical even though $Q_1$ may not be.


\section{Operations on Wavepackets}

\subsection{Pointwise Operations}

\subsubsection{Addition} % includes adding a constant 
The sum of two wavepackets is in general not another wavepacket (with some exceptions below), so more complex functions must be represented by a sum
\[ f(x) = \sum_n e^{a_nx^2+b_nx+c_n}. \]
A constant $k$ can be represented as $e^{\ln k}$ where $\ln k$ is a suitable version of the complex natural logarithm, so setting $a=b=0$, $c=\ln k$ gives a wavepacket that can be added to the sum.

\subsubsection{Simplification of Sums}
If two wavepackets have equal $a$ and $b$, then they can be combined:
\[ e^{ax^2+bx+c_1}+e^{ax^2+bx+c_2} = e^{ax^2+bx}\left(e^{c_1}+e^{c_2}\right) = e^{ax^2+bx+\ln(e^{c_1}+e^{c_2})}. \]
If $e^{c_1}+e^{c_2}=0$ then the two wavepackets cancel.

\subsubsection{Scalar Multiplication} % includes negation
\[ ke^{ax^2+bx+c} = e^{ax^2+bx+(c+\ln k)} \]
\[ -e^{ax^2+bx+c} = e^{ax^2+bx+(c+\pi i)} \]

\subsubsection{Multiplication of Two Wavepackets}
\[ e^{a_1x^2+b_1x+c_1}e^{a_2x^2+b_2x+c_2} = e^{(a_1+a_2)x^2+(b_1+b_2)x+(c_1+c_2)} \]

\subsubsection{Powers}
For real integer $k$,
\[ \left(e^{ax^2+bx+c}\right)^k = e^{kax^2+kbx+kc}. \]

\subsubsection{Modulus and Phase}\label{-modpha}
For real $x$,
\[ \left|e^{ax^2+bx+c}\right| = e^{(\re a)x^2+(\re b)x+\re c}; \]
\[ \arg e^{ax^2+bx+c} = (\im a)x^2+(\im b)x+\im c. \]

\subsubsection{Complex Conjugate}
For real $x$,
\[ \overline{e^{ax^2+bx+c}} = e^{\bar{a}x^2+\bar{b}x+\bar{c}}. \]

\subsubsection{Real and Imaginary Parts}
For real $x$,
\[ \re e^{ax^2+bx+c} = e^{ax^2+bx+(c-\ln 2)} + e^{\bar{a}x^2+\bar{b}x+(\bar{c}-\ln 2)}; \]
\[ \im e^{ax^2+bx+c} = e^{ax^2+bx+(c-\ln 2-\frac\pi2 i)} + e^{\bar{a}x^2+\bar{b}x+(\bar{c}-\ln 2+\frac\pi2 i)}. \]

\subsection{Geometrical Operations}
\subsubsection{Centre}\label{-centre}
Every wavepacket with a finite integral has a single point where its modulus is a maximum.  This is
\[ \arg\max \left|e^{ax^2+bx+c}\right| 
 = \arg\max e^{(\re a)x^2+(\re b)x+\re c}
 = \arg\max ((\re a)x^2+(\re b)x+\re c).
\]
At the maximum, the derivative $2(\re a)x+\re b$ will be equal to zero, thus $x = \frac{-\re b}{2 \re a}$.

In multiple dimensions we need to find $\arg\max (\x^T(\re A)\x+(\re\bv)\cdot\x+\re c)$.  Assuming that $A=A^T$, the vector gradient of this is $2(\re A)\x+\re\bv$.  This is zero when $\x=-\frac12(\re A)^{-1}\re\bv$.

\subsubsection{Translation}
Displacing a wavepacket a distance $s$ gives
\[ e^{a(x-s)^2+b(x-s)+c} = e^{ax^2+(-2as+b)x+(as^2-bs+c)}. \]
For multiple dimensions and a vector displacement $\sv$, this is
\[ e^{(\x-\sv)^TA(\x-\sv) + \bv\cdot(\x-\sv) + c} = e^{\x^TA\x + (-2A\sv+\bv)\cdot\x + (\sv^TA\sv-\bv\cdot\sv+c)}, \]
assuming $A=A^T$.

\subsubsection{Rotation}

\subsection{Integral Operations}

\subsubsection{Integral Over All Space}\label{-int}
\newcommand{\intR}[1]{\int^\infty_{-\infty} #1 \dx}
This can be computed by using the identity $\intR{e^{-\pi x^2}}=1$.
\begin{eqnarray*}
\intR{e^{ax^2+bx+c}} &=& \intR{e^{a\left(x^2+\frac{b}{a}x+\frac{c}{a}\right)}} \\
&=& \intR{e^{a\left(\left(x+\frac{b}{2a}\right)^2-\frac{b^2}{4a^2}+\frac{c}{a}\right)}} \\ 
&=& e^{c-\frac{b^2}{4a}} \intR{e^{a\left(x+\frac{b}{2a}\right)^2}} \\ 
&=& e^{c-\frac{b^2}{4a}} \intR{e^{ax^2}} \\ 
&=& e^{c-\frac{b^2}{4a}} \sqrt{\frac{-\pi}{a}} \intR{e^{-\pi x^2}} \\ 
&=& e^{c-\frac{b^2}{4a}} \sqrt{\frac{-\pi}{a}}.
\end{eqnarray*}

In the case of complex coefficients, the identity $\intR{e^{ax^2}}=\sqrt{\frac{-\pi}{a}}$ holds as long as $\re a<0$.  The shift of origin $\intR{e^{a\left(x+\frac{b}{2a}\right)^2}} =  \intR{e^{ax^2}}$ is not intuitively true if $b/a$ is not real.  However, a contour integration argument shows the infinite line of integration can also be moved in the imaginary direction without affecting the integral so long as it does not cross a pole (there are no poles in $e^{Q(x)}$) and the integral remains convergent.

The case for multiple dimensions is analogous but slightly more complicated; three dimensions will be used here as an example.  The basic identity works in multiple dimensions:
\newcommand{\intRy}[1]{\int^\infty_{-\infty} #1 \dy}
\newcommand{\intRz}[1]{\int^\infty_{-\infty} #1 \dz}
\newcommand{\intRV}[1]{\int_{\mathbb{R}^3} #1 \dV} % needs \usepackage{amssymb} for \mathbb
\[ \intR{e^{-\pi x^2}}=\intRy{e^{-\pi y^2}}=\intRz{e^{-\pi z^2}}=1 \]
\[ \Rightarrow \qquad \intR{e^{-\pi x^2}}\intRy{e^{-\pi y^2}}\intRz{e^{-\pi z^2}}=1 \]
\[ \Rightarrow \qquad \intRV{e^{-\pi x^2}e^{-\pi y^2}e^{-\pi z^2}}=\intRV{e^{-\pi \x\cdot\x}}=1. \]
\[ \Rightarrow \qquad \intRV{e^{-\x\cdot\x}}=\pi^{3/2}, \]
where the final step is done by a scaling of $\sqrt\pi$ in each axis.

The matrix $A$ in the exponent $\x^TA\x$ can be symmetric without loss of generality, thus $-A$ is also symmetric.  Symmetric matrices always factorise as $-A=B^TB$.  This means that 
\[ \intRV{e^{\x^TA\x}}=\intRV{e^{-\x^TB^TB\x}}=\intRV{e^{-(B\x)^TB\x}}. \]
If $A$ was real, then $B$ is real and this is equal to 
\[ \frac1{\det B}\intRV{e^{-\x^T\x}} = \frac{\pi^{3/2}}{\sqrt{\det(-A)}}
= \sqrt{\frac{\pi^3}{\det(-A)}} =\sqrt{\frac{(-\pi)^3}{\det A}}, \]
where we have used $\det(-A)=\det{B^T}\det{B}=(\det{B})^2$ and the fact that $\det(-A)=(-1)^3\det A$ where $3$ is the number of dimensions of space.  The general version of this result incorporates the value $\sqrt{\frac{-\pi}{a}}$ for one dimension.  Subtle point: both $B$ and $-B$ work in the factorisation of $-A$, so we choose the one with positive determinant that will not invert the sign of the volume element and also ensures $\det B$ is the positive root $\sqrt{\det(-A)}$.

If $A$ has complex entries, similar matrix factorisations are possible but it is not intuitively clear that arguments about scaling the volume element with the determinant work.  However, it is probably correct to rely on the analytic continuation of the real result into the complex plane.

The final step to integrating the full exponent $\x^TA\x+\bv\cdot\x+c$ is `completing the square' in multiple dimensions.  Assuming $A=A^T$, consider a displacement $\dv$ to the quadratic term:
\[ (\x+\dv)^TA(\x+\dv) = \x^TA\x+2\dv^TA\x+\dv^TA\dv \]
and note that the linear term is reconstructed if $2\dv^TA\x=\bv^T\x$, which is true if $\dv=\frac12A^{-1}\bv$.  Now we have
\[ \x^TA\x+\bv\cdot\x+c =  (\x+\dv)^TA(\x+\dv) - \dv^TA\dv + c \]
and
\begin{eqnarray*}
\intRV{e^{\x^TA\x+\bv\cdot\x+c}} &=& \intRV{e^{(\x+\dv)^TA(\x+\dv) - \dv^TA\dv + c}} \\
&=& e^{c-\dv^TA\dv} \intRV{e^{(\x+\dv)^TA(\x+\dv)}} \\
&=& e^{c-\frac14\bv^TA^{-1}\bv} \intRV{e^{\x^TA\x}} \\
&=& e^{c-\frac14\bv^TA^{-1}\bv} \sqrt{\frac{(-\pi)^3}{\det A}}.
\end{eqnarray*}
The complex shifts of the range of integration are allowed as before and note that this formula incorporates the one-dimensional result (if 3 is replaced by 1).

\subsubsection{Inner Product}
\newcommand{\innerprod}[2]{\left(#1,#2\right)}
The standard inner product $\innerprod{f}{g}$ between complex functions is the integral $\int \bar f g$ over the whole space.  In one dimension,
\begin{eqnarray*}
 \innerprod{e^{a_1x^2+b_1x+c_1}}{e^{a_2x^2+b_2x+c_2}} &=& \intR{e^{(\bar a_1+a_2)x^2+(\bar b_1+b_2)x+(\bar c_1+c_2)}} \\
&=&  e^{\bar c_1+c_2-\frac{(\bar b_1+b_2)^2}{4(\bar a_1+a_2)}} \sqrt{\frac{-\pi}{\bar a_1+a_2}}. 
\end{eqnarray*}
In multiple dimensions, the analogous substitutions are made into the multidimensional integral-over-all-space formula.

\subsubsection{$L_2$ Norm}
\newcommand{\norm}[1]{||#1||}
The `Euclidean' or $L_2$ norm is defined via $\norm{f}=\sqrt{\innerprod{f}{f}}$ so that
\begin{eqnarray*}
 \norm{e^{ax^2+bx+c}} &=& \sqrt{\intR{e^{(\bar a+a)x^2+(\bar b+b)x+(\bar c+c)}}} \\
&=&  \sqrt{e^{2\re c-\frac{(2\re b)^2}{8\re a}} \sqrt{\frac{-\pi}{2 \re a}}} \\
&=& e^{\re c-\frac{(\re b)^2}{4\re a}} \sqrt[4]{\frac{-\pi}{2 \re a}}.
\end{eqnarray*}
The fact that $\innerprod{f}{f}=\int \bar f f=\int |f|^2$ means that $\norm{f}^2$ corresponds to the integral of the probability density for quantum wavefunctions.

\subsubsection{Fourier Transform}\label{-ft}
If the Fourier transform is defined by $\tilde f(\omega) = \frac1{\sqrt{2\pi}}\intR{f(x)e^{-i\omega x}}$,
%with inverse $f(x) = \frac1{\sqrt{2\pi}}\int^\infty_{-\infty}f(\omega)e^{i\omega x}\,\mathrm d\omega$,
then the transform of a 1D wavepacket, making use of the integral formulae in section \ref{-int}, is:
\begin{eqnarray*}
 \tilde f(\omega) &=& \frac1{\sqrt{2\pi}}\intR{e^{ax^2+bx+c}e^{-i\omega x}} \\
&=& \frac1{\sqrt{2\pi}}\intR{e^{ax^2+(b-i\omega)x+c}} \\
&=& \frac1{\sqrt{2\pi}}e^{c-\frac{(b-i\omega)^2}{4a}} \sqrt{\frac{-\pi}{a}} \\
&=& e^{c-\frac{b^2-2ib\omega-\omega^2}{4a}} \sqrt{\frac{-1}{2a}} \\
&=& e^{\frac{\omega^2}{4a}+\frac{ib}{2a}\omega+(c-\frac{b^2}{4a}-\frac12\ln(-2a))}.
\end{eqnarray*}
This has the form of a wavepacket in $\omega$ with parameters
\[ \tilde a = \frac1{4a}, \qquad 
\tilde b = \frac{ib}{2a}, \qquad
\tilde c = c-\frac{b^2}{4a}-\frac12\ln(-2a). \]
Applying these formulae twice gives $\tilde{\tilde a}=a$, $\tilde{\tilde b}=-b$ and $\tilde{\tilde c}=c$, consistent with the fact that the transform of the transform is $f(-x)$.  The inverse transform formula also gives a way to express a wavepacket as a superposition of waves:
\[ e^{ax^2+bx+c} = f(x) = \frac1{\sqrt{2\pi}}\int^\infty_{-\infty} \tilde f(\omega)e^{i\omega x} \,\mathrm d\omega= \frac1{\sqrt{2\pi}}\int^\infty_{-\infty}e^{\tilde a\omega^2+\tilde b\omega+\tilde c}e^{i\omega x} \,\mathrm d\omega. \]

\newcommand{\vomega}{\boldsymbol\omega}
\newcommand{\intRW}[1]{\int_{\mathbb{R}^3} #1 \,\mathrm d\vomega} % needs \usepackage{amssymb} for \mathbb
In three dimensions (although it generalises to any number), the Fourier transform of $f(\x)$ is done on each variable in turn, giving  $\tilde f(\vomega) = \frac1{(2\pi)^{3/2}}\intRV{f(\x)e^{-i\vomega\cdot\x}}$.  The transform of a 3D wavepacket is then:
\begin{eqnarray*}
 \tilde f(\vomega) &=& \frac1{(2\pi)^{3/2}}\intRV{e^{\x^TA\x + \bv\cdot\x + c}e^{-i\vomega\cdot\x}} \\
&=& \frac1{(2\pi)^{3/2}}\intRV{e^{\x^TA\x + (\bv-i\vomega)\cdot\x + c}} \\
&=& \frac1{(2\pi)^{3/2}} e^{c-\frac14(\bv-i\vomega)^TA^{-1}(\bv-i\vomega)} \sqrt{\frac{(-\pi)^3}{\det A}} \\
&=& e^{c+\frac14\vomega^TA^{-1}\vomega+\frac{i}2\bv^TA^{-1}\vomega-\frac14\bv^TA^{-1}\bv}  \sqrt{\frac{(-\frac12)^3}{\det A}} \\
&=& e^{\vomega^T\frac14A^{-1}\vomega + (\frac{i}2A^{-1}\bv)\cdot\vomega + (c-\frac14\bv^TA^{-1}\bv-\frac12\ln((-2)^3\det A))}.
\end{eqnarray*}
This has the form of a wavepacket in $\vomega$ with parameters
\[ \tilde A = \frac14A^{-1}, \qquad 
\tilde \bv = \frac{i}2A^{-1}\bv, \qquad
\tilde c = c-\frac14\bv^TA^{-1}\bv-\frac12\ln((-2)^3\det A). \]

\subsubsection{Wigner Function}
In quantum mechanics, the Wigner function of a wavefunction $\psi(x)$ is the following (always real) function of position $x$ and momentum $p$:
\[ W(x,p) = \frac1{\pi\hbar}\int_{-\infty}^{\infty}\bar\psi(x+y)\psi(x-y)e^{2ipy/\hbar} \dy. \]
If $\psi$ is a 1D wavepacket, this becomes
\begin{eqnarray*}
W(x,p) &=& \frac1{\pi\hbar}\int_{-\infty}^{\infty} e^{\bar a(x+y)^2+\bar b(x+y)+\bar c} e^{a(x-y)^2+b(x-y)+c} e^{2ipy/\hbar} \dy \\
&=& \frac1{\pi\hbar} e^{2(\re a)x^2+2(\re b)x+2\re c} \int_{-\infty}^{\infty} e^{2(\re a)y^2 + (-4i(\im a)x-2i\im b+2ip/\hbar)y} \dy \\
&=& \frac1{\pi\hbar} e^{2(\re a)x^2+2(\re b)x+2\re c} e^{-(-4i(\im a)x-2i\im b+2ip/\hbar)^2/(8\re a)} \sqrt{\frac{-\pi}{2\re a}} \\
&=& \frac1{\pi\hbar} e^{2(\re a)x^2+2(\re b)x+2\re c} e^{(-2(\im a)x-\im b+p/\hbar)^2/(2\re a)} \sqrt{\frac{-\pi}{2\re a}}.
\end{eqnarray*}
This is a positive 2D Gaussian in $(x,p)$ space, with the $x$ dependency (first exponential term) actually being equal to $|\psi(x)|^2$.  Recalling that $\arg\psi(x) = (\im a)x^2+(\im b)x+\im c$ and thus $\dd{\arg\psi(x)}{x} = 2(\im a)x+\im b$, this can also be written
\[ W(x,p) = \frac1{\pi\hbar} |\psi(x)|^2 e^{(p/\hbar-\dd{\arg\psi(x)}{x})^2/(2\re a)} \sqrt{\frac{-\pi}{2\re a}}. \]

\subsection{Splitting/Slicing}\label{-split}
The absolute value of a wavepacket (see section \ref{-modpha}) has a Gaussian profile, which we want to split into smaller Gaussians that add up to approximately the original.  The `phase' part $e^{i\im Q(x)}$ can be factored out during this whole process.  Additionally, a linear transformation in space and scaling in amplitude can make the absolute value equal to $e^{-x^2}$ (or $e^{-|\x|^2}$ in multiple dimensions).

Consider an approximation 
\[ e^{-x^2} \simeq \sum_{n=-\infty}^{\infty} c_n e^{-((x-na)/b)^2}. \]
This is a sum of Gaussians regularly spaced by $a$ and having width $b<1$ times the original.  The left-hand side integrates to $\sqrt{\pi}$ and each term on the right integrates to $c_nb\sqrt{\pi}$.  In the limit where $a$ is small, in the region near $x=0$, there are many Gaussians spaced by $a$ adding up to a nearly constant value of $1$, so we want $c_nb\sqrt{\pi}\simeq a$ near $x=0$.  That makes one choice for the coefficients
\[ c_n = \frac{a}{b\sqrt{\pi}} e^{-(na)^2}, \]
which simply modulates this by the overall envelope of the Gaussian we are trying to approximate.  For larger $b$, there will be some spreading in the overall variance, which can be counteracted with
\[ c_n = \frac{a}{b\sqrt{\pi(1-b^2)}} e^{-(na)^2/(1-b^2)}. \]

\subsection{Merging}
Although the sum of two (or more) wavepackets is usually not another wavepacket, it may sometimes be desired to merge them if a single wavepacket is a good enough approximation to the sum.


\section{Wavepackets in Relativistic Wave Equations}
When special relativity is taken into account, the set of wavepackets studied so far is no longer closed under time evolution.  However, they may still be used as initial conditions and arbitrary time evolution in the free wave equation can be calculated by Fourier methods.

\subsection{Example: Klein--Gordon Equation (non-quantum)}
A relativistic scalar particle field $\phi(x,t)$ propagates according the the wave equation:
\[ \frac1{c^2}\ddot\phi = \left(\nabla^2-\frac{m^2c^2}{\hbar^2}\right)\phi, \]
noting that $\phi$ is not a wavefunction this time, it is the shape of a single excitation of the scalar field operator, which is itself quantised.  Substituting a plane wave solution $\phi=e^{i(\Omega t+\vomega\cdot\x)}$ gives
\[ \frac{-1}{c^2}\Omega^2 = -|\vomega|^2-\frac{m^2c^2}{\hbar^2}
\qquad \Rightarrow \qquad \Omega = \pm c\sqrt{|\vomega|^2+\frac{m^2c^2}{\hbar^2}}. \]
There are two solutions because the differential equation is second order in time, so the initial conditions are two pieces of information: $\phi$ and $\dot\phi$.  The function $\Omega(\vomega)$ is defined to mean the positive square root in the expression above.  One of these solutions corresponds to the opposite charge anti-particle, which appears to have negative `energy' if the energy operator $i\hbar\partial_t$ from the Schr\"odinger equation is used on it.  As $i\hbar\partial_t\phi=i\hbar i\Omega\phi=-\hbar\Omega\phi$, the positive $\Omega$ solution is the anti-particle, labelled by $a$ below and the negative $\Omega$ solution is the particle, labelled $b$.  The combined solution can be expressed in Fourier space:
\begin{eqnarray*}
 \phi &=& \frac1{(2\pi)^{3/2}}\intRW{e^{i\vomega\cdot\x}(a(\vomega)e^{i\Omega(\vomega) t}+b(\vomega)e^{-i\Omega(\vomega) t})}, \\
 \dot\phi &=& \frac1{(2\pi)^{3/2}}\intRW{e^{i\vomega\cdot\x}i\Omega(\vomega)(a(\vomega)e^{i\Omega(\vomega) t}-b(\vomega)e^{-i\Omega(\vomega) t})} \\
 \phi_{t=0} &=& \frac1{(2\pi)^{3/2}}\intRW{e^{i\vomega\cdot\x}(a(\vomega)+b(\vomega))}, \\
 \dot\phi_{t=0} &=& \frac1{(2\pi)^{3/2}}\intRW{e^{i\vomega\cdot\x}i\Omega(\vomega)(a(\vomega)-b(\vomega))} \\
\Rightarrow \qquad
\tilde\phi_{t=0} &=& a(\vomega)+b(\vomega), \\
\tilde{\dot\phi}_{t=0} &=& i\Omega(\vomega)(a(\vomega)-b(\vomega)) \\
\Rightarrow \qquad
a(\vomega) &=& \frac12\left(\tilde\phi_{t=0}+\frac{\tilde{\dot\phi}_{t=0}}{i\Omega(\vomega)}\right), \\
b(\vomega) &=& \frac12\left(\tilde\phi_{t=0}-\frac{\tilde{\dot\phi}_{t=0}}{i\Omega(\vomega)}\right).
\end{eqnarray*}
With these coefficients known in terms of the initial conditions, the Fourier representation at any later time can be written:
\begin{eqnarray*}
\tilde\phi &=& a(\vomega)e^{i\Omega(\vomega) t}+b(\vomega)e^{-i\Omega(\vomega) t} \\
\Rightarrow \qquad 
\tilde\phi &=& \frac12\left(\tilde\phi_{t=0}+\frac{\tilde{\dot\phi}_{t=0}}{i\Omega(\vomega)}\right)e^{i\Omega(\vomega) t}+\frac12\left(\tilde\phi_{t=0}-\frac{\tilde{\dot\phi}_{t=0}}{i\Omega(\vomega)}\right)e^{-i\Omega(\vomega) t} \\
&=& \frac12(e^{i\Omega(\vomega) t}+e^{-i\Omega(\vomega) t})\tilde\phi_{t=0} + \frac1{2i\Omega(\vomega)}(e^{i\Omega(\vomega) t}-e^{-i\Omega(\vomega) t})\tilde{\dot\phi}_{t=0} \\
&=& \cos(\Omega(\vomega) t)\tilde\phi_{t=0} + \frac{\sin(\Omega(\vomega) t)}{\Omega(\vomega)}\tilde{\dot\phi}_{t=0}.
\end{eqnarray*}
If $a(\vomega)$ and $b(\vomega)$ can be expressed as 3D wavepackets as defined in section \ref{-multid} (or small sums of these), the time evolution problem reduces to calculating the Fourier transform of:
\[ \tilde\phi = e^{\vomega^T\tilde A\vomega + \tilde\bv\cdot\vomega + \tilde c}e^{i\Omega(\vomega) t}, \]
showing only the $a(\vomega)$ part since $-t$ can be substituted for $t$ in the $b(\vomega)$ terms.

\subsubsection{Splitting Scheme}
The time evolution is not a simple wavepacket because $\Omega(\vomega)$ is not a quadratic function of $\vomega$, but one approach is to approximate it piecewise by quadratic functions in local areas of frequency space.  If enough sufficiently-overlapping wavepackets are used with the correct normalisation, this can be an arbitrarily good approximation (see section \ref{-split}).  The value of $\Omega$ at an offset from a central frequency $\vomega_0$ can be written:
\begin{eqnarray*}
\Omega(\vomega_0+\vomega) &=& c\left(|\vomega_0+\vomega|^2+\frac{m^2c^2}{\hbar^2}\right)^{1/2} \\
&=& c\left(\frac{m^2c^2}{\hbar^2}+|\vomega_0|^2+2\vomega_0\cdot\vomega+|\vomega|^2\right)^{1/2} \\
&=& c\left(\frac1{c^2}\Omega(\vomega_0)^2+2\vomega_0\cdot\vomega+|\vomega|^2\right)^{1/2} \\
&=& \left(\Omega(\vomega_0)^2+2c^2\vomega_0\cdot\vomega+c^2|\vomega|^2\right)^{1/2} \\
&=& \Omega(\vomega_0)\left(1+\frac{2c^2}{\Omega(\vomega_0)^2}\vomega_0\cdot\vomega+\frac{c^2}{\Omega(\vomega_0)^2}|\vomega|^2\right)^{1/2} \\
&=& \Omega(\vomega_0)\left(1+2k^2\vomega_0\cdot\vomega+k^2|\vomega|^2\right)^{1/2},
\end{eqnarray*}
where $k=\frac{c}{\Omega(\vomega_0)}$ has been defined for brevity.  Using the Taylor expansion $(1+x)^{1/2}=1+\frac12x-\frac18x^2+\frac1{16}x^3+O(x^4)$, the square root becomes
\[
1
+\frac12\left(2k^2\vomega_0\cdot\vomega+k^2|\vomega|^2\right)
-\frac18\left(4k^4(\vomega_0\cdot\vomega)^2+4k^4(\vomega_0\cdot\vomega)|\vomega|^2\right)
+\frac1{16}\left(8k^6(\vomega_0\cdot\vomega)^3\right)
+O(|\vomega|^4)
\]
\[ = 
1
+k^2\vomega_0\cdot\vomega+\frac{k^2}2|\vomega|^2
-\frac{k^4}2(\vomega_0\cdot\vomega)^2-\frac{k^4}2(\vomega_0\cdot\vomega)|\vomega|^2
+\frac{k^6}2(\vomega_0\cdot\vomega)^3
+O(|\vomega|^4).
\]
Renaming $\vomega$ above to $\vomega_1$, the quadratic part can be incorporated into the existing wavepacket as a function of $\vomega=\vomega_0+\vomega_1$:
\begin{eqnarray*}
 \tilde\phi &=& e^{\vomega^T\tilde A\vomega + \tilde\bv\cdot\vomega + \tilde c + i\Omega(\vomega) t} \\
&\simeq& e^{\vomega^T\tilde A\vomega + \tilde\bv\cdot\vomega + \tilde c + i\Omega(\vomega_0)
\left(1+k^2\vomega_0\cdot\vomega_1+\frac{k^2}2|\vomega_1|^2-\frac{k^4}2(\vomega_0\cdot\vomega_1)^2\right) t}.
\end{eqnarray*}
Substituting $\vomega_1=\vomega-\vomega_0$ makes the expression in brackets equal to:
\[ 1+k^2\vomega_0\cdot(\vomega-\vomega_0)+\frac{k^2}2|\vomega-\vomega_0|^2-\frac{k^4}2(\vomega_0\cdot(\vomega-\vomega_0))^2 \]
\[ = 1
+k^2\vomega_0\cdot\vomega-k^2|\vomega_0|^2
+\frac{k^2}2\left(|\vomega|^2-2\vomega_0\cdot\vomega+|\vomega_0|^2\right)
-\frac{k^4}2(\vomega_0\cdot\vomega-|\vomega_0|^2)^2 \]
\[ = 1
-\frac{k^2}2|\vomega_0|^2
+\frac{k^2}2|\vomega|^2
-\frac{k^4}2|\vomega_0|^4+k^4|\vomega_0|^2(\vomega_0\cdot\vomega)-\frac{k^4}2(\vomega_0\cdot\vomega)^2, \]
so that the (approximate) time-evolved parameters of the wavepacket are:
\begin{eqnarray*}
\tilde A(t) &=& \tilde A+i\Omega(\vomega_0)t\left(\frac{k^2}2I-\frac{k^4}2\vomega_0\vomega_0^T\right), \\
\tilde\bv(t) &=& \tilde\bv+i\Omega(\vomega_0)t\left(k^4|\vomega_0|^2\vomega_0\right), \\
\tilde c(t) &=& \tilde c+i\Omega(\vomega_0)t\left(1-\frac{k^2}2|\vomega_0|^2-\frac{k^4}2|\vomega_0|^4\right).
\end{eqnarray*}
Note again that $k=\frac{c}{\Omega(\vomega_0)}$ and these may be turned back into real space with the inverse of the formula in section \ref{-ft}:
\[ A = \frac14\tilde A^{-1}, \qquad 
 \bv = \frac{-i}2\tilde A^{-1}\tilde\bv, \qquad
 c = \tilde c-\frac14\tilde\bv^T\tilde A^{-1}\tilde \bv-\frac12\ln((-2)^3\det \tilde A). \]

\subsubsection{Error Bound}
The error of the quadratic approximation to $\Omega(\vomega_0+\vomega)$ used in the previous section is:
\[ E=\Omega(\vomega_0)\left(-\frac{k^4}2(\vomega_0\cdot\vomega)|\vomega|^2+\frac{k^6}2(\vomega_0\cdot\vomega)^3\right)+O(|\vomega|^4). \]
Assuming that the leading order cubic contribution (denoted $E_3$) dominates, an upper bound on the error can be obtained using the triangle inequality:
\begin{eqnarray*}
|E_3| &\le& \Omega(\vomega_0)\left(\left|\frac{k^4}2(\vomega_0\cdot\vomega)|\vomega|^2\right|+\left|\frac{k^6}2(\vomega_0\cdot\vomega)^3\right|\right) \\
&\le& \Omega(\vomega_0)\left(\frac{k^4}2|\vomega_0||\vomega|^3+\frac{k^6}2|\vomega_0|^3|\vomega|^3\right)
\end{eqnarray*}
% Could try this again with direction dependence (don't use tri-ineq so early)

\subsubsection{Evolution in Real Space}
The previous sections have assumed the wavepacket parameters will be evolved in Fourier space and only converted back to real space at the end.  However, for some applications the wavepackets may be specified in real space to start with.  The time derivatives of the parameters in real and Fourier space are related by:
\[ \dot A = \frac14\partial_t(\tilde A^{-1}), \qquad 
 \dot\bv = \frac{-i}2\partial_t(\tilde A^{-1}\tilde\bv), \qquad
 \dot c = \dot{\tilde c}-\frac14\partial_t(\tilde\bv^T\tilde A^{-1}\tilde \bv)-\frac12\partial_t\ln((-2)^3\det \tilde A). \]
To evaluate this requires the two formulae $\partial_t(A^{-1})=-A^{-1}\dot AA^{-1}$ and $\partial_t\det A=\det A\tr(A^{-1}\dot A)$, as well as the transform formulae such as $\tilde A^{-1} = 4A$ and $\tilde A^{-1}\tilde\bv = 2i\bv$ to get the right-hand side mostly in terms of real space again.
\[ \dot A = \frac{-1}4\tilde A^{-1}\dot{\tilde A}\tilde A^{-1}
 = -4A\dot{\tilde A}A, \]
\[ \dot\bv = \frac{-i}2\left(\partial_t(\tilde A^{-1})\tilde\bv+\tilde A^{-1}\dot{\tilde\bv}\right) 
 = \frac{i}2\tilde A^{-1}\dot{\tilde A}\tilde A^{-1}\tilde\bv-\frac{i}2\tilde A^{-1}\dot{\tilde\bv}
 = -4A\dot{\tilde A}\bv-2iA\dot{\tilde\bv}, \]
\[ \dot c = \dot{\tilde c}-\frac14\partial_t(\tilde\bv^T\tilde A^{-1}\tilde \bv)-\frac12\partial_t\ln((-2)^3\det \tilde A) \]
\[ = \dot{\tilde c}-\frac14\tilde\bv^T\partial_t(\tilde A^{-1})\tilde \bv-\frac12\dot{\tilde\bv}^T\tilde A^{-1}\tilde \bv-\frac1{2(-2)^3\det \tilde A}(-2)^3\partial_t\det \tilde A \]
\[ = \dot{\tilde c}+\frac14\tilde\bv^T\tilde A^{-1}\dot{\tilde A}\tilde A^{-1}\tilde \bv-i\dot{\tilde\bv}^T\bv-\frac{\det \tilde A\tr(\tilde A^{-1}\dot{\tilde A})}{2\det \tilde A} \]
\[ = \dot{\tilde c}-\bv^T\dot{\tilde A}\bv-i\dot{\tilde\bv}^T\bv-2\tr(A\dot{\tilde A}). \]
The specific time derivatives of the Fourier-space parameters for the Klein--Gordon equation come from the previous section:
\begin{eqnarray*}
\dot{\tilde A} &=& i\Omega(\vomega_0)\left(\frac{k^2}2I-\frac{k^4}2\vomega_0\vomega_0^T\right), \\
\dot{\tilde\bv} &=& i\Omega(\vomega_0)\left(k^4|\vomega_0|^2\vomega_0\right), \\
\dot{\tilde c} &=& i\Omega(\vomega_0)\left(1-\frac{k^2}2|\vomega_0|^2-\frac{k^4}2|\vomega_0|^4\right).
\end{eqnarray*}
The only ambiguity is how to choose $\vomega_0$, which is the frequency about which the quadratic approximation is most accurate.  A natural choice is the centre (as defined in section \ref{-centre}) of the Fourier-space wavepacket
\[ \vomega_0 = -\frac12(\re \tilde A)^{-1}\re\tilde\bv. \]

%%%
\comment{Using derivative formulae from a previous section, a 3D wavepacket  $f=e^{\x^TA\x + \bv\cdot\x + c}$ satisfies
\begin{eqnarray*}
\nabla^2f = \partial_i\partial_i f &=& (2a_{ii}+(2a_{ij}x_j+b_i)(2a_{ik}x_k+b_i))f \\
&=& (2\tr A + |2A\x+\bv|^2)f,
\end{eqnarray*}
where repeated indices are implicitly summed over.}

\comment{Since this equation is second-order in time, it can be rewritten as a first-order equation acting on a 2-component field:
\[ \pdd{}{t}\vectwo{\phi}{\dot\phi} = \vectwo{\dot\phi}{\left(c^2\nabla^2-\frac{m^2c^4}{\hbar^2}\right)\phi}. \]}

\subsection{Example: Dirac Equation (non-quantum)}
The spinor field $\psi_\alpha(\x,t)$ of electrons and positrons satisfies the equation $i\hbar\gamma^\mu\partial_\mu\psi-mc\psi=0$, using the Einstein summation convention and the standard gamma matrices of particle physics.  This notation suppresses the spinor indices $\alpha=1...4$, each spinor component being a complex number, so the full equation can be written $i\hbar\gamma^\mu_{\alpha\beta}\partial_\mu\psi_\beta-mc\psi_\alpha = 0$. 

Substituting a plane wave solution $\psi=\alpha e^{i(\Omega t+\vomega\cdot\x)}$ with $\alpha$ a constant spinor gives
\begin{eqnarray*}
\partial_0 \psi &=& i\Omega \psi \qquad /**/\\ % Probable error here: one of these derivatives probably has a factor of c in it
\qquad \partial_i \psi &=& i\omega_i \psi, 
\end{eqnarray*}
so that 
\[ -\hbar\gamma^0\Omega\psi-\hbar\gamma^i\omega_i\psi-mc\psi=0. \]
The exponential term in $\psi$ is never zero, so dividing it out both sides (and changing sign) leaves:
\[ (\hbar\gamma^0\Omega+\hbar\gamma^i\omega_i+mc)\alpha=0. \]
The term in brackets is a 4$\times$4 complex matrix that acts on spinors (like the gamma matrices).  A solution with $\alpha=0$ would be zero everywhere, so nontrivial solutions require that
\[ \det(\hbar\Omega\gamma^0+\hbar\omega_i\gamma^i+mcI_4) = 0 \]
because the matrix must send some nonzero spinor $\alpha$ to zero.

Written out in full (using the Dirac representation for $\gamma^\mu$), the determinant is
\[
\left|\begin{array}{cccc}
mc+\hbar\Omega&0&\hbar\omega_3&\hbar\omega_1-i\hbar\omega_2\\
0&mc+\hbar\Omega&\hbar\omega_1+i\hbar\omega_2&-\hbar\omega_3\\
-\hbar\omega_3&-\hbar\omega_1+i\hbar\omega_2&mc-\hbar\Omega&0\\
-\hbar\omega_1-i\hbar\omega_2&\hbar\omega_3&0&mc-\hbar\Omega
\end{array}\right|
\]
\[ =((mc)^2-(\hbar\Omega)^2+(\hbar\omega_1)^2+(\hbar\omega_2)^2+(\hbar\omega_3)^2)^2 = (m^2c^2+\hbar^2(|\vomega|^2-\Omega^2))^2. \]
So the angular frequency $\Omega$ must satisfy
\[ m^2c^2+\hbar^2(|\vomega|^2-\Omega^2)=0
\qquad \Rightarrow \qquad
\Omega^2=|\vomega|^2+\frac{m^2c^2}{\hbar^2}. \]


\section{Wavepackets and Quantum Fields}

\subsection{Fermionic Fields}
\newcommand{\y}{\mathbf y}
\newcommand{\vac}{|0\rangle}
\newcommand{\dthreex}{\,\mathrm d^3 \x}
\newcommand{\dthreey}{\,\mathrm d^3 \y}
Fermionic field operators satsify the canonical anticommutation relations (CAR)
\[ \{\psi_\alpha(\x),\psi_\beta^\dag(\y)\}=\delta_{\alpha\beta}\delta^3(\x-\y) 
\qquad \mbox{and} \qquad
\{\psi_\alpha(\x),\psi_\beta(\y)\}=\{\psi_\alpha^\dag(\x),\psi_\beta^\dag(\y)\}=0, \]
where $\{A,B\}=AB+BA$ and spinor indices are shown, since most fermions are spinors.  The zero state also satisfies $\psi_\alpha(\x)\vac=0$.  States containing particles are created by applying the $\psi^\dag$ (creation-type) operators to $\vac$.  However, a particle existing at only one point is not differentiable, so define the function-shaped creation operator (and its Hermitian conjugate)
\[ \psi_\alpha^\dag[f] = \int f(\x) \psi_\alpha^\dag(\x) \dthreex,
\qquad
\psi_\alpha[f] = \int \bar f(\x) \psi_\alpha(\x) \dthreex. \]
A single application of this gives the state that shall be labelled $|f_\alpha\rangle=\psi_\alpha^\dag[f]\vac$.  The normalisation of this state works like this:
\begin{eqnarray*}
 \langle f_\alpha|f_\alpha\rangle &=& \langle 0| \int \bar{f}(\y) \psi_\alpha(\y) \dthreey \int f(\x) \psi_\alpha^\dag(\x) \dthreex\,\vac \\
&=& \langle 0| \iint f(\x) \bar{f}(\y) \psi_\alpha(\y) \psi_\alpha^\dag(\x) \dthreex \dthreey\,\vac \\ 
&=& \langle 0| \iint f(\x) \bar{f}(\y) \left(\delta_{\alpha\alpha}\delta^3(\y-\x) - \psi_\alpha^\dag(\x) \psi_\alpha(\y)\right) \dthreex \dthreey\,\vac \\ 
&=& \langle 0| \iint f(\x) \bar{f}(\y) \delta^3(\y-\x) \dthreex \dthreey\,\vac \\ 
&=& \int |f(\x)|^2 \dthreex, \\ 
\end{eqnarray*}
using the rearrangement $AB=\{A,B\}-BA$ with the CAR value, followed by the fact that $\psi$ states annihilate on $\vac$, so the term ending in $\psi$ can be dropped, leaving only a scalar in the integral and finally $\langle 0\vac=1$.  Thus $|f_\alpha\rangle$ is normalised if the norm squared of $f$ integrates to one, which will be assumed from now on.

It will be useful to know the commutation properties of the wavepacket creation operators:
\begin{eqnarray*}
\{\psi_\alpha[f],\psi_\beta^\dag[g]\}
&=&  \iint \bar f(\x) g(\y) \{\psi_\alpha(\x),\psi_\beta^\dag(\y)\} \dthreex \dthreey \\ 
&=&  \iint \bar f(\x) g(\y) \delta_{\alpha\beta}\delta^3(\x-\y) \dthreex \dthreey \\ 
&=&  \delta_{\alpha\beta} \int \bar f(\x) g(\x) \dthreex = \delta_{\alpha\beta} (f,g), \\ 
\{\psi_\alpha[f],\psi_\beta[g]\} = \{\psi_\alpha^\dag[f],\psi_\beta^\dag[g]\} &=& 0,
\end{eqnarray*}
where the inner product on functions $(f,g)=\int \bar f(\x) g(\x) \dthreex$ has been used.  This gives a faster way of calculating
\[ \langle f_\alpha|f_\alpha\rangle
= \langle 0| \psi_\alpha[f] \psi_\alpha^\dag[f] \vac
= \langle 0| \left( \delta_{\alpha\alpha} (f,f) -\psi_\alpha^\dag[f] \psi_\alpha[f] \right) \vac
= (f,f). \]
If a two-particle state is written $|f_\alpha,g_\beta\rangle= \psi_\beta^\dag[g] \psi_\alpha^\dag[f] \vac$, then
\begin{eqnarray*}
\langle f_\alpha,g_\beta|f_\alpha,g_\beta\rangle 
&=& \langle 0| \psi_\alpha[f]\psi_\beta[g] \psi_\beta^\dag[g]\psi_\alpha^\dag[f] \vac \\
&=& \langle 0| \psi_\alpha[f]\left(\delta_{\beta\beta}(g,g) -\psi_\beta^\dag[g] \psi_\beta[g] \right)\psi_\alpha^\dag[f] \vac \\
&=& (f,f)(g,g) -  \langle 0| \psi_\alpha[f]\psi_\beta^\dag[g] \psi_\beta[g]\psi_\alpha^\dag[f] \vac \\
&=& (f,f)(g,g) -  \langle 0| \psi_\alpha[f]\psi_\beta^\dag[g] \left( \delta_{\beta\alpha}(g,f) - \psi_\alpha^\dag[f]\psi_\beta[g] \right) \vac \\
&=& (f,f)(g,g) - \delta_{\beta\alpha}(g,f) \langle 0| \psi_\alpha[f]\psi_\beta^\dag[g] \vac \\
&=& (f,f)(g,g) - \delta_{\beta\alpha}(g,f) \delta_{\alpha\beta}(f,g) \\
&=& (f,f)(g,g) - \delta_{\alpha\beta} |(f,g)|^2. \\
\end{eqnarray*}
If $f$ and $g$ are normalised, then this is $1 - \delta_{\alpha\beta} |(f,g)|^2$, so the two particle state is normalised unless $\alpha=\beta$ and $(f,g)\ne 0$.  In such `overlapping' cases, the amplitude is reduced: in the case where $f=g$, it is reduced to zero.

\end{document}