{'fraction_non_alphanumeric': 0.0373458278784538, 'fraction_numerical': 0.01319506649211052, 'mean_word_length': 4.476792452830189, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability of line-node semimetals with strong Coulomb interactions and properties of the symmetry-broken state

We employ diagrammatic Monte Carlo simulations to establish criteria for the stability of line-node semimetals in the presence of Coulomb interactions. Our results indicate a phase transition to a chiral insulating state that occurs at a finite interaction threshold which we determine. We also compute the Landau levels for out-of-plane and in-plane magnetic fields in the symmetric and symmetry-broken phases. We find that the magnetic field couples to the chiral order parameter, implying that this degree of freedom can be manipulated in situ in experiments. Finally, we check the existence of edge states in the symmetry-broken phase. On the system's boundary, we note that the metallic"drum-head"states that exist in the symmetric phase are gapped out. However, the symmetry-broken phase permits topological defects in the macroscopic order parameter in the form of domain walls, which host metallic"interface states."These consist of line-like gap-closings that occur on the two-dimensional interfaces.

Topological semimetals exhibit a band structure that is gapped everywhere except at a few nodal points, where the bands meet. The topological nature of these protects them against perturbations and also gives rise to unique phenomena like the chiral anomaly-where nodal points act as sources and sinks of a spontaneous current-and Fermi arcs, which are direct manifestations of the bulk topology [1]. Applications of topological semimetals are thus far mainly as building blocks of information technology [2][3][4]. Several of these materials exhibit a large magnetoresistance effect [5] that may be exploited in magnetic field sensors [6], and spintronic devices [7].
In the search for new and technologically useful semimetals, symmetry-protected topological phases have emerged as an important platform that dramatically expands the types of nodal features realized in a material. Typically, this involves a combination of point group symmetries and explicitly broken time-reversal or inversion symmetry. The result is a wide class of band touching points that carry more than unit topological charge [8,9] and may also involve multiple bands [10,11] or have to be classified as line-nodes [12][13][14][15][16][17][18].
While the explicit reliance on symmetry significantly widens the scope of topological and semimetallic materials, it also has implications for the role of correlations in these: In tight-binding models of bilayer graphene, the dispersion is quadratic around the nodal points, yet the system develops a nematic instability at an infinitesimal interaction [19]. For multiple-charge Weyl nodes and line-node semimetals with contact interaction, renormalization group theory indicates instabilities at a finite threshold [20,21]. This is also true for single-layer graphene, where actual material parameters are situated close to a chiral symmetry-breaking regime [22]. Thus, semimetallic phases that depend explicitly on symmetry may be susceptible to correlation effects that destroy the underlying symmetry.
Reliably predicting the parameter regimes where symmetry-protected topological phases remain stable faces several delicate problems. The absence of screening means that interactions are effectively long-ranged, leading to infrared divergencies. Furthermore, competition between correlation effects originating from different length scales is very likely. For example, in graphene, the long-range part of the interaction drives the system towards an asymptotically free Dirac liquid with divergent Fermi velocity as T → 0, without causing any instabilities [23]. Thus, the renormalization of the dispersion occurring primarily at the infrared end increases the effective kinetic energy relative to the short-range part of the interaction, which is believed to drive the phase transition. Therefore, an accurate solution to this class of problems requires that all length scales are treated on an equal footing.
A type of symmetry-protected topological phase that has attracted considerable interest is the nodal-line semimetal. Predictions of this state has been made in TlTaSe 2 [15], CaAgP [16] and Ca 3 P 2 [17] and ZrSiS [24] based on reflection symmetry. In PbTaSe 2 it has also been confirmed by angle resolved photo emission spectroscopy [14]. Recently, the observations of strongly renormalized transport properties and Fermi velocity-as compared to DFT calculations-in ZrSiS was interpreted as an indication of a strongly correlated linenode semimetal [25].
In this work, we employ diagrammatic Monte Carlo simulations [26] to establish quantitative criteria for the stability of line-node semimetals in the presence of long-range interactions and also characterize the symmetry-broken phase which occurs for sufficiently strong coupling. We find evidence for a chiral insulator that supports metallic interfaces on domain walls that interpolate between different signs of the order parameter and can be manipulated in situ via an external field.

MODEL
We consider the case of a single nodal line with a bare Fermi velocity of v 0 f running along the z-axis Here, λ is a fictitious screening length introduced to regularize the series, and we are thus principally interested in the limit λ → ∞. We consider a cylindrical domain given by where Λ is the ultra violet cutoff. Because of a scale invariance associated with the linear dispersion, the only relevant length scale in the low-temperature limit is the ratio of the UV cutoff and the inverse screening length Λ/λ −1 . To see this, we may choose a temperature and energy scale where the temperature is unity by rewriting the partition function z(β, H) = z(1, βH). This gives a bare Greens function where we have exploited the linearity of H 0 in k. Diagrammatic corrections to the Greens function take the form where N is the expansion order. If we introduce a change of scale k = βk we obtain For a screened Coulomb interaction in D = 3 we find The UV cutoff changes scales as Λ → βΛ, giving which is characterized by the ratio Λ/λ −1 in the limit β → ∞.
In the perturbative regime, the nodal line (1) is protected by a symmetry due to being odd under an orthonormal map k → −k. The implication of this symmetry is that on the k zaxis, the Greens function must have a pole at zero energy as long as the series expansion remains convergent [27]. Correspondingly, destroying the semimetallic phase requires breaking this symmetry. In a diagrammatic framework, this phase transition can be identified via a divergent susceptibility with respect to a symmetry-breaking perturbation.

CONTACT INTERACTION
For contact interaction, the self-consistent Fock theory can be solved analytically due to translation invariance in momentum space. Specifically, the self-energy satisfies the relation Here, it should be noted that at the level of Fock theory, Σ is independent of frequency, and thus de facto takes the form of a correction to the effective dispersion. Furthermore, contact interaction does not renormalize the Fermi velocity since H 0 (k) is an odd function. Since the self energy is translation invariant, it must therefore take the form Σ(k) = ∆σ z . Inserting this self-energy in (9) and summing over frequency, we obtain where the integral of H 0 over k vanishes. Thus, Eq. (10) provides a self-consistent equation for ∆ as a function of the coupling strength, whose solutions will provide the gap parameter in this regime. Solutions for which ∆ is finite correspond to a symmetry-broken state, while the symmetric phase is characterized by a vanishing gap. In the low-temperature limit, and for a cylindrical domain (3) with Λ = 1, the integral (10) provides an algebraic expression for the gap of the form where we have introduced η = α(2π) −2 . Equation (11) predicts a critical coupling strength η c /v 0 f = 2, see Fig. 1. Above this threshold, the gap is given by The onset of chiral a phase at a finite interaction strength is consistent with results from renormalization group theory for contact interaction [21].

SIMULATIONS
To check the stability of the semimetallic phase for longrange interactions, we employ diagrammatic Monte Carlo, which is a computational protocol based on the stochastic sampling of the diagrammatic expansion. Within this framework, the space of connected diagrams for the self-energy is sampled through a metropolis type random walk [26,[28][29][30]. The Greens function is then obtained via Dyson's equation [31] Here, we use a sampling protocol based on the worm algorithm as described in [32]. We employ a bold scheme where the expansion is conducted in dressed Greens functions while retaining only skeleton graphs. Thus, at order N = 1, the solution corresponds to self-consistent Fock theory. We do not employ bold interactions lines since this is expected to have little advantage for a semimetallic system. Following the scaling relation (8) we can without loss of generality set Λ = 1. This gives a volume of the momentum space of 2π. We then rewrite the integral over k as which defines a rescaled interaction parameter η and a set of units where the integral over momenta is of measure unity. We parameterize the temperature and scale in terms of a variable γ so that The limit γ → ∞ thus corresponds to zero temperature and a divergent ratio Λ/λ −1 . An observable that is convergent in this limit should correspondingly be a function of η/v 0 f . To obtain a self-consistent solution for the model (1-2), we consider a starting guess for the frequency-independent selfenergy Σ 0 (k), which in turn provides a corresponding Greens function G 0 (ω, k). A stochastic summation of the expansion in V gives a new self-energy Σ 1 (ω, k) which is subsequently used. This scheme is repeated until relevant observables have converged. Near the phase transition, this typically requires several hundred iterations. We have used two starting configurations for the self-energy, featuring extremely small or relatively large symmetry-breaking terms, respectively. For most parameter regimes, these result in identical solutions. However, at low temperatures and for a coupling strength that is slightly larger than the critical coupling, we observe a family of very fragile meta-stable symmetric solutions that likely result from competition between the symmetric and antisymmetric parts of the self-energy.
To track the onset of a symmetry-broken phase, we define the chiral order parameter as follows. First, we note that the frequency-independent part of the self-energy can be written The chiral symmetry-breaking is generated by the zcomponent, prompting us to define an order parameter of the form Since we consider a straight nodal line, we assume a solution that is translation invariant in the z−direction. Furthermore, we assume a symmetry of the self-energy where R z φ represents a rotation around the z−axis by φ. The results from the diagrammatic Monte Carlo simulations are summarized in Fig. 2. For most parameter ranges, the order parameter scales approximately as O ∼ (η/v 0 f ) 2 , prompting us to plot the square root. The solutions correspond to different values of γ, which controls the model parameters according to (15). As we progressively decrease the temperature and increase the screening length, the order parameter saturates to a single line which depends only on η/v 0 f , indicating that the chiral order exhibits a well-defined IR limit at zero temperature. The presented data corresponds to a first and second-order expansion.
In Fig. 3 we see the critical coupling strength as a function of γ, extracted from the data presented in Fig. 2. The critical point saturates to η c /v 0 f = 0.45 ± 0.1. The correction from first to second order falls within the error bars, indicating that this problem is well captured by self-consistent Fock theory. This is consistent with previous applications of diagrammatic techniques to semimetallic systems: In Weyl semimetals, the correction to the Greens function is almost entirely contained in the frequency-independent part of the self-energy, leading to the emergence of virtually free fermions [33]. In graphene, at least the long-range part of the interaction drives the system towards an asymptotically free Dirac liquid [23], while for short-range interactions, the convergence of the series has been demonstrated analytically up to a finite threshold [34]. Most likely, this results from the exponential suppression of diagram topologies that involve excitation of the background in semimetals.

LANDAU LEVELS AND MAGNETIC RESPONSE
To compute the Landau levels arising when the line-node is placed in a magnetic field, we consider a dispersion of the form (1) and take the Fermi velocity to be unity. This gives with energy bands Thus, the nodal line runs along the z-axis (0, 0, k z ).  γ=2  γ=3  γ=4  γ=5  γ=6  γ=7   γ=0  γ=1  γ=2  γ=3  γ=4  γ=5  γ=6  γ=7 sents a rotation around the z axis by . rom the diagrammatic Monte Carlo simularized in Fig. 1. Setting the COUPLING TO MAGNETIC FIELD

SURFACE STATES
Line nodes: DiagMC-Method Strong interactions: [28][29][30][31] endent part of the self energy can be written metry-breaking is generated by the zmpting us to define an order parameter of the er an infinitely long nodal line, we assume a translation invariant in the z direction. Fursume a symmetry of the self energy sents a rotation around the z axis by . rom the diagrammatic Monte Carlo simularized in Fig. 1. Setting the . Note that we display the square root of the order parameter since it is approximately quadratic in ⌘ over most parameter ranges. Here, the screening parameter is 1 = 10 2 ⇥ 2 , while v 0 f = 10 2 ⇥ 2 so that the solutions correspond to progressively lower temperatures and longer screening lengths. For larger values of the solutions collapse onto a single line indicating that the chiral order remains convergent in the infrared limit.

COUPLING TO MAGNETIC FIELD SURFACE STATES
Line nodes: DiagMC-Method Strong interactions: [28][29][30][31] 3 is is expected to have m. e can without loss of me of the momentum al over k as rameter ⌘ and a set of is of measure unity. d scale in terms of a zero temperature and e that is convergent in unction of ⌘/v 0 f . or the model (1-2), we ency-independent self corresponding Greens ation of the expansion which is subsequently vant observables have this typically requires ed two starting configxtremely small or relrespectively. For most tical solutions. Howupling strength that is g, we observe a famic solutions that likely mmetric and antisymly (see also discussion oken phase, we define First, we note that the energy can be written . enerated by the zorder parameter of the ).
rameter ⌘ and a set of is of measure unity. d scale in terms of a zero temperature and e that is convergent in unction of ⌘/v 0 f . or the model (1-2), we ency-independent self corresponding Greens ation of the expansion which is subsequently vant observables have this typically requires ed two starting configxtremely small or relrespectively. For most tical solutions. Howupling strength that is g, we observe a famic solutions that likely mmetric and antisym- . Note that we display the square root of the order parameter since it is approximately quadratic in η over most parameter ranges. Here, the screening parameter is λ −1 = 10 −2 × 2 −γ , while βv 0 f = 10 2 × 2 γ so that the solutions correspond to progressively lower temperatures and longer screening lengths. For larger values of γ, the solutions collapse onto a single line indicating that the chiral order remains convergent in the infrared limit. To couple the system to an external magnetic field, we introduce the displacement of the momentum k by the vector potential A: First, we consider a magnetic field along the z-direction, i.e., B = Bẑ. Working in the axial gauge, the vector potential reads A = (−By/2, Bx/2, 0), so that we can define the ladder operators in terms of the new momenta escaled interaction parameter ⌘ and a set of ntegral over momenta is of measure unity. ize the temperature and scale in terms of a t thus corresponds to zero temperature and ⇤/ 1 . An observable that is convergent in correspondingly be a function of ⌘/v 0 f . f-consistent solution for the model (1-2), we g guess for the frequency-independent self hich in turn provides a corresponding Greens ). A stochastic summation of the expansion self energy ⌃ 1 (!, k) which is subsequently e is repeated until relevant observables have the phase transition, this typically requires terations. We have used two starting configelf energy, featuring extremely small or relmetry-breaking terms respectively. For most es these result in identical solutions. Howperatures and for a coupling strength that is an the critical coupling, we observe a fame meta-stable symmetric solutions that likely etition between the symmetric and antisyme self energy respectively (see also discussion nset of a symmetry-broken phase, we define arameter as follows. First, we note that the ndent part of the self energy can be written metry-breaking is generated by the zpting us to define an order parameter of the r an infinitely long nodal line, we assume a Plotting the square root of the orde As we increase the screening leng perature, the order parameter saturat lapse. Here, = 10 2 ⇥2 , so that the ( = 7) corresponds to = 12800, w ture in units of the bare Fermi velocit

SURFACE STAT
Line nodes: DiagMC-Method Strong interactio Emergent free femrion [32] Symmetry: [25] Graphene: [21]. Bilayer graphene: [17] scale in terms of a = 2 ˜ 1 . (21) ro temperature and hat is convergent in ction of ⌘/v 0 f . the model (1-2), we cy-independent self rresponding Greens on of the expansion ich is subsequently nt observables have s typically requires two starting configemely small or relpectively. For most al solutions. Howling strength that is we observe a famsolutions that likely etric and antisym-(see also discussion n phase, we define st, we note that the rgy can be written (22) erated by the zer parameter of the (23) l line, we assume a z direction. Furlf energy z axis by . onte Carlo simulat parameter ranges,  Line nodes: DiagMC-Method Strong interactions: [28][29][30][31] Emergent free femrion [32] Symmetry: [25] Graphene: [21]. Bilayer graphene: [17] The Landau levels can be easily found from the eigenequation HΦ = EΦ. For Φ = (|φ 1 , |φ 2 ) T , we obtain the two equations √ 2Ba|φ 2 = E|φ 1 , and, by inserting the first equation into the second, we arrive at which describes a harmonic oscillator with where c is a normalization constant. Then, for |φ 1 we have Therefore, for a magnetic field in the z direction, the Landau levels are given by with eigenstates The eigenstates appear as a spectrum of bands that are flat in all directions. The gap between these is controlled by the external field, as shown in Fig. 4. In contrast to Weyl semimetals, line nodes do not give rise to a chiral anomaly [35] in the presence of a magnetic field. Instead, a single band remains at the Fermi level because the dispersion is independent of k z .
For an in-plane field of the form B = Bx, we may choose a vector potential of the form A = (0, −Bz/2, By/2). The ladder operators can then be constructed as Expressed in this language, the Hamiltonian (19) takes the form From the eigenequation HΦ = EΦ, we obtain for |φ 1 and |φ 2 As before, we can take the first equation and plug it into the second. By using the commutator [a, a † ] = 1, we finally arrive at Solving the equation (35) is complicated by the presence of terms of the form ∼ a 2 and ∼ (a † ) 2 , which render it anharmonic so that standard recipes for extracting the Landau levels are not applicable. For this reason, we adopt the Bargmann representation [36][37][38], which has been widely used in this scenario. Notably, this technique was applied to an anharmonic oscillator with a quartic potential [39] and the twomode squeeze harmonic oscillator and the kth-order harmonic generation [40]. In this representation, the ladder operators are related to a complex variable z according to whilst the wave function is a holomorphic function of z only, namely, Expressed in this formalism, Eq. (35) takes the form where we have used the notation ϕ 2 = dϕ2 dz and ϕ 2 = d 2 ϕ2 dz 2 . The different quantum states correspond to solutions of this equation for corresponding quantum numbers, such as polynomials of degree n. In principle, it is possible to extract a solution in the form of a power series in z, though this turns out to be highly inefficient. Thus, we instead introduce a reparameterization of ϕ 2 (z) given by The differential equation then takes the form The trivial solution consisting of a combination of two exponentials does not correspond to the Landau levels, and the energy E still appears as an arbitrary constant. To extract the nontrivial solutions, we need to introduce a change of variables of the form which finally transforms the equation for ψ 2 into Now, we may find solutions of this equation as a polynomial in ρ of degree n by considering Plugging this into the equation we obtain an expression in terms of the coefficients of the expansion f i Figure 5. Landau levels for the line-node semimetal in the presence of an external field B =x (B = 1) as given by Eq. (45). In contrast to the case of a magnetic field along theẑ direction in Fig. 4, the levels do exhibit a momentum dependence. The two energy levels that are closest to the Fermi surface cross in the origin.
Since our assumption of a polynomial of degree n implies f n = 0, it trivially follows that n 2 + ω 2 = 0, giving Landau levels with an energy together with f i = 0, ∀ i = n. The function ψ 2 then reads where c is a normalization constant. Returning to the original wave function component ϕ 2 we get The ϕ 1 component may be written as In contrast to the Weyl semimetals, the states corresponding to n = 0 can be treated on an equal footing with the remaining levels. This quantum number gives two different states related by E 1 = −E 2 that cross the Fermi level at k x = 0, as depicted in Fig. 5. Finally, we will review these results once a symmetry breaking term of the form ∆ σ z is included in the Hamiltonian. This gives (49) Hence, the energy bands are now given by so that the system is an insulator. The symmetry-breaking perturbation introduced in the system has now gapped out the nodal line, suggesting that the Landau levels will form away from the Fermi level irrespectively of the orientation of the magnetic field. For the case of an external field B = Bẑ with ladder operators defined according to Eq. (22), the Hamiltonian reads Following the same approach as above and considering the eigenequation HΦ = EΦ, we obtain a harmonic oscillatorlike equation for |φ 2 with the solution where c is a normalization. Comparing with Eq. (29), we see that the effect of the perturbation ∆ is to introduce a displacement of the Landau levels. On the other hand, for |φ 1 we have This scenario is slightly different from the unperturbed case, and the value n = 0 needs to be considered separately since |φ 1 = 0 and |φ 2 = |0 . As a result, the energy for n = 0 is given by E 0 = −∆. In conclusion, we thus find E n,± = ± 2Bn + ∆ 2 , ∀ n = 0.
(54) Fig. 6 shows the Landau levels (54) for B = 1 as a function of the symmetry breaking parameter ∆. Besides the aforementioned displacement of the energy levels for n = 0, the unperturbed E 0 = 0 state is shifted relative to the Fermi level depending on ∆ in a similar manner to the k z -dependence of the Landau levels appearing in Weyl semimetals. Thus, for a symmetry-broken state, the magnetic field lifts the degeneracy between the two chiralities, implying that this degree of freedom can be manipulated by an external field in experiments.
Finally, we consider the case of an in-plane magnetic field B = Bx in the symmetry-broken phase. Expressed in the ladder operators defined in Eq. (31), the Hamiltonian takes the form As before, the eigenvalue equation gives rise to an anharmonic problem, meaning that we have to rely on the Bargmann representation. For |φ 2 the solution is given by Eq. (35), except for a shift in energy given by E 2 → E 2 − ∆ 2 . Hence, we can use the same protocol as above for the in-plane field with the ω parameter accordingly modified to  Figure 6. Landau levels for a symmetry-broken line-node semimetal in the presence of an external field B =ẑ (B = 1) as a function of the symmetry breaking parameter ∆. The applied field explicitly breaks the symmetry between the chiralities ∆ and −∆ respectively, indicating that the chiral order can be manipulated via a magnetic field.
Landau levels for the symmetry-broken line-node semimetal in the presence of an external field B =x (B = 1).
As the system becomes gapped, the crossing of the Landau levels shown in Fig. 5 is lifted.
The Landau levels are given by with the wave function components where c is a normalization constant.
In the symmetry-broken phase, the system becomes gapped, with Landau levels situated at finite energies, as seen in Fig. 7. The levels closest to the Fermi surface attain a gap of |∆| which is thus independent of the chirality.

INTERFACE STATES
In the symmetric phase, the line-node semimetals exhibit metallic drum-head surface states [15], which are generalizations of the Fermi arcs that occur in the Weyl semimetals [1]. These states are stabilized by a combination of topology and symmetry in the sense that the states are topologically protected in a subspace generated by the symmetry. Once this symmetry is spontaneously broken, the lines are gapped out, and the edge states are no longer protected.
However, the symmetry-broken phase permits domain walls that interpolate between regions of different chirality, on which the symmetry-breaking term ∆ changes sign. This opens the possibility for metallic interface states that are bound to these topological defects. To model this scenario, we consider a domain wall described by where ∆ + > 0 and ∆ − < 0 are constants. To identify the interface states, we apply an analytical approach based on trial functions that has been applied to Fermi arcs within Weyl semimetals in semi-infinite systems [41,42]. Since ∆(y) is translation invariant in thex andẑ directions but not alongŷ, it follows that k x and k z are good quantum numbers while k y is not. Hence, we conduct the substitution k y → −i∂ y , which transforms the perturbed Hamiltonian into H(k x , −i∂ y , k z , y) = k x σ x − i∂ y σ y + ∆(y)σ z . (61) Next, we introduce a trial wave function of the form Ψ(x, y, z) = ψ λ |x, z = ψ 1 ψ 2 e λy |x, z .
Therefore, our problem is reduced to the eigenequation with a continuity condition at Ψ(y = 0). The secular equation, det|H(k x , −i∂ y , k z , y) − E| = 0, may be used to find the possible values of λ, namely, Requiring the wave function to vanish at y → ±∞, we need to separate the two regions of different chirality into Ψ + for y > 0 and Ψ − for y < 0. Then, we obtain Ψ + = c + ψ λ+ |x, z = c + ψ + 1 ψ + 2 e λ+y |x, z ,  One may note that, after the introduction of ∆(j), the system is gapped in every layer but the one corresponding to j 0 , where ∆(j 0 ) = 0. It should therefore be expected that the interface states fall off as where c 0 is a normalization constant. The red curve in Fig.  9 shows a fit of the form (86) with c 0 as a free parameter, revealing that there is an excellent agreement. The result for an even number of layers was found to be indistinguishable from the case of an odd number implying that the metallic interface states do not depend on details of the domain wall.

SUMMARY
In conclusion, we have examined the stability of line-node semimetals in the presence of Coulomb interactions and found a chiral instability occurring at a finite interaction strength. The chiral order parameter exhibits a well-defined behavior in the limit of an infinite screening length despite the presence of infrared divergencies in this problem. By computing the Landau levels, we observe that an out-of-plane magnetic field couples to the chiral order parameter, implying that this degree of freedom can be controlled in experiments. While the drum-head edge states associated with line-node semimetals vanish in the chiral phase, we observe metallic interface-states in this regime, which exist on domain walls interpolating between regions of different chirality. These domain walls could conceivably be trapped on a sample with a concave geometry in experiments.
This work was supported by the Swedish Research Council (VR) through grant 2018-03882 and Stiftelsen Olle Engkvist via grant 204-0185. Computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre in Linköping, Sweden. J. C. would like to thank Lars Fritz for important input and discussions.