We write for the THz field vector potential \({\mathbf {A}}_T=A_{T0}(t){\mathbf {A}}_{T1}\) with, \({\mathbf {A}}_{T1}({\mathbf {r}},t)= {\mathbf {A}}_{T1}({\mathbf {r}},t+T)\) and \(T=2\pi /\omega\). To follow the evolution of the SC \(A_{T0}\) is ramped up at the rate \(a_T = 16.59 \text {T}\mu\)m/ps until a maximum value \(A_{T0}=59.733 \,\text {T}\mu\)m is reached. In the “Methods” section an analysis is performed of the dominant frequencies contained in the ramped THz field. Of interest are the maxima and minima of \(n_{\mathrm {s}}\) which are displayed in Fig. 2a. The order parameter amplitude oscillates as the field develops in time. For sufficiently strong fields superconductivity can be suppressed. The SC state does not break down homogeneously in the entire sample but respects the local distribution of currents, as set by the vector potential and subject to the boundary conditions (cf. Eq. (4)). This can be seen in Fig. 2b. It should be mentioned that we used relatively strong fields in our calculations. Currently, the reported local maximum amplitude of the B-field are in the range of \(B_\mathrm {T,max}\approx 0.1 T\). The maximum field values in our calculation were on the order of Tesla (cf. Fig. 3).
We note that \(\Psi\) is mostly affected by the eddy currents since their magnitude is several orders higher than that of \(\mathbf{j}_{\mathrm {s}}\). For very high field strengths the normal currents start to destroy superconductivity at the center of the SC. Also the currents flowing around the edges become relevant and lead to the formation of several superconducting islands in the sample. The breakdown of superconductivity at the center of the sample is not only a consequence of the strong current flow. Figure 4a evidences that the vector potential \({\mathbf {A}}_T\) changes its spatial orientation rapidly at the disc center. e.g., the x component changes from positive to negative on a length scale that is comparable to the GL coherence length. The superconducting phase field tends to follow the external potential. The steep phase gradient created in this way also affects the order parameter density and leads to a gradual suppression of the superconductivity at the disc center. In turn, the reduction of \(\Psi\) affects the supercurrents since their magnitude is directly proportional to \(n_{\mathrm {s}}\). In contrast, the normal current distribution is only slightly affected by the breakdown of superconductivity. The distribution of the total current is shown in Fig. 4b. The overall features of the current align with its relation Eq. (4) to the vector potential (Fig. 4a) and the imposed boundary conditions forcing the current density to flow around the cylinder which leads to the formation of local whirls in the current density.
One may expect the current flow to be governed by the spatial variation of the applied vector potentials. The time evolution of the phase field \(\theta\), and the order parameter amplitude \(n_\mathrm {s}\) show, however, a non-trivial time evolution. Technically, this is because both quantities are calculated from coupled and nonlinear partial differential equations. On the other hand, the dynamics of \(\Psi\) happen on a finite time scale and the order parameter can not adapt instantaneously to the external vector potential. These effects alter the current flow in the disc leading to a quite involved time evolution.
To study the SC behavior during switching off the field, we monitor the SC dynamic starting from the initial value \(A_{T0}=59.733\,\)Tμm and decreasing it to zero (after \(t=0.18\,\)ps) at a rate of \(a_{T} = 331.8\,\)Tμm/ps. The order parameter is strongly suppressed when the field is switched off.
The field was ramped down with a rate that is much higher than the initial ramping rate. Therefore, the normal currents generated by the time variation of \({\mathbf {A}}\) are higher in this case and reduce the superconducting order parameter. In this process the phase \(\theta _{\mathrm {s}}\) develops a steep gradient which leads to delicate supercurrent distributions (see Fig. 5a,b), which is mostly in-plane. The time evolution of the phase suggests that also vortex-antivortex (VAV) pairs are created in the rapid decrease of \(|\Psi |^2\) (not shown here). They appear as little knots in the phase field close to the border of the sample. But since no additional magnetic field is applied the VAV-pairs are highly unstable and annihilate eventually. For the stabilization, a static magnetic field is needed which is by itself not capable of creating the VAV pairs but it preserves them, once created by \({\mathbf {A}}_T\). The normal current also displays a more involved time evolution than in the previous case. A number of magnetic whirls along the sample edges are generated. However, these whirls are not equivalent to the conventional superconducting vortices since there is no phase singularity and no normal conducting vortex core.
To illustrate the role of the orbital and vectorial texture of \({\mathbf {A}}_T\), we investigate how the superconductor reacts to a plane-wave THz field (no spatial textures). This amounts to using the vector potential
$$\begin{aligned} {\mathbf {A}}_T(x,t) = -a_Tt\sin \left( \frac{2\pi }{\lambda }\left( ct-x \right) \right) \mathbf{e}_y+\nabla u \end{aligned}$$
(9)
as an input for the first time dependent GL-equation. The function \(u({\mathbf {r}}, t)\) is chosen such that \(\nabla \cdot \mathbf{A}_\mathrm {T} = 0\) and \({\mathbf {A}}_\mathrm {T}\cdot {\mathbf {n}}=0\) holds at all times. The corresponding magnetic field has only components in the z direction. To reduce the numerical efforts, the calculations are performed for a circular SC of radius \(R = 250\,\)nm located in the plane \(z = 100\,\)nm. Furthermore, we neglect the SC self field \(\mathbf {\mathbf {A}}_\mathrm {s}\) which is allowed in the cases where \(\lambda _\mathrm {eff}\gg \xi _\mathrm {GL}\)41. As before, we choose the frequency of the potential to be \(f=10\,\)THz with a corresponding wavelength \(\lambda = 30 \,\) \(\mu\)m, i.e \(\lambda \gg R\). Thus, the space variation of \(B_\mathrm {T} = \nabla \times {\mathbf {A}}_\mathrm {T}\) is negligible and the external field is basically homogeneous in space. The prefactor in Eq. (9) has the value \(a_\mathrm {T} = 0.41\,\)Tμm/ps, meaning that after \(t=5\,\)ps the external magnetic field has reached a value of \(B_\mathrm {T}=0.39\,\)T. The time evolution of the system is shown in Fig. 6. As in the previous calculations, the effect of \(A_T\) on the superconductor consists of gradual suppression of the order parameter. However, for the plane-wave vector potential, the suppression of \(\Psi\) happens solely around the sample edges, whereas the center remains superconducting, highlighting so the role of finite-size effects. With an increasing external field strength, the normal conducting ring domain grows until the material becomes completely normal conducting again. An interesting observation is that the phase of the order parameter is only slightly affected by the external field. In contrast, the field created by the spintronic THz emitters results in a strong modulation of the order parameter phase. Also, the normal conducting domains at intermediate field strengths have a more complicated structure.
In light of the above observations, the question arises of how an ensemble of superconducting vortices reacts to structured THz field pulses. To address this issue, an additional external magnetic field \({\mathbf {B}}_\mathrm {e} = B_\mathrm {e} {\mathbf {e}}_\mathrm {z}\) is applied to the sample with \(B_\mathrm {e} = 0.17\,\)T. For such a field, the equilibrium state of the system consists of a circular arrangement of 8 vortices (cf. Fig. 7a)a. Due to finite-size effects, the vortices experience a strong confining potential which is brought up by the Meissner currents flowing around the sample edges. As a consequence, the vortex lattice reflects the symmetry of the sample. In the next step, we investigate how the state under finite \(\mathbf{B}_\mathrm {e}\) is affected by THz field pulses of different topology, amplitude, and frequency. We start by comparing the vortex behavior under a conventional plane-wave THz field pulse with the vortex dynamics steered by the STEs described above. In Fig. 7a, one can see that the application of the plane-wave field leads to a stronger order parameter suppression compared to the STE-pulse, despite a weaker field-magnitude. This is because the fields lead to different current distributions in the sample, which in turn affect the order parameter differently. Another observation is made when comparing the right panels of Fig. 7a,b): The plane-wave pulse is able to drive the vortices out of the sample for sufficiently high field amplitude. In contrast, the STE-pulse preserves the topology of the phase field and the vortex lattice experiences a structural change without vortex creation or annihilation. It should be noted that the qualitative dynamic behavior of the vortex ensemble is found to be independent of the rate at which the external field is increased.
The frequency of \(f=10\,\)THz is substantially higher than the gap frequency of Nb \(f_{Nb}=2\Delta _{\mathrm {Nb}}/h=701\,\)GHz and the radiation field has, in addition, a pair breaking effect on the superconducting condensate. In this respect, cuprate superconductors with gaps on the order of \(\Delta \sim 20\,\)meV would be better-suited42. To avoid effects related to pair breaking, we investigated how the dynamics of the system change if the pulse-frequency is reduced down to \(f=500\,\)GHz (cf. Fig. 7c). The general behavior of \(n_\mathrm {s}\) is the same as in the case of \(f=10\,\)THz. However, the maximum value of \(n_s\) exhibits stronger oscillations since the order parameter is offered more time to adapt to the more slowly varying external field. The vortex lattice is also stronger affected by the pulse with \(f=500\,\)GHz. In particular, the vortex cores are notably deformed even at moderate field strength values. In addition, we find that the low-frequency pulse has a certain tendency to make the vortices clump together (see panel two in Fig. 7c). This tendency becomes even more clear in the high field regime, where we have essentially four giant vortices instead of 8 distinct vortex cores. Furthermore, the phase field \(\theta\) indicates the temporal appearance of vortex-antivortex pairs. As in previous calculations, the observed vortex dynamics are essentially independent of the rate at which the external field increases, as long as the frequency is kept constant.
To highlight the role of geometric confinement in our system we performed the calculations discussed above on samples of different geometry and spatial extension. The applied STE-pulse has a frequency \(f=10\,\)THz and a ramping factor \(a_T = 1.659\)Tμm/ps. The sample is a square with side length \(a = 250\,\)nm or a superconducting disc of radius \(r = 1\, \upmu\)m (cf. Fig. 8).
A consequence of the modified geometry is inferred from the upper panels on the right, where the vortices form a simple cubic lattice in the square sample. The current flow around the edges of the square creates an energetic barrier that prevents an immediate relaxation of the vortex lattice into a hexagonal arrangement. However, if this vortex ensemble is subject to the THz pulse, one can see a very similar behavior as in the case of the THz field-driven disc with \(r=250\,\)nm. With increasing the amplitude of the magnetic field, the vortex cores grow in size until most of the sample becomes normal conducting. But also, in the square sample, the THz field pulse is not changing the topology of the order parameter phase-field \(\theta\) and the vortex number is preserved (not shown here). The situation is different if we return to the superconducting disc with an increased radius of \(r=1\, \upmu\)m (cf. lower panels on the right in Fig. 8). This system is less affected by geometric confinement and the vortex lattice has adopted a hexagonal symmetry. An important difference from the previous systems is the higher sensitivity of the order parameter to the applied THz field. Already at a maximum value of \(B_\mathrm {T}=0.16\,\)T, one can observe the formation of four almost normal conducting domains. At \(B_\mathrm {T}=0.23\,\)T, the system is mostly normal conducting. Due to the presence of the fully developed vortex lattice, the larger system size and the lesser importance of the boundary currents, the total current distribution in the disc differs from that of the previous systems. That also means the currents induced by the THz field pulse have a different impact on the order parameter development. For even larger systems, the finite spatial extent of the radiation field will become an important factor as well, and superconductivity is only suppressed at the center of the disc. The normal conducting domains induced by the time-varying magnetic field and the local heating of the sample then might provide another way to detect and characterize the properties of the pulse. Particularly, the formation of stable vortex clusters of characteristic shape inside the normal domains might serve as a fingerprint for the specifics of the radiation field43.
To explore the effect of the THz fields on transport, we inject an electric transport current into the superconductor in the presence of the STE generated THz field pulse. The idea is to simulate a system which is similar to a superconducting single-photon detector. It is known that in such systems the absorption of an incoming photon creates a local hot spot, which eventually leads to the breakdown of superconductivity in the entire sample44. The transition into the resistive state is usually accompanied by the appearance of highly mobile vortex-antivortex pairs. Here, we investigate if such a transient dynamic state can also be created with the aid of a structured THz field pulse. The external current is introduced into our system by solving the continuity equation
$$\begin{aligned} 0 = \nabla \cdot {\mathbf {j}}_\mathrm {s} – \sigma \Delta \phi \end{aligned}$$
(10)
with \(\sigma = 38\,(\mu \Omega \mathrm {m})^{-1}\) and boundary conditions
$$\begin{aligned}&\Psi = 0, ~~ \text {on} ~~ \partial \Omega _e \end{aligned}$$
(11)
$$\begin{aligned}&- \sigma \nabla \phi \cdot {\mathbf {N}}_e = j_e, ~~ \text {on} ~~ \partial \Omega _e \end{aligned}$$
(12)
where \(\Omega _\mathrm {v}\) and \(\Omega _\mathrm {e}\) are the SC/vacuum and SC/electrode boundaries, and \({\mathbf {N}}_\mathrm {v}\) and \(\mathbf{N}_\mathrm {e}\) are the corresponding surface normals45. In our case the external current flows along the y-direction and the superconductor/electrode boundaries are located at \(y=\pm \,250\,\)nm. In a first step the external bias is quasi-adiabatically increased until the critical current \(j_\mathrm {c}\) is found. In our case the system is observed to transit into the resistive state for \(j_\mathrm {e}=j_\mathrm {c} = 0.45\,\)GA\(/\mathrm {cm}^2\). Next, the external current is fixed to \(j_\mathrm {e}=0.95j_\mathrm {c}\) and the system is allowed to relax into the equilibrium fully. Finally, the current-driven superconductor is subjected to the STE emitted field pulse with a frequency \(f=10\,\)THz and a ramping factor \(a_T = 0.166\)Tμm/ps. In the simulations, the value of the transport current is chosen relatively high, which we expect to render the superconductor more susceptible to the applied THz radiation. In a real experimental setup, the applied radiation field may heat the sample locally, providing an additional driving factor for the breakdown of superconductivity, an effect which is not taken into account here.
Fig. 9 shows the general trend of the superconducting dynamics in a square sample with a side length \(a=250\,\) nm under the application of the STE generated field pulse with \(f=10\,\)THz and \(a_T = 0.166\)Tμm/ps. The field is ramped up slower than in the previous calculations and the order parameter is only gradually affected at first. Once the field amplitude is high enough, the system transits quickly into the normal conducting state. The transition by itself was not observed to be accompanied by the emergence of vortex-antivortex pairs, as in the case of the classic single-photon detector. Instead, the breakdown of superconductivity starts at the superconductor/electrode boundaries, where the order parameter is already suppressed. The magnetic field of the THz field pulse drives the growth of the normal conducting domains until the entire material becomes normal conducting again. This behavior can also be observed for THz field pulses of lower frequency and ramping rate (not shown here). An important factor is the relatively small size of the system. In larger superconductors, the THz field might affect the order parameter evolution differently and vortex formation might become possible. This might be especially interesting if we take the local heating of the sample into account. In this way, the combination of a transverse current flow and a quench of the locally induced hot spot could lead to the formation of vortex avalanches which are characteristics for the applied THz field pulse46.