System energy

Introduction

The energy of the nuclear system consists of a few terms. Those terms will be presented below. Constants before these terms are dependent on decided functionals. More info about functionals one can find here.

This page is organized as follows. We divide it into sections. One section corresponds to one energy term. First, we describe it we provide a mathematical formula for that term. At the end of each energy section, we will present a pseudo-code that LISE is using during simulations.

Technically LISE is calculating each term for each point of the box. All of the terms are directly dependent on the densities. More info about densities here.

Notation

Two sumation will be often use in this page. Here we present the meaning of them

\begin{align} \rho_{0} &= \int_V{\left(\rho_n + \rho_p\right)}\\ \rho_{1} &= \int_V{\left(\rho_n - \rho_p\right)}, \end{align}

where \rho_n and \rho_p, corresponds to density of neutrons and protons respectively.

Inside the code sections, we will call them rho_0 and rho_1 respectively, but yet inside the code it is a one dimensional array that contains a sum or a difference of the densities.

This notation will apply also for other densities such as kinetic \tau_{0,1}, anomalous \nu_{0,1} and currents j_{0,1}

Kinetic energy

Kinetic energy is a term devoted to the motion of the nuclei.

\begin{align} E_{kin} &= \tfrac{\hbar c}{m_p + m_n} \int_V \tau_n + \tau_p \\ &= \tfrac{\hbar c}{m_p + m_n} \int_V \tau_0, \end{align}

where mp and mn are the mass of protons and neutrons respectively, τ is kinetic density.

e_kin = 0.0
for ixyz in nxyz:
    e_kin += hbarc2m * (density.tau_n[ixyz] + density.tau_p[ixyz]) * dxyz = 
                hbarc2m * tau_0[ixyz] * dxyz

Normal density energy

todo: fill explanation

\begin{align} E_{\rho} &= \int_V\ C^\rho_0 \left(\rho_n + \rho_p\right) + \int_V\ C^\rho_1 \left(\rho_n - \rho_p\right) \\ &= \int_V\ (C^\rho_0 \rho_{0} + C^\rho_1 \rho_1), \end{align}

where

\begin{align} C^\rho_0 &= \tfrac{3}{8}t_0\\ C^\rho_1 &= -\tfrac{1}{4}\left(x_0+\tfrac{1}{2}\right) \end{align} 

e_rho = 0.0
for ixyz in nxyz:
    e_rho += (c_rho_0 * (density.rho_n[ixyz] + density.rho_p[ixyz]) +
                c_rho_1 * (density.rho_n[ixyz] - density.rho_p[ixyz])) 
            * dxyz
            = (c_rho_0 * rho_0[ixyz] + c_rho_1 * rho_1[ixyz]) * dxyz

Normal density kinetic energy todo: naming

todo: fill explanation

\begin{align} E_{\rho\tau} &= C^{\rho\tau}_0 \int_V\ (\rho_n + \rho_p)(\tau_n + \tau_p) + C^{\rho\tau}_1 \int_V\ (\rho_n + \rho_p)(\tau_n - \tau_p) \\ &= C^{\rho\tau}_0 \int_V\ \rho_0\tau_0 + C^{\rho\tau}_1 \int_V\ \rho_1\tau_1, \end{align}

where

\begin{align} C^{\rho\tau}_0 &= \tfrac{3}{16}t_1 + \tfrac{1}{4}t_2\left(x_2 + \tfrac{5}{4}\right)\\ C^{\rho\tau}_1 &= -\tfrac{1}{8} \left[ t_1 \left( x_1 + \tfrac{1}{2} \right) - t_2 \left( x_2 + \tfrac{1}{2} \right) \right] \end{align}

and \tau is described here

e_rhotau = 0.0
for ixyz in nxyz:
    e_rhotau += (c_rhotau_0 * 
                    ((density.rho_n[ixyz] + density.rho_p[ixyz]) * (density.tau_n[ixyz] + density.tau_p[ixyz])) +
                    c_rhotau_1 * 
                    ((density.rho_n[ixyz] - density.rho_p[ixyz]) * (density.tau_n[ixyz] - density.tau_p[ixyz])))
                * dxyz
                = (c_rhotau_0 * (density.rho_0[ixyz] * density.tau_0[ixyz]) +
                    c_rhotau_1 * (density.rho_1[ixyz] * density.tau_1[ixyz]))
                * dxyz

Gamma energy

todo: fill explanation

\begin{align} E_\gamma &= C^\gamma_0 \int_V\ (\rho_n+\rho_p)^{\gamma + 2} + C^{\gamma}_1 \int_V\ (\rho_n +\rho_p)^{\gamma} (\rho_n - \rho_p)^{2} \\ &= C^\gamma_0 \int_V\ \rho_0^{\gamma + 2} + C^{\gamma}_1 \int_V\ \rho_0^{\gamma} \rho_1^{2}, \end{align}

where

\begin{align} C^\gamma_0 &= \tfrac{t_3}{16}\\ C^\gamma_1 &= -\tfrac{1}{24}\left(x_3+\tfrac{1}{2}\right)t_3 \end{align} 

e_gamma = 0.0
for ixyz in nxyz:
    e_gamma += (c_gamma_0 * (density.rho_n[ixyz] + density.rho_p)**(gamma+2) + 
                c_gamma_1 * ((density.rho_n[ixyz] + density.rho_p)[ixyz]**gamma * (density.rho_n[ixyz] - density.rho_p)[ixyz]**2))) * dxyz
                = (c_gamma_0 * rho_0[ixyz]**(gamma+2) + c_gamma_1 * rho_0[ixyz]**gamma * rho_1[ixyz]**2) * dxyz

Spin-orbit energy

todo: fill explanation

\begin{align} E_{so} &= C^{\nabla J}_0 \int_V\ (\rho_n + \rho_p)(\nabla J_n + \nabla J_p) + C^{\nabla J}_1 \int_V\ (\rho_n - \rho_p) (\nabla J_n - \nabla J_p) \\ &= C^{\nabla J}_0 \int_V\ \rho_0 \nabla J_0 + C^{\nabla J}_1 \int_V\ \rho_1 \nabla J_1, \end{align}

where

\begin{align} C^{\nabla J}_0 &= -\tfrac{3}{4} W_0\\ C^{\nabla J}_1 &= -\tfrac{1}{4} W_0 \end{align} 

e_so = 0.0
for ixyz in nxyz:
    e_so += (c_divjj_0 * ((density.rho_n[ixyz] + density.rho_p[ixyz]) * (density.divjj_n[ixyz] + density.divjj_p[ixyz])) +
                c_divjj_1 * ((density.rho_n[ixyz] - density.rho_p[ixyz]) * (density.divjj_n[ixyz] - density.divjj_p[ixyz])))
                * dxyz
            = (c_divjj_0 * (density.rho_0[ixyz] * density.divjj_0[ixyz]) +
                c_divjj_1 * (density.rho_1[ixyz] * density.divjj_1[ixyz]))
e_so *= dxyz

Density laplacian energy

todo: fill explanation

\begin{align} E_{\nabla\rho} &= C^\nabla_0 \int_V\ (\rho_n+\rho_p) \nabla (\rho_n+\rho_p) + C^\nabla_1 \int_V\ (\rho_n-\rho_p) \nabla (\rho_n-\rho_p) \\ &= C^\nabla_0 \int_V\ \rho_0 \nabla \rho_0 + C^\nabla_1 \int_V\ \rho_1 \nabla \rho_1, \end{align}

where

\begin{align} C^\nabla_0 &= -\tfrac{1}{16}\left[\tfrac{9}{4}t_1 - t_2\left(x_2+\tfrac{5}{4}\right)\right]\\ C^\nabla_1 &= -\tfrac{1}{32}\left[3t_1 - \left(x_1+\tfrac{1}{2}\right) - t_2\left(x_2+\tfrac{1}{2}\right)\right] \end{align} 

e_laprho = 0.0
for ixyz in nxyz:
    e_laprho += (c_laprho_0 * ((density.rho_n[ixyz] + density.rho_p[ixyz]) * (density.lap_rho_n[ixyz] + density.lap_rho_p[ixyz]))
                + c_laprho_1 * ((density.rho_n[ixyz] - density.rho_n[ixyz]) * (density.lap_rho_n[ixyz] - density.lap_rho_p[ixyz])))
                * dxyz
                = (c_laprho_0 * (density.rho_n[ixyz] * density.lap_rho_0[ixyz])
                + c_laprho_1 * (density.rho_1[ixyz] * density.lap_rho_1[ixyz]))
e_laprho *= dxyz

Pairing energy

todo: fill explanation

E_{pairing, q} = - \int_V\ \operatorname{Re}(\Delta_q)\ \nu^*_q

where q stands for protons or neutrons

e_pair= 0.0
for ixyz in nxyz:
    e_pair_q -= creal(potentials.delta_q[ixyz]) * conj(density.nu_q[ixyz]) * dxyz

where pairing energy for protons and neutrons are calculatede separately and are indicated by the _q in the code.

Flow energy

todo: fill explanation

E_{flow,q} = \int_V\ (j_{x,q}^2 + j_{y,q}^2 + j_{z,q}^2) / \rho_q

Flow energy is calculated only for the lattice points in which \rho is greater than a threshold value. In the current version of the code it is set to be 10^{-7}

e_flow_q = 0.0
for ixyz in nxyz:
    if density.rho_q[ixyz] > 1.e-7:
        e_flow_q += (pow(density.current_x_q[ixyz],2) + pow(density.current_y_q[ixyz],2) + pow(density.current_z_p[ixyz],2))
                        / density.rho_q[ixyz]
e_flow_q *= dxyz

Current energy

todo: fill explanation

\begin{align} E_j = &C^j_0 \int_V\ (j_{x,n} + j_{x,p})^2 + (j_{y,n} + j_{y,p})^2 + (j_{z,n} + j_{z,p})^2 \\ + &C^j_1 \int_V\ (j_{x,n} - j_{x,p})^2 + (j_{y,n} - j_{y,p})^2 + (j_{z,n} - j_{z,p})^2 \\ + &C^{\nabla j}_0 \int_V\ (s_{x,n} + s_{x,p})(sc_{x,n} + sc_{x,p}) + (s_{y,n} + s_{y,p})(sc_{y,n} + sc_{y,p}) + (s_{z,n} + s_{z,p})(sc_{z,n} + sc_{z,p}) \\ + &C^{\nabla j}_0 \int_V\ (s_{x,n} - s_{x,p})(sc_{x,n} - sc_{x,p}) + (s_{y,n} - s_{y,p})(sc_{y,n} - sc_{y,p}) + (s_{z,n} - s_{z,p})(sc_{z,n} - sc_{z,p}) \end{align}

where s stands for spin density and sc for spin current. Moreover constants can be described as follows

\begin{align} C^j_0 &= - C^{\rho\tau}_0 = -\tfrac{3}{16}t_1 + \tfrac{1}{4}t_2\left(x_2 + \tfrac{5}{4}\right)\\ C^j_1 &= - C^{\rho\tau}_1 = \tfrac{1}{8} \left[ t_1 \left( x_1 + \tfrac{1}{2} \right) - t_2 \left( x_2 + \tfrac{1}{2} \right) \right] \\ C^{\nabla j}_0 &= C^{\nabla J}_0 = -\tfrac{3}{4} W_0\\ C^{\nabla j}_1 &= C^{\nabla J}_1 = -\tfrac{1}{4} W_0 \end{align} 

e_j = 0.0
for ixyz in nxyz:
    e_j += c_j_0 * (pow(density.j_n_x[ixyz] + density.j_p_x[ixyz], 2) + 
                    pow(density.j_n_y[ixyz] + density.j_p_y[ixyz], 2) + 
                    pow(density.j_n_z[ixyz] + density.j_p_z[ixyz], 2)) +
           c_j_1 * (pow(density.j_n_x[ixyz] - density.j_p_x[ixyz], 2) + 
                    pow(density.j_n_y[ixyz] - density.j_p_y[ixyz], 2) + 
                    pow(density.j_n_z[ixyz] - density.j_p_z[ixyz], 2)) +
           c_divj_0 * ((density.s_n_x[ixyz] + density.s_p_x[ixyz])(density.sc_n_x[ixyz] + density.sc_p_x[ixyz]) + 
                       (density.s_n_y[ixyz] + density.s_p_y[ixyz])(density.sc_n_y[ixyz] + density.sc_p_y[ixyz]) + 
                       (density.s_n_z[ixyz] + density.s_p_z[ixyz])(density.sc_n_z[ixyz] + density.sc_p_z[ixyz])) +
           c_divj_1 * ((density.s_n_x[ixyz] - density.s_p_x[ixyz])(density.sc_n_x[ixyz] - density.sc_p_x[ixyz]) + 
                       (density.s_n_y[ixyz] - density.s_p_y[ixyz])(density.sc_n_y[ixyz] - density.sc_p_y[ixyz]) + 
                       (density.s_n_z[ixyz] - density.s_p_z[ixyz])(density.sc_n_z[ixyz] - density.sc_p_z[ixyz]))
         = c_j_0 * (pow(density.j_x_0[ixyz], 2) + pow(density.j_y_0[ixyz], 2) + pow(density.j_z_0[ixyz], 2)) +
           c_j_1 * (pow(density.j_x_1[ixyz], 2) + pow(density.j_y_1[ixyz], 2) + pow(density.j_z_1[ixyz], 2)) +
           c_divj_0 * ((density.s_x_0[ixyz]*density.sc_x_0[ixyz]) + (density.s_y_0[ixyz]*density.sc_y_0[ixyz])) + (density.s_z_0[ixyz]*density.sc_z_0[ixyz]) +
           c_divj_1 * ((density.s_x_1[ixyz]*density.sc_x_1[ixyz]) + (density.s_y_1[ixyz]*density.sc_y_1[ixyz])) + (density.s_z_1[ixyz]*density.sc_z_1[ixyz])
e_j *= dxyz