# Investigating the superconducting state of 2 H – NbS 2 as seen by the vortex lattice

A. Alshemi, E. Campillo, E. M. Forgan, R. Cubitt, M. Abdel-Hafiez, E. Blackburn

► To cite this version:

A. Alshemi, E. Campillo, E. M. Forgan, R. Cubitt, M. Abdel-Hafiez, et al. Investigating the superconducting state of 2 H – NbS 2 as seen by the vortex lattice. Physical Review Research, 2024, 6 (3), pp.033218. 10.1103/PhysRevResearch.6.033218. hal-05520013

HAL Id: hal-05520013

https://hal.science/hal-05520013v1

Submitted on 26 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Investigating the superconducting state of 2H-NbS₂ as seen by the vortex lattice

A. Alshemi ◦, ¹,* E. Campillo ◦, ¹ E. M. Forgan ◦, ² R. Cubitt ◦, ³ M. Abdel-Hafiez, ⁴,⁵ and E. Blackburn ◦, ¹

¹Division of Synchrotron Radiation Research, Department of Physics, Lund University, SE-22100 Lund, Sweden

²School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham, B15 2TT, United Kingdom

³Institut Laue Langevin, 71 Avenue des Martyrs, F-38000 Grenoble cedex 9, France

⁴Department of Applied Physics and Astronomy, University of Sharjah, P.O. Box 27272 Sharjah, United Arab Emirates

⁵Department of Physics and Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden

<figure><img src="image_2.png" /></figure>

(Received 30 April 2024; accepted 25 July 2024; published 26 August 2024)

2H-NbS₂ is a classic example of an anisotropic multiband superconductor, with significant recent work focusing on the interesting responses seen when high magnetic fields are applied precisely parallel to the hexagonal niobium planes. It is often contrasted with its sister compound 2H-NbSe₂ because they have similar onset temperatures for superconductivity, but 2H-NbS₂ has no charge density wave whereas in 2H-NbSe₂ the charge density wave order couples strongly to the superconductivity. Using small-angle neutron scattering, a bulk-sensitive probe, we have studied the vortex lattice and how it responds to the underlying superconducting anisotropy. This is done by controlling the orientation of the field with respect to the Nb planes. The superconducting anisotropy, Γ<sub>ac</sub> = 7.07 ± 0.2, is found to be field independent over the range measured (0.15 to 1.25 T), and the magnetic field distribution as a function of the applied magnetic field is found to be in excellent quantitative agreement with anisotropic London theory modified with a core-size cutoff correction, providing the first complete validation of this model. We find values of λ<sub>ab</sub> = 141.9 ± 1.5 nm for the in-plane London penetration depth, and λ<sub>c</sub> ~ 1 µm for the out-of-plane response. The field-independence indicates that we are primarily sampling the larger of the two gaps generating the superconductivity in this material.

DOI: 10.1103/PhysRevResearch.6.033218

I. INTRODUCTION

Recent studies of the superconducting phase diagram of 2H-NbS₂ [1] and 2H-NbSe₂ [2] have revealed that the superconducting state in these sister materials is very sensitive to the orientation of the magnetic field with respect to the basal plane. Different types of spatially textured superconductivity [3,4] have been conjectured at high magnetic fields in both materials.

These two compounds belong to the larger family of transition metal dichalcogenides (TMDCs). This is a class of highly two-dimensional (2D) materials, with hexagonal planes of transition metals weakly coupled along the **c** axis, with members displaying multiple types of electronic order, including charge density wave (CDW) order, Mott-insulating behavior and superconductivity [5]. 2H-NbSe₂ is a classic example of the interaction between a charge density wave state and superconductivity, and the family 2H-MX₂ (M = Nb, Ti, Ta; X = S, Se) all have similar electronic band structures in the normal state. When they become superconducting, typically a larger gap develops on Fermi surface sheets with a more 2D

nature, and the smaller gap appears on Fermi surface sheets with a more three-dimensional (3D) character. Within this family, 2H-NbS₂ stands out because it is the only one in which CDW order has not been seen in the bulk (see Table I); in the literature there is some dispute if it exists in the monolayer 1H-NbS₂ [6,7]. This means that comparing 2H-NbS₂ and 2H-NbSe₂ is a clean way to check out how the CDW affects the physics observed in these materials.

Both materials are considered to be excellent examples of two-band superconductors with two different s-wave gaps, from scanning tunneling spectroscopy [11], specific heat [12,13], and Andreev reflection [14] studies. For both materials, the Fermi surfaces are similar, with three types of Fermi surface sheet. The Fermi surface sheets that arise from the Nb 4d bands and exhibit superconductivity are cylinders centered around the **Γ** and **K** points in the Brillouin zone. They have different levels of corrugation leading to more 3D character in those around the **Γ** point. There is also a (non-superconducting) smaller pancakelike sheet at **Γ** associated with the chalcogen. This has been reported by many independent groups; a nice description is given by Noat et al. [15]. Where there are two superconducting sheets, one may expect different gap magnitudes and superconducting anisotropies, although the effective vortex core radii (related to the coherence lengths) are expected to lock together in most cases [16,17].

In this context, we use the vortex lattice (VL) that develops in the superconducting state to probe the superconducting response, with a particular eye on the effects of anisotropy. In 2H-NbSe₂, the vortex lattice has been studied using a

<footnote>*Contact author: ahmed.alshemi@sljus.lu.se</footnote>

<footnote>Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI. Funded
by Bibsam.</footnote>
<table><caption>TABLE I. Transition temperatures for the charge density wave and superconducting states in selected transition metal dichalcogenides. Values are taken from Refs. [8–10].</caption><thead><tr><th>TMDC</th><th><i>T</i><sub>CDW</sub> (K)</th><th><i>T</i><sub>c</sub> (K)</th></tr></thead><tbody><tr><td>2H-NbSe<sub>2</sub></td><td>33.5</td><td>7.3</td></tr><tr><td>2H-NbS<sub>2</sub></td><td>None</td><td>5.5</td></tr><tr><td>2H-TaSe<sub>2</sub></td><td>122.3</td><td>0.15</td></tr><tr><td>2H-TaS<sub>2</sub></td><td>78</td><td>0.8</td></tr></tbody></table>

variety of methods. Two results stand out. First, the individual vortex cores have a six-fold star-shaped structure [11], reflecting the symmetry of the CDW order. Second, the vortex lattice consists of one hexagonal domain; to a first approximation the vortices lie parallel to the external field [see Fig. 1(b)]. On rotating the magnetic field toward the basal plane, this domain distorts, reflecting the underlying effective mass anisotropy. However, the unit cell vectors of this domain do not change direction and are, in fact, pinned to the crystallographic *a*\* axis [18,19]. This unexpected observation

indicates that in 2H-NbSe₂ the orientation of the VL is not in agreement with the predictions of anisotropic London theory, as seen in, for example, YBa₂Cu₃O₇-δ [20], Sr₂RuO₄ [21], and KFe₂As₂ [22].

In contrast, there are relatively few direct studies of the vortex lattice in 2H-NbS₂, with a scanning tunneling microscopy and spectroscopy study by Guillamon *et al.* [11] confirming that a well-ordered vortex lattice can be seen at the sample surface, and that the vortex core has a standard circular shape.

Here, we present a study of the vortex lattice as measured deep inside the superconducting state at 1.5 K, extracting information on the superconducting anisotropy, the penetration depth λ, and the coherence length ξ. This is done by neutron diffraction from the magnetic field distribution associated with the vortex lattice, whereby the Fourier components of this periodic distribution can be extracted from the resulting Bragg reflections. A schematic of the experimental setup is shown in Fig. 1. This is done in the small-angle scattering regime because the inter-vortex distances are much larger than available neutron wavelengths. To probe the anisotropy of the system, we create vortex lattices with different angles Ω

<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematic diagram of the experimental setup. (a) Typical small-angle neutron scattering (SANS) instrument setup for diffraction by a superconducting vortex lattice. The neutrons pass through a velocity selector, which sets their average wavelength (usually with a full-width half-maximum (FWHM) spread ~10%). The beam is well collimated, with a long evacuated flight path before and after the sample to minimize air scattering. The transmitted main beam is caught on a neutron-absorbing beamstop to allow the Bragg reflections, which are scattered at small angles, to be visible on a 2D multidetector. The sample crystal structure is represented here, while the actual crystal orientation used during the experiment is shown in (c). (b) Sketch of the “parallel” field geometry used in the experiment; here the individual vortices are aligned parallel to the applied field, and then by rotating the magnet and sample together through angles ω or φ, the Bragg condition for the vortex lattice can be met. When the data obtained at different angles are summed together, an image such as that shown on the multidetector will be seen. (c) Experimental geometry. The coordinate system is defined with the *z* direction along **B**, and the components of observed VL scattering vectors are denoted as *q*<sub>x</sub> and *q*<sub>y</sub>. The applied magnetic field **B** is rotated away from the *c*-direction by an angle Ω, and the longitudinal and transverse components of the field modulation are denoted by *b*<sub>*z*</sub> and *b*<sub>*x*</sub> &amp; *b*<sub>*y*</sub>, respectively.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. SANS diffraction patterns of the VL in 2H-NbS₂ as a function of field rotation angle Ω. We show the results of measurements made at 1.5 K for seven different field angles, Ω, at either 0.15 T or 0.25 T, depending on the closeness of the Bragg spots to the direct beam. The images are sums over φ rocking scans about the horizontal axis perpendicular to the incoming beam, minus backgrounds. The direct beam has been masked off in the software. The white and green dashed lines in panels (a), (b), and (c) represent the opening angles of α and β for Domain I (●●●) and Domain II (■■■), respectively. The white dashed ellipse depicted in panel (d) lies on top of the Bragg reflections from the distorted VL. It has a major-to-minor axis ratio Γ<sub>VL</sub> and the reciprocal-space area of the ellipse is given by πq<sub>Δ</sub><sup>2</sup>, with q<sub>Δ</sub> = 2π(2B/√3Φ<sub>0</sub>)<sup>1/2</sup>.</figcaption></figure>

between the applied magnetic field and the **c** axis. Full details of the experimental protocols are given in Appendix A. Figure 2 shows a range of diffraction patterns collected at *T* = 1.5 K and magnetic fields ranging from 0.15 T to 0.25 T, for a range of values of *Ω*, running from 0° to 87.5°, prepared using the field cooling process described in Appendix A. It is clear that rotating the field from being parallel to **c** axis (*Ω* = 0°) to nearly being in the basal plane (*Ω* = 87.5°) distorts the hexagonal vortex lattice, indicating the anisotropic nature of 2H-NbS₂.

## II. RESULTS AND DISCUSSION

### A. Vortex lattice structure

When the field (*B* = 0.15 T) is parallel to **c** axis (*Ω* = 0°), the diffraction pattern is perfectly hexagonal [Fig. 2(a)], with the Bragg peaks appearing at *q*<sub>Δ</sub> = 5.6(1) × 10<sup>-3</sup> Å<sup>-1</sup>, as compared with the ideal value for a hexagonal VL *q*<sub>Δ</sub> = 2π(2B/√3Φ<sub>0</sub>)<sup>1/2</sup> = 5.7 × 10<sup>-3</sup> Å<sup>-1</sup> at *B* = 0.15 T. Φ<sub>0</sub> = *h*/2*e is the flux quantum. Because of the sixfold crystal symmetry, one of two hexagonal domains separated by 30° might be expected to be energetically favored. We observe that one domain (Hexagonal Domain I [●●●]) is dominant, with trace amounts of Hexagonal Domain II [■■■]. Prior to changing Ω, the diffraction spots of Domain I lie parallel to the **a**/*b* axes, indicating that, in real space, the vortex lattice planes are perpendicular to the Nb nearest neighbor directions. When Ω = 0°, the *q*<sub>x</sub> and *q*<sub>y</sub> directions on the detector are probing directions at right angles within the basal plane, and we expect the superconducting parameters (λ and ξ) to be essentially isotropic around this plane.

At Ω = 30° at 0.15 T, the ratio of the two domains has shifted, making Domain II easier to see, and we can also see the VL begin to distort due to the difference in the superconducting properties in- and out-of-plane. The Bragg peaks move further away from the beam center along the horizontal axis and closer along the vertical axis. The same distortion applies to both domains, as can be seen by considering that all of the Bragg spots lie on the same ellipse [Figs. 2(b) and 2(d)].

As the field is rotated further, Domain II becomes the dominant form and, indeed, is the only domain seen for Ω > 70°. As the distortion increases, the top and bottom reflections become more prominent. The side reflections also move out of the window captured by the 2D detector. The regions in which the different domains are seen as a function of magnetic field and Ω are shown in Fig. 3. No intermediate structures are seen at any point. These observations suggest that the transition between Domains I and II is first order in character, as reported in CaAlSi [23].

As Ω is increased, *q*<sub>y</sub> remains in the same orientation in the basal plane of the crystal lattice, but *q*<sub>x</sub> starts to mix basal plane and out-of-plane components. A VL structure with Bragg peaks pinned to the *q*<sub>y</sub> direction should then be preferred, as observed in, for example, YBa<sub>2</sub>Cu<sub>3</sub>O<sub>7-δ</sub> [20], Sr<sub>2</sub>RuO<sub>4</sub> [21], and KFe<sub>2</sub>As<sub>2</sub> [22]. However, as noted above, this was *not* observed in 2H-NbS<sub>2</sub>'s sister compound, 2H-NbSe<sub>2</sub>, presumably due to a stabilizing effect associated with the CDW formation [19].

### B. Superconducting anisotropy

The vortex lattice becomes increasingly anisotropic as the field is rotated toward the basal plane, because the London
<figure><img src="image_5.png" /></figure>

FIG. 4. Spin-splitting of the VL Bragg reflections due to spin-flip scattering. Vortex lattice diffraction patterns as a function of rotation angle (Ω) Vortex lattice rocking curves at 0.45 T and 1.5 K. Each rocking curve is fitted by two Lorentzian. The VL anisotropy increases as Ω comes close to the basal plane, as indicated by the dashed white lines in the upper panels (a–c). Lower panels (d–f) represent the rocking curves corresponding to the diffraction patterns in the upper panels. The rocking curves show the scattered intensity distribution plotted as a function of the tilt angle deviation (φ + Δφ) relative to the rocking curve center which has slight zero offset at φ = 0.09°. Two distinct peaks, indicative of Zeeman splitting from transverse field modulation (spin flip), the top Bragg reflection is represented in red, while the bottom Bragg reflection is depicted in purple. Both sets of peaks are fitted with two Lorentzians, as shown by the dashed lines, delineating the peak intensities and width for the top and bottom spots as a function of rotation angle (Ω). The scattering geometry triangles are depicted in panel (e) (with angles exaggerated for clarity) for the two different SF processes: Spin-up to spin-down (left-hand side) and spin-down to spin-up (right-hand side), The scattering angle (2θ = 2φ₀) is the same in both cases, but different tilt angles are required to satisfy the Bragg condition φ₂/₃ (kinetic energy loss) and φ₁/₄ (kinetic energy gain) for the top and bottom Bragg spots, respectively.

average of all methods at each condition is shown in Fig. 5. By
fitting to Eq. (1), the superconducting anisotropy is found to
be Γ<sub>ac</sub> = 7.07 ± 0.2.

<figure><img src="image_6.png" /></figure>

FIG. 5. Field-independent vortex lattice anisotropy. The VL anisotropy measured at 1.5 K as a function of the applied field and the angle between the field and the c axis (Ω). The dashed line shows the VL anisotropy calculated using Eq. (1), with Γ<sub>ac</sub> = 7.07; the solid line is for Γ<sub>ac</sub> = ∞.

For anisotropic superconductors, like 2H-NbS₂, the vari-
ations in the penetration depth ultimately arise due to
differences in the Fermi velocities and effective masses of
the carriers in the different directions. Γ<sub>ac</sub> can therefore
also be observed in the ratios of the superconducting coher-
ence length, ξ<sub>ab</sub>/ξ<sub>c</sub>, and the upper critical fields, B<sub>c2</sub><sup>ab</sup>/B<sub>c2</sub><sup>c</sup>.
The upper critical fields have been measured by transport
at 2 K, giving anisotropies ranging from 6.5 to 8.1 [27,28],
and a temperature-dependent study of the heat capacity by
Kačmarčík et al. [12] finds a value of 7 for temperatures
0.3 T<sub>c</sub> < T < T<sub>c</sub>. While this value is in agreement with ours,
Leroux et al. estimated a value of 11 from extrapolations of the
penetration depth to 0 K [29], while Cho et al. [24] used torque
magnetometry to trace out the upper critical fields, combined
with extremely precise in-plane angular alignment. From this
latter work, the in-plane upper critical field has an unusual
temperature dependence. Using the critical fields measured at
1.5 K gives B<sub>c2</sub><sup>ab</sup>/B<sub>c2</sub><sup>c</sup> = 7.33. However, this upper critical field
appears to be significantly lower than the orbital upper critical
field estimated from measurements close to T<sub>c</sub>, which would
give Γ<sub>ac</sub> ~ 15. Bi et al. [30] have pointed out that surface
superconductivity (up to B<sub>c3</sub> = 1.695B<sub>c2</sub>) may be playing a
role here.

In multiband superconductors, the different bands may
have different anisotropies, leading to field- and temperature-
dependent superconducting anisotropies; a classic example is MgB₂ [31]. 2H-NbS₂ is widely accepted to be a two-band superconductor, but we do not observe any field dependence in our superconducting anisotropy, as was the case in the heat capacity studies of Kačmarčík *et al.* [12]. To engage with this question further, we first need to extract more information about the characteristic superconducting lengths from our data.

### C. Integrated intensity and the field-dependent form factor

For a structurally two-dimensional, well-ordered VL, the spatial variation of the magnetic field in the mixed state can be written as a Fourier series with components at different momentum transfers. The Fourier component associated with the first-order Bragg reflections from the diffracted lattice is referred to as the “form factor” and can be related to the integrated intensity $I(q_i)$ of a Bragg reflection from VL domain i via the Christen formula [32]:

$$I(q_i) = 2\pi V_i S \left(\frac{\gamma}{4}\right)^2 \frac{\lambda_n^2}{\Phi_0^2 q_i \cos(\zeta)} |F(q_i)|^2. \quad (4)$$

Here, $V_i$ is the volume of the sample occupied by the VL domain i (with $V_1 + V_{II} = V$, the total sample volume). S is the incident neutron flux density (extracted through a measurement of the direct beam with known aperture size), $\lambda_n$ is the neutron wavelength, $\gamma$ is the gyromagnetic ratio of the neutron, $\Phi_0$ is the flux quantum, $q_i$ is the magnitude of the scattering vector for the relevant Bragg spot in the diffraction pattern, and $\cos(\zeta)$ is the Lorentz factor, where $\zeta$ is the angle between the reciprocal lattice vector $q_i$ and the direction that is at right angles to the rocking axis [33]. For each domain, the relevant integrated intensities of individual Bragg peaks (see Appendix A) are averaged to give $I(q_i)$ for that domain; this is related to the form factor via Eq. (4).

In an isotropic superconductor, the flux lines lie parallel to the applied field, with the screening supercurrents perfectly perpendicular to the direction of the average field, so that the local fields are all parallel to the applied field. The form factor of the vortex lattice therefore contains only components parallel to the field [b_z in Fig. 1(c)]. In 2H-NbS₂, this case applies if the field is perfectly parallel to the c axis. As soon as the field rotates away from this, transverse field components (b_{x,y}) will develop, leading to the spin-split scattering discussed above. This happens because the supercurrents tend to flow within the “easy” basal plane. Whereas the flux lines still follow the average field direction, the form factor arises from spatially varying contributions from both the longitudinal and transverse field components. This has been fully described using anisotropic London theory by Kogan [34] for uniaxial crystal systems like 2H-NbS₂, and expanded to account for the effective mass anisotropy of the carriers by Thiemann *et al.* [26]. Kealey *et al.* applied this practically to the biaxial superconductor YBa₂Cu₃O₇-δ [20], correcting some misprints in Ref. [26], but were only able to find qualitative agreement, as was the case in studies on Sr₂RuO₄ [21,35] and KFe₂As₂ [22]. This may be due to difficulties in handling field-dependent superconducting anisotropies. In 2H-NbS₂, we have a fixed value of $\Gamma_{ac}$, and we are able to validate the model of

Thiemann *et al.* over almost the entire angular range at multiple fields.

In the mixed state, the magnetic field distribution $B(r)$ can be decomposed into a Fourier series over the set of reciprocal space wave vectors $q$:

$$B(r) = \sum_{q} b(q) \exp(iq \cdot r). \quad (5)$$

The average field is in the z direction, parallel to the applied field, and the transverse fields $b_x(r)$ and $b_y(r)$ are in the xy plane. For a given Bragg reflection at $q = (q_x, q_y)$, the field components are [26]

$$b_x = \frac{(\lambda^2 m_{xz} q_y^2) B}{d}, \quad (6)$$

$$b_y = \frac{(-\lambda^2 m_{xz} q_x q_y) B}{d}, \quad (7)$$

$$b_z = \frac{(1 + \lambda^2 m_{zz} q^2) B}{d}, \quad (8)$$

where

$$d = (1 + \lambda^2 m_{yy} q_x^2 + \lambda^2 m_{xx} q_y^2)(1 + \lambda^2 m_{zz} q^2) - \lambda^4 m_{xz} q^2 q_y^2. \quad (9)$$

Here B is the average field and $\lambda = (\lambda_{ab}^2 \lambda_c)^{1/3}$ is the geometric mean of the penetration depths in the ab plane and along the c axis, and $m_{ij}$ are the tensor components of the effective mass of the charge carriers, referred to the x, y, and z axes defined in Fig. 1(c). The normalized effective masses along the unit cell axes are $m_a$, $m_b$, and $m_c$. For 2H-NbS₂, $m_a = m_b < m_c$, so $m_a^2 m_c = 1$. The superconducting anisotropy can be quantified as the ratio of the normalized effective masses, $\Gamma_{ac} = (m_c/m_a)^{1/2}$, so we can rewrite the effective mass components in the vortex frame into functions of $\Gamma_{ac}$ and $\Omega$:

$$m_{xx} = \Gamma_{ac}^{-2/3} \cos^2 \Omega + \Gamma_{ac}^{4/3} \sin^2 \Omega, \quad (10)$$

$$m_{yy} = \Gamma_{ac}^{-2/3}, \quad (11)$$

$$m_{zz} = \Gamma_{ac}^{-2/3} \sin^2 \Omega + \Gamma_{ac}^{4/3} \cos^2 \Omega, \quad (12)$$

$$m_{xz} = (\Gamma_{ac}^{-2/3} - \Gamma_{ac}^{4/3}) \sin \Omega \cos \Omega. \quad (13)$$

With these expressions we can then calculate the theoretical form factor, $F(q)$, for a given Bragg reflection, including a Gaussian cutoff term to account for the finite size of the vortex core [36]:

$$F(q) = (b_x^2 + b_y^2 + b_z^2)^{1/2} \exp[-c(q_x^2 \xi_{\perp}^2 + q_y^2 \xi_{ab}^2)]. \quad (14)$$

$\xi_{\perp} = [\xi_{ab}^2 \cos^2 \Omega + \xi_c^2 \sin^2 \Omega]$ represents the core width along the $q_x$ direction. Here, $\xi_{\perp}$ is $\xi_c$ when $\Omega = 90^{\circ}$, i.e., $B \parallel$ crystal planes. $\xi_{ab}$ represents the in-plane coherence length while $\xi_c$ is the coherence length for the c-axis direction. The constant c is a core cutoff parameter; a quantitative comparison of this model and the numerical solution of the Eilenberger equations indicates that the most suitable value is $c = 0.44$ [37].

The theoretical form factors depend on sample properties $\lambda$, $\xi_{ab}$, and $\xi_c$. We fitted the field-dependent data separately
<figure><img src="image_7.png" /><figcaption>FIG. 6. Field dependence of the vortex lattice form factor for different field-sample angle $\Omega$. Each panel shows the vortex lattice form factor as a function of the applied magnetic field $B$, measured at 1.5 K for a particular value of $\Omega$, ranging from 30° to 87.5°. For the upper panels, scattering from two distinct hexagonal vortex lattice domains was present, denoted as Domain I (●●) and Domain II (■■), while for the lower panels, the scattering arises only from Domain II. The data points represent experimental measurements of $|F(q_i)|^2$ with their respective error bars (see Supplemental Material Note 1 [48]). The dashed black lines are the results from fits to Eq. (14), using theoretical values for the scattering vectors $q$ of the Bragg peaks assuming $\Gamma_{ac} = 7$. Hexagons are depicted for each $\Omega$ to visualize the vortex lattice, constructed based on the actual anisotropy ratio $\Gamma_{VL}(\Omega)$.</figcaption></figure>

for different values of $\Omega$. The fits are very insensitive to the value of $\xi_c$, as it only influences the Domain I spots, and even then its contribution varies as $\xi_c^2 \sin^2 \Omega \ll \xi_{ab}^2 \cos^2 \Omega$, and so has little effect. Accordingly, we fixed this parameter as $\xi_c = \xi_{ab}/\Gamma_{ac}$, and let $\lambda$ and $\xi_{ab}$ vary as a fittable parameters. Where two VL domains were present, their data were fitted simultaneously, with one additional (field-dependent) fitting parameter representing the fractions of the sample volume occupied by the two domains: $V_I/V$ and $V_{II}/V$. In Fig. 6, we show the VL form factor variation as a function of field at $T = 1.5$ K for six values of $\Omega$. At each field, we distinguish between the form factor values obtained from different types of Bragg spots, as these different spots will have different amounts of longitudinal and transverse field, so that the Domain I spots typically have lower form factors than Domain II (illustrated more completely in Appendix C). Each value represents the average taken over equivalent spots for a given domain. The $(V_i/V)|F(q_i)|^2$ were calculated from the experimental data using Eq. (4) and then fitted to Eq. (14) using two different methods.

In the first approach, the intensities were fitted using the theoretical values for the scattering vector $q$ associated with the peak given that $\Gamma_{ac} = 7$. For Domain I(●●), $q_x = (q_Δ/2)(\Gamma_{VL})^{1/2}$ and $q_y = (q_Δ/2)(3/\Gamma_{VL})^{1/2}$, while for Domain II(■■) $q_x = 0$ and $q_y = (q_Δ)/(\Gamma_{VL})^{1/2}$. In the second approach, we used the experimentally measured values

for $q_x$ and $q_y$. These two approaches yielded statistically indistinguishable results, given the experimental errors (an example is shown in Appendix D). This consistency increases confidence in the robustness of our findings. In Fig. 6 only the first approach is presented, as smooth fit lines can be calculated from Eq. (14) using the theoretical Bragg spot positions.

The fits show excellent agreement with the description of the anisotropic superconductor given by Thiemann et al. [26], indicating that the tendency for supercurrents to flow in the ab-plane can be described in this way. The values obtained for $\lambda$ and $\xi$ are given in Table II.

<table><caption>TABLE II. Fitted penetration depth and coherence length values at each angle $\Omega$. These were obtained by fitting the field dependence of the form factor data measured at 1.5 K with the anisotropic London model with core-size correction factor described in Eq. (14).</caption><thead><tr><th>Ω (°)</th><th>λ<sub>GM</sub> (nm)</th><th>λ<sub>ab</sub> (nm)</th><th>λ<sub>c</sub> (nm)</th><th>ξ<sub>ab,eff</sub> (nm)</th></tr></thead><tbody><tr><td>30</td><td>277 ± 5</td><td>145 ± 3</td><td>1022 ± 28</td><td>19.5 ± 0.5</td></tr><tr><td>60</td><td>275 ± 3</td><td>144 ± 2</td><td>1014 ± 22</td><td>20.1 ± 0.2</td></tr><tr><td>70</td><td>270 ± 2</td><td>141 ± 1</td><td>995 ± 21</td><td>23.2 ± 0.4</td></tr><tr><td>80</td><td>273 ± 7</td><td>142 ± 4</td><td>1006 ± 32</td><td>34 ± 1</td></tr><tr><td>85</td><td>272 ± 10</td><td>142 ± 5</td><td>1001 ± 40</td><td>34 ± 1.5</td></tr><tr><td>87.5</td><td>273 ± 3</td><td>142 ± 1</td><td>1006 ± 21</td><td>24.6 ± 0.4</td></tr></tbody></table>
Considering $\xi$ first, ideally, the in-plane coherence length $\xi_{ab}$ should represent the core size of the vortices, in which case it can be related to the upper critical field if the superconductivity is destroyed by orbital overlap, using the Ginzburg-Landau expression $B_{c2}^{\parallel c}(0\,K) = \Phi_0/2\pi\xi_{ab}^2$. By extrapolating to 0 K, its value can be estimated from the zero-temperature upper critical field $B_{c2}^{\parallel c}$. This upper critical field is well documented for NbS$_2$, and for our sample, it is 1.8 T [14]. If this is the orbital limit, then $\xi_{ab}(0\,K) \sim 13.5$ nm. Our fit results vary between 19.5 and 34 nm in size for this parameter, and so we consider that our measured value is an effective value, rather than the intrinsic value. This could be because of flux pinning affecting the regularity of the vortex lattice. This can be represented by a 'static Debye-Waller factor' which arises from static disorder within the vortex lattice which can include local wiggling of the vortices [38], or zigzagging of the vortices between basal planes at higher values of $\Omega$. While we could not resolve $\xi_c$ in our fits, we note that using $B_{c2}^{\perp c}(0\,K) = \Phi_0/2\pi\xi_{ab}\xi_c$, if we take the orbital value for the upper critical field of 24 T, as calculated by Cho *et al.* [24] based on the Werthamer-Helfand-Hohenberg model, this gives $\xi_c \sim 1$ nm.

We now turn to the penetration depth, where the geometric mean at all values of $\Omega$ is found to be angle independent, as expected, with an average value of $272.3 \pm 1.3$ nm. From this, $\lambda_{ab} = 141.9 \pm 1.5$ nm and $\lambda_c = 1003 \pm 20$ nm. From the literature, there is sparse agreement: $\lambda_{ab}$ values between 83 nm [29] and 590 nm [39] are reported. On our samples, a value for $\lambda_{ab}$ of 131 nm was obtained by Majumdar *et al.* [14].

Recently, Kogan *et al.* [40] established a relationship between the zero-temperature penetration depth $\lambda(0)$ and the slope of the penetration depth $\lambda^{-2}(T)$ near $T_c$, using a calculation analogous to the Helfand-Werthamer relationship between the zero-temperature upper critical magnetic field and its slope at $T_c$. With further thermodynamic information, they obtained

$$ \lambda^2(0) \approx \left| \left( \frac{dH_{c2}}{dT} \right)_{T_c} \right| \frac{1}{T_c \gamma}, \quad (15) $$

where $(dH_{c2}/dT)_{T_c}$ is the slope of the upper critical magnetic field $H_{c2}$ with respect to temperature ($T$) at $T_c$, $T_c$ is the critical temperature of the superconductor in question and $\gamma$ is the specific heat coefficient per unit volume. Kogan *et al.'s* model has been developed for isotropic s-wave superconductivity with nonmagnetic scattering, and was successfully cross-checked against experimental values for $\lambda$ in $V_3Si$ and $Nb_3Sn$.

Using this model we have calculated $\lambda_{ab}(0) \approx 142.9$ nm for $2H-NbS_2$, taking as inputs $T_c = 5.5$ K, $(dH_{c2}^c/dT)_{T_c} \approx -0.25 \times 10^4$ Oe/K [12,24,28], and $\gamma = 0.494 \times 10^4$ erg/cm$^3$K$^2$ [41]. This is in excellent agreement with our $\lambda_{ab}$ (1.5 K) as extracted from the fit.

On the face of it, this is surprising, as $2H-NbS_2$ is an anisotropic two-band superconductor, which also makes the field-independent superconducting anisotropy unexpected. However, as Kogan *et al.* noted, their results should still be applicable if the order parameter is constant over a Fermi surface of any given shape. Indeed, Kogan *et al.* tested their

approach successfully on the two-band superconductor $MgB_2$, calculating $\lambda(0) \approx 176$ nm, as compared to reported values of 180–185 nm.

We therefore suggest that over the field range explored here at 1.5 K, there are no major changes in the contributions from the two gaps, with the dominant response coming from the large-gap band in $2H-NbS_2$, which has predominantly 2D character, as proposed by Kačmarčík *et al.* [12]. Interestingly, Noat *et al.* [15] have suggested that, in reality, only the large gap sheets (the cylinders centered on the K points) are intrinsically superconducting and that the other gap develops parasitically by coupling through the pancakelike sheet coming from the chalcogen p-bands.

### III. CONCLUSIONS

We have measured neutron diffraction by the vortex lattice in $2H-NbS_2$ as a function of field angle and field magnitude at 1.5 K, and find that in fields up to 1.25 T, the intensity of the observed diffraction from the vortex lattice can be described extremely well by the Thiemann model using anisotropic London theory [26], with the addition of a core-correction factor to account for the finite size of the vortex core. This forms the first full validation of this model. This process is aided by our clear observation that the superconducting anisotropy is constant and field-independent at $\Gamma_{ac} \sim 7$. From this, we extract values for the London penetration depth of $\lambda_{ab} = 141.9 \pm 1.5$ nm and $\lambda_c \sim 1$ µm. For $\lambda_{ab}$, this experimental result fits with the recent model developed by Kogan *et al.* [40].

The dominant Fermi sheets in these conditions are expected to be the cylinders around the K point in reciprocal space, which hold most of the density of states at the Fermi surface, and are highly two-dimensional. Coupling between these sheets and the more three-dimensional cylinders around the $\Gamma$ point, therefore, controls the overall superconducting anisotropy. The interband coupling for $2H-NbS_2$ has been found by Noat *et al.* [15] to be half that in $2H-NbSe_2$, from fitting tunneling spectroscopy data. It is therefore not surprising that $2H-NbS_2$ shows a stronger 2D character with a superconducting anisotropy of $7.07 \pm 0.2$, as compared to the value of $3.2 \pm 0.2$ measured by Gammel *et al.* [19] for $2H-NbSe_2$ using the same method presented here.

The experimental data collected at the Institut Laue-Langevin is publicly available [47]. Processed experimental data are available upon request to the corresponding authors.

### ACKNOWLEDGMENTS

A.A. and E.B. acknowledge support from the Crafoordska Stiftelsen (Grant No. 20190930) and the Swedish Research Council (Vetenskapsrådet) under Project No. 201804704. A.A., E.M.F., and E.B. thank the Institut Laue-Langevin for travel support to attend the experiment. The authors gratefully acknowledge the Institut Laue-Langevin for the allocated beamtime.

The authors have no conflicts to disclose.
## APPENDIX A: METHODS

Thin platelets of 2H-NbS₂ samples were grown using the chemical vapor transport technique, yielding high-quality single crystals with optically flat surfaces on the macroscopic scale. A detailed description of this process for growing TMD crystals is given by Chareev *et al.* [42]. The sample quality was confirmed by specific heat measurements at zero field, wherein a sharp jump centered at the superconducting transition temperature *T*<sub>c</sub> = 5.5 K is seen. This is considered to be the most reliable way to check sample quality [43], and matches well with the behavior reported for other samples grown in the same laboratory [1].

The 2H-NbS₂ polytype is hexagonal, with space group *P*6₃/mmc and lattice parameters *a* = *b* ~ 3.31 Å and *c* ~ 11.86 Å [43]. The platelets grow with the *c* axis normal to the platelet surface. A mosaic of nine coaligned crystals was made by mounting the platelets onto two aluminum sheets, so that all the crystal faces lay within an area 10 × 10 mm². The mosaic had a total crystal thickness of 80 µm with a total mass of 42.3 mg. The in-plane alignment (along the *a*\* and *b*\* axes) was carried out using an optical microscope, with the clearly visible crystal hexagonal facets serving as points of reference. Overall, the misalignment of the *a*\*/*b*\* axes was measured to be less than 1°.

To diffract off the vortex lattice, small angle neutron scattering measurements were carried out at the D33 instrument [44] at the Institut Laue-Langevin, Grenoble, France. The incident neutrons had wavelength $\lambda_n = 10$ Å, collimation of 10.3 m and a wavelength spread $\Delta\lambda_n/\lambda_n = 10\%$ full-width half-maximum (FWHM). The scattering patterns were collected on a (256 × 128 pixels) two-dimensional position-sensitive multidetector placed 10.0435 m after the sample, which was mounted in a 9 T horizontal-field cryomagnet. The crystals were aligned with the *a*\* axis vertical and the *b*\* axis horizontal. A schematic of the arrangement of a typical SANS instrument when used for VL studies is shown in Fig. 1(a).

To create the vortex lattice, a magnetic field was applied. As shown in Fig. 1(b), the magnetic field was applied essentially parallel to the incoming neutron beam. The advantage of this parallel field geometry is that by rocking the sample and magnet together through small rocking angles about axes perpendicular to the neutron beam, denoted $\phi$ and $\omega$, the diffraction conditions for all of the VL Bragg reflections can be met to give diffraction peaks on the 2D detector. In the work presented here, the sample and magnet were rocked about the horizontal axis perpendicular to the beam direction, i.e., the $\phi$ angle. At $\Omega = 0°$, the sample and magnet were also rocked along the $\Omega$ angle.

To study the vortex lattice through the entire bulk of the sample, a relatively well-ordered vortex lattice needs to be prepared. To do this, the sample was field cooled from above $T_c$ to $T = 1.5$ K while oscillating the magnitude of the applied magnetic field by $\pm 1\%$ about the desired value. This procedure improves the orientational ordering of the VL, by keeping the vortices away from local pinning potentials [45]. This was done for magnetic fields over the range 0.15 to 1.25 T.

To probe the superconducting anisotropy of 2H-NbS₂, we also rotated the magnetic field away from the *c* axis toward

the *b* direction by an angle $\Omega$. These rotations were all carried out above $T_c$. The experimental geometry is shown schematically in Fig. 1(c). As $\Omega$ is changed, the profile of the sample with respect to the beam changes, and so different sample apertures can be used at different angles; ideally only the sample volume should be illuminated by the neutron beam. For $\Omega = 0°, 30°, 80°, \text{and } 85°$, a 12 mm diameter circular aperture was used. For $\Omega = 60°$ and $70°$, a 7H × 10V mm² rectangular aperture was used. For $\Omega = 87.5°$, the aperture size was 3H × 10V mm².

Background data were collected in the normal state at 6.5 K, using the same rotation and tilt angles as those of the “foreground” measurements at 1.5 K. The background was then subtracted from the foreground, leaving only the vortex lattice signal. The analysis was done using the software program GRASP [46]. Experimentally, the integrated intensity of a Bragg peak from the VL is determined by counting the number of neutrons detected in the region of the detector where the peak appears as a function of the rocking angle. For monitoring the form factor, only $\phi$ rocking scans were used, so the Bragg condition for spots close to the horizontal axis is not fully met. This means that for Domain I, the integrated intensities of the four diagonal Bragg spots, represented by $\odot\odot$, are used, whereas for Domain II, only the top and bottom Bragg spots, $\odot$, are included. For the form factor calculations, the correct region of the detector to include for each peak was tested iteratively to maximize the signal-to-noise ratio after the background subtraction. The resulting rocking curves were fitted to a Lorentzian function with the background set to zero to give $I(q_i)$. For $\Omega \le 70°$, each Bragg spot gave one maximum in the rocking curve. At higher angles, the spin-splitting effect led to two maxima in each rocking curve. The integrated intensity in these cases consists of the sum of the intensities obtained from the two peaks, since each represents half of the incident neutron flux, corresponding to one direction of the neutron spin.

## APPENDIX B: DERIVATION OF EXPRESSION FOR VORTEX LATTICE ANISOTROPY $\Gamma_{VL}$ IN TERMS OF INTRINSIC SUPERCONDUCTING ANISOTROPY $\Gamma_{ac}$ AND ANGLE $\Omega$ BETWEEN B AND C AXIS

In the isotropic case where **B** $\parallel$ **c**, the nearest neighbor flux lattice vectors form a regular hexagon, sitting on a circle centered at the origin. When the direction of **B** is tilted away from the **c** axis, all *x* components of the flux lattice vectors are multiplied by a factor $\gamma$ ($\gamma < 1$ for planar materials such as 2H-NbS₂), while the *y* components are divided by $\gamma$ (Eq. (7) in Ref. [25]). The flux lattice vectors then sit on an ellipse centered at the origin with an axial ratio $\Gamma_{VL} = 1/\gamma^2$. We note that for a given value of *B*, the area of the ellipse is the same as that of the circle, due to flux quantization [25].

The first-order reciprocal lattice vectors (the measured diffraction spots) will be similarly distorted and sit on an ellipse centered about the origin that is also described by $\Gamma_{VL}$, although that ellipse will be rotated by 90° about the **B** direction.
Following Eq. (9) from Campbell et al. [25],

$$ \frac{1}{\gamma^2} = \left( \frac{m_3}{m_{zz}} \right)^{1/2} = \Gamma_{VL}. \qquad (B1) $$

In our notation, Eq. (2) in Ref. [25] can then be reformulated as

$$ m_{zz} = m_1 \sin^2 \Omega + m_3 \cos^2 \Omega, \qquad (B2) $$

which allows us to express $\Gamma_{VL}$ as

$$ \Gamma_{VL} = \left( \frac{m_3}{m_1 \sin^2 \Omega + m_3 \cos^2 \Omega} \right)^{\frac{1}{2}}. \qquad (B3) $$

Now, by introducing the parameter to describe the intrinsic superconducting anisotropy, $\Gamma_{ac}$, defined as $\Gamma_{ac} = (m_3/m_1)^{1/2}$, the final relationship linking the vortex lattice anisotropy with the intrinsic superconductor anisotropy is

$$ \Gamma_{VL} = \frac{\Gamma_{ac}}{\sqrt{\sin^2 \Omega + \Gamma_{ac}^2 \cos^2 \Omega}}. \qquad (B4) $$

## APPENDIX C: LONGITUDINAL AND TRANSVERSE MAGNETIC FIELD COMPONENTS OF THE FLUX LINE LATTICE IN THE ANISOTROPIC 2H-NbS$_2$

As indicated in Fig. 6, the Bragg reflections associated with Domain I have noticeably lower form factors than those for Domain II. This is due to the fact that there are two effects of the angular rotation on the VL form factor that bring this about. First, the penetration depth is changing as $\Omega$ changes. Hence, the form factors of the observed Bragg spots on the left and right in the SANS patterns for Domain I and the diagonal spots for Domain II both decrease due to the longer due to the longer penetration depth for currents in the c direction. However, this effect is insignificant for the top and bottom spots, whose form factor is primarily influenced by currents in the basal plane and remains relatively stable with small field rotations away from the c-axis. The second effect comes from $B_{c2}$, as the field is rotated away from the c axis, the upper critical field $B_{c2}$ increases since it is larger when the field is closer to the basal plane. Consequently, a given magnetic field strength represents a smaller fraction of $B_{c2}$, reducing the core overlap effects in the superconductor. This reduction in core overlap is uniform for all diffraction spots in the SANS pattern. The combined result of these changes is an increase in the form factor for the top and bottom spots, making them more intense than the left and right spots, which diminish in intensity upon rotation away from the c-axis. This specific outcome is notably observed at rotation angles of 80°, 85°, and 87.5° as shown in Fig. 2. In Fig. 7 we show how the form factor in the form of transverse and longitudinal field components $b_x$, $b_y$, and $b_z$ is changing as a function of both rotation angle ($\Omega$) and applied magnetic field ($B$). The calculations involve the diagonal spots of Domain I(●) and only the top and bottom spots of Domain II(●) using $\Gamma_{ac} = 7$ and Eqs. (6)–(8).

<figure><img src="image_8.png" /><figcaption>FIG. 7. Magnetic field components $b_x$, $b_y$, and $b_z$ generating the form factor from the vortex lattice at the Bragg peaks ● from Domain I and ● from Domain II, as a function of angle $\Omega$ taken at 1.5 K for fields 0.15 and 1.25 T.</figcaption></figure>

## APPENDIX D: FIELD-DEPENDENT FORM FACTOR FIT USING THE EXPERIMENTAL VALUES OF $q_x$ AND $q_y$ (APPROACH II)

In Fig. 6, we illustrated the fit of the form factor data using the theoretical definition for $q_x$ and $q_y$ (Approach I). Here, we show our Approach II by fitting the data with the experimentally measured values for $q_x$ and $q_y$. As shown in Fig. 8, we show the fit for angles $\Omega = 80°, 85°, \text{ and } 87.5°$ as an example. The red data points are the results from the fitting process using Approach II for each angle $\Omega$. The black
<figure><img src="image_9.png" /><figcaption>FIG. 8. Vortex lattice form factor for angles $\Omega = 80^{\circ}$, $85^{\circ}$, and $87.5^{\circ}$ as a function of the field at 1.5 K.</figcaption></figure>

dashed lines are from the Approach I fit shown in Fig. 6. This
shows that both fit procedures are robust and consistent, and

the output fit from both approaches yields the same $\lambda_{\text{GM}}$, $\xi_{ab}$, and $\xi_c$ within the experimental error.

[1] C. Cho, J. Lyu, C. Y. Ng, J. J. He, K. T. Lo, D. Chareev, T. A. Abdel-Baset, M. Abdel-Hafiez, and R. Lortz, Evidence for the Fulde-Ferrell-Larkin-Ovchinnikov in bulk NbS₂, Nat. Commun. **12**, 3676 (2021).

[2] P. Wan, O. Zheliuk, N. F. Q. Yuan, X. Peng, L. Zhang, M. Liang, U. Zeitler, S. Wiedmann, N. E. Hussey, T. T. M. Palstra, and J. Ye, Orbital Fulde-Ferrell-Larkin-Ovchinnikov state in an Ising superconductor, Nature (London) **619**, 46 (2023).

[3] P. Fulde and R. A. Ferrell, Superconductivity in a strong spin-exchange field, Phys. Rev. **135**, A550 (1964).

[4] A. Larkin and Y. N. Ovchinnikov, Inhomogeneous state of superconductors, Sov. Phys. JETP **20**, 762 (1965).

[5] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, 2D transition metal dichalcogenides, Nat. Rev. Mater. **2**, 17033 (2017).

[6] T. Knispel, J. Berges, A. Schobert, E. G. C. P. van Loon, W. Jolie, T. Wehling, T. Michely, and J. Fischer, Unconventional charge-density-wave gap in monolayer NbS₂, Nano Lett. **24**, 1045 (2024).

[7] D. Lin, S. Li, J. Wen, H. Berger, L. Forró, H. Zhou, S. Jia, T. Taniguchi, K. Watanabe, X. Xi, and M. S. Bahramy, Patterns and driving forces of dimensionality-dependent charge density waves in 2H-type transition metal dichalcogenides, Nat. Commun. **11**, 2406 (2020).

[8] D. E. Moncton, J. D. Axe, and F. J. DiSalvo, Neutron scattering study of the charge-density wave transitions in 2H-TaSe₂ and 2H-NbSe₂, Phys. Rev. B **16**, 801 (1977).

[9] J. A. Wilson, F. J. Di Salvo, and S. Mahajan, Charge-density waves and superlattices in the metallic layered transition metal dichalcogenides, Adv. Phys. **24**, 117 (1975).

[10] M. Abdel-Hafiez, X.-M. Zhao, A. A. Kordyuk, Y.-W. Fang, B. Pan, Z. He, C.-G. Duan, J. Zhao, and X.-J. Chen, Enhancement of superconductivity under pressure and the magnetic phase diagram of tantalum disulfide single crystals, Sci. Rep. **6**, 31824 (2016).

[11] I. Guillamon, H. Suderow, S. Vieira, L. Cario, P. Diener, and P. Rodière, Superconducting density of states and vortex cores of 2H-NbS₂, Phys. Rev. Lett. **101**, 166407 (2008).

[12] J. Kačmarčík, Z. Pribulová, C. Marcenat, T. Klein, P. Rodière, L. Cario, and P. Samuely, Specific heat measurements of

a superconducting NbS₂ single crystal in an external magnetic field: Energy gap structure, Phys. Rev. B **82**, 014518 (2010).

[13] N. Kobayashi, K. Noto, and Y. Muto, Thermodynamic properties of the layered superconductor 2H-NbSe₂, J. Low Temp. Phys. **27**, 217 (1977).

[14] A. Majumdar, D. VanGennep, J. Brisbois, D. Chareev, A. V. Sadakov, A. S. Usoltsev, M. Mito, A. V. Silhanek, T. Sarkar, A. Hassan, O. Karis, R. Ahuja, and M. Abdel-Hafiez, Interplay of charge density wave and multiband superconductivity in layered quasi-two-dimensional materials: The case of 2H-NbS₂ and 2H-NbSe₂, Phys. Rev. Mater. **4**, 084005 (2020).

[15] Y. Noat, J. A. Silva-Guillén, T. Cren, V. Cherkez, C. Brun, S. Pons, F. Debontridder, D. Roditchev, W. Sacks, L. Cario, P. Ordejón, A. García, and E. Canadell, Quasiparticle spectra of 2H - NbSe₂: Two-band superconductivity and the role of tunneling selectivity, Phys. Rev. B **92**, 134510 (2015).

[16] M. Ichioka, V. G. Kogan, and J. Schmalian, Locking of length scales in two-band superconductors, Phys. Rev. B **95**, 064512 (2017).

[17] M. Silaev and E. Babaev, Microscopic theory of type-1.5 superconductivity in multiband systems, Phys. Rev. B **84**, 094515 (2011).

[18] H. F. Hess, C. A. Murray, and J. V. Waszczak, Flux lattice and vortex structure in 2H-NbSe₂ in inclined fields, Phys. Rev. B **50**, 16528 (1994).

[19] P. L. Gammel, D. A. Huse, R. N. Kleiman, B. Batlogg, C. S. Oglesby, E. Bucher, D. J. Bishop, T. E. Mason, and K. Mortensen, Small angle neutron scattering study of the magnetic flux-line lattice in single crystal 2H-NbSe₂, Phys. Rev. Lett. **72**, 278 (1994).

[20] P. G. Kealey, D. Charalambous, E. M. Forgan, S. L. Lee, S. T. Johnson, P. Schleger, R. Cubitt, D. M. Paul, C. M. Aegerter, S. Tajima, and A. Rykov, Transverse-field components of the flux-line lattice in the anisotropic superconductor YBa₂Cu₃O₇₋₈, Phys. Rev. B **64**, 174501 (2001).

[21] C. Rastovski, C. D. Dewhurst, W. J. Gannon, D. C. Peets, H. Takatsu, Y. Maeno, M. Ichioka, K. Machida, and M. R. Eskildsen, Anisotropy of the superconducting state in Sr₂RuO₄, Phys. Rev. Lett. **111**, 087003 (2013).
[22] S. J. Kuhn, H. Kawano-Furukawa, E. Jellyman, R. Riyat, E. M. Forgan, M. Ono, K. Kihou, C. H. Lee, F. Hardy, P. Adelmann, T. Wolf, C. Meingast, J. Gavilano, and M. R. Eskildsen, Si-multaneous evidence for pauli paramagnetic effects and multi-band superconductivity in KFe₂As₂ by small-angle neutron scattering studies of the vortex lattice, Phys. Rev. B **93**, 104527 (2016).

[23] P. K. Biswas, M. R. Lees, G. Balakrishnan, D. Q. Liao, D. S. Keeble, J. L. Gavilano, N. Egetenmeyer, C. D. Dewhurst, and D. M. Paul, First-order reorientation transition of the flux-line lattice in CaAlSi, Phys. Rev. Lett. **108**, 077001 (2012).

[24] Chang-woo Cho, C. Y. Ng, C. H. Wong, M. Abdel-Hafiez, A. N. Vasiliev, D. A. Chareev, A. G. Lebed, and R. Lortz, Competition between orbital effects, pauli limiting, and Fulde-Ferrell-Larkin-Ovchinnikov states in 2D transition metal dichalcogenide superconductors, New J. Phys. **24**, 083001 (2022).

[25] L. J. Campbell, M. M. Doria, and V. G. Kogan, Vortex lattice structures in uniaxial superconductors, Phys. Rev. B **38**, 2439 (1988).

[26] S. L. Thiemann, Z. Radović, and V. G. Kogan, Field structure of vortex lattices in uniaxial superconductors, Phys. Rev. B **39**, 11406 (1989).

[27] K. Onabe, M. Naito, and S. Tanaka, Anisotropy of upper critical field in superconducting 2H-NbS₂, J. Phys. Soc. Jpn. **45**, 50 (1978).

[28] D. Huang, H. Nakamura, K. Küster, U. Wedig, N. B. M. Schröter, V. N. Strocov, U. Starke, and H. Takagi, Probing the interlayer coupling in 2H-NbS₂ via soft x-ray angle-resolved photoemission spectroscopy, Phys. Rev. B **105**, 245145 (2022).

[29] M. Leroux, P. Rodière, L. Cario, and T. Klein, Anisotropy and temperature dependence of the first critical field in 2H-NbS₂, Phys. B: Condens. Matter **407**, 1813 (2012).

[30] X. Bi, Z. Li, J. Huang, F. Qin, C. Zhang, Z. Xu, L. Zhou, M. Tang, C. Qiu, P. Tang, T. Ideue, T. Nojima, Y. Iwasa, and H. Yuan, Orbital-selective two-dimensional superconductivity in 2H-NbS₂, Phys. Rev. Res. **4**, 013188 (2022).

[31] R. Cubitt, M. R. Eskildsen, C. D. Dewhurst, J. Jun, S. M. Kazakov, and J. Karpinski, Effects of two-band superconductivity on the flux-line lattice in magnesium diboride, Phys. Rev. Lett. **91**, 047002 (2003).

[32] D. K. Christen, F. Tasset, S. Spooner, and H. A. Mook, Study of the intermediate mixed state of niobium by small-angle neutron scattering, Phys. Rev. B **15**, 4506 (1977).

[33] G. L. Squires, *Introduction to the Theory of Thermal Neutron Scattering*, 3rd ed. (Cambridge University Press, Cambridge, UK, 2012)

[34] V. G. Kogan, London approach to anisotropic type-II superconductors, Phys. Rev. B **24**, 1572 (1981).

[35] S. J. Kuhn, W. Morgenlander, E. R. Louden, C. Rastovski, W. J. Gannon, H. Takatsu, D. C. Peets, Y. Maeno, C. D. Dewhurst, J. Gavilano, and M. R. Eskildsen, Anisotropy and multiband superconductivity in Sr₂RuO₄ determined by small-angle neutron scattering studies of the vortex lattice, Phys. Rev. B **96**, 174507 (2017).

[36] A. Yaouanc, P. Dalmas de Réotier, and E. H. Brandt, Effect of the vortex core on the magnetic field in hard superconductors, Phys. Rev. B **55**, 11107 (1997).

[37] E. Campillo, M. Bartkowiak, R. Riyat, E. Jellyman, A. S. Cameron, A. T. Holmes, O. Prokhenko, W.-D. Stein, A. Erb, E. M. Forgan, and E. Blackburn, Deviations from the extended london model at high magnetic fields in YBa₂Cu₃O₇, Phys. Rev. B **105**, 184508 (2022).

[38] M. Tinkham, *Introduction to Superconductivity* (McGraw-Hill, New York, NY, 1996)

[39] A. Tulapurkar, A. Grover, S. Ramakrishnan, A. Niazi, and A. Rastogi, Peak effect studies in 2H-NbS₂, Phys. B: Condens. Matter **312-313**, 118 (2002).

[40] V. G. Kogan, M. A. Tanatar, and R. Prozorov, London penetration depth at zero temperature and near the superconducting transition, Phys. Rev. B **101**, 104510 (2020).

[41] Y. Hamaue and R. Aoki, Effects of organic intercalation on lattice vibrations and superconducting properties of 2H-NbS₂, J. Phys. Soc. Jpn. **55**, 1327 (1986).

[42] D. A. Chareev, P. Evstigneeva, D. Phuyal, G. J. Man, H. Rensmo, A. N. Vasiliev, and M. Abdel-Hafiez, Growth of transition-metal dichalcogenides by solvent evaporation technique, Cryst. Growth Des. **20**, 6930 (2020).

[43] C. Witteveen, K. Gornicka, J. Chang, M. Månsson, T. Klimczuk, and F. O. von Rohr, Polytypism and superconductivity in the NbS₂ system, Dalton Trans. **50**, 3216 (2021).

[44] C. D. Dewhurst, I. Grillo, D. Honecker, M. Bonnaud, M. Jacques, C. Amrouni, A. Perillo-Marcone, G. Manzin, and R. Cubitt, The small-angle neutron scattering instrument D33 at the Institut Laue-Langevin, J. Appl. Crystallogr. **49**, 1 (2016).

[45] M. Marziali Bermúdez, M. R. Eskildsen, M. Bartkowiak, G. Nagy, V. Bekeris, and G. Pasquini, Dynamic reorganization of vortex matter into partially disordered lattices, Phys. Rev. Lett. **115**, 067001 (2015).

[46] C. D. Dewhurst, Graphical reduction and analysis small-angle neutron scattering program: GRASP, J. Appl. Crystallogr. **56**, 1595 (2023).

[47] A. Alshemi, E. Blackburn, E. Campillo, R. Cubitt, and E. Forgan, Investigation of the vortex lattice in NbS₂—A potential FFLO candidate, Institut Laue-Langevin (ILL) (2021), https://doi.ill.fr/10.5291/ILL-DATA.5-42-558.

[48] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.6.033218 for experimental form factor values used in the fitting shown in Fig. 6.