Exact restoration of Galilei invariance in density functional calculations with quantum Monte Carlo

Within the HF approach, the nuclear wave function is assumed to be a Slater determinant $|\Phi\rangle$ formed by a set of single-particle wave functions $\chi_{\alpha}$, where $\alpha$ stands for the spherical quantum numbers:

In the above equation, $\mathcal{A}$ is the antisymmetrization operator and $|X\rangle=|x_{1},\ldots,x_{A}\rangle$, and $x_{i}=\{{\bf r}_{i},\sigma_{i},\tau_{i}\}$ are generalized coordinates representing position, spin, and isospin of the $i$th nucleon.

The nuclear mean-field is determined by finding the Slater determinant $|\Phi\rangle$ that minimizes the energy

which is a functional of the single-particle density

written in terms of the creation and annihilation operators $a^{\dagger}_{\alpha},a_{\alpha}$.

Following [11], the solution of such a minimization defines the HF single-particle basis $\phi_{a}$ and the HF average potential

where $h$ depends on the density and obeys the relation

The HF solution provides a single-particle basis in which both $h$ and $\rho$ are diagonal.

In order to determine the mean-field wave functions, they are expanded in a spherical Woods-Saxon basis $\chi_{\alpha}$

so that

The self-consistent HF equations can then be written as an eigenvalue problem (see [11], Eq. (5.38))

where the matrix elements $h_{\alpha\alpha^{\prime}}$ are given by

and $t$ and $\bar{v}$ denote, respectively, the matrix elements of the kinetic energy and the antisymmetrized matrix elements of the two-body interaction. These equations must be solved by iteration, assuming an initial choice for the transformation coefficients $D_{\alpha}$. In the case of density-dependent interactions, the matrix element $\bar{v}$ must be recalculated at each iteration, and furthermore the rearrangement term $\partial\bar{v}/\partial\rho$ must be taken into account. The Woods-Saxon basis used to compute the matrix elements of the two-body interaction ${\bar{v}}$ is obtained by solving the Schrödinger equation for a Woods-Saxon potential of standard form [34] in a box of $20\,\rm fm$, using an energy cutoff of $400\,\rm MeV$ [35]. The use of a Woods-Saxon basis is particularly convenient in the case of weakly-bound systems, but it is equivalent to the use of a harmonic oscillator basis for the case of the well bound nuclei considered here.

According to the fitting procedure adopted to establish the parameters of the Gogny interaction [19], HF calculations are performed correcting for the violation of translation invariance by subtracting the kinetic CM contribution $T_{\rm CM}$ from the Hamiltonian, where

and the mass of the nucleus is given by $M=A\,m_{N}$. The first term, a one-body operator, can be dealt with by simply rescaling the total kinetic energy:

On the other hand, the last term of equation (16) is a two-body operator, and it has to be treated on the same footing as the two-body components of the potential of equation (2). We remark that in the fitting protocol of other effective interactions, only the one-body part of the kinetic energy correction is taken into account.

In spite of this subtraction procedure, the resulting wave functions still violate the basic translational symmetry. While symmetry breaking in nuclear physics has been the subject of an immense literature, in the case of translational invariance (see [36, 37, 38, 39] for recent studies) there are surprisingly few quantitative estimates of the errors due to the violation of this symmetry in HF calculations. This is mostly due to the fact that mean-field theory is mainly applied to medium-heavy nuclei, for which this violation, which is expected to scale as $1/A$, is small. A remarkable exception is represented by a series of papers by K.W. Schmid, who calculated the binding energies and radii of few closed-shell nuclei, including \isotope[4]He, \isotope[16]O, and \isotope[40]Ca, fully restoring the Galilean invariance in the case of density-independent interactions by making use of the analytic properties of harmonic oscillator configurations [15]. It should be noted that the density-dependent part of the Gogny interaction $v^{\rm ddp}_{ij}$, that depends on the density at the center of mass of the nucleon pair ${\bf R}_{ij}$, explicitly violates Galilean invariance. Schmid modified this dependence using the density $\rho^{\alpha}({\bf R}_{ij}-{\bf R}_{\rm CM})$ instead of $\rho^{\alpha}({\bf R}_{ij})$ (see also the discussion about internal density in [40]). This procedure restores Galilean invariance, but it introduces a $A$-body force that makes HF calculations much more difficult to perform. On the other hand, this modification is suitable for QMC calculations, as discussed in the next section.

In the VMC method [21], provided a trial wave function $\Psi_{T}$, the expectation value of the Hamiltonian is given by

where $E_{0}$ is the energy of the true ground state with the same quantum numbers as $\Psi_{T}$, and the rightmost equality is valid only if the wave function is the exact ground-state wave function $\Psi_{0}$. The energy expectation value of equation (18) typically depends on the quality of the employed wave function. In the VMC method, one minimizes $E_{V}$ with respect to changes in the variational parameters, in order to obtain $\Psi_{T}$ as close as possible to $\Psi_{0}$.

In this work we employ a trial wave function of the form

where $\langle X|\Phi\rangle$ is the same Slater determinant of equation (7). In our VMC calculations, the single-particle states are taken to be

where $R_{nl}(r)$ is the radial function, $Y_{ll_{z}}$ is the spherical harmonic, and $Y_{ss_{z}}(\sigma)$ and $Y_{tt_{z}}(\tau)$ are the complex spinors describing the spin and isospin of the single-particle state.

As described in [23], the radial functions $R_{nl}(r)$, are expressed as a sum of cubic splines, characterized by a smooth first derivative and a continuous second derivative. The large number of variational parameters involved in the construction of such radial components allows enough flexibility to obtain very accurate trial wave functions, with optimal variational energies very close to those calculated by performing the imaginary-time propagation (see [21] and references therein for more details).

Note that, using the wave function of equation (19) and inserting a completeness over the generalized coordinate $X$, the expectation value of a generic operator $\mathcal{O}$ can be expressed as

where $P(X)=|\Psi_{T}(X)|^{2}$ can be interpreted as a probability distribution of points $\{X\}$ in a multidimensional space, and

is referred to as the local expectation value. Equation (21) actually corresponds to a multidimensional integral that can be calculated using Monte Carlo sampling. According to the Metropolis algorithm [41], a number of configurations $X_{i}$ are sampled from the probability distribution $P(X)$, and the expectation value of the operator $\mathcal{O}$ is calculated as

where $\mathcal{M}$ is the number of sampled configurations. Details on the sampling procedure and Monte Carlo statistical errors can be found, e.g. in [42].

Note that, in order to remove CM contributions, the nuclear wave function and the observables are calculated in the CM rest frame (COMF in the following) instead than in a fixed rest frame (FRF in the following), subtracting the CM position from all the spatial coordinates:

By introducing the “intrinsic” coordinate ${\bf r}^{\rm int}_{i}={\bf r}_{i}-{\bf R}_{\rm CM}$, the above procedure is equivalent to replacing the radial function and the spherical harmonic of equation (20):

Each single-particle orbital depends now upon the coordinates of all the nucleons, in a way that resembles backflow correlations [43]. Hence, although the trial wave function is constituted by a single Slater determinant, subtracting the CM coordinate brings about some degree of correlations among the nucleons, whose importance decreases with $A$. As a consequence, the determination of the radial wave functions that minimize the expectation value of the energy becomes a numerically challenging problem, because it involves the calculation of non-factorizable integrals in $3A$ dimensions. However, this problem lends itself to a solution based on QMC techniques, following the iterative procedure outlined in the next section.

In order to optimize the radial components of the trial wave function, in this work we adopt the linear method (LM) [22], that was applied for the first time in a nuclear quantum Monte Carlo calculation in [23]. Let us first define the normalized trial variational state

as a function of the $\mathcal{N}_{p}$ variational parameters ${\bf p}=\{p_{1},\dots,p_{\mathcal{N}_{p}}\}$.

Within the LM, at each optimization step one performs a first-order expansion around the current set of variational parameters ${\bf p}^{0}$

where $|\bar{\Psi}_{T}^{0}({\bf p}^{0})\rangle\equiv|\Psi_{T}({\bf p}^{0})\rangle$, and for $i>0$

In the latter equation the first derivative with respect to the $i$th parameter is given by

while the overlap matrix is defined as

By using the normalization freedom we can impose $\langle\bar{\Psi}_{T}^{0}({\bf p}^{0})|\bar{\Psi}_{T}^{0}({\bf p}^{0})\rangle=1$, so that the derivatives of $|\bar{\Psi}_{T}({\bf p})\rangle$ are orthogonal to $|\bar{\Psi}_{T}^{0}({\bf p}^{0})\rangle$

The eigenvalue equation for the Hamiltonian in the basis formed by the $(\mathcal{N}_{p}+1)$-dimensional basis $\Big{\{}|\bar{\Psi}_{T}^{0}({\bf p}^{0})\rangle,\dots,|\bar{\Psi}^{\mathcal{N}_{p}}_{T}({\bf p}^{0})\rangle\Big{\}}$ reads

Multiplying the latter equation by $\langle\bar{\Psi}^{i}_{T}({\bf p}^{0})|$ yields to the generalized eigenvalue equation

where the Hamiltonian and overlap matrix elements are

The linear method consists of solving equation (33) for the lowest eigenvalue and associated eigenvector $\Delta\bar{{\bf p}}$. It has to be noted that equation (33) can alternatively be obtained by minimizing the energy expectation value on the linear wave function

with respect to changes in the variational parameters [22].

The expressions of the above matrix elements for real variational parameters can be found in [44]. Here we report their expressions for the complex case [45], needed in nuclear QMC calculations. Inserting a completeness over the generalized coordinate $X$, the Hamiltonian and overlap matrix elements read (for brevity the dependence on ${\bf p}^{0}$ is understood)

By making explicit the definition of $|\bar{\Psi}^{i}_{T}\rangle$ one obtains

where the matrix element $\langle X|H|\Psi^{i}_{T}\rangle$ can be expressed in terms of the the local energy $E_{L}(X)=\langle X|H|\Psi_{T}\rangle/\langle X|\Psi_{T}\rangle$ (see equation (22)) and its derivative $E_{L}^{i}(X)$:

The derivatives of the wave function and of the local energy are numerically evaluated using the three-point stencil rule by shifting the $i-$parameter by a small quantity $\delta_{p}$: $p_{i}\to p_{i}\pm\delta_{p}$.

The evaluation of $S_{0i}$ requires inserting an additional completeness relation

where the rightmost term is intended as in equation (21).

Collecting the above results, the Hamiltonian and overlap matrix elements can then be expressed as an average over finite Monte Carlo samples, as in equation (23):

The hermiticity of the matrix $\bar{H}$ is recovered only in the limit of an infinite Monte Carlo sample. However, even for a finite sample, we refrain from symmetrizing $H_{ij}$, as employing non Hermitian estimators leads to a stronger zero-variance principle [46], and to smaller statistical errors on a finite sample than using its symmetrized analog [44]. It has been shown that the Monte Carlo statistical uncertainty is largely reduced by the fact that the above matrix elements are written in terms of covariances [44]. Nevertheless, for a finite sample size, the matrix $\bar{H}$ can be ill-conditioned, preventing a stable solution of the eigenvalue problem. In order to stabilize the algorithm, we add a small positive constant $\varepsilon$ to the diagonal matrix elements of $\bar{H}$ except for the first one, $\bar{H}_{ij}\to\bar{H}_{ij}+\varepsilon\,(1-\delta_{i0})\,\delta_{ij}$. This reduces the length of $\Delta\bar{{\bf p}}$ and it rotates it toward the steepest-descent direction.

Strong nonlinearities in the variational parameters potentially make $|\bar{\Psi}_{T}^{\text{lin}}({\bf p})\rangle$ significantly different from $|\bar{\Psi}_{T}({\bf p}^{0}+\Delta{\bf p})\rangle$. To alleviate this problem we employ the heuristic procedure of [23]. For a given value of $\varepsilon$, equation (33) is solved. If the linear variation of the trial state for ${\bf p}={\bf p}^{0}+\Delta{{\bf p}}$ is small,

a short correlated run is performed in which the energy expectation value is estimated along the full variation of the trial state for a set of possible values of $\varepsilon$ (typically $50$ values are considered). The optimal $\varepsilon$ is the one corresponding to the lowest eigenvalue, provided that

Note that, in latter expression, at variance to equation (41), the full trial state instead of its linearized approximation appears in the numerator. This additional constraint suppresses the potential instabilities caused by the nonlinear dependence of the trial state on the variational parameters. When using the “standard” version of the LM, there were instances in which, despite the variation of the linear trial state being well below the threshold of equation (41), the full trial state fluctuated significantly more, preventing the convergence of the minimization algorithm. We found that choosing $\delta=0.2$ guarantees an ideal compromise between the convergence-rate and the stability of the algorithm.

When using the LM in conjunction with a density dependent interaction, one has to take into account that the parameters will now become themselves density dependent, as a consequence of the presence of the term $v^{\rm ddp}_{ij}$. This needs to be implemented in the calculation of the derivative of the local energy $E_{L}^{i}(X)$. To this aim, for each iteration of the LM, we first perform a preliminary Monte Carlo run, in which we compute the $2\,\mathcal{N}_{p}+1$ single-nucleon densities corresponding to the current set of variational parameters plus those with $p_{i}\to p_{i}+\delta_{p}$ and $p_{i}\to p_{i}-\delta_{p}$. We have found that this procedure reduces the benefit of using the non-hermitian estimator of $H_{ij}$ of equation (40). Indeed, explicitly symmetrizing $H_{ij}$ by means of $H_{ij}\to H_{ij}+H^{*}_{ji}$ does not significantly alter the convergence of the algorithm.

We first carried out a comparison between HF and VMC approaches by optimizing the single-particle orbitals of \isotope[4]He without introducing correlations in the VMC wave function. In both methods, we applied the same CM corrections, consisting in the subtraction of one- and two-body kinetic terms as in equation (16). Hence, for this particular comparison, in the VMC calculation we refrain from subtracting the CM coordinates as in equation (24). The convergence of the LM, shown in figure 1, is fast. Already after 10 optimization steps, the LM converges to the HF energy.

In table 2 we compare HF and LM results in \isotope[4]He and \isotope[16]O for the total energy $E$, the kinetic energy $T$, the expectation value of the different contributions of the potential of equation (2), and the point-nucleon radius. The latter is defined as

where ${\cal N}$ is the number of protons or neutrons,

is the projector operator onto protons or neutrons, and ${\bf r}_{i}$ is the spatial rescaled coordinate defined in equation (24).

The agreement between HF and LM results is good, for both energies and radii. The reason for the small discrepancies is mainly due to the fact that in our VMC calculations we use the regularized version of equation (6) rather than a pure contact density-dependent interaction. Note that, in \isotope[16]O, the deviation between the experimental energy ($-127.6\,\rm MeV$) and the HF result ($-129.1\,\rm MeV$) is $\approx 2\,\rm MeV$. The HF value is consistent with the result of [49] ($129.6\,\rm MeV$) once taking into account that the spin-orbit term of the Gogny interaction, not included in the present calculations, gives a contribution of about $0.7\,\rm MeV$ to the binding energy of \isotope[16]O (its contribution to the point-nucleon radius is negligible).

In figures 2 and 3 we show the total HF and LM one-body point-nucleon density (protons plus neutrons) in \isotope[4]He and \isotope[16]O, respectively. Similarly to the point-nucleon radius, the one-body point-nucleon density is defined as

where $\mathcal{P}_{N_{i}}$ is the projector operator of equation (44) and $\rho_{N}$ integrates to the number of nucleons. In QMC calculations the above expression is implemented by constructing the histogram of the nucleon radial coordinates $r_{i}$ multiplied by the expectation value of the projector operator, normalized to the number Monte Carlo samples and the volume element centered in $r_{i}$

where $dr$ is the size of the bin. Note that, since the nucleon coordinates are sampled from a probability distribution constructed using the trial wave function, and the latter is also used to calculate the expectation value of the projection operator, the density of equation (45) is always constructed for the actual wave function employed in the Monte Carlo run (correlated or uncorrelated).

The agreement between the HF and LM densities is excellent, proving once again the accuracy of the LM in optimizing the radial components of the wave function.

To gauge CM contamination in HF calculations, analogously to what done in [50], we performed VMC calculations for \isotope[4]He and \isotope[16]O in the COMF, i.e. subtracting the CM coordinate as in equation (24). Furthermore, we modify the density-dependent term of the Gogny potential, which is manifestly not Galilean invariant [15], by changing the definition of $v^{\rm ddp}_{ij}$ in equation (5) to

In this way, Galilean invariance can be completely restored.

As a first step, we used the single-particle orbitals optimized for the comparative study between HF and LM, performing a PAV calculation. Results in the COMF are reported in the first column of tables 3 and 4, to be compared to those obtained in table 2. As one could have expected, the point-nucleon radii are much smaller when the CM motion is removed from the wave function. Correcting for the one- and two-body terms of equation (16) takes care for most of the CM effects in the kinetic energy in the FRF. Note that the expectation values of these two terms exactly cancel in the COMF, as the CM kinetic energy vanishes in this reference frame. Although CM contamination mainly affects the expectation value of $v_{\rm ddp}$, the expectation values of $v_{4}+v_{\rm C}$ do also change from the FRF to the COMF. This might be somewhat surprising, as these terms only depend on the relative coordinate of the nucleon pair.

In order to clarify this point, let us consider the Hamiltonian [40]

where the first term only depends on ${\bf R}_{\rm CM}$, while the second term is a function of the Jacobi coordinates. The latter, defined as

are invariant under the CM subtraction of equation (24). The ground-state wave function corresponding to the Hamiltonian of equation (48) factorizes as

Hence, the expectation value of any intrinsic operator $\mathcal{O}_{\rm int}({\boldsymbol{\xi}}_{1},\dots,{\boldsymbol{\xi}}_{A-1})$ is independent of the CM coordinates

as the factors $\int d{\bf R}_{\rm CM}|\Psi_{\rm CM}|^{2}$ from the numerator and the denominator simplify. However, since the density-dependent Hamiltonian of the Gogny interaction cannot be written as in equation (48), the factorization of equation (50) does not apply and the expectation value of the intrinsic operator $v_{4}+v_{\rm C}$ might change when computed in the FRF or in the COMF.

As a second step, we performed a VAP calculation by optimizing the single-particle orbitals of \isotope[4]He and \isotope[16]O with the LM in the COMF. To compensate for the increased density-dependent term, the optimization procedure broadens the single-particle orbitals. As a consequence, in the results shown in the second column of tables 3 and 4, the point-nucleon radii are larger than those obtained in the FRF, see table table 2. Consistently, the expectation values of the kinetic energy and of the separate potential terms are also smaller (in absolute values). These effects are, as expected, more prominent in \isotope[4]He, although both nuclei are appreciably under-bound. The effect of CM projection on the binding energies of ⁴He and \isotope[16]O can be obtained comparing the second column of tables 3 and 4 with the HF results reported in table 2 and in figure 1. It accounts for $10.2\,\rm MeV$ in ⁴He and $9.6\,\rm MeV$ in \isotope[16]O. These values are in good agreement with those obtained by Schmid (see figure 7 in [15]).

The set of parameters of the Gogny force we have employed so far are meaningful only if the system is described in the FRF. However, it is possible, in principle, to refit such parameters in order to reproduce observables in the COMF. We leave this task to a future work. Here, we limit ourselves to reducing the strength of the density-dependent term of the Gogny interaction by renormalizing the coefficient $\tau_{0}$ of equation (5) ($1390.6$ MeV) so as to reproduce the experimental binding energies in the case of \isotope[4]He and \isotope[16]O. Since CM effects in $v^{\rm ddp}_{ij}$ are larger in lighter nuclei, the quenching is stronger in \isotope[4]He ($\tau_{0}^{*}=1192.4\,\rm MeV$) than in \isotope[16]O ($\tau_{0}^{*}=1357.5\,\rm MeV$), consistently with the findings of [15]. The binding energies and radii of \isotope[4]He and \isotope[16]O obtained employing the renormalized value of $\tau_{0}^{*}$, denoted as VAP^∗, are listed in the third column of tables 3 and 4.

The quenching of the repulsive density-dependent term of equation (2) implies a reduction of the point-nucleon radii, which is reflected in the longitudinal elastic form factors (charge form factors), as shown in figures 4 and 5. The charge form factor is expressed as the ground-state expectation value of the one-body charge operator [51]

where $\tilde{\rho}_{N}(q)$ is the Fourier transform of the one-body point-nucleon density defined in equation (45), and $Q^{2}_{\rm el}={\bf q}^{2}-\omega_{\rm el}^{2}$ is the four-momentum squared, with $\omega_{\rm el}=\sqrt{q^{2}+m_{A}^{2}}-m_{A}$ the energy transfer corresponding to the elastic peak, $m_{A}$ being the mass of the target nucleus. $G_{E}^{N}(Q^{2})$ are the nucleon electric form factors, for which we adopt Kelly’s parametrization [52]. The above expression is derived ignoring small spin-orbit contributions in the one-body current.

In \isotope[4]He, subtracting the CM contributions and using the quenched $\tau_{0}^{*}$ coefficient significantly improves the agreement with experimental data. As discussed in [21, 61, 62], meson-exchange currents are needed to shift the peaks of the longitudinal elastic form factor to lower values of the momentum transfer and achieve agreement with experiment. On the other hand, the fact that experimental data at small $q$ are underestimated, and consequently the proton radius is overpredicted, suggests that the simple quenching of $\tau_{0}$ does not suffice to adequately describe \isotope[4]He.

CM spuriosities have less effect on the longitudinal form factor of \isotope[16]O. Even in this case, however, working in the COMF significantly improves the agreement with experiment in the whole momentum-transfer region. In particular, our VAP^∗ calculations reproduce very well both the first and the second diffraction peaks, despite the Gogny D1S interaction is not expected to be reliable in this kinematical regime.

Finally we remark that the trends observed in figures 4 and 5, and in particular the differences between the momentum dependence of LM and VAP results, are qualitatively similar to those calculated in [17], although quantitative statements are not possible since a simplified interaction was adopted in [17].

We complete our analyses by studying a heavier nucleus, namely \isotope[40]Ca. In this case, in order to save computing time, we did not perform LM calculations and we used instead the single-particle orbitals obtained from HF as a starting point for obtaining the PAV, VAP, and VAP^∗ results. As shown in table 5, the main differences between the HF and PAV binding energies originate from the density-dependent contribution, whose increase can be directly related to the smaller PAV radii. On the other hand, the HF and PAV expectation values of $v_{4}+v_{C}$ are very similar, reflecting the fact that the factorization of equation (51) is approximately fullfilled for ⁴⁰Ca. Consistently with \isotope[4]He and \isotope[16]O, in order to reduce the contribution of $v_{ddp}$, the LM in the COMF broadens the single-particle orbitals, so that the VAP radii are larger than the PAV ones. In order for VAP^∗ results to match the experimental binding energy, a milder quenching with respect to the original value, $\tau_{0}^{*}=1390.0\,\rm MeV$, is needed in the case of \isotope[40]Ca ($\tau_{0}^{*}=1378.1\,\rm MeV$) in comparison to \isotope[4]He ($\tau_{0}^{*}=1192.4\,\rm MeV$) and \isotope[16]O $(\tau_{0}^{*}=1357.5\,\rm MeV)$. As expected, this procedure yields to smaller point radii, which are in better agreement with the experimental values. A better reproduction of experimental data for the longitudinal form factor, displayed in figure 6, is also observed, with the first, second, and third diffraction peaks correctly reproduced.