### Background  
1. C-field theory
 * Projected fields
 * Projected Gross-Pitaevskii Equation
2. Spectral-Galerkin method
 * $2M$ Rule for quadrature integration
 * Integrating general nonlinearities

### C-field definition
\begin{equation}
\psi(x,t)=\sum_n c_n(t)\phi_n(x),
\end{equation}
where the modes $\phi_n(x)$ are a complete and orthonormal set:
\begin{equation}
\int dx\;\phi^*_n(x)\phi_m(x)=\delta_{nm},\\
\mathbb{1}=\sum_n|\phi_n\rangle\langle\phi_n|.
\end{equation}
Completeness then gives a construction for the Dirac delta:
\begin{equation}
\delta(x-x')=\sum_n\phi_n(x)\phi_n^*(x').
\end{equation}

### Defining the C-region
To define the C-region, we introduce the projection operator
\begin{equation}
\hat{\cal P}=\bar{\sum_n}|\phi_n\rangle\langle \phi_n|
\end{equation}

where the bar notation refers to the cutoff defining the $C$-region. The quantum number is restricted to $n\leq M-1$, giving a total of $M$ modes in the 1D $C$-region, indexed by quantum numbers $n\in C=\{ 0,1,\dots, M-1\}$.

The Dirac delta function for the $C$-region is 
\begin{equation}
\langle x|\hat{\cal P}|x'\rangle=\bar{\sum_n}\phi_n(x)\phi_n^*(x')\equiv\delta(x,x').
\end{equation}

### Projected $\delta$-function
The projected delta function $\delta(x,x')$ plays a fundamental role in  c-field theory. It forms the kernel of the projector acting on functions in  position space:
\begin{equation}
{\cal P}f(x)\equiv \int dx'\;\delta(x,x')f(x')
\end{equation}
For any function residing in $C$, projection is the identity. Equivalently, the action of $\delta(x,x')$ is equivalent to the defining property of $\delta(x-x')$:
\begin{equation}
\int dx'\;\delta(x,x')f(x')=f(x).
\end{equation}
However, when considering $\delta(x,x')$, the order of arguments is significant.

## Projected  Gross-Pitaevskii equation
### Introduction to PGPE
Given a basis as above, one start from the GPE
\begin{equation}
i\hbar\frac{\partial \psi}{\partial t}=\left(-\frac{\hbar^2\nabla^2}{2m}+V(x)+g|\psi|^2\right)\psi
\end{equation}
and substitute 
\begin{equation}
\psi=\bar{\sum_n}c_n(t)\phi_n(x),
\end{equation}
and project onto $\phi_n(x)$.

## Spectral basis
We express the equation in the spectral basis 
\begin{equation}
(\hat T+\hat V)\phi_n(x)=\epsilon_n\phi_n(x).
\end{equation}
We arrive at the equation
\begin{equation}
i\hbar\frac{\partial c_n}{\partial t}=\epsilon_nc_n+g\int dx\;\phi_n^*(x)|\psi(x)|^2\psi(x)
\end{equation}
In the harmonic oscillator basis, our modes take the form 
\begin{equation}
\phi_n(x)=e^{-x^2/2}P_n(x)
\end{equation}
where $P_n(x)$ contains the polynomial part of degree $n$ and normalization factors. 

There are two obvious ways to express the integral.
#### 1. Direct substitution of basis decomposition
\begin{equation}
I_n\equiv \int dx\;\phi_n^*(x)|\psi(x)|^2\psi(x)
\end{equation}
Firstly, we can substitute the wavefunction:
\begin{equation}
I_n=\bar{\sum}_{mpq}c_m^*c_pc_q\int dx\;\phi_n^*(x)\phi_m^*(x)\phi_p(x)\phi_q(x),
\end{equation}
a triple summation over a set of known matrix elements for each mode $n$. For each $n$ this triple summation  requires $O(M^3)$ arithmetic operations. The total cost of evaluating the PGPE right hand side via summation is $O(M^4)$. This approach is numerically expensive, and scales poorly with $M$.

####  2. Quadrature method
This computational effort can be reduced via an alternate approach using Gaussian quadrature. 

>#### Philosophy:   
We will transform onto a *particular* (in general non-uniform) spatial grid that will guarantee that the non-linear term is projected faithfully onto $\phi_n(x)$, *for all modes in the C-field*.

For the integral at hand, the highest order term can be found by considering the cutoff mode for each wavefunction $\psi(x)$, and can be written as
\begin{equation}
I_{M-1}=\int dx\;e^{-2x^2}R_{4(M-1)}(x),
\end{equation}
where $R$ expresses the polynomial part.

### The $2M$ rule
As we will see below, for this $C$ region consisting of $M$ modes, the nonlinear integral may be evaluated exactly by mapping the coefficients $c_n$ to a spatial grid of $2M$ points, corresponding to a Gaussian quadrature rule of order  $2M$. Thus we can state our general rule, for a cubic nonlinearity of Gross-Pitaevskii form:

>#### The 2M rule:
A field represented by $M$ modes may be evolved at working precision using a particular grid of $2M$ spatial points.

The Fourier spectral method is a special case for which the wieght function becomes unity, and for which the $2M$ rule achieves spatial sampling at the *Nyquist frequency*.

### Gauss-Hermite quadrature
For any integral of the form
\begin{equation}
I=\int_{-\infty}^{\infty} dx \;e^{-x^2}Q_{2N-1}(x),
\end{equation}
with $Q_{2N-1}(x)$ restricted to maximum polynomial degree $2N-1$
\begin{equation}
Q_{2N-1}(x)=a_0+a_1x+\dots+a_{2N-1}x^{2N-1},
\end{equation}
there exists a quadrature rule  of order $N$, involving $N$ roots $x_k$, and $N$ weights $w_k$, $k\in 1,\dots, N$, that will evaluated all such integrals exactly. 

From a linear algebra point of view this result is not so surprising: The Gaussian-weighted integral of an arbitrary polynomial of degree $2N-1$ involves $2N$ uknown coefficients, and $2N$ known integrals. 

Given $2N$ free parameters, here given by the roots and weights of the quadrature rule, the $2N$ unknowns may be solved for exactly, by solving a system of simultaneous linear equations. 

The integral may be evaluated "exactly" as
\begin{equation}
I=\sum_k w_kQ_{2N-1}(x_k),
\end{equation}
where the meaning of  "exact" is that it is accurate to machine precision; typically the accuracy will be of order 1e-16. 

### Projected time evolution: work step
#### Algorithm
>1. cast the integral in the above form $I$, 
2. Transform to the quadrature grid $c_n\to \psi(x_k)$
2. Evaluate $Q_{2N-1}(x_k)$ for the resulting polynomial part, and 
3. Evaluate the sum weighted by $w_k$, giving the projection of $|\psi|^2\psi\to c_n'$

Step 1. is done by the change of variables $x=y/\sqrt{2}$

we find
\begin{align}
I_{M-1}&=\int_{-\infty}^{\infty} dy\;e^{-y^2}\left[\frac{R_{4(M-1)}(y/\sqrt{2})}{\sqrt{2}}\right],
\end{align}
where the terms in parentheses correspond to $Q_{2N-1}(x)$ defined above. 

According to the rules of Gaussian quadrature, a rule of order $N$ will integrate all polynomials up to and including degree $2N-1$. A rule of order $2M$ will thus evaluate all terms up to order $4M-1$. This is the lowest order rule that will guarantee exactness.
Hence
>**All projected integrals in the PGPE for $M$ modes are exact within a Gaussian quadrature of order $2M$**.

A rule of order $n$ will integrate polynomials up to and including degree $2n-1$. Thus a rule of order $JM$ will integrate terms up to order $2JM-1$, giving numerically exact evalution of the nonlinear term for all modes in $C$. The variable transformation is 
$$ x=\frac{y}{\sqrt{J}}$$
giving the integral 
$$I_{M-1}=\int_{-\infty}^{\infty} dy\;e^{-y^2} \frac{1}{\sqrt{J}}R_{2J(M-1)}(y/\sqrt{J}),$$

All powers of $\psi(x)$, i.e. both odd and even, are accommodated by putting $2J\to K$, to give
$$
I_{M-1}\equiv \int dx\;\phi_n^*(x)\psi(x)^{K-2}\psi(x),
$$
to give the form
$$
I_{M-1}=\int_{-\infty}^{\infty} dx \;e^{-Kx^2/2}R_{K(M-1)}(x).
$$
A rule of order $KM/2$ will integrate all terms up to order $KM-1>K(M-1)$. The variable transformation is 
$$x=\sqrt{\frac{2}{K}}y,$$ 
giving the integral
$$I_{M-1}=\int_{-\infty}^{\infty} dy\;e^{-y^2} \sqrt{\frac{2}{K}}R_{K(M-1)}\left(\sqrt{\frac{2}{K}}y\right) ,$$

### General Anisotropic form
The most general formulation we require must accommodate different oscillator fequencies in different spatial directions:
$$
I_{M-1}\equiv \int dx\;\phi_n^*(x)\psi(x)^{K-2}\psi(x),
$$
to give the form
$$
I_{M-1}=\int_{-\infty}^{\infty} dx \;e^{-K\omega x^2/2}R_{K(M-1)}(x).
$$
A rule of order $KM/2$ will integrate all terms with polynomial degree up to $KM-1>K(M-1)$. The variable transformation is 
$$x=\sqrt{\frac{2}{K\omega}}y,$$ 
giving the integral
$$I_{M-1}=\int_{-\infty}^{\infty} dy\;e^{-y^2} \sqrt{\frac{2}{K\omega}}R_{K(M-1)}\left(\sqrt{\frac{2}{K\omega}}y\right).$$

In the package `ProjectedGPE.jl`, the function that sets up the roots `x`, weights `w`, and transformation matrix `T` to evaluate integrals of this form is 
`nfieldtrans.jl`

You should be able to run this example to compute wavefunction norm:


```julia
using ProjectedGPE;
```


```julia
M = 20
x,w,T = nfieldtrans("Hermite",M,2) #2 field product in the integral to be evaluated
c = randn(M) + im*randn(M);
ψ = T*c
N = sum(w.*abs(ψ).^2)
```




50.816025132662084



and compare with the direct summation of coefficients:


```julia
Nsum = sum(abs(c).^2)
```




50.81602513266172




```julia
abs(N-Nsum)
```




3.623767952376511e-13



Integrals of higher order field products are computed with comparable accuracy.


```julia

```

e58f12cc2469fbe8f8ce4f9d1911d2d27604e1bf
ipynb
Jupyter Notebook
examples/quadrature methods.ipynb
AshtonSBradley/ProjectedGPE.jl
16623c1e00bbbae73e7448bd9b38f5b7e97979e9
MIT
2017-06-05T09:44:59.000Z
2021-01-15T17:59:48.000Z
examples/quadrature methods.ipynb
AshtonSBradley/ProjectedGPE.jl
16623c1e00bbbae73e7448bd9b38f5b7e97979e9
MIT
2018-08-30T05:46:58.000Z
2018-08-30T05:46:58.000Z
examples/.ipynb_checkpoints/projected GPE example-checkpoint.ipynb
AshtonSBradley/ProjectedGPE.jl
16623c1e00bbbae73e7448bd9b38f5b7e97979e9
MIT
2017-06-05T09:53:03.000Z
2017-06-05T09:53:03.000Z
Qwen/Qwen-72B
1. YES
2. YES