Electronic polymers and soft-matter-like broken symmetries in underdoped cuprates

Empirical evidence in heavy fermion, pnictide, and other systems suggests that unconventional superconductivity appears associated to some form of real-space electronic order. For the cuprates, despite several proposals, the emergence of order in the phase diagram between the commensurate antiferromagnetic state and the superconducting state is not well understood. Here we show that in this regime doped holes assemble in"electronic polymers."Within a Monte Carlo study we find, that in clean systems by lowering the temperature the polymer melt condenses first in a smectic state and then in a Wigner crystal both with the addition of inversion symmetry breaking. Disorder blurs the positional order leaving a robust inversion symmetry breaking and a nematic order, accompanied by vector chiral spin order and with the persistence of a thermodynamic transition. Such electronic phases, whose properties are reminiscent of soft-matter physics, produce charge and spin responses in good accord with experiments.

Intensity of the charge peak where c i,σ (c † i,σ ) destroys (creates) an electron with spin σ at the site i, and n i, Mapping to an effective XY model. Since the textures we considered are long-ranged and the cluster has periodic boundary conditions, the energy of each configuration is strongly dependent on the system size which is limited to small clusters (16 × 16) in the Hubbard model. One expects, however, that the short-range part of the interaction, which is due to the overlap of the hole-wave functions within the core of the segment, is well converged in our cluster. Fortunately, because the lowest energy textures we find are planar, the long-range part of the interaction can be well reproduced with an effective XY model which allows to extrapolate the results to very large clusters.
In the effective XY model the localized holes are replaced by vacant sites and an effective frustrating interaction is introduced between the spins in sites nearest neighbors to the holes. A similar approach was used in Ref. 4 although in that paper the minimal topological dipole corresponds to one hole instead of two holes as in our case. In fact in our GA computations we find that the one hole state is a collinear spin polaron and only for two holes or more the GA yields the dipolar distortion.
Supplementary Figure 1(b) shows an example with the green segments representing the effective interactions across the chain of VA pairs.
The effective hamiltonian we use reads, where J ij = J for nearest-neighbor sites if i and j are not vacant sites, J ij = J > 0 across each VA pair (which plays the role of an elementary topological dipole) and J ij = 0 elsewhere. J is adjusted to match the stiffness of the Hubbard model in the GA and J is adjusted so that the dipole strength of the two models agree as explained below.
In order to obtain the dipole strength p 0 of a dipole oriented along the (1, 1) direction, we compute the average spin current in the directionû = 1 √ 2 (1, −1) by averaging over a line parallel to the dipole moment. The spin current can be shown from the mapping to a Coulomb gas model to be directly proportional to the dipole strength, Then we adjust J of the XY model in such a way to match the dipole strength of the Hubbard model computed in the same way. This ensures that the long-range part of the energy of the XY model coincides with the long-range part of the energy of the Hubbard model for the same texture.
Supplementary Figure 2 shows the configurations considered to obtain the effective interaction among the topological charges once the mapping to the Coulomb gas is done.
In order to extrapolate to large clusters we assume that the energy of each configuration for the cluster of size L can be written as the sum of short-range (sr) terms plus long-range (lr) terms, E Hub,L×L = E Hub,sr,L×L + E lr,L×L .
By construction E lr,L×L coincides with the same term in the effective XY model, therefore we can write E XY,L×L = E XY,sr,L×L + E lr,L×L .
The energy of the Hubbard model for a large cluster of size L is obtained assuming that the short-range terms do not depend on the system size according to Next according to the Coulomb gas mapping 5 we assume that the long-range part of the energy can be represented by topological charges of strength ±k.
Examining the GA textures we find that topological charges are centered close to (but not exactly at) lattice sites. Since in the Monte Carlo computation only the dipole strength p = ±kl plays a role we use the freedom on k and l to force the topological charges to be located exactly at lattice sites. k was determined in such a way that when a dipole p is formed with charges on the diagonal of a plaquette of length l = √ 2a, the dipole moment p 0 = ±kl reproduces that of the GA and XY computations yielding k = 0.8.
Next, we decompose E Hub,L ×L for N h holes in the different contributions, Here, E AF,L ×L is the energy of the antiferromagnetic solution in a cluster of the same size, r ij is the distance between topological charges i and j, η i = ± is the sign of the topological charges, 0 is a self-energy term independent of the distance among the topological charges. The last two terms represent the long-and short-range parts of the interaction. Our remaining task is to determine the short-range interaction. Supplementary Figure 3 shows the total potential V (r) between two topological charges of the same sign (a) and of different sign (b), as a function of their distance r, including short-and long-range parts (green lines) and compared with the long-range part alone (red lines).
We have assigned high energy values (10t) to the terms δV +− r=1 and δV ++ r=2 , since these configurations cannot be stabilized even as metastable states. This effectively has the effect of eliminating these configurations from the Monte Carlo computations, so that the computations will not depend on the precise value of these constants.
Notice that the effective interaction between topological charges should not be taken literally.
Indeed if two holes are separated at a long distance the Coulomb gas energy is logarithmically large because it assumes the vortex and antivortex configuration to persists indefinitely. Instead for large distances the GA solution converges to two collinear spin polarons which have a smaller energy.
Thus for those configurations the effective interaction breaks down. The interaction makes sense only if the VA pairs are at short distances and form dipoles. Thus in the Monte Carlo computations we restrict the temperature range in such a way that this condition is met.

Supplementary Note 2: Treatment of quenched disorder due to dopants
In the effective Coulomb gas model for the Monte Carlo simulations, we take into account the disorder due to out-of-plane ions. Since each V or A is associated with a hole carrier, we associate to each V or A also a real positive Coulomb charge. These hole charges feel the disordered potential out of the plane generated by the doping out-of-plane ions.
We use as Supplementary Figure 4 shows the potential field generated by a given configuration of disorder in the plane where the holes move,  show that for Q ion /Q rep = 0.125, the low temperature width of the peaks in the spin structure factor is very similar to the experimental structure factor (see Refs. [9][10][11]. It also show that, as a function of temperature, the structure factor evolves from incommensurate peaks at low temperatures to a commensurate peak above the ferronematic ordering temperature reported in the main article (45 K). The same data are plotted convoluted with a Gaussian in panel (a) of Figure 3 of the main article.

Supplementary Note 5: Nature of phase transitions
The study of phase transitions with Monte Carlo calculations on finite size clusters poses obvious problems. In particular it is of physical relevance to distinguish whether a vanishing order parameter is small because the symmetry of the phase is not broken or because the Monte Carlo evolution is "ergodically" exploring and spending equal times in regions with opposite values of the order parameters. To this purpose, to find the stablest phases at a given temperature we use the method explained in a previous work 6 , which monitors and records the Monte Carlo evolution of the order parameter(s). In our present case, we characterize Monte Carlo generated configurations by the instantaneous value of a three-dimensional (3D) order parameter defined as the vector (φ, P (1,1) , P (1,−1) ).
Each different phase is characterized by a "cloud" of points in different regions of the 3D space. In order to characterize these different phases, we compute the 3D density distribution of points at each temperature of the simulation. The points are expected to be distributed with a Boltzmann weight exp [−F (φ, P (1,1) , P (1,−1) )/(k B T )], so that higher density of points corresponds to the lower free energy F and indicates the stablest phase. This distribution can be visualized by plotting isosurfaces with a given number of points per unit volume. As an example we report in Supplementary Figure 9 this type of analysis for a low temperature state and an high temperature state, in a clean system.
Since at low temperature (Supplementary Figure 9(a)) the highest density of points occurs inside the two regions around φ ∼ +1 and P (1,1) ∼ ±1, and the two regions around φ ∼ −1, and P (1,−1) ∼ ±1, this phase is identified as ferrosmectic or ferronematic depending on the presence of charge ordering, which can be checked a posteriori from the charge structure factor. At high temperature [ Supplementary Figure 9(b)] the densest cloud of points occurs in correspondence of the zero value of the three order parameters indicating that at this temperature the system is in the disordered phase.
The phase transition can be visualized more easily by plotting 2D cuts of the 3D histogram, i.e. we plot only the plane P (1,−1) = constant, for φ > 0, or P (1,−1) = constant, for φ < 0. The constant is determined in such a way that the plane contains the bin with the maximum number of points of all the 3D histogram which usually yields a value for the constant nearly zero. So, if φ > 0, we take the plane P (1,−1) = 0 and plot the 3D bins in this plane as a 2D histogram, whose axes are φ and P (1,1) . In the case of Supplementary Figure 9(a) this yields a 2D histogram with maxima at φ ∼ 1 and P (1,1) ∼ ±1 while in the disordered case a single maximum appears at the origin (2D histograms not shown).
Using this procedure Supplementary Figure 10 shows that the phase transition from the disordered phase to the ferronematic phase remains sharp in the presence of realistic disorder. Indeed, in a narrow temperature range the maximum of the distribution shifts from finite values of polarization and nematicity (a) to the disordered phase characterized by a maximum of the distribution close to the origin (b).