First-principles estimation of the superconducting transition temperature of a metallic hydrogen liquid

In the conventional Eliashberg theory, $T_{c}$ is determined by solving a set of linearized Eliashberg equations. Liu et al. show rigorously that the same set of the equations still hold for liquids and general systems, provided that relevant parameters are properly determined. The set of the linearized Eliashberg equations reads

	$\displaystyle\rho\Delta_{n}$	$\displaystyle=$	$\displaystyle\sum_{n^{\prime}}\left[\lambda(n^{\prime}-n)-\mu^{\star}-\frac{\hbar\beta}{\pi}|\tilde{\omega}(n)|\delta_{nn^{\prime}}\right]\Delta_{n^{\prime}},\quad$		(1)
	$\displaystyle\tilde{\omega}(n)$	$\displaystyle=$	$\displaystyle\frac{\pi}{\hbar\beta}\left(2n+1+\lambda(0)+2\sum_{m=1}^{n}\lambda(m)\right),n\geq 0\quad$		(2)

and $|\tilde{\omega}(-n)|=|\tilde{\omega}(n-1)|$, where $n$ and $n^{\prime}$ are integers indexing the Matsubara frequencies of Fermions $\omega_{n}=(2n+1)\pi/\hbar\beta$, $\beta\equiv 1/k_{B}T$ and $T$ is the temperature. $\{\lambda(m),\,m\in Z\}$ is a set of interaction parameters characterizing the effective pairing interaction, and $\mu^{*}$ is the empirical Morel-Anderson pseudopotential introduced to take account of the exchange-correlation effect of electrons [24]. Solving Eq. (1) yields a set of eigenvalues $\rho$ and eigenvectors $\Delta_{n}$. For normal states, all eigenvalues are negative. The emergence of a non-negative eigenvalue indicates that the system enters into the superconducting state. One can thus estimate $T_{c}$ by solving the equations for a number of temperatures to see when the maximal eigenvalue changes sign. We note that $\Delta_{n}$ here is not related to the physical gap.

Two comments are in order. First, in calculations for solids using the standard approach (see, e.g., Ref. 7), it is common to use the empirical McMillan formula and estimate $T_{c}$ from a small number of parameters evaluated at the zero temperature including the mass enhancement factor $\lambda\equiv\lambda(0)$, an average phonon frequency and $\mu^{*}$. It simplifies the calculations. Unfortunately, such an approach does not work well for liquids because (a) unlike harmonic solids, $\lambda$ and the average phonon frequency of a liquid are temperature-dependent; (b) the empirical formulas are tuned for solids, and do not fit well with liquids which have qualitatively different density spectral functions due to the presence of low-frequency diffusive modes [25, 26]. Second, the choice of parameter $\mu^{*}$ is empirical. For metallic hydrogen systems, we use $\mu^{*}\approx 0.089$ as suggested in Ref. 7. The empirical parameter introduces an uncertainty in predicting $T_{c}$. Ideally, one should base calculations on the DFT for superconductors [27, 28, 29], which provides a rigorous framework for treating exchange-correlation effects in superconductors. This is not yet implemented in our current calculation.

To apply the Eliashberg equations, the interaction parameters $\lambda(m)$ are the central quantities to be determined. They are related to the effective pairing interaction $W$ by

$\displaystyle\lambda(n-n^{\prime})=-\sum_{\bm{k}^{\prime}}W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})\delta(\tilde{\epsilon}_{\bm{k}^{\prime}}-\epsilon_{F}),$

(3)

where $W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})$ denotes the matrix element of the effective pairing interaction, which is a function of the quasi-wave-vector transfer $\bm{k}-\bm{k}^{\prime}$ and the frequency transfer $\omega_{n}-\omega_{n^{\prime}}$, $\tilde{\epsilon}_{\bm{k}^{\prime}}$ is the renormalized quasi-electron dispersion, and $\epsilon_{F}$ is the Fermi energy. To apply the formula, it is necessary to determine the effective pairing interaction as well as the renormalized quasi-electron dispersion.

To determine these quantities in a liquid, Liu et al. propose a stochastic path integral approach. It is based on the classical isomorphism [30], i.e., mapping quantum ions to classical ring polymers by applying the imaginary-time path integral formalism [31]. The spatial configuration of the classical ensemble of ring polymers is specified by a set of $\tau$-dependent coordinates $\{\bm{R}_{a}(\tau)\}$, where $\tau\in[0,\hbar\beta)$ is the imaginary time or interpreted as a variable to parameterize the rings. In a PIMD simulation, the isomorphic ensemble is simulated by using the classical molecular dynamics. It provides a statistical sampling of the spatial configurations of the ring polymers.

Quantities related to electrons can be evaluated in the isomorphic ensemble. Electrons in the ensemble experience a $\tau$-dependent Kohn Sham potential $V_{\mathrm{KS}}(\tau)$ which is self-consistently determined by the ion configuration $\{\bm{R}_{a}(\tau)\}$ at the specific imaginary-time $\tau$. As a result, the single-particle Green’s function of electrons with respect to a given configuration of ring polymers, in the context of the Kohn-Sham theory, can be determined by

$$\left[-\frac{\partial}{\partial\tau}-\frac{\hat{H}(\tau)-\epsilon_{F}\hat{\mathbb{I}}}{\hbar}\right]\hat{\mathcal{G}}(\tau,\tau^{\prime})=\delta(\tau-\tau^{\prime})\hat{\mathbb{I}},$$

(4)

where we write the equation in an operator form, $\hat{H}(\tau)\equiv\hat{K}+\hat{V}_{\mathrm{KS}}(\tau)$ with $\hat{K}$ being the operators of the kinetic energy. The Green’s function is not a physical Green’s function for electrons in a liquid. To get the physical one, we need to further conduct an ensemble average:

$\displaystyle\mathcal{G}_{\mathrm{phy.}}\equiv\bar{\mathcal{G}}(\tau-\tau^{\prime})=\langle\mathcal{G}(\tau,\tau^{\prime})\rangle_{C},$

(5)

where $\langle\dots\rangle_{C}$ denotes the average over the ring-polymer configurations. The averaged (physical) Green’s function $\bar{\mathcal{G}}(\tau-\tau^{\prime})$ defines the dispersion and lifetime of a quasi-electron renormalized by the fluctuating ionic field. The quasi-electron dispersion $\tilde{\epsilon}_{\bm{k}^{\prime}}$ required by Eq. (3) can then be inferred from it.

For a given spatial configuration of ring polymers, we can determine the electron-ion scattering $\mathcal{T}$-matrix. The $\mathcal{T}$-matrix is related to the Green’s function according to the Lippmann-Schwinger equation

$\displaystyle\hat{\mathcal{T}}=\hat{\mathcal{V}}+\frac{1}{\hbar}\hat{\mathcal{V}}\hat{\bar{\mathcal{G}}}\hat{\mathcal{T}},$

(6)

where $\hat{\mathcal{V}}=\hat{V}_{\mathrm{KS}}(\tau)-\hat{\bar{\Sigma}}$ is the effective scattering potential for electrons, and $\hat{\bar{\Sigma}}$ is the self-energy inferred from $\hat{\bar{\mathcal{G}}}$.

Liu et al. establish a rigorous relation between the effective pairing interaction and the $\mathcal{T}$-matrices. One first determines the pair scattering amplitude $\Gamma$, i.e., the fluctuation of the $\mathcal{T}$-matrices:

$\displaystyle\Gamma_{11^{\prime}}=-\beta\langle|\mathcal{T}_{11^{\prime}}|^{2}\rangle_{C},$

(7)

where $\mathcal{T}_{11^{\prime}}$ denotes the matrix element of the $\mathcal{T}$-matrix, and we adopt a short-hand notation by using numbers to represent the states of electrons: $1\equiv(\bm{k},\omega_{n})$, $\bar{1}\equiv(-\bm{k},-\omega_{n})$, and similarly for $1^{\prime}$ and $\bar{1}^{\prime}$. Then, the effective pairing interaction can be obtained by solving the Bethe-Salpeter equation:

$\displaystyle W_{11^{\prime}}=\Gamma_{11^{\prime}}+\frac{1}{\hbar^{2}\beta}\sum_{2}\Gamma_{12}|\bar{\mathcal{G}}_{2}|^{2}W_{21^{\prime}}.$

(8)

The obtained $W_{11^{\prime}}$ then enters into Eq. (3), in which we express it explicitly as $W_{11^{\prime}}\equiv W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})$. We note that in general $W_{11^{\prime}}$ depends on $1$ and $1^{\prime}$ independently. For EPC, it is a good approximation to regard $W_{11^{\prime}}$ as a function of $\omega_{n}-\omega_{n^{\prime}}$ as long as $\hbar|\omega_{n(n^{\prime})}|$ is much smaller than the Fermi energy. Moreover, since liquids are isotropic and only matrix elements with both $\bm{k}$ and $\bm{k}^{\prime}$ residing on the Fermi surface are needed for evaluating Eq. (3), $W_{11^{\prime}}$ can be written as a function of $\bm{k}-\bm{k}^{\prime}$.

To develop an algorithm for estimating $T_{c}$ based on the the rigorous formalism shown in the last subsection, Liu et al. employ a number of approximations. In this subsection, we will discuss these approximations, and identify the one to be improved.

The most important approximation which makes an implementation possible is the quasi-static approximation. It is necessary because solving Eqs. (4) and (22) exactly is deemed to be impossible in practice as it requires a time resolution $\sim\hbar/\epsilon_{F}$ or a corresponding large cutoff frequency $\sim\epsilon_{F}/\hbar$. Fortunately, one can exploit the fact that ions move much slower than electrons. It means that the ionic field only has significant components for frequencies $\nu_{m}\equiv 2\pi m/\hbar\beta\lesssim\omega_{\mathrm{ph}}$, where $\omega_{\mathrm{ph}}$ is a characteristic frequency of phonons. In this case, the equation can be solved by applying the quasi-static approximation. This is to choose a large $n$ such that $\omega_{n}\gg\omega_{\mathrm{ph}}$, and solve Eq. (4) as if $\hat{H}(\tau)$ is time-independent:

$\displaystyle\hat{\mathcal{G}}_{n}(\tau)=\hbar\left[\left(\mathrm{i}\hbar\omega_{n}+\epsilon_{F}\right)\hat{\mathbb{I}}-\hat{H}(\tau)\right]^{-1},$

(9)

The ensemble average $\hat{\bar{\mathcal{G}}}_{n}=\langle\hat{\mathcal{G}}_{n}(\tau)\rangle_{C}$ will be independent of $\tau$ and a good approximation for the Fourier transform of $\bar{\mathcal{G}}(\tau-\tau^{\prime})$ at the frequency $\omega_{n}$ [1]. The $\mathcal{T}$-matrix can then be calculated in the time-domain in the same manner, and Fourier transformed back to the frequency domain. See Sec. III.2 for details. For determining the effective pairing potential, we only need to solve for a specific $n$ with $\omega_{\mathrm{ph}}\ll\omega_{n}\ll\epsilon_{F}/\hbar$. This is because $W_{11^{\prime}}\equiv W(\omega_{n}-\omega_{n^{\prime}},\bm{k}-\bm{k}^{\prime})$ is assumed to be a function of $\omega_{n}-\omega_{n}^{\prime}\equiv\nu_{m}$.

More approximations are involved in PIMD simulations. It is necessary to simulate a system in a finite size supercell, and discretize the imaginary time to a finite number of beads. In Ref. 1, it is found that the discretization truncates the high-frequency components of the density correlation function of ions and the effective pairing interaction. It results in an incorrect asymptotic behavior of $\lambda(m)$ at large $m$. To address the issue, Liu et al. develop an oversampling scheme to interpolate the density correlation function in the temporal domain, and obtain an over-sampled density correlation function. See Ref. 1 for details. To get the effective pairing interaction, they make use of the empirical relation:

$$W(\nu_{m},\bm{q})=|M(\nu_{m},\bm{q})|^{2}\chi(\nu_{m},\bm{q}),$$

(10)

where $\chi(\nu_{m},\bm{q})$ is the density correlation function of ions, and $M(\nu_{m},\bm{q})$ could be regarded as the EPC matrix element as defined in the Eliashberg theory. It is found numerically that the relation holds well for metallic hydrogen liquids, and $|M(\nu_{m},\bm{q})|^{2}$ can be obtained from a linear regression. By substituting the density correlation function in Eq. (10) with the over-sampled one, one can obtain an over-sampled effective pairing interaction. It is found that $\lambda(m)$ determined from the over-sampled effective pairing interaction has a correct asymptotic behavior at large $m$.

Lastly, Liu et al. adopt the LSA for evaluating the ionic field. This is to approximate the Fourier transform of the ionic potential as

$$V_{\mathrm{ion}}(\bm{q},\tau)=\frac{e^{2}}{\varepsilon_{\mathrm{et}}(q)q^{2}}\rho_{\mathrm{i}}(\bm{q},\tau),$$

(11)

where $\rho_{\mathrm{i}}(\bm{q},\tau)=\sum_{a}\exp[-\mathrm{i}\bm{q}\cdot\bm{R}_{a}(\tau)]$ is the Fourier transform of the density distribution of ions, and $\varepsilon_{\mathrm{et}}(q)$ is the static electron-test charge dielectric function [32] which has a form similar to that of the random phase approximation but with a local field correction factor determined by Ichimaru et al. [33]. It is derived from the jellium model which models ions as a uniform positive charge background, an approximation which obviously breaks down in core regions close to ions. While it is reasonable to expect that such an approximation could work well for metallic hydrogen in which ions (protons) have small core volumes and electrons are shared by all ions, it becomes problematic when other atoms are involved and bonding between atoms is covalent in nature.

Among the three approximations, the LSA is the one that severely limits the applicability of Liu et al.’s implementation for systems other than metallic hydrogen. This will be the focus of this paper aiming to improve the implementation. The effect of the finite simulation supercell and the finite number of beads should be assessed to make sure that calculations converge. This will also be tested. Finally, it is the quasi-static approximation which makes an implementation feasible for EPC. The approximation is valid because phonons have an energy scale much smaller than that of electrons, a feature unique to conventional EPC superconductors.

In the PAW method, electrons are described by PS wavefunctions. A transformation operator is introduced to map a PS wavefunction $\tilde{\Psi}$ to an AE wave function $\Psi$:

$\displaystyle\Psi=\hat{T}\tilde{\Psi}\equiv\tilde{\Psi}+\sum_{a}\hat{T}^{a}|\tilde{\Psi}\rangle,$

(12)

where $\hat{T}^{a}$ is an operator transforming the core state of the ion $a$ from the PS space to the AE space, and is nonzero only within an augmentation sphere surrounding the ion. The radius of the sphere is chosen such that there is no overlap between spheres associated with different ions. Inside a sphere, the wave functions are expanded in the basis of partial waves:

$\displaystyle\hat{T}^{a}|\tilde{\Psi}\rangle=|\Psi^{a}\rangle-|\tilde{\Psi}^{a}\rangle=\sum_{i}\left(|\phi_{i}^{a}\rangle-|\tilde{\phi}_{i}^{a}\rangle\right)\langle\tilde{p}_{i}^{a}|\tilde{\Psi}\rangle,\quad$

(13)

where $i$ indexes the partial wave states with angular momentum quantum numbers $(l_{i}m_{i})$ and at a given reference energy $\epsilon_{i}$. The AE partial wavefunctions $\langle\bm{r}|\phi_{i}^{a}\rangle\equiv\phi_{l_{i}\epsilon_{i}}^{a}(r)Y_{l_{i}m_{i}}(\hat{\bm{r}})$ at the radial coordinate $r$ and along the direction $\hat{\bm{r}}$ are obtained by solving a spherical AE Schrodinger equation at the reference energy in the presence of an ion at the center. The corresponding PS partial wavefunctions $\langle\bm{r}|\tilde{\phi}_{i}^{a}\rangle$ are chosen such that they are smooth and coincide with the AE partial wavefunctions outside the sphere. $\langle\tilde{p}_{i}^{a}|$ are a set of projectors which are biorthogonal to the PS partial wave states:

$\displaystyle\langle\tilde{p}_{j}^{a}|\tilde{\phi}_{i}^{a}\rangle=\delta_{i,j}.$

(14)

With a complete set of the basis, the closure relation in each sphere can be written as:

$\displaystyle\sum_{i}|\tilde{\phi}_{i}^{a}\rangle\langle\tilde{p}_{i}^{a}|=\hat{\mathbb{I}}^{a}.$

(15)

In the PS basis, the Kohn-Sham (KS) equation is transformed to

	$\displaystyle\hat{\tilde{H}}\tilde{\Psi}_{n}$	$\displaystyle=$	$\displaystyle\epsilon_{n}\hat{S}\tilde{\Psi}_{n},$		(16a)
	$\displaystyle\hat{\tilde{H}}$	$\displaystyle=$	$\displaystyle\hat{T}^{\dagger}\hat{H}\hat{T},$		(16b)
	$\displaystyle\hat{S}$	$\displaystyle=$	$\displaystyle\hat{T}^{\dagger}\hat{T},$		(16c)

where $\hat{H}$ and $\hat{\tilde{H}}$ denote the Hamiltonians in the AE and PS basis, respectively. This is a generalized eigenvalue problem, and the orthogonality condition between PS eigen-wavefunctions becomes

$\displaystyle\langle\Psi_{m}|\Psi_{n}\rangle=\langle\tilde{\Psi}_{m}|\hat{S}|\tilde{\Psi}_{n}\rangle=\delta_{m,n}.$

(17)

We note that in DFT calculations based on the PAW method, the PS Hamiltonian is directly derived from a total energy functional expressed in terms of the PS wavefunctions. It is not necessary to construct the AE Hamiltonian first and then transform it to the PS Hamiltonian by using Eq. (16b).

With the preparation, we now proceed to derive formulae for determining the AE Green’s function and the $\mathcal{T}$-matrix by using the PAW method.

We seek for an expression of the Green’s function appropriate for the PAW method. Under the quasi-static approximation, the AE Green’s function is determined by Eq. (9). However, it depends on the AE Hamiltonian $\hat{H}(\tau)$. We can transform it to a form that only depends on quantities in the PS basis. To see that, we start from the spectral representation of a general Green’s function $\hat{\mathcal{G}}(\epsilon)=\hbar(\epsilon\hat{\mathbb{I}}-\hat{H})^{-1}$:

	$\displaystyle\hat{\mathcal{G}}(\epsilon)$	$\displaystyle=$	$\displaystyle\hbar\sum_{n}\frac{|\Psi_{n}\rangle\langle\Psi_{n}|}{\epsilon-\epsilon_{n}}$		(18)
		$\displaystyle=$	$\displaystyle\hat{T}\left[\hbar\sum_{n}\frac{|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|}{\epsilon-\epsilon_{n}}\right]\hat{T}^{\dagger},$

where $|\Psi_{n}\rangle$ and $|\tilde{\Psi}_{n}\rangle$ are the AE and PS eigen-wavefunctions, respectively, and $\epsilon_{n}$ is the corresponding eigen-energy.

We can convert the expression to a more convenient form. By using the orthonormality condition Eq. (17), it is easy to verify the closure relation

$\displaystyle\sum_{n}\hat{S}|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|=\hat{\mathbb{I}}.$

(19)

We then have

	$\displaystyle\left(\epsilon\hat{S}-\hat{\tilde{H}}\right)^{-1}=$		$\displaystyle\sum_{n}\left(\epsilon\hat{S}-\hat{\tilde{H}}\right)^{-1}\hat{S}|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|$
	$\displaystyle=$		$\displaystyle\sum_{n}\frac{|\tilde{\Psi}_{n}\rangle\langle\tilde{\Psi}_{n}|}{\epsilon-\epsilon_{n}},$		(20)

where we make use of the eigen-equation Eq. (16a). By substituting the equality into Eq. (18) and setting $\epsilon=\mathrm{i}\hbar\omega_{n}+\epsilon_{F}$, we obtain the formula of determining the AE Green’s function in the PAW method:

$$\hat{\mathcal{G}}_{n}(\tau)=\hbar\hat{T}(\tau)\left[(\mathrm{i}\hbar\omega_{n}+\epsilon_{F})\hat{S}(\tau)-\hat{\tilde{H}}(\tau)\right]^{-1}\hat{T}^{\dagger}(\tau).\quad$$

(21)

We note that $\hat{\tilde{H}}$, $\hat{S}$ and $\hat{T}$ all depend on the configuration of ions, therefore are time-dependent.

With the Green’s function, the scattering $\mathcal{T}$-matrix can be determined. To avoid solving the Lippmann-Schwinger equation in the AE space, we apply an identity

$$\hat{\mathcal{T}}_{n}(\tau)=\hbar\hat{\bar{\mathcal{G}}}_{n}^{-1}\hat{\mathcal{G}}(\tau)\hat{\bar{\mathcal{G}}}_{n}^{-1}-\hbar\hat{\bar{\mathcal{G}}}_{n}^{-1}$$

(22)

to calculate the $\mathcal{T}$-matirx in time domain. It is easy to prove that the $\mathcal{T}$-matirx calculated in this way satisfies Eq. (6). We note that in a liquid, because of the spatial translation symmetry, the average Green’s function is diagonal in the plane-wave basis. As a result, the matrix inverse and multiplication can be carried out in a truncated Hilbert space spanned by the basis without introducing errors.

The matrix elements of $\mathcal{T}$-matrix in frequency domain can then be determined by applying the Fourier transform:

$\displaystyle\hat{\mathcal{T}}_{\omega_{n}+\nu_{m},\omega_{n}}\approx\frac{1}{\hbar\beta}\int_{0}^{\hbar\beta}d\tau\hat{\mathcal{T}}_{n}(\tau)\mathrm{e}^{\mathrm{i}\nu_{m}\tau}.$

(23)

The matrix elements are all that are needed for evaluating the pair scattering amplitude and the effective pairing interaction (see Sec. II.1). We note that for the purpose only the matrix elements with wave-vectors close to the Fermi surface are needed. Therefore, a large energy cutoff, which is required for recovering the AE wavefunctions from the PS wavefunctions, is not needed for our calculation.

To apply the formulae derived in the last subsection, it is necessary to determine the transformation operator $\hat{T}$ and the overlap matrix $\hat{S}$.

The transformation operator transforms a PS state to an AE state. From Eq. (12) and Eq. (13), the operator is

$$\hat{T}=\hat{\mathbb{I}}+{\sum_{j,a}}^{\prime}\left(|\phi_{j}^{a}\rangle-|\tilde{\phi}_{j}^{a}\rangle\right)\langle\tilde{p}_{j}^{a}|,$$

(24)

where the summation is over ions and partial waves, with the prime indicating explicitly that the summation is always truncated to a small number of partial waves in real calculations. For liquids, the basis states in both the AE space and the PS space are the plane waves. A matrix element in the plane wave basis is

$$\langle\psi_{\bm{k}}|\hat{T}|\tilde{\psi}_{\bm{k}^{\prime}}\rangle=\delta_{\bm{k}\bm{k}^{\prime}}\\ +{\sum_{i,a}}^{\prime}\left(\langle\psi_{\bm{k}}|\phi_{i}^{a}\rangle-\langle\psi_{\bm{k}}|\tilde{\phi}_{i}^{a}\rangle\right)\langle\tilde{p}_{i}^{a}|\tilde{\psi}_{\bm{k}^{\prime}}\rangle,$$

(25)

where $|\psi_{\bm{k}}\rangle$ and $|\tilde{\psi}_{\bm{k}^{\prime}}\rangle$ are the plane waves in the AE and PS space, respectively. The overlap matrix element between a plane wave and a partial wave can be evaluated by

$$\langle\psi_{\bm{k}}|\phi_{i}^{a}\rangle=\frac{4\pi}{\sqrt{V}}\mathrm{e}^{-\mathrm{i}\bm{k}\cdot\bm{R}_{a}}Y_{l_{i}m_{i}}(\hat{\bm{k}})\\ \times\int_{0}^{r_{c}}r^{2}\mathrm{d}r(-\mathrm{i})^{l_{i}}j_{l_{i}}(kr)\phi_{l_{i},m_{i}}^{a}(r),$$

(26)

and similarly for other matrix elements, where $Y_{lm}$ denotes the spherical harmonics, $j_{l}(kr)$ is the spherical Bessel function, $r_{c}$ denotes the cutoff radius of the augmentation sphere in defining $\hat{T}^{a}$, and $\bm{R}_{a}$ is the position of the ion $a$.

The overlap matrix $\hat{S}=\hat{T}^{\dagger}\hat{T}$ can also be determined by using the partial waves. From Eq. (24), we have

$$\hat{S}=-\hat{\mathbb{I}}+\hat{T}+\hat{T}^{\dagger}\\ +{\sum_{ij,a}}^{\prime}\Ket{\tilde{p}_{i}^{a}}\Braket{\phi_{i}^{a}-\tilde{\phi}_{i}^{a}}{\phi_{j}^{a}-\tilde{\phi}_{j}^{a}}\Bra{\tilde{p}_{j}^{a}},$$

(27)

where we assume that augmentation spheres associating with different ions do not overlap. We can rewrite it as

$$\hat{S}=\hat{S}^{\prime}+\left(\hat{T}-\hat{T}^{\prime}+\mathrm{h.c.}\right),$$

(28)

where

	$\displaystyle\hat{S}^{\prime}$	$\displaystyle=$	$\displaystyle\hat{\mathbb{I}}+{\sum_{ij,a}}^{\prime}|\tilde{p}_{i}^{a}\rangle\left(\langle\phi_{i}^{a}|\phi_{j}^{a}\rangle-\langle\tilde{\phi}_{i}^{a}|\tilde{\phi}_{j}^{a}\rangle\right)\langle\tilde{p}_{j}^{a}|,$		(29)
	$\displaystyle\hat{T}^{\prime}$	$\displaystyle=$	$\displaystyle\hat{\mathbb{I}}+{\sum_{ij,a}}^{\prime}|\tilde{p}_{j}^{a}\rangle\left(\langle\tilde{\phi}_{j}^{a}|\phi_{i}^{a}\rangle-\langle\tilde{\phi}_{j}^{a}|\tilde{\phi}_{i}^{a}\rangle\right)\langle\tilde{p}_{i}^{a}|.$		(30)

We note that $\hat{S}^{\prime}$ is the form commonly adopted in literatures. Actually, by using the closure relation Eq. (15), it is easy to see that $\hat{T}^{\prime}$ is equal to $\hat{T}$, and therefore $\hat{S}=\hat{S}^{\prime}$ if the summation is over a complete set of partial waves. However, since real calculations always use just a small set of partial waves instead of a complete set, the correction due to $\hat{T}-\hat{T}^{\prime}$ could be non-negligible. For the determination of the AE Green’s function and the $\mathcal{T}$-matrix, it is important to maintain the consistency between $\hat{T}$ and $\hat{S}$ since they both appear in Eq. (21). For the reason, we use Eq. (28) instead of Eq. (29) for calculating the overlap matrix.

We employ the GPU version of VASP [34, 35] for calculating electron structures. The radius of the augmentation spheres is set to $0.52$ Å, which may result in overlaps of the spheres in some ionic configurations but saves computation time. We test a smaller radius of $0.43$ Å and find that differences in calculation results are negligible. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [36] is used in all calculations. A Monkhorst-Pack k-point mesh of spacing not larger than $2\pi\times 0.05$ Å^-1 is used to sample the Brillouin zone, and an energy cutoff of $600$ eV is employed to expand the wavefunctions, yielding a convergence in energy better than $\sim 2\,\mathrm{meV/atom}$.

To determine the Green’s function and the scattering $\mathcal{T}$-matrix by using Eqs. (21, 22), we need to obtain the pseudo Hamiltonian $\hat{\tilde{H}}$, the transformation matrix $\hat{T}$ and the overlap matrix $\hat{S}$ from the electron structure calculations. Since VASP does not calculate or output these quantities directly, we modify it to output intermediate quantities including the local part of the effective KS potential $\tilde{v}_{\mathrm{loc}}$, the coefficients of the non-local part of the effective KS potential $(\tilde{v}^{\prime})_{ij}^{a}$ as well as the projector matrix elements $\langle\tilde{\psi}_{\bm{k}}|\tilde{p}_{i}^{a}\rangle$. With them, the matrix elements $\langle\tilde{\psi}_{\bm{k}}|\hat{\tilde{H}}|\tilde{\psi}_{\bm{k}^{\prime}}\rangle$ of the pseudo Hamiltonian can be determined by using the definition

$\displaystyle\hat{\tilde{H}}=\hat{K}+\tilde{v}_{\mathrm{loc}}+\sum_{ij,a}|\tilde{p}_{i}^{a}\rangle(\tilde{v}^{\prime})_{ij}^{a}\langle\tilde{p}_{j}^{a}|.$

(31)

The matrix elements of the transformation operators $\hat{T}$, $\hat{T}^{\prime}$ and the overlap matrix $\hat{S}^{\prime}$ can also be determined by applying Eq. (24), Eq. (30), and Eq. (29), respectively. It is then straightforward to determine $\hat{S}$ by applying Eq. (28).

We can then proceed to carry out the analysis outlined in Sec. II by following the procedure detailed in Ref. 1. Figure 1 shows the flowchart of our implementation.

To test our first-principles implementation based on PAW method, we apply it to metallic hydrogen liquids. The ion configurations from the constant-volume PIMD simulations of 200 hydrogen atoms in Ref. [1] are re-used in the current study. In this subsection, we summarize our calculation results.

Estimated $T_{c}$’s as well as relevant parameters of metallic hydrogen liquids under pressures from $0.5$ to $1.5$ TPa are summarized in Table 1. For comparison, the LSA results are also shown. We find that the results coincide well with the LSA results. It confirms Liu et al.’s prediction that the superconducting transition temperature of metallic hydrogen is higher than its predicted melting temperature [1]. We note that our prediction of $T_{c}$ is only valid for the liquid phase. $T_{c}$ for the solid phase of metallic hydrogen had been thoroughly studied in Ref. [7].

In Fig. 2, we show the comparison of the pair scattering amplitude $\Gamma_{m}(q)$ and the effective pairing interaction $W_{m}(q,\nu)$ between the two implementations for $P=700$ GPa and $T=450$ K. We find that the results yielded by the two implementations are very close. We also show the interaction parameters $\lambda(n)$ in Fig. 2(c). To determine $\lambda(n)$, we recast the matrix element $W_{11^{\prime}}$ as a function of $q=|\bm{k}-\bm{k}^{\prime}|$ for $0.8k_{F}\leq|\bm{k}|,|\bm{k}^{\prime}|\leq 1.14k_{F}$, and evaluate Eq. (3) by integrating over $q$ [5]. It is evident that the results also coincide well. It is not surprising that the LSA can work well in metallic hydrogen systems because hydrogen has no inner-core electron and the effective KS potential is dominated by its local part which can be well approximated by the LSA with an appropriate dielectric function. For more general systems like hydrides, the LSA is obviously inadequate, and our implementation based on the standard first-principles approach will be useful.

We test the convergence of our calculation by varying the size of the simulation supercell and the number of beads discretizing the imaginary time. In Ref. 1 and the calculation shown above, we use a supercell with 200 hydrogen atoms and discretize the imaginary time to 24 beads. In Fig. 3(a), we show the dependence of estimated $T_{c}$ on the number of atoms. We find that $T_{c}$ slightly decreases with the increasing simulation size, but the difference are within the error bars. In Fig. 3(b), we show the dependence of estimated $T_{c}$ on the number of beads. We find that only a relatively small number of beads ($>16$) is needed to get a converged result. The oversampling procedure discussed in Sec. II.2 is effective in eliminating the discretization errors.