Orthogonality of Spherical Bessel function J

4 Comments

The orthogonality of Bessel function for finite range is

\displaystyle \int_{0}^{a} J_{\nu}\left( \alpha_{\nu n} \frac{x}{a} \right) J_{\nu}\left( \alpha_{\nu m} \frac{x}{a} \right) x dx = \frac{a^2}{2} J_{\nu+1}^2(\alpha_{\nu n}) \delta_{nm}

where \alpha_{\nu n} is the n-th root of the Bessel function J_{\nu}. The above Bessel function is scaled to length a with n-th node. For example, a = 1, \nu = 0

Simply the integration by setting k_n = \alpha_{\nu n} / a

\displaystyle \int_{0}^{a} J_{\nu} (k_n x)J_{\nu}(k_m x) x dx = \frac{a^2}{2} J_{\nu+1}^2(\alpha_{\nu n}) \delta_{nm}


The relation between Riccati-Bessel function, spherical Bessel function, and Bessel function are:

\displaystyle \hat{j_l}(x) = x j_l(x), ~~ j_l(x) = \sqrt{\frac{\pi}{2x}} J_{l+1/2}(x)

The orthogonal condition of the spherical Bessel function should be

\displaystyle \int_{0}^{a} j_l( \kappa_n x) j_l(\kappa_m x) x^2 dx , ~~ \kappa_n = \frac{\alpha_{l+1/2, n}}{a}

From the Bessel Function orthogonality,

\displaystyle \int_{0}^{a} J_{l+1/2} (\kappa_n x)J_{l+1/2}(\kappa_m x) x dx = \frac{2}{\pi} \kappa_n \int_0^a x^2 j_l(\kappa_n x) j_l(\kappa_m x) dx

\displaystyle \int_0^a j_l (\kappa_n x) j_l( \kappa_m x) x^2 dx = \frac{a^3}{2} j_{l+1}^2(\alpha_{l+1/2, n} ) \delta_{nm}

In term of Riccati-Bessel function, \hat{j_l}(kx) = kx j_{l}(kx) ,

\displaystyle \int_0^a \hat{j_l} (\kappa_n x) \hat{j_l}( \kappa_m x) dx = \kappa_n^2 \frac{a^3}{2} j_{l+1}^2(\alpha_{l+1/2, n} ) \delta_{nm} = \frac{a}{2} \hat{j_{l+1}}(\alpha_{l+1/2, n}) \delta_{nm}

Using the orthogonality, the spherical Bessel function can be basis that span the whole space with f(a) = 0 . We can also normalized the spherical Bessel function. For Riccati-Bessel function, it span the whole space with f(0)=0, f(a)=0. Finite-range function expands as Bessel function is called Fourier-Bessel series.


The infinite range orthogonality of Bessel function J is (according to wiki and many reference, why?)

\displaystyle \int_0^\infty J_\alpha( u x) J_\alpha( vx) x dx = \frac{1}{u}\delta(u-v), ~~ \alpha > -\frac{1}{2}

in term of spherical Bessel function

\displaystyle \int_0^\infty j_{\alpha}(u x) j_\alpha(vx) x^2 dx = \frac{\pi}{2u^2} \delta(u-v), ~~ \alpha > -1

and using \hat{j_l}(kx) = kx j_{l}(kx) , in term of Riccati-Bessel function is

\displaystyle \int_0^\infty \hat{j_{\alpha}}(u x) \hat{j_\alpha}(vx) dx = \frac{\pi}{2} \delta(u-v), ~~ \alpha > -1

For example, the spherical Bessel J of order zero,

\displaystyle j_0(x) = \frac{\sin(x)}{x} ,

Let’s integrate it

\displaystyle \int_0^\infty \frac{\sin(x)}{x} \frac{\sin(x)}{x} x^2 dx = \int_0^\infty \sin^2(x) dx

which is infinite, and the right-hand side is a delta function \delta(0) = \infty .

The problem with the above formulae is, that the orthogonal condition is OK. but it should be also normal, i.e. orthonormal, to be a normalizable basis.


The Hankel Transform is

\displaystyle F_\nu(k) = \int_0^\infty f(r) J_\nu(kr) r dr , ~~ \nu > -\frac{1}{2}, k \geq 0

and this can change to spherical-Bessel function as,

\displaystyle \psi (k) = \sqrt{\frac{2}{\pi}}\int_0^\infty \psi(r) j_l(kr) r^2 dr

The foundation for the transformation is the orthogonality of the Bessel function. Substitute f(r) as Bessel function.

\displaystyle \int_0^\infty J_0(h r) J_0(kr) r dr =\frac{1}{h} \delta (k-h)

So, what do I miss?

For example, J_\nu(x) is oscillating around zero, the square of it is always positive, and with a weighting x, it approach to \sin^2(x) with a phase shift,

\displaystyle J_\nu^2(x) x \xrightarrow{x\rightarrow \infty} \frac{2}{\pi} \sin^2(x + \theta_\nu)

And the integrate of \sin^2(x) is \pi/2 . Thus, the integral of Bessel function with weight x is an infinite sum for between each zeros, and each terms is ~1. How this infinite sum becomes 1?


Recall the plane wave expansion,

\displaystyle \exp(i \vec{k}\cdot\vec{r}) = \sum_{l=0}^\infty (2l+1) i^l j_l(kr) P_l(\hat{k}\cdot\hat{r})

As we show it before, it is a special case of the Hankel transform, or belong the general Fourier Transform.

Wave function in momentum space

Leave a comment

The wave function often calculated in spatial coordinate. However, in experimental point of view, the momentum distribution can be extracted directly from the experimental data.

The conversion between momentum space and position space is the Fourier transform

\displaystyle \phi(\vec{k}) = \frac{1}{\sqrt{2\pi}^3} \int Exp\left(-i \vec{k}\cdot \vec{r} \right) \phi(\vec{r}) d\vec{r}

Using the plane wave expansion

\displaystyle \exp(i k\cdot r) = \sum_{l=0}^\infty (2l+1) i^l j_l(kr) P_l(\hat{k}\cdot\hat{r})

or

\displaystyle \exp(i k\cdot r) = 4\pi \sum_{l=0}^\infty \sum_{m=-l}^{l} i^l j_l(kr) Y_{lm}(\Omega_k) Y_{lm}^{*}(\Omega_r)

Thus,

\displaystyle \phi(\vec{k}) = \frac{4\pi}{\sqrt{2\pi}^3} \sum_{l=0}^\infty (-i)^l \sum_{m=-l}^{l} \int j_l(k r) Y_{lm}(\Omega_k) Y_{lm}^*(\Omega) \phi(\vec{r}) r^2 dr d\Omega

where j_l (x) is spherical Bessel function. Usually, due to conservation of angular momentum, the angular part can be separated from the spatial part. Let assume the wave function in position space is

\phi(\vec{r}) = \psi(r) Y_{l_r m_r}(\Omega)

\phi(\vec{k}) = \psi(k) Y_{l_k m_k}(\Omega_k)

Then we have

\displaystyle \psi(k) = \frac{4\pi}{\sqrt{(2\pi)^3}} \int j_l(k r) \psi(r) r^2 dr

\displaystyle \phi(\vec{k}) = \psi(k) (-i)^l Y_{lm}(\Omega_k)

where l = l_r = l_k, m=m_r = m_k , due to the orthogonality of spherical harmonics.

For s, p, d, f-state, the spherical Bessel function is

\displaystyle j_0(x) = \frac{\sin(x)}{x}

\displaystyle j_1(x) = \frac{\sin(x)}{x^2} - \frac{\cos(x)}{x}

\displaystyle j_2(x) = \sin(x)\frac{2-x^2}{x^3} - \cos(x)\frac{3}{x^2}

\displaystyle j_3(x) = \sin(x)\frac{15-6x^2}{x^4} - \cos(x)\frac{15-x^2}{x^3}

For Hydrogen-like wave function, the non-normalized momentum distribution is

\displaystyle \psi_{10}(k) =  \frac{4 Z^{5/2}}{(k^2 + Z^2)^2}

\displaystyle \psi_{20}(k) = \frac{16 \sqrt{2}Z^{5/2}(4k^2-Z^2)}{(4k^2 + Z^2)^3}

\displaystyle \psi_{21}(k) =\sqrt{\frac{2}{3}}\frac{64 k Z^{7/2}}{(4k^2 + Z^2)^3}

\displaystyle \psi_{30}(k) =  \frac{36 \sqrt{3} Z^{5/2} (81k^3-30k^2Z^2+Z^4)}{(9k^2 + Z^2)^4}

\displaystyle \psi_{31}(k) = \frac{144 \sqrt{6} Z^{7/2} (9k^3-kZ^2)}{(9k^2 + Z^2)^4}

Here is the plot for momentum distribution (\psi_{nl}(k))^2 k^2 .

k_1.PNG

k_2.PNG

k_3.PNG

It is interesting that, the number of node decrease with higher angular momentum. But be-aware that it is only in atomic case, not a universal.

The higher the principle quantum number, the smaller of the spread of momentum. This is because, the spread of position wave function getting larger, and the uncertainty in momentum space will be smaller. This is a universal principle.

We also plot the Hydrogen radial function in here \psi(r)^2 x^2 , for reference,

r_1.PNG

r_2.PNG

r_3.PNG

 

 

A little summery of what I am doing

Leave a comment

My current position is developing a computational program that can measure the system of polarized target automatically and repeatedly. The program needs to connect with the microwave generator, voltage supply, power meter, and oscilloscope. That is purely technical and nothing special.

Later, after the program started to collecting data, another program is needed to analysis the data. Although there are many program that can do analysis, but those program are not so easy to use, in both control and display. So, I make an analysis program. The program is simply read the 2-D time-domain signal, and then determine some parameters of a specific function that fit the signal. However, the signal could be very noisy, so FFTW and wavelet analysis were implemented. That’s why the wavelet analysis appeared in this blog.

After that, the sample of polarized target is made from various pentacene derivatives, that no known energy level exist. One way to find out the energy levels of the singlet excited states is to measure the absorption spectrum. However, the energy of the triplet state is difficult to measure. And the energy level of the triplet state is critical for Dynamic Nuclear Polarization to be happened. Thus, one way to find out is doing computation chemistry.

My last chemistry class was like 15 years ago, when I was junior high school. But the basic of computation chemistry is solving the Schrodinger equation. That is what theoretical nuclear physicists do! The variation method, the Hartree-fock method, I heard it and somehow know it, but never do it with my own hand and computer. That is why I revisit the Hartree-fock method, and found out my previous understanding is so naive.

In the course of studying Hartree-Fock method, and one problem is evaluating the overlap integral

\displaystyle \int \psi_a(r_1) \psi_b(r_2) \frac{1}{r_{12}} \psi_c(r_1) \psi_d(r_2) dr_1dr_2

For a 3-D system, the integral involves product of multiple spherical harmonics. That is really troublesome. Therefore, I move to study the spherical harmonics, and the related rotational invariant, Wigner D-matrix, Clebcsh-Gordon series and Fourier series.  The spherical harmonics arises from solving the Laplace equation in spherical coordinate. A general theory of the solution of Lapalace equation involves Legendre polynomial, which is a special case for Hypergeometric function.  And a very interesting connection is that the elliptical function of the 1st and 2nd kind are also two special cases of hypergeometric function, that, the solution of Laplace equation for elliptical boundary condition is elliptical functions. That connects spherical harmonics and elliptical function! Wow!

I am now a bit off-track, that I am very interesting on the function of all functions. We know that there are many elementary functions, such as sin, cos, and Log. And even more kind of special functions, such as

  • Hermite — solving 1-D harmonic oscillator
  • Laguerre — the radial function of hydrogen
  • Legendre — the solution of the “\theta” of  Laplace equation
  • Gamma — a continuation of factorial
  • Bessel — solution of 3-D infinite square well
  • Elliptic — magnetic field of a solenoid
  • Dirac delta
  • Gaussian or Error function

For discrete argument

  • Clebsch-Gordon
  • Factoral
  • Binomial

As far as I know, the Hypergemetric function is like “mother of functions”, although not all special functions can be expressed as it. I am driven by curiosity, so, I am not sure where I will go. For instance, the transformation of Hypergeometric function is very interesting.

So, for now, as an ending, I found one article is very interesting. The History and Future of Special Functions, by Stephen Wolfram.

Review on rotation

Leave a comment

The rotation of a vector in a vector space can be done by either rotating the basis vector or the coordinate of the vector. Here, we always use fixed basis for rotation.

For a rigid body, its rotation can be accomplished using Euler rotation, or rotation around an axis.

Whenever a transform preserves the norm of the vector, it is a unitary transform. Rotation preserves the norm and it is a unitary transform, can it can be represented by a unitary matrix. As a unitary matrix, the eigen states are an convenient basis for the vector space.

We will start from 2-D space. Within the 2-D space, we discuss about rotation started by vector and then function. The vector function does not explicitly discussed, but it was touched when discussing on functions. In the course, the eigen state is a key concept, as it is a convenient basis. We skipped the discussion for 3-D space, the connection between 2-D and 3-D space was already discussed in previous post. At the end, we take about direct product space.


In 2-D space. A 2-D vector is rotated by a transform R, and the representation matrix of R has eigen value

\exp(\pm i \omega)

and eigenvector

\displaystyle \hat{e}_\pm = \mp \frac{ \hat{e}_x \pm i \hat{e}_y}{\sqrt{2}}

If all vector expand as a linear combination of the eigen vector, then the rotation can be done by simply multiplying the eigen value.

Now, for a 2-D function, the rotation is done by changing of coordinate. However, The functional space is also a vector space, such that

  1. a* f_1 + b* f_2 still in the space,
  2. exist of  unit and inverse of addition,
  3. the norm can be defined on a suitable domain by \int |f(x,y)|^2 dxdy

For example, the two functions \phi_1(x,y) = x, \phi_2(x,y) = y , the rotation can be done by a rotational matrix,

\displaystyle R = \begin{pmatrix} \cos(\omega) & -\sin(\omega) \\ \sin(\omega) & \cos(\omega) \end{pmatrix}

And, the product x^2, y^2, xy also from a basis. And the rotation on this new basis was induced from the original rotation.

\displaystyle R_2 = \begin{pmatrix} c^2 & s^2 & -2cs \\ s^2 & c^2 & 2cs \\ cs & -cs & c^2 - s^2 \end{pmatrix}

where c = \cos(\omega), s = \sin(\omega) . The space becomes “3-dimensional” because xy = yx, otherwise, it will becomes “4-dimensional”.

The 2-D function can also be expressed in polar coordinate, f(r, \theta) , and further decomposed into g(r) h(\theta) .


How can we find the eigen function for the angular part?

One way is using an operator that commutes with rotation, so that the eigen function of the operator is also the eigen function of the rotation. an example is the Laplacian.

The eigen function for the 2-D Lapacian is the Fourier series.

Therefore, if we can express the function into a polynomial of r^n (\exp(i n \theta)  , \exp(-i n \theta)) , the rotation of the function is simply multiplied by the rotation matrix.

The eigen function is

\displaystyle \phi_{nm}(\theta) = e^{i m \theta}, m = \pm

The D-matrix of rotation (D for Darstellung, representation in German)  \omega is

D^n_{mm'}(\omega) = \delta_{mm'} e^{i m \omega}

The delta function of m, m' indicates that a rotation does not mix the spaces. The transformation of the eigen function is

\displaystyle \phi_{nm}(\theta') = \sum_{nm} \phi_{nm'}(\theta) D^n_{m'm}(\omega)

for example,

f(x,y) = x^2 + k y^2

write in polar coordinate

\displaystyle f(r, \theta) = r^2 (\cos^2(\theta) + k \sin^2(\theta)) = \frac{r^2}{4} \sum_{nm} a_{nm} \phi_{nm}(\theta)

where a_0 = 2 + 2k, a_{2+} = a_{2-} = 1-a, a_{other} = 0.

The rotation is

\displaystyle f(r, \theta' = \theta + \omega ) = \frac{r^2}{4} \sum_{nm} a_{nm} \phi_{nm}(\theta) D^n_{mm}(\omega)  = \frac{r^2}{4} \sum_{nm} a_{nm} \phi_{nm}(\theta + \omega)

If we write the rotated function in Cartesian form,

f(x',y') = x'^2 + k y'^2 = (c^2 + k s^2)x^2 + (s^2 + k c^2)y^2 + 2(k-1) c s x y

where c = \cos(\omega), s = \sin(\omega) .


In 3-D space, the same logic still applicable.

The spherical harmonics Y_{lm} serves as the basis for eigenvalue of l(l+1), eigen spaces for difference l are orthogonal. This is an extension of the 2-D eigen function \exp(\pm n i \theta) .

A 3-D function can be expressed in spherical harmonics, and the rotation is simple multiplied with the Wigner D-matrix.


On above, we show an example of higher order rotation induced by product space. I called it the induced space (I am not sure it is the correct name or not), because the space is the same, but the order is higher.

For two particles system, the direct product space is formed by the product of the basis from two distinct space (could be identical space).

Capture.PNG

Some common direct product spaces are

  • combining two spins
  • combining two orbital angular momentum
  • two particles system

No matter induced space or direct product space, there structure are very similar. In 3-D rotation, the two spaces and the direct product space is related by the Clebsch-Gordon coefficient. While in 2-D rotation, we can see from the above discussion, the coefficient is simply 1.

Lets use 2-D space to show the “induced product” space. For order n=1, which is the primary base that contains only x, y.

For n=2, the space has x^2, y^2, xy, but the linear combination x^2 + y^2 is unchanged after rotation. Thus, the size of the space reduced 3-1 = 2.

For n = 3, the space has x^3, y^3, x^2y, xy^3 , this time, the linear combinations x^3 + xy^2 = x(x^2+y^2) behave like x and y^3 + x^2y behave like y, thus the size of the space reduce to 4 - 2 = 2.

For higher order, the total combination of x^ay^b, a+b = n is C^{n+1}_1 = n+1 , and we can find n-1 repeated combinations, thus the size of the irreducible space of order n is always 2.

For 3-D space, the size of combination of x^ay^bz^c, a + b+ c = n is C^{n+2}_2 = (n+1)(n+2)/2 . We can find n(n-1)/2 repeated combination, thus, the size of the irreducible  space of order n is always 2n+1.

Hard & soft thresholding

Leave a comment

In usual Fourier transform (FT), the filter is cut-off certain frequency.

This trick is also suitable for wavelet transform (WT). However, there could be some “features” located in high frequency scale (or octave) , a simply cut-off would remove these features.

If the signal to noise level is large, that means the noise has smaller amplitude than that the signal, we can use hard or soft thresholding, which zero any coefficient, which is after the FT or WT,  less then a threshold.

Lets X be the coefficient. The hard thresholding is

Y=\begin{cases} 0, & |X| <\sigma \\ X, & \mbox{else} \end{cases}

ht

The soft thresholding is

Y = \begin{cases} 0, & |X| < \sigma \\ sign(X) f(|X|, \sigma), & \mbox{else} \end{cases}

A popular function

\displaystyle f(x, \sigma) = \frac{x - \sigma}{ X_{max} - \sigma } X_{max}

st.PNG

or

\displaystyle f(x, \sigma) = x - \sigma

st2.PNG

Wavelet Analysis or MRA

2 Comments

Although the Fourier transform is a very powerful tool for data analysis, it has some limit due to lack of time information. From physics point of view, any time-data should live in time-frequency space. Since the Fourier transform has very narrow frequency resolution, according to  uncertainty principle, the time resolution will be very large, therefore, no time information can be given by Fourier transform.

Usually, such limitation would not be a problem. However, when analysis musics, long term performance of a device, or seismic survey, time information is very crucial.

To over come this difficulty, there a short-time Fourier transform (STFT) was developed. The idea is the applied a time-window (a piecewise uniform function, or Gaussian) on the data first, then FT. By applying the time-window on difference time of the data (or shifting the window), we can get the time information. However, since the frequency range of the time-window  always covers the low frequency, this means the high frequency  signal is hard to extract.

To improve the STFT, the time-window can be scaled (usually by 2). When the time window is shrink by factor of 2, the frequency range is expanded by factor of 2. If we can subtract the frequency ranges for the time-window and the shrink-time-window, the high frequency range is isolated.

To be more clear, let say the time-window function be

\phi_{[0,1)}(t) = 1 , 0 \leq t < 1

its FT is

\hat{\phi}(\omega) = sinc(\pi \omega)

Lets also define a dilation operator

Df(t) = \sqrt{2} f(2t)

the factor \sqrt{2} is for normalization.

The FT of D\phi(t) has smaller frequency range, like the following graph.

Capture.PNG

We can subtract the orange and blue curve to get the green curve. Then FT back the green curve to get the high-frequency time-window.

We can see that, we can repeat the dilation, or anti-dilation infinite time. Because of this, we can drop the FT basis Exp(-2\pi i t \omega), only use the low-pass time-window to see the low-frequency behaviour of the data, and use the high-pass time-window to see the high-frequency behaviour of the data. Now, we stepped into the Multi-resolution analysis (MRA).

In MRA, the low-pass time-window is called scaling function \phi(t) , and the high-pass time-window is called wavelet \psi(t).

Since the scaling function is craetd by dilation, it has the property

\phi(t) = \sum_{k} g_{0}(k) \phi(2t-k)

where k is integer. This means the vector space span by {\phi(t-k)}_{k}=V_0 is a subspace of the dilated space DV_0 =V_1. The dilation can be go one forever, so that the whole frequency domain will be covered by V_{\infty}.

Also, the space span by the wavelet, {\psi(t-k)}=W_0, is also a subspace of V_1. Thus, we can imagine the structure of MRA is:

Capture.PNG

Therefore, any function f(t) can also be expressed into the wavelet spaces. i.e.

f(t) = \sum_{j,k} w_{j,k} 2^{j/2}\psi(2^j t - k)

where j, k are integers.

I know this introduction is very rough, but it gives a relatively smooth transition from FT to WT (wavelet transform), when compare to the available material on the web.

 

 

 

 

Qt with FFTW

2 Comments

To use Qt with FFTW (Fast Fourier Transfrom in the West) in windows.

  1. download the FFTW pre-compiled *.dll file
    http://www.fftw.org/install/windows.html
  2. extract the zip. Copy the fftw3.h and libfftw3-3.dll and put them into your project folder, the folder same as *.pro file.
  3. in *.pro, add
    LIBS += "$$PWD/libfftw3-3.dll"
  4. in QT creator, on the “projects” column, right click on the project name and select “Add Existing Files…”, choose the fftw.h
  5. in main.cpp, add
    #include "fftw.h"
  6. got to the build directory, usually named as “Build-…..”, you will see “debug” or “release” folders. copy the libfftw3-3.dll and paste it in the “debug” or “release” folder depends on you are debug or release.

The step 6 is very important, without so, the program will not even started, while the compilation is completely normal.

[Pol. p target] laser duty 20% and 30%

Leave a comment

This were what we doen today.

the peaks of Fourier spectrum position seem depend on crystal angle.

[Pol. p target] laser duty 10%

Leave a comment

we did 10% laser duty. the result is no surprise.

since we employed Fourier Analysis, so, the Fourier spectrum has 2 peaks, separated by 60 kHz. the central frequency is 12.6 MHz.

Discrete Fourier Transform (shift of the spectrum)

Leave a comment

as the last post shown the formula of discrete Fourier transform:

a_n = \sum_{i=0}^m { x_i e^{i (n-1) 2\pi (i-1) /m}}

where (m+1) is the number of data-point, n is from 1 to (m+1). if n is bigger then m+1, it is a_n is start to repeat. for example, n = m+1+r

a_{m+1+r} = a_{r}

Therefore,

a_{m+1-r} = a_{-r}

Thus, we can rearrange the a_n from -(m+1)/2  to (m+1)/2. This, set the center be zero-frequency. for example, if m+1 = 1000, then we can have n from 1 to 1000, or -499 to 0 and 1 to 500. since n=1 mean frequency = 0, we can shift n by 1 unit and get n from -500 to 499.

this it very useful if out signal is actually from 2 signals, 1 is raw signal and the other is a reference signal. the signal we get is by subtracting these 2. thus, when the raw signal is very close to the reference signal, the frequency is very small and the peak of Fourier spectrum will go to the edge and hardly recognized. However, by shift the spectrum and set the a_0 at the center. we can see How the raw signal is different from the reference signal.

 

Older Entries