Imaging anisotropic vortex dynamics in FeSe
Irene P. Zhang Johanna C. Palmstrom Hilary Noad Logan Bishop-Van Horn Yusuke Iguchi Zheng Cui John R. Kirtley Ian R. Fisher Kathryn A. Moler











Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA


Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA


Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA

















Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
Department of Physics, Stanford University, Stanford, CA 94305, USA


Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA


Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA


Geballe Laboratory for Advanced Materials, Stanford University, Stanford, CA 94305, USA

Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
Geballe Laboratory for Advanced Materials, Stanford University, Stanford, CA 94305, USA

Department of Applied Physics, Stanford University, Stanford, CA 94305, USA
Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
Geballe Laboratory for Advanced Materials, Stanford University, Stanford, CA 94305, USA



Strong vortex pinning in FeSe could be useful for technological applications and could provide clues about the  coexistence of superconductivity and nematicity. To characterize the pinning of individual, isolated vortices, we simultaneously apply a local magnetic field and image the vortex motion with scanning SQUID susceptibility. We find that the pinning is highly anisotropic: the vortices move easily along directions that are parallel to the orientations of twin domain walls and pin strongly in a perpendicular direction. These results are consistent with a scenario in which the anisotropy arises from vortex pinning on domain walls and quantify the dynamics of individual vortex pinning in FeSe.



                              
Imaging anisotropic vortex dynamics in FeSe
    Kathryn A. Moler
    December 30, 2023
===========================================










§ INTRODUCTION

FeSe is a particularly simple layered iron based superconductor (Fe SC), with an un-strained critical temperature of about 8 K <cit.>, an increase in this critical temperature with pressure to 36.7 K at  8 GPa <cit.>, and with a report of interface-induced high-temperature superconductivity above 50 K in single unit-cell films on SrTiO_3 <cit.>. In general, iron-based superconductors exhibit complex interplay between superconducting, nematic, and magnetic orders, and undergo tetragonal-to-orthorhombic phase transitions close to magnetic ordering transitions <cit.>. Unlike other Fe SCs, FeSe does not magnetically order, providing an opportunity to study the superconductivity-nematicity relationship without the added complexity of an ordered magnetic state. The structural transition of FeSe occurs at 90 K and is understood to be driven by electronic nematic order <cit.>. In the orthorhombic state, FeSe and other Fe SCs form domains separated by twin boundaries (TBs). The superconducting pairing mechanism has been discussed in terms of spin-fluctuation pairing <cit.> but the nodal character of the gap remains controversial. Some superfluid density <cit.>, thermal conductivity <cit.>, and tunneling spectroscopy measurements <cit.> are consistent with line nodes in the orbital component of the superconducting order parameter; however, recent thermal conductivity <cit.>, STM <cit.>, and London penetration depth <cit.> studies suggest that FeSe is fully gapped but with deep gap minima. The relationship between nematic order and superconductivity is likewise controversial, with NMR studies suggesting that nematic order competes with superconductivity <cit.> and heat capacity and thermal expansion studies suggest that it enhances superconductivity <cit.>. 
 
Although vortex pinning in superconductors is of great practical as well as fundamental interest, its mechanism is still poorly understood  <cit.>. Recently it has been possible to study this pinning directly by imaging individual vortices while manipulating them.  Auslaender et al. <cit.>  dragged vortices in a cuprate superconductor over distances of a few microns using a magnetic force microscope tip. They found an enhanced response of the vortex to pulling when the tip was oscillated transversely. They also found enhanced vortex pinning anisotropy, which they attributed to clustering of oxygen vacancies in their sample. Later work from Shapira et al. used MFM to drag vortices along twin boundaries in  YBa2Cu3Ot-δ and demonstrated that the vortices moved in a series jumps, consistent with power-law behavior <cit.>. Kalisky et al. <cit.> showed that vortices, when dragged by a scanning SQUID microscope, avoided crossing twin boundaries in underdoped Ba(Fe_1-xCo_x)_2As_2. Embon et al. <cit.> used a SQUID on a tip to image the movement over a few tens of nanometers of vortices driven by  applied supercurrents in a thin Pb film. They were able to map out anisotropic and spatially inhomogeneous pinning forces on the vortices, which they attributed to multiple overlapping pinning sites. 
 
The vortex pinning properties of FeSe have attracted interest <cit.> and make it a potential competitor to high-T_c cuprates for high field applications <cit.>. Critical current density studies of vortex pinning in FeSe have found that it is dominated by strong point-like pinning <cit.>, while STM studies have shown that vortices preferentially pin on TBs in FeSe, where the superfluid density is reduced <cit.>. In this paper we present scanning SQUID magnetometry and susceptibility images of vortices trapped in single crystals of FeSe. The susceptibility images show structures that we attribute to motion of the pinned vortices driven by the magnetic fields applied by the field coil integrated into our SQUID susceptometers. Analysis of our data using a simple model is consistent with a quadratic dependence of the restoring force on displacement along the direction of the TBs, with a much larger restoring force in the orthogonal direction. The anisotropy in the in-plane restoring forces can be as large as a factor of 20. 




§ MATERIALS AND METHODS




We used scanning SQUID microscopy to image magnetometry and susceptibility 
in bulk single crystal FeSe. Single crystal FeSe samples were grown by 
chemical vapor transport following the growth procedure outlined in Ref. boehmer2016growth and exfoliated with Kapton tape and silver paint (Dupont 4929N) as in Ref. massat2017 to achieve a surface flat enough 
to be scanned using our susceptometers. The superconducting transition and structural transition temperatures of the batch of samples used were found to be 8.2 K and 88 K respectively. The bulk superconducting transition temperature was extracted from magnetic-susceptibility measurements taken using the vibrating sample mount option of the MPMS 3 from Quantum Design. The structural transition temperature was determined from resistivity measurements on a free standing crystal taken using the MPMS 3 paired with a Linear Research Model LR-700 AC resistance bridge. The susceptibility and resistivity data are shown in Fig. <ref>.

Our SQUID susceptometers contain two Nb pickup loop/field coil pairs 
arranged in a gradiometric layout <cit.>. The pickup loop and 
field coil are covered by Nb shielding so that flux passes only through the 
loop and not through the gaps between the leads. The inner radius of the 
pickup loop was 0.3 um, resulting in sub-micron spatial resolution. As the 
susceptometer scans across the surface of the sample, we record the magnetic
flux passing through the pickup loop. The dc signal is recorded as 
magnetometry and is reported in units of the flux quantum Φ_0 = h/2e. 
We use an SR830 lock-in amplifier to pass an ac current through the field 
coil, creating a local magnetic field, and record the ac flux at that 
frequency to measure susceptibility. The gradiometric design cancels out the
flux due to the field coil so that the pickup loop only measures the 
magnetic response of the sample. The susceptibility is normalized by the 
lock-in amplifier current and is reported in units of Φ_0/A. This 
design allows us to image local magnetic fields and susceptibility at the 
surface of a sample simultaneously.



§ RESULTS



 §.§ Imaging vortex motion


In four measured samples, we found two-lobed features ("butterflies") in the susceptibility 
images, accompanying superconducting vortices. Not all vortices had these 
“butterflies" and their brightness varied from vortex to 
vortex. Increasing the sample temperature increases the brightness of the 
lobes up to Tc, after which both the vortex and accompanying butterflies 
disappear. The butterflies orient along one of two directions which are 
perpendicular to each other. Further, the relative brightness of the two 
lobes is consistent within butterflies of the same orientation, and the 
dimmer lobes are located in areas opposite to the location of the field coil
shielding. The brightness of the butterflies was varied from vortex to 
vortex, ranging from around one tenth of a Φ_0/A to a few Φ_0/A. 
Besides wiggling the vortex, the field coil could also push the vortex to a different pinning site. Not all vortices had butterflies, likely because they
were pinned too strongly to show a signal in susceptibility. Due to
unevennesss in the sample surface, it was difficult to control for the
height and angle of the pickup loop/field coil pair relative to the surface
area being scanned.
We used x-ray diffraction to determine the tetragonal crystal 
axes in Samples 1 and 2 and found that a line cutting through the two lobes orients the butterflies either along or perpedicular to the tetragonal [100] direction. This orientation is 
consistent with the butterflies being aligned with TBs. Fig. <ref> shows two representative susceptibility butterflies from Sample
1, along with the corresponding magnetometry image showing the vortex for 
the first butterfly. A background is subtracted from all susceptibility 
images by fitting a plane to the four corners of the scan. 



In Fig. <ref>, we show a series of susceptibility scans for three vortices with butterflies at variouss temperatures, as well as a temperature series showing the development of stripes in diamagnetic susceptibility close to T_c. Near the superconducting transition, the vortex becomes easier to move and the features of the butterflies become sharper. In one sample we also observed striped variations in susceptibility run along the tetragonal [100] direction.



 §.§ Modeling vortex motion



We model our susceptibility images by 1) calculating the magnetic fields inside the superconductor at the vortex position due to the applied currents through the field coil, 2) calculating the motion of the vortex in response to these fields using a simple model with an anisotropic pinning potential, and then 3) calculating the change in flux through the susceptometer due to the vortex motion. The SQUID susceptibility is given by the response flux Φ divided by the field coil current I.


  §.§.§ Applied fields


Consider a geometry in which a scanning SQUID susceptometer, composed of niobium films with penetration depth λ_ Nb, with layout in the pickup loop/field coil region as illustrated if Fig. <ref>b, is assumed oriented parallel to the sample surface, has a spacing z_0 between the surface of the susceptometer and the surface of the sample, and is in the half-space z>0 (region 1). The superconducting sample, with penetration depth λ, is in the half-space z<0 (region 2). The fields generated by the susceptometer are calculated following Ref. brandt2005thin as described for our sensors in Ref. KirtleyRevSciInstrum16. For a thin superconducting film with thickness t < 
λ_Nb, one can define a stream function g by J = ẑ×∇ g, where J is the sheet current <cit.>. One interpretation of the stream function is that it defines a density of magnetization in the ẑ direction, or equivalently a collection of small current loops in the plane. Once the stream functions g_j,l are known for all of the grid points j in all of the superconducting layers l in the susceptometer, the source potential at any point i in the half-space z>0 but outside of the susceptometer superconducting layers is given by <cit.>

    φ_s(r⃗_i,z_i)=-∑_l∑_j w g_j,l/4πz_i-z_j,l/((z_i-z_j,l)^2+ρ_i,j,l^2)^3/2,

where w is the area of the pixels used, and ρ_i,j,l=√((x_i-x_j,l)^2+(y_i-y_j,l)^2). For what follows we calculate the source potential φ_s for all r⃗_i and z=0.

Once the source fields are known, we match boundary conditions at z=0. This is done by expanding the scalar magnetic potential φ_1(r⃗,z) outside the sample (z>0) and the magnetic field H⃗_⃗2⃗(r⃗,z) inside the sample (z<0) in Fourier series that are constructed to satisfy Maxwell's equations for z>0 and London's equation for z<0 <cit.>.


    φ_1(r⃗,z)= 1/(2π)^2∫ d^2κ(φ_s(κ⃗)e^κ z+φ_r(κ⃗)e^-κ z)e^i κ⃗·r⃗



    H⃗_2(r⃗,z) =1/(2π)^2∫ d^2 κ h⃗_2(κ⃗) e^q z e^i κ⃗·r⃗,

where κ⃗=κ_x x̂+κ_y ŷ, κ = √(κ_x^2+κ_y^2), q^2=κ^2+1/λ^2 and H⃗_1 = ∇⃗φ_1. Using the boundary conditions B⃗·ẑ and H⃗×ẑ continuous at z=0, as well as ∇⃗·B⃗_⃗2⃗=0 and B⃗=μ_0 H⃗ results in <cit.>


    h_2,z(κ⃗) = 2κ^2/κ+qφ_s(κ⃗)



    h⃗_2,||(κ⃗)=2 i κ⃗ q/κ+qφ_s(κ⃗)

The magnetic field H⃗_2(r⃗,z) penetrating the superconductor is calculated taking the Fourier transform of h⃗_2(κ⃗) (Eq. <ref>). 



  §.§.§ Vortex motion

The additional magnetic energy δ u of the vortex per unit length due to the susceptometer applied field is given by

    δ u(r⃗) = 1/2∫ d^2r'  H⃗_2(r⃗') ·B⃗_V(r⃗-r⃗')

where B⃗_V(r⃗'-r⃗) is the field of the vortex centered at x,y. If we assume the vortex is oriented in the ẑ direction and neglect the finite extent of the vortex fields, this becomes simply

    δ u(r⃗) = Φ_0 H_2,z(r⃗)/2

The total additional energy δ U of the vortex due to the field penetration into the superconductor is then

    δ U(x,y) = Φ_0/2(2π)^2∫_0^-∞ dz  ∫ d^2κ  h_2,z(κ⃗) e^q z e^i κ⃗·r⃗

The integration over z can be done analytically, resulting in

    δ U(x,y) = Φ_0/2(2π)^2∫ d^2κ  h_2,z(κ⃗) e^i κ⃗·r⃗/q

The force on the vortex due to the susceptometer fields is the gradient of this extra potential:

    F⃗_ SQUID(x,y) = ∇⃗(δ U(x,y))


We find that the force is reduced when the field coil shield is above all
or part of the vortex, meaning that when comparing scans to the 
corresponding susceptometer layout, the areas opposite to the field coil 
shield will show reduced signal.

Because we do not know the exact form or mechanism of the vortex pinning potential in FeSe, we 
use a simple quadratic model with spring constants k_w and k_s  (k_s ≥ k_w) associated with orthogonal axes ŵ and ŝ rotated by an angle θ in the ab plane relative to the scan axes x̂ and ŷ . When the susceptometer scans in the xy plane relative to the vortex position, the susceptometer applied fields pull the vortex towards (or away) from the pickup loop/field coil. Since we modulate the current through the field coil at about 1 kHz, we assume that the vortex response is much faster than the applied force, and therefore that the displacement of the vortex d⃗r⃗ = dw ŵ+ds ŝ from its equilibrium position can be calculated from the balances of forces condition

    F⃗_ SQUID = k_s ds  ŝ + k_w dw  ŵ.



  §.§.§ Response flux due to vortex motion


The expected ac flux from the vortex motion is calculated using the 
gradient of the dc flux through the susceptometer pickup loop due to the vortex:

    Φ_ac = dΦ_dc/dx dx + dΦ_dc/dy dy

The flux due to a vortex is calculated following the procedure in Ref. KirtleyRevSciInstrum16; 
the derivatives in Eq. <ref> are taken numerically, and the susceptibility is calculated by dividing Φ_ac by the applied field coil current.





  §.§.§ Results of the model



Example predictions for the ac flux due to vortex motion with the susceptometer geometry 
in Fig. <ref> is shown in Fig. <ref> for various assumed values of k along the x- and y-axes of the images, with z_0=2 μm. The resulting 
shapes for the isotropic cases (k_s = k_w) are similar to an incomplete 
torus, while the shapes for the anisotropic cases (k_s ≠ k_w) are 
distinctly lobed. The weaker axis (k_w) spring constant determines the intensity of the susceptibility signal,
while the stronger axis spring constant (k_s) largely determines the shape. As the 
strength of the strong axis spring constant is increased, the lobes become more axially symmetric, and the dark region in the center fades. Consistent with the force
profile, the lobes located across from the field coil shielding are dimmer 
than their partners. Crucially, we find that the apparent axis of the 
anisotropy is not qualitatively changed by the asymmetric susceptometer geometry.



Using reduced chi-square fitting, we calculated the optimal spring constant 
in the weak axis (“k_w”) direction and a lower bound for the spring 
constant in the strong axis direction (“k_s”) for 12 butterflies from two samples. We find 
that in all cases the signal along the hard axis is so weak compared to the noise level 
that the optimal k_s value approaches infinity. In Fig. <ref> we show the fitting process for the butterflies shown in Fig. <ref>. We fit k_w and k_s separately by taking cuts along the weak and strong axes, respectively. Full images for the data, model, and difference between data and model are shown in Fig. <ref> The difference in brightness between model and experiment in the dimmer lobe can be explained by the limited spatial extent of the field coil shielding in the model SQUID geometry.



Taking the optimal k_w
value and the lower limit for k_s results in a lower limit for the 
ratio k_s/k_w, which we use to characterize the anisotropy of vortex motion. The range of k_s/k_w varies from around 4 to over 20, and three vortices show ratios above 20. The results are summarized in Table <ref>, along with the temperature at which the data were taken.



Another way to infer the vortex dynamics in these samples is to directly 
extract the vortex displacement from scans. This approach avoids the 
assumptions of a toy model potential entirely. Providing that the direction 
of the vortex motion is known, we can use Equation <ref> and the gradient from the magnetometry image to obtain the maximum vortex displacement. In this analysis we 
exclude the values close to the vortex center and far from the vortex, as
the magnetometry gradient is very small in these areas and amplifies any 
noise in the susceptibility. This is the case along the lines where F_s = 
0 and F_w = 0, respectively. Furthermore, since the motion in the ŝ
direction is severely limited, we can assume the vortex moves mostly in the 
ŵ direction. The motion along lines through F_s = 0 and with constant but 
small F_w is therefore mostly radial. In Fig. <ref> we show line cuts of the 
experimental vortex displacement and calculated susceptometer force along the weak and strong axes of 
a butterfly from Sample 1. Because the signal along the strong axis is 
dominated by the background, we subtracted a linear fit from susceptibility 
line cuts along this axis. The absolute magnitude of the vortex displacement is at least an order of 
magnitude smaller in the strong-direction compared to the weak-direction. 
Furthermore, while the displacement along the weak axis tracks with 
susceptometer force, the displacement along the strong axis appears to be 
largely independent of the force applied by the field coil. A linear fit to 
the weak axis data is plotted in Fig 1(b). The effective k_w in this case 
is 1.0 ± 0.06 × 10^-8 N/m, about 2.5 times the optimal k_w found by
fitting the data to the quadratic model. These results 
are consistent with a much larger spring constant along the strong axis than along
the weak axis, and also consistent with a linear 
force-displacement relation in the weak direction (Fig. <ref>(b)).




§ DISCUSSION


The ac current through the field coil generates an oscillating local magnetic field and superconducting screening currents
that pull and push the vortex towards 
and away from the center of the pickup loop/field coil.  


As a result, as long as the vortex response to the applied force is faster than the time variation of the field coil current, the flux through the pickup loop is in-phase with 
the field coil current, as
reflected in the real part of the susceptibility image. This scenario is rather general to 
scanning SQUID measurements with simultaneous magnetometry and 
susceptibility imaging and does not rely on the particulars of the FeSe 
samples, except that in FeSe the motion is both large and highly 
anisotropic.

We can first rule out that the motion is isotropic but that the observed lobed shape is 
simply due to the SQUID susceptometer geometry based on the simulations in Fig 
2. The expected ac signal in the isotropic case is proportional to the 
susceptometer force, which is mostly radially symmetric except for shielded 
regions. The force is small when the pickup loop/field coil center is above 
the vortex, increases to a maximum at a distance equal to the field coil 
radius and then decreases again as the field coil moves away from the 
vortex. This force profile results in the “doughnut" shape discussed 
earlier. The effect of increasing the spring potential along one axis is to 
weaken the effect of the force along that axis so that as k becomes very 
large the signal along that axis becomes indiscernible, producing lobed 
features. 

As noted in the introduction, previous MFM and scanning SQUID microscopy measurements have measured vortex
dynamics in superconductors<cit.>. 








In contrast to the Auslaender et al. and Kalisky et al. studies, in
which the vortices were dragged to new locations on the sample, in our case and for small displacements in the Embon et al. <cit.> study the vortex returns to its original location and is instead oscillated in place. In the Embon et al. study the displacements are typically a few tens of nanometers, the forces are pico-Newtons, and the spring constants are ∼ 10^-5 N/m, whereas in the present study the displacements are typically one micron, the forces are typically femto-newtons, and the weak spring constants are typically 10^-8 N/m, 3 orders of magnitude weaker in FeSe than in Pb.

A 2012 STM study <cit.> found that in thin film FeSe, superconductivity is suppressed along twin domain boundaries and that vortices tend to pin on these boundaries. The cause of this suppression is thought to be due to the increased height of the Se atoms along the boundary. We did not observe chains of vortices along any direction; however, it is possible our vortex density was not high enough to observe such patterns. In three samples for which we were able to determine the direction of the crystalline axes, the vortex motion was highly anisotropic along the TB directions, suggesting that it is easier to pull vortices along the TBs than acorss them. Further evidence of this can be found in the striped variations in susceptibility signal close to T_c which run along the same direction as the weak axis of the butterflies. Given that the only other known symmetry-breaking through lattice orthorhombicity, for which the crystal axes are oriented at 45 degrees to the butterfly and stripe direction, the only reasonable conclusion is that the vortices are pinned on TBs and the variations susceptibility reflect the suppressed superconductivity on the boundaries.



§ CONCLUSION


In conclusion, scanning SQUID susceptibility can be used to image the motion of vortices under the influence of locally applied magnetic fields. We have applied this method in the instance of vortices in FeSe. Detailed calculations of the magnetic fields generated by the susceptometer, combined with a simple model for the pinning forces on the vortex, show that these pinning forces can be highly anisotropic and is consistent with a quadratic dependence of the restoring force on displacement. We calculate an effective spring constant for the weak and strong axes and show that this is consistent with vortex pinning that is strong across TBs and weaker along them.


We thank J. Straquadine and E. Rosenberg, and C. Watson for experimental assistance and P. Massat for feedback on the manuscript. This work was supported by the Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, under Contract No. DE-AC02-76SF00515


apsrev4-2

