Nonlinear interaction of external perturbations in warm dense matter

We present extensive new ab initio path integral Monte Carlo (PIMC) results for an electron gas at warm dense matter conditions that is subject to multiple harmonic perturbations. In addition to the previously investigated nonlinear effects at the original wave number (Dornheim et al., Phys. Rev. Lett., 2020, 125, 085001) and the excitation of higher harmonics (Dornheim et al., Phys. Rev. Res., 2021, 3, 033231), the presence of multiple external potentials leads to mode‐coupling effects, which constitute the dominant nonlinear effect and lead to a substantially more complicated density response compared to linear response theory. One possibility to estimate mode‐coupling effects from a PIMC simulation of the unperturbed system is given in terms of generalized imaginary‐time correlation functions that have been recently introduced by Dornheim et al. (J. Chem. Phys., 2021, 155, 054110). In addition, we extend our previous analytical theory of the nonlinear density response of the electron gas in terms of the static local field correction (Dornheim et al., 2020, Phys. Rev. Lett., 125, 235001), which allows for a highly accurate description of the PIMC results with negligible computational cost.

This interest in the properties of WDM has led to a spark of new developments in the QMC simulation of the UEG at finite temperature, [17][18][19][20][21][22][23][24][25] which culminated in the first accurate parametrization of the exchange-correlation (XC) free energy f xc covering the entire relevant range of densities and temperatures. [8,[26][27][28] These new XC-functionals have already been utilized for thermal DFT simulations of WDM, [29,30] which has further substantiated the impact of thermal effects on material properties for some parameters.
In addition to f xc , accurate data for the warm dense UEG have been presented for a multitude of quantities including the static structure factor S(k), [25,[31][32][33] the momentum distribution function n(k) [34][35][36][37][38] the dynamic structure factor S(k, ), [39][40][41] and other dynamic properties like the conductivity. [42,43] Of particular importance for many practical applications [44][45][46][47] in WDM theory is the response of the UEG to an external perturbation. [2] Typically, such effects are described within linear response theory (LRT), which presupposes a simple, linear relation between response and perturbation and, thus, leads to a drastic simplification of the underlying theory. Therefore, accurate data for the linear response of the UEG based on QMC simulations are available both in the electronic ground state [5,48,49] and at finite temperature. [50][51][52][53][54] Very recently, Dornheim et al. [32,55,56] have gone beyond the assumption of LRT by carrying out extensive ab initio path integral Monte Carlo (PIMC) simulations of a harmonically perturbed electron gas. [48,[57][58][59] First and foremost, this has allowed for the first time to unambiguously check the validity range of LRT. Indeed, it has been reported that nonlinear effects can be important in some situations that are of relevance for state-of-the-art WDM experiments, for example, using free electron lasers. [60] In addition, these PIMC studies have allowed to obtain the first data for various nonlinear density response functions both at the wave vector of the original perturbation [32] and its integer harmonics. [55] This, in turn, has allowed the same group to present an accurate analytical theory for the description of the nonlinear electronic density response in terms of the effectively static local field correction (LFC) G(k), which is available as a readily usable analytical representation. [54] Further new results into this direction include the exploration of the straightforward relation between nonlinear effects and higher order correlation functions known from many-body theory [61] and the computation of the nonlinear density response based on imaginary-time correlation functions defined with respect of the unperturbed system. [62] Yet, all of these works concentrated on the case of a single harmonic perturbation. This is reasonable as long as the superposition principle holds as it allows to decompose any perturbation into a sum of independent spectral components. However, in case of a strong excitation that requires one to go beyond LRT, the superposition principle will fail. Any additional excitation will then not simply add up but the combined effect of the excitations may be more complicated. This is a serious problem for applications because many realistic field-matter interaction experiments utilize combinations of different exciting fields associated with different frequencies or harmonics, for example, refs. [45,63]. Examples are four wave mixing or Raman spectroscopy, including Stokes and Anti-Stokes Raman scattering, stimulated Raman spectroscopy or coherent anti-Stokes Raman spectroscopy.
In many WDM experiments, the situation is similar as often combinations of signals are being used. An example are pump-probe experiments. Also, the plasmon signal in X-ray Thomson scattering (XRTS) of WDM can be viewed as a combination of X-ray and plasmon oscillations. Naturally, this is of no consequence for weak fields, where the superposition principle holds. However, this drastically changes once (at least) one of the involved fields has a high amplitude and nonlinear effects become significant. It should be expected that under these conditions interesting and non-trivial mode-coupling effects will occur.
In the present work, we rigorously study these nonlinear mode-coupling effects by carrying out extensive new PIMC simulations of a warm dense electron gas that is subject to two harmonic perturbations. This generic case should be sufficient to understand the most important nonlinear effects. Remarkably, we find that mode-coupling constitutes the dominant nonlinear effect for weak to moderate values of the perturbation amplitude A. In addition, we practically demonstrate that these effects can be estimated from the imaginary-time structure of the unperturbed system. Finally, we extend our earlier analytical theory of the nonlinear density response of the UEG [55] and find excellent agreement between our LFC-based theory and the exact PIMC results.
For completeness, we note that the mode coupling effects observed in our paper should not be confused with mode coupling theory (MCT). MCT means linear superposition of field modes, and this term is used extensively in electromagnetism, nanophotonics and other fields which are based on Maxwell's equations the linearity of which implies the superposition principle of individual modes. Deviations occur only in the case that the field penetrates a nonlinear medium. This is, in fact, the case of WDM.
The paper is organized as follows: In Section 2, we introduce the relevant theoretical background including the model system (Section 2.1) and different approaches to the density response (Section 2.2). Section 3 is devoted to the discussion of our new simulation results, and the paper is concluded by a brief summary and outlook in Section 4.

THEORY
We assume Hartree atomic units throughout this work. Furthermore, we restrict ourselves to a fully unpolarized (paramagnetic) system where the number of spin-up and -down electrons are equal, that is, N ↑ = N ↓ = N/2.

Model system
The UEG is typically characterized by two dimensionless parameters, which are both of the order of unity in the WDM regime [64] : (a) the density parameter (also known as Wigner-Seitz radius) r s = r∕a B , where r and a B are the average inter-particle distance and first Bohr radius, and (b) the degeneracy temperature = k B T/E F , with E F = k F 2 /2, and the Fermi wave number being defined as From a physical perspective, these two parameters allow for a straight forward interpretation: for ≫ 1, quantum effects are negligible and the system attains the classical limit, whereas ≪ 1 indicates the electronic ground state, where quantum degeneracy effects predominate; regarding the Wigner-Seitz radius, r s ≪ 1 indicates the weak coupling regime, and the UEG actually converges towards the ideal Fermi gas in the limit r s → 0. In contrast, r s ≫ 1 indicates strong coupling, where the UEG first transforms into an electron liquid [51,53] and subsequently forms a Wigner crystal. [65] Thus, the density parameter plays the role of an effective coupling parameter for the UEG, which can be easily seen by considering the kinetic and interaction contributions K and W in their respective lowest order: the Hartree-Fock kinetic energy is given by K = a HF ∕r 2 s , and the mean-field part of the interaction scales as W = b/r s . The usual definition of a coupling parameter thus directly results in Γ = W/K ∼ r s .
Let us consider a Hamiltonian of the formĤ whereĤ UEG is the standard Hamiltonian of the unperturbed UEG [8] and the external potential V(r) is comprised of multiple harmonic perturbations, which immediately leads to the Fourier transform In practice, we simulate the inhomogeneous electron gas that is governed by Equation (2) using the direct PIMC method [66,67] without any nodal constraints. Therefore, our simulations are computationally involved due to the fermion sign problem, [68][69][70] but exact within the given statistical uncertainty.

Density response
The induced density contains then a variety of combinations of higher harmonics of the incoming perturbing potential ind (q) = (q) Specifically, (q) denotes the usual static limit of the linear response function, [2,71] Y (q 1 , q 2 ) the generalized quadratic response function, and Z(q 1 , q 2 , q 3 ) the generalized response function in cubic order; see ref. [55] for a more detailed derivation.
In addition, we note that the generalized quadratic density response function (see ref. [62] for details) is connected to the imaginary-time structure of the system by the relation Here we define the density operator̃( to be not normalized, and r l, denotes the position of particle l at an imaginary time ∈ [0, ]. Equation (6) thus directly implies that all quadratic terms of the nonlinear density response (including mode-coupling effects, see below) can be obtained from a single simulation of the unperturbed UEG. Moreover, a similar relation exists for the cubic response function Z(q 1 , q 2 , q 3 ), but it is not used in the present work. While Equation (6) is exact, it is highly desirable to have an accurate theory for the density response that can be evaluated without the need for a computationally expensive PIMC simulation. To this end, we follow the considerations from ref. [55], and give an (approximate) mean-field expression for the quadratic density response, taking the form. [72] ] .
Moreover, the screening terms in the denominator of Equation (8) can be improved by including electronic XC-effects in the form of the effectively static LFC G(k), which leads to ] .
In practice, we use the analytical representation of G (k; r s , ) from ref. [54].

RESULTS
All PIMC simulation results that are presented in this work have been obtained using a canonical adaption [73] of the worm algorithm by Boninsegni et al. [74,75] . Further, we use a primitive factorization of the density matrix, and the convergence with the number of imaginary-time steps P has been carefully checked. Let us start this investigation by considering the full wave-number dependence of the density response for the case of two harmonic perturbations at r s = 2 and = 1 shown in Figure 1. We note that this parameter combination corresponds to a metallic density (e.g., aluminium [76] ) in the WDM regime, and has frequently been studied in previous   Response for a single perturbation with q = q 1 (q = q 2 ), taken from ref. [55] investigations [32,55,56,61,62,77,78] of the nonlinear electronic density response. The left panel has been obtained for the case of equal perturbation amplitudes, A 1 = A 2 = A, and for the two wave numbers q 1 = (1, 0, 0) T 2 /L, and q 2 = (2, 0, 0) T 2 /L. The yellow triangles depict the prediction from LRT, which is given by two independent signals at q 1 and q 2 ; any interplay between multiple perturbations is completely neglected.
The black squares have been obtained from our PIMC simulations with a comparably small perturbation amplitude A = 0.02. In this case, LRT is relatively accurate and the density response is close to zero for all k ≠ q i , although small yet significant signals can be seen at k = 3q 1 and k = 2q 2 . The latter effect is substantially increased for A = 0.1, as can be seen from the red circles in the same plot. Finally, the green crosses have been obtained for the case of a strong perturbation, A = 0.5, and the depicted spectrum of the density response strongly deviates from LRT and exhibits at least five different signals.
To understand these nontrivial findings, we have to evaluate the different contributions to the mode-coupled density response given in Equation (5). Firstly, the first line corresponds to the usual LRT contributions at the original wave numbers q i , which are always present. The second part corresponds to all quadratic contributions, which, in turn, can be further sub-divided into two distinct categories: (i) the case i = j gives the quadratic density response at the second harmonic of the respective perturbation. This, too, is not a mode-coupling effect and has been previously reported in refs. [55,62]; (ii) the case i ≠ j entails the mode-coupling signal, and predicts a quadratic density response for ±k = q 1 ± q 2 . For the present case, this results in signals at k = q 1 (which comes in addition to the LRT term) and at k = 3q 1 . In fact, the latter constitutes the dominant nonlinear signal at these parameters, which, too, can be easily understood from Equation (5). Following the notation from ref. [55], the density response at the second harmonic of a perturbation at wave number q is given by the quadratic density response function (2) (q) = Y (k − q, q) k,2q . The mode-coupling response at k = 3q 1 is instead given by (2) and thus entails two evaluations of the generalized quadratic response function, namely Y (q 2 , q 1 ) and Y (q 1 , q 2 ). The last effect from the left panel of Figure 1 that requires an explanation is the signal at k = 5q 1 that emerges for strong perturbations, A = 0.5. Evidently, this has to be a mode-coupling effect, as it is not at an integer harmonic of q 2 and would correspond to the fifth harmonic of q 1 , which is negligible for such parameters. [55] More specifically, it corresponds to a cubic contribution to Equation (5) at k = 2q 2 + q 1 that is described by the generalized cubic density response function Z.
Let us next consider the right panel of Figure 1, where we compare the wave number dependence of the density response from different scenarios. Again, the yellow triangles depict the linear response and have been included as a reference. The red circles correspond to the doubly perturbed simulation for A = 0.1, that is, the same as the red circles in the left panel. Furthermore, the black squares and green crosses have been obtained from PIMC simulations with only a single perturbation at q 1 and q 2 , respectively, and are taken from ref. [55]. For the black squares, the signal at the original perturbation can hardly be distinguished from the LRT prediction, whereas the corresponding red circle exhibits a substantially reduced response. Naturally, the latter is a direct consequence of the quadratic mode coupling, as q 1 = q 2 − q 1 in this case. In addition, we find a weak signal at k = 2q 1 for the black circles, which is the incipient quadratic response at the second harmonic. The green crosses exhibit a similar behaviour, although the signal at k = q 2 is not as strongly reduced for the double perturbation compared to the other points. This makes sense, as the main deviation of the red circles from LRT is given by the aforementioned contribution from the second harmonic of q 1 , which is smaller in magnitude than the mode-coupling effect at k = q 1 .
In addition, we find a signal at the second harmonic of q 2 , which is similarly pronounced for both the green crosses and the red circles. Finally, we stress that the signal at k = 3q 1 is only present in the simulation of two external perturbations, as it is expected. We thus briefly summarize the following interim conclusions: mode-coupling (i) typically constitutes the strongest nonlinear effect as it is comprised of multiple contributions at the same wave number k (cf. Equation 10) and (ii) leads to signals at wave numbers where, otherwise, the response would vanish.
To provide a more quantitative assessment of the effects due to mode-coupling, we investigate the perturbation amplitude dependence of the density response in Figure 2. In particular, the left panel corresponds to the case of k = q 1 , that is, at one of the original perturbations, and the red circles show our new PIMC data that have been obtained for the case of two perturbations for different A. As a reference, we also include PIMC data from ref. [55] that have been obtained for a single perturbation at q 1 as the green crosses. The solid yellow line shows the prediction from LRT, which is in good agreement to both PIMC data sets for very small A. Strikingly, LRT holds for substantially larger A in the case of a single perturbation compared to the more complicated behaviour due to mode-coupling between q 1 and q 2 . In fact, the dominant nonlinear effect in the green crosses is given by the cubic density response at the first harmonic, which is depicted by the dashed green line, and is in excellent agreement to the data points. In stark contrast, the first nonlinear term in the red circles is quadratic in the perturbation amplitude A and determined by the generalized response function Y (−q 1 , q 2 ) = Y (q 2 , −q 1 ). The latter can be conveniently estimated by performing a PIMC simulation of the unperturbed electron gas and evaluating Equation (6); see ref. [62] for a detailed description of this procedure. The final result for both LRT and the thus obtained quadratic mode-coupling term is depicted by the dash-dotted black curve, which is in excellent agreement to the independent data points for A ≤ 0.2. For larger A, cubic terms in Equation (5) start to become noticeable, as it is expected.
The right panel of Figure 2 shows the same information, but for k = q 1 + q 2 . In this case, the green crosses exhibit a very weak response that is well described by the cubic response of the third harmonic of q 1 , that is, the dashed green curve. In stark contrast, the red circles are larger by an order of magnitude and exhibit a parabolic behaviour (cf. Equation 10) in leading order that is again well described by the theoretical curve that we obtain by evaluating Equation (6). In addition, the solid blue curve has been obtained by evaluating Equation (9) using as input the static LFC by Dornheim et al. [33] Evidently, the resulting parabola is in excellent agreement with the exact A → 0 limit predicted by the ITCF formalism. This, in turn, further corroborates the high value of the LFC for the description of nonlinear effects that has been reported in a previous work. [55] A further interesting topic is the behaviour of the density response when the perturbation amplitudes of the two perturbations are not equal. Such a case is investigated in Figure 3, where we consider the same wave vectors q 1 and q 2 as before. More specifically, we show the density response at k = q 1 , and the x-axis shows the perturbation amplitude A 1 = A; the perturbation amplitude A 2 = 0.5 is being kept constant. The resulting data are shown as the black squares, and, approximately, exhibit a linear behaviour with A. At the same time, the pre-factor substantially differs from the LRT prediction (solid yellow), and from the simulation results obtained for only a single perturbation (green crosses). In other words, we find an apparently linear effect that is not described by LRT.
Yet, this seeming contradiction is resolved by the second line in Equation (5), which predicts a mode-coupling contribution at k = q 1 that is proportional to A 1 × A 2 . Since A 2 is being kept constant for the results shown in Figure 3, the quadratic mode-coupling term manifests with a linear dependence on A 1 = A, which resolves this apparent conundrum. The corresponding estimation combining both LRT and the quadratic term that has been obtained from the generalized ITCF (Equation 6) is depicted as the dash-dotted black line in Figure 3, and qualitatively captures the correct trend. The comparably small residual differences between data and theory are due to the fact that at A 2 = 0.5, even the quadratic LFC ITCF

F I G U R E 2
Dependence of the density response at k = q 1 (left) and k = q 1 + q 2 (right) on the perturbation amplitude A for the same conditions as in Figure 1. Red circles: PIMC results for the density response at k; dash-dotted black: A → 0 limit predicted by the ITCF, Equation (6); green crosses: PIMC results for the density response at k from a simulation with a single harmonic potential at q 1 taken from ref. [55]; dashed green: corresponding description based on the cubic density response; solid yellow: LRT; solid blue: LFC-based description of the A → 0 limit from Equation (9 perturbation at q 2 = (2, 0, 0) T 2 /L is applied with A 2 = 0.5 being fixed. Black squares: PIMC data for the density response at k = q 1 ; dash-dotted black: ITCF prediction, Equation (6); dotted red: linear fit for A ≤ 0.2; green crosses: PIMC results for a single external perturbation at q 1 taken from ref. [55]; solid yellow: LRT prediction for the latter case; dashed green: corresponding LRT + cubic response description is not fully sufficient; see Figure 2 above. Finally, the dotted red line has been obtained from a linear fit to the PIMC data for A ≤ 0.2, and, therefore, depicts the correct nonlinear result for the linear response in this case. We stress that this finding might be of considerable practical relevance, as it directly implies that the presence of a second external perturbation with a large amplitude leads to a break-down of LRT for the first perturbation even in the limit of A → 0. For example, an external ionic potential might be expanded into such harmonics, and the corresponding coefficients might be substantial. Let us next consider the impact of the direction of the external perturbation. To this end, we investigate the wave number dependence of the density response again for r s = 2 and = 1 in Figure 4. Specifically, the black squares and green crosses have been obtained for the case of single perturbations along q 1 and q 2 , respectively, with A = 0.1. In addition, the blue diamonds show the signal for the case of a double perturbation at both q 1 and q 2 with A 1 = A 2 = A. Finally, the red circles have been obtained from a new simulation with two perturbations along q 1 and q 3 = (0, 2, 0) T 2 /L. We note that the k-vectors shown in Figure 4 have been selected along the x-direction insofar as this is possible. Therefore, the red circles are indistinguishable from the black squares as no mode-coupling affects the density response in this direction.
This notion is further validated in Figure 5, where we show the A-dependence of the density response at k = q 1 . The blue diamonds show the PIMC data for the simulation with two perturbations along q 1 and q 2 , which are well described by LRT and the quadratic mode-coupling response (dash-dotted black curve) as it has been explained before. The green crosses are the simulation results for only a single perturbation at q 1 , which follow the combination of LRT and the cubic response at the first harmonic. Finally, the red circles have been obtained from the simulations with perturbations at q 1 and q 3 . Evidently, the response at k = q 1 is hardly affected by the presence of the second perturbation along another direction, and the circles closely follow the green crosses for small to medium perturbations. For completeness, we mention that the small deviations at large perturbation amplitudes A are a consequence of cubic mode-coupling terms, cf. Equation (5).
Let us conclude our study of mode-coupling effects in the density response of the warm dense UEG with a more complicated combination of perturbation wave vectors. This situation is investigated in Figure 6, where we consider the case of q 1 = (2, 0, 0) T 2 /L and q 2 = (3, 0, 0) T 2 /L, again with A 1 = A 2 = A. In particular, the qL/2 are given by prime numbers, so that no quadratic mode-coupling contribution coincides with the original wave numbers or their integer harmonics. As usual, the yellow triangles depict the prediction from LRT, and the black squares, red circles, and green crosses show PIMC results for A = 0.05, A = 0.1, and A = 0.5, respectively. We note that we have divided the actual density response by so that the individual data points give an intuitive measure for the respective contribution of each wave number to the total density response of the system. Evidently, this choice of the q i results in a number of signals at different wave vectors k even in the case of a rather weak perturbation, A = 0.05. More specifically, there appear four quadratic signals (the second harmonics of q 1 and q 2 , and the mode-coupling signals at k = q 2 ± q 1 ) in addition to the two main peaks at the original wave vectors of the external perturbation. In the regime of large perturbations, A = 0.5, the situation becomes even more complicated as additional signals occur due to cubic effects. We thus conclude that the density response to a superposition of multiple external perturbations becomes dramatically more complicated beyond the LRT regime.

SUMMARY AND DISCUSSION
In summary, we have analysed nonlinear mode-coupling effects in the density response of an electron gas in the WDM regime. This has been achieved on the basis of extensive new ab initio PIMC results that have been obtained by applying multiple external harmonic perturbations at the same time. First and foremost, we note that mode-coupling effects, while being absent in LRT, constitute the dominant nonlinear effect for weak to moderate perturbation amplitudes. In addition, the presence of a strong, constant perturbation at q 2 leads to a substantial modification of the linear response at a different wave vector q 1 even in the limit of A → 0, which has considerable implications for a non-UEG in an external potential, for example due to positively charged ions. Another important application will be pump-probe experiments. Our results indicate that the effect of even a very weak probe signal cannot be correctly described within the common linear transport and optics schemes, in case of a high pump amplitude. From a theoretical perspective, we have shown that mode-coupling effects can be estimated numerically exactly from a PIMC simulation of the unperturbed UEG by computing the generalized imaginary-time correlation functions that have recently been introduced by Dornheim et al. [62] This leads to a substantial reduction of the computational effort, as, in principle, the full information about mode-coupling between all wave vectors can be obtained from a single simulation.
In addition, we have extended our earlier analytical theory [55] for the nonlinear density response in terms of the static LFC, [33,50,54] and find excellent agreement between theory and simulations with negligible computational cost.
There exists a large variety of nonlinear response phenomena in WDM. An interesting question is to inquire about the physical origin of the nonlinearity. For example, the Coulomb interaction between the electrons automatically gives rise to nonlinearities: the 1/r distance dependence of the Coulomb potential transforms a harmonic perturbation of any particle into a nonlinear response of its neighbours. A familiar example is the generation of high harmonics of a laser field in a gas. [79] Therefore, our results will be particularly important for situations where the plasma is strongly correlated. While we have studied only moderately coupled WDM (r s = 2), the consequences will be even more significant when r s is increased. At strong correlations a number of additional nonlinear effects should be expected. For example, modes propagating in different directions may couple if the isotropy of a nonideal system is broken, for example, in the presence of a magnetic field. Such an effect is known from strongly couple classical plasmas where it leads to nontrivial coupling of transport processes in different directions, including diffusion [80] and heat transport. [81] Similar nonlinear effects that are mediated by strong Coulomb interaction should also be expected in WDM.
We expect that our new results will open up many avenues for future research in the field of WDM theory and beyond. First and foremost, we mention the potential utility of nonlinear effects in the electronic density response as a method of diagnostics in experiments. Specifically, the interpretation of current XRTS experiments is solely based on LRT, [45] and the inference of important plasma parameters such as the electronic temperature T e is notoriously difficult. [44] In this context, Moldabekov et al. [82] have recently suggested that the application of an external perturbation to generate inhomogeneous WDM samples on purpose leads to a modified experimental signal, that more strongly depends on T e . Naturally, the accurate theoretical description in this case will have to consider the mode-coupling between this external potential and the XRTS probe, and the current LFC-based theory is uniquely suited for this endeavour. A second potential application of our work is given by the possibility to use the nonlinear density response as a probe for three-and even four-body correlation functions known from many-body theory, [61] which might give unprecedented insights into the physical mechanisms of WDM. Finally, we mention that a general theory of the nonlinear electronic density response can be directly incorporated into many theoretical approached that have hitherto, often by necessity, been limited to LRT. Prominent examples include the construction of effectively screened ionic potentials [83][84][85] and the electronic stopping power. [47] Specifically, going beyond LRT is a crucial step for understanding and exploring the effective ion-ion attraction in the media with strongly correlated electrons, [86,87] which remains to be an unsolved problem. [88,89] Additionally, the effect of electronic strong correlations on a projectile energy dissipation (stopping power) is one of the long-standing problems in WDM, [90,91] which can be addressed using the non-LRT. [72] ACKNOWLEDGMENT This work was partly funded by the Center of Advanced Systems Understanding (CASUS) which is financed by Germany's Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament, and by the Deutsche Forschungsgemeinschaft (DFG) via project BO1366/15. The PIMC calculations were carried out at the Norddeutscher Verbund für Hochund Höchstleistungsrechnen (HLRN) under grant shp00026, on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden, on the cluster hemera at Helmholtz-Zentrum Dresden-Rossendorf (HZDR), and at the computing center (Rechenzentrum) of Kiel university. Open Access funding enabled and organized by Projekt DEAL.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.