A Scalable Translationally Invariant Variational Theory of Ab Initio Polarons

This paper introduces a scalable, translationally invariant variational theory for ab initio polarons that combines momentum-projected wavefunctions with low-rank kernel factorization to accurately model carrier behavior across coupling regimes in the thermodynamic limit, revealing significant biases in existing diagrammatic Monte Carlo results for strong-coupling hole polarons in LiF.

The Gist

The Big Picture: The "Dressed" Electron

Imagine an electron moving through a solid crystal (like a piece of salt or a semiconductor) as a person walking through a crowded dance floor.

• The Electron: The person walking.
• The Lattice: The crowd of people (atoms) dancing.
• The Polaron: When the walker moves, they bump into people, causing the crowd to shuffle and rearrange around them. The walker is now "dressed" in a cloud of moving people. This combined package (walker + crowd) is called a polaron.

Scientists have long wanted to calculate exactly how heavy this "dressed" package is and how fast it can move. However, doing this math is incredibly difficult because the crowd is huge, and the interactions are complex.

The Problem: The "Supercell" Trap

Previous methods to solve this problem had two main flaws:

1. They were too slow: To get accurate answers, scientists had to simulate a tiny, artificial chunk of the material (a "supercell") and repeat it over and over. This is like trying to understand how a whole city moves traffic by only studying one single block. It's computationally expensive and often inaccurate.
2. They were biased: Some methods worked well if the walker was moving slowly (weak coupling), while others worked well if the walker was stuck in a deep hole created by the crowd (strong coupling). No single method could handle both situations accurately without breaking the math.

The Solution: A New "Scalable" Theory

The authors (Baumgarten, Wu, Jiang, and Lee) introduced a new mathematical framework that solves these problems. Think of their approach as a new way to simulate the dance floor that doesn't require building a fake city block.

1. The "Momentum Projected" Wavefunction (The Magic Mirror)
Imagine you have a photo of a person standing still in a crowd (a localized state). In the old methods, you had to pick a specific spot for the person, which broke the symmetry of the room.
The authors use a trick called momentum projection. Imagine taking that photo of the person and creating a "ghostly superposition" where the person is simultaneously standing in every possible spot on the dance floor at once. This restores the natural symmetry of the crystal. It allows the math to describe a polaron that is either stuck in one spot (strong coupling) or zooming freely across the whole room (weak coupling) using the same set of rules.

2. The "Low-Rank Factorization" (The Compression Trick)
The math behind electron-crowd interactions usually involves a massive spreadsheet of numbers that gets too big to handle as the simulation gets bigger.
The authors used a technique called low-rank factorization.

• Analogy: Imagine you have a 10,000-page instruction manual on how the crowd reacts. Instead of reading every single page, you realize that 99% of the instructions are just variations of the same 50 core rules.
• By compressing the data into these "core rules" (singular vectors), they reduced the computational cost. Instead of the time needed growing quadratically (getting much slower as the grid gets bigger), it now grows almost linearly. This means they can simulate a massive, dense crowd (a dense grid of points) on a standard computer without waiting years for the result.

What They Found (The Benchmarks)

They tested their new method on four different materials: Lithium Fluoride (LiF), and two types of Titanium Dioxide (Anatase and Rutile).

• The "Gold Standard" Check: They compared their results against a method called DiagMC (Diagrammatic Monte Carlo), which is considered a very accurate, unbiased benchmark.
• The Surprise:
  • For weak-coupling cases (like the electron in LiF), their new method matched DiagMC perfectly.
  • For strong-coupling cases (like the hole in LiF), their new method agreed with other reliable methods (VMC), but disagreed significantly with the published DiagMC results.
  • The Conclusion: The authors suggest that the DiagMC results for the strong-coupling LiF hole were likely biased or inaccurate due to sampling errors. Their new method, being "translationally invariant" (symmetric), seems to be the more reliable truth in these tough scenarios.

Real-World Visualization

The paper didn't just calculate numbers; they visualized the "shape" of the polaron.

• LiF Electron: The polaron is a large, fluffy cloud spreading out evenly in all directions (isotropic).
• Rutile Electron: The polaron is a tight, compact ball.
• Anatase Electron: The polaron is a flat, pancake-like shape (anisotropic), spreading out in two dimensions but staying thin in the third.

Summary

This paper presents a new, faster, and more accurate way to calculate how electrons interact with the atoms they move through.

1. It's Scalable: It can handle huge, realistic simulations without needing supercomputers to run for centuries.
2. It's Universal: It works for both "free" electrons and "stuck" electrons.
3. It's Corrective: It revealed that a previous "gold standard" calculation might have been wrong for certain difficult cases, offering a more trustworthy path forward for understanding materials.

In short, they built a better, faster, and more symmetrical lens to see how electrons move through the solid world.

Technical Summary

Technical Summary: A Scalable Translationally Invariant Variational Theory of Ab Initio Polarons

Problem Statement
Polarons, charge carriers dressed by lattice vibrations, are central to condensed-matter phenomena. However, a comprehensive ab initio description that remains translationally invariant, applicable across all electron-phonon coupling regimes (weak to strong), and computationally tractable on dense Brillouin-zone meshes has been lacking. Existing perturbative and Green's-function methods struggle in strong-coupling regimes, while traditional variational approaches either favor strong coupling (e.g., localized adiabatic states) or become prohibitively expensive for thermodynamic-limit (TDL) extrapolations due to unfavorable scaling with $k$-point sampling. Furthermore, first-principles diagrammatic Monte Carlo (DiagMC), while unbiased in principle, faces significant sampling difficulties in strongly interacting regimes.

Methodology
The authors introduce a scalable, translationally invariant variational framework that bridges the gap between weak and strong coupling without resorting to supercells. The methodology consists of two primary components:

1. Momentum-Projected Toyozawa-Type Wavefunction:
   The approach starts with the Landau–Pekar (LP) product state, which describes a localized electron coupled to a coherent state of lattice deformation (Davydov D2 ansatz). While accurate for strong coupling, the LP state breaks translational symmetry. To restore this symmetry and capture delocalized carriers, the authors apply a Peierls–Yoccoz (PY) projection operator to the D2 state. This constructs a momentum eigenstate (labeled by crystal momentum $K$) that retains the localized physics of the adiabatic state while recovering the delocalized character required in the weak-coupling limit. This results in the delocalized D2 (dD2) wavefunction.

2. Low-Rank Factorization for Scalability:
   The direct evaluation of the variational energy for the dD2 wavefunction typically scales as $O(N_k^2 N_b^2 N_{mod})$, where $N_k$ is the number of $k$-points, $N_b$ the number of bands, and $N_{mod}$ the number of phonon modes. This quadratic scaling in $N_k$ makes dense sampling infeasible. The authors overcome this by employing a recently introduced low-rank representation of the short-range electron-phonon (eph) coupling kernel. By expressing the coupling matrix elements $g_{mn\nu}(k, q)$ in terms of singular vectors and Wannier-to-Bloch transformation matrices, the contraction over $k$ and $q$ can be reorganized. This allows the dominant computational bottleneck to be evaluated using Fast Fourier Transforms (FFT), reducing the scaling to near-linear, $O(N_c N_k \log N_k N_b^2 N_{mod})$, where $N_c$ is the number of retained singular vectors.

Key Contributions

• Unified Framework: The method provides a single variational ansatz capable of describing polarons from the weak-coupling (delocalized) to the strong-coupling (self-trapped) regimes without switching formalisms.
• Computational Efficiency: The integration of low-rank factorization with momentum projection enables near-linear scaling with Brillouin-zone sampling, making TDL extrapolations feasible for realistic materials.
• Systematic Improvability: The dD2 wavefunction serves as a foundation for systematically improvable hierarchies, such as adding perturbative corrections (e.g., CSPT2) on top of the mean-field solution.

Results
The theory was benchmarked against the Fröhlich model and first-principles calculations for LiF, anatase TiO$_2$, and rutile TiO$_2$.

• Fröhlich Model: The dD2 method successfully recovers the correct weak-coupling ($E \propto \alpha$) and strong-coupling ($E \propto \alpha^2$) limits, outperforming the unprojected LP state and matching the accuracy of all-coupling Green's function methods.
• LiF (Strong Coupling): For the hole polaron in LiF, the dD2 and D2 methods yield a binding energy of 1.933 eV, in excellent agreement with Variational Monte Carlo (VMC) and CSPT2. Notably, these results are significantly lower than the previously published DiagMC value (2.26 eV), suggesting a systematic bias in the DiagMC results for this strong-coupling case.
• LiF (Weak Coupling): For the electron polaron in LiF, the momentum-projected dD2 yields a binding energy of -0.395 eV, closely matching the DiagMC value (-0.408 eV) and VMC. In contrast, the unprojected D2 underbinds significantly, failing to capture the delocalized nature of the carrier.
• TiO$_2$ Polymorphs:
  • In rutile, the electron polaron is relatively localized; dD2 improves upon D2, but CSPT2 (based on an undisplaced reference) yields the lowest energy, indicating a regime where low-phonon excitations dominate.
  • In anatase, the electron polaron is highly delocalized. D2 predicts a negligible binding energy, whereas dD2 stabilizes the polaron with a binding energy of -0.138 eV, consistent with the need for translational invariance in weak-to-intermediate coupling.
• Band Structures and Extent: The method provides direct access to TDL band structures and real-space correlation functions. For the LiF electron polaron, the dD2 band structure agrees well with DiagMC near the $\Gamma$ point but lies systematically lower at finite momenta, further indicating potential inaccuracies in the DiagMC sampling or Hamiltonian approximation in those sectors. Real-space analysis reveals distinct anisotropy and size differences among the polarons, with LiF hole polarons being the most localized and anatase electron polarons being the most extended.

Significance
The paper establishes momentum-projected variational wavefunctions as a powerful, systematically improvable route for studying polarons in real materials at the thermodynamic limit. By resolving the computational bottleneck of dense Brillouin-zone sampling, the method offers a scalable alternative to existing approaches. The comparison with DiagMC highlights potential systematic errors in stochastic sampling for strong-coupling regimes and inaccuracies in Hamiltonian approximations used in previous studies. The work demonstrates that a translationally invariant mean-field approach, when combined with low-rank factorization, can accurately capture both self-trapped and delocalized carriers, providing a robust foundation for predictive polaron physics.