A phenomenological theory of superconductor diodes

In presence of SOC and a Zeeman field, a generalized GL free energy of a SC can be written as

$$\begin{equation}F=\int \mathrm{d}\boldsymbol{q}\left\{(\alpha +\gamma {\boldsymbol{q}}^{2}+{\gamma }^{\prime }{\boldsymbol{q}}^{4}+{\eta }_{\boldsymbol{q}})\vert {\psi }_{\boldsymbol{q}}{\vert }^{2}+\frac{1}{2}(\beta +{\beta }_{2}{\boldsymbol{q}}^{2}+h{\boldsymbol{\beta }}_{1}\cdot \boldsymbol{q})\vert {\psi }_{\boldsymbol{q}}{\vert }^{4}\right\},\end{equation} \tag{ 1 }$$

where ψ _q is the order parameter in the reciprocal space. The parameters α, β and γ are conventional GL coefficients. The terms ${\eta }_{\boldsymbol{q}}=h{\sum }_{lmn}{\kappa }_{lmn}{q}_{x}^{l}{q}_{y}^{n}{q}_{z}^{m}$ (l + m + n being an odd integer) and h β ₁ ⋅ q , originating from SOC and external Zeeman field h, break both inversion $(\mathcal{P})$ and time-reversal $(\mathcal{T})$ symmetries and lead to magnetochiral anisotropy. It is assumed that h ≪ T_c (we omit the Bohr magneton μ_B, the Boltzmann constant k_B and the reduced Planck constant ℏ throughout this paper). The higher order terms corresponding to γ' and β₂ are also included for reasons that will be clear later.

The structure of the coupling constants κ_lmn and β ₁ is directly related to the symmetries of the system and thus depends on the form of the SOC. For example, in a system with continuous rotational symmetry ${\mathcal{C}}_{\infty }$ (rotation axis along $\hat{\boldsymbol{z}}$), the group representation leads to ${\eta }_{\boldsymbol{q}}=({a}_{0}+{a}_{2}\vert \boldsymbol{q}{\vert }^{2})(\boldsymbol{h}\cdot \boldsymbol{q})+({b}_{0}+{b}_{2}\vert \boldsymbol{q}{\vert }^{2})(\boldsymbol{h}\times \boldsymbol{q})\cdot \hat{\boldsymbol{z}}$, up to the linear order in h and the third order in q . The h ⋅ q term breaks all mirror symmetries, while $(\boldsymbol{h}\times \boldsymbol{q})\cdot \hat{\boldsymbol{z}}$ breaks ${\mathcal{M}}_{z}$, the mirror symmetry in the z-direction. When ${\mathcal{C}}_{\infty }$ is reduced to a discrete rotation ${\mathcal{C}}_{n}$, a form of ${\eta }_{\boldsymbol{q}}=({c}_{1}\vert {\boldsymbol{q}}_{{\Vert}}\vert +{c}_{3}\vert {\boldsymbol{q}}_{{\Vert}}{\vert }^{3}){h}_{z}\,\mathrm{cos}(n{\theta }_{{\boldsymbol{q}}_{{\Vert}}})$ is allowed, which preserves ${\mathcal{M}}_{z}$.

In general, even if both $\mathcal{P}$ and $\mathcal{T}$ are broken, nonreciprocal effects are not necessarily expected. To see that, let us define the inversion operators of each dimension, ${\mathcal{P}}_{x},{\mathcal{P}}_{y}$ and ${\mathcal{P}}_{z}$, so that the corresponding symmetry invariance requires ${\mathcal{P}}_{x/y/z}H({k}_{x/y/z}){\mathcal{P}}_{x/y/z}^{-1}=H(-{k}_{x/y/z})$ respectively, where H( k ) is the Hamiltonian of the system. Obviously $\mathcal{P}={\mathcal{P}}_{x}{\mathcal{P}}_{y}{\mathcal{P}}_{z}$ is broken. However, since breaking ${\mathcal{P}}_{x}$ is necessary for any possible difference between the currents along the ±x directions, a term such as $h{q}_{x}^{2}{q}_{y}$, although breaking $\mathcal{P}$ and $\mathcal{T}$, would not cause such a nonreciprocity. This means that the direction of the magnetic field to induce nonreciprocal effects (in ±x-directions) needs to be determined by symmetries—it should break all possible ${\mathcal{P}}_{x}$ of the Hamiltonian. This will be illustrated with the example discussed in the later part of this paper.

Consider the case where the magnitude of the SC order parameter is uniform and it only varies in its phase along the x-direction, i.e. ψ( r ) = |ψ|e^iϕ(x). This assumption is valid as long as we are dealing with SCs of thicknesses much less than the coherence length. In this simplified case, the free energy becomes

$$\begin{equation}F=\int \mathrm{d}q\left\{[\alpha +\gamma {q}^{2}+{\gamma }^{\prime }{q}^{4}+q(h{\kappa }_{1}+h{\kappa }_{3}{q}^{2})]\vert {\psi }_{q}{\vert }^{2}+\frac{1}{2}(\beta +h{\beta }_{1}q+{\beta }_{2}{q}^{2})\vert {\psi }_{q}{\vert }^{4}\right\},\end{equation} \tag{ 2 }$$

with q = ∂_x ϕ(x). The order parameter may be multi-component in general. In that case, ψ denotes a certain linear combination of these components which minimizes the energy, and the internal structure of the order parameter does not affect our discussion. Also note that the magnetic field appears as a scalar since it is assumed in the proper direction to be determined by the specific form of the SOC in a concrete model.

The terms in equation (2) do not affect the coherence length since the spatial variation is assumed in the phase ϕ only. As for the κ₁ term, there should be a linear derivative term with respect to the magnitude |ψ| correspondingly. However, it vanishes after integrated over space. Higher-order derivative terms of |ψ| could modify the coherence length, but which is a small effect when we consider a weak magnetic field.

The supercurrent along the x-direction is

$$\begin{equation}I=-2e\left[(2\gamma q+4{\gamma }^{\prime }{q}^{3}+h{\kappa }_{1}+3h{\kappa }_{3}{q}^{2})\vert {\psi }_{q}{\vert }^{2}+\frac{1}{2}(h{\beta }_{1}+2{\beta }_{2}q)\vert {\psi }_{q}{\vert }^{4}\right],\end{equation} \tag{ 3 }$$

|ψ_q | is obtained by minimizing F, which leads to $\vert {\psi }_{q}{\vert }^{2}=\frac{\vert \alpha \vert }{\beta }f(\tilde{q}),$ with

$$\begin{equation}f(\tilde{q})=\frac{1-{\tilde{q}}^{2}-{\tilde{\gamma }}^{\prime }{\tilde{q}}^{4}-{\tilde{\kappa }}_{1}\tilde{q}-{\tilde{\kappa }}_{3}{\tilde{q}}^{3}}{1+{\tilde{\beta }}_{1}\tilde{q}+{\tilde{\beta }}_{2}{\tilde{q}}^{2}}.\end{equation} \tag{ 4 }$$

We have introduced the dimensionless variables $\tilde{q}\equiv q\sqrt{\gamma /\vert \alpha \vert }$, ${\tilde{\gamma }}^{\prime }\equiv {\gamma }^{\prime }\vert \alpha \vert /{\gamma }^{2}$, ${\tilde{\beta }}_{1}\equiv h{\beta }_{1}{\beta }^{-1}\sqrt{\vert \alpha \vert /\gamma }$, ${\tilde{\beta }}_{2}\equiv {\beta }_{2}{\beta }^{-1}\vert \alpha \vert /\gamma$, and ${\tilde{\kappa }}_{n}\equiv h{\kappa }_{n}\vert \alpha {\vert }^{n/2-1}{\gamma }^{-n/2}$. Substitution of |ψ_q | into equation (3) yields

$$\begin{equation}I=-2e\frac{\sqrt{\vert \alpha {\vert }^{3}\gamma }}{\beta }\left[(2\tilde{q}+4{\tilde{\gamma }}^{\prime }{\tilde{q}}^{3}+{\tilde{\kappa }}_{1}+3{\tilde{\kappa }}_{3}{\tilde{q}}^{2})f(\tilde{q})+\frac{1}{2}({\tilde{\beta }}_{1}+2{\tilde{\beta }}_{2}\tilde{q})f{(\tilde{q})}^{2}\right].\end{equation} \tag{ 5 }$$

The current I as a function of $\tilde{q}$ has a maximum I_c+ and a minimum −I_c−, I_c± being the critical currents along the positive and negative directions respectively. It should be noted that the supercurrent I is nonzero when q vanishes, as can be seen in equation (5). This indicates that the ground state is not with zero q. Instead, the value of ground-state q is determined by I = 0 which leads to q = q₀ ≠ 0. This kind of finite-momentum pairing usually accompanies the SC diode effect. However, while a q-linear term in the free energy, i.e. κ₁ ≠ 0 in equation (2), is enough to induce finite-momentum pairing, the SC diode effect requires more, as will be clear soon.

When h = 0, the maximum and minimum are readily found at ${\tilde{q}}_{\pm }\to \pm 1/\sqrt{3}$, and thus ${q}_{\pm }\to \sqrt{\vert \alpha \vert /3\gamma }\sim \sqrt{{\epsilon}}$, where ≡ 1 − T/T_c. With nonzero but small h, the variables ${\tilde{\gamma }}^{\prime },{\tilde{\kappa }}_{3}$ and ${\tilde{\beta }}_{1,2}$ are much smaller than unity and the solutions ${\tilde{q}}_{\pm }$ are only slightly shifted. The extrema can be obtained by expansion, which leads to the critical currents (up to the first order in $h\sqrt{{\epsilon}}$)

$$\begin{equation}{I}_{\text{c}\pm }\approx \frac{8e}{3\sqrt{3}}\frac{\sqrt{\vert \alpha {\vert }^{3}\gamma }}{\beta }\left(1\pm \frac{Q}{2}\right),\end{equation} \tag{ 6 }$$

where we defined the diode quality parameter

$$\begin{equation}Q\equiv \frac{{I}_{\text{c}+}-{I}_{\text{c}-}}{({I}_{\text{c}+}+{I}_{\text{c}-})/2}=\frac{1}{2\sqrt{3}}(2{\tilde{\beta }}_{1}+4{\tilde{\kappa }}_{1}{\tilde{\beta }}_{2}-4{\tilde{\kappa }}_{3}+5{\tilde{\kappa }}_{1}{\tilde{\gamma }}^{\prime }).\end{equation} \tag{ 7 }$$

To see whether a given term in equation (2) is important to the SC diode effect up to the lowest orders in h and , one may count the exponent of in it. Since ${\tilde{\kappa }}_{1}\sim {{\epsilon}}^{-1/2},\,\,{\tilde{\kappa }}_{3}\sim {{\epsilon}}^{1/2},\,\,{\tilde{\gamma }}^{\prime }\sim {\epsilon},\,\,{\tilde{\beta }}_{1}\sim {{\epsilon}}^{1/2}$ and ${\tilde{\beta }}_{2}\sim {\epsilon}$, all the terms in equation (7) are linear in $h\sqrt{{\epsilon}}$. On the other hand, one can show that all terms contributing to Q up to the order $\sim h\sqrt{{\epsilon}}$ have been included in equations (1) and (2). Thus, all the terms in equations (1) and (2) are important while other higher order terms can be neglected. From equation (7), it is clear that the q-linear term in the kinetic energy of Cooper pairs, i.e., the ${\tilde{\kappa }}_{1}$ term in η _q , would not change the critical current alone, because it only shifts the positions of the maximum and the minimum of equation (5) while keeping their values unchanged. (For this reason, the divergence of ${\tilde{\kappa }}_{1}$ at → 0 does not cause problems.)

For example of SC diode effects, let us consider two-dimensional SCs with Rashba SOC. The normal Hamiltonian can be written as

$$\begin{equation}{H}^{\text{R}}(\boldsymbol{k})=\left(\frac{{k}^{2}}{2m}-\mu \right){\sigma }_{0}+{\lambda }_{\text{R}}({k}_{x}{\sigma }_{y}-{k}_{y}{\sigma }_{x})+\boldsymbol{h}\cdot \boldsymbol{\sigma },\end{equation} \tag{ 8 }$$

where λ_R is the Rashba SOC strength, k = (k_x , k_y ) is the electron wave vector, h is the magnetic field, μ is the chemical potential, and σ_x,y are Pauli matrices. Two x-inverting symmetries, ${\mathcal{P}}_{x1}={\sigma }_{x}$ and ${\mathcal{P}}_{x2}={\sigma }_{z}$, are preserved when h = 0. Their effects on magnetic fields in different directions are shown in table 1. To break a symmetry, the Zeeman term must be odd under the symmetry operation. Table 1 shows that only a Zeeman field along the y-direction breaks both ${\mathcal{P}}_{x1}$ and ${\mathcal{P}}_{x2}$. According to our previous symmetry analysis, a nonreciprocity in the ±x-directions is expected only if h_y ≠ 0.

Table 1. Symmetry operations on the Zeeman fields along the three directions in Rashba systems. A plus (minus) sign means that the Zeeman term is even (odd) under the symmetry operation.

	${\mathcal{P}}_{x1}$	${\mathcal{P}}_{x2}$
hx σx	+	−
hy σy	−	−
hz σz	−	+

The Rashba SOC and the Zeeman field result in the following term of the GL free energy (up to the linear order in h ),

$$\begin{equation}\delta {F}^{\text{R}}=\int \mathrm{d}\boldsymbol{q}({q}_{x}{h}_{y}-{q}_{y}{h}_{x})\left({\kappa }_{1}^{\text{R}}+{\kappa }_{3}^{\text{R}}\vert \boldsymbol{q}{\vert }^{2}+\frac{{\beta }_{1}^{\text{R}}}{2}\vert {\psi }_{\boldsymbol{q}}{\vert }^{2}\right)\vert {\psi }_{\boldsymbol{q}}{\vert }^{2}.\end{equation} \tag{ 9 }$$

Thus, if a magnetic field along the y-direction, h = (0, h_y , 0), is applied, the critical currents along the ±x-direction will be different, as previously obtained in equations (6) and (7) and consistent with the symmetry analysis.

Assuming $\vert \boldsymbol{h}\vert \ll {T}_{\text{c}}\ll {E}_{\text{R}}=\frac{1}{2}m{\lambda }_{\text{R}}^{2}$ and treating the problem in the band basis, one may neglect the inter-band terms and consider only the intra-band pairing Δ. With this simplification, the GL coefficients in equation (2) can be obtained and the resulting Q-parameter is

$$\begin{equation}{Q}^{\text{R}}=\frac{2.7{\lambda }_{\text{R}}}{\vert {\lambda }_{\text{R}}\vert }\frac{h\sqrt{{\epsilon}}}{{T}_{\text{c}}}\times \left\{\begin{array}{cc}\hfill {(1+\tilde{\mu })}^{-1/2},\hfill & \hfill (\tilde{\mu } > 0)\hfill \\ \hfill \frac{8}{7}+\frac{16}{21}\tilde{\mu }+{(1+\tilde{\mu })}^{1/2},\hfill & \hfill (-1< \tilde{\mu }< 0)\hfill \end{array}\right.,\end{equation} \tag{ 10 }$$

where $\tilde{\mu }\equiv \mu /{E}_{\text{R}}$. (Note that $\mu /{E}_{\text{R}}={(2\mu /{\lambda }_{\text{R}}{k}_{\text{F}})}^{2}$ where αk_F is the Rashba splitting energy at the Fermi surface.) Q^R as a function of the chemical potential μ is shown in figure 1. The parameter Q^R has its maximum at the band crossing point μ = 0 and decreases as the Fermi level moves away either towards the band edge μ = −E_R or towards the limit μ ≫ E_R. At μ = 0, there exists a discontinuity due to the flip of the helicity of the spin-momentum locking. Note that the discontinuity appears also because we took the limit T_c/E_R → 0 and neglect the inter-band pairing. The calculation is done in the band basis assuming a constant pairing breaking energy near the Fermi surface, which is true when both μ + E_R ≫ Δ and |μ| ≫ Δ are satisfied. Near μ = 0, moreover, the smallness of the Fermi wave vector k_F, compared to the Cooper pair wave vector | q |, invalidates the series expansion over q /k_F for the GL theory. Thus, the discontinuity shall be smoothed out when T_c/E_R is not infinitesimal. And our GL theory calculations do not apply near μ = 0 or μ = −E_R.

Zoom In Zoom Out Reset image size

The quality parameter Q^R may also be obtained using a self-consistent Bogoliubov–de Gennes mean-field Hamiltonian

$$\begin{equation}{\hat{H}}_{\text{BdG}}=\sum\limits _{\boldsymbol{k}}{H}_{ij}^{\text{R}}(\boldsymbol{k}){\psi }_{i}^{{\dagger}}(\boldsymbol{k}){\psi }_{i}(\boldsymbol{k})+{\Delta}{\psi }_{{\uparrow}}^{{\dagger}}\left(\frac{\boldsymbol{q}}{2}+\boldsymbol{k}\right){\psi }_{{\downarrow}}^{{\dagger}}\left(\frac{\boldsymbol{q}}{2}-\boldsymbol{k}\right)+\text{h.c.},\end{equation} \tag{ 11 }$$

where i, j =  ↑↓ are matrix indexes in the spin space. This method applies to wider parameter regions although it is feasible only numerically. Note that the pairing gap Δ depends on the wave vector q since it is determined by minimizing the free energy $F(\boldsymbol{q})=-T\,{\sum }_{n,\boldsymbol{k}}\mathrm{ln}(1+{e}^{-{{\epsilon}}_{n}/T})$, where _n are the eigenvalues of ${\hat{H}}_{\text{BdG}}$. For a given q , the corresponding supercurrent is I_x ( q ) = 2e∂F/∂q_x . The critical currents I_c+ and I_c− are obtained by finding the maximums of I_x ( q ) and −I_x ( q ) respectively. The diode quality parameter Q^R, defined in equation (7), as a function of μ is shown in figure 2, which has qualitatively the same features as those of figure 1 obtained with the generalized GL method. The discontinuity at μ = 0 becomes smooth since T_c/E_R is not so small. In the large μ limit, Q^R ∼ μ^1/2, as shown by the log-scale plot in the inset of figure 2.

Zoom In Zoom Out Reset image size

The temperature dependence of Q^R is shown in figure 3. It gradually increases as T is lowered from T_c, consistent with the prediction, ${Q}^{\text{R}}\sim \sqrt{{T}_{\text{c}}-T}$, by the generalized GL theory. However, as T further decreases, Q^R starts to increase dramatically. (Results with temperatures near zero cannot be obtained here due to a numerical convergence problem.)

Zoom In Zoom Out Reset image size

Both analytical and numerical calculations show that the SC diode effect in two-dimensional Rashba systems reaches its maximum at the band crossing point. This suggests that stronger experimental signals may be achieved by tuning the chemical potential closer to zero by, for example, gating, as well as by increasing the magnetic field or decreasing the temperature.

The critical currents, or the extrema of equation (5), are calculated perturbatively. We first find the zero-order solutions by assuming ${\tilde{\kappa }}_{1}={\tilde{\kappa }}_{3}={\tilde{\gamma }}^{\prime }={\tilde{\beta }}_{1}={\tilde{\beta }}_{2}=0$. They are found at

$$\begin{equation}{\tilde{q}}_{0\pm }=\pm 1/\sqrt{3}.\end{equation} \tag{ 12 }$$

With nonzero ${\tilde{\kappa }}_{1},{\tilde{\kappa }}_{3},{\tilde{\gamma }}^{\prime },{\tilde{\beta }}_{1}$ and ${\tilde{\beta }}_{2}$, $I(\tilde{q})$ can be expanded around ${\tilde{q}}_{0\pm }$ up to the order of ${\tilde{q}}^{2}$. Keeping only the lowest-order terms in ${\tilde{\kappa }}_{1},{\tilde{\kappa }}_{3},{\tilde{\gamma }}^{\prime },{\tilde{\beta }}_{1},{\tilde{\beta }}_{2}$, one find the extrema of $I(\tilde{q})$ near ${\tilde{q}}_{0\pm }$ as given in equations (6) and (7).

The GL coefficients are obtained in a standard way by applying perturbation method to the Bogoliubov–de Gennes mean-field Hamiltonian,

$$\begin{equation}{H}_{\mathrm{B}\mathrm{d}\mathrm{G}}(\boldsymbol{k})=\left(\begin{matrix}\hfill {H}^{\text{R}}(\boldsymbol{q}/2+\boldsymbol{k})\hfill & \hfill {\hat{{\Delta}}}_{q}\hfill \\ \hfill {\hat{{\Delta}}}_{q}^{{\dagger}}\hfill & \hfill -{H}^{\text{R}}{(\boldsymbol{q}/2-\boldsymbol{k})}^{\ast }\hfill \end{matrix}\right),\end{equation} \tag{ 13 }$$

where H^R( k ) is the normal Hamiltonian defined in equation (8) and ${\hat{{\Delta}}}_{q}={{\Delta}}_{q}i{\sigma }_{y}$ is the SC pairing term. The free energy (q-integrand) up to ${{\Delta}}_{q}^{2}$ is calculated by

$$\begin{equation}{f}^{(2)}(q)=\frac{\vert {{\Delta}}_{q}{\vert }^{2}}{g}-\frac{T}{4{\pi }^{2}}\sum\limits _{\boldsymbol{k},n}\mathrm{Tr}[G({\omega }_{n},\boldsymbol{q}/2+\boldsymbol{k}){{\Delta}}_{q}{G}^{\text{T}}(-{\omega }_{n},\boldsymbol{q}/2-\boldsymbol{k}){{\Delta}}_{q}^{{\dagger}}],\end{equation} \tag{ 14 }$$

and the fourth order term is

$$\begin{equation}{f}^{(4)}(q)=-\frac{T}{8{\pi }^{2}}\sum\limits _{\boldsymbol{k},n}\mathrm{Tr}[{(G({\omega }_{n},\boldsymbol{q}/2+\boldsymbol{k}){{\Delta}}_{q}{G}^{\text{T}}(-{\omega }_{n},\boldsymbol{q}/2-\boldsymbol{k}){{\Delta}}_{q}^{{\dagger}})}^{2}].\end{equation} \tag{ 15 }$$

g > 0 is the on-site attractive interaction strength, which is to be determined self-consistently for a given T_c.

In a Rashba SC, we neglect the inter-band pairing and get

$$\begin{align}\hfill {f}^{(2)}(q)& =\frac{\vert {{\Delta}}_{q}{\vert }^{2}}{g}-\frac{T\vert {{\Delta}}_{q}{\vert }^{2}}{4{\pi }^{2}}\sum\limits _{n,\boldsymbol{k},\pm }\frac{1}{i{\omega }_{n}-{\xi }_{\pm }(\boldsymbol{k})-{\delta }_{\pm }(\boldsymbol{k})}\frac{1}{i{\omega }_{n}+{\xi }_{\pm }(\boldsymbol{k})-{\delta }_{\pm }(\boldsymbol{k})},\hfill \\ \hfill & =\frac{\vert {{\Delta}}_{q}{\vert }^{2}}{g}-\frac{T\vert {{\Delta}}_{q}{\vert }^{2}}{4{\pi }^{2}}\sum\limits _{n,\pm }\int \mathrm{d}{\xi }_{\pm }{\nu }_{\pm }\frac{1}{i{\omega }_{n}-{\xi }_{\pm }-{\delta }_{\pm }}\frac{1}{i{\omega }_{n}+{\xi }_{\pm }-{\delta }_{\pm }},\hfill \end{align} \tag{ 16 }$$

where ${\xi }_{\pm }(\boldsymbol{k})=\frac{1}{2}({\xi }_{\pm }^{e}-{\xi }_{\pm }^{h})$ and ${\delta }_{\pm }(\boldsymbol{k})=\frac{1}{2}({\xi }_{\pm }^{e}+{\xi }_{\pm }^{h})$, with ${\xi }_{\pm }^{e}(\boldsymbol{k})$ and ${\xi }_{\pm }^{h}(\boldsymbol{k})$ being the eigen-values of H^R( q /2 + k ) and −H^R( q /2 − k )* respectively. The pair-breaking energies δ_±( k ) contains contributions from both the Cooper pair wave vector q and the Zeeman field $\boldsymbol{h}={h}_{y}\hat{\boldsymbol{y}}$. In the second line, we changed the summation over k into the integral over the energy ξ_± by introducing the densities of states ν_±. Assuming Δ_q is small, δ_± and ν_± may be treated as a constants, and thus

$$\begin{align}\hfill {f}^{(2)}(q)& =\frac{\vert {{\Delta}}_{q}{\vert }^{2}}{g}-\frac{T\vert {{\Delta}}_{q}{\vert }^{2}}{4{\pi }^{2}}\sum\limits _{\pm }\mathfrak{R}\left[\sum\limits _{n=0}^{\infty }\frac{2\pi }{{\omega }_{n}+i\delta }\right]\hfill \\ \hfill & =\frac{\vert {{\Delta}}_{q}{\vert }^{2}}{g}-\frac{T\vert {{\Delta}}_{q}{\vert }^{2}}{4{\pi }^{2}}\sum\limits _{\pm }\mathfrak{R}\left[\sum\limits _{n=0}^{\infty }\frac{2\pi }{{\omega }_{n}}\sum\limits _{l=0}^{\infty }{\left(\frac{-i{\delta }_{\pm }}{{\omega }_{n}}\right)}^{l}\right]\qquad \hfill \end{align} \tag{ 17 }$$

$$\begin{align}\hfill & \qquad =\frac{\vert {{\Delta}}_{q}{\vert }^{2}}{g}-\frac{T\vert {{\Delta}}_{q}{\vert }^{2}}{4{\pi }^{2}}\sum\limits _{\pm }\sum\limits _{n=0}^{\infty }\frac{2\pi }{{\omega }_{n}}\sum\limits _{l=0}^{\infty }{(-1)}^{l}{\left(\frac{{\delta }_{\pm }}{{\omega }_{n}}\right)}^{2l}\hfill \\ \hfill & \qquad =\frac{\vert {{\Delta}}_{q}{\vert }^{2}}{g}-\frac{T\vert {{\Delta}}_{q}{\vert }^{2}}{4{\pi }^{2}}\sum\limits _{\pm }\left[\frac{{T}_{\text{c}}-T}{{T}_{\text{c}}^{2}}-{\delta }_{\pm }^{2}\frac{7\zeta (3)}{4{\pi }^{2}{T}^{3}}+{\delta }_{\pm }^{4}\frac{31\zeta (5)}{16{\pi }^{4}{T}^{5}}\right]+O({\delta }_{\pm }^{5}).\hfill \end{align} \tag{ 18 }$$

The calculation of f  ⁽⁴⁾(q) follows a similar procedure.

Keeping the terms in f  ⁽²⁾(q) and f  ⁽⁴⁾(q) up to the fourth order in q and to the first order in h_y , we obtain the GL coefficients as follows.

When μ > 0,

$$\begin{equation}{\alpha }_{+}^{\text{R}}=\frac{m(T-{T}_{\text{c}})}{\pi {T}_{\text{c}}}, \end{equation} \tag{ 19 }$$

$$\begin{equation}{\gamma }_{+}^{\text{R}}=\frac{7\zeta (3)}{16{\pi }^{3}}\frac{{E}_{\text{R}}+\mu }{{T}_{\text{c}}^{2}}, \end{equation} \tag{ 20 }$$

$$\begin{equation}{\kappa }_{1+}^{\text{R}}=-\frac{7\zeta (3)}{4{\pi }^{3}}\frac{{h}_{y}m{\lambda }_{\text{R}}\left(1+\frac{1}{2}\mu /{E}_{\text{R}}\right)}{{T}_{\text{c}}^{2}},\qquad \end{equation} \tag{ 21 }$$

$$\begin{equation}{\kappa }_{3+}^{\text{R}}=\frac{93\zeta (5)}{64{\pi }^{5}}\frac{{h}_{y}{\lambda }_{\text{R}}{E}_{\text{R}}}{{T}_{\text{c}}^{4}}\left(1+\frac{1}{2}\mu /{E}_{\text{R}}\right)(1+\mu /{E}_{\text{R}}),\!\!\end{equation} \tag{ 22 }$$

$$\begin{equation}{\gamma }_{+}^{\prime \text{R}}=-\frac{93\zeta (5)}{512{\pi }^{5}}\frac{{({E}_{\text{R}}+\mu )}^{2}}{m{T}_{\text{c}}^{4}}, \end{equation} \tag{ 23 }$$

$$\begin{equation}{\beta }_{+}^{\text{R}}=\frac{7\zeta (3)}{16{\pi }^{3}}\frac{m}{{T}_{\text{c}}^{2}}, \end{equation} \tag{ 24 }$$

$$\begin{equation}{\beta }_{1+}^{\text{R}}=\frac{93\zeta (5){h}_{y}}{32{\pi }^{5}}\frac{2{E}_{\text{R}}+\mu }{{T}_{\text{c}}^{4}{\lambda }_{\text{R}}}, \end{equation} \tag{ 25 }$$

$$\begin{equation}{\beta }_{2+}^{\text{R}}=-\frac{93\zeta (5)}{128{\pi }^{5}}\frac{{E}_{\text{R}}+\mu }{{T}_{\text{c}}^{4}}. \end{equation} \tag{ 26 }$$

When μ < 0,

$$\begin{equation}{\alpha }_{-}^{\text{R}}=\frac{{\alpha }_{+}^{\text{R}}}{\sqrt{1+\mu /{E}_{\text{R}}}}, \end{equation} \tag{ 27 }$$

$$\begin{equation}{\gamma }_{-}^{\text{R}}=\frac{7\zeta (3)}{16{\pi }^{3}}\frac{\sqrt{{E}_{\text{R}}}\sqrt{{E}_{\text{R}}+\mu }}{{T}_{\text{c}}^{2}}, \end{equation} \tag{ 28 }$$

$$\begin{equation}{\kappa }_{1-}^{\text{R}}=-\frac{7\zeta (3)}{4{\pi }^{3}}\frac{{h}_{y}m{\lambda }_{\text{R}}\sqrt{1+\mu /{E}_{\text{R}}}}{{T}_{\text{c}}^{2}},\qquad \end{equation} \tag{ 29 }$$

$$\begin{equation}{\kappa }_{3-}^{\text{R}}=\frac{93\zeta (5)}{64{\pi }^{5}}\frac{{h}_{y}{\lambda }_{\text{R}}{E}_{\text{R}}}{{T}_{\text{c}}^{4}}{(1+\mu /{E}_{\text{R}})}^{3/2}, \end{equation} \tag{ 30 }$$

$$\begin{equation}{\gamma }_{-}^{\prime \text{R}}=\frac{{\gamma }_{+}^{\prime \text{R}}}{\sqrt{1+\mu /{E}_{\text{R}}}}, \end{equation} \tag{ 31 }$$

$$\begin{equation}{\beta }_{-}^{\text{R}}=\frac{{\beta }_{+}^{\text{R}}}{\sqrt{1+\mu /{E}_{\text{R}}}}, \end{equation} \tag{ 32 }$$

$$\begin{equation}{\beta }_{1-}^{\text{R}}=\frac{93\zeta (5){h}_{y}}{32{\pi }^{5}}\frac{{(\sqrt{{E}_{\text{R}}}+\sqrt{{E}_{\text{R}}+\mu })}^{2}}{{T}_{\text{c}}^{4}{\lambda }_{\text{R}}},\qquad \end{equation} \tag{ 33 }$$

$$\begin{equation}{\beta }_{2-}^{\text{R}}=-\frac{93\zeta (5)}{128{\pi }^{5}}\frac{\sqrt{{E}_{\text{R}}+\mu }\left(\sqrt{{E}_{\text{R}}+\mu }+\sqrt{{E}_{\text{R}}}\right)}{{T}_{\text{c}}^{4}}.\quad \end{equation} \tag{ 34 }$$

ζ(x) is the Riemann zeta function and ${E}_{\text{R}}=m{\lambda }_{\text{R}}^{2}/2$.

The expression of Q^R in equation (10) is obtained by substituting the above results into equation (7). The contributions of the four terms are of the same order of magnitude (see supplementary (https://stacks.iop.org/NJP/24/053014/mmedia)) and thus none of them can be neglected. It is also found that the discontinuity of Q^R at μ = 0 comes from β₁ and β₂, i.e. from the terms quartic in Δ_q .