Phase-channel dynamics reveal the role of impurities and screening in a quasi-one-dimensional charge-density wave system

Charge density waves (CDWs), i.e. the periodic spatial modulation of coupled electronic and lattice density, are ubiquitous in low-dimensional conductors and have taken on renewed relevance due their role in state-of-the-art materials, e.g. high-T c superconductors, topological insulators and low-dimensional carbon. As CDWs are described by a complex order parameter to represent both the amplitude and phase, they are formally analogous to BCS superconductors and spin-waves, providing a prototype of collective phenomena for the further development of field theories and ab-initio calculations of complex solids. The low-energy excitations are mixed electron-phonon quanta which ideally separate into an amplitude and phase channel, and provide a sensitive probe of the ground state and non-equilibrium dynamics, including ultrafast photoinduced phase transitions. While recent studies of the amplitude modes have brought substantial progress aided by a phenomenological Ginzburg-Landau framework, we focus here on the phase modes using ultrafast terahertz spectroscopy. Experiments on K0.3MoO3 provide a more complete picture, and reveal a high sensitivity to interactions with impurities and screening effects from photogenerated carriers, both of which can be accounted for by generalizations of the model. Moreover, our considerations emphasize the need to revisit the treatment of inherent electronic damping in quantum-mechanical CDW theories.


Lock-in
B. Derivation of approximate relation between differential reflectivity and conductivity for an exponential excitation profile The complex field reflection coefficient for a sample where the excitation depth is small compared to the probe wavelength (as is the case for blue bronze with near-infrared excitation and THz probing) deviates significantly from the Fresnel formula for a homogeneous medium [1,2]. For a pump-induced change in the sample complex permittivity ∆ε r (ω,z) = ∆ε r (ω)e −z D ex (where D ex = 1 α ex is the excitation depth) the reflectivity can still be expressed in the form r = (1 − n ′ ) (1 + n ′ ) where n ′ = n + ∆ε r X β (ξ ), n = √ ε r is the background (unpumped) refractive index, and X β (ξ ) = I β +1 (ξ ) I β (ξ ) is the ratio of modified Bessel functions with β = 2iωnD ex c and ξ = 2iω ∆ε r D ex c. For D ex ≪ λ THz , one can approximate the expression for X β (ξ ) to first-order as: X β (ξ ) → 1 2 ξ (β + 1). Substituting this into the expression for r, noting that the ground-state reflectivity r 0 = (1 − n) (1 + n), one can readily derive (assuming ∆n = n ′ − n ≪ n): where ∆r = r − r 0 , k 0 = ω c, and ∆σ = iε 0 ω∆ε r is the pump-induced change in complex conductivity. The validity of these approximations was thoroughly checked in numerical tests with parameters appropriate to the experiments presented in the main paper. Note that if one were instead to use a simple Fresnel reflection treatment assuming homogeneous excitation vs. depth, this introduces an imaginary unit in the coefficient relating ∆r r 0 and ∆σ , such that real and imaginary parts would be erroneously exchanged in calculating ∆σ from the data.   (Fig. 3)). Experimental data (left two columns) and Lorentzian-Drude fits (right two columns); real and imaginary parts of ∆r r 0 as indicated.

D. Application of the time-dependent Ginzburg-Landau treatment
Here we consider generalizations of the time-dependent Ginzburg-Landau (TDGL) model prescribed in [3,4]. For completeness, we first briefly summarize this, in order to clarify the notation and modifications that follow. We consider the potential function, in terms of the complex electronic order parameter∆ = ∆e iϕ = ∆ 1 +i∆ 2 and complex bare phonon coordinatesξ n = ξ n e iχn = ξ n1 +iξ n2 (n = 1...N) (where all coordinates refer to the complex envelope amplitudes of the q = 2k F components) as per: is the Mexican hat potential, U ξ n = 1 2 ω 2 0n ξ 2 n represents the elastic energy stored in the bare phonon mode n with frequency ω 0n , and U c = −m n (∆ 1 ξ n1 + ∆ 2 ξ n2 ) = −m n ∆ ⋅ ξ n cos(ϕ − χ n ) is the linear coupling term (summations over n are left implicit). This has the equilibrium solution is the renormalized critical temperature. Taking ϕ 0 = χ n0 = 0 as the equilibrium phase (at first, arbitrarily, as there is no pinning term and hence U depends only on (ϕ − χ n )), and calculating the Hessian matrix of the potential about∆ 0 = ∆ 0 yields the linearized equations of motion [3,4]: where∆ 1 ≈ ∆ − ∆ 0 and∆ 2 ≈ ∆ 0 ϕ represent the amplitude and phase deviations from equilibrium (likewise forξ n1 ,ξ n2 ), and we have added phenomenological damping constants γ 1,2 for∆ 1,2 . Note that while in a classical treatment one has γ 1 = γ 2 = γ [5], we allow here for γ 1 ≠ γ 2 to account for different quasi-particle scattering channels. Moreover, the variables are normalized to unit mass and we do not include any inherent damping for the bare phonon modes (which did not assist in fitting the phase-phonon data). The equations for the amplitude (Raman-active) and phase (IR-active) channels are decoupled, and converge to the same set of solutions for T → T c (corresponding to the T -independent phase channel).
In the overdamped limit for the EOP (as adopted in [3,4]), i.e. ∂ 2 t ∆ j ≪ γ∂ t ∆ j (which is equivalent to neglecting the inertial mass), Eq.s 3(a-b) become: where κ j = γ −1 j is the scattering time. As per [3,4], one can solve these linear equations with the ansatz ∝ e λt for the eigenvalues λ k = −Γ k 2 + iΩ 0k (and corresponding eigenvectors, which reflect the relative contribution of the EOP and each phonon to the mode), which yields either N + 2 (general damping, Eq. 3) or N + 1 (overdamped limit, Eq. 4) modes for each channel (after discounting for the complex-conjugate solutions with Ω 0k < 0). For each channel, one obtains a set of N finite-frequency modes which, at least for the parameters used for K 0.3 MoO 3 in [3,4], each contain a dominant contribution from one of the bare modes, and hence have frequencies Ω 0k only moderately shifted from a given bare phonon at ω 0n (although in general, stronger mixing and avoided crossings can result [4]). In the overdamped limit, the remaining eigenvalues have I(λ ) = 0, and also R(λ ) = 0 for the phase channel (and amplitude channel for T → T c ). The latter was associated in [3,4] with a true soft mode of the system [6]. As the authors asserted, this indeed represents a distinct physical picture to that from earlier quantum-mechanical (QM) treatments, as discussed below.   [3] and as the basis for the IR-active data measured in the present paper.
In order to provide a starting point for the analysis of our phase-phonon data vs T , we first re-fitted the model curves for the Raman band data ( Fig. 5(a,c) in main paper), i.e. the frequencies and bandwidths extracted from [3] (note that only the fitted values of the parameters κm 2 n were given there). This yielded a fitted set of parameters as shown in Table I, where the values of κm 2 n are in agreement with those given in [3].

Phenomenological impurity pinning/scattering vs. temperature
As discussed in the main paper, the nominal TDGL model summarized above predicts temperatureindependent band parameters for the phase-phonons, whereas a clear T -dependence is observed in the experimental THz data ( Fig. 5(b,d), main paper). Hence several physical generalizations/modifications of the model were investigated, in order to reconcile this T -dependence, which would still retain the good agreement with the literature amplitude-phonon data. This included (i) general damping (c.f. Eq. 3, i.e. including the underdamped regime); (ii) a temperature dependence for β -although this does not actually affect the T -dependence of the phonon data, rather only the magnitude of the equilibrium coordinates and renormalization of the critical temperature T c0 → T c . (iii) Our efforts to generalize the TDGL potential to include gradient terms of the complex order parameter (as suggested in [3]), in particular the term ∂ z∆ 2 [7,8], did not yield any coupling between the amplitude-and phase-channels (at least to first-order, once reducing the spatial equations to the q = 2k F -components∆ andξ n ). (iv) A finite pinning potential for the CDW phase [8][9][10], which was omitted in [3,4] presumably as it should not strongly affect the amplitude channel and hence was not necessary to treat the behavior of the Raman-active bands. (v) T -dependent damping γ 2 for the EOP phase (∆ 2 ). In this case, one could consider two opposite trends: i.e. γ 2 (T ) increasing with T , e.g. due to friction (scattering) with an increasing number of normal carriers thermally excited across the CDW gap [11]; γ 2 (T ) decreasing with T , e.g. due to the loss of thermal carriers which may screen impurity/Coulomb interactions between charges in the CDW condensate and suppress the associated scattering processes [12,13].
We found that including effects (iv) and (v) (with γ 2 decreasing with T ) allowed one to reproduce the (semi-quantitative) behavior of the phase-phonon band data vs. T , i.e. with band A stiffening (and crossing its bare mode frequency ν (0) 0A at intermediate T ) while bands B,C soften as T → T c , and a significant reduction in the bandwidth of band A at low T . These effects were implemented as follows. For the impurity pinning, we treat a single strong-pinning center (which could be the consolidated average of several impurity potentials in each local region) and add to the function U (Eq. 2) a harmonic potential U i = 1 2 Ω 2 i (T )∆ 2 2 [10,14] (taking the origin of the impurity at ϕ = 0). As both impurity pinning and scattering should depend on the (equilibrium) amplitude of the EOP, we took both the restoring force Ω 2 i (T ) and damping γ 2 (T ) to depend on ∆ 0 , i.e. Ω 2 i (T ) = Ω 2 i (0)δ n 0 (T ) and γ 2 (T ) = γ + γ i (T ) with γ i (T ) γ = g i δ n 0 (T ), where δ 0 (T ) ≡ ∆ 0 (T ) ∆ 0 (0) = (1 − T T c ) 1 2 is the normalized EOP amplitude. Tests showed that the value n = 2 for the exponent was required to obtain reasonable agreement with the data (the physical significance of this is discussed in the main paper). The best fit of the phase-phonon data corresponded to the values Ω i (0) 2π = 9.3 THz and g i = 1.9 (using the same values for the other parameters above, and hence main-taining the fit to the literature amplitude-phonon data). The separate and combined effect of each term is illustrated in Fig. 3 (c.f. Fig. 5 in the main paper). In particular, one sees that the additional damping causes the eigenvalues to approach those of the bare modes as T decreases. This reduction in mixing between the (increasingly damped) EOP and the bare phonons can also be understood by considering the mechanical analogy of coupling of oscillators to an overdamped "bath". Also, the pinning potential is seen to play the dominant role for band A in terms the crossing of ν

EOP-phonon coupling with both coherent and incoherent contributions
As mentioned in the main text, the present results and those from the literature [15] indicate that only the phase-phonon bands are strongly affected (blue-shifted) by photoexcited free carriers. The issue was discussed in how to reconcile such a situation within the TDGL model. As asserted in the main paper, one explanation involves a loss of equivalence in the EOP-phonon coupling m n for the amplitude and phase channels, due to the coupling losing a strict coherent dependence on the relative phases. Here we describe how this can be included in the TDGL model quantitatively.
The coupling term U c between the EOP and each bare phonon is replaced by where 0 < η < 1 represents the loss of coherence. Physically, this could arise in the photoexcited state due to local fluctuations, whereby the coexistence of the EOP and phonon amplitude still yields a stabilization of the local energy density, only that the strict local phase relation is suppressed. Using Eq. 5 in Eq. 2, one still obtains the same result for the equilibrium given above (∆ 0 ,ξ 0n ). The Hessian about the equilibrium is however altered, leading to the equations of motion: Comparison with Eq.s 3 shows that the amplitude channel is unaffected (Eq.s 6(a,c)), while the effective coupling for the phase channel is reduced (essentially by a factor (1 − η), although the first term in Eq. 6(b) corresponds rather to m n → m n √ 1 − η). Hence, this modification predicts that the phase-phonons should be shifted toward their bare frequencies upon photoexcitation, as per the dominant trend seen in the data (see Fig. 5 (b) in the main paper).

E. Comparison of the TDGL model with quantum mechanical models
It is instructive to compare the qualitative predictions of the TDGL model with those from earlier quantum-mechanical treatments (QM) in the literature (based on the Fröhlich Hamiltonian). The extension of the original LRA model [16] to N>1 bare phonons [17] provides a closed-form solution for the Q = 0 amplitude-and phase-phonon spectra for T → 0. This yields N renormalized modes for each channel, with 0n (as well as the single-particle gap near 2∆ in the conductivity spectrum). For typical values of the e-ph coupling coefficients, one has: (i) the lowest amplitude mode ("amplitudon") is most strongly red-shifted; (ii) the lowest phase mode ("phason") is driven to ν = 0 (or close to it, if pinning is included [17]) leaving only N−1 renormalized phase-phonons at finite frequency. This contrasts to the predictions of the TDGL model, where, as detailed above, one has N finite-frequency modes also for the phase channel in the overdamped limit (and an additional mode at ν = 0). Hence in the nominal QM theory, we would have consider a re-assessment of our assignment for band A (the lowest-frequency phase-phonon), and, for that matter, the basis of the TDGL model. We note that subsequent extensions of the QM model to include long-range Coulomb terms [18,19] predict that the phason resonance can actually break up into both acoustic and optical resonances, the latter of which could provide an alternate assignment for our band A (although the T -dependence of the example results in [18] do not correlate with the frequency/bandwidth behavior of our band here). Interestingly, if one reduces the phase damping γ 2 in the TDGL model to the underdamped regime, we find that one obtains a set of phase-phonons qualitatively consistent with the QM model (i.e. with an additional phase-mode at ν = 0 and with non-zero damping), although this does not readily aid in interpreting our results here. Conversely, one could consider that the interactions included to calculate the QM spectral response functions do not adequately predict the phase damping, and should rather converge to the TDGL result if this were incorporated. However, in this case, the QM theory would no longer predict the phason at ν ∼ 0, where a band with a peak at ν ∼ 100 GHz has indeed been assigned in experimental studies [20]. Clearly this issue deserves further theoretical investigation, as we assert in the main paper.

F. Simulations of transient broadening of phonon response around delay zero
As presented in the main paper, some of the initial spectral features in the reflectivity spectra ∆r r 0 (i.e. for τ ≲ 2 ps) could not be fitted using the quasi-stationary Drude-Lorentz dispersion model, where a form of transient broadening occurs for the Lorentzian bands. This phenomena is primarily due to the frequencymixing which occurs when the THz probe pulse interacts with a system whose polarization response changes on the time scale of the THz cycle period. While these effects have been studied theoretically for a Drude response (e.g. [21]), we did not find any exposition for the perturbation of (narrow) Lorentzian bands, where the temporal response (free-induction decay) is longer and hence temporal-spectral effects can be more severe. Here we present the simulation results of such OP-TP reflection experiments, for the case of a single Lorentzian band which undergoes a blue-shift in its instantaneous response on a short time scale, to demonstrate how the experimentally observed broadening (as identified in the main paper) manifests.
While the finite-difference time-domain (FDTD) simulations are a proven method to simulate pumpprobe experiments with rapidly changing response functions [22,23], these can be time-consuming, especially for batch runs (i.e. vs. pump-probe delay τ). Instead, here we develop and apply a model for a homogeneously excited thin-film between air (refractive index n 1 ) and substrate (n 3 ) which allows timeintegration of only a small number of field/polarization values. The treatment is based on applying the time-domain boundary conditions for the fields at both interfaces and expanding the field development to first-order in the sample depth (z ∈ [0,L]). The resulting coupled equations for the reflected (E r ) and transmitted (E t ) fields, in terms of the known incident field (E i ) are where each P{E n } is the co-integrated time-domain solution to the corresponding polarization equation where Γ(t), ω 0 (t) and S(t) are the time-dependent damping, resonant frequency and strength of the electronic oscillator (whose temporal development is initiated by a pump pulse centered at t = −τ). These equations are readily generalized to multiple Lorentzian bands and the addition of an instantaneous response due to ε br . Their validity was checked by comparison of numerical data for a time-stationary medium with the analytic formula for the reflected/transmitted fields. To simulate the experimental data, Eq.s 7-8 are integrated for a large set of τ-values. The resulting raw data (integrated vs. laboratory time t lab ) are then interpolated onto the experimental time base t (which is skewed relative to t lab due to the use of delay stage in the THz emitter beam path [24]), and after time-windowing and Fourier transformation ∆r(ν,τ) r 0 (ν) is calculated. As a representative simulation for comparison with the experimental 2D signals ∆r(ν,τ) r 0 (ν) (see Fig. 3(c,d) in main paper), we consider the case of a single Lorentzian band, whose frequency and strength are perturbed by excitation with a Gaussian pulse with duration 150 fs, as shown in Fig. 4(a). Also shown for comparison is the "fictitious" simulated data ∆r r 0 obtained if one uses a stationary response for each "frozen" value of the Lorentzian band parameters (which would be the case for a much slower transition in the parameters). Clearly one observes the dynamic broadening artifact in the simulated data (especially during the first ∼2 ps), which bears a close resemblance to those in the experimental data (both for the real and imaginary parts of ∆r r 0 ). Note that this is an inherent effect in light-matter interaction close to the uncertainty limit (i.e. during rapid evolution compared to the oscillation period, spectral resonances must become broadened), and accentuates that OP-TP transients on a sub-ps time scale cannot be interpreted directly in terms of a quasi-stationary picture. A more rigorous approach is to consider the response of the perturbed system rather in the dual-frequency domain [24], although this blurs the ability to analyse/interpret the time-evolution of the spectral response (which becomes well-defined again for longer time scales).