Briefly explained, this code is provided to solve specifically the deuteron equations and get its eigenfunctions and eigenvalue.
In order to find the results, a particular model for the potential is implemented. However, another option could be applied, making use of the user's preferred election.
With the aim of better understanding what it is done in this code, let us introduce the main theoretical concepts.
The deuteron is the only 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 main properties:
- Attractive: forms a bound state
- Short range: of order of
$\sim 1fm$ - Non-central: it has a tensor component connecting different states
- Hard core: there is a range in which the potential is completely repulsive, so that the nucleons do not collapse
- Conserves parity
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 on two main quantities:
- Central potential
$V_C$ - Tensorial potential
$V_T$
Both can also be described by diverse phenomenological models. One option (the one we will use) is to define them through Yukawa potential wells:
where
Now, let us derive the form of the Hamiltonian using some of the deuteron's properties. The deuteron has total isospin
These results lead us from the total Hamiltonian
into two radial equations describing the system which, after inserting the discussed characteristics and quantum numbers, have the next form:
These are exactly the equations that our code will try to solve! So, we will get the two wavefunctions
Also, we must take into account the initial and boundary conditions we expect for this kind of system. Exposed in a simple form:
Apart from that, we will be able to also calculate the electric quadrupole moment Q. This term indicates that the charge distribution is not spherically symmetric. This is why we must take into account the d-wave
The value of Q is obtained through its definition:
which can be further developed by introducing the operator form as:
where
This is exactly the formula used in the code to obtain the result of the quarupole moment, which will also work as a proof of the correctness of the model used.
Let us describe the way in which this code works and what is each module used for.
Firstly, we will make a brief description of each module following slightly an order of usefulness, so that the explanations can be better understood:
-
generate_data: here we insert the values for our model. We will implement the initial conditions and boundary conditions for both
$u_s$ and$u_d$ . We insert two different combinations, because we will use a combination of two independent solutions to reach the final result obeying all conditions (later explained). Also the range of the radius will be present, and the initial guess for the energy eigenvalue. All the data will be saved in deuteron_values. - get_values: we can get all values implemented using this module. It takes all the different data from the data files.
- potentials: contains the functions describing both potentials we have defined in the theory.
- plotVs: we get the image of both potentials to better imagine the behaviour and interaction between particles.
- equations: this file contains the form of the system we must solve. We take the 2 second order differential equations and create a system of 4 first order differential equations.
- find_results: we define the root functions that we will use in root finding processes implemented in other modules to find the values of A, B, C and D constants and also the E eigenvalue.
- get_solutions: using functions in find_results, we apply root finding processes to reach first the proper values of A, B, C and D giving the continuity of the eigenfunctions in a selected midpoint. Later, with the obtained values, we apply again a root finding function to reach the appropiate value for E energy. All the results will be saved in values_ABCDE and they will be extracted through get_values to be used in other files.
- get_graph: we write a function which solves the equations with the previously obtained correct values.
-
plotWFs: with this module we finally reach the form of the
$u_s$ and$u_d$ wavefunctions. We make use of get_graph to solve the integration. The two functions are plotted vs.$r (fm)$ (the distance between both nucleons). -
calculate_properties: finally we calculate the probability of the system being in the
$u_d$ D-wave state (after normalizing the total wavefunction). We also reach the Q electric quadrupole moment.
Let us see what procedure we must follow to develop the code.
- We introduce the model we want for the potential and the desired initial and boundary conditions. In our case, we used the fact that
$u_s(r) \sim r$ for small r values and$u_d(r) \sim r^3$ . This way, we had a method to implement the initial conditions for two independent solutions which would obey these properties in the initial point. For the final point, we used the fact that both eigenfunctions follow the form$u \sim e^{-kr}$ for large r values. So we have an idea of how the wavefunctions will behave in the final value. - We implement the equations we must solve. As we are working with second order differential equations, it is convenient to convert each on a system on two first order ODEs. Later, this systems will be ready to be solved by many integration methods which work for first order differential equations.
- Once we have inserted the system, we must follow a specific procedure to get all conditions obeyed. Instead of implementing the shooting method directly from the starting point to the final one, we use the shooting method to a fitting point. This means that we will make the solutions be continuous in a midpoint R rather than trying to make the functions obey the final boundary conditions (small changes make the functions diverge due to the exponential behaviour for large r values). So, we solve the 1st order ODEs using an integration method implemented by SciPy and get the functions from the initial (ending) point to the midpoint.
- Having this situation in mind, since we are facing two coupled second order equations, we will have two linearly independent solutions that are regular in the origin and obey our initial conditions:
$u_{s_A}, u_{d_A}$ and$u_{s_B}, u_{d_B}$ . The same happesn in the final point with the boundary conditions. Consequently, we have:
So we must find the values for A, B, C and D which obey:
This is precisely defined in find_results.py, where we describe 4 root functions in error_ABCD to find these constants. Later, we make use of the still non-defined continuity condition:
- We get the roots for the equations by applying root-finding methods from SciPy both for a system of equations (to get A, B, C and D constants) and for a single function (to get the energy E) and save the results.
- We use the obtained values in plotWFs to plot both wavefunctions using get_graph to solve again the system through an integration method, but this time with the saved appropriate values.
- With the saved results we can now normalize the final wavefunction
$\psi_{final}= \alpha u_s(r) + \beta u_d(r)$ . After that, we calculate the Q electric quadrupole moment and the probability of being in the$L=2$ D-state.
In order to get the whole procedure done, the libraries required will be:
If we want to get the final results of this code straightforwardly, we just must be sure that all the libraries needed have been correctly installed. After that, we can download the whole code and run module plotWFs to reach the graph of the wavefunctions. We should get one picture like
Also, to obtain the values for Q and