Coulomb Energy in Nuclear Unit

Leave a comment

Coulomb energy between 2 charge q_1, q_2 in SI unit is

\displaystyle U_c = \frac{1}{4\pi\epsilon_0} \frac{q_1 q_2}{r}  [J]

\displaystyle U_c = \frac{e^2}{4\pi\epsilon_0} \frac{Z_1 Z_2}{r [m]} [J \cdot m] = \frac{Z_1 Z_2}{r [m]} (\frac{e^2}{4\pi \epsilon_0}) [J \cdot m] = \frac{Z_1 Z_2}{r [m]} 2.30708 \times 10^{-28} [J \cdot m]

We need to convert the SI unit into nuclear unit:

1 [J] = \frac{10^{-6}}{e} [MeV]

1 [m] = 10^{15} [fm]

Since the unit of the r depends on the nominator, the change of the unit does not need to be compensated.

\displaystyle U_c = \frac{Z_1 Z_2}{r [fm]} 1.43996 [MeV \cdot fm]

Therefore, a simple expression

\displaystyle U_c = \frac{e^2}{r} Z_1 Z_2

where e^2 = 1.44 [MeV\cdot fm]

Other useful quantities are:

  • \hbar c = 197.327 [MeV\cdot fm]
  • e^2/\hbar c = 1/137.036
  • \hbar = 6.58212 [MeV\cdot s]
  • c = 2.99792458\times 10^{23} [fm/s]

Qt with FFTW

Leave a comment

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.

Electron spin resonance in power sample

Leave a comment

The topic was discussed in a paper “Dynamic nuclear polarization by photoexcited-triplet electron spins in polycrystalline samples” by K. Takeda, K. Takegoshi, and T. Terao. I am going to give some supplementary on how the calculation can be done.


In a spin-1 system, the spin operator S_x, S_y, S_z and the rotation operator are easy to calculate.

The zero-field splitting Hamiltonian in the ordinary basis of S_z is

\displaystyle H_{zfs} = \begin{pmatrix} D/3 & 0 & E \\ 0 & -2D/3 & 0 \\ E & 0 & D/3 \end{pmatrix}

The eigen vectors in the paper are

\displaystyle \left|X\right> = \frac{1}{\sqrt{2}}( \left|-\right> - \left|+\right> )

\displaystyle \left|Y\right> = \frac{i}{\sqrt{2}}( \left|-\right> + \left|+\right> )

\displaystyle \left|Z\right> = \left|0\right>

with eigen values X = D/3 -E, Y = D/3 +E, and Z = -2D/3 . The convention that the \left| Y \right> has an imaginary number i is not necessary. In the paper, there are typos on the eigen values.

In other way, the zero-field splitting Hamiltonian can be diagonalized into

\displaystyle H_{zfs} =  P\cdot H_{XYZ} \cdot P^{-1}

\displaystyle H_{XYZ} = \begin{pmatrix} X & 0 & 0 \\ 0 & Y & 0 \\ 0 & 0 & Z \end{pmatrix}

\displaystyle P = \frac{1}{\sqrt{2}} \begin{pmatrix} -1 & i & 0 \\ 0 & 0 & \sqrt{2} \\ 1 & i & 0 \end{pmatrix}

where P is the transformation changing from ordinary basis to eigen basis.


Given that the population of the eigen basis is w = ( w_X, w_Y, w_Z) , the density matrix is

\displaystyle \rho_{XYZ} = \sum_{i = X, Y, Z} \left| i \right>\left< i \right| w_i

The corresponding population in ordinary basis is the diagonal elements of

\rho = P \cdot \rho_{XYZ} \cdot P^{-1}

In the paper, it gives w = (0.76, 0.16, 0.08). In order to get the population when the molecule X-axis is aligned with the Z-axis Lab frame, the \rho has to be rotated by 90 degree, and we have the population w' = ( 0.12, 0.76, 0.12) .


The Zeeman Hamiltonian is H_B = \omega S_Z , where \omega = \gamma_e  B  \hbar . In the paper, the author fix the crystal and rotating the external magnetic field. To express the H_B into eigen basis. We simply apply a transformation as the density matrix. But first we need to rotate the magnetic field. The rotation is simple, as the direction of the magnetic field is (\theta, \phi) . Thus, the Zeeman Hamiltonian is H_B = \omega \hat{n} \cdot \vec{S} .

The transformation is

H_{B_XYZ} = P^{-1} \cdot H_B \cdot P

again, the total Hamiltonian is difference from the author.

\displaystyle H = \begin{pmatrix} X & -i \omega cos(\theta) & i \omega cos(\phi) sin(\theta) \\ i \omega cos(\theta) & Y & -i\omega sin(\phi) sin(\theta) \\ -i\omega cos(\phi) sin(\theta) & i \omega sin(\phi) sin(\theta) & Z \end{pmatrix}

The difference is because the normal vector \hat{n} = ( sin(\phi) sin(\theta), cos(\phi) sin(\theta), cos(\theta) ) is used in the paper.

After that, the calculation for the power sample is the same.


Instead of using integration, we can generated the power using isotropic distribution. I generated 10k molecule, the result of 9.65 GHz microwave is

power_ESR.PNG

where the y axis is the magnetic field in Tesla. The result is similar to the paper. I think that their calculation somehow changed the eigen value, so that the eigen vectors were matching with the eigen energy. And the notation of the normal vector and the constant factor of the eigen vectors does not matter.

Phase shift of elastics scattering

Leave a comment

I found that the derivation of most “google result” is not clear enough. So here is my derivation. Before process, people may need to review the pervious post.

Most people start from,

u_l = A( \hat{h}_l^- - S_l \hat{h}_l^+ )

that annoying me because it is somehow not “natural”. Why there is a “minus” sign? Why the \hat{h}_l^- is the first term? For my self, a more natural way is,

u_l = a \hat{h}_l^+ + b \hat{h}_l^-

where a, b are complex numbers, but that is still not so natural, because in numerical calculation, for simplicity, there is no complex number, we only have,

u_l = \alpha \hat{j}_l + \beta \hat{n}_l

The first term is alway there as it is the free solution and bounded at r = 0. the second term is caused by the potential.


The goal is to find a solution take the form

\displaystyle \psi = A \left( e^{i \vec{k} \cdot \vec{r}} + f(\theta) \frac{e^{ikr}}{r} \right)

where the first term is free wave and the second term is scattered wave. The solution for elastics scattering is

\displaystyle \psi = \sum C_l P_l (\cos\theta) \frac{u_l}{kr} \rightarrow \sum C_l P_l(\cos\theta) (\alpha \hat{j}_l + \beta \hat{n}_l)

we used the substitution,

\displaystyle R_l(r) = \frac{u_l(\rho)}{\rho},  \rho = kr.

The radial function can be solved using Rungu-Kutta method on the equation,

\displaystyle \frac{d^2}{d\rho^2} u_l = \frac{2 m_\mu}{\hbar^2} (V-E) u_l + \frac{l(l+1)}{\rho^2}

and the solution of u_l at far away is,

u_l \rightarrow  \alpha \hat{j}_l + \beta \hat{n}_l.

the arrow means r \rightarrow \infty. So, the problem is how to rewrite the solution. In the way, we will see how the phase shift or the S-matrix was found.


The free solution is the spherical wave,

\displaystyle e^{i \vec{k} \cdot \vec{r}} = \sum_l (2l+1) i^l P_l(\cos\theta) j_l(kr)

The spherical Bessel function j_l(kr) cna be express as Heankel function

h_l^{\pm} = n_l \pm i j_l \rightarrow e^{\pm i (kr - l \frac{\pi}{2})}

The + sign is outgoing wave.


\displaystyle u_l \rightarrow (\alpha \hat{j}_l + \beta \hat{n}_l)

\displaystyle = \frac{\alpha}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \frac{\beta}{2}(\hat{h}_l^{+} + \hat{h}_l^{-})

\displaystyle = \frac{\alpha + i \beta}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \beta \hat{h}_l^{+}

\displaystyle = (\alpha - i \beta ) \left( \frac{\hat{h}_l^+ - \hat{h}_l^-}{2i} + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

\displaystyle = (\alpha - i \beta ) \left( \hat{j}_l + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

Since the u_l should be normalized, we can se \alpha = \cos \delta and \beta = \sin\delta.

\displaystyle \frac{\beta}{\alpha - i \beta } = \sin(\delta) e^{i\delta}

We put u_l back

\displaystyle \psi \rightarrow \sum_l C_l P_l (cos\theta)(\alpha - i \beta ) \left( \hat{j}_l + \sin(\delta) e^{i\delta} \hat{h}_l^+\right)

By setting

\displaystyle C_l = A i^l \frac{2l+1}{\alpha - i \beta} ,

we have the first term is the free wave function. In the second term, \hat{h}_l^+ \rightarrow e^{i(kr - l \frac{\pi}{2}}) / kr . Notice that

e^{i l \frac{\pi}{2}} = i^{-l}

That cancel the i^l term in C_l . And we have

 \displaystyle f(\theta) = \sum (2l+1) P_l (\cos\theta) \frac{\sin(\delta) e^{i\delta}}{k}


some people will write the u_l as \hat{h}_l^{\pm} and the S-matrix,

\displaystyle u_l = \frac{\alpha + i \beta} {2i} \hat{h}_l^+ - \frac{\alpha - i \beta}{2i} \hat{h}_l^-

\displaystyle = -\frac{\alpha - i \beta}{2i} \left( \hat{h}_l^- - \frac{\alpha + i \beta}{\alpha - i \beta} \hat{h}_l^+ \right)

\displaystyle = A' (\hat{h}_l^- - S_l \hat{h}_l^+)

where

\displaystyle S_l =\frac{\alpha + i \beta}{\alpha - i \beta} = e^{2i\delta} .

Remember that this is the S-matrix for elastics scattering.

 

Cross section IV

Leave a comment

A more microscopic view of the cross section is

\displaystyle  \frac{dN}{dS} \sigma \rho dt dA= dn,

where N is the number of incident particle, S is the area of the beam, \sigma is the total cross section, \rho is the particle density in the target, A is the area of the target, t is the thickness of the target, and n is the particle detected.

The cross section is a constant, so, after integration on the target area and thickness,

\displaystyle \sigma = \frac{n}{N \rho t} .

If the beam is not uniform, dN/dS = f(S) is a function of S. The integration has to be careful.

 

A reminder of Exponential function

Leave a comment

Suddenly found that exponential function is very wonderful. It is the eigen-function of differential operator. because of that, it appears in almost everywhere as differential operator constantly appears.

The definition of the exponential e is,

\displaystyle e = \lim_{n\to\infty} \left(1+\frac{1}{n} \right)^n

this is the 1% growth continually in a unit of time. It is easy to see that

\displaystyle e^x =  \lim_{n\to\infty} \left(1+\frac{x}{n} \right)^n

e^x e^y = e^{x+y}

expanding,

\displaystyle \left(1+\frac{x}{n} \right)^n = \sum_{r=0}^{n} C_r^n \left(\frac{x}{n} \right)^r  = \sum_{r=0}^{n} \left(\prod_{k=0}^{r-1} \frac{n-k}{n} \right)\frac{x^r}{r!},

Take n\to\infty,

\displaystyle e^x = \sum_{r=0}^{\infty} \frac{x^r}{r!}

It is easy to see that

\displaystyle \frac{d}{dx} e^x = e^x

 

Levenberg-Marquardt Algorithm

Leave a comment

In pervious post, we shows the Gauss-Newton method for fitting non-linear function. The disadvantage of that method is that the inverse matrix could be ill-defined. This makes the method unstable.

Back to the basic, we want to minimize the sum of square of residual (SSR). The SSR is,

SSR(\beta) = (Y - f(\beta))^T\cdot (Y-f(\beta))

The derivative on \beta,

\frac{d}{d\beta} SSR(\beta) = -2 (Y-f(\beta))^T \cdot \nabla f(\beta)

Many literatures denote \nabla f = J, which is the Jacobian. The second derivative of f is Hessian matrix H = \nabla^2 f \sim J^T\cdot J.

The Gradient Descent method is that ,

h = \beta - \beta_0 = \alpha J^T \cdot (Y - f(\beta_0))

where \alpha is a step size. The gradient descent changes the SSR using the steepest path. The step size \alpha has to be adjusted. The simplest way to adjust is testing the \delta = SSR(\beta_0 + h) - SSR(\beta_0). If \delta < 0 , the \alpha increases, else decreases. This method is slow but stable. It is slow because of finding the \alpha . It is stable because the method is always computable.


Thus, we have 2 methods, Gradient Descent is stable and slow, Gauss-Newton method is unstable but fast. Levenberg-Marquardt Algorithm combined this 2 methods so that it is stable and fast by solving,

(J^T \cdot J + \lambda I) h = J^T \cdot (Y - f)

where \lambda is an adjustable parameter. When \lambda >> 1 , the J^T\cdot J is neglectable and the method becomes Gradient Descent with small \alpha. When the \lambda << 1, the method becomes Gauss-Newton method.

Usually, the \lambda_0 is small. The Gauss-Newton method is very good near the minimum of SSR, while Gradient Descent is better far away.

When the \delta < 0, \lambda_{i+1} = \lambda_i / 10 , else \lambda_{i+1} = \lambda_i * 10. I don’t know the exact reason for this setting. In fact if you set oppositely, the method is still work in most cases.

The method add \lambda I on the J^T\cdot J , the inverse is always well-define. Therefore, this method is stable.

Older Entries Newer Entries