A Time-Dependent Ginzburg-Landau Framework for Sample-Specific Simulation of Superconductors for SRF Applications

Ginzburg-Landau (GL) theory is one of the oldest theories of superconductivity, and it remains relevant today owing to its relative simplicity and direct physical insights into the electrodynamic response of superconductors under static applied fields and currents [52]. The time-dependent Ginzburg-Landau (TDGL) equations were originally proposed by Schmid [53] in 1966 and Gor’kov and Eliashberg [54] derived them rigorously from BCS theory later in 1968. The TDGL equations (in Gaussian units) are given by:

	$\displaystyle\Gamma\left(\frac{\partial\psi}{\partial t}+\frac{ie_{s}\phi}{\hbar}\psi\right)+\frac{1}{2m_{s}}\left(-i\hbar\nabla-\frac{e_{s}}{c}\mathbf{A}\right)^{2}\psi+\alpha\psi+\beta|\psi|^{2}\psi=0$		(1)
	$\displaystyle\frac{4\pi\sigma_{n}}{c}\left(\frac{1}{c}\frac{\partial\mathbf{A}}{\partial t}+\nabla\phi\right)+\nabla\times\nabla\times\mathbf{A}-\frac{2\pi ie_{s}\hbar}{m_{s}c}\left(\psi^{*}\nabla\psi-\psi\nabla\psi^{*}\right)+\frac{4\pi e^{2}_{s}}{m_{s}c^{2}}|\psi|^{2}\mathbf{A}=0.$		(2)

These equations are solved for the complex superconducting order parameter, $\psi$, and the magnetic vector potential, $\mathbf{A}$. The magnitude squared of $\psi$ is proportional to the density of superconducting electrons. The parameters $\alpha$ and $\beta$ are phenomenological, and were originally introduced as coefficients of the series expansion of the Ginzburg-Landau free energy. Additionally, $\phi$ is the scalar potential; $\sigma_{n}$ is the normal electron conductivity; $\Gamma$ is the phenomenological relaxation rate of $\psi$. Furthermore, $e_{s}=2e$ and $m_{s}=2m_{e}$ represent the total charge and total effective mass of a Cooper pair, respectively. The TDGL equations are subject to boundary conditions

	$\displaystyle\left(i\hbar\nabla\psi+\frac{e_{s}}{c}\mathbf{A}\psi\right)\cdot\mathbf{n}=0$		(3)
	$\displaystyle\left(\nabla\times\mathbf{A}\right)\times\mathbf{n}=\mathbf{H}_{a}\times\mathbf{n}$		(4)
	$\displaystyle\left(\nabla\phi+\frac{1}{c}\frac{\partial\mathbf{A}}{\partial t}\right)\cdot\mathbf{n}=0,$		(5)

where $\mathbf{n}$ is the outward normal vector to the boundary surface and $\mathbf{H}_{a}$ is the applied magnetic field. Eq. 3 ensures no current flows out of the superconducting domain, and noting that $E=-\nabla\phi-\frac{1}{c}\frac{\partial\mathbf{A}}{\partial t}$, Eqs. 4 and 5 are electromagnetic interface conditions with an applied magnetic field.

The parameters $\alpha$, $\beta$, and $\Gamma$ were originally introduced into the theory as phenomenological, temperature-dependent constants [55]. It is worth noting that $\alpha<0$ corresponds to the superconducting state whereas $\alpha\geq 0$ corresponds to the normal state; $\beta$ is strictly positive regardless of the system’s state. The TDGL equations can also be derived from microscopic theory using the time-dependent Gor’Kov equations [54]. A useful consequence of this derivation is that it allows the TDGL parameters to be directly related to experimentally observable properties of the material in question. The material dependencies are given by Ref. 56:

	$\displaystyle\alpha(\nu(0),T_{c},T)$	$\displaystyle=-\nu(0)\left(\frac{1-\frac{T^{2}}{T_{c}^{2}}}{1+\frac{T^{2}}{T_{c}^{2}}}\right)$		(6)
		$\displaystyle\approx-\nu(0)\left(1-\frac{T}{T_{c}}\right)$
	$\displaystyle\beta(\nu(0),T_{c},T)$	$\displaystyle=\frac{7\zeta(3)\nu(0)}{8\pi^{2}T_{c}^{2}}\left(\frac{1}{1+\frac{T^{2}}{T_{c}^{2}}}\right)^{2}$		(7)
		$\displaystyle\approx\frac{7\zeta(3)\nu(0)}{8\pi^{2}T_{c}^{2}}$
	$\displaystyle\Gamma(\nu(0),T_{c})$	$\displaystyle=\frac{\nu(0)\pi\hbar}{8T_{c}},$		(8)

where $\nu(0)$ is the density of states at the Fermi-level, $T_{c}$ is the critical temperature, $T$ is the temperature, and $\zeta(x)$ is the Riemann zeta function.

A major contribution of this paper is to introduce a framework for modeling sample-specific features of superconducting materials by connecting ab initio calculations of the material’s properties to experimental characterizations of the material. Eqs. 6-8 determine how the parameters of TDGL depend on the underlying material properties. We allow these parameters to vary spatially to capture the sample-specific features observed from experimental characterizations.

Additionally, when solving the TDGL equations numerically, it is standard to normalize all the parameters of the model in order to obtain dimensionless quantities. In order to satisfy the above requirements, the steps are as follows: Pick a reference value for each of $\alpha$, $\beta$, $\Gamma$, and $\sigma_{n}$, these will most often just be the corresponding values for the bulk material. Label these reference values $\alpha_{0}$, $\beta_{0}$, $\Gamma_{0}$, and $\sigma_{n0}$. Define dimensionless spatially varying functions, $a(\mathbf{r})$, $b(\mathbf{r})$, $\gamma(\mathbf{r})$, and $s(\mathbf{r})$, relative to their reference values. Apply the following transformations to Eqs. 1 and 2:

	$\displaystyle\alpha$	$\displaystyle\longrightarrow\alpha_{0}a(\mathbf{r})$		(9)
	$\displaystyle\beta$	$\displaystyle\longrightarrow\beta_{0}b(\mathbf{r})$		(10)
	$\displaystyle\Gamma$	$\displaystyle\longrightarrow\Gamma_{0}\gamma(\mathbf{r})$		(11)
	$\displaystyle\sigma_{n}$	$\displaystyle\longrightarrow\sigma_{n0}s(\mathbf{r}).$		(12)

Then proceed with standard nondimensionalization procedures for TDGL (see the Appendix for more details). The advantage of this approach is that the nondimensionalization procedures, when used on these transformed equations, leave behind only the dimensionless functions which capture the spatial variation of the sample material in natural units. The resulting equations are

	$\displaystyle\gamma\left(\frac{\partial\psi}{\partial t}+i\kappa_{0}\phi\psi\right)+\left(\frac{-i}{\kappa_{0}}\nabla-\mathbf{A}\right)^{2}\psi-a\psi+b|\psi|^{2}\psi=0$		(13)
	$\displaystyle\frac{s}{u_{0}}\left(\frac{\partial\mathbf{A}}{\partial t}+\nabla\phi\right)+\nabla\times\nabla\times\mathbf{A}+\frac{i}{2\kappa_{0}}\left(\psi^{*}\nabla\psi-\psi\nabla\psi^{*}\right)+|\psi|^{2}\mathbf{A}=0,$		(14)

where $\kappa_{0}=\frac{\lambda_{0}}{\xi_{0}}$ is the Ginzburg-Landau parameter of the reference material. The quantity $\lambda_{0}=\sqrt{\frac{m_{s}c^{2}\beta_{0}}{4\pi e^{2}_{s}|\alpha_{0}|}}$ is the penetration depth of the reference material, and $\xi_{0}=\sqrt{\frac{\hbar^{2}}{2m_{s}|\alpha_{0}|}}$ is its coherence length. The parameter $u_{0}=\frac{\tau_{\psi_{0}}}{\tau_{j_{0}}}$ is a ratio of characteristic time scales in the reference material, where $\tau_{\psi_{0}}=\frac{\Gamma_{0}}{|\alpha_{0}|}$ is the characteristic relaxation time of $\psi$ in the reference material and $\tau_{j_{0}}=\frac{\sigma_{n0}m_{s}\beta_{0}}{e_{s}^{2}|\alpha_{0}|}$ is the characteristic relaxation time of the current. We have also inserted a minus in front of $a$, which is just a convention to make positive values of $a$ correspond to the superconducting state (Note: this is the opposite of how $\alpha$ is usually interpreted within Ginzburg-Landau theory, however it is standard to make this change when performing nondimensionalization). The boundary conditions become:

	$\displaystyle\left(\frac{i}{\kappa_{0}}\nabla\psi+\mathbf{A}\psi\right)\cdot n=0$		(15)
	$\displaystyle\left(\nabla\times\mathbf{A}\right)\times n=\mathbf{H}_{a}\times n$		(16)
	$\displaystyle\left(\nabla\phi+\frac{\partial\mathbf{A}}{\partial t}\right)\cdot n=0.$		(17)

It should be noted that $\gamma(\mathbf{r})$ and $s(\mathbf{r})$ allow the local characteristic time relaxations to vary, which only affects the dynamics of the solutions. Where the time dynamics are relevant these parameters cannot be ignored; however, in many cases, the primary interest of TDGL calculations is in determining the energetic stability and critical fields, such as the superheating field. In these cases, the time dynamics are not relevant and $\gamma$ and $s$ can be safely set to arbitrary constant values, such as unity.

In the previous section, we have shown how to introduce spatial variation to the TDGL parameters. The process of calculating the values of $a(\mathbf{r})$, $b(\mathbf{r})$, and $\gamma(\mathbf{r})$ is where we bring in experiment and ab initio theory. From Eqs. 6-8, we know these parameters depend on well-defined microscopic quantities, namely $\nu(0)$ and $T_{c}$. These quantities can be calculated using density-functional theory (DFT). Local densities of states are straightforward to compute in DFT, providing local values for $\nu(0)$ and superconducting quantities such as $T_{c}$ are calculated by applying Eliashberg theory within a DFT framework and directly calculating electron-phonon coupling from first principles. Experiment can then give information about the compositions of sample materials, and DFT calculations can determine the $\nu(0)$ and $T_{c}$ associated with these compositions. Using these values, $a(\mathbf{r})$, $b(\mathbf{r})$, and $\gamma(\mathbf{r})$ are calculated from Eqs. 6-8, and the material geometries from the experimental results determine the spatial variation.

In this paper, we will demonstrate our framework on Sn-deficient islands. These defects have been studied extensively [35, 36, 40, 41] so it is straightforward to find estimates of the general size of these islands and their material compositions in the literature. DFT has also been used to calculate $\nu(0)$ and $T_{c}$ for Nb-Sn systems of different concentrations. In this paper we use a commonly observed Sn concentration of 18% and take the DFT values from Ref. 46 to construct $a(\mathbf{r})$, $b(\mathbf{r})$, and $\gamma(\mathbf{r})$. This paper will primarily focus on the results of our TDGL calculations, for further information regarding the details and results from the DFT calculations used in this paper, we refer the interested reader to Refs. 45 and 46.

In this paper we solve the TDGL equations in 3D via a finite element method proposed by Gao [57]. We choose the popular temporal gauge which sets the scalar potential $\phi=0$ (see Du [58] for a more detailed overview of gauge choices for TDGL). Taking the curl of Eq. 14 allows us to directly solve for the magnetic field and we get a system of 3 equations for $\psi$, $\mathbf{A}$, and $\mathbf{H}$:

	$\displaystyle\frac{\partial\psi}{\partial t}+\left(\frac{i}{\kappa_{0}}\nabla+\mathbf{A}\right)^{2}\psi-a\psi+b|\psi|^{2}\psi=0$		(18)
	$\displaystyle\frac{1}{u_{0}}\frac{\partial\mathbf{A}}{\partial t}+\nabla\times\mathbf{H}+\frac{i}{2\kappa_{0}}\left(\psi^{*}\nabla\psi-\psi\nabla\psi^{*}\right)+|\psi|^{2}\mathbf{A}=0$		(19)
	$\displaystyle\frac{1}{u_{0}}\frac{\partial\mathbf{H}}{\partial t}+\nabla\times\nabla\times\mathbf{H}+\frac{i}{\kappa_{0}}\nabla\psi^{*}\times\nabla\psi+|\psi|^{2}\mathbf{H}+\nabla|\psi|^{2}\times\mathbf{A}=0.$		(20)

These equations are subject to the boundary conditions:

	$\displaystyle\left(\frac{i}{\kappa_{0}}\nabla\psi+\mathbf{A}\psi\right)\cdot\mathbf{n}=0$		(21)
	$\displaystyle\mathbf{H}\times\mathbf{n}=\mathbf{H}_{a}\times\mathbf{n},$		(22)

where $\mathbf{n}$ is still the outward normal unit vector and $\mathbf{H}_{a}$ is the applied field. We are then able to solve the system by discretizing the test and trial functions with Lagrange finite elements, and using a backward difference scheme for the time derivatives. All computations in this paper were done using the open source software FEniCS [59].

When simulating SRF materials, dissipation is often a physical quantity of interest. Under TDGL, a formula for dissipation can be derived by considering the time derivative of the free energy and the free energy current flux density. A more detailed derivation is found in Ref. 56, but we quote the final result here:

$$D=2\Gamma\left|\left(\frac{\partial\psi}{\partial t}+\frac{ie_{s}\phi\psi}{\hbar}\right)\right|^{2}+\sigma_{n}\mathbf{E}^{2}.$$

(23)

This quantity is a power density, with the first term corresponding to the superconducting dissipation arising from the relaxation of the order parameter. The second term is the dissipation of normal currents which are largest near the surface where magnetic field can still appreciably penetrate.

It is worth considering how this expression for the dissipation density is related to existing theories of RF power loss and surface resistance. The first term in Eq. 23 is associated with the dissipation due to the rate of change of $\psi$. This term is typically small, except in the vortex state where it becomes the dominant source of dissipation. A dissipation of this form is similar to that proposed by Tinkham [60], who suggested the vortex dissipation should be proportional to $\left(\frac{\partial\psi}{\partial t}\right)^{2}$. The additional term within the parenthesis in Eq. 23 is a result of the gauge invariance of TDGL.

The second term in Eq. 23 can be directly related to the phenomenological “two-fluid model,” which was first proposed by Gorter and Casimir [61] in 1934 and was applied by London [62] later that year to calculate the power loss of a superconductor. The two-fluid model approximates the electrons of the system as consisting of two non-interacting ‘fluids’: the superconducting electrons, in the form of cooper pairs which carry lossless supercurrent, and the normal electrons, which exist as thermally excited quasiparticles that produce typical dissipative currents. Under the two-fluid model, the normal fluid losses produce dissipation of the form [63, 64]

$$P\propto\sigma_{n}E^{2},$$

which is identical to the second term of Eq. 23. For AC applied currents, the electric field is proportional to the frequency and magnetic field, $E\propto\omega H$, meaning that overall the power loss will be of the form

$$P\sim\omega^{2}H^{2}.$$

(24)

It is also well known that within RF cavities, the power loss is given by

$$P=\frac{1}{2}\oint_{A}R_{s}|H(r)|^{2}dA,$$

(25)

where $R_{s}$ is the surface resistance of the cavity walls and $A$ is the surface area. Comparing Eqs. 24 and 25, we have that the second term of Eq. 23 directly leads to a surface resistance proportional to the square of the frequency, $R_{s}\propto\omega^{2}$. In the late 1950s, Mattis and Bardeen [65] as well as Abrikosov, Gor’Kov, and Khalatnikov [66] independently derived the now well-known “BCS resistance.” Under the low frequency and low field limit, the BCS resistance reduces to the form [67, 68]

$$R_{bcs}\simeq\frac{\omega^{2}}{T}e^{\frac{-\Delta}{k_{B}T}}$$

(26)

where $\Delta=1.76k_{B}T_{c}$ is the superconducting energy gap [56]. We thus see that our calculated expression for the surface resistance is in good agreement with the $\omega^{2}$ dependence of the BCS prediction in the low frequency and field limit.

One of the most important quantities that can be directly estimated from the dissipation is the cavity quality factor, Q. Quality factor is given by

$$Q=\frac{2\pi E}{\Delta E},$$

(27)

where $E$ is the energy stored in the cavity and $\Delta E$ is the energy dissipated in the cavity walls each RF period. These quantities (working in SI units for this section) can be expressed as integrals:

	$\displaystyle E=\frac{1}{2}\mu_{0}\int_{V}dV\mathbf{H}^{2}$		(28)
	$\displaystyle\Delta E=\int_{0}^{T}dt\int_{V_{surf}}dV_{surf}D,$		(29)

where $V$ is the cavity volume, $T$ is the RF period, and $D$ is given by Equation 23. $V_{surf}$ is the volume in the first few penetration depths of the cavity surface where essentially all of the dissipation occurs. TDGL simulation outputs are unit-free, so it is helpful to pull constants with units out of these integrals, leaving behind dimensionless functions which can be calculated from TDGL solutions. We start by expressing Equation 28 in cylindrical coordinates:

$$E=\frac{1}{2}\mu_{0}\int rdr\int d\phi\int dz\mathbf{H}^{2}$$

We then define some dimensionless quantities:

	$\displaystyle\tilde{r}$	$\displaystyle=\frac{r}{R}$		(30)
	$\displaystyle\tilde{z}$	$\displaystyle=\frac{z}{L}$		(31)
	$\displaystyle\tilde{\mathbf{H}}$	$\displaystyle=\frac{\mathbf{H}}{H_{a}}$		(32)

Where $R$ is the maximum radius of the cavity, $L$ is the length of the cavity in the axial direction, and $H_{a}$ is the maximum value of the applied field at the surface of the cavity during an RF period. These quantities allow the definition of a unit-less integral that only depends on the cavity geometry:

$$I_{H}\equiv\int\tilde{r}d\tilde{r}\int d\tilde{z}\tilde{\mathbf{H^{2}}}$$

(33)

Using these definitions with Equation 28 and assuming that $\mathbf{H}$ has azimuthal symmetry results in

$$E=\pi\mu_{0}H_{a}^{2}LR^{2}I_{H}.$$

(34)

Turning to the dissipated energy integral, suppose all of the dissipation occurs within a distance $d$ below the cavity surface, where $d<<R$. This allows the cylindrical integral to be converted into cartesian coordinates, with the azimuthal direction becoming the new $x$ direction, the axial direction becoming the new $y$ direction, and the radial direction becoming the new $z$ direction. Doing this, we have

$$\Delta E=\int_{0}^{T}dt\int_{0}^{2\pi R}dx\int_{0}^{L}dy\int_{0}^{d}dzD.$$

(35)

When calculating this from simulation outputs, the integral is necessarily calculated over a small region of the overall cavity surface. Let $L_{x}$ and $L_{y}$ be the simulation domain size in the x and y directions respectively, and let $N$ be the total number of simulation areas needed to fully partition the cavity surface. Then the dissipation integral becomes

$$\Delta E=N\int_{0}^{T}dt\int_{0}^{L_{x}}dx\int_{0}^{L_{y}}dy\int_{0}^{d}dzD,$$

(36)

and $N$ can be approximated as

$$N=\frac{2\pi RL}{L_{x}L_{y}}.$$

(37)

Contuining as before we again define dimensionless coordinates:

	$\displaystyle\tilde{x}$	$\displaystyle=\frac{x}{\lambda}$		(38)
	$\displaystyle\tilde{y}$	$\displaystyle=\frac{y}{\lambda}$		(39)
	$\displaystyle\tilde{z}$	$\displaystyle=\frac{z}{\lambda}$		(40)
	$\displaystyle\tilde{t}$	$\displaystyle=\frac{T_{sim}}{T}t$		(41)

where $\lambda$ is the penetration depth and $T_{sim}$ is the period in units of simulation time. These convert the integral to

$$\Delta E=\lambda^{3}\frac{T}{T_{sim}}N\int_{0}^{T_{sim}}d\tilde{t}\int_{0}^{\frac{L_{x}}{\lambda}}d\tilde{x}\int_{0}^{\frac{L_{y}}{\lambda}}d\tilde{y}\int_{0}^{\frac{d}{\lambda}}d\tilde{z}D.$$

(42)

Additionally, under the temporal gauge ($\phi=0$), Equation 23 can be expressed as

$$D=2\mu_{0}H_{c}^{2}\frac{T_{sim}}{T}\left(\left|\frac{\partial\tilde{\psi}}{\partial\tilde{t}}\right|^{2}+\sigma_{n}\mu_{0}\lambda^{2}\frac{T_{sim}}{T}\left(\frac{\partial\tilde{\mathbf{A}}}{\partial\tilde{t}}\right)^{2}\right),$$

(43)

where $\tilde{\psi}$ and $\tilde{\mathbf{A}}$ are the unit-free versions of the vector potential and order parameter that are solved for with Eqs. 13 and 14 (a derivation of Equation 43 can be found in the Appendix). Finally, we define some more dimensionless integrals over the TDGL solutions:

	$\displaystyle I_{\psi}$	$\displaystyle\equiv\int_{0}^{T_{sim}}d\tilde{t}\int_{0}^{\frac{L_{x}}{\lambda}}d\tilde{x}\int_{0}^{\frac{L_{y}}{\lambda}}d\tilde{y}\int_{0}^{\frac{d}{\lambda}}d\tilde{z}\left|\frac{\partial\tilde{\psi}}{\partial\tau}\right|^{2}$		(44)
	$\displaystyle I_{A}$	$\displaystyle\equiv\int_{0}^{T_{sim}}d\tilde{t}\int_{0}^{\frac{L_{x}}{\lambda}}d\tilde{x}\int_{0}^{\frac{L_{y}}{\lambda}}d\tilde{y}\int_{0}^{\frac{d}{\lambda}}d\tilde{z}\left(\frac{\partial\tilde{\mathbf{A}}}{\partial\tau}\right)^{2}$		(45)

Combining everything and noting that $\omega=\frac{2\pi}{T}$, we get

$$\Delta E=2\mu_{0}H_{c}^{2}\lambda^{3}\frac{2\pi RL}{L_{x}L_{y}}\left(I_{\psi}+\omega\frac{\sigma_{n}\mu_{0}\lambda^{2}T_{sim}}{2\pi}I_{A}\right).$$

(46)

Now using Equations 27, 34, and 46 we get an expression for the quality factor,

$$Q=\frac{\tilde{H_{a}^{2}}RL_{x}L_{y}I_{H}}{2\lambda^{3}\left(I_{\psi}+\omega\frac{\sigma_{n}\mu_{0}\lambda^{2}T_{sim}}{2\pi}I_{A}\right)},$$

(47)

where $\tilde{H_{a}}\equiv\frac{H_{a}}{\sqrt{2}H_{c}}$ is the applied field in simulation units. It is common to express the quality factor as

$$Q=\frac{G}{R_{s}},$$

(48)

where $R_{s}$ is the cavity surface resistance and $G$ is a geometric factor that depends only on quantities which are determined by the cavity geometry. We can define these quantities under the framework we have presented as

	$\displaystyle G$	$\displaystyle=\frac{1}{2}\mu_{0}\omega RI_{H}$		(49)
	$\displaystyle R_{s}$	$\displaystyle=\frac{\mu_{0}\omega\lambda^{3}}{\tilde{H_{a}^{2}}L_{x}L_{y}}\left(I_{\psi}+\omega\frac{\sigma_{n}\mu_{0}\lambda^{2}T_{sim}}{2\pi}I_{A}\right).$		(50)

For a typical 1.3 GHz 9-cell Nb TESLA cavity, $G=270$ $\Omega$ [69], so in practice we can just use this value or other known values of $G$, and only calculate $R_{s}$ from Equation 50.

The $\sigma_{n}$ in the above equations is the conductivity of only the normal electrons. A simple way to estimate this value is to consider the Drude model conductivity [70],

$$\sigma=\frac{ne^{2}\tau}{m},$$

(51)

where $e$ and $m$ are the electron charge and mass, $n$ is the number density of electrons, and $\tau$ is their mean free collision time. In the superconducting state, we can make the usual assumption from the two-fluid model that $n=n_{s}+n_{n}$ where $n_{s}$ and $n_{n}$ are the number densities of the superconducting and normal electrons respectively. Under TDGL, $n_{s}=\frac{\alpha}{\beta}|\tilde{\psi}|^{2}$ and so the normal conductivity is

$$n_{n}=n-\frac{\alpha}{\beta}|\tilde{\psi}|^{2}.$$

(52)

Using this along with Equation 51, the normal electron conductivity in the Meissner state can be estimated with

$$\sigma_{n}=\frac{ne^{2}\tau}{m}\left(1-\frac{1}{n}\frac{\alpha}{\beta}|\tilde{\psi}|^{2}\right).$$

(53)

In the above calculations, it was assumed that the cavity surface is partitioned into small fractional areas, and the dissipated energy is calculated over one of these areas, and then multiplied by the total number of them. If some defect was present in the simulation domain, it would mean that a $Q$ calculation based purely on that value would be assuming that such a defect is uniformly distributed over the surface of the domain. In Equation 50, the only thing that changes when simulating a different material or different defect is the value of the quantity in parenthesis, every other part of the process for calculating $Q$ remains the same. This means that we can calculate a more realistic cavity surface resistance by calculating the $R_{s}$ with Equation 50 for a few different simulations, and then taking a weighted average of these values. Let $R_{i}$ be the surface resistance of the $i$th simulation, and let $p_{i}$ be the percentage of the fractional areas partitioning the cavity surface which are represented by this simulation. Then

$$R_{tot}=\sum_{i}p_{i}R_{i},$$

(54)

where $\sum_{i}p_{i}=1$. The simplest application for Equation 54 is to perform two simulations, one baseline simulation with no defects or material inhomogeneity, and then another simulation containing some defect of interest. Choosing a value $p$ to represent the percentage of fractional areas containing the defect and then following through with the rest of the quality factor calculation provides a simple way to estimate $Q$ for different average defect densities by simply changing the value of $p$.

The value of Q calculated from TDGL outputs will typically be underestimated at low field. This is because of the assumption of gapless superconductivity, which results in higher surface resistances than is predicted with the BCS surface resistance (Equation 26). Despite this, our approach still often predicts quality factors within an order of magnitude of the experimental values. Additionally, the relative behavior of Q at different applied fields, especially when averaging the impact of multiple kinds of defects as described above, qualitatively captures effects such as high field Q-slope.

We start by demonstrating the calculation of quality factor using Eqs. 48 and 50. We simulate a 1.3 GHz Nb cavity ($G=270$ $\Omega$), which typically have quality factors in the range of $10^{10}-10^{11}$. Our simulation domain for these simulations is a cube with side length $20\lambda$ and we use a sinusoidal applied field with a period of $2000\tau_{\psi}$. Since $\tau_{\psi}=\frac{\Gamma}{\alpha}$, we can use Eqs. 6 and 8, and $\nu(0)=90$ states/(eV nm³) (value obtain from private correspondence) to get that for Nb at 2K, $\tau_{\psi}=3.7\times 10^{-13}$ s. As a result, a period of $2000\tau_{\psi}$ corresponds to a frequency of $\sim$1.3 GHz. Our results for this simulation are shown in Figure 2. We see that the low field quality factor is $\sim$$2\times 10^{9}$, this is only a little below the expected range of $10^{10}-10^{11}$, which can be easily explained by the assumption of gapless superconductivity under TDGL, which causes $Q$ to be underestimated. There is a sharp drop in $Q$ once the applied field crosses the superheating field of niobium (the black dotted line) and vortices start to nucleate.

We can also plot the unitless dissipation values that are used to calculate $Q$. Figure 3 depicts the dissipation values for the two terms in parentheses in Equation 46. We find that, as expected, the term due to the dissipation of the electric field ($\frac{\sigma_{n}\mu_{0}\lambda^{2}}{\tau_{\psi}}I_{A}$) is larger than the term due to order parameter relaxation ($I_{\psi}$) by about two orders of magnitude. The $I_{\psi}$ term grows faster than the $I_{A}$ term and the two are equal at an applied field of $\frac{H_{a}}{\sqrt{2}H_{c}}=0.6$, and at higher fields, the $I_{\psi}$ term remains dominant, particularly above the superheating field where it grows much faster due to vortex nucleation and motion. These results demonstrate how the methods in this paper can be used to calculate the dissipation for ideal cases like uniform niobium. We will use them again for the case of Sn-deficient islands in Section III.2.

There has been a considerable amount of work both experimentally and theoretically understanding the various defects present in Nb₃Sn [34, 35, 36, 37, 38, 39, 19, 35, 36, 40, 41]. For the purposes of this paper, we will focus on Sn-deficient islands, small regions within Nb₃Sn grains that have lower than expected Sn concentrations. To demonstrate the use of our sample specific TDGL framework, we will estimate the potential impact to SRF performance by determining the vortex penetration field, $H_{vort}$. This field is the lowest field at which a vortex nucleates, it is a generalization of the superheating field for surfaces containing defects, and it is equal to the superheating field in the limit of a completely uniform flat surface.

In order to simulate Sn-deficient islands, we need estimates of the properties of the material within these regions. We choose a Sn concentration of 18 at.% for our islands. Using this value as well as DFT calculations for Nb₃Sn from Sitaraman et al. [46], we estimate the critical temperature and density of states. Using these values, we construct $a(\mathbf{r})$, $b(\mathbf{r})$, and $\gamma(\mathbf{r})$; the values of these quantities in the bulk as well as in the island are reported in Table 1. The simulation geometry is the same as depicted in Fig. 1, with the domain a cube with side length $10\lambda$. The applied field for these simulations was a constant value.

Figure 4 depicts vortex nucleation occuring at $H_{vort}$ for a particular island. We find that varying the distance of the island from the surface between $0.1$ and $2$ penetration depths (i.e. between $\sim$10-200 nm for Nb₃Sn), we observe a reduction in $H_{vort}$ by as much as $\sim$70% down from the bulk value of the superheating field for islands very near the surface, as shown in Fig. 5. We did this for several different island sizes, and found that the severity of this drop in $H_{vort}$ increases with increasing island size. As such, we can conclude that large islands within 200 nm of the cavity surface are a potential cause for concern when it comes to SRF performance, particularly ones very close to the surface.

Using these solutions, we can also calculate the dissipation due to the Sn-deficient islands, and estimate their quality factors. To do so, we change our applied field to be sinusoidal with a period of 1000$\tau_{\psi_{0}}$, which corresponds to a frequency on the order of $\sim$5 GHz. Most Nb₃Sn cavities use a frequency lower than this, but these simulations already had immense computational cost, so using a $5\times$ longer period would make them prohibitively expensive. When running TDGL simulations, to get good solution convergence, the mesh elements must be smaller than a coherence length (and ideally much smaller). This means simulating Nb₃Sn is extremely challenging, as it has $\kappa=$ $\sim$$26$, meaning the domain both has to be large enough to capture at least a few penetration depths, while also being refined enough to have elements that are a few times smaller than a coherence length. When determining $H_{vort}$, the mesh size (and therefore computational cost) can be reduced by only significantly refining the mesh right near the surface, where vortices will nucleate. Once a vortex begins to nucleate, we know that the applied field is above $H_{vort}$ and thus the rest of that simulation is no longer relevant. However, calculating dissipation requires enough resolution throughout the mesh for a vortex to be able to nucleate, and then move around to dissipate energy, and with a sinusoidal applied field, the vortices can partially nucleate, and then get pushed back out of the superconductor as the field changes. The result of this is that simulations to estimate the dissipation and quality factor of Sn-deficient islands are much more computationally expensive than estimating $H_{vort}$, and so we will only be simulating the smallest island size at a distance from the surface of $0.5\lambda$ and we decreased our domain size to a $5\lambda\times 2.5\lambda\times 2.5\lambda$ rectangular cuboid with the long direction oriented in the same direction as the applied field (and thus the initial orientation of any vortices that nucleate).

Figure 6 shows the results for our quality factors. We ran two simulations, a baseline simulation with no island and a Sn-deficient island with the size described above. Then using Equation 54, we estimated the effective total cavity surface resistance for a number of different island densities. The percentage of island coverage in this figure is calculated by determining the fraction of the total simulation surface area occupied by the island projected into the XY plane and multiplying it by the weighting on the surface resistance from the island simulations used in Equation 54. the 0% island coverage curve (in blue) is the quality factor calculated only from the simulation without an island, and the 6.3% curve (in pink) is the quality factor calculated only from the simulation with an island. We see that above $H_{vort}$ for this island size and distance from the surface vortices begin nucleating, with basically all of the deviations between different island coverages happening only above this value. For cavities with $\geq 1.5\%$ of their surface containing an island, the drop in quality factor starts to become very significant, with even small increases in island density leading to significant drops in quality factor. The noise in the calculations here is largely due to numerical variations in how many vortices enter, and how much they move around. The temperature of the simulation is held constant at the Nb₃Sn operating temperature of 4.2K, so it is likely that these vortices would realistically lead to heating and subsequent cavity quench once the field is much higher than $H_{vort}$. The low field magnitude of the quality factor in these simulations is $\sim$$10^{8}$. This is lower than the expected value of $\sim$$10^{10}$, partly because of being at a higher frequency than 1.3 GHz, and also because the embedded assumption of gapless superconductivity in TDGL means that it will typically underestimate $Q$ (as discussed in Section II.5). Given infinite computational time (or a lower $\kappa$ material), this kind of calculation could be repeated for a variety of island sizes and distances, and applying the defect averaging approach described in Section II.5, the resulting $Q$ plots would look increasingly like the high field Q slope plots taken experimentally from real cavities.