Effective Hamiltonian based Monte Carlo for the BCS to BEC crossover in the attractive Hubbard model

We present an effective Hamiltonian based real-space approach for studying the weak-coupling BCS to the strong-coupling Bose-Einstein condensate crossover in the two-dimensional attractive Hubbard model at finite temperatures. We introduce and justify an effective classical Hamiltonian to describe the thermal fluctuations of the relevant auxiliary fields. Our results for $T_c$ and phase diagrams compare very well with those obtained from more sophisticated and CPU-intensive numerical methods. We demonstrate that the method works in the presence of disorder and can be a powerful tool for a real-space description of the effect of disorder on superconductivity. From a combined analysis of the superconducting order parameter, the distribution of auxiliary fields, and the quasiparticle density of states, we identify the regions of metallic, insulating, superconducting, and pseudogapped behavior. Our finding of the importance of phase fluctuations for the pseudogap behavior is consistent with the conclusions drawn from recent experiments on NbN superconductors. The method can be generalized to study superconductors with nontrivial order-parameter symmetries by identifying the relevant auxiliary variables.

Figure 1

A schematic picture describing the method to isolate the contribution of a single rotor pair to the total energy for the effective classical Hamiltonian. (a) $i$ th rotor is oriented away from the otherwise phase coherent arrangement of rotors. The double lines connecting site $i$ to all other sites indicate the pairs that contribute to the change in energy due to change in the orientation of the $i$ th rotor. (b) $j$ th rotor is rotated by an angle $\theta$, and the single lines indicate the pairs contributing to change in energy. (c) Both $i$ th and $j$ th rotors are rotated by an angle $\theta$. Note that in this case the pair $ij$ does not contribute to the change in energy.

Figure 2

(a)–(f) Change in energy as a function of the orientation angle between a single nearest-neighbor pair of rotors for different values of $U$. Symbols are the results of numerical calculations and the solid line in each panel is a fit to the functional form $J(1-cos\theta )+K(1-{cos}^2\theta )$. The best-fit values of $J$ and $K$ are indicated in the figure.

Figure 3

(a)–(f) Change in energy as a function of the change in magnitude of a single rotor. Symbols are the numerical data and the solid line shows a fit to the functional form $\delta E=k\left(\right|\mathrm{\Delta }|-|{\mathrm{\Delta }}_0{\left|\right)}^2$, with the values of best-fit parameter noted in the figure.

Figure 4

Parameters of $H_{\mathrm{cl}}$ as extracted from the change in energy about the mean-field value. (a) Coupling constant $J$ (open circles) as a function of $U$ calculated by assuming a cosine form for the change in energy and using only $\theta =0$ and $\pi$. Filled squares show a comparison with $J$ obtained from the fits shown in Fig. 2. Green solid line is the phase stiffness calculated as expectation value of the kinetic-energy operator. The dashed blue line represents the $1/U$ behavior valid in the large-$U$ limit. (b) The stiffness $k$ to change in local amplitude $|{\mathrm{\Delta }}_i|$, as a function of $U$ calculated by assuming a $k\left(\right|\mathrm{\Delta }|-|{\mathrm{\Delta }}_0{\left|\right)}^2$ form and using only two points from the data (open circles). Filled squares are the values obtained from the best fits shown in Fig. 2. Dashed blue line represents $k=U$.

Figure 5

(a) The temperature dependence of the superconducting order parameter ${\mathrm{\Delta }}_{\mathrm{op}}$ normalized by its low-temperature value, ${\mathrm{\Delta }}_{\mathrm{op}}\left(0\right)$, for different values of $U$. Inset shows the variation of ${\mathrm{\Delta }}_{\mathrm{op}}\left(0\right)$ with $U$. (b) Vortex $n_v$ and antivortex $n_a$ density as function of temperature. (c) Spectral gap ${\mathrm{\Delta }}_g$ as a function of temperature for different $U$. (d) Circles show the transition temperature $T_c$, as inferred from the inflection point in the $T$ dependence of the order parameter, as a function of $U$. Squares and triangles mark, respectively, the expected variations of $T_c$ in the small-$U$ and large-$U$ limits.

Figure 6

Variation of the pair-pair correlation function $P_S$ with inverse temperature for different system sizes. The results are obtained for $U=4t$ and $n\approx 0.5$.

Figure 7

The distribution $P\left(\right|\mathrm{\Delta }\left|\right)$ of magnitudes $\left|\mathrm{\Delta }\right|$ at different temperatures, $0.1T_c, 1.0T_c, 2.0T_c$, for (a) $U=1.5$, (b) $U=3$, and (c) $U=10$. (d) The ratio of the variance to the mean value of $P\left(\right|\mathrm{\Delta }\left|\right)$ as function of $U$. This ratio decreases with increasing $U$ highlighting the importance of amplitude fluctuations at weak $U$ values.

Figure 8

The distribution $P\left(D\right)$ of the nn phase correlators $\text{cos}({\varphi }_i-{\varphi }_j)$ (see text) at different temperatures, $0.1T_c, 0.5T_c, 1.0T_c, 2.0T_c$, for (a) $U=1.5$, (b) $U=3$, and (c) $U=10$. (d) The maximum value, $P_{\text{max}}$, of $P\left(D\right)$ as a function of $T$ for various $U$ values.

Figure 9

Real-space plots for two different temperatures $T\sim 0.1T_c$ (first column) and $T\sim T_c$ (second column) for $U/t$ = 1.5. (a) and (b) show the nn phase correlations $D_{ij}=cos({\varphi }_i-{\varphi }_j)$, and (c) and (d) show the amplitude variables $|{\mathrm{\Delta }}_i|$.

Figure 10

Real-space plots for two different temperatures $T\sim 0.1T_c$ (first column) and $T\sim T_c$ (second column) for $U/t$ = 16. (a) and (b) show the nn phase correlations $D_{ij}=cos({\varphi }_i-{\varphi }_j)$, and (c) and (d) show the amplitude variables $|{\mathrm{\Delta }}_i|$.

Figure 11

Variation of quasiparticle density of states, $N\left(\omega \right)$, with temperatures for coupling strengths (a) $U=1.5$, (b) $U=3$, (c) $U=6$, and (d) $U=10$. For small $U$, the spectral gap vanishes at $T_c$ as expected in the BCS regime.

Figure 12

(a) Temperature variation of coherence length $\left(\xi \right)$ for different $U$. (b) Variation of $\xi$ as function of $U$ at low temperatures. The filled symbols represent the Monte Carlo data, and the other three data sets are for the hypothetical auxiliary field configurations corresponding to fluctuations in magnitude of ${\mathrm{\Delta }}_i$, phase of ${\mathrm{\Delta }}_i$, and both phase and magnitude of ${\mathrm{\Delta }}_i$ (see text).

Figure 13

(a) The $T-U$ phase diagram showing the superconducting, the normal metal, the nonsuperconducting gapped, and the pseudogapped phases. Symbols are the data points obtained from simulations and the dashed lines are guides to eye. (b) The ratio of spectral gap to ordering temperature, $\frac{{\mathrm{\Delta }}_g}{T_c}$, as a function of $U$. The ratio is shown for the gap at $T=0$ and that at $T\sim T_c$.

Figure 14

(a)–(f) Change in energy as a function of the orientation angle between a single pair of rotors for different values of $U$ and $V$. The scatter of points is due to the inequivalence of nearest-neighbor pairs. In order to show the variation for different pairs on the same scale, we have normalized the variation in energy for each pair by its maximum value, hence all the points lie between zero and one. The solid line in each panel is the function $(1-cos\theta )/2$.

Figure 15

Probability distributions for coupling $J$ for (a) $U=4$ and (c) $U=10$ and probability distributions for amplitude stiffness $k$ for (b) $U=4$ and (d) $U=10$ at various disorders.

Figure 16

Temperature dependence of the superconducting order parameter normalized by its low-temperature value ${\mathrm{\Delta }}_{\mathrm{op}}/{\mathrm{\Delta }}_{\mathrm{op}}\left(0\right)$ for (a) $U=4$ and (b) $U=6$. Inset in (b) shows the variation of ${\mathrm{\Delta }}_{\mathrm{op}}\left(0\right)$ with disorder strength $V$. The phase diagram in the temperature-disorder plane showing superconducting, nonsuperconducting gapped, and pseudogapped regimes for (c) $U=4$ and (d) $U=6$.