Introduction

In metals or semimetals, application of a strong magnetic field (B) at low temperatures (T) leads to the quantization of electronic states into discrete Landau Levels with distinct energies, resulting in periodic oscillations of the density of states near the Fermi surface, dubbed quantum oscillation (QO)1. These oscillations were experimentally observed in magnetization and resistivity measurements in bismuth during the 1930s, known as de Haas–van Alphen (dHvA) oscillations and Shubnikov–de Haas (SdH) oscillations, respectively2,3. The oscillatory period is traditionally represented in an inverse magnetic field plot. The corresponding frequency (Fr) is proportional to the extremal cross-section of the Fermi surface (AF) perpendicular to the applied field, as dictated by the Onsager relation Fr = (φ0/2π2)AF, where φ0 = h/2e is the magnetic flux quantum. This characteristic establishes quantum oscillation as one of the most powerful tools for investigating the geometry of Fermi surface.

In recent years, nodal-line semimetals that feature Dirac band crossings as lines or closed loops in momentum space have attracted significant interest as new topological systems4,5,6,7,8,9. Among the experimentally verified nodal-line semimetals, ZrSiS stands out with its strong linear Dirac dispersion spanning an energy range of ~2 eV, indicating the presence of nearly massless transport fermions6,7. It is also notable for its intricate Fermi surface configuration and up to six fundamental frequencies: F(δ1) ~ 8.5 T, F(δ2) ~ 16.5 T, F(γ1) ~ 22 T, F(γ2) ~ 92 T, F(α) ~ 240 T, F(β)~419 T. Accordingly, clear QOs have been detected in various physical properties, including resistivity10,11,12, magnetization13, thermopower14, and magnetic torque15,16. Notably, the α and β pockets located on the Z-R-A plane can additionally exhibit the magnetic breakdown phenomena, leading to QO frequencies up to ~8000 T. However, no single method has been successful in fully capturing all these frequencies at temperatures above 2 K and in magnetic fields below 14 T. Moreover, even under high magnetic fields up to 31 T and temperatures as low as 0.4 K, certain frequencies like F(β) remained absent in SdH oscillations, which was proposed to be related to specific scattering mechanism16. Thus, discrepancies among different studies have only provided fragmented insights into the complete picture of the Fermi surface.

Typically, the magnitude of QOs at a specific frequency increases with declining temperatures and increasing magnetic fields. This can be ascribed to the thermal damping factor RT = αTμ/[Bsinh(αTμ/B)] and the field-related damping factor RD = exp(-αTDμ/B) in the Lifshitz-Kosevich (LK) formula, respectively1. Here, α = (2π2kBm0)/(ℏe) is a constant, μ = m*/m0 represents the relative effective mass, m0 denotes the mass of a free electron, and TD is the Dingle temperature. Consequently, to extract comprehensive information about the Fermi surface, investigations of QOs often require extreme conditions including ultralow temperatures (<1 K) and high magnetic fields (>20 T). However, achieving these conditions is challenging with standard cryogenic magnets, which typically operate at relatively moderate temperatures (>1.8 K) and within a limited field range (<14 T). Under these conditions, each measurement technique may only reveal a subset of QO frequencies, like in ZrSiS11,12,13,14. Currently, there is no single measurement capable of capturing all its quantum oscillation frequencies under general experimental conditions. Such a comprehensive measurement is highly desirable for exploring the complex and non-trivial Fermiology in semimetals. It would enable us to investigate complete Fermi surfaces in a single attempt, in case the sample orientation varies slightly when combining different methodologies. Therefore, this approach is particularly crucial for samples with Fermi surfaces that are both complex and sensitive to tilting angles.

Recently, we adopted a composite magnetoelectric (ME) technique as a highly sensitive tool for investigating magnetic materials17,18,19. A composite laminate magnetoelectric configuration is used, stemming from renewed interest in room temperature multiferroics from the year 200020. In this method, a thin sample is affixed to a piezoelectric transducer using silver epoxy to form a bilayer composite magnetoelectric laminate. The piezoelectric layer is a thin 0.7Pb(Mg1/3Nb2/3)O3-0.3PbTiO3 (PMN-PT) [001]-cut single crystal with a record high piezoelectric coefficient, which also became feasible around the year 200021,22. This configuration facilitates the conversion of the in-plane magnetostriction of the specimen into an electrical signal across the piezo-layer. The magnetostriction (λ) measures the response of sample geometry in external magnetic fields, i.e., λ(B) = ΔL(B)/L0, where L0 is the sample length in zero magnetic field, and ΔL(B) = L(B) – L0 is the length change induced by magnetic field. Unlike traditional dc measurements, the ac ME method employs the ac-lock-in technique, which achieves an enhanced signal-to-noise ratio and high resolution by applying an ac-driven magnetic field (Bac). This technique has demonstrated high accuracy in characterizing magnetic skyrmion phases in materials such as MnSi17 and Co7Zn8Mn518, and in unveiling ferromagnetic spin-spin correlations in quasi-2D magnets, even within paramagnetic phases19. Hence, the ac ME method has the potential for detecting weak QOs in the magnetostriction of semimetals under moderate field and temperature conditions.

Theoretically, the ac ME technique can achieve superior resolution in detecting QOs of magnetostriction in semimetals with high frequencies. The real part of the ac voltage signal Vx is expected to be proportional to the piezomagnetic coefficient μ0(∂λ)/∂B (where μ0 is the permeability of vacuum)23. This derivative (∂λ/∂B) of the oscillatory part of magnetostriction can be estimated to be:

$${V}_{x}\propto \frac{\partial \widetilde{{\rm{\lambda }}}}{\partial B} \sim -\sum _{r}\frac{2\pi {F}_{r}}{{B}^{2}}\widetilde{{\rm{\lambda }}}\left({\varphi }_{r}+\frac{\pi }{2}\right)$$
(1)

Where Fr is the oscillatory frequency and r is the index for it, φr is the corresponding oscillatory phase. This relation displays a phase shift of π/2 comparing to the original oscillatory magnetostriction, with an amplification factor of 2πFrB-2. The B-2 part of this factor can significantly enhance the low-field oscillatory signals. The 2πFr part of the factor provides a global enhancement to the oscillations, thus, high-frequency oscillations would be amplified further. On the other hand, the magnetostriction itself increases with decreasing Fr, i.e., with decreasing the size of the charge-carrier group24. Notably, for high frequencies above 100 T, the onset field for observing QOs may be significantly reduced since the amplification factor becomes even larger in lower magnetic fields.

In this study, we conducted an investigation on the nodal-line semimetal ZrSiS. We compared the ac composite ME method with four other techniques, including in-plane resistivity (ρ), out-of-plane magnetization (M), out-of-plane magnetostriction λ, and in-plane thermopower (S), using the same thin ZrSiS single crystal. The measurements were carried out at temperatures down to 2 K and in magnetic fields up to 13 T. Our findings demonstrate that the ac ME method can effectively reveal all six fundamental frequencies while other methods cannot. Besides, results of the ME method display a series of magnetic breakdown frequencies around 8000 T with a low onset magnetic field of 7.6 T. And the effective mass values of δ1 and α bands can be determined from temperature-dependent data, showing nearly zero values. These results highlight the potential of the ME technique for exploring Dirac and Weyl physics as well as Fermi surface properties in other metals, under moderate experimental conditions.

Results

Basic properties of ZrSiS

ZrSiS belongs to the WHM (W = Zr, Hf; H = Si, Ge, Sn, Sb; M = S, Se, Te) material system and crystallizes in the PbFCl structure type with a tetragonal P4/nmm space group (No. 129)25, as depicted in Fig. 1a. The Si square net is situated in the ab-plane and is enclosed between Zr-S layers. X-ray diffraction analysis of a single crystal (inset of Fig. 1b) confirms the high quality of the ZrSiS samples. The diffraction pattern in Fig. 1b exhibits well-defined sharp (00 L) peaks, and the calculated c-axis lattice constant is 8.060 Å using the Nelson-Riley method, consistent with previous reports26,27.

Fig. 1: Structure and basic sample characterizations of ZrSiS single crystal.
figure 1

a Schematic crystal structure of ZrSiS with P4/nmm space group. b Single crystal X-ray diffraction pattern measured at room temperature. Inset: A photograph of ZrSiS single crystal. c Temperature-dependent resistivity under external magnetic fields of 0 and 8 T, applied along the c-axis.

Full size image

Figure 1c illustrates the temperature dependence of resistivity along the a-axis for two magnetic fields: 0 T and 8 T, where the field is aligned parallel to the c-axis. In the absence of a magnetic field, the resistivity shows typical metallic behavior, gradually decreasing as the temperature drops. The high quality of the crystal is evident from the calculated residual resistivity ratio of approximately 50. Conversely, in 8 T, a clear transition occurs, leading to a pronounced insulating behavior below 100 K. This field-induced crossover can be linked to the emergence of a band gap in topological semimetals, a phenomenon consistent with prior research11.

Quantum oscillations in magnetic properties: magnetostriction and magnetization

As mentioned in various studies, the QOs are expected to be noticeable under high magnetic fields and low temperatures in different magnetic properties10,11,12,13,14,15,16. QOs in magnetostriction λ is a fundamental thermodynamic phenomenon that can directly reflect strain dependence of the Fermi surface. As depicted in Fig. 2a, field-dependent magnetostriction λ(B) for both B and λ along the c-axis, exceeding 10-5 at 13 T, is measured with a capacitive dilatometer at 1.6 K. An oscillatory component with a minor amplitude (10-7-10-8) and low frequency is shown in the inset of Fig. 2a by subtracting a dominating linear background. The onset field of the oscillations is about 1.5 T. Other higher frequency oscillations reported in previous studies are not found in λ(B) due to the thin crystal thickness (0.36 mm) and large background signal. In contrast, field-dependent magnetization M(B) (Fig. 2b) exhibits both low- and high-frequency oscillations. A similar dominant linear background is also identified while the onset field of the oscillations is reduced to 1 T, indicating better sensitivity in terms of QO. As a result, these two non-transport methods yield considerably less information within limited field and temperature condition compared to the known exotic features found in QO of ZrSiS appearing in high magnetic fields10,15,16.

Fig. 2: Field dependent physical properties at ~2 K with B // c-axis.
figure 2

a Field-dependent magnetostriction along the c-axis (L0 = 0.36 mm). Inset of (a) displays the oscillatory part with the background subtracted. b Field-dependent magnetization. d Field-dependent resistivity. e Field-dependent thermopower. Inset of (e) shows ultra-high-frequency oscillations in the field range above 12 T. c Real part of the ac voltage signal induced by in-plane deformation of the sample, which is adhered to PMN-PT. The deformation is stimulated by a sinusoidal magnetic field Bac (left panel of inset of (c)). Right panel of inset of (c) shows ultra-high frequency oscillations in the field range above 12 T.

Full size image

Quantum oscillations of ZrSiS in in-plane resistivity and thermopower

To show the QOs detected by transport properties, we measured the field-dependent resistivity ρ and thermopower S on the same crystal used in magnetic property measurements, as depicted in Figs. 2d, e, respectively. Field-dependent resistivity ρ(B) (Fig. 2d) presents strong high-frequency and weak low-frequency oscillations. The dominant background of the ρ curve follows a B2-like profile, which is common in metals. The substantial background raises its discernable onset field of oscillations to over 4 T. Conversely, field-dependent thermopower S(B) at 1.6 K reveals more complex oscillations (Fig. 2e). Its background is much smaller than that of other methods. Consequently, its oscillation part shows an onset field of the oscillations as low as 1.5 T. The 8000 T QO, which has a relatively weak magnitude, appears around 8.7 T. The overall field-dependent profile of S(B) is consistent with previous studies14.

Quantum oscillations of ZrSiS with composite magnetoelectric configuration

We found that the composite ME method exhibits a superior sensitivity compared to the above four methods. To employ the composite ME method on ZrSiS, we mechanically bonded the single crystal with a thin PMN-PT to produce a ZrSiS/PMN-PT laminate, as depicted in the inset of Fig. 2c. The real part of ME signal Vx as a function of the dc magnetic field B is measured at 2 K, as shown in Fig. 2c. QOs with various frequencies are displayed, revealing more details than the above four methods. Vx has a weak linear background, leading to a dominant oscillation signal as the magnetic field increases. The onset field of the oscillations drops to 0.5 T, which suggests that the oscillatory component of magnetostriction far exceeds the noise level. Additionally, Vx exhibits oscillations around 8000 T, as shown in the inset of Fig. 2c, similar to that in thermopower (inset of Fig. 2e). However, a lower starting field of 7.6 T is observed in the ME technique. Based on these features, the ME technique appears to be an excellent method for characterizing QOs.

Feature comparison and FFT analyses of QOs in ZrSiS

In this study, QOs measured by different methods are quantitatively analyzed. Fast Fourier transformation (FFT) analyses are employed to extract the oscillatory frequencies and amplitudes for each method, which are presented in Fig. 3k-o. The oscillations of various physical properties are plotted as a function of 1/B (inverse magnetic field) with the background subtracted. The entire field region is divided into two sections: high field region (0.07 < 1/B < 0.2 T-1, Fig. 3a–e) and low field region (0.2 < 1/B < 0.5 T–1, Fig. 3f–j). For 1/B > 0.5 T–1, the oscillation is too weak to be analyzed. The six fundamental frequencies are recognized by FFT in different field ranges in our ac ME method (Supplementary Note 1 and Supplementary Figure 1), and marked as dashed lines in the FFT plots. The discrepancies between different methods are revealed in the comparison. Oscillation frequencies of F(δ1) ~ 8 T and F(δ2) ~ 16 T are common among all the methods, except in resistivity where these two frequencies are not resolvable, as shown in Fig. 3n. F(γ1) ~ 22 T appears to exist only in the oscillations of magnetostriction and ac ME signals (Vx), but the FFT peak in Vx is more pronounced than that in λ, with its onset field (~2.85 T) indicated in the low field region (Fig. 3g). F(γ2) ~ 91 T is difficult to discern in any oscillation data directly, while only the FFT result of Vx (Fig. 3i) shows a discrete peak at this frequency. For the fundamental frequencies above 100 T, F(α) ~240 T oscillations are common in these methods except in λ (Fig. 3l–o). F(β) ~420 T only appears in Vx, with the smallest amplitude in the FFT profile (Fig. 3l) compared to all the other fundamental frequencies.

Fig. 3: 1/B plots and FFT analyses of the oscillations of five physical properties.
figure 3

a-e 1/B plots ranging from 0.07 T–1 to 0.25 T–1 for five physical properties. Scale bar length in (b) is 1/22.5 T–1, which implies a single period of F(γ1) oscillations. f–j 1/B plots for the same properties ranging from 0.25 T–1 to 0.8 T–1. With a scale bar length of 1/8.5 T–1 in h, one single period of F(δ1) oscillations can be revealed. k–o FFT analyses for each of the five physical properties. The analyses highlight six fundamental oscillatory frequencies in ZrSiS, indicated by six dotted lines.

Full size image

It is noted that, in the high-field region (Fig. 3b) of Vx, the oscillations exhibit a complex pattern, indicating that more frequencies are included. The FFT analysis suggests a series of oscillations in the frequency range of 60 to 600 T (65, 183, and 593 T), which are combinations between F(α) and F(β) orbits due to the magnetic breakdown phenomenon16. Only the 593 T frequency can be found in thermopower data as well (Fig. 3o). All these fundamental and breakdown frequencies below 1 kT are summarized in Table 1 and the frequencies above 1 kT will be discussed below. From Table 1, a combination of magnetostriction, resistivity, magnetization, and thermopower is insufficient to capture all the fundamental frequencies, and F(γ2) is still missing. In contrast, the ME method demonstrates an ultra-high-sensitive performance, independently revealing all the allowed frequencies reported previously10,11,12,13,14,15,16.

Table 1 Main oscillatory frequencies in FFT analyses for ZrSiS.
Full size table

Temperature-dependent QOs of V x

We noticed that the effective mass m* can be estimated by analyzing the temperature-dependent amplitude of ME QOs in FFT according to the Lifshitz-Kosevich formula28. Therefore, the field-dependent ME signals at selected temperatures are measured and depicted in Fig. 4a. A clear damping effect of oscillation amplitudes is revealed when the temperature increases from 2 to 20 K. The corresponding FFT amplitudes of the featured frequencies also reflect this tendency, as shown in Fig. 4b. F(δ1) and F(α) are chosen to fit the thermal damping factor RT in the conventional LK formula, as shown in Fig. 4c. Their temperature-dependent behaviors exhibit metallic fermion, and the effective masses are extracted to be m*(δ1) = 0.0357 m0 and m*(α) = 0.0451 m0. Such small effective masses indicate nearly massless fermions for these two orbits, consistent with previous dHvA oscillation studies13.

Fig. 4: Field-dependent oscillations of Vx at selected temperatures and FFT analyses.
figure 4

a Background-subtracted field-dependent oscillations of Vx at selected temperatures with B // c. b FFT spectrum (5–1000 T) over a magnetic field range of 1–13 T. c Effective mass fitting of δ1 and α orbitals, revealing conventional metallic behavior consistent with the Lifshitz-Kosevich theory. d FFT spectrum (7.5–11.5 kT) across a magnetic field range of 1–13 T. e Effective mass fitting for orbitals III and IV.

Full size image

A series of ultra-high frequencies can also be observed in ME signals (Fig. 4d), which have been reported by magnetic torque and resistivity up to 35 T and down to 0.34 K10,15,16. The most pronounced peaks appear in a region of around 8.5 kT, labeled with six Roman numbers I to VI. The differences in frequency between these distinct peaks are ~F(α) except the one between peaks IV and V, which is ~F(β) −F(α). Two more peaks in a frequency region of around 11 kT are also detected but with much reduced amplitudes, and the difference between them is ~F(α). The gap in frequency between these two regions is indicated with a frequency of ~2330 T, whose origin is still controversial in the literatures10,15,16. We observed that the effective masses of orbits III and IV are fitted to be m*(III) = 0.124 m0 and m*(IV) = 0.159 m0, which hold the lowest magnitudes among these breakdown orbits (I to VI). This observation is based on their ability to persist at higher temperatures and their relatively higher amplitude compared to others. Notably, these effective masses are several times larger than those of the fundamental orbits. This significant increase in effective mass is indicative of a longer orbital length, which is a characteristic feature of giant breakdown orbits. These apparent features are overall consistent with the latest dHvA oscillation measurements16, which further validate the high sensitivity of our method under moderate field and temperature conditions.

Discusion

It is noticeably that there are several pronounced peaks between 20 T to 80 T in Fig. 3o (thermopower), while they are small in other methods. We interpret these peaks as harmonics of the δ2 frequency (~16 T), aligning with the analysis by Marcin Matusiak et al. 14, but we diverge from their interpretation, seeing the ~48 T frequency not as a new fundamental orbit but only as a third harmonic of δ2. The heightened amplitude of these peaks can be explained by the spin damping factor Rs in LK formula:

$${R}_{s}=\cos \frac{p\pi g{m}^{* }}{2{m}_{0}}.$$
(2)

Where p is harmonic index, g is the Landé g-factor. Adjusting Rs, can significantly alter the oscillatory amplitude of harmonic frequencies, sometimes even surpassing the fundamental frequency. Moreover, this effect varies between different measurement techniques. For instance, in MoSi2 with magnetic field orientation B // [001], the dHvA oscillations29 show higher amplitudes for the second harmonics compared to the fundamental frequencies. In contrast, the SdH oscillations30 under the same conditions show much weaker second harmonics. Thus, in the ac composite ME technique and the other four methods, the smaller peaks in the 20–80 T range can be attributed to this method-varying effect of enhanced harmonic oscillations.

We need to point out the limitation of our method. It would not be expected to work well for all types of materials. It tends to perform exceptionally well in non-magnetic metals, especially those that are multi-band and exhibit higher magnetostriction24,31. In metals and semimetals with magnetic orderings, the dominant magnetostriction signal from magnetic phases can overshadow the quantum oscillatory signals. This limitation also applies to type-II superconductors, where the vortex phase below the upper critical field can produce an overwhelmingly large signal32, potentially masking the quantum oscillation signals in our method.

In summary, we have demonstrated that the ac composite ME method possesses high sensitivity for studying the complex QOs in the Dirac semimetal ZrSiS. In comparison to four other typical methods with the magnetic field and temperature limited to below 13 T and above 2 K, respectively, the QOs in the ME signal provide the most complete configuration of Fermi surfaces in ZrSiS with the lowest onset field ~0.5 T, whereas the combination of all four methods still lacks some of the allowed orbitals. In particular, FFT result of the ME oscillations yields a series of clearly separated peaks above 7.5 kT. Meanwhile, the effective masses of some orbitals can be obtained by studying the temperature-damping effect of the ME QOs. The outstanding performance of the ME method may stem from the additional factor introduced by ∂λ/∂B. This factor shows significant enhancement when compared to directly measured magnetostriction. In terms of experimental techniques, the use of ac-driven and lock-in technique effectively reduces noise levels. As a result, our method enhances the sensitivity of probing quantum oscillations in magnetostriction. This increased sensitivity makes it practical for studying Dirac Fermions and other topological semimetals under moderate magnetic field and temperature conditions. This could pave the way for further investigations into Fermi surface behavior based on magnetostriction oscillations.

Methods

Sample preparation

Single crystals of ZrSiS were synthesized by the chemical vapor transport method. The stoichiometric mixture of Zr (Alfa Aesar, 99.9%), Si (Alfa Aesar, 99.9%) and S (Alfa Aesar, 99.9%) powders were sealed in an evacuated quartz tube with I2 being used as transport agent. The sealed quartz tubes were put in a two-zone tube furnace. The hot side is about 1020 oC and the cold side is 920 oC, and dwelled for two weeks. Plate-like shape single crystals with shinning surfaces were obtained. The size of the crystal is 2.72 × 2.4 × 0.36 mm3.

Transport measurements

Temperature dependent resistivity was measured using a four-probe method. All the measurements were done in a Quantum Design Physical Properties Measurement System (PPMS-9T) for 1.8 K < T < 400 K and B < 9 T.

Structural characterization

X-ray diffraction measurement of single crystal was performed by the PANalytical X’Pert diffractometer using the Cu Kα radiation (λ = 0.15406 nm) at room temperature.

Magnetization measurements

Magnetization measurements were conducted on oriented crystals in steady magnetic fields up to 14 T and at a temperature of 2 K. Either a commercial superconducting quantum interference device (SQUID) or a vibrating-sample-magnetometry from Quantum Design is used.

Striction measurements

High resolution magnetostriction measurements were achieved by using a capacitive dilatometer made of CuBe. To precisely measure the capacitance change, an Andeen Hagerling 2550 A (f = 1 kHz) capacitance bridge with driven voltage of 15 V was used.

Composite ME measurements

The setup for this method is illustrated in Fig. 2e. A PMN-PT [001]-cut single crystal, typically 0.2 mm thick, serves as a piezoelectric transducer to detect magnetostriction in the sample. Silver epoxy (H20E, EPO-TEK) is applied to function both as electrodes and as a binder. This configuration effectively converts the in-plane magnetostriction of the sample into a voltage signal with high conversion efficiency. To mitigate noise interference from our liquid-Helium-free cryostat and electrical sources, an ac lock-in technique is employed. A constant ac magnetic field (Bac, approximately 0.5 Oe at a frequency of 211 Hz) is applied in conjunction with the static field B. This is achieved using a custom coil encircling the sample. The resulting piezoelectrically-induced ac voltage, with a resolution of approximately 0.1 μV, is captured using a lock-in amplifier (OE1022, SYSU Scientific Instruments). It is crucial to maintain a consistent relative positioning between the sample and the coil to avoid signal discrepancies.

Thermopower measurements

The thermopower S was measured along the a-axis of a single crystal with the magnetic field parallel to the c-axis. The sample is bridged over two oxygen-free copper blocks to generate a temperature difference. When measuring the quantum oscillations, heating power was held constant to maintain a temperature difference of ~0.5 K for the sample.

FFT analysis

These analyses were performed using the rectangle window function to prevent the modification of oscillatory amplitude. Appropriate field regions were chosen to ensure signal-to-noise ratio, with suitable zero padding to ensure resolution of the X-axis (frequency) of the results.