Repository files navigation

Briefly explained, this code is provided to solve specifically the most simplified case of the deuteron equation(s) and get its eigenfunction and eigenvalue.

In order to find the results, a particular model for the potential of the interaction is implemented. However, another option could be applied making use of the user's preferred election! 😄 Also, the process could be adjusted to another similar system for which the eigenfunction and eigenvalue want to be reached.

With the aim of better understanding what it is done in this project, let us introduce the main theoretical concepts.

The deuteron is theonly bound state of two nucleons and its bound is weak, therefore it only exists in the ground state. It is formed by one proton and one neutron and the nuclear force between them has the next properties: attractive (to form the bound state), short range (of order of$\sim 1fm$ ) and non-central (so it has a tensor component).

All these assumptions have been made through the obtained experimental data. Furthermore, the potential describing the interaction between both particles has been contructed more precisely to fit the actual results. It can be widely analysed and it will show many different components. We instead will focus just on themain quantity, avoiding the tentsorial part and other constituents. We will then have:

• Central potential$V_C$

This can also be described by different models. One option, the one we will use, is defining it through asquare well potential:

$$ V(r < r_0)=-V_0, \hspace{0.5cm} V(r \geq r_0)=0.$$

where$V_0$ and$r_0$ are determined again from experimental data. This will be the exact form of the central potential component$V_C$.

Now, let us derive the form of the Hamiltonian using deuteron's properties. The deuteron has total isospin$T=0$ and spin parity$J^{\pi}=1^+$. Since the spin of both proton and neutron is$s_n=s_p=1/2$, the total spin could in principle be$S=1,0$. Apart from that, in analogy with the Hydrogen atom, it is reasonable to assume that the ground state of the deuteron also has zero orbital angular momentum, so:$L=0$. At this point, due to the asymmetry of the whole wavefunction in a system of two fermions, we have$S=1$ for the total spin.

These results lead us from the total Hamiltonian

$$ H=-\frac{\hbar^2}{2M}\frac{d^2}{dr^2}r+\frac{\hbar^2}{M}\frac{L(L+1)}{r^2}+V_C(r)$$

to the simplified Schrödinger equation for the deuteron:

$$ -\frac{\hbar^2}{2M}\frac{d^2u_S}{dr^2}r+V_C(r)u_s(r)=Eu_S(r)$$

This is the equation that the code will solve! Therefore, we willget the deuteron's wavefunction$u_S$ and energy (by experimental data we know it is:$E=-2.225 MeV$, which is the same as its binding energy but with a negative sign).

As already said, this is actually the most simple model for this system. In fact, the potential would have a tensorial part, which generates a mixture between the eigenfunction$u_S$ for the s-wave ($L=0$) and$u_D$ for the d-wave ($L=2$). This could be the next step to improve this code!

To use this code properly, the next steps must be followed:

1. Be sure that all libraries needed are installed. For this process the next Python libraries are required:

• NumPy
• SciPy
• Matplotlib
• JSONLines
• Seaborn

2. Check that all the data we want to implement for the deuteron is correct ingenerate_data. 🔴 In case that any change is made, please run that file so that the new values will be saved invalues_deuteron and the code will use those new components.
3. If the potential model wants to be changed, look at the filepotential and modify carefully the function. If not, just continue to step 4.
4. We run therun main file and we get the desired eigenfunction$u_s$ and eigenvalue$E$ for deuteron plotted along a specified range of r! 👏

If the values and models are correctly applied, the execution will show a plot like this:

The code is divided in the next files:

• Ingenerate_data we insert the data we want to apply in the potential and the initial and boundary conditions. This will generate a data file calledvalues_deuteron in JSON lines format and it will allow us to read the data and work with it on an easy way.
• Inget_values we extract the data saved invalues_deuteron to later insert it straightfordwardly in the rest of the files.
• Inpotential we have the function defining the potential which will be put to model the deuteron interaction. This is my simplified implementation, but it could be replaced by the user's choice taking care of the form of the function, the inputs, outputs and the dimension of the parameters to be performed.
• Inschro_equation we express the second order differential equation in a system of two first order differential equations through the functionrad_equation. This system of 1st order ODEs can be solved with many methods fromscipy.integrate.
• Infind_energy we take the just explainedrad_equation function and insert it in thesolve_ivp SciPy method. After that, we compare the obtained eigenfunction's value in the last point of the radius and the actual boundary condition we are looking for. By doing so, we get the difference between those two values as the output of theerror_E function defined in this file.
• Inrun we finally take theerror_E function and apply it to a root finding process to obtain the energy value actually obeying the boundary condition. In this case the chosen method is called 'secant', which is defined inside theroot_scalar function of SciPy. After getting the proper eigenvalue, we can solve again our system described inschro_equation but this time with the correct E value. At this point, the code will just show the obtained result for the energy and the eigenfunction in a plot. We achieved the final result!
• Intests we test that each function of each file of the project works as it should and has the expected outputs.