Introduction

The low-temperature CDW phase in organic and inorganic solids1 serves as an important prototype for a broad range of collective phenomena involving strongly coupled degrees of freedom and spontaneously symmetry-broken ground states2. Recently, the role of CDWs in low-dimensional systems has become more generally recognized, as they are found to arise in high-T c superconductors3 and topological insulators4 as competing excitations, as well as in graphene- and carbon-nanotube-based materials5, 6. The low-energy excitations of the incommensurate CDW ground state (GS) correspond to distortions of the coupled electron and lattice density modulation (at q = 2k F), which ideally decouple into a (Raman-active) amplitude channel and (infrared-active, Goldstone) phase channel, and can be exploited as sensitive spectroscopic probes of the CDW physics. Although many studies focused on a single coupled phonon1, 7, which leads to an “amplitudon” and “phason” (the latter a soft mode at frequency ν ~ 0), in general multiple phonons can couple to the CDW electrons yielding a set of amplitude- and phase-“phonons”8,9,10,11,12. While quantum-mechanical (QM) descriptions were developed at an early stage8, 9 which predict such modes (based on fluctuations of the Fröhlich Hamiltonian7), recently a phenomenological time-dependent Ginzburg-Landau (TDGL) model was employed to account for both the temperature dependence of the modes and non-equilibrium dynamics upon photoexcitation11,12,13. The model describes the Q = 0 distortions of the CDW via classical equations of motion for the complex coordinates of the electronic order parameter (EOP, \(\tilde{{\rm{\Delta }}}\)) with linear coupling to the bare phonons (\({\tilde{\xi }}_{n}\), Fig. 1). This provides a versatile framework to include effects such as time-dependent perturbations of the potentials and higher spatial dimensions12, 13, which should find more general applicability to other collective phenomena and photo-induced phase transitions.

Figure 1
figure 1

Representation of quasi-1D charge density wave and structure in K0.3MoO3. (a) Schematic of phase-space potentials (left) and spatial modulation (right) associated with electronic CDW (order parameter, EOP) and a single phonon, with wavevectors q = 2k F and complex amplitudes \(\tilde{{\rm{\Delta }}}\) and \(\tilde{\xi }\), respectively. (b) Structure of K0.3MoO3 (T = 100 K43) including CDW distortion (exaggerated by a factor 20 for clarity), assuming q ≈ [0,0.75,0.5]. Crystal is truncated by the experimental (\(20\bar{1}\)) sample surface plane (containing [010] and [102] directions); Mo-O octahedra color-coded by inequivalent Mo atoms in the undistorted C2/m lattice.

Full size image

In such optical-pump optical-probe (OP-OP) experiments11,12,13 one specifically probes the amplitude channel via coherent oscillations in the interband permittivity. While the phase channel has been probed in certain CDW materials with terahertz (THz) and infrared spectroscopy10 including non-equilibrium dynamics14, 15, these studies did not interpret the results in terms of the QM or TDGL models. However, for further development and testing of the theoretical description, the phase channel plays a crucial role: involving the translation of the CDW condensate, it is more directly related to the non-linear transport16 and should be more sensitive to effects such as impurities, intra- and inter-chain interactions17,18,19,20. Moreover, a close inspection of the nominal QM8, 9 and TDGL11, 12 models reveals qualitative differences in the predictions specifically for the phase modes, even for T ≈ T c (as addressed below). Hence, further experimental studies of the phase-channel and theory development are required to reach a complete and unified description of the CDW physics.

Towards this goal, in the present paper we target the quasi-1D conductor blue bronze (K0.3MoO3)21, which forms an incommensurate CDW below T c = 183 K, and whose low-energy spectrum has been intensely studied via neutron scattering22, 23, Raman24, 25, and far-infrared10, 26 spectroscopy. Moreover, the non-equilibrium response following femtosecond optical excitation has been studied recently using optical12, 13, 27,28,29, photoelectron30 and electron-diffraction31 probes, including application of the TDGL model to the amplitude modes11,12,13. For the phase modes, while previous GS (far-)IR measurements covered a wide frequency range10, 26, these were confined to low temperature (≤10 K). Moreover, to our knowledge, no studies of the non-equilibrium response of the phase-phonons have been reported (although a signature assigned to the phason below 50 GHz was resolved in OP-OP measurements32). Here we apply optical-pump THz-probe spectroscopy (OP-TP), which probes the phase-phonons via the broadband, transient complex conductivity, and reconcile both the non-equilibrium dynamics and T-dependence of the bands with a generalized TDGL model.

Results

Non-equilibrium dynamics

The OP-TP experiments were performed with both optical-pump (hν ex = 1.6 eV) and THz-probe pulses (time delay τ) normally incident on the sample (Fig. 2(a)), and the reflected THz field was measured with time-domain sampling (see Methods section) In Fig. 2(b) we show examples of the reflected reference pulse E 0(t) and photoinduced differential signal ΔE(t, τ) = E(t, τ) − E 0(t) at τ = 0 (with the THz field polarized parallel to the crystal b-axis) for T = 50 K and an excitation fluence of F ex = 550 μJ cm−2. The full 2D signal built up from a sequence of measurements vs. τ is shown in Fig. 2(c). One sees that the temporal onset of ΔE is shifted somewhat from the main peak of E 0 and exhibits a ringing signature composed of multiple beating frequencies due to a perturbation of the GS phase-phonons. The corresponding differential reflectivity spectra Δr(ν, τ)/r 0(ν) are presented in Fig. 3(a–d). Accounting for the short excitation depth D ex = 1/α ex compared to the THz probe wavelengths, one can show that the complex conductivity change Δσ(ν, τ) is essentially in-phase with Δr/r 0 (see Supplemental Information), and hence spectral features in Δr/r 0 correlate closely with those in Δσ (with minor distortions due to a pre-factor containing the GS permittivity). The spectra in Fig. 3(a) contain a sequence of derivative-shape features in the range \(\mathop{ > }\limits_{ \tilde {}}\)1.5 THz for Re{Δr/r 0} (τ = 2 ps), whose polarity implies a photoinduced blue-shift of existing bands. The most dominant feature (with a zero-crossing at 1.78 THz) occurs about a frequency where a weak spectral modulation was indeed observed previously10 but not assigned as a phase-phonon. From additional GS THz measurements and the analysis below, we deduce that the signatures arise from three main bands (labeled A–C from hereon, as shown in Fig. 3(a)), where bands (B, C) correspond to those previously found in FTIR measurements10. The weak signature at low frequency can be fit well with a Drude model (see below), which we assign to the generation of mobile carriers. From the full 2D data (Fig. 3(c,d)) one sees that all features persist for delays beyond 50 ps, after initial relaxation in the first few ps.

Figure 2
figure 2

Time-domain field response to optical excitation. (a) Experimental OP-TP reflection geometry. (b) Detected temporal electric fields (T = 50 K, pump fluence F ex = 550 μJ cm−2): Reference field E 0(t) (reflected signal without optical excitation); differential field change ΔE(t, τ) for pump-probe delay τ = 0. (c) Full 2D map of ΔE(t, τ). Note that the data is plotted with two different scales for the ranges τ ∈ [−2, 5] ps and τ ∈ [5, 70] ps, respectively.

Full size image
Figure 3
figure 3

Time-resolved phase-phonon spectra and model analysis. Differential field reflectivity spectra Δr/r 0 ∝ Δσ/(ε r − 1) at T = 50 K (left graphs Re{Δr/r 0}; right graphs Im{Δr/r 0}). (a,b) For pump-probe delay τ = 2 ps, including both experimental data (points) and fitted model (see text). The GS frequencies of the three main bands (labeled A–C) from the OP-TP analysis are included as vertical lines. (c,d) Full experimental data Δr(ν, τ)/r 0(ν) obtained from Fourier analysis of each row in Fig. 2(c); (e,f) corresponding model fits.

Full size image

In order to disentangle the band contributions and extract parameter kinetics vs. τ, we fit the full complex data based on photoinduced modification of three Lorentzian bands and a Drude contribution (see Methods section). The resulting fitted spectra are shown in Fig. 3(e,f) (and also as curves in Fig. 3(a,b) for τ = 2 ps), which reproduce the experimental data well, except for the additional broadening seen in the first 2 ps. Simulations for a single perturbed phonon show that this transient broadening can be attributed to the response of a time-non-stationary medium33 due to sub-ps dynamics (see Supplemental Information), and hence cannot be reproduced using the quasi-stationary spectral analysis here. In Fig. 4(a–d) we plot selected fitted band parameters as a function of delay τ, where the blue-shifts for all bands A–C are clearly evident (ν 0n (τ), Fig. 4(a)). For band A, this is accompanied with a transient reduction in band strength (S A, Fig. 4(b)), and additional broadening (after initial development in the first few ps, Fig. 4(c)). As expected from inspection of the low-frequency region in Fig. 3(c), the fitted Drude plasma frequency rises to a fairly constant plateau at ν pl = 2 THz (Fig. 4(d)) with a scattering time τ s ~ 90 fs. Based on the low-frequency value of \({\varepsilon }_{{\rm{r}}}(\nu \to 0) \sim 100\), ν pl corresponds to an excitation density \({N}_{{\rm{ex}}}/({m}_{{\rm{eff}}}/{m}_{0}) \sim 5\times {10}^{18}\,{{\rm{cm}}}^{-3}\). As these photoexcited carriers may well populate orbitals distinct from those of the CDW condensate, we do not necessarily associate the effective mass here with that of the CDW (\({m}_{{\rm{eff}}}^{{\rm{CDW}}} \sim 300\cdot {m}_{0}\) 21).

Figure 4
figure 4

Parameter kinetics for phase-phonon bands and free-carrier response. Results of fitting Δr/r 0 data in Fig. 3 (T = 50 K) vs. pump-probe delay τ. (a) Band frequencies ν 0n for the three bands n = A,B,C. GS values obtained during global fit indicated by dashed horizontal lines. (b,c) Band strength and bandwidth (FWHM) for Band A, respectively. (d) Fitted Drude plasma frequency ν pl for free-carrier contribution. Note that the fluctuations at later delays arise due to small drifts and noise during the measurement run, and any structure should rather be attributed to the non-monotonic order of delays used for the measurement run. (e) Experimental kinetics of Im{Δr/r 0} at two selected frequencies, as indicated, and global bi-exponential fit (solid curves). (f) Corresponding residual from fitting in (e). Representative confidence intervals in the fitted frequency shifts and plasma frequency indicated by the error bars in (a,d), respectively (see Fig. 5(b,d) for errors in the GS parameters).

Full size image

To extract kinetic time scales, we fit the experimental data directly, i.e. the kinetics of Im{Δr/r 0} for ν = 1.76 and 2.28 THz (the two peaks in Fig. 3(b)) as shown in Fig. 4(e). A global bi-exponential fit yielded time constants of τ 1 = 5.1 ps and τ 2 = 60 ps, in addition to a small but significant, long-lived offset. The values τ 1,2 also reproduce the band-shift kinetics well (solid curve in Fig. 4(a) for band A). The fit residuals (Fig. 4(f)) exhibit roughly two oscillation cycles with period ~2 s followed by a weaker oscillation with period ~10 ps. This additional slow modulation of the system, which is in-phase for both bands, presumably results from coherent charge/lattice oscillations polarized perpendicular to the surface. As discussed below, we attribute the blue-shifts to a reduction of the e-ph coupling due to the presence of the free carriers (i.e. a reduction of the red-shift of the dressed phonons). As the plasma frequency ν pl(τ) shows no such decay kinetics, this would imply that after the initial reaction of the phonons to the free charges, the system can reconfigure in such a way that the e-ph interactions can be partially reestablished. This process may again involve relaxation of excitation distributions perpendicular to the surface. Compared to OP-OP studies11, 28 of the amplitude-phonons (and interband electronic contributions), our value of τ 1 is consistent with a component in the range 5–10 ps11, 12, 29 which was assigned to a second stage of CDW recovery, whereas the initial sub-ps CDW recovery (on the same time scale as the period of the THz waves here) would rather contribute to the transient broadening, as seen in our data. A persistent offset in the (electronic) OP-OP signals was attributed tentatively to residual cooling28, although our results here indicate that a long-lived free-carrier population may well contribute.

Temperature dependence and extended Ginzburg-Landau model

We turn now to the temperature dependence of the phase-phonons, obtained from measurements as per those above for T = 50, 80, 130 and 160 K (see Supplemental Information for full data sets), which yielded the fitted band frequencies ν 0n and widths Δν n shown in Fig. 5 (right) for both the GS, and excited state for τ = 2 ps. One sees a clear softening for band B and C, but rather a stiffening for band A, with increasing T, accompanied with significant broadening for band A. The GS data for bands B, C from ref. 10 at T = 6 K are also included, and are reasonably consistent with an extrapolation of our parameters (recall that band A was not assigned in ref. 10).

Figure 5
figure 5

Temperature dependence of amplitude- and phase-phonon band parameters. Amplitude- (left graphs) and phase-phonons (right) for T < T c; center frequencies ν 0n (top graphs) and bandwidths Δν n (FWHM, bottom), for bands n = A (red), B (magenta) and C (blue). Literature data from neutron scattering (\(\blacktriangle\) 23), Raman scattering (\(\{\backslash triangledown\}\) 24; \(\diamond \) 25), optical-pump optical-probe reflectivity (\(\blacksquare\) 11, 12), and far-infrared FTIR measurements (○10, T = 6 K). THz TDS data (present work) from global analysis of OP-TP data (for T = 50, 80, 130 and 160 K) for both GS (●) and τ = 2 ps after excitation (☆), including confidence intervals for the GS parameters. Also included are model curves from the TDGL model, with parameters fit to the amplitude phonons in ref. 11, both without (solid) and with (dashed) pinning/scattering contributions from impurities (with Ωi(0)/2π = 9.3 THz and g i = 1.9, see text; note that this only affects the predictions for the phase channel, so both curves are superposed for the amplitude channel).

Full size image

Literature data for the amplitude modes are also shown in Fig. 5 (left), comprising GS neutron- and Raman-spectroscopy for band A, and OP-OP data for all three bands (sources cited in caption), which are reasonably consistent (except for the deviation between neutron and OP-OP data for ν 0A approaching T c). Also included are the theoretical curves based on the TDGL model adopted from refs 11, 12, which demonstrates the very good agreement obtained by those authors (even for \(T\ll {T}_{{\rm{c}}}\)). However, the nominal TDGL model predicts T-independent phase-phonon bands (which coincide with the amplitude modes for T → T c; solid horizontal lines in Fig. 5(b,d)), which clearly is not consistent with our results away from T c. While any GL treatment2 is only a first-order expansion about T ≈ T c, it is not clear why it should break down here specifically for the phase channel. Hence we considered possible T-dependent effects which should only affect the phase channel (Δ2, ξ n2) in the TDGL equations. A straightforward generalization which reproduces the observed trends for all bands is based on impurity effects (see Supplemental Information), which dominantly affect the CDW phase channel via pinning17, 34, 35 and scattering18, 36. We incorporated these via a pinning potential \({U}_{{\rm{i}}}=-\,\tfrac{1}{2}{{\rm{\Omega }}}_{{\rm{i}}}^{2}(T){{\rm{\Delta }}}_{2}^{2}\) and damping term γ 2 → γ + γ i(T), taken to depend on the relative equilibrium EOP amplitude δ 0(T) = (1 − T/T c)1/2 as per \({{\rm{\Omega }}}_{{\rm{i}}}^{2}(T)={{\rm{\Omega }}}_{{\rm{i}}}^{2}(0){\delta }_{0}^{n}(T)\) and \({\gamma }_{i}(T)/\gamma ={g}_{{\rm{i}}}{\delta }_{0}^{n}(T)\), where the exponent n = 2 was required to obtain reasonable agreement with the data. As can be seen in Fig. 5(b,d), the generalized model reproduces the qualitative trends in the phase-phonon frequencies ν 0n (T) (including the crossing of ν 0A through \({\nu }_{0{\rm{A}}}^{(0)}\) at T ~ 120 K) and significantly improves the predicted trend for Δν A(T), while maintaining all other TDGL parameters used for the amplitude-phonons. The exponent n = 2 deviates from the value n = 1 which is commonly taken for the potential of a single impurity: \({U}_{{\rm{i}}}\propto {{\rm{\Omega }}}_{{\rm{i}}}^{2}\propto {\delta }_{0}\) 17, 34, 35. However, QM treatments have predicted impurity pinning forces which grow faster than Δ0(T)37, 38 due to competition between the CDW and formation of Friedel oscillations about each impurity. Also, a divergence in phase damping with decreasing T was deduced from non-linear transport experiments16 and was discussed in terms of screening39 and impurity-induced phase deformations of the CDW, which enhance phase scattering as compensating thermal carriers are frozen out (which dominates over the damping due to the corresponding CDW-carrier scattering19 which would instead grow with T).

Considering the photoinduced band changes, one can see that at least qualitatively, the observed blue-shifts are consistent with a reduction of the couplings m n between EOP and bare phonons11 (shifting the frequencies ν 0n in the direction of \({\nu }_{0n}^{(0)}\)), which we attribute to the presence of the long-lived free carriers. The fact that ν 0A shifts even above \({\nu }_{0{\rm{A}}}^{(0)}\) for T < 160 K (Fig. 5(b)) may be due to an additional partial suppression of the impurity pinning (recall that the solid lines are the predicted frequencies in the absence of pinning). This raises the question whether the amplitude-phonons in OP-OP measurements are also blue-shifted11,12,13. A study vs. fluence F ex 29, however, revealed a small red-shift at high fluence. One mechanism which could account for this different behavior involves relaxing the strict phase coherence of the coupling between EOP and bare phonons after excitation,

$${U}_{n}\to -\,{m}_{n}{\rm{\Delta }}\cdot {\xi }_{n}[(1-\eta )\,\cos \,(\phi -\chi )+\eta ],\,(0\le \eta \le 1)$$

e.g. due to the presence of inhomogeneous distortions of the CDW. This modification only reduces the effective coupling for the phase channel (see Supplemental Information), consistent with the experimental findings.

Discussion

While the extended TDGL model allows us to assign and account for all three phase-phonons here, it is important to note that this contrasts with the predictions of the nominal QM theory (even for T ≈ T c)8, 9, where one of the renormalized modes is driven to a near-zero-frequency “phason” (which was experimentally assigned10 at ν ~ 0.1 THz for K0.3MoO3). This discrepancy manifests much more clearly for the phase channel, as all amplitude-phonons are predicted to remain at finite frequencies in both models. A key difference appears to lie in the treatment of electronic damping in the QM model (see Supplemental Information), which in the phenomenological TDGL treatment is dictated by free parameters chosen to agree with experiment (corresponding to an overdamped EOP11, 12). This motivates further theory development to include additional interactions in deriving the spectral response function. A first step was taken in this direction with extensions of the QM model to include long-range Coulomb effects19, 20, which were shown to strongly affect the phase channel, causing spectral weight to be transferred from the phason to an additional finite-frequency resonance.

This study demonstrates how coherent THz spectroscopy can provide unique insight into the underlying CDW physics, such as impurity and screening effects which predominantly affect the phase channel. The inclusion of these effects within a generalized TDGL model can account for the T-dependence and photoinduced dynamics observed here for K0.3MoO3. Given that the TDGL model provides a highly practical framework for interpretation of (non-equilibrium) experiments11,12,13, we strongly advocate further investigations into its applicability and correspondence to microscopic QM treatments.

Methods

The single crystals of K0.3MoO3 were grown according to the temperature-gradient flux method40, reaching facial dimensions up to 4 × 5 mm2, and the crystal structure and orientation were confirmed by x-ray diffraction on a small single crystal. The crystal surfaces were sufficiently flat/smooth to be used in the experiments as-grown.

The terahertz time-domain spectroscopy system used for the optical-pump THz-probe experiments is shown in the Supplemental Information. The system is based on a 1-kHz Ti:Al2O3 amplifier laser (Clark-MXR CPA-2101, λ ex = 775 nm, pulse duration ~150 fs FWHM). The THz pulses were generated by a large-area InAs surface emitter, while a 0.5-mm-thick 〈110〉-cut ZnTe is used for electrooptic detection (frequency range ~0.2–3 THz). The blue bronze samples were mounted in a liquid-helium cryostat (Oxford Microstat) and could be oriented with the conductive b-axis either parallel or perpendicular to the p-polarized THz probe field. A thin polypropylene film (50-μm-thick, Goodfellow PP301350) was used as the window to minimize any distortion of the THz beam (and two-photon absorption effects for the optical excitation pulse). In order to realize a normal-incidence reflection geometry, a two-way beamsplitter (0.5-mm-thick high-resistivity Si wafer) is placed in the THz beam before the sample. The optical excitation beam for the sample passes through a small hole in the sample focusing mirror (off-axis paraboloidal mirror, effective focal distance f eff = 101.6 mm), brought to a weak focus on the sample (with a spot diameter of ~2 mm) using a telescope. While measurements were performed with the THz field polarized both parallel and perpendicular to the b-axis, the phase-phonon signals presented here were observed only in the former case. Also, measurements were made with fluence in the range F ex = 140–720 μJ cm−2, although we present here only data for F ex = 550 μJ cm−2 which yielded a favorable compromise between signal strength and saturation effects29. The mechanical delay stage for THz detection was placed in the emitter arm, such that the time-variable t corresponds to a fixed pump-probe delay τ for each point in the acquired THz waveform41. Besides yielding a 2D pump-probe signal ΔE(t, τ) = E(t, τ) − E 0(t) with a synchronous zero-delay for all t, this also rejects any pump-only emission signals in the resulting spectral data ΔE(ν, τ) (which were also observed). A dual-chopping scheme with two lock-in amplifiers was employed to synchronously measure ΔE and E 0 42, which was essential to track phase drifts in the reflected THz field due to unavoidable longitudinal movement of the cryostat (of some 10 μm) during measurements.

The fit model for each complex field reflectivity spectrum Δr(ν, τ)/r 0(ν) in the OP-TP data was based on assuming three Lorentzian bands and a Drude contribution (as well as a frequency-independent background permittivity ε br to account for higher-lying resonances), whose GS parameters (frequency ν  0n , width Δν n and strength S  n ) are modified by the optical excitation to \({\nu }_{0n}^{^{\prime} }(\tau )\), \({\rm{\Delta }}{\nu }_{n}^{^{\prime} }(\tau )\), \({S}_{n}^{^{\prime} }(\tau )\) and \({\varepsilon }_{{\rm{br}}}^{^{\prime} }(\tau )\). By definition, ν 0,D ≡ \({\nu }_{0,{\rm{D}}}^{^{\prime} }\) ≡ 0 for the Drude contribution, and we assumed no free carriers in the GS (S D = 0). Note that we also tested the use of Fano lineshapes (as was applied in the vicinity of a single band in an OP-TP study on 1T-TaS2 14), especially as this could be appropriate for the proposed interaction with free carriers. However, the Fano interaction (which essentially only mixes the real and imaginary parts of the band conductivity) could not reproduce such large shifts, especially for band A. Due to baseline/resolution issues in the GS THz-TDS results, we included global GS band parameters in the fit procedure (the fitted GS band strengths were constrained to remain comparable to those from GS THz-TDS data). In order to implement a robust and practical fitting scheme including the GS parameters, we first selected spectra from a small number of delays (τ n  = 2, 10 ps), and performed the fitting of all parameters. Thereafter, the GS parameters were held fixed and each spectrum was fitted sequentially for the excited-state parameters. At each iteration, we first used an approximate expression for Δr/r 0 valid for \({D}_{{\rm{ex}}}\ll {\lambda }_{{\rm{THz}}}\) (Eq. 1 in Supplemental Information) to estimate the model differential field reflectivity. While this equation is nonlinear in the GS complex permittivity, it is linear in the differential change. This allowed to determine the ES band strengths \({S}_{n}^{^{\prime} }(\tau )\) via generalized regression, which relieved the fitting algorithm of a fraction of the search parameters and facilitated the fitting procedure. Then the exact expression for Δr/r 0 was calculated and used to compute the current misfit. To avoid stagnation, the nonlinear optimization was based on Covariance Matrix Adaptation Evolution Strategy (CMA-ES), followed by local refinement with the Nelder-Mead simplex method. Confidence intervals in the fitted parameters were obtained from the numerical Hessian of the misfit function, corresponding to ±1σ. Note that the excitation depth D ex was also fitted (as the nominal value from linear measurements could not account for the experimental data, see main paper). However, after the optimum value was determined with the data for T = 50 K, it was kept fixed for all other temperatures T = 80, 130, 160 K.

Data availability

The raw data and analysis results of the current study are available from the corresponding author on reasonable request.