Nuclear Charge Radii of $^{10,11}$B

TL;DR

This paper reports the first laser spectroscopy measurement of nuclear charge radii for boron isotopes $^{10,11}$B, combining high-precision calculations and measurements to improve accuracy and provide insights into nuclear structure.

Contribution

It introduces a novel laser spectroscopy method combined with ab initio calculations to determine charge radii of boron isotopes with unprecedented precision.

Findings

Measured charge radius difference: -0.49(12) fm^2

Results support a cluster structure in boron nuclei

Benchmark data for nuclear structure models

Abstract

The first determination of the nuclear charge radius by laser spectroscopy for a five-electron system is reported. This is achieved by combining high-accuracy ab initio mass-shift calculations and a high-resolution measurement of the isotope shift in the $2s^2 2p\, ^2\mathrm{P}_{1/2} \rightarrow 2s^2 3s\, ^2\mathrm{S}_{1/2}$ ground state transition in boron atoms. Accuracy is increased by orders of magnitude for the stable isotopes $^{10,11}$B and the results are used to extract their difference in the mean-square charge radius $\langle r^2_\mathrm{c}\rangle^{11} - \langle r^2_\mathrm{c}\rangle^{10} = -0.49\,(12)\,\mathrm{fm}^2$. The result is qualitatively explained by a possible cluster structure of the boron nuclei and quantitatively used as a benchmark for new ab initio nuclear structure calculations using the no-core shell model and Green's function Monte Carlo approaches.

Tables

Table 1. Table 1: Components of the isotope shift δ ​ ν IS = ν ​ ( B 10 ) − ν ​ ( B 11 ) 𝛿 subscript 𝜈 IS 𝜈 superscript B 10 𝜈 superscript B 11 \delta\nu_{\mathrm{IS}}=\nu\left({}^{10}\mathrm{B}\right)-\nu\left({}^{11}\mathrm{B}\right) in the 2 ​ p 2 ​ P 1 / 2 → 3 ​ s 2 ​ S 1 / 2 → 2 superscript 𝑝 2 subscript 𝑃 1 2 3 superscript 𝑠 2 subscript 𝑆 1 2 2p\,^{2}P_{1/2}\rightarrow 3s\,^{2}S_{1/2} transition of the boron atom.

Component	Value (MHz)
$m{\alpha }^2$	$5027.27$±	$0.03$
$m{\alpha }^4$	$-4.78$±	$0.07$
$m{\alpha }^5$	$0.572$±	$0.005$
$m{\alpha }^6$	$-0.058$±	$0.002$
Total	$5023.00$±	$0.08$
$C\left[\mathrm{MHz}/{\mathrm{fm}}^2\right]$	$16.91$±	$0.09$

Table 2. Table 2: Isotope shift δ ​ ν IS = ν ​ ( B 10 ) − ν ​ ( B 11 ) 𝛿 subscript 𝜈 IS 𝜈 superscript B 10 𝜈 superscript B 11 \delta\nu_{\mathrm{IS}}=\nu\left({}^{10}\mathrm{B}\right)-\nu\left({}^{11}\mathrm{B}\right) in the 2 ​ p 2 ​ P 1 / 2 → 3 ​ s 2 ​ S 1 / 2 → 2 superscript 𝑝 2 subscript P 1 2 3 superscript 𝑠 2 subscript S 1 2 2p\,^{2}\mathrm{P}_{1/2}\rightarrow 3s\,^{2}\mathrm{S}_{1/2} transition of atomic boron. The calculated mass shift listed in Tab.  1 is used to extract the field shift that is related to the finite nuclear size.

	Value (MHz)		Reference
Isotope Shift	$5250$±	$360$	Mrozowski (1939)
	$4110$±	$360$	Burke (1955)
	$5220$±	$150$	Johansson et al. (1993)
	$5031.3$±	$2.0$	this work, Exp.
Mass Shift	$5023.00$±	$0.08$	this work, Theory
Field Shift	$8.8$±	$2.0$	extracted

Equations

$$E(\alpha )$$

$$\Psi (\mathbf{r},\mathbf{\sigma })=\hat{\mathcal{A}}\left({\Xi }_{S,M_S}(\vec{\sigma })\sum _{l=1}^K{c}_l{\psi }_l(\vec{r})\right).$$

$${\psi }_l(\vec{r})=$$

$${\vec{\psi }}_l(\vec{r})=$$

$${\nu }^{(n)}(X\to Y)$$

$${\delta }^{(n)}{\nu }_{\mathrm{MS}}$$

$$E_{\mathrm{fs}}^{(4)}$$

$$E_{\mathrm{fs},\log }^{(6)}=-(Z\alpha {)}^2\ln (Z\alpha m⟨r_{\mathrm{c}}^2⟩)E_{\mathrm{fs}}^{(4)}.$$

$$2s^{2}2p{}^2{\mathrm{P}}_{1/2}\underset{\text{249.75nm}}{\overset{{\lambda }_1}{\to }}2s^{2}3s{}^2{\mathrm{S}}_{1/2}\underset{<\text{271nm}}{\overset{{\lambda }_2}{\to }}{\mathrm{B}}^{+}+{\mathrm{e}}^{-}$$

$$\delta {\nu }_{\mathrm{IS}}={\nu }_{\mathrm{SP}}{(}^{10}\mathrm{B})-{\nu }_{\mathrm{SP}}{(}^{11}\mathrm{B})-\Delta {\nu }_{\mathrm{DP}}{(}^{11}\mathrm{B})/2\cdot r$$

Full text

Introduction—The lightest elements play an exceptional role for the advancement of nuclear and atomic physics: Only here theoretical approaches are sufficiently advanced to calculate both electronic and the nuclear structure from first principles. Laser spectroscopy provides unique benchmarks to test and further advance those models of the fundamental structure of nature. For hydrogen-like systems, atomic theory is sufficiently accurate to calculate transition frequencies including quantum electrodynamic (QED) corrections to such a precision that the mean-square nuclear charge radius $\left\langle r_{\mathrm{c}}^{2}\right\rangle$ can be extracted. This has been demonstrated for hydrogen Udem et al. (1997), muonic hydrogen Pohl et al. (2010), and muonic deuterium Pohl et al. (2016). But already for two-electron systems this is so far not feasible even though first progress towards this goal has been reported Yerokhin et al. (2018). However, calculating the mass-dependent isotope shift $\delta\nu_{\mathrm{MS}}^{A,A^{\prime}}$ between two isotopes $A$ and $A^{\prime}$ in an optical transition has the advantage that all mass-independent contributions and their related uncertainties cancel. This allows to isolate the small portion of the isotope shift that is caused by the change in the mean-square nuclear charge radius $\delta\left\langle r_{\mathrm{c}}^{2}\right\rangle$ between isotopes. A few experiments have already utilized this technique to obtain nuclear charge radii of stable and short-lived isotopes of He Wang et al. (2004); Mueller et al. (2007), Li Ewald et al. (2004); Sánchez et al. (2006), and Be+ Nörtershäuser et al. (2009); Krieger et al. (2012), based on corresponding mass-shift calculations in two- and three-electron systems Yan (2001); Yan et al. (2008); Pachucki and Komasa (2003); Puchalski and Pachucki (2008).

Here, we report the first application of this technique to the five-electron system of atomic boron. We present high-precision calculations of the mass shift and the field-shift factor required to extract the difference in mean-square nuclear charge radius $\delta\left\langle r_{\mathrm{c}}^{2}\right\rangle$ between the two stable boron isotopes 10,11B and a measurement of the isotope shift using resonance ionization mass spectrometry on a thermal atomic beam. The results are compared to new ab initio nuclear structure calculations.

Isotope Shift Calculations— The approach to extract $\delta\left\langle r_{\mathrm{c}}^{2}\right\rangle$ between 10B and 11B is similar to our previous work Lu et al. (2013): The isotope shift is composed of the mass shift (MS) and the field shift (FS) $\delta\nu_{\mathrm{IS}}=\delta\nu_{\mathrm{MS}}+\delta\nu_{\mathrm{FS}}$, where the latter $\delta\nu_{\mathrm{FS}}=C\,\delta\left\langle r_{\mathrm{c}}^{2}\right\rangle$ contains the information on the charge radius difference. To obtain $\delta\nu_{\mathrm{MS}}$ with the required high accuracy we calculate the shift due to the finite nuclear mass in powers of the fine structure constant $\alpha$. Namely, the atomic energy levels are considered as a function of $\alpha$, which is expanded in a power series

[TABLE]

These expansion coefficients are calculated including finite nuclear mass effects but neglecting the nuclear spin, which gives rise to a hyperfine splitting but does not shift the energy levels in first order.

The leading term $E^{(2)}$ is the eigenvalue of the non-relativistic Hamiltonian in the center-of-mass system. This Hamiltonian is used to obtain the non-relativistic wave function $\Psi$ using the variational approach, where it is expressed in the form of $K$-term linear combinations of the five-electron basis functions $\psi_{l}({\vec{r}})$

[TABLE]

The operator $\hat{\mathcal{A}}$ ensures the antisymmetry of the total wave function with respect to the exchange of the electrons. The $\Xi_{S,M_{S}}(\vec{\sigma})$ is an $n$-electron spin eigenfunction with the quantum numbers $S$ and $M_{S}$, and $\vec{\sigma}$ and ${\vec{r}}$ are the $n$-electron vectors in spin and coordinate space. The spatial basis functions are the five-electron explicitly correlated Gaussian (ECG) functions of S and P symmetry, respectively:

[TABLE]

with ${\vec{r}}_{a}$ being the coordinate of the $a$-th particle (electrons and the nucleus). The linear parameters $c_{l}$ are obtained by the standard inverse iteration method. The nonlinear parameters $c_{ab}$ are determined variationally for each basis function in an extensive optimization of the non-relativistic energy $E^{(2)}$ with progressively doubled size from $K=1024$ to $K=8192$ terms Puchalski et al. (2015).

The obtained non-relativistic wave functions are then used for further perturbative calculations of relativistic and QED contributions using the respective Hamiltonians. Finite nuclear mass terms of the non-relativistic energy are treated up to $(m/M)^{2}$ and those of higher order in $\alpha$ up to $m/M$. Accordingly, $E^{(4)}=\langle\Psi|H^{(4)}|\Psi\rangle$ is calculated using the relativistic Hamiltonian $H^{(4)}$ and the non-relativistic wave function $\Psi$. The similar calculation for the leading QED correction $E^{(5)}=\langle\Psi|H^{(5)}|\Psi\rangle$ can be found in the supplemental material. At present, a complete numerical evaluation of $E^{(6)}$ correction for a five-electron system is unfeasible, as the full calculation of $E^{(6)}$ has been performed only for one- and two-electron systems Patkóš et al. (2016, 2017). Based on this experience, the $E^{(6)}$ term is estimated using its dominating contribution built of the leading one-electron terms which are proportional to the contact term $\delta^{3}(r_{a})$. No higher order terms are needed here, since the related uncertainty is much smaller than the one from $E^{(4)}$. Following the convention introduced for $E^{(n)}$, formulas for related contributions to the transition energy $\nu$ between the atomic states $X$ and $Y$ and to the mass shift between the isotopes can be written as

[TABLE]

and are summarized in Tab. 1. To determine the field shift factor C we consider the leading $r_{c}$-dependent correction

[TABLE]

and the logarithmic relativistic correction to the wave function at the origin

[TABLE]

Our recommended value for the constant $C$ is obtained as a sum of two components $C=C^{(4)}+C^{(6)}_{\mathrm{log}}$ and an arithmetic average over the isotopes $C=(C_{11}+C_{10})/2$. The remaining term $(C_{11}-C_{10})\,(\langle r_{\mathrm{c}}^{2}\rangle^{A}+\langle r_{\mathrm{c}}^{2}\rangle^{A^{\prime}})/2$ is as small as the neglected relativistic $\alpha^{4}\,m^{3}/M^{2}$ corrections and is therefore also neglected. The result is given in Tab. 1 with an uncertainty estimated as 50% of the logarithmic correction. More details on this calculation can be found in the supplemental material.

Experiment—Resonance ionization spectroscopy was performed in the two-step ionization scheme

[TABLE]

of neutral boron on a thermal atomic beam. The experimental setup is shown in Fig. 1. Stable boron atoms are emitted from a graphite tube filled with amorphous boron powder which is heated to approximately 2500 K. The beam traverses first an electron-impact ionization segment used to monitor the boron production rate in between the laser spectroscopic measurements and then enters the laser interaction region. The apertures along this path limit the angular spread of the atomic beam to 2 mrad.

For laser excitation and ionization two single-mode continuous-wave laser systems are used: A titanium-sapphire (Ti:Sa) laser is generating light output at approximately 1000 nm and is frequency-doubled twice to generate the resonance transition wavelength of $\lambda_{1}=249.75$ nm. The Ti:Sa frequency is scanned stepwise during the experiment, and stabilized against short term fluctuations via an optical reference cavity at each measurement point. A frequency comb is employed to monitor long-term drifts of the reference cavity. The output power of the resonant laser is kept below 1 mW to avoid power broadening and to limit ionization by the resonant laser light. For non-resonant ionization, a commercial cw frequency-doubled Nd:YAG Laser at 532 nm is externally frequency-doubled to $\lambda_{2}=266$ nm providing a minimum of 500 mW power. The resonant and the non-resonant laser beam are overlapped before being intersected perpendicularly with the boron beam.

At resonance, the tunable laser excites boron atoms to the $3s$ state from where they are ionized by absorption of another photon from either the resonant laser or the ionization laser. The generated ions are guided electrostatically into a quadrupole mass spectrometer (QMS) for mass separation and single-ion detection at low background using a channeltron detector. The data acquisition system records the number of detected ions as a function of the frequency of the resonant laser.

After the combined laser beams intersect the atomic beam once, they are back-reflected from a mirror outside the vacuum chamber and overlapped with a precision better than 1 mrad to the incoming laser beams.

To obtain a single spectrum, the laser is scanned and stabilized in steps of 8 to 10 MHz covering the atomic resonance which spans 1.6 GHz and 2.4 GHz for 10B and 11B, respectively. For each frequency step a measurement with (DP = double pass) and without (SP = single pass) retroreflected beam is taken by switching a remote-controlled shutter placed in front of the reflecting mirror. The dwell time in each spectrum is adapted between 1 and 10 s according to ionization laser power, resonant laser power and boron source performance so that a count of more than 200 photoions per step is achieved at a peak. After recording a resonance signal for 11B, which has 80% natural abundance, laser and QMS are tuned to settings for 10B and a spectrum of this isotope is recorded under identical experimental conditions but with increased dwell time in SP configuration to improve statistics at the lower abundance. The intersection angle $\alpha$ of atom and laser beam is then deliberately varied around $\pi/2$ before another set of spectra is recorded.

Single-pass spectra are fitted using Voigt profiles for the individual hyperfine components and shared widths for all hyperfine components. We take into account the natural Lorentzian linewidth of $40.1$ MHz and Gaussian broadening effects originating from the thermal atom beam. In the 11B DP spectra, the transition is probed and recorded twice with a small shift in frequency $\Delta\nu_{\mathrm{DP}}^{11}$ since the intersection angle has never been exactly $\pi/2$ and thus the rest-frame transition frequency for the two laser beam orientations is Doppler-shifted in opposite directions. To reduce the number of degrees of freedom in the fit with back-reflection, the hyperfine factors are fixed to their literature value for the $2p$ state and the values extracted from a single-pass fit for the $3s$ level. It is assumed that the peaks from both directions each share the same width but have different intensity due to power losses and differences in beam diameter after retro-reflection. Since the single-pass spectrum was recorded simultaneously using the shutter system, one of the transition centroids for the double-pass spectra was fixed to the value from the single-pass. With these constraints, the model for the double-pass spectra converged with good agreement. An example dataset with the respective fit results is shown in Fig. 2. The spectra of 11B and 10B are taken under the same ambient and spatial conditions, the transition frequency centroids are nevertheless Doppler-shifted with different magnitudes due to their difference in the mass-dependent mean velocity. The size of this effect can be determined from the double-pass spectrum of 11B and scaled for the thermal velocities: since the atoms are emitted from a thermal source, both isotopes inherit the same average kinetic energy, resulting in a velocity ratio of $r=v(^{11}\mathrm{B})/v(^{10}\mathrm{B})=\sqrt{m(^{10}\mathrm{B})/m(^{11}\mathrm{B})}.$ The isotope shift between the stable isotopes can therefore be calculated as

[TABLE]

where $\nu_{\mathrm{SP}}(^{A}\mathrm{B})$ is the respective Doppler-shifted transition centroid from the single-pass spectra of 11B and 10B. With this approach, the determination of the isotope shift depends mainly on the parameters extracted from single-pass spectra, which are determined with much smaller systematic and statistical uncertainties.

Combining all data yields a value for the isotope shift of -5031.3 (2.0) MHz. Table 2 lists previous measurements of the isotope shift as well as our theoretical and experimental results. None of the previous publications denotes the isotope shift explicitly but report transition frequencies of the two isotopes with respective uncertainties. The influence of correlated uncertainties in the calculation of $\delta\nu_{\mathrm{IS}}$ is therefore unknown and Gaussian error propagation is used. Our value agrees within $1.3\sigma$ with the most recent value obtained by Johansson Johansson et al. (1993) but has about two orders of magnitude higher accuracy. It is close to the calculated mass shift of $-5023.00(8)$ MHz and the difference of $8.8\,(2.0)$ MHz can be attributed to the finite size effect. This corresponds to a change in mean-square nuclear charge radius of $\delta\langle r_{\mathrm{c}}^{2}\rangle=-0.49\,(12)\,$fm2 taking the calculated field shift factor $C=16.91\,(9)\,$MHz/fm2 into account. At first glance it might seem counterintuitive that the charge radius of 10B is larger than that in 11B, but it is well in line with similar results in the lithium and beryllium chains Ewald et al. (2004); Sánchez et al. (2006); Nörtershäuser et al. (2009); Krieger et al. (2012). There, a minimum of the charge radius has been observed at $N=6$ and was attributed to the cluster structure of the lighter nuclei. Similar arguments can also be used for boron: 9B is unbound since it consists of two $\alpha$-clusters and an additional proton that does not support binding of the two $\alpha$’s like the neutron does in the case of 9Be. 10B can be considered to be of $\alpha+\alpha+d$ structure and is expected to have a rather large charge radius and the decrease towards 11B would be analog to the decrease of $\delta\langle r_{\mathrm{c}}^{2}\rangle=-0.734\,(40)\,$fm2 from 6Li ($\alpha+d$) to 7Li ($\alpha+t$) Nörtershäuser et al. (2011), which is indeed of similar size. In Fig. 3 the charge radius of 11B is plotted versus the charge radius of 10B. Apart from these more empirical arguments we compare the extracted change of the charge radius to new ab initio calculations of the stable boron isotopes in the following.

Nuclear Structure Theory—We employ two state-of-the-art ab initio nuclear structure methods, the no-core shell model (NCSM) and the Green’s function Monte Carlo (GFMC) approach, to compute the charge radii of 10B and 11B.

For the NCSM we use different two-nucleon (NN) and three-nucleon (3N) interactions from chiral effective field theory: (a) the N2LO-SAT interaction with NN and 3N interaction at N2LO with flow parameter $\alpha=0.08\,\text{fm}^{4}$ Ekström et al. (2015); (b) the NN interaction at N3LO by Entem and Machleidt Entem and Machleidt (2003) supplemented with a 3N interaction at N2LO with local regulator and reduced cutoff ($\Lambda_{\text{3N}}=400\,\text{MeV/c}$, $\alpha=0.08\,\text{fm}^{4}$) that has been widely used in the past years Roth et al. (2012); Hergert et al. (2013); Gebrerufael et al. (2017); (c) the Entem-Machleidt NN interaction at N3LO with a new 3N interaction at N2LO with nonlocal regulator ($\Lambda=500\,\text{MeV/c}$, $c_{D}=0.8$, $\alpha=0.12\,\text{fm}^{4}$); and (d) the recent NN interaction at N4LO by Entem, Machleidt, and Nosyk Entem et al. (2017) plus a 3N interaction at N2LO with nonlocal regulator ($\Lambda=500\,\text{MeV/c}$, $c_{D}=-1.8$, $\alpha=0.16\,\text{fm}^{4}$) Hüther et al. . Only (a) uses information beyond the few-body sector and explicit constraints on nuclear radii to determine the low-energy constants, in all other cases the NN interaction is fitted exclusively to two-nucleon scattering data and the 3N interaction to the triton binding energy, the triton $\beta$-decay half-life, or properties of the $\alpha$ particle. For all interactions we employ a consistent similarity renormalization group evolution up to the three-body level for the Hamiltonian and up to the two-body level for the radius operator. We have confirmed that the impact of variations of the flow parameter $\alpha$ on the radii is much smaller than the model-space convergence uncertainties. For each interaction large-scale NCSM calculations are performed for model spaces from $N_{\max}=2$ to $10$ using harmonic oscillator frequencies $\hbar\Omega=12,13,...,18$ MeV. To extract the nominal value and uncertainty for the point-proton radius, we first identify the $\hbar\Omega$-value that provides the most stable radius as function of $N_{\max}$ and then use the neighboring $\hbar\Omega$-values and the residual $N_{\max}$-dependence to estimate the many-body uncertainties.

Green’s function Monte Carlo (GFMC) uses imaginary-time projection techniques to solve for nuclear ground- and excited-state energies at the 1% accuracy level for a given Hamiltonian Carlson et al. (2015). Here we employ the AV18+IL7 Hamiltonian, containing the Argonne $v_{18}$ NN potential Wiringa et al. (1995) and Illinois-7 3N potential Pieper et al. (2001); Pieper (2008). The four parameters characterizing this 3N potential were fit to low-lying nuclear spectra in the mass range A=3-10 and give an excellent reproduction of approximately 100 ground- and excited-state energies up to A=12. GFMC calculations successfully predicted the charge radii of 6,8He isotopes, while Li and Be isotope radii tend to be a little smaller than current experimental values Lu et al. (2013).

While our calculated binding energies for 10B agrees well with the experimental value, we slightly underbind 11B. To remedy this shortcoming of the AV18+IL7 interaction, when computing the radius of 11B we slightly quench the phenomenological repulsive term of the 3N force to reproduce the experimental binding energy.

The results for charge radii and ground state energies from both NCSM and GFMC calculations are summarized in Tab. SI of the supplemental material. As in Ref. Lonardoni et al. (2018) and Lu et al. (2013), the charge radii are derived from the estimates for the point-proton radii taking into account the finite-size of the nucleons and the Darwin-Foldy correction. The one-body spin-orbit correction of Ref. Lu et al. (2013) has been estimated in a variational Monte Carlo calculation; it is found to be about five times smaller than the Darwin-Foldy term in 11B, and to vanish in 10B, so it is neglected here. Two-body terms in the charge operator have also been neglected in this analysis, as their contribution was proven to be small in the GFMC calculations of the 12C charge form factor Lovato et al. (2013).

Conclusion and Outlook—

In Fig. 3 the charge radius of 11B is plotted versus that of 10B. Our newly obtained value for the difference in mean-square charge radius $\delta\langle r_{\mathrm{c}}^{2}\rangle=-0.49\,(12)\,\mathrm{fm}^{2}$ is depicted as dotted line with grey shaded uncertainty area. It improves the experimental precision from electron scattering (f) by at least a factor of five. Our value is used as benchmark for ab initio nuclear structure theory. Independent of the applied many-body method and interaction, both NCSM (a-d) and GFMC (e) predict values which are slightly smaller, but are consistent with the experimental value for $\delta\langle r_{c}^{2}\rangle$ within uncertainties. Absolute values of the radii vary with different interactions, which is partly correlated to a variation of the ground-state energies of the two isotopes; for the NCSM calculations with interaction (a) both isotopes are overbound by about 4 MeV compared to experiment, while for interactions (b,c,d) they are all underbound by about 4 MeV. For GFMC the ground-state energies are reproduced by construction and tweaking the 3N potential to match the experimental binding energies.

In summary, the value obtained from high precision laser spectroscopy and novel atomic theory calculations for the 5e- system show good agreement with ab initio nuclear structure calculations. The advances made in all three fields of physics open the window to explore a wider area of light nuclei to be investigated with high precision experiments that will provide valuable benchmarks for understanding and refining our nuclear models.

Acknowledgement—We wish to acknowledge our special thanks to the late Steven C. Pieper, who contributed so much to the progress of nuclear quantum Monte Carlo.

This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 279384907 - SFB 1245 and by the U.S. DOE, Office of Science, Office of Nuclear Physics, under contract DE-AC02-06CH11357. MP and KP acknowledge support from National Science Center (Poland) Grant No. 2014/15/B/ST4/05022. FS and BM acknowledge support from HGS-HIRe. Part of the NCSM calculations were performed on the LICHTENBERG high-performance computer at TU Darmstadt. Under an award of computer time provided by the INCITE program, this research used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory.