Finite-density phase diagram of a (1+1)-d non-abelian lattice gauge theory with tensor networks

We investigate the finite-density phase diagram of a non-abelian SU(2) lattice gauge theory in (1+1)-dimensions using tensor network methods. We numerically characterise the phase diagram as a function of the matter filling and of the matter-field coupling, identifying different phases, some of them appearing only at finite densities. For weak matter-field coupling we find a meson BCS liquid phase, which is confirmed by second-order analytical perturbation theory. At unit filling and for strong coupling, the system undergoes a phase transition to a charge density wave of single-site (spin-0) mesons via spontaneous chiral symmetry breaking. At finite densities, the chiral symmetry is restored almost everywhere, and the meson BCS liquid becomes a simple liquid at strong couplings, with the exception of filling two-thirds, where a charge density wave of mesons spreading over neighbouring sites appears. Finally, we identify two tri-critical points between the chiral and the two liquid phases which are compatible with a $SU(2)_2$ Wess-Zumino-Novikov-Witten model. Here we do not perform the continuum limit but we explicitly address the global $U(1)$ charge conservation symmetry.


Introduction
Lattice gauge theories lie at the heart of the description of nature. They provide the theoretical ground to formulate theories ranging from the Standard Model in high energy physics [1,2] to effective condensed matter models [3,4], while being able to to capture superconducting and topological states of matter. Despite the tremendous success of such theories, our capability of non-perturbative analysis -which mostly relies on numerical techniques -has been severely limited to zero-density scenarios. Indeed, elucidating the properties of gauge theories at finite density is notoriously challenging [5,6] as, differently from their zero-density counterparts where Monte Carlo simulations have been tremendously successful [7,8,9], importance sampling is invalidated due to the so called sign problem. Similarly, quantum simulators of lattice gauge theories [10,11,12,13,14,15,16], despite being extremely promising platforms, are still in the early stage of their development [17,18].
In this work, we investigate the zerotemperature, finite-chemical-potential phase diagram of a SU(2) lattice gauge theory in (1+1)dimensions using Tensor Network (TN) methods [19,20,21,22,23]. The finite-density phase diagrams of non-Abelian lattice gauge theories have been suggested to host paradigmatic sce-narios for strongly correlated quantum matter. In condensed matter settings, SU(2) gauge theories naturally emerge as a low-energy description of the Hubbard model [3,24], where superconductivity at finite doping represents the equivalent of meson formation at finite chemical potential. In high-energy physics, the large density regime of quantum chromodynamics (QCD) has attracted a lot of attention especially in view of possible exotic superconducting phases which may appear in regimes close to the interior of neutron stars [25,26]. Here we focus on (1+1)-dimensional non-Abelian models, and specifically on the SU(2) Yang-Mills theory cast as a quantum link model [27,28,29,30], with compact gauge fields coupled to Kogut-Susskind fermions [31]. Firstly, we show how these models support Bardeen-Cooper-Schrieffer (BCS) liquid states of matter pairs (mesons), in a wide range of Hamiltonian parameters, and that the stability of this phase is preserved while increasing the particle density. Since there are no bare interactions between our fermionic fields, this 'colorlike-superfluidity' is solely driven by the interactions mediated by the gauge fields, very much akin to suggested scenarios in SU(2) gauge theories for the Hubbard model [3]. Secondly, we analyse the strong-coupling regime, where we found that chiral symmetry is spontaneously broken, but is immediately restored at finite densities, where a liquid phase with gapless single particle excitations is present. Finally, we track the behaviour of the entanglement entropy at the transition between BCS liquid and the chiralsymmetry-broken phase. Our results signal that the transition is of second order, and is compatible with a SU (2) 2 Wess-Zumino-Novikov-Witten model [32,33,34], which has been put forward as a key scenario for gauge theories with SU(2) gauge invariance [35].
In the last two decades, TNs, and in particular the matrix-product state methods applied here, have been widely employed to numerically tackle strongly correlated, low-dimensional many-body problems [36,37,19,20,38,22]. These wavefunction based methods do not rely on importance sampling, and thus are unaffected by the sign problem, so they can be applied to zero-and finite-chemical potential regimes on equal footing. Moreover, they can easily give information about quantum correlations, such as, e.g., entan-glement entropies, which have recently attracted an increasing attention in the context of gauge theories [39,40,41]. As such, both in terms of regimes of applicability and observables, TNs provide a computational tool which is complementary to Monte Carlo methods. Recently, TN methods have been applied to gauge theories, with investigations ranging from the phase diagrams of Abelian theories in (1+1) and (2+1)-d, to the real time dynamics of U(1) and SU(2) models [42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59]. Here, we focus instead on the effects of finite density of fermions in a non-Abelian lattice gauge theory in the quantum link model formalism, which is very convenient for TN implementations due to its finite-dimensional Hilbert space on each link. Finally, we employ a recently developed TN structure which allows to directly work on the gauge invariant Hilbert space and to achieve the necessary numerical efficiency [45,47].
From a theoretical viewpoint, our study shows how tensor methods allow to tackle questions, such as the presence of superfluid phases at finite density and the restoration of chiral symmetry in non-Abelian gauge theories, which play a fundamental role in considerably more complex theories such as QCD [60,61]. The analysis presented here can be extended to other gauge models, such as SU(3), more complex geometries (ladders, cylinders), while the continuum limit can be studied along the lines of recent TN results [49].

Results
Quantum link formalism for SU(2)− We characterise the phase diagram of a SU(2) lattice gauge theory in (1+1)-d, with spin-1 2 fermionic matter field c [M ] j,s coupled to non-abelian gauge fields U j,j+1;s,s within the quantum link formulation of lattice gauge theories [29,30,62,44,45,15]. One of the main features of quantum link formulation is that the local Hilbert space of a quantum link is finite dimensional and as such it has a finite (fermionic) representation. In this language, gauge (tensor) fields read as compound operators U j,j+1;s,s = c s,s are the entries (row s, column s ) of the Pauli matrix σ (µ) in direction µ ∈ {x, y, z}. Within such quantum link formulation, the group of left (resp. right) transformations applies solely on the rishon mode L (R) of the pair, and this guarantees that the two sets of generators, and thus the two groups, mutually commute [J Notice that, due to the fact that the fermionic representation of the link operators U j,j+1;s,s and J [τ ](µ) j,j+1 are bilinear operators, there is a local conservation of the total number of fermions within a link [29,30,62,44,45,15]. Hence, we fix the total rishon population at every link (L j ,R j+1 ), on which the effective many-body dynamics will indeed depend. From now on, we set n [L] j,↑ + n [L] j,↓ + n [R] j+1,↑ + n [R] j+1,↓ = 2 fermions per link. Within this selection rule, every single quantum link degree of freedom has effective dimension 6, and, as discussed in Ref. [30] and reminded below, the microscopic model has a SU(2) gauge symmetry.
The local or gauge symmetry allows us to impose a SU(2)-equivalent Gauss' law at every site, which restricts the space of 'physical' states |Ψ phys ; in this case, the Gauss' law translates into the fact the total spin around every site j must equal to a spin-0 (see Method's subsection).
Model − The Hamiltonian of model of interest encodes the microscopical dynamics of a Yang-Mills theory on a lattice, and is composed by three terms The link operators play the role of parallel transporters, and appear in the coupling term between gauge fields and matter, as [31] H coupl = t . . . describes the pure-gauge field free energy written in terms of the fermion occupation n j,s , where g 1 = g 0 3/8. If the quantum dynamics is ruled only by H free and H coupl there is an additional, accidental, local conservation of the number of fermions s=↑,↓ n j,s around every site j. As a consequence, the resulting gauge symmetry of the model is SU(2) × U(1) = U(2) [30]. In order to reduce the symmetry of the quantum link model form U(2) to SU(2) we break the additional U(1) gauge symmetry by adding the real part of the determinant of each link matrix to the Hamiltonian, which allows the total population of every site to fluctuate. This extra component of the dynamics moves a fermion pair, which is a total spin-0 due to anti-commutation rules, from the left L j mode to the right R j+1 model of a quantum link, or vice versa. It preserves the total spin on every site since it carries zero spin current, as well as the total population on the quantum link, keeping the symmetries previously discussed intact but the total fermion number over site j and site j + 1. Finally, from now on we set = 1 and = 5, while g 1 ≡ 1 sets the scale of energies. break spontaneously translational invariance in the thermodynamical limit. The first non-critical phase (type A insulator) appears only at matter filling f M = 1, where f M is the total number of matter fermions divided by the number of sites. This phase exhibits CDW order at wave-vector k = π, i.e. it has effective periodicity of 2 composite sites. The second non-critical phase (type B insulator) appears at filling f M = 2/3 instead and it show an effective 3 sites periodicity (wavevector is k = 2π/3). These locally ordered phases arise at large matter-field coupling t: In the following, we report two independent methods to identify the transition point First of all, we pinpoint these phases by investigating the bipartite entanglement quantified by the Von-Neumann entropy S = Tr[ρ log ρ ] of the reduced density matrix ρ = Tr +1..L [|Ψ 0 Ψ 0 |] for sites {1.. } with < L, where |Ψ 0 is the ground state we find, given g 1 , , t and f M . The scaling of bipartite entanglement is a discriminating property between critical and non-critical theories in 1 + 1 dimensions [63]: While for critical systems S S 0 = c 6 log(L sin(π /L)) + c diverges logarithmically with , with c being the central charge of the corresponding conformal field theory (and c a nonuniversal constant), for gapped systems S saturates to a finite value at the thermodynamical limit. Moreover, since we are performing simulations of fermions at fixed filling f M , one has to take into account Fermi oscillations as discussed thoroughly in Refs. [64,65]. In this extended framework, S follows a behaviour according to where the Fermi wave-vector k F enters in the oscillatory contribution. We then use the theoretical expression (5) to fit the simulated bipartite entanglement S and estimate the central charge c, and the Fermi momentum k F of the theory. In Fig. 2 (side panels) we plot several bipartite entanglement profiles for a strong coupling theory at various matter filling fractions 0 < f M ≤ 1.
Fermi oscillations are clearly visible, and allow us to estimate the Fermi wave-vector close to Fig. 2. The two approximate scalings cross at f M = 2/3, where, indeed, the fitted central charge c has a sudden drop as shown in Fig. 2 (central panel), and is compatible with a gapped phase (insulator B). A second drop is found at f M = 1, another gapped phase (insulator A), while elsewhere the central charge is compatible to a constant c = 1 signaling a Luttinger liquid phase (see next section).
The CDW order, emerging in each of the gapped phases, is signalled by the local order pa- To avoid potential problems of symmetry restoration due to finite-size simulations, we extract the CDW order parameter from the static structure factor, at wave-vector k, of the matter density-density correlations [66,67] Fig. 3 we show the typical behaviour of the relevant order parameter ζ k , at filling f M = 1 3 ) in the left (right) panel as a function of the coupling parameter t. To identify the CDW ordered phase we compute ζ k for different system lengths L and couplings t and check when the order parameter survives at the thermodynamical limit lim L→∞ ζ k = 0. From numerical extrapolations (see the dot-dashed curves in Fig. 3), it appears that the CDW order is established for large t values, specifically t > t A c (resp. t > t B c ) for f M = 1 (f M = 2/3). The precise identification of the transition points t A c and t B c between the liquid and ordered phases is performed in the next section.
Transition points − We can employ the CDW order parameter data to estimate the transition values t A c and t B c , using various approaches. A first approach is the steepest slope method, which locates the transition point by taking the thermodynamical limit extrapolation lim L→∞ t of the parameter value t where the differential ∂ζ k /∂t| t is maximal. This method is known to provide robust results as long as the critical exponent β is smaller than 1 [66]. We estimated β A 0.25 ± 0.08 and β B 0.13 ± 0.05 for the transition into insulator A and B respectively, via a finite-size scaling procedure (not shown) and thus we can safely apply this procedure, which results in t A c = 12 ± 1 and t B c = 35 ± 2 as shown in Fig. 3 (insets).
Another independent approach to locate the critical points is to analyse the central charge of the model as a function of the system parameters, obtained from the entanglement entropy S fitted via Eq. (5), as reported in Fig. 4. As expected, as the thermodynamical limit is approached, a discontinuity is highlighted in the fitted c val- Additionally, it appears that the transition point is described by a critical theory different from the Luttinger phases. Indeed, the spike survives at increasing system length L and leads to a two-sided discontinuity in c as a function of t for increasing . Specifically, while  Fig. 4). The extrapolated values c (A) and c (B) have each been obtained via two independent methods, providing compatible results: The first method is the one previously discussed, and expressed by Eq. (5) (bullet points in the insets of Fig. 4). The second method (boxes points) is based on a postprocessing on the entanglement entropy profile, in order to remove the Fermi oscillations and recover a simple conformal trend S 0 , whose c can be easily fitted. This post-processing is somewhat technical and extensively discussed in Ref. [68]. Our results suggest that eventual logarithmic corrections [69,70,71], if present, do not visibly affect the quality of our fits. (see examples shown in Fig. 2, obtained via the first method). Moreover, by extrapolating the position t * of the peak itself at the thermodynamical limit (upper insets in Fig. 4), we can locate the transitions points t A c and t B c . Indeed, this procedure provides results compatible with the previous method t A c = 12 ± 1 and t B c = 33 ± 2. The values we obtained for the central charge are compatible with the SU(2) 2 WZNW model, which is an expected critical behaviour for gauge theories [35]. (1) exhibits a gapless phase with central charge c = 1, is already a strong signature that such phase is a band conductor (albeit an interacting one) and its long-range properties at the thermodynamical limit are captured by the Luttinger liquid paradigm [72,73]. In this section we aim to confirm the existence of a mesonic BCS phase, where the system behaves like a Luttinger liquid of bound fermion pairs (forming a spin singlet together) rather than a liquid of single fermion quasiparticles (simple metal). This is predicted for small coupling t from a perturbation theory analysis presented in the methods section. However, even for strong t we also observe emergence of BCS behaviour at low fillings: In Fig. 2, looking at the Fermi oscillations in the entanglement entropy, we clearly detect the shift from a scenario at f M 1 where the number of oscillations is equal to the number of matter fermion holes with respect to insulator A (k f ≈ π(1 − f M ), single fermionic quasiparticle liquid), to a scenario at f M 0 where the number of oscillations is equal to half the matter fermions in the system (k f ≈ π 2 f M , meson BCS liquid) [74].
We confirm this picture by studying the quasi- j,↑ of the superfluid mesons. Since particle number conservation cannot be explicitly broken in (1+1)-d due to the Mermin-Wagner-Hohenberg theorem [75,76], true long-range order cannot be established. Therefore, we investigate the correlation length ξ of the meson superfluid order, defined either as [77] where C is the bulk average of −η e − /ξ . We found that the two methods provide compatible results: henceforth, we report those obtained via the first method. Typical results of the correlation length are shown in Fig. 5, where we specifically consider the filling values which allow for insulating phases (f M = 2 3 , 1): The correlation length grows linearly with the system size L inside the whole gapless region t < t A c (t < t B c ), while it stays finite and scales as ξ ∝ (t−t A,B c ) ν A,B in the insulating phase. That is, the meson BCS order survives until the insulating phase is established. An analysis of the criti-cal exponent ν when approaching the BCS phase from the insulators is carried out in Fig. 6. For other filling fractions we also observed linear divergence of the correlation length, and especially at high matter filling f M 1 this occurs only at small coupling t <t, despite the central charge remaining constant c = 1. This is, again, a clear signature of the transition from the meson BCS liquid into a single fermion quasiparticle liquid, as we previously argued from the entanglement entropies.
The final phase diagram summarizing all our findings is sketched in Fig. 1.

Discussion
In this manuscript we applied gauge invariant tensor networks methods to study the finite density phase diagram of a non-Abelian lattice gauge theory. This model was built to include the basic ingredients of the Yang-Mills theory for a SU(2) gauge symmetry, following a Kogut-Susskind formulation, and it has been recast in a quantum link formalism. By means of TN numerical ansatz for quantum link models, we numerically investigated the phase diagram at the thermodynamical limit. At small matter-field coupling t we confirmed the existence of a meson BCS conductor, a Luttinger liquid phase of matter fermion pairs, which agrees with an analytical prediction based on second-order degenerate perturbation theory. We observed that the meson BCS may not survive at larger coupling t, depending on the matter filling fraction f M , eventually in favor of a simple liquid phase. Moreover, chiral symmetry breaking at strong coupling leads to the presence of two insulating phases, at specific matter filling fractions f M = 2 3 and f M = 1: Both these insulators exhibit charge density wave order, respectively at k = 2 3 π and k = π. Finally, using entanglement entropies, we characterised in detail the transition points t showing that they are compatible with a SU (2) 2 Wess-Zumino-Novikov-Witten critical theory.
Our work paves the way towards nonperturbative studies of non-Abelian gauge theories in finite-density regimes where MonteCarlo simulations are biased by the sign problem. Since the local Hilbert space dimension only depends weakly on the symmetry group, it is possible to extend the present study to SU(3) models, which promise to host an even richer pairing scenario. In this context, as we have shown above, while low-dimensional models differ significantly to their higher dimensional counterparts in manyrespects, they still host several common phenomena, including confinement and the restoration of chiral symmetry at finite chemical potential. The methods applied here can immediately be extended to ladder systems, and potentially to (2+1)-d lattices, providing a complementary route to recent approaches based on projected entangled paired states [78,79]. In particular, the possibility of directly accessing entanglement properties enables the exploration of entanglement properties in gauge theories, a topic which has recently surged to attention [39,40] and where key questions, such as the role of sub-leading corrections to the area law, promise to give tremendous insights on both critical and topological phenomena [80,81].
In this study we were not interested in the continuum limit of the theory, but studied the zero temperature regime of the regularised lattice version. In this sense, our results are directly related to possible realisation in cold atom systems, where the lattice theory can be implemented [17,10,12]. We foresee a possible study of the continuum limit of the theory carried out via, for instance, dimensional reduction. This method has already given successful results when studying gauge theories with tensor networks, and showed that continuum limit extrapolations can be carried out even for moderately small dimension of the truncated gauge fields [51,49,82,83]. In the future, we plan to investigate this aspect within the quantum link formulation we used here.

Methods
Gauge symmetry − By definition, gauge invariance implies that H commutes with the local generators J 34 for f M = 2/3. Insets: Correlation length rescaled by the total length L of the system, thus ξ/L. In the gapless phases this quantity tends to a finite nonzero value at the thermodynamical limit, which signals quasi long-range order of meson superfluidity.
and a right group or rotations This suggests that the compound gauge transformation generator J j at site j can be naturally decomposed into right-gauge, matter, and leftgauge components (11) Its commutation relations with the gauge field, which read respectively j,s +c [τ ] † j,s ), which is a Z 2 group, since Π 2 = 1.
Boundaries − As we simulate a onedimensional lattice gauge system in Open Boundary Conditions (OBC), the rishon modes at the boundaries remain uncoupled, precisely c 1,s and c [L] L,s . The quantum state over these modes is left invariant by the quantum dynamics. Nevertheless, the actual space of 'physical' states is affected by the boundary conditions. Throughout our simulations, we set specific Von-Neumann boundary conditions such that there is no SU(2)charge at the boundaries, or equivalently, no free-field energy density at the boundaries, e.g. n with g 1 = g 0 3 8 . This model is exactly solvable in the case of zero coupling t = 0: Indeed, when H coupl is absent, the remaining Hamiltonian contributions act trivially (as the identity operator 1) upon the matter modes, while every quantum link (L j , R j+1 ) does not interact with any other link. Therefore, at zero temperature, every quantum link is in the state Such a quantum link state, which is also graphically sketched in Fig. 1 (bottom-right panel), minimises H free and H break separately, and its energy density is − . The full ground state will thus be a tensor product state where rishon modes are in |Link 0 state and matter models are in any arbitrary state which satisfies the Gauss' law: namely, either in |0 j,M or in |2 j,M . This ambiguity produces a ground state degeneracy of 2 L configurations, each of them with total ground energy E 0 = −(L − 1) . Such degenerate ground space states can be classified as spin-1 2 chain states, according to where {m 1 ..m L } is any binary string and σ − j = c [M ] j,↓ c [M ] j,↑ obey the standard Lie algebra of Pauli lowering operators: † creates a pair of matter fermions at site j, which we will call 'meson' in analogy to high-energy physics, and possesses hardcore boson statistics.
For small coupling regimes, t g 2 1 , , it is possible to develop a degenerate perturbation theory analysis.
First order perturbation theory leads to no contribution whatsoever.
Second order degenerate perturbation theory allows us to construct an effective Hamiltonian H eff over the unperturbed ground space, such that Ψ m where the sum runs over the excited states |ε of the unperturbed (t = 0) Hamiltonian. After a lengthy but straightforward algebra, the resulting H eff can be rewritten in the spin-1 2 language introduced before as an antiferromagnetic Heisenberg model As the Heisenberg model is exactly solvable by Bethe ansatz, a gapless phase is realised in the limit t → 0. Finally, exploiting the mapping into a Hubbard-like Model, we expect -as confirmed by our numerical analysis -a phase of superfluid mesons for small matter-field coupling, which can be equivalently interpreted as a BCS conductor (since mesons are matter fermion pairs).
Gauge invariant tensor networks − We simulate numerically the model of Eq. (1) in equilibrium at zero temperature,i.e. in the ground state, via variational algorithm based on tensor network ansatz states for lattice gauge models introduced in Ref. [45]. More precisely, we adopt a Matrix Product State (MPS) ansatz [19,20,84,22,36,85], with additional embedding of both lattice gauge and global symmetries. Each simulation is performed on OBC, for finite lattice size L. Numerical convergence is tested as a function of the various numerical approximation parameters (e.g. the MPS bond link dimension χ, see later on). The thermodynamical limit is then estimated by extrapolations to L → ∞ of the various observable quantities. In conclusion, the equilibrium properties of Hamiltonian (1) are studied as a function of parameters g 1 , t, and finally the matter filling f M = N M /L, which we can control in a canonical fashion.
Hereafter we review and discuss the simulation details and the approximations involved. We first introduce the convention for classifying and labelling the local (canonical) basis states: we consider the three-mode compound site j of the lattice (precisely, modes R j , M j and L j ) as a single site j for our simulation. Within this space, we construct the effective local basis {|b j }, which includes all and only the d = 14 invariant states |b j which satisfy the Gauss' law |J j | 2 |b j = 0, where |J j | 2 = µ=x,y,z J (µ)2 j . Starting from this construction of the local basis space, we build a quantum many-body (QMB) tailored wave-function ansatz based on the Matrix Product State formulation (17) with OBC, where the variational information is conveniently stored in rank-3 tensors A [j] . The variational complexity of ansatz (17) is determined by χ, which plays the role of a refinement parameter, and is often referred to as bond link dimension.
Finally, we encode the additional symmetries within this variational picture: the global particle conservation, which translates into the selection rule j n [M ] j = N M , and the link symmetry, whose selection rule reads n j+1 for every j. This is accomplished by associating quantum numbers to the auxiliary indices α j [86,87,45], thus considering α j = (q j , j , t j ), where q tracks the matter particles number, j carries information on the link symmetry at link (L j , R j+1 ), and t j contains all the residual information (symmetry degeneracy index). Using this formalism one obtains where the first Kronecker delta fixes the matter conservation symmetry (once we set q 0 = 0 and q L = N M ), and the other deltas fix the link symmetries. The residual tensors R [j] contains the variational parameters to be optimised.
To find the ground state of Hamiltonian (1) within the variational space given by Eqs. (17) and (18), we adopt an annealing scheme (imaginary time evolution) implemented via the Time Evolved Block Decimation (TEBD) method [38,88], evolving an initial (random) state |Ψ 0 in imaginary time until convergence is reached. Numerical resources and precision − For the system lengths that we considered (20 ≤ L ≤ 150) we employed a MPS bond link dimension up to χ ∼ 400. We observed that such dimensions were sufficient to provide a robust phase analysis. Specifically, we could reach precisions of the order ∼ 10 −11 in the entanglement spectra for gapless phases (as shown in Fig. 7), while we easily hit machine precision (λ min 10 −32 ) for gapped phases. In turn, this is expected to produce a precision of the order of √ λ min ∼ 10 −5 or better in the expectation values of local observables and correlators, even in the gapless phases. We confirmed this estimation by comparing simulations carried out at different bond links: Fig. 8 shows such a comparison, between ξ = 250 and ξ = 300, and highlights that the errors of the quantities involved are of the order ∼ 10 −5 or below.

Acknowledgements
The authors acknowledge fruitful and stimulating discussions with Peter Zoller and Uwe-Jens Wiese. SM gratefully acknowledges the support of the DFG via a Heisenberg fellowship, and the periodic visits at IQOQI during which part of this work has been developed. We acknowledge support from the EU via the RYSQ, SIQS and SCALEQIT projects, the ERC Synergy Grant UQUAM, the DFG via the SFB/TRR21 and SFB FoQuS (FWF Project No. F4016-N23), and the Baden-Württemberg Stiftung via Eliteprogramm for Postdocs. We also acknowledge financial support from Spanish MINECO FIS2015-69983-P, UPV/EHU Project No. EHUA15/17, UPV/EHU UFI 11/55.