Quantum simulations of nuclei and nuclear pasta with the multiresolution adaptive numerical environment for scientific simulations

Background: Neutron star and supernova matter at densities just below the nuclear matter saturation density is expected to form a lattice of exotic shapes. These so-called nuclear pasta phases are caused by Coulomb frustration. Their elastic and transport properties are believed to play an important role for thermal and magnetic field evolution, rotation, and oscillation of neutron stars. Furthermore, they can impact neutrino opacities in core-collapse supernovae.

Purpose: In this work, we present proof-of-principle three-dimensional (3D) Skyrme Hartree-Fock (SHF) simulations of nuclear pasta with the Multi-resolution ADaptive Numerical Environment for Scientific Simulations (MADNESS).

Methods: We perform benchmark studies of ${}^{16}\mathrm{O}, {}^{208}\mathrm{Pb}$, and ${}^{238}\mathrm{U}$ nuclear ground states and calculate binding energies via 3D SHF simulations. Results are compared with experimentally measured binding energies as well as with theoretically predicted values from an established SHF code. The nuclear pasta simulation is initialized in the so-called waffle geometry as obtained by the Indiana University Molecular Dynamics (IUMD) code. The size of the unit cell is 24 fm with an average density of about $\rho =0.05{\mathrm{fm}}^{-3}$, proton fraction of $Y_p=0.3$, and temperature of $T=0$ MeV.

Results: Our calculations reproduce the binding energies and shapes of light and heavy nuclei with different geometries. For the pasta simulation, we find that the final geometry is very similar to the initial waffle state. We compare calculations with and without spin-orbit forces. We find that while subtle differences are present, the pasta phase remains in the waffle geometry.

Conclusions: Within the MADNESS framework, we can successfully perform calculations of inhomogeneous nuclear matter. By using pasta configurations from IUMD it is possible to explore different geometries and test the impact of self-consistent calculations on the latter.

Figure 1

(a) Number density profiles of the ${}^{16}\mathrm{O}$ nucleus ground state, as obtained with our MADNESS code, taken along the $x, y$, and $z$ axes. The profiles show that the resulting nucleus is spherically symmetric. (b) Number density, kinetic density $\tau$, and Laplacian of the number density $\mathrm{\Delta }\rho$ along the $x$ axis.

Figure 2

Evolution of the maximum error and binding energy with iterations of the ${}^{16}\mathrm{O}$ calculation.

Figure 3

(a) Number density profiles of the ${}^{208}\mathrm{Pb}$ along the $x, y$, and $z$ axes in the initial state (subscript $i$) and at convergence (subscript $f$). Panel (b) shows the number density, kinetic density $\tau$ and Laplacian of the number density $\mathrm{\Delta }\rho$ at iteration 280.

Figure 4

Zoom of Fig. 3 showing the $\rho , \tau$, and $\mathrm{\Delta }\rho$ profiles along $x$ for the ${}^{208}\mathrm{Pb}$ simulation (I) with smoothing and $\epsilon ={10}^{-5}$ (marked by a $☆$) and $\epsilon ={10}^{-7}$ without smoothing. While the oscillation patterns of the three quantities are similar with and without smoothing, small differences are visible in the amplitudes.

Figure 5

Evolution of the maximum error and binding energy with iterations of the ${}^{208}\mathrm{Pb}$ calculation.

Figure 6

(a) Number density profiles of the ${}^{238}\mathrm{U}$ along the $x, y$, and $z$ axis at convergence (subscript $f$) and the beginning of the iterations (subscript $i$). Panel (b) shows the number density, kinetic density $\tau$, and Laplacian of the number density $\mathrm{\Delta }\rho$ at iteration 2800.

Figure 7

Zoom of Fig. 6 showing the $\rho , \tau$, and $\mathrm{\Delta }\rho$ profiles along $x$ for the ${}^{238}\mathrm{U}$ simulation (I) with smoothing and $\epsilon ={10}^{-5}$ (marked by a $☆$) and $\epsilon ={10}^{-7}$ without smoothing.

Figure 8

Evolution of the maximum wave function change $\delta \psi$ and binding energy $E_{\mathrm{bind}}$ with iterations for the ${}^{238}\mathrm{U}$ nucleus.

Figure 9

Evolution of the maximum wave function change $\delta \psi$ and binding energy $E_{\mathrm{bind}}$ with iterations for the ${}^{238}\mathrm{U}$ nucleus for simulation setup (II), $\epsilon ={10}^{-7}$, and no Gaussian smoothing.

Figure 10

Density isosurface of the ${}^{238}\mathrm{U}$ nucleus for $\rho =0.08{\mathrm{fm}}^3$ for simulation (II) with $\epsilon ={10}^{-7}$ and no Gaussian smoothing. Panel (a) shows the orientation of the nucleus in the simulation space at iteration 4000. Panel (b) shows the nucleus at iteration $19400$. The black lines mark the symmetry axis.

Figure 11

Positions of neutrons (blue) and protons (gray) in the converged IUMD pasta simulation. Panels show different orientations of the simulations space. Sphere sizes are for visualization only.

Figure 12

Isosurfaces of initial total baryon number density $\rho \left(\vec{r}\right)$ of M-SHF. Nucleon wave functions are Gaussians folded around coordinates from IUMD shown in Fig. 11. Panels correspond to orientations of the simulations space as in Fig. 11. See text for details.

Figure 13

Evolution of the maximum wave-function change $\delta \psi$ and binding energy $E_{\mathrm{bind}}$ with iterations for the MD M-SHF pasta simulation.

Figure 14

Isosurfaces of the total baryon density $\rho \left(\vec{r}\right)$ of the converged M-SHF simulation at iteration 8000. Panels are as in Fig. 12.

Figure 15

(a) Number density profiles of the nuclear pasta configuration along the $x, y$, and $z$ axes at iteration 4600 and in the final, i.e., converged, state at iteration 8000. (b) Comparison of the $x$ profiles of the total number density $\rho$, kinetic density $\tau$, and Laplacian of the number density $\mathrm{\Delta }\rho$ in the final, i.e., converged, state of the nuclear pasta and at iteration 4600.

Figure 16

Evolution of the maximum wave-function change $\delta \psi$ and binding energy $E_{\mathrm{bind}}$ with iterations for the MD M-SHF pasta simulation.

Figure 17

Isosurfaces of $\rho \left(\vec{r}\right)$ as in Fig. 14. The orientation corresponds to the top plate of the converged waffle phase. Panel (a) shows the result of the M-SHF simulation without spin-orbit and (b) with spin-orbit interactions. Panel (c) is the converged Sky3D simulation with spin-orbit.

Figure 18

Isosurfaces of $\rho \left(\vec{r}\right)$ as in Fig. 14. The orientation corresponds to the bottom plate of the converged waffle phase. Panel (a) shows the result of the M-SHF simulation without spin-orbit and (b) with spin-orbit interactions. Panel (c) is the converged Sky3D simulation with spin-orbit.