Universality of continuous phase transitions on random Voronoi graphs

TL;DR

This study investigates phase transitions of the Ising model and contact process on 2D random Voronoi graphs, finding that quenched disorder does not alter their universal critical behavior, unlike in other lattice types.

Contribution

It demonstrates that phase transitions on random Voronoi graphs belong to the same universality classes as on regular lattices, highlighting the irrelevance of quenched disorder in this context.

Findings

Both models exhibit critical behavior consistent with their clean universality classes.

Critical points are determined with high precision.

First correction-to-scaling exponent for the Ising model is reported.

Abstract

The Voronoi construction is ubiquitous across the natural sciences and engineering. In statistical mechanics, though, critical phenomena have so far been only investigated on the Delaunay triangulation, the dual of a Voronoi graph. In this paper we set to fill this gap by studying the two most prominent systems of classical statistical mechanics, the equilibrium spin-1/2 Ising model and the non-equilibrium contact process, on two-dimensional random Voronoi graphs. Particular motivation comes from the fact that these graphs have vertices of constant coordination number, making it possible to isolate topological effects of quenched disorder from node-intrinsic coordination number disorder. Using large-scale numerical simulations and finite-size-scaling techniques, we are able to demonstrate that both systems belong to their respective clean universality classes. Therefore, quenched disorder introduced by the randomness of the lattice is irrelevant and does not influence the character of the phase transitions. We report the critical points of both models to considerable precision and, for the Ising model, also the first correction-to-scaling exponent.

Figures

Tables

Table 1. Table 1: Results of the simultaneous fits according to Eq. ( 7 ).

$L_{\text{min}}$	${\left(\xi /L\right)}^{\ast }$	$U_4^{\ast }$	$U_6^{\ast }$	$\omega$	${\chi }^2$/d.o.f
16	0.9078(2)	0.61067(2)	1.4563(2)	1.36(2)	12.6
18	0.9070(2)	0.61066(2)	1.4564(2)	1.43(2)	7.2
20	0.9066(3)	0.61066(2)	1.4564(2)	1.47(3)	5.8
24	0.9062(3)	0.61065(2)	1.4564(2)	1.53(4)	5.0
32	0.9058(3)	0.61065(3)	1.4563(2)	1.59(7)	4.7
40	0.9060(5)	0.61067(4)	1.4561(3)	1.50(12)	4.1
48	0.9063(7)	0.61069(5)	1.4559(4)	1.37(18)	4.5

Table 2. Table 2: Number of disorder replicas and cluster updates per EMCS for the Ising simulations at criticality. All systems were simulated at T = 1.47205 𝑇 1.47205 T=1.47205 .

	$L$	$N_r$	$n_{\mathrm{Wolff}}$
	16, 18, 20, 24	100 000	10
	32, 36, 40, 48	100 000	12
	64, 80, 96	100 000	15
	128, 192, 256	35 000	17
	384, 512	15 000	23
	768	5000	27
	1024	4000	29

Equations

$$\mathcal{H}=-\sum _{⟨i,j⟩}{J}_{ij}{s}_i{s}_j+\sum _i{h}_i{s}_i,s_i=\pm 1$$

$$[\mathcal{O}{]}_{\text{avg}}\equiv \frac{1}{N_r}\sum _{i=1}^{N_r}{\mathcal{O}}_i$$

$$M$$

$$E$$

$$\chi$$

$$\xi =\frac{1}{2\sin (k_{\text{min}}/2)}\sqrt{\frac{[⟨|{\mathcal{F}}^2(0)|⟩{]}_{\text{avg}}}{[⟨|{\mathcal{F}}^2({\mathbf{k}}_{\text{min}})|⟩{]}_{\text{avg}}}-1},$$

$$\mathcal{F}(\mathbf{k})=\sum _j{s}_j\exp (i{\mathbf{kx}}_j),$$

$$\displaystyle U_{4}=\Bigg{[}1-\frac{\langle m^{4}\rangle}{3\langle m^{2}\rangle^{2}}\Bigg{]}_{\text{avg}},\quad U_{6}=\Bigg{[}\frac{\langle m^{6}\rangle}{\langle m^{2}\rangle^{3}}\Bigg{]}_{\text{avg}}.$$

$$R{|}_{Q_{\xi }=s}=R^{\ast }+a_R{L}^{-\omega },$$

$$\omega =1.54(16)$$

$$(\xi /L{)}^{\ast }$$

$$U_4^{\ast }$$

$$U_6^{\ast }$$

$$T{|}_{Q_{\xi }=s}=T_c+a{L}^{-\omega -1/\nu },$$

$$T_c=1.4720633(31).$$

$$[⟨m⟩{]}_{\text{avg}}$$

$$\chi$$

$$U_4$$

$$x\equiv (T-T_c)L^{1/\nu }.$$

$$p_c=0.649788(1)$$

$$\rho (L,t)=t^{-\delta }\tilde{\rho }\left(t/L^z\right),$$

$$\rho (\Delta ,t)=t^{-\delta }\hat{\rho }\left(\Delta t^{1/{\nu }_{\parallel }}\right),$$

Full text

I Introduction

Many location optimization problems can be approached through area-of-influence considerations. A simple example is that of several fire stations distributed over a large city. Rather naturally, the area of responsibility attributed to a particular fire station should include those buildings which lie closer to it than to any other station. The resulting tessellation of the city map defines the so-called Voronoi graph (VG). In Fig. 1, the red dots would denote fire stations, with the corresponding VG being depicted in green. Due to its conceptual simplicity, Voronoi constructions can be found in a large number of applications spanning all fields of physical sciences, including climate modeling Ringler et al. (2008); Ju et al. (2011), crystal structure Wigner and Seitz (1933), cosmology Weygaert et al. (1989); Ramella et al. (2001), microbiology Poupon (2004), and growth processes William and Mehl (1939), as well as optimization problems Okabe and Suzuki (1997), game theory Ahn et al. (2004); Cheong et al. (2004); Dürr and Thang (2007), artificial intelligence Forbus et al. (2002); Gavrilova (2008) and, recently, also in the field of machine learning Ward et al. (2017), among others. Moreover, numerous generalizations have been defined, such as weighted graphs, and Voronoi graphs on spherical and general curved surfaces, Na et al. (2002); Boissonnat et al. (2006), as well as for fuzzy point sets Jooyandeh et al. (2009) or metrics other then Euclidean 111In a more realistic fire station example one might consider using the Manhattan distance rather than the Euclidean one.

In this paper, the VG is used as a discrete spatial lattice on top of which the physical system is evolved. During the evolution the lattice remains static, i.e., we consider quenched disorder Binder and Young (1986). As the graph is constructed from a randomly distributed set of points it is said to present topological randomness, as opposed to conventional disorder, as e.g., in the case of randomly diluted regular structures.

Quenched disorder in general constitutes a relevant perturbation on the phase transition, i.e., it may change the universal properties of the model Nishimori and Ortiz (2010). The first major result on disorder relevance was brought forward by Harris in 1974 Harris (1974), who argued that disorder is relevant when $d\nu<2$, where $d$ is the dimension of the system and $\nu$ denotes the correlation length exponent. Using Harris’ criterion, numerical results on 2D diluted regular lattices can be entirely explained, including the rise of non-conventional activated scaling and strong Griffiths effects for the contact process (CP) Vojta and Dickison (2005); Vojta et al. (2009); Wada and de Oliveira (2017), as well as ambiguities related to strong logarithmic corrections for the Ising model Kim and Patrascioiu (1994); Ziegler (1994); Kühn (1994); Dotsenko and Dotsenko (1981); Jug (1983); Shankar (1987); Martins and Plascak (2007); Fytas and Theodorakis (2013); Zhu et al. (2015), which represents the marginal case, as in two dimensions $\nu=1$. Despite its success in describing the effects of uncorrelated disorder, the Harris criterion fails to explain certain results for the Ising model and CP on two-dimensional Delaunay triangulations Janke et al. (1993, 1994); Lima et al. (2000); Janke and Villanova (2002); de Oliveira et al. (2008, 2016), where both systems retain their clean universal properties. In order to explain these results, a generalization has been proposed Barghathi and Vojta (2014), which attributes the non-relevance to strong spatial anti-correlations in the coordination numbers of the lattice nodes. Also for this new criterion, however, violations have been found Schrauth et al. (2018a). Moreover, as the argument specifically relies on fluctuations in the local coordination numbers, this criterion is silent about random structures with constant coordination number, such as the VG, where each site has exactly three neighbors (compare Fig. 1).

Therefore, finding a complete criterion for the influence of topological disorder on continuous phase transitions still remains an unsolved puzzle, to which the study of critical phenomena on constant coordination lattices can provide valuable new pieces. One such puzzle piece is given in the current work, where numerical simulations of the Ising model and the CP on two-dimensional VGs show clean universal behavior for both systems, revealing the perturbations by quenched topological randomness to be irrelevant for these systems.

The paper is organized as follows. In Sec. II we review the construction rules of the VG. In Sec. III and IV we present and analyze numerical simulations of the Ising model and the CP, respectively. Finally, in Sec .V we present our concluding remarks.

II Voronoi Graph

The Voronoi diagram is a partition of the plane into cells, generated by a set of points $P=\{p_{1},p_{2},\dotsc,p_{N}\}$ such that for each cell corresponding to the point $p_{i}$, every point $q$ in that cell is closer to $p_{i}$ than to any other point $p_{j}$, i.e., $d(q,p_{i})<d(q,p_{j})\,\,\forall p_{j}\neq p_{i}$. It is the geometric dual of the Delaunay triangulation, which can be constructed by connecting points corresponding to adjacent Voronoi cells. The VG is defined by taking as sites the corners of the cells and, as edges, the boundaries between the cells. Therefore the new set of points $P^{\prime}={p^{\prime}_{1},p^{\prime}_{2},\dotsc,p^{\prime}_{2N}}$ is twice as large as the original set that defined the cells. This can be easily seen from the Euler characteristic of a finite graph, which is defined as $\chi=N-E+F$, where $N$, $E$ and $F$ are the number of vertices, edges and faces. For the periodic boundary conditions used here $\chi=0$ holds. In a Delaunay triangulation, $E=3N$ as the average coordination number is exactly $q=6$ and any edge is shared by two triangles. Therefore we end up with $F=2N$ faces in the triangulation and hence $2N$ points in the VG due to the duality property. The location of the VG sites is given by the center of the circumcircle of the corresponding triangle in the Delaunay triangulation. A sample of a periodic VG is shown in Fig. 1. It can be easily seen that all cells have convex shape and that every site has exactly three neighbors. This latter feature, the absence of coordination number fluctuations, constitutes a major motivation of this study. We build VGs from Delaunay triangulations, and for constructing the initial Delaunay triangulations we use the CGAL library Kruithof (2017).

III Ising Model

The classical Ising model Ising (1925) is defined by the Hamiltonian

[TABLE]

where $s_{i}$ are discrete spins on the lattice sites, $J_{ij}$ denotes the coupling between nearest neighbors $\langle i,j\rangle$ and $h_{i}$ is the external field at site $s_{i}$. For equilibrium lattice models, quenched disorder can be introduced in a variety of ways. For fixed ferromagnetic coupling $J_{ij}=J>0$ and random external field, the system is called random field Ising model (RFIM) and has been investigated thoroughly over the last decades Nattermann (1998). In contrast, for vanishing external field but randomly distributed (anti)ferromagnetic bonds, the system shows the behavior of a spin glass (Binder and Young, 1986, and references therein). In the present study, quenched disorder is introduced as topological randomness encoded in the implicit connectivity of the Voronoi graph. We therefore fix all couplings to $J_{ij}=1$ at vanishing external field.

For the simulation of the Ising model at criticality state-of-the-art importance-sampling Monte Carlo methods are employed, using cluster, as well as local update algorithms. In particular, we use the algorithm proposed by Wolff Wolff (1989), which significantly reduces the critical slowing down and is furthermore straightforwardly applicable to disordered lattices. Although the cluster updates preserve ergodicity, we add local Metropolis updates Metropolis et al. (1953) in order to make sure that the short-wavelength modes are properly thermalized.

In the study of disordered systems, it is necessary to average physical observables over many different, independent disorder realizations of the system, also called replicas. The quenched averages over $N_{r}$ replicas are performed at the level of (extensive) observables, rather than at the level of the partition function Binder and Young (1986). Denoting quenched averages as

[TABLE]

and thermal averages as $\langle...\rangle$, we define magnetization, energy and susceptibility as

[TABLE]

respectively, where $m$ and $e$ denote the magnetization and energy per site. Furthermore, the two-point finite-size correlation function is given by

[TABLE]

with the Fourier transform of the magnetization being defined by

[TABLE]

where $\mathbf{x}_{j}$ denotes the spatial coordinate of site $j$ and $\mathbf{k}_{\text{min}}=(2\pi/L,0)$ represents the smallest non-zero wave vector in the finite lattice. The time series of the quantities $m$, $e$ and $\mathcal{F}^{2}(\mathbf{k}_{\text{min}})$ is recorded during the Monte Carlo run, which allows to compute all relevant observables subsequently. Finally, the fourth- and sixth-order magnetic cumulant, also called Binder ratios Binder (1981), are given by

[TABLE]

These quantities together with the correlation length in unit of the linear system size share the property of being invariant under RG transformation, i.e., they tend towards well-defined fixed point values as the system approaches the critical point. As a first step in our analysis we determine these values for $R\equiv\{\xi/L,U_{4},U_{6}\}$ as well as the eigenvalue $\omega$ of the first irrelevant RG operator using the quotient method Ballesteros et al. (1996, 1997, 1998). This finite-size scaling (FSS) technique allows for great statistical accuracy and does not require a precise knowledge of the critical temperature. Specifically, the scaling functions $R$ are evaluated at crossing points of $\xi/L$ curves, where they are expected to scale as

[TABLE]

neglecting additional sub-leading correction terms. Here, $\omega$ is the leading correction-to-scaling exponent, $R^{*}$ denotes the value of the RG-invariant scaling functions at the fixed point and the amplitudes $a_{R}$ of the leading corrections depend on the respective scaling functions $R$. Furthermore, in Eq. (7) the quotients are defined as $Q_{O}=O(sL,T)/O(L,T)$, which means that the observable $O$ – in our case $\xi$ – is measured at a temperature for which the correlation length in units of the linear lattice size is the same for the pair $(L,sL)$. We fix $s=2$.

Using histogram reweighting techniques Ferrenberg and Swendsen (1988, 1989) we determine the crossing points of $\xi/L$ precisely. The reweighting procedure is performed for every disorder replica individually and the curves are averaged afterwards. Up to $10^{5}$ disorder realizations are used for the smallest lattices and at least 4000 replicas for the largest ones. Every replica is initially prepared in a cold spin configuration and is thermalized using 1000 Elementary Monte-Carlo step (EMCS). We checked for a proper thermalization by also performing simulations starting from a hot configuration, which gives identical results within numerical precision. In our update procedure, one EMCS consists of a full Metropolis lattice sweep and several single-cluster updates. Since the average cluster size $\langle|C|\rangle$ in a $d$-dimensional system at criticality scales as $L^{d-\gamma/\nu}$, we increase the number of cluster updates with lattice size according to $L^{0.25}$ in order to keep the fraction of flipped spins approximately independent of the lattice size Janke and Villanova (2002). A detailed list of replica configurations and cluster steps can be found in Tab. 2. The simulations reported in this section took about half a million CPU hours on an Intel Xeon E5-2697 v3 processor.

In order to obtain the fixed point phenomenological couplings and the correction exponent, we perform simultaneous fits of Eq. (7) for all three couplings $\{\xi/L,U_{4},U_{6}\}$ with joint $\omega$ and for different $L_{\text{min}}$, i.e., discarding the smallest lattices in the fits. As the quotients $Q_{O}$ are naturally correlated in pairs $(L,2L)$, we implement a fitting procedure that optimizes a generalized $\chi^{2}$, including the full self-covariance information, as proposed in Ref. Gordillo-Guerrero and Ruiz-Lorenzo (2007). Uncertainties for the fit parameters are obtained by bootstrap resampling methods Efron and Tibshirani (1994)222We construct 250 bootstrap samples of the full data set by averaging the observables over $N_{r}(L)$ randomly drawn disorder replicas rather than performing a simple average where every replica is considered exactly once. For every of those bootstrap samples the fits are performed, resulting in 250 estimates for the fit variables. Averages and standard deviations of these estimates are reported as the final results.. The results are shown in Tab. 1. As $L_{\text{min}}$ is increased, the fit results show slight systematic trends, caused by higher-order corrections. Above $L_{\text{min}}\approx 24$ the values saturate and the $\chi^{2}$/d.o.f of the fit does not improve further. We therefore use, as our final estimates, the averages for $L_{\text{min}}=24,32,40$ and adopt a rather conservative uncertainty which includes the fluctuations among the single estimates as well as their individual uncertainties. This yields, as our final results

[TABLE]

and

[TABLE]

Comparing our estimates for the critical couplings with reference values of the Ising model on a regular square lattice, which are known exactly, up to small uncertainties from numerical integration, $(\xi/L)^{*}=0.90505\dotso$, $U_{4}^{*}=0.61069\dotso$, $U_{6}^{*}=1.45565\dotso$ Kamieniarz and Blote (1993); Salas and Sokal (2000) we find that even though these quantities are only considered universal in a limited sense (they weakly depend on certain geometrical characteristics of the system Selke and Shchur (2005, 2009); Malakis et al. (2014)) they compare remarkably well, giving a first indication that the Ising model on a VG stays in the universality class of the clean model. Also our result for the correction exponent is noticeably smaller than the reference value on a square lattice, $\omega=2$ Blöte and den Nijs (1988), though not particularly small in absolute numbers, which explains why corrections to scaling turn out to be relatively weak in the scaling collapses described below.

Essential for computing scaling collapses is a precise knowledge of the location of the critical point, which depends on the details of the lattice structure and is therefore in general not known in advance. In the framework of quotient-FSS the critical temperature can be obtained using infinite-volume extrapolations, as the crossing points are expected to scale according to Binder (1981)

[TABLE]

where higher-order terms have been neglected as in Eq. (7), and we adopt the clean exponent $\nu=1$. We perform four series of fits, where in the first three the correction exponent is fixed to our previous estimate (8), plus and minus its uncertainty, i.e., $\omega=1.38,1.54,1.70$. In the last series of fits $\omega$ is treated as a free parameter. The results are displayed in Fig. 2 and listed in detail in the Supplementary Material, including $\chi^{2}$/d.o.f values. It can be seen that for $L_{\text{min}}\gtrsim 40$ all four fits are compatible within their error bars. As our final estimate we take the average of the fixed-$\omega$ fits for $L_{\text{min}}=64$, obtaining

[TABLE]

In the next step of the analysis we simulate the Ising model on Voronoi graphs of size $L=24,32,\dotsc,384$ for several temperatures in the vicinity of the critical point, using at least $1000$ disorder replicas for every lattice size and temperature. Similar to the precision simulations directly at criticality reported above, we start from cold configurations and perform 2500 measurements after a thermalization time of 500 EMCS. The number of cluster updates per EMCS is reduced by about a factor of five with respect to the values reported in Tab. 2.

Finite-size scaling theory predicts that magnetization, susceptibility and Binder ratio scale according to

[TABLE]

where $\beta$, $\gamma$ and $\nu$ are critical exponents of the model and the universal scaling functions $f$ have the argument

[TABLE]

These equations describe the scaling behavior to first order. Corrections of higher order are expected to become irrelevant for large system sizes. In Fig. 3 we show the scaling collapse plots, fixing all critical exponents to their exactly known values ($\nu=1$, $\beta=1/8$, $\gamma=7/4$) and $T_{c}$ to the estimate (11). As can be seen, a flawless collapses for all three scaling functions is obtained even for small lattices, which shows that the Ising model on a 2D random Voronoi graphs belongs to the universality class of the clean 2D Ising model.

IV Directed Percolation

The most prominent family of non-equilibrium phase transitions arguably is the so-called directed percolation (DP) universality class. According to the Janssen-Grassberger DP conjecture, any system featuring a fluctuating phase and a unique absorbing state with scalar order parameter, no additional symmetries and only local interactions falls into this class Janssen (1981); Grassberger (1982). On a lattice, the DP universality class is realized by certain reaction-diffusion schemes, such as the contact process (CP). In this model every site can be in either of two states, active or inactive. In the epidemic language, often used in this context, one refers to infected and recovered states. The system is evolved as a Markov process, where in each time step, one random active site is picked. With probability $p$ it infects a random neighbor whereas with probability $1-p$ the particle spontaneously recovers and is removed from the active cluster. Time is then incremented by $1/N_{a}$, where $N_{a}$ denotes the size of the active cluster before the update. As soon as the system enters the so-called absorbing state where every site is inactive, the dynamics terminates.

In order to show that the CP on a random VG shows clean universal behavior we conduct numerical simulations as described in Ref. Schrauth et al. (2018a), in total amounting to approximately 300 000 CPU hours on an Intel Xeon E5-2697 v3 processor. As a first step, from seed simulations we determine the critical point, by rescaling the cluster size and survival probability according to their expected power law behavior, $N_{a}(t)\sim t^{\theta}$ and $P(t)\sim t^{-\delta}$, where $\theta$ and $\delta$ denote critical exponents. In total, we use $10^{5}$ independent disorder realizations of linear size $L=2048$ with periodic boundary conditions and performed $10^{4}$ seed runs on each of them. Using reference values from Ref. Dickman (1999), $\theta=0.2293(4)$ and $\delta=0.4523(10)$, we obtain the critical probability

[TABLE]

by determining the asymptotically constant curve, as shown in Fig. 4. The uncertainty is determined from curves which noticeably bend away from horizontal behavior. Note that using larger lattices and longer simulation times would not significantly increase the precision of this estimate, as the analysis is limited by the uncertainties of the reference values.

Once the critical probability is known, we perform decay simulations starting from a initially fully occupied lattice for different system sizes precisely at criticality and monitor the density $\rho(t)$ of active sites until the system reaches the absorbing state. This allows us to obtain the exponents $\delta$ and $z$ via data collapses according to

[TABLE]

where $\tilde{\rho}$ denotes a universal scaling function. In a second set of simulations we perform decay simulations for lattices of fixed size $L=1024$ in the vicinity of the critical point, which gives us the exponents $\delta$ and $\nu_{\parallel}=z\nu$ by means of the relation

[TABLE]

where $\Delta=p-p_{c}$ is the distance from criticality and $\hat{\rho}$ a scaling function. Both scaling collapses, which turn out flawless, are shown in Fig. 5, where we used the reference values, $\delta=0.4523(10)$, $z=1.7674(6)$ and $\nu_{\parallel}=1.292(4)$ from Ref. Dickman (1999). In the top panel, all curves are averages over 1400 disorder realizations with 5 runs per realization, whereas in the bottom panel we used 250 disorder replicas and 5 runs per realization. The insets show the respective non-rescaled density as a function of time. This confirms the critical exponent values used and provides compelling evidence that CP on the random VG belongs to the clean DP universality class.

V Conclusion

Critical phenomena on Voronoi graphs have, to the best of our knowledge, not yet been investigated, focusing instead on lattices such as its dual, the Delaunay triangulation. In order to correct for this omission, we conducted large-scale numerical simulations of the classical Ising model and contact process on two-dimensional VGs constructed from randomly distributed sites. We establish reference values for the critical points of both models and, for the Ising model, also obtain the first correction-to-scaling exponent. Furthermore, using finite-size-scaling techniques we show that both systems display clean universal exponents at criticality, i.e., we reveal that the VG disorder is – in the sense of the RG – an irrelevant perturbation to their phase transitions. Although we only analyzed two particular models, this result has implications for other classes of transitions as well. From the RG perspective, the correlation length exponent $\nu$ is directly related to the relevance of quenched disorder Harris (1974); Weinrib and Halperin (1983); Kinzel (1985); Harris (2016); Barghathi and Vojta (2014). For instance, the phase transition of regular (isotropic) percolation can be predicted to also remain unchanged on a Voronoi graph, since its exponent, $\nu=4/3$, is larger than that of both models considered. Moreover, our results are especially relevant for the search for a general disorder relevance criterion. The Voronoi graph has constant coordination number, similar to the constant coordination (CC) lattice we have recently introduced Schrauth et al. (2018b) and refined Schrauth and Portela (2019). Studies of both the Ising and DP phase transition on the CC lattice Schrauth et al. (2018b, 2019) found disorder to be probably marginal in the Ising case and clearly relevant for the DP universality class. The contrast of this result with the irrelevance of the VG disorder shows that the absence of coordination number fluctuations is non-predictive of disorder relevance. Thus, coordination number fluctuations do not play a central role in determining the influence of quenched topological disorder on continuous phase transitions and a different direction, such as considering a measure of connectivity Schrauth et al. (2019), should be explored in the search for a general relevance criterion.

Acknowledgements.

We thank H. Hinrichsen, F. Goth and F. P. Toldin for helpful discussions. M.S. thanks the Studienstiftung des deutschen Volkes for financial support. This work is part of the DFG research project Hi 744/9-1. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de).