# Investigating the superconducting state of 2$H$-NbS2 as seen by the vortex
lattice

A. Alshemi<EMAIL_ADDRESS>Division of Synchrotron Radiation
Research, Department of Physics, Lund University, SE-22100 Lund, Sweden  E.
Campillo  Division of Synchrotron Radiation Research, Department of Physics,
Lund University, SE-22100 Lund, Sweden  E. M. Forgan  School of Physics and
Astronomy, University of Birmingham, Edgbaston, Birmingham, B15 2TT, UK  R.
Cubitt  Institut Laue Langevin, 71 Avenue des Martyrs, F-38000 Grenoble cedex
9, France  M. Abdel-Hafiez  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  E. Blackburn  Division of Synchrotron Radiation Research,
Department of Physics, Lund University, SE-22100 Lund, Sweden

###### Abstract

2$H$-NbS2 is a classic example of an anisotropic multi-band superconductor,
with significant recent work focussing 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 2$H$-NbSe2 because
they have similar onset temperatures for superconductivity, but 2$H$-NbS2 has
no charge density wave whereas in 2$H$-NbSe2 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, $\Gamma_{ac}=7.07\pm 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 cut-off correction, providing the first complete validation of
this model. We find values of $\lambda_{ab}=141.9\pm 1.5$ nm for the in-plane
London penetration depth, and $\lambda_{c}\sim$ 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.

superconductivity, Vortex lattice, flux lines

## Introduction

Recent studies of the superconducting phase diagram of 2$H$-NbS2 Cho  _et al._
(2021) and 2$H$-NbSe2 Wan  _et al._ (2023) 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 Fulde and Ferrell (1964); Larkin
and Ovchinnikov (1965) 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 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 behaviour and superconductivity
Manzeli  _et al._ (2017). 2$H$-NbSe2 is a classic example of the interaction
between a charge density wave state and superconductivity, and the family
2$H$-MX2 (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 3D character. Within
this family, 2$H$-NbS2 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 1$H$-NbS2 Knispel  _et al._ (2024);
Lin  _et al._ (2020). This means that comparing 2$H$-NbS2 and 2$H$-NbSe2 is a
clean way to check out how the CDW affects the physics observed in these
materials.

Figure 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 $\sim 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. (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 $\Omega$ or $\phi$, 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 multi-detector 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_{x}$ and $q_{y}$. The applied magnetic field B is rotated away from the 2$H$-NbS2 plane by an angle $\Omega$ and the longitudinal and transverse components of the field modulation are denoted by $b_{z}$ and $b_{x}~{}\&~{}b_{y}$, respectively. TMDC |  $T_{\mathrm{CDW}}$ (K) |  $T_{c}$ (K)
---|---|---
2$H$-NbSe2 |  33.5 |  7.3
2$H$-NbS2 |  none |  5.5
2$H$-TaSe2 |  122.3 |  0.15
2$H$-TaS2 |  78 |  0.8

Table 1: Transition temperatures for the charge density wave and
superconducting states in selected transition metal dichalcogenides. Values
are taken from Refs. Moncton  _et al._ (1977); Wilson  _et al._ (1975); Abdel-
Hafiez  _et al._ (2016).

Both materials are considered to be excellent examples of two-band
superconductors with two different $s$-wave gaps, from scanning tunnelling
spectroscopy Guillamon  _et al._ (2008), specific heat Kačmarčík  _et al._
(2010); Kobayashi  _et al._ (1977) and Andreev reflection Majumdar  _et al._
(2020) 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
4$d$ bands and exhibit superconductivity are cylinders centred around the
$\Gamma$ and K points in the Brillouin zone. They have different levels of
corrugation leading to more 3D character in those around the $\Gamma$ point.
There is also a (non-superconducting) smaller pancake-like sheet at $\Gamma$
associated with the chalcogen. This has been reported by many independent
groups; a nice description is given by Noat _et al._ Noat  _et al._ (2015).
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 Ichioka  _et al._ (2017).

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 2$H$-NbSe2, the vortex lattice has been
studied using a variety of methods. Two results stand out. Firstly, the
individual vortex cores have a six-fold star-shaped structure Guillamon  _et
al._ (2008), reflecting the symmetry of the CDW order. Secondly, the vortex
lattice consists of one hexagonal domain; to a first approximation the
vortices lie parallel to the external field (see Fig. 1b). On rotating the
magnetic field towards 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 $\textbf{a}^{*}$ axis Hess  _et al._ (1994); Gammel  _et al._
(1994). This unexpected observation indicates that in 2$H$-NbSe2 the
orientation of the VL is not in agreement with the predictions of anisotropic
London theory, as seen in, for example, YBa2Cu3O7-δ Kealey  _et al._ (2001),
Sr2RuO4 Rastovski  _et al._ (2013) and KFe2As2 Kuhn  _et al._ (2016).

In contrast, there are relatively few direct studies of the vortex lattice in
2$H$-NbS2, with a scanning tunneling microscopy and spectroscopy study by
Guillamon et al. Guillamon  _et al._ (2008) confirming that a well-ordered
vortex lattice can be seen at the sample surface, and that the vortex core has
a standard circular shape.

Figure 2: SANS diffraction pattern of the VL in 2$H$-NbS2 as a function of
field rotation angle $\boldsymbol{\Omega}$. We show the results of
measurements made at 1.5 K for seven different field angles, $\Omega$, 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 $\phi$ rocking scans about the
horizontal axis perpendicular to the incoming beam, minus backgrounds. The
direct beam has been masked off in software. The white and green dashed lines
in panels (a), (b) and (c) represent the opening angles of $\alpha$ and
$\beta$ 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 $\Gamma_{\mathrm{VL}}$ and
the reciprocal-space area of the ellipse is given by $\pi
q_{\vartriangle}^{2}$, with
$q_{\vartriangle}=2\pi\left(2B/\sqrt{3}\Phi_{0}\right)^{1/2}$.

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 $\lambda$ and the coherence length $\xi$.
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 $\Omega$ between the
applied magnetic field and the $c$ axis. Full details of the experimental
protocols are given in the Methods. Fig. 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 $\Omega$, running from 0∘ to 87.5∘, prepared
using the field cooling process described in the Methods. It is clear that
rotating the field from being parallel to c axis ($\Omega=0^{\circ}$) to
nearly being in the basal plane ($\Omega=87.5^{\circ}$) distorts the hexagonal
vortex lattice, indicating the anisotropic nature of 2$H$-NbS2.

## Results and Discussion

### Vortex lattice structure

When the field ($B$ = 0.15 T) is parallel to c axis ($\Omega=0^{\circ}$), the
diffraction pattern is perfectly hexagonal (Fig. 2a), with the Bragg peaks
appearing at $q_{\vartriangle}=5.6(1)\times 10^{-3}$ Å-1, as compared with the
ideal value for a hexagonal VL
$q_{\vartriangle}=2\pi\left(2B/\sqrt{3}\Phi_{0}\right)^{1/2}=5.7\times
10^{-3}$Å-1 at $B$ = 0.15 T. $\Phi_{0}=h/2e$ is the flux quantum. Because of
the six-fold crystal symmetry, one of two hexagonal domains separated by 30∘
might be expected to be energetically favoured. We observe that one domain
(Hexagonal Domain I [  ]) is dominant, with trace amounts of Hexagonal Domain
II [  ]. Prior to changing $\Omega$, the diffraction spots of Domain I lie
parallel to the $\textbf{a}/\textbf{b}$ axes, indicating that, in real space,
the vortex lattice planes are perpendicular to the Nb nearest neighbour
directions. When $\Omega=0^{\circ}$, the $q_{x}$ and $q_{y}$ directions on the
detector are probing directions at right angles within the basal plane, and we
expect the superconducting parameters ($\lambda$ and $\xi$) to be essentially
isotropic around this plane.

At $\Omega=30^{\circ}$ 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 centre 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 (Fig. 2b,d).

As the field is rotated further, Domain II becomes the dominant form, and
indeed is the only domain seen for $\Omega>70^{\circ}$. 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 $\Omega$ 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 Biswas  _et al._
(2012).

As $\Omega$ is increased, $q_{y}$ remains in the same orientation in the basal
plane of the crystal lattice, but $q_{x}$ starts to mix basal plane and out-
of-plane components. A VL structure with Bragg peaks pinned to the $q_{y}$
direction should then be preferred, as observed in, for example, YBa2Cu3O7-δ
Kealey  _et al._ (2001), Sr2RuO4 Rastovski  _et al._ (2013) and KFe2As2 Kuhn
_et al._ (2016). However, as noted above, this was _not_ observed in
2$H$-NbS2’s sister compound, 2$H$-NbSe2, presumably due to a stabilising
effect associated with the CDW formation Gammel  _et al._ (1994).

Figure 3: The B-$\Omega$ phase diagram of the VL in 2$H$-NbS2 at $T$ = 1.5 K.
The shaded areas are represented by a gradient of green color, transitioning
from dark green to pale green, signifying the change in VL structure as the
direction of the magnetic field B is changed relative to the crystallographic
$ab$-plane. Hexagon symbols represent distinct vortex lattice domains, with
dark blue hexagons for domain I (  ), light green for domain II (  ), and a
teal dodecagon indicating regions with mixed domains (I $\&$ II). The dashed
black lines define the transitions between these domains as a function of $B$
and $\Omega$. The upper critical field ($B_{\text{c2}}$) line is the
parameterisation from the ratio of $B_{c2}^{\text{orb
}}(\theta)/B_{c2\|ab}^{\text{orb }}$ using values from Ref. woo Cho  _et al._
(2022). $B_{c1\|c}$ is estimated to be $\sim$ 30 mT Majumdar  _et al._ (2020).
Figure 4: Spin-splitting of the VL Bragg reflections due to spin-flip
scattering. Vortex lattice diffraction patterns as a function of rotation
angle ($\Omega$) 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
$\Omega$ 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 ($\phi+\Delta\phi)$ relative to the rocking curve
center which has slight zero offset at $\phi=0.09^{\circ}$. 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
($\Omega$). The scattering geometry triangles are depicted in (e) (with angles
exaggerated for clarity) for the two different SF processes: Spin-up to spin-
down (L.H.S) and spin-down to spin-up (R.H.S), The scattering angle (
$2\theta=2\phi_{0}$) is the same in both cases, but different tilt angles are
required to satisfy the Bragg condition $\phi_{2/3}$ (kinetic energy loss) and
$\phi_{1/4}$ (kinetic energy gain) for the top and bottom Bragg spots,
respectively.

### Superconducting anisotropy

The vortex lattice becomes increasingly anisotropic as the field is rotated
toward the basal plane, because the London penetration depth is different in
the basal plane and along the c axis. This leads to the distortion of the
diffraction pattern because the screening currents no longer follow circular
paths close to the vortex core. 2$H$-NbS2 falls into the class of uniaxial
anisotropic superconductors, which have a distortion of the hexagonal VL as a
function of $\Omega$. The anisotropy of the vortex lattice,
$\Gamma_{\mathrm{VL}}$ is defined as the ratio of the semi-major and semi-
minor axes of the ellipse circumscribing the diffraction spots, as given by
Campbell  _et al._ (1988)

$\Gamma_{\mathrm{VL}}=\frac{\Gamma_{ac}}{\sqrt{\sin^{2}\Omega+(\Gamma_{ac}\;\cos\Omega)^{2}}},$
(1)

where $\Gamma_{ac}$ corresponds to the $a-c$ anisotropy of the penetration
depth. To evaluate $\Gamma_{ac}$, we therefore need to measure
$\Gamma_{\mathrm{VL}}$ as a function of $\Omega$. There are multiple ways to
extract this information, and different ways have to be used at different
angles, primarily because not all of the spots are measured at higher angles
(Fig. 2).

For $\Omega\leq 70^{\circ}$, the following methods were used:

  * •

The area of the Brillouin zone associated with the VL is fixed for a given
value of field, and can be calculated directly for the perfect hexagon. Given
this, the position in reciprocal space of the spots, or, equivalently, the
opening angles $\alpha$ and $\beta$ (marked in Fig. 2a-2c) can be used to
calculate $\Gamma_{\mathrm{VL}}$, although the two domains require slightly
different treatment. For Domain I (  ),
$\Gamma_{\mathrm{VL}}=\sqrt{3}~{}[\tan(\alpha/2)]$, or by taking the $q_{y}$
value for the spots  , $\Gamma_{\mathrm{VL}}=3/4[q_{\vartriangle}/q_{y}]^{2}$.
For Domain II (  ), we have
$\Gamma_{\mathrm{VL}}=[q_{\vartriangle}/q_{y}]^{2}$ where $q_{y}$ is measured
for the spots  , or $\Gamma_{\mathrm{VL}}=1/(\sqrt{3}~{}[\tan(\beta/2)])$.

  * •

Fitting an ellipse that meets the area constraint ($A_{ellipse}=\pi
q_{\vartriangle}^{2}=8\pi^{3}B/\sqrt{3\Phi_{0}}$, shown in Fig. 2d) to the six
Bragg spot positions gives the semi-major and semi-minor axes, and hence
$\Gamma_{\mathrm{VL}}$ directly. This method cannot be used for
$\Omega>70^{\circ}$ as only the top and bottom spots are visible.

For $\Omega>70^{\circ}$, only the top and bottom Bragg spots of Domain II are
visible, so the first of these methods is employed; in this case the position
of the spots is $q_{\mathrm{VL}}=q_{y}$. However, in anisotropic
superconductors when the vortices are tilted away from a principal axis, the
field distribution associated with the vortices develop transverse field
components that vary with $\Omega$ Thiemann  _et al._ (1989). (We later give
expressions for these in Eqs. 6, 7, 8) Even if the near-horizontal spots at
large angles did fall on the detector, they would have almost no intensity.
This because the transverse components (labelled $b_{x}$ & $b_{y}$ in Fig. 1c)
are close to zero at small $q_{y}$ and the longitudinal component $b_{z}$ also
falls off at large $\Omega$. For the vertical spots, the transverse field
component $b_{x}$ dominates and flips the spin of the scattering neutrons.
Neutrons after spin flipping parallel or anti-parallel to the applied magnetic
field will have slightly different Zeeman energies, leading to a change in the
kinetic energy and hence velocity of the scattered neutron. This results in
the Bragg spots splitting into two separate peaks in the $\phi$ rocking scans
with maxima at $\phi=\phi_{B}\pm\Delta\phi$, where $\phi_{B}$ is the Bragg
angle expected for elastic scattering from the distorted vortex lattice
Rastovski  _et al._ (2013); Kuhn  _et al._ (2016); Kealey  _et al._ (2001),
and $\Delta\phi$ is the magnitude of the spin splitting of $\phi$. Examples of
the data are given in Fig. 4, together with schematic illustrations in Fig. 4e
of the effect of the energy change on the scattering process.

This spin-splitting of the peaks, 2$\Delta\phi$ becomes more pronounced as
$\Omega$ increases, because it is a function of $\Gamma_{\mathrm{VL}}$. This
arises because
$2\Delta\phi=\left(2k_{0}/q_{\mathrm{VL}}\right)\left(\Delta\varepsilon/\varepsilon_{0}\right)$
where $q_{\mathrm{VL}}$ is the magnitude of the scattering vector along the
minor axis i.e. $q_{\mathrm{VL}}=q_{\mathrm{y}}$ as depicted in Fig. 4b,
$\Delta\varepsilon=\gamma\mu_{N}B$ and
$\varepsilon_{0}=\hbar^{2}k_{0}^{2}/2m_{n}$, the neutron gyromagnetic ratio
$\gamma=1.913$, the nuclear magneton $\mu_{N}=e\hbar/2m_{n}$ and $m_{n}$ is
the neutron mass. This gives $\Delta\phi=C\,\Gamma_{\mathrm{VL}}\phi_{B}$
where $C=\gamma\sqrt{3}/4\pi=0.2635$. We therefore have four peaks, two
associated with the spot seen in the upper half of the detector in Fig. 4,
with centers $\phi_{1/2}$ and two associated with the lower spot,
$\phi_{3/4}$. Their centres are at angles:

$\displaystyle\phi_{1/2}$ $\displaystyle=\phi_{B}\pm C\ \Gamma_{\mathrm{VL}}\
\phi_{B}$ (2) $\displaystyle\phi_{3/4}$ $\displaystyle=-\phi_{B}\mp C\
\Gamma_{\mathrm{VL}}\ \phi_{B}$ (3)

$\Gamma_{\mathrm{VL}}$ can then be extracted from the rocking curves by a
simultaneous fit of the four peaks. The individual peaks are treated as
Lorentzians with a common integrated intensity and width under each field and
$\Omega$ condition. The relations between the peak centres are fixed by the
equations above, giving two outputs: $\Gamma_{\mathrm{VL}}$ and the zero error
in the $\phi$ motor positioning (found to be 0.09∘).

All of these methods have been used to evaluate $\Gamma_{\mathrm{VL}}$ where
possible, and they all agree within experimental error. The values obtained at
each $\Omega$ are field independent. The weighted average of all methods at
each condition is shown in Fig. 5. By fitting to Eqn. 1, the superconducting
anisotropy is found to be $\Gamma_{ac}=7.07\pm 0.2$.

For anisotropic superconductors, like 2$H$-NbS2, the variations in the
penetration depth ultimately arise due to differences in the Fermi velocities
and effective masses of the carriers in the different directions.
$\Gamma_{ac}$ can therefore also be observed in the ratios of the
superconducting coherence length, $\xi_{\mathrm{ab}}/\xi_{\mathrm{c}}$, and
the upper critical fields, $B_{c2}^{ab}/B_{c2}^{c}$. The upper critical fields
have been measured by transport at 2 K, giving anisotropies ranging from 6.5
to 8.1 Onabe  _et al._ (1978); Huang  _et al._ (2022), and a temperature
dependent study of the heat capacity by Kačmarčík _et al._ Kačmarčík  _et al._
(2010) finds a value of 7 for temperatures $0.3~{}T_{c}<T<T_{c}$. 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 Leroux  _et al._ (2012), while
Cho _et al._ woo Cho  _et al._ (2022) 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_{c2}^{ab}/B_{c2}^{c}$ = 7.33. However, this upper critical field
appears to be significantly lower than the orbital upper critical field
estimated from measurements close to $T_{c}$, which would give $\Gamma_{ac}$
$\sim$ 15\. Bi _et al._ Bi  _et al._ (2022) have pointed out that surface
superconductivity (up to $B_{c3}=1.695B_{c2}$) may be playing a role here.

In multi-band superconductors, the different bands may have different
anisotropies, leading to field- and temperature-dependent superconducting
anisotropies; a classic example is MgB2 Cubitt  _et al._ (2003). 2$H$-NbS2 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._ Kačmarčík  _et al._ (2010). To
engage with this question further, we first need to extract more information
about the characteristic superconducting lengths from our data.

Figure 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 ($\Omega$). The dashed line shows the VL anisotropy
calculated using Eqn. 1, with $\Gamma_{ac}=7.07$; the solid line is for
$\Gamma_{ac}=\infty$.

### 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(\textbf{q}_{i})$ of a Bragg reflection from VL domain $i$ via the Christen
formula Christen  _et al._ (1977):

$I(\textbf{q}_{i})=2\pi
V_{i}S\big{(}\frac{\gamma}{4}\big{)}^{2}\frac{\lambda_{n}^{2}}{\Phi_{0}^{2}q_{i}\cos(\zeta)}|F(\textbf{q}_{i})|^{2}.$
(4)

Here, $V_{i}$ is the volume of the sample occupied by the VL domain $i$ (with
$V_{\textrm{I}}+V_{\textrm{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, $\textbf{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 $\textbf{q}_{i}$ and the direction that is at right angles to the
rocking axis Squires (2012). For each domain, the relevant integrated
intensities of individual Bragg peaks (see Methods) are averaged to give
$I(\textbf{q}_{i})$ for that domain; this is related to the form factor via
Eqn. 4.

Figure 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^{\circ}$ to $87.5^{\circ}$. 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(\textbf{q}_{i})|^{2}$ with their
respective error bars (see Supplementary Note 1). The dashed black lines are
the results from fits to Eqn. 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_{\text{VL}}$($\Omega$).

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. 1c). In 2$H$-NbS2, 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 Kogan (1981) for uniaxial crystal systems like
2$H$-NbS2, and expanded to account for the effective mass anisotropy of the
carriers by Thiemann _et al._ Thiemann  _et al._ (1989). Kealey _et al._
applied this practically to the biaxial superconductor YBa2Cu3O7-δ Kealey  _et
al._ (2001), correcting some misprints in Ref. Thiemann  _et al._ (1989), but
were only able to find qualitative agreement, as was the case in studies on
Sr2RuO4 Rastovski  _et al._ (2013); Kuhn  _et al._ (2017) and KFe2As2 Kuhn
_et al._ (2016). This may be due to difficulties in handling field-dependent
superconducting anisotropies. In 2$H$-NbS2, 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 $\mathbf{B}(\mathbf{r})$
can be decomposed into a Fourier series over the set of reciprocal space wave
vectors $\mathbf{q}$:

$\mathbf{B}(\mathbf{r})=\sum_{\mathbf{q}}\mathbf{b}(\mathbf{q})\exp(i\mathbf{q}\cdot\mathbf{r})$
(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 $\bf{q}$ = $(q_{x},q_{y})$, the field components are
Thiemann  _et al._ (1989):

$\displaystyle b_{x}$ $\displaystyle=\frac{(\lambda^{2}m_{xz}q_{y}^{2})B}{d}$
(6) $\displaystyle b_{y}$
$\displaystyle=\frac{(-\lambda^{2}m_{xz}q_{x}q_{y})B}{d}$ (7) $\displaystyle
b_{z}$ $\displaystyle=\frac{(1+\lambda^{2}m_{zz}q^{2})B}{d}$ (8)

where

$\displaystyle d=$
$\displaystyle\left(1+\lambda^{2}m_{yy}q_{x}^{2}+\lambda^{2}m_{xx}q_{y}^{2}\right)\left(1+\lambda^{2}m_{zz}q^{2}\right)-\lambda^{4}m_{xz}^{2}q^{2}q_{y}^{2}.$
(9)

Here $B$ is the average field and
$\lambda=\left(\lambda_{ab}^{2}\lambda_{c}\right)^{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. 1c. The normalized
effective masses along the unit cell axes are $m_{a}$, $m_{b}$, and $m_{c}$.
For 2$H$-NbS2, $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$:

$\displaystyle m_{xx}$
$\displaystyle=\Gamma_{ac}^{-2/3}\cos^{2}\Omega+\Gamma_{ac}^{4/3}\sin^{2}\Omega,$
(10) $\displaystyle m_{yy}$ $\displaystyle=\Gamma_{ac}^{-2/3},$ (11)
$\displaystyle m_{zz}$
$\displaystyle=\Gamma_{ac}^{-2/3}\sin^{2}\Omega+\Gamma_{ac}^{4/3}\cos^{2}\Omega,$
(12) $\displaystyle m_{xz}$
$\displaystyle=\left(\Gamma_{ac}^{-2/3}-\Gamma_{ac}^{4/3}\right)\sin\Omega~{}\cos\Omega.$
(13)

With these expressions we can then calculate the theoretical form factor,
$F(\mathbf{q})$, for a given Bragg reflection, including a Gaussian cutoff
term to account for the finite size of the vortex core Yaouanc  _et al._
(1997):

$\displaystyle
F(\mathbf{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})].$
(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∘ 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
Campillo  _et al._ (2022).

$\boldsymbol{\Omega}$ (∘) |  $\lambda_{\textbf{GM}}$ (nm) |  $\lambda_{\textbf{ab}}$ (nm) |  $\lambda_{\textbf{c}}$ (nm) |  $\xi_{\textbf{ab,eff}}$ (nm)
---|---|---|---|---
30 |  277 $\pm$ 5 |  145 $\pm$ 3 |  1022 $\pm$ 28 |  19.5 $\pm$ 0.5
60 |  275 $\pm$ 3 |  144 $\pm$ 2 |  1014 $\pm$ 22 |  20.1 $\pm$ 0.2
70 |  270 $\pm$ 2 |  141 $\pm$ 1 |  995 $\pm$ 21 |  23.2 $\pm$ 0.4
80 |  273 $\pm$ 7 |  142 $\pm$ 4 |  1006 $\pm$ 32 |  34 $\pm$ 1
85 |  272 $\pm$ 10 |  142 $\pm$ 5 |  1001 $\pm$ 40 |  34 $\pm$ 1.5
87.5 |  273 $\pm$ 3 |  142 $\pm$ 1 |  1006 $\pm$ 21 |  24.6 $\pm$ 0.4
Table 2: 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 Eqn. 14.

The theoretical form factors depend on sample properties $\lambda$, $\xi_{ab}$
and $\xi_{c}$. We fitted the field-dependent data separately 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^{2}_{c}\sin^{2}\Omega<<\xi^{2}_{ab}\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_{\textrm{I}}/V$ and $V_{\textrm{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 A). Each value represents the average taken over
equivalent spots for a given domain. The $(V_{i}/V)|F(\mathbf{q}_{i})|^{2}$
were calculated from the experimental data using Eqn. 4 and then fitted to
Eqn. 14 using two different methods.

In the first approach, the intensities were fitted using the theoretical
values for the scattering vector $\mathbf{q}$ associated with the peak given
that $\Gamma_{ac}=7$. For Domain I (  ),
$q_{x}=(q_{\vartriangle}/2)({\Gamma_{\mathrm{VL}})^{1/2}}$ and
$q_{y}=(q_{\vartriangle}/2)({3/\Gamma_{\mathrm{VL}}})^{1/2}$, while for Domain
II (  ) $q_{x}=0$ and
$q_{y}=(q_{\vartriangle})/({\Gamma_{\mathrm{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 B). 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 Eqn.
14 using the theoretical Bragg spot positions.

The fits show excellent agreement with the description of the anisotropic
superconductor given by Thiemann _et al._ Thiemann  _et al._ (1989),
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.

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 NbS2, and for our
sample, it is 1.8 T Majumdar  _et al._ (2020). 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 Tinkham (1996), 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._ woo Cho  _et al._ (2022) 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 Leroux  _et al._ (2012) and 590
nm Tulapurkar  _et al._ (2002) are reported. On our samples, a value for
$\lambda_{ab}$ of 131 nm was obtained by Majumdar _et al._ Majumdar  _et al._
(2020).

Recently, Kogan _et al._ Kogan  _et al._ (2020) established a relationship
between the zero-temperature penetration depth $\lambda(0)$ and the slope of
the penetration depth $\lambda^{-2}(T)$ near $T_{\rm c}$, using a calculation
analogous to the Helfand-Werthamer relationship between the zero-temperature
upper critical magnetic field and its slope at $T_{\rm 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},$
(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_{\rm c}$, $T_{\rm 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 non-magnetic
scattering, and was successfully cross-checked against experimental values for
$\lambda$ in V3Si and Nb3Sn.

Using this model we have calculated $\lambda_{ab}$(0) $\approx 142.9$ nm for
2$H$-NbS2, taking as inputs $T_{\rm c}=5.5~{}K$,
$(dH^{c}_{c2}/dT)_{T=T_{c}}\approx-0.25\times 10^{4}\ \text{Oe/K}$ Kačmarčík
_et al._ (2010); Huang  _et al._ (2022); woo Cho  _et al._ (2022), and
$\gamma=0.494\times 10^{4}\ \text{erg/cm}^{3}\text{K}^{2}$ Hamaue and Aoki
(1986). 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 2$H$-NbS2 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 MgB2, 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 2$H$-NbS2, which has predominantly
2D character, as proposed by Kačmarčík _et al._ Kačmarčík  _et al._ (2010).
Interestingly, Noat _et al._ Noat  _et al._ (2015) have suggested that, in
reality, only the large gap sheets (the cylinders centred on the K points) are
intrinsically superconducting, and that the other gap develops parasitically
by coupling through the pancake-like sheet coming from the chalcogen
$p$-bands.

## Conclusions

We have measured neutron diffraction by the vortex lattice in 2$H$-NbS2 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 Thiemann  _et al._ (1989), 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._ Kogan  _et al._ (2020).

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 2$H$-NbS2 has been found by Noat _et al._ Noat  _et al._ (2015)
to be half that in 2$H$-NbSe2, from fitting tunnelling spectroscopy data. It
is therefore not surprising that 2$H$-NbS2 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._ Gammel  _et al._ (1994) for
2$H$-NbSe2 using the same method presented here.

## Methods

Thin platelets of 2$H$-NbS2 samples were grown using the chemical vapour
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._ Chareev  _et al._ (2020).
The sample quality was confirmed by specific heat measurements at zero field,
wherein a sharp jump centered at the superconducting transition temperature
$T_{\rm c}$ = 5.5 K is seen. This is considered to be the most reliable way to
check sample quality Witteveen  _et al._ (2021), and matches well with the
behaviour reported for other samples grown in the same laboratory Cho  _et
al._ (2021).

The 2$H$-NbS2 polytype is hexagonal, with space group $P6_{3}/mmc$ and lattice
parameters $a=b\sim 3.31$ Å and $c\sim 11.86$ Å Witteveen  _et al._ (2021).
The platelets grow with the c axis normal to the platelet surface. A mosaic of
nine co-aligned crystals was made by mounting the platelets onto two aluminum
sheets, so that all the crystal faces lay within an area $10\times
10\,\text{mm}^{2}$. The mosaic had a total crystal thickness of 80 µm with a
total mass of $42.3\,\text{mg}$. The in-plane alignment (along the
$\textbf{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 $\textbf{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 Dewhurst  _et al._ (2016)
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$\times$128 pixels) two-
dimensional position-sensitive multi-detector placed 10.0435 m after the
sample, which was mounted in a 9 T horizontal-field cryomagnet. The crystals
were aligned with the $\textbf{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^{\circ}$, 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_{\rm 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 Marziali
Bermúdez  _et al._ (2015). This was done for magnetic fields over the range
0.15 T to 1.25 T.

To probe the superconducting anisotropy of 2$H$-NbS2, we also rotated the
magnetic field away from the c axis towards the b direction by an angle
$\Omega$. These rotations were all carried out above $T_{\rm 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^{\circ}$, 30∘, 80∘ and 85∘, a 12 mm diameter circular aperture was
used. For $\Omega=60^{\circ}$ and 70∘, a 7H $\times$ 10V mm2 rectangular
aperture was used. For $\Omega=87.5^{\circ}$, the aperture size was 3H
$\times$ 10V mm2.

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
Dewhurst (2023). 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  , are used, whereas for Domain II, only the top and
bottom Bragg spots,  , are included. For the form factor calculations, the
correct region of the detector to include for each peak was tested iteratively
to maximise 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(\textbf{q}_{i})$. For $\Omega\leq
70^{\circ}$, 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.

## Acknowledgements

AA and EB acknowledge support from the Crafoord Foundation (no. 20190930) and
the Swedish Research Council under Project No. 201804704. AA, EMF and EB thank
the Institut Laue-Langevin for travel support to attend the experiment. The
authors gratefully acknowledge the Institut Laue-Langevin for the allocated
beamtime.

## Competing interests

The authors have no conflicts to disclose.

## Materials and correspondence

The experimental data collected at the Institut Laue-Langevin is publicly
available Alshemi  _et al._ (2021). Processed experimental data are available
upon request to the corresponding authors.

## APPENDIX A: Longitudinal and Transverse magnetic field components of the
flux line lattice in the anisotropic 2$H$-NbS2

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. Firstly, the penetration depth is changing as $\Omega$
changes. Hence, the form factor of the observed Bragg spots on the left and
right in the SANS patterns for both Domain I and Domain II (diagonal spots)
decreases 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_{\text{c2}}$, as the field is rotated away from the c
axis, the upper critical field $B_{\text{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_{\text{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^{\circ}$, 85∘, and $87.5^{\circ}$ 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 Eqns. 6-8.

Figure 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
T and 1.25 T.

## APPENDIX B: 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^{\circ}$, 85∘ and $87.5^{\circ}$ as an example. The red data points are
the results from the fitting process using Approach II for each angle
$\Omega$. The black 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_{GM}$, $\xi_{ab}$, and
$\xi_{c}$ within the experimental error.

Figure 8:  Vortex lattice form factor for angles $\Omega=80^{\circ}$, 85∘ and
$87.5^{\circ}$ as a function of the field at 1.5 K.

## References

  * Cho  _et al._ (2021) 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, Nature Communications 12, 3676 (2021).
  * Wan  _et al._ (2023) 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, Nature 619, 46 (2023).
  * Fulde and Ferrell (1964) P. Fulde and R. A. Ferrell, Physical Review 135, A550 (1964).
  * Larkin and Ovchinnikov (1965) A. Larkin and Y. N. Ovchinnikov, Soviet Physics-JETP 20, 762 (1965).
  * Manzeli  _et al._ (2017) S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nature Reviews Materials 2, 17033 (2017).
  * Knispel  _et al._ (2024) T. Knispel, J. Berges, A. Schobert, E. G. C. P. van Loon, W. Jolie, T. Wehling, T. Michely, and J. Fischer, Nano Letters 24, 1045 (2024), pMID: 38232959, https://doi.org/10.1021/acs.nanolett.3c02787 .
  * Lin  _et al._ (2020) D. Lin, S. Li, J. Wen, H. Berger, L. Forró, H. Zhou, S. Jia, T. Taniguchi, K. Watanabe, X. Xi, and M. S. Bahramy, Nature Communications 11, 2406 (2020).
  * Moncton  _et al._ (1977) D. E. Moncton, J. D. Axe, and F. J. DiSalvo, Phys. Rev. B 16, 801 (1977).
  * Wilson  _et al._ (1975) J. A. Wilson, F. J. Di Salvo, and S. Mahajan, Advances in Physics 24, 117 (1975).
  * Abdel-Hafiez  _et al._ (2016) 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, Sci. Rep. 6, 31824 (2016).
  * Guillamon  _et al._ (2008) I. Guillamon, H. Suderow, S. Vieira, L. Cario, P. Diener, and P. Rodière, Phys. Rev. Lett. 101, 166407 (2008).
  * Kačmarčík  _et al._ (2010) J. Kačmarčík, Z. Pribulová, C. Marcenat, T. Klein, P. Rodière, L. Cario, and P. Samuely, Phys. Rev. B 82, 014518 (2010).
  * Kobayashi  _et al._ (1977) N. Kobayashi, K. Noto, and Y. Muto, Journal of Low Temperature Physics 27, 217 (1977).
  * Majumdar  _et al._ (2020) 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, Phys. Rev. Mater. 4, 084005 (2020).
  * Noat  _et al._ (2015) 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, Phys. Rev. B 92, 134510 (2015).
  * Ichioka  _et al._ (2017) M. Ichioka, V. G. Kogan, and J. Schmalian, Phys. Rev. B 95, 064512 (2017).
  * Hess  _et al._ (1994) H. F. Hess, C. A. Murray, and J. V. Waszczak, Phys. Rev. B 50, 16528 (1994).
  * Gammel  _et al._ (1994) 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, Phys. Rev. Lett. 72, 278 (1994).
  * Kealey  _et al._ (2001) 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, Phys. Rev. B 64, 174501 (2001).
  * Rastovski  _et al._ (2013) C. Rastovski, C. D. Dewhurst, W. J. Gannon, D. C. Peets, H. Takatsu, Y. Maeno, M. Ichioka, K. Machida, and M. R. Eskildsen, Phys. Rev. Lett. 111, 087003 (2013).
  * Kuhn  _et al._ (2016) 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, Phys. Rev. B 93, 104527 (2016).
  * Biswas  _et al._ (2012) 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, Phys. Rev. Lett. 108, 077001 (2012).
  * woo Cho  _et al._ (2022) C. woo Cho, C. Y. Ng, C. H. Wong, M. Abdel-Hafiez, A. N. Vasiliev, D. A. Chareev, A. G. Lebed, and R. Lortz, New Journal of Physics 24, 083001 (2022).
  * Campbell  _et al._ (1988) L. J. Campbell, M. M. Doria, and V. G. Kogan, Phys. Rev. B 38, 2439 (1988).
  * Thiemann  _et al._ (1989) S. L. Thiemann, Z. Radović, and V. G. Kogan, Phys. Rev. B 39, 11406 (1989).
  * Onabe  _et al._ (1978) K. Onabe, M. Naito, and S. Tanaka, J. Phys. Soc. Jpn. 45, 50 (1978).
  * Huang  _et al._ (2022) D. Huang, H. Nakamura, K. Küster, U. Wedig, N. B. M. Schröter, V. N. Strocov, U. Starke, and H. Takagi, Phys. Rev. B 105, 245145 (2022).
  * Leroux  _et al._ (2012) M. Leroux, P. Rodière, L. Cario, and T. Klein, Physica B 407, 1813 (2012).
  * Bi  _et al._ (2022) 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, Phys. Rev. Res. 4, 013188 (2022).
  * Cubitt  _et al._ (2003) R. Cubitt, M. R. Eskildsen, C. D. Dewhurst, J. Jun, S. M. Kazakov, and J. Karpinski, Phys. Rev. Lett. 91, 047002 (2003).
  * Christen  _et al._ (1977) D. K. Christen, F. Tasset, S. Spooner, and H. A. Mook, Phys. Rev. B 15, 4506 (1977).
  * Squires (2012) G. L. Squires,  _Introduction to the Theory of Thermal Neutron Scattering_, 3rd ed. (Cambridge University Press, 2012).
  * Kogan (1981) V. G. Kogan, Phys. Rev. B 24, 1572 (1981).
  * Kuhn  _et al._ (2017) 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, Phys. Rev. B 96, 174507 (2017).
  * Yaouanc  _et al._ (1997) A. Yaouanc, P. Dalmas de Réotier, and E. H. Brandt, Phys. Rev. B 55, 11107 (1997).
  * Campillo  _et al._ (2022) E. Campillo, M. Bartkowiak, R. Riyat, E. Jellyman, A. S. Cameron, A. T. Holmes, O. Prokhnenko, W.-D. Stein, A. Erb, E. M. Forgan, and E. Blackburn, Phys. Rev. B 105, 184508 (2022).
  * Tinkham (1996) M. Tinkham,  _Introduction to Superconductivity_ (McGraw-Hill, New York, 1996).
  * Tulapurkar  _et al._ (2002) A. Tulapurkar, A. Grover, S. Ramakrishnan, A. Niazi, and A. Rastogi, Physica B: Condensed Matter 312-313, 118 (2002), the International Conference on Strongly Correlated Electron Systems.
  * Kogan  _et al._ (2020) V. G. Kogan, M. A. Tanatar, and R. Prozorov, Phys. Rev. B 101, 104510 (2020).
  * Hamaue and Aoki (1986) Y. Hamaue and R. Aoki, Journal of the Physical Society of Japan 55, 1327 (1986).
  * Chareev  _et al._ (2020) D. A. Chareev, P. Evstigneeva, D. Phuyal, G. J. Man, H. Rensmo, A. N. Vasiliev, and M. Abdel-Hafiez, Crystal Growth & Design 20, 6930 (2020).
  * Witteveen  _et al._ (2021) C. Witteveen, K. Gornicka, J. Chang, M. Månsson, T. Klimczuk, and F. O. von Rohr, Dalton Transactions 50, 3216 (2021).
  * Dewhurst  _et al._ (2016) C. D. Dewhurst, I. Grillo, D. Honecker, M. Bonnaud, M. Jacques, C. Amrouni, A. Perillo-Marcone, G. Manzin, and R. Cubitt, Journal of Applied Crystallography 49, 1 (2016).
  * Marziali Bermúdez  _et al._ (2015) M. Marziali Bermúdez, M. R. Eskildsen, M. Bartkowiak, G. Nagy, V. Bekeris, and G. Pasquini, Phys. Rev. Lett. 115, 067001 (2015).
  * Dewhurst (2023) C. D. Dewhurst, Journal of Applied Crystallography 56, 1595 (2023).
  * Alshemi  _et al._ (2021) A. Alshemi, E. Blackburn, E. Campillo, R. Cubitt, E. M. Forgan, L. Shen, and A. Stellhorn, doi:10.5291/ILL-DATA.5-42-558 (2021).