Amplitude dynamics of charge density wave in LaTe$_3$: theoretical description of pump-probe experiments

TL;DR

This paper develops a theoretical model combining Ginzburg-Landau equations and a three-temperature approach to describe the ultrafast dynamics of charge density waves in LaTe$_3$ following photoexcitation, matching experimental observations.

Contribution

It introduces a coupled dynamical model for CDW amplitude evolution that integrates electronic and lattice temperature effects, providing a comprehensive description of pump-probe experiment results.

Findings

Model accurately reproduces experimental CDW amplitude dynamics.

Good agreement with ultrafast electron diffraction data.

Provides insights into the interplay of electronic and lattice responses.

Abstract

We formulate a dynamical model to describe a photo-induced charge density wave (CDW) quench transition and apply it to recent multi-probe experiments on LaTe$_3$ [A. Zong et al., Nat. Phys. 15, 27 (2019)]. Our approach relies on coupled time-dependent Ginzburg-Landau equations tracking two order parameters that represent the modulations of the electronic density and the ionic positions. We aim at describing the amplitude of the order parameters under the assumption that they are homogeneous in space. This description is supplemented by a three-temperature model, which treats separately the electronic temperature, temperature of the lattice phonons with stronger couplings to the electronic subsystem, and temperature of all other phonons. The broad scope of available data for LaTe$_3$ and similar materials as well as the synergy between different time-resolved spectroscopies allow us to extract model parameters. The resulting calculations are in good agreement with ultra-fast electron diffraction experiments, reproducing qualitative and quantitative features of the CDW amplitude evolution during the initial few picoseconds after photoexcitation.

Tables

Table 1. Table 1: Parameters used for numerical simulations of Eqs. ( 12 ) and ( 15 ), and Eqs. ( 16 – 18 ).

TDGL sector			Temperature-evolution sector
Parameter	Value	Physical meaning	Parameter	Value	Physical meaning
${\omega }_0/\left(2\pi \right)$	3.1 THz	Unrenormalized frequency	$c_0$	1.1 mJ/mol$\cdot$K2	Coefficient in the low-temperature
		of the CDW phonon			electronic heat capacity, Eq. (23)
${\tau }_0$	20 fs	Electronic density	$c$	4 mJ/mol$\cdot$K2	Coefficient in the high-temperature
		relaxation time, Eq. (15)			electronic heat capacity, Eq. (24)
$\zeta$	1.1	Critical temperature renormalization	$C_{\mathrm{tot}}$	99.7 J/mol$\cdot$K	Total heat capacity of the lattice
		parameter, Eq. (6)
${\gamma }_y$	$0.04$	Damping parameter for	$\kappa$	0.2	Fraction of “hot phonons” with
		the CDW phonon			respect to all phonons
$T_c$	670 K	Temperature of the CDW	$G_{\mathrm{eL}}$	5.5 J/ps$\cdot$K$\cdot$mol	Energy exchange rate between
		phase transition			electrons and “hot phonons”
			$G_{\mathrm{LL}}$	7.25 J/ps$\cdot$K$\cdot$mol	Energy exchange rate between “hot
					phonons” and the rest of the lattice
			${\tau }_{\mathrm{DW}}$	2.2 ps	Time constant for the second stage
					of the temperature relaxation, Eq. (27)
			$T_{\mathrm{env}}$	300K	Equilibrium temperature of the entire
					system before the laser pulse

Equations

$$\delta {\rho }_{\mathrm{e}}(\mathbf{r})=A\exp (i\mathbf{Q}\cdot \mathbf{r})+c.c.,$$

$$\delta {\mathbf{r}}_n=\left[iu\exp (i\mathbf{Q}\cdot {\mathbf{r}}_n)+c.c.\right]{\hat{e}}_{\mathrm{CDW}}$$

$$\mathcal{F}[A,u]=-a|A{|}^2+\frac{b}{2}|A{|}^4-\eta (A{u}^{\ast }+A^{\ast }u)+\mathcal{K}|u{|}^2,$$

$$a=\alpha (T_0-T),$$

$$T_{\mathrm{c}}=T_0+\frac{{\eta }^2}{\alpha \mathcal{K}}.$$

$$T_0=T_{\mathrm{c}}(1-\zeta ),$$

$$\zeta =\frac{{\eta }^2}{\alpha T_{\mathrm{c}}\mathcal{K}}$$

$$A_{\mathrm{eq}}=\sqrt{\frac{\alpha }{b}(T_{\mathrm{c}}-T)},u_{\mathrm{eq}}=\frac{\eta }{\mathcal{K}}{A}_{\mathrm{eq}}.$$

$$x=\frac{A}{A_{\mathrm{eq}}(T=0)},y=\frac{u}{u_{\mathrm{eq}}(T=0)}.$$

$$x_{\mathrm{eq}}=y_{\mathrm{eq}}=\sqrt{\Theta },$$

$$\Gamma \frac{dA}{dt}=-\frac{\partial \mathcal{F}}{\partial A^{\ast }}=aA-b|A{|}^{2}A+\eta u,$$

$${\tau }_0\frac{dx}{dt}-\Theta x+|x{|}^{2}x+\zeta (x-y)=0,$$

$$\Theta \bar{x}-{\bar{x}}^3+\zeta (y-\bar{x})=0.$$

$$f_u=-\frac{d\mathcal{F}}{d{u}^{\ast }}=-\mathcal{K}u+\eta A.$$

$$\frac{1}{{\omega }_0^2}\frac{d^{2}y}{d{t}^2}+\frac{{\gamma }_y}{{\omega }_0}\frac{dy}{dt}+(y-x)=0,$$

$$C_{\mathrm{e}}(T_{\mathrm{e}})\frac{d{T}_{\mathrm{e}}}{dt}=-G_{\mathrm{eL}}(T_{\mathrm{e}}-T_{L2}),$$

$$C_{L1}\frac{d{T}_{L1}}{dt}=-G_{\mathrm{LL}}(T_{L1}-T_{L2}),$$

$$C_{L2}\frac{d{T}_{L2}}{dt}=-G_{\mathrm{eL}}(T_{L2}-T_{\mathrm{e}})-G_{\mathrm{LL}}(T_{L2}-T_{L1}).$$

$$\hslash {\omega }_{\gamma }\mathcal{V}F=\underset{T_{\mathrm{env}}}{\overset{T_{\mathrm{e}}(0)}{\int }}{C}_{\mathrm{e}}(T)dT.$$

$$\frac{{\omega }_{\mathrm{AM}}}{2\pi }={\nu }_{\mathrm{AM}}=2.2\mathrm{THz},{\gamma }_{\mathrm{AM}}=1.26\times 10^{12}{\mathrm{s}}^{-1}.$$

$${\tau }_0\gtrsim \frac{1}{2\Delta }\approx 1\mathrm{fs},$$

$${\omega }_{\mathrm{AM}}(T=0)={\lambda }_{\mathrm{CDW}}^{1/2}{\omega }_0,$$

$$C_{\mathrm{e}}^0(T_{\mathrm{e}})=c_0{T}_{\mathrm{e}},$$

$$C_{\mathrm{e}}(T_{\mathrm{e}})=\{\begin{array}{ll}{c}_0{T}_{\mathrm{e}},& \text{if }{T}_{\mathrm{e}}<T_{\mathrm{env}},\\ c_0{T}_{\mathrm{env}}+(T_{\mathrm{e}}-T_{\mathrm{env}})c,& \text{if }{T}_{\mathrm{env}}<T_{\mathrm{e}},\end{array}$$

$$C_{L1}=(1-\kappa )C_{\mathrm{tot}},C_{L2}=\kappa C_{\mathrm{tot}},$$

$$C_{\mathrm{e}}(T_{\mathrm{e}})\frac{d{T}_{\mathrm{e}}}{dt}\approx -G_{\mathrm{eL}}{T}_{\mathrm{e}},$$

$${\tau }_{\mathrm{DW}}=\kappa (1-\kappa )C_{\mathrm{tot}}/G_{\mathrm{LL}}.$$

$$T_{\mathrm{e}}({\tau }_{\mathrm{e}})-T_{\mathrm{env}}=\frac{T_{\mathrm{e}}(0)-T_{\mathrm{env}}}{e},e=2.718\dots .$$

$$T_{\mathrm{e}}({\tau }_{\mathrm{c}})=T_{\mathrm{c}},$$

$$1\times 10^{20}{\mathrm{cm}}^{-3}\lesssim F_{\mathrm{c}}<4\times 10^{20}{\mathrm{cm}}^{-3}.$$

$$g(t)=\frac{1}{\sqrt{2\pi }w}\exp \left(-\frac{t^2}{2w^2}\right)$$

Full text

I Introduction

Dynamics of phase transitions associated with spontaneous symmetry breaking remains an interesting subject both theoretically and experimentally. Thanks to the advances in time-resolved pump-probe techniques, it is now possible [​​ 1; 2; 3; 4; 5; 6; 7; 8; 9] to perturb an ordered state and then monitor its fast non-adiabatic recovery. For strong perturbations, one can observe a passage through an ordering transition, register the emergence of ordered phases, and measure time evolution of diverse system parameters with a subpicosecond resolution. The responses of ordered phases, such as superconducting phase [​​ 1; 6; 10; 11], spin-density-wave [​​ 2; 7; 9] and charge-density-wave [​​ 12; 3; 13; 14; 15; 16; 17; 4; 5; 18; 19; 20; 21; 22; 8; 23; 24; 25; 26] phases, have been studied this way.

The focus of the present work is on the non-equilibrium dynamics across a CDW transition. Despite long and thorough scrutiny [​​ 27; 28], the CDW state continues to generate ample amount of research activity motivated by interesting many-body physics (collective transport phenomena [​​ 28; 29; 30; 31], non-mean field critical exponents [​​ 32; 33; 34], exotic metastable ‘hidden’ states [​​ 35; 36; 37; 38]), and large number of experimentally available model systems. In particular, one can mention the “classical” CDW materials, such as [​​ 39; 40; 41; 42] NbSe3, NbSe2, TaS2, blue bronzes K0.30MoO3 and Rb0.30MoO3. CDW phase was also observed and actively investigated in the family of rare-earth tritellurides $R$Te3 [​​ 43; 44; 45; 46; 47; 33; 48; 3; 49; 5; 4; 50; 51; 52; 24].

Recently, we reported [​​ 8] results of an experimental multi-probe study of a photo-induced CDW transition in a member of the rare-earth tritelluride family LaTe3. In these experiments, the post-pump relaxation was monitored with the help of three different time-resolved techniques: ultra-fast electron diffraction (UED), transient reflectivity, and time- and angle-resolved photo-emission spectroscopy (tr-ARPES). These measurements delivered a wealth of complementary information about both electronic and lattice degrees of freedom, and, at the same time, highlighted the need for advancing a quantitative theoretical description of the far-from-equilibrium dynamics in CDW materials and beyond. In the present paper, the experiments of Ref. 8 will serve as the primary testing ground for a rather general theoretical approach.

The CDW order is characterized by the amplitude and the phase of charge modulations. Both can strongly fluctuate in time and space. The amplitude fluctuations are sometimes referred to as “the Higgs modes”, while the phase fluctuations are associated with the appearance of phasons and topological defects, e.g., dislocations. The experiments of Ref. 8 produced evidence that the relaxation of the phase is significantly slower than that of the amplitude. The slowness of the phase relaxation was interpreted as being due to the presence of topological defects. Modeling the latter, however, is rather expensive computationally, because it requires space-resolved simulations. It is, therefore, reasonable to ask first what part of the observed phenomenology can be accounted for by the amplitude relaxation only. This is what we do in the present work. The phase relaxation is to be investigated elsewhere.

In what follows, we develop a theoretical framework for describing the time evolution of space-averaged CDW amplitude in response to a strong quench induced by a femtosecond laser pulse. In the course of the ensuing nonequilibrium evolution, electrons relax much faster than the lattice. Therefore, the modulation of the electronic density and the modulation of the lattice should be treated as two distinct components of the CDW order. We do this using time-dependent Ginzburg-Landau (TDGL) equations [​​ 53; 3; 4; 54; 55; 56; 24; 23; 57] with two order parameters [​​ 54; 55]. We also approximate energy redistribution between different degrees of freedom using the so-called “three-temperature model” [​​ 58; 59; 60; 61; 62; 63], which assigns separate temperatures to (i) electrons, (ii) phonons efficiently coupled to the electronic subsystem, and (iii) all other phonons. Both TDGL equation with two order parameters and the three temperature model were previously considered in the above-mentioned references but not in combination with each other. The former was also tested only in the small-oscillations regime [​​ 54; 55], which characterizes the response of the system to a weak quench.

Fusing together the two-order-parameter TDGL equations and the three-temperature model into a single capable formalism is the main agenda of this paper. Our analysis indicates that the combination of the above two ingredients constitutes a quantitatively adequate and yet efficient theoretical framework for treating strong photo-induced quenches in CDW materials. This framework should also be applicable to other materials, where electronic order is coupled to the lattice.

The paper is organized as follows. In Sec. II we present the theoretical formalism. The values of parameters for the resulting equations are fixed in Sec. III. Numerical simulations are compared with experimental data in Sec. IV. Emergent time scales are summarized in Sec. V. Section VI contains discussion. Finally, the conclusions are presented in Sec. VII. Technically involved derivations are relegated to the Appendices.

II Theoretical formalism

II.1 Preliminary qualitative considerations

In an experiment, a laser pumping pulse initially excites mostly electronic degrees of freedom, while keeping the lattice unaffected. The ensuing internal thermalization of the electronic subsystem is much faster than that of the lattice. It occurs via electron-electron interactions. Given our assumption of a strong laser pulse, the resulting electronic temperature is significantly higher than the initial temperature of the system – possibly higher than the CDW transition temperature. To be specific, let us focus on the latter possibility. In such a case, if the system were completely in equilibrium, then electronic density modulations would quickly disappear. However, since the ions had no time to respond yet, the lattice modulation associated with CDW remains intact, which imposes an external periodic potential on the electronic subsystem. Therefore, once the electronic density relaxes, the electronic order parameter starts tracking the lattice order parameter. During the subsequent evolution, the lattice, on the one hand, experiences a diminished push from electrons to assume a modulated structure; hence the amplitude of the lattice modulation gradually decreases. This happens in an oscillatory way due to the motion of heavy ions near their equilibrium positions. On the other hand, the electronic subsystem, whose heat capacity is much smaller than that of the lattice, rapidly loses energy to the lattice. As a result, the electronic temperature also decreases and eventually falls below the temperature of the CDW phase transition. Once this happens, both the electronic and the lattice orders start recovering, while the electronic temperature continues to decrease until it reaches the temperature of the lattice.

Below we develop a theoretical model capturing the above non-equilibrium evolution. It describes the electronic and the lattice CDW amplitudes via the TDGL formalism. For a non-equilibrium state, we employ the so-called three temperature model, where the electronic and the two lattice subsystems are are assumed to thermalize internally to quasi-equilibrium distributions characterized by different temperatures. We then use rate equations to describe the energy flow between these subsystems.

II.2 TDGL sector

Among theoretical tools [​​ 64; 65; 27; 66; 67; 68; 69; 70; 71; 72; 73] capturing the dynamics of an order parameter, the TDGL equation [​​ 53; 3; 4; 54; 55; 56; 24; 23; 57] is one of the most popular. Despite the known issues with its microscopic justification [​​ 74], TDGL formalism remains in wide use due to its simplicity and intuitive appeal. Below, we introduce the ingredients of this formalism adapted to a setting with two order parameters.

II.2.1 Static Landau functional

The CDW state is characterized by both the modulation of the electronic density

[TABLE]

and the displacements of ions

[TABLE]

from high-symmetry lattice positions ${\bf r}_{n}$. Here $\bf Q$ is the CDW wave vector; $\hat{e}_{\rm CDW}$ is the unit vector along the CDW displacements. Quantities $A$ and $u$ describe complex electronic and lattice order parameters. In Appendix A, we discuss the connection between the lattice order and the diffraction measurements.

To describe equilibrium properties of the CDW system, one can introduce a Landau-type functional [​​ 75], which depends either exclusively on $A$, or exclusively on $u$, because in equilibrium $A$ is proportional to $u$ (for sufficiently small $A$ and $u$). However, to account for the non-equilibrium properties of the CDW during photo-induced transition, we need both $A$ and $u$. We thus consider the Landau functional of the following form [​​ 54; 55]:

[TABLE]

where $a$, $b$, $\eta$ and ${\cal K}$ are the expansion parameters. The first two terms in Eq. (3) are of purely electronic origin. The last term corresponds to the elastic lattice energy, which increases if the ions are shifted from their most symmetric positions. Finally, the term proportional to $\eta$ describes the electron-lattice coupling – often the main driving force behind the CDW transition.

Below we assume that parameters of the Landau functional are temperature independent, except for

[TABLE]

where $\alpha$ is a positive proportionality coefficient, $T$ is the temperature, and $T_{0}$ is the “bare” transition temperature for a hypothetical situation of vanishing electron-lattice interaction. (In principle, the parameter $T_{0}$ is not a physical temperature but simply a parameter characterizing the structure of the Landau functional; hence, it can be negative.) Due to finite coupling between $A$ and $u$, the actual transition into the ordered phase occurs at the critical temperature

[TABLE]

For LaTe3, $T_{\rm c}\approx 670$K [​​ 8]. As for $T_{0}$, it can be estimated as

[TABLE]

where

[TABLE]

characterizes the strength of the electron-phonon interaction. For the parameters chosen in Sec. III to represent LaTe3, we obtained $T_{0}=-67$ K.

Minimizing the functional (3) for $T<T_{\rm c}$, one finds equilibrium values of the order parameters:

[TABLE]

For our calculations, it is convenient to work with the dimensionless quantities:

[TABLE]

For $T\leq T_{\rm c}$, equilibrium values of $x$ and $y$ are

[TABLE]

where ${\Theta}=(T_{\rm c}-T)/T_{\rm c}$.

II.2.2 Time-dependent equations

The next step is to generalize the static Landau theory to non-equilibrium situations. We will describe the dynamics of the electronic degree of freedom $A$ as

[TABLE]

where $\Gamma$ is a damping parameter. In dimensionless variables (9), Eq. (11) reads

[TABLE]

where $\tau_{0}=\frac{\Gamma}{\alpha T_{\rm c}}$ is the electronic density relaxation time.

From the viewpoint of the true microscopic kinetics, Eq. (12) is a crude approximation. However, for the purposes of the present work, this approximation should be sufficient given that we are mainly interested in the dynamics of variable $y$, which unfolds on the time scale much longer than $\tau_{0}$. This longer time scale is associated with the motion of ions, which are much heavier and, thus, much slower than electrons. It is only important for us that, on fast time scale of $\tau_{0}$, variable $x$ relaxes to the “instantaneous equilibrium” value $\bar{x}$ given by the real-valued root of the equation

[TABLE]

To model the evolution of $u$, we keep in mind that $u$ is associated with displacements of heavy ions, which can be viewed, approximately, as classical objects. The forces acting on the ions are

[TABLE]

The term proportional to ${\cal K}$ in Eq. (14) describes elastic force that pulls ions back to their high-symmetry positions. The term proportional to $\eta$ originates from the interaction with the modulated electron density (1). By adding damping, we arrive at a classical equation of motion for $u$. In the rescaled variables (9), it reads

[TABLE]

where $\omega_{0}=\sqrt{{\cal K}/m}$ is unrenormalized frequency of a phonon mode actively involved in the CDW formation (we call this mode “the CDW phonon”); $m$ is the ion mass parameter; $\gamma_{y}$ describes damping. Equations (12) and (15) constitute the desired TDGL equations in a dimensionless form.

Previously TDGL formalism [54; 55] with two order parameters was applied to describe small deviations from equilibrium in a CDW material. However, in the present work we deal with a far-from-equilibrium response, which requires an additional ingredient in the theory. This additional ingredient is the three-temperature model introduced in the next subsection.

II.3 Three-temperature model

Temperature $T$ that appears in Eq. (4) is the electronic temperature which, from now on, we denote as $T_{\rm e}$. We need this change of notation to distinguish the electronic temperature from the temperature (or temperatures) of the lattice. The TDGL equation (12) thus depends on $T_{\rm e}$.

Let us expand on the qualitative scenario given in Section II.1 using the numbers specific for the experiment of Ref. 8. According to that scenario, all photons from the pumping laser pulse are mostly absorbed by the electronic subsystem. Right after the laser pulse, the electronic subsystem is far from equilibrium. However, it quickly reaches quasistationary thermal state with temperature $T_{\rm e}(0)\sim 1000$ K, which significantly exceeds the initial lattice temperature of about 300 K. Lattice phonons, whose heat capacity is much larger than that of the electrons, then act as a heat sink gradually absorbing energy of the hot electrons. Let us emphasize that, during this process, it can happen that phonon distribution function becomes highly non-thermal.

The simplest approach to capture the above dynamics would be to introduce a two-temperature model [​​ 76; 77; 78] tracking the electronic and the lattice temperatures. However, given all available experimental data – including (i) heat capacity measurements [​​ 46] in LaTe3, (ii) UED data (in particular, the time dependence of the crystal Bragg peaks intensity [​​ 8]), (iii) short-time transient-reflectivity dynamics [​​ 8], and (iv) tr-ARPES data [​​ 8], which allows to estimate electronic heat capacity (see Sec. III below) – we were unable to adequately reproduce the observed behavior using the two-temperature model.

We attribute the above discrepancy to an intrinsic limitation of the two-temperature model. Namely, such a model describes the entire lattice by a single temperature, while, in reality, not all phonons are thermalized simultaneously. Instead, we expect that the electronic energy is preferentially absorbed by a smaller subset of the phonon ensemble. Such a selective coupling was discussed for rare-earth tritellurides in Refs. 79; 80 and, in a broader context, in Refs. 61; 81; 82; 63. In particular, the authors of Ref. 60 argued that layered structure makes preferential electron coupling to a subset of phonon modes more likely. For LaTe3, we can identify two factors that increase the selectivity of the electron-phonon coupling. First, electronic states forming the Fermi surface are located in the tellurium-only layers DiMasi et al. (1995) and, hence, are weakly coupled to the lattice degrees of freedom localized in the $R$Te layers. Second, at sufficiently low excitation energies of the electronic quasiparticles, the momentum conservation law implies that the probability of emission (or absorption) of a phonon with momentum ${\bf q}$ is proportional to $\sum_{\bf k}\delta(\varepsilon_{\bf k}-\varepsilon_{\rm F})\delta(\varepsilon_{{\bf k}+{\bf q}}-\varepsilon_{\rm F})$, where $\varepsilon_{\rm F}$ is the Fermi energy, and $\varepsilon_{\bf k}$ is the electron dispersion. In the rare-earth tritellurides, due to proximity to the Fermi surface nesting condition, the latter sum is a strongly non-uniform function of the momentum ${\bf q}$, see Eq. (2) and Fig. 6a of Ref. 83. This non-uniformity imposes additional selection criteria on the phonon modes participating in the electron scattering.

In order to account for the selectivity of the electron-phonon coupling, we split the lattice phonons into two subgroups: (i) the “hot phonons”, the ones strongly coupled to the electrons, and (ii) “cold phonons”, i.e. the rest of the lattice modes. The hot phonons are to be characterized by temperature $T_{\rm L2}$, while cold phonons by temperature $T_{\rm L1}$. To simplify the model, we neglect the direct energy transfer between the electronic subsystem and the cold phonons.

The post-pulse temperature dynamics is then described by the three-temperature model [​​ 58; 59]:

[TABLE]

Here $C_{\rm L1}$ and $C_{\rm L2}$ are the heat capacities of the cold phonons and of the hot phonons, respectively; $C_{\rm e}(T_{\rm e})$ is the temperature-dependent electronic heat capacity whose functional form is specified in subsection III.2. Parameters $G_{\rm eL}$ and $G_{\rm LL}$ describe energy exchange rates.

The initial conditions for Eqs. (16-18) are to be chosen as follows. For the lattice, $T_{\rm L1}(0)=T_{\rm L2}(0)=T_{\rm env}$, where $T_{\rm env}$, is the pre-pulse equilibrium temperature of all three subsystems. For electrons, $T_{\rm e}(0)$ is defined as the temperature right after the laser pulse and the ensuing fast electronic self-thermalization. The value $T_{\rm e}(0)$ is, therefore, a function of absorbed electromagnetic energy per mole, which is, in turn, proportional to (i) the photoexcitation density $F$ (the number of absorbed photons per unit volume), (ii) the pump photon energy $\hbar\omega_{\gamma}$, and (iii) the molar volume ${\cal V}$ of LaTe3. The energy balance condition then gives the equation determining $T_{\rm e}(0)$:

[TABLE]

III Choice of parameters for $\text{LaTe}_{3}$

Overall, our formalism includes five equations: (12) and (15) for the TDGL sector, and (16–18) for the temperature-evolution sector. To perform simulations, one needs to select specific values for the model parameters: $\omega_{0}$, $\tau_{0}$, $\zeta$, and $\gamma_{y}$ in the TDGL sector; and $C_{\rm L1}$, $C_{\rm L2}$, $G_{\rm eL}$, and $G_{\rm LL}$ in the temperature-evolution sector. One also needs a concrete functional form of the temperature dependence of the electronic heat capacity $C_{\rm e}(T_{\rm e})$.

Since the total number of free parameters is large, extracting their values through an indiscriminate fitting might produce misleading results. To circumvent this issue, we split the whole task into several steps to be presented in the following subsections. In each step, only a small number of the unknowns are fixed. This approach relies on the availability of a broad array of experimental results for the rare-earth tritellurides $R$Te3, in particular LaTe3 [​​ 8]. For convenience of the readers, the final values of all parameters are collected in Table 1.

III.1 TDGL parameters

Here we fix the TDGL parameters by matching the frequency $\omega_{\rm AM}$ and the damping constant $\gamma_{\rm AM}$ of the CDW amplitude mode (AM) obtained theoretically with the values measured experimentally. Theoretically, we apply small-oscillations formalism [​​ 54; 55] to the TDGL equations (12) and (15), see Appendix B for detailed derivations. On the experimental side, the phonon spectrum in $R$Te3 and its temperature dependence were reported in many works [​​ 47; 49; 79; 80; 8].

From the average frequency and the frequency width of the measured transient reflectivity oscillations, we know that [​​ 8], at $T=300$ K,

[TABLE]

We use the above values to define a complex parameter $\lambda=-\gamma_{\rm AM}+i\omega_{\rm AM}$, which is then substituted into Eqs. (59, 58) describing the small-oscillations eigenvalue problem. This leads to two constraints on parameters $\omega_{0}$, $\tau_{0}$, $\zeta$, and $\gamma_{y}$ expressed by Eqs. (62) and (63). Thereby, we reduce the total number of adjustable TDGL parameters from four to two. For our analysis, it is convenient to treat Eqs. (62) and (63) as definitions of the implicit functions $\gamma_{y}(\omega_{0},\tau_{0})$ and $\zeta(\omega_{0},\tau_{0})$, which specify the dependence of $\zeta$ and $\gamma_{y}$ on $\omega_{0}$ and $\tau_{0}$. Thus, once $\omega_{0}$ and $\tau_{0}$ are fixed, the TDGL sector contains no unknown parameters.

Estimating $\omega_{0}$ and $\tau_{0}$, one must be mindful of several relevant theoretical and experimental restrictions. The first of them is the physical requirement $\gamma_{y}(\omega_{0},\tau_{0})\geq 0$. It limits the allowed space for $\omega_{0}$ and $\tau_{0}$ to the region left and below the red curve in Fig. 1.

The next restriction is related to whether the AM frequency $\omega_{\rm AM}$ softens to zero close to the CDW transition temperature $T_{\rm c}$ or not. In Refs. 54; 55, the former regime is called “adiabatic” and the latter “non-adiabatic” (see examples in Appendix B). Although the behavior of the AM in LaTe3 near $T_{\rm c}$ is not accessible experimentally, we rely on the reported universality of the AM characteristics for several members of the $R$Te3 family [​​ 48]. (The most noticeable aspect of this universality is the same low-temperature value of the AM frequency $\nu_{\rm AM}\approx 2.2$ THz.) Specifically, experiments suggest that the AM in TbTe3 [see Fig. 3(b,c) in Ref. 79] and DyTe3 [see Fig. 8(a,c) in Ref. 80] softens to zero close to the transition temperature. Thus, we assume the adiabatic regime for LaTe3 as well. As shown in Fig. 1, such an assumption further confines $\omega_{0}$ and $\tau_{0}$ to the region above the dashed line. Together with the previous constraint, this implies that $\tau_{0}\lesssim 30$ fs.

Now we note that the Heisenberg uncertainty principle suggests that

[TABLE]

where $2\Delta\approx 700$ meV is the CDW gap at 300 K [​​ 51]. Given the above constraints, we assign $\tau_{0}=20$ fs. The dynamics of the CDW order parameters obtained from our numerical simulations is not very sensitive to the specific choice of $\tau_{0}$ as long as we are interested in time scales much longer than $\tau_{0}$.

Finally, to set $\omega_{0}$, we use the relation

[TABLE]

where $\lambda_{\rm CDW}$ is the electron-phonon coupling constant [​​ 84] responsible for the CDW instability. It is assumed both for the rare-earth tritellurides [​49] and for a broader class of materials [see, e.g., discussion after Eq. (5) in Ref. 27, and Table 3.1 of Ref. 84] and that $\lambda_{\rm CDW}\approx 0.5$. When this value together with $\omega_{\rm AM}$ given by Eq. (20) is substituted into Eq. (22), we obtain an estimate $\omega_{0}/(2\pi)\approx 3.1$ THz. Here we assumed that $\omega_{\rm AM}$ is the same at $T=0$ and at $T=T_{\rm env}$. Indeed, since $T_{\rm env}$ is significantly below $T_{c}$, it is permissible to treat $\omega_{\rm AM}$ as being temperature-independent in the range $T\leq T_{\rm env}$.

Once $\tau_{0}$ and $\omega_{0}$ are determined, both $\gamma_{y}=\gamma_{y}(\omega_{0},\tau_{0})$ and $\zeta=\zeta(\omega_{0},\tau_{0})$ are obtained. The final values are summarized in Table 1.

III.2 Three-temperature model parameters

In the context of LaTe3 experiments of Ref. 8, we assume that the equilibrium pre-pulse temperature of all three subsystems is $T_{\rm env}=300$ K. For each value of the photoexcitation density $F$, the electronic temperature right after the initial self-thermalization of the electronic subsystem $T_{\rm e}(0)$ is to be calculated using Eq. (19) with $\hbar\omega_{\gamma}=1.19\,{\rm eV}=1.9\times 10^{-19}\,{\rm J}$ and ${\cal V}=76.8\,{\rm cm^{3}/mol}$.

Let us now turn to the temperature-dependence of $C_{\rm e}(T_{\rm e})$. At sufficiently low temperatures, electronic heat capacity is a linear function of temperature

[TABLE]

where [​​ 46] $c_{0}=1.1$ mJ/mol K2 for LaTe3. In the proximity to and above $T_{\rm c}=670$ K, we do not expect Eq. (23) to remain valid. It would imply that the electronic temperature following the maximum intensity laser pulse reaches the value $\approx 4000$ K, while the analysis of our tr-ARPES data presented in Appendix C.2 reveals that $T_{\rm e}\lesssim 2000$ K. The deviation from the linear temperature dependence (23) is also expected on the basis of purely theoretical reasoning outlined in Appendix C.1. Following that reasoning, we approximate $C_{\rm e}(T_{\rm e})$ by a piecewise linear ansatz:

[TABLE]

where $c=4$ mJ$/$mol K2 is a parameter extracted from tr-ARPES experiments in Appendix C.1 .

Next we turn to lattice heat capacities $C_{\rm L1}$ and $C_{\rm L2}$. Since these parameters are associated with two complementary groups of phonons, we express

[TABLE]

where $C_{\rm tot}$ is the total heat capacity of the lattice, and $\kappa$ is the coefficient determining the fraction of the phonon modes contributing to $C_{\rm L2}$. We approximate $C_{\rm tot}$ by the Dulong-Petit value $99.7$ J/(mol K) for LaTe3, which is permissible in the temperature range of interest and consistent with experiment [​​ 46]. The value of $\kappa$ is fixed to be equal to $0.2$ in subsection IV.2 on the basis of our model fitting to the UED Bragg peak intensities.

A typical post-pulse time evolution of the three temperatures calculated on basis of Eqs. (16-18) is shown in Fig. 2. Here we assume that the initial rise of the electronic temperature occurs on a very short time scale, which we approximate as instantaneous. The remaining evolution can be divided into two stages. During the first stage, $T_{\rm e}$ relaxes to $T_{\rm L2}$ on the time scale of 1–2 ps. The second stage unfolds for $t\gtrsim 1.5$ ps, where the common temperature of the electrons and the hot phonons ($T_{\rm e}\approx T_{\rm L2}$) approaches $T_{\rm L1}$.

For sufficiently strong laser pulses, such that $T_{\rm e}(0)\gg T_{\rm L1,L2}$, the first stage can be accurately described by the approximate equation

[TABLE]

which is governed by a single parameter $G_{\rm eL}$. Its value can be estimated by assuming that, at high excitation densities, the initial decay of the transient reflectivity, measured in Ref. 8, is controlled by $T_{\rm e}(t)$. This way, we obtain $G_{\rm eL}=5.5$ J/(ps$\cdot$K$\cdot$mol), see subsection IV.1 for further details.

During the second stage, the temperature relaxation process is exponential, characterized by the time constant

[TABLE]

Here, following the notation of Ref. 8, we use the subscript ‘DW’, which stands for ‘Debye-Waller’, because the above time constant controls the evolution of the Bragg peak intensity in the late-time regime. Expression (27) can be derived with the help of Eqs. (17) and (18) in the limit $T_{\rm e}=T_{\rm L2}$ (corresponding to $t>1.5$ ps in Fig. 2). From the measured relaxation of the Bragg peak intensity [​​ 8], we have $\tau_{\rm DW}=2.2$ ps. Thereby, Eq. (27) defines $G_{\rm LL}=7.25$ J/(ps$\cdot$K$\cdot$mol).

IV Comparison of the experiments with the numerical simulations

In this section, we present the results of the simulations and compare them with experiments of Ref. 8. The simulation parameters are given in Table 1.

IV.1 Comparison to the short-time transient reflectivity

measurements

Laser pulse initially excites electronic degrees of freedom, which, in turn, excite the lattice. Both the electrons and the lattice contribute to the change in the transient reflectivity signal. We expect that the monotonically decaying part of the transient reflectivity (measured transient reflectivity with the oscillating contribution from the amplitude mode subtracted [​​ 8]) tracks the dynamics of the electronic temperature $T_{\rm e}(t)$.

The computed time dependence $T_{\rm e}(t)$ is shown in Fig. 3 for different excitation densities $F$. The crossover between the first rapid stage to the second slow stage is clearly seen. For longer times $T_{\rm e}(t)$ approaches some photoexcitation-density-dependent base temperature, which only slightly exceeds $T_{\rm env}$ for all excitation densities used.

To compare the computed temperature evolution with experiment, we introduce time $\tau_{\rm e}$ by condition

[TABLE]

It characterizes the time scale of the electronic temperature cooling down during the first rapid stage. In Fig. 4, we compare $\tau_{\rm e}$ with the relaxation time $\tau_{\rm R}$, extracted from the transient reflectivity experiment [​​ 8]. Given the simplicity of our model, the agreement between the theory and the experiment is reasonable. It was attained by adjusting the parameter $G_{\rm eL}$ in Eq. (26), while other parameters affecting the latter equation were fixed as described in subsection III.2.

IV.2 Dynamics of the order parameters

IV.2.1 Melting of the CDW order

In Fig. 5, the plots of $x^{2}(t)$ and $y^{2}(t)$ illustrate the typical dynamics of the electronic and the lattice CDW order parameters after the arrival of a laser pulse. For low excitation densities, such as that of Fig. 5a, the laser pulse does not completely destroys the CDW order – it only excites damped AM oscillations around the equilibrium values of the order parameters. For stronger pulses, as in panels (b–d), both $x(t)$ and $y(t)$ cross zero, which, despite the lack of equilibrium, indicates the proximity to the melting of the CDW order.

To investigate the onset of the CDW melting, it is useful to look back at the corresponding temperature evolution $T_{\rm e}(t)$ shown in Fig. 3, where we see that $T_{\rm e}(t)$ stays above $T_{\rm c}$ for a finite time $\tau_{\rm c}(F)$. This time is defined by the condition

[TABLE]

which is indicated in Fig. 3 by the dash-dotted horizontal line. Once $F$ is larger than a certain threshold, $\tau_{\rm c}$ grows monotonically with $F$, as shown in Fig. 4.

In Fig. 5(b–d), the time $t=\tau_{\rm c}$ is indicated by vertical dashed lines. We see that $\tau_{\rm c}$ indeed gives the correct estimate of the time when the CDW order ceases to decay and starts to recover. We also note that for panel (d), where $\tau_{\rm c}\gtrsim 1$ ps, the order parameters demonstrate multiple passages through zero with decreasing amplitude. After the oscillations fade, the order parameters remain suppressed for about 0.5 ps.

In general, the notion of melting in the course of a non-equilibrium evolution is not sharply defined. Here we adopt the criterion that the CDW order undergoes melting when the electronic and the lattice order parameters do not simply cross zero but rather approach zero in a damped oscillatory (or non-oscillatory) fashion. From such a perspective, Fig. 5d represents the melting behavior, while Fig. 5b does not, and Fig. 5c is the border case. On the basis of the above analysis, we conclude that the critical excitation density $F_{\rm c}$, defined as the lowest border for melting, satisfies the following inequality

[TABLE]

Experimentally [​​ 8] $F_{\rm c}\sim 2.0\times 10^{20}\,{\rm cm}^{-3}$, in agreement with the above constraints.

All plots in Fig. 5 exhibit prominent oscillatory behavior of the order parameters. At lower photoexcitation densities, such as in Fig. 5a, the oscillations are clearly related to the appearance [​​11] of the AM observed in the transient reflectivity experiments. However, the experiment indicates significant reduction of the oscillation amplitude for $F\gtrsim 2\times 10^{20}$ cm*-3*. This discrepancy can be attributed to our assumption that the order parameters are homogeneous in space, while in the real system, the spatial configuration of the order parameters following the melting and the subsequent re-emergence of the CDW is likely strongly inhomogeneous due to the appearance of topological defects in the order parameter texture. As a result of this inhomogeneity, the system has relatively small coherent CDW domains of varying size with different size-dependent frequencies $\omega_{0}$, which, in turn, leads to the strong dephasing of the oscillations, once the signal is averaged over the entire sample.

Further analyzing oscillations in Fig. 5d, we observe that the frequency of transient oscillations for $t\lesssim 1.5$ ps is twice the AM frequency $\omega_{\rm AM}$. Such a doubling occurs because of the interplay of two factors: (i) in Fig. 5 we plot $x^{2}(t)$ and $y^{2}(t)$ instead of $x(t)$ and $y(t)$ and (ii) the order parameters oscillate near $x=y=0$. The experiments of Refs. 85; 57 indicate that such a frequency doubling may, actually, occur in real systems.

As for the UED experiments of Ref. 8, they have insufficient time resolution to detect the order parameters oscillations. To represent the experimental observations, we use the quantity $(y^{2}*g)(t)$, where the asterisk denotes time convolution, and

[TABLE]

is the Gaussian filter, with parameter $w$ representing the time resolution of the experiment. For the UED experiment of Ref. 8, $w=0.42$ ps. The results of the convolution are shown in Fig. 6 for different $F$’s. We can see that the oscillations present in all panels of Fig. 5 were smeared out by the filter.

IV.2.2 Two kinds of diffraction peaks in UED experiments

The UED experiments [​​ 8] observe two kinds of diffraction peaks associated either with the underlying crystal structure of LaTe3 or with the CDW order, see Fig. 7. The measurements of Ref. 8 were done in the higher-order Brillouin zones, which implies that the measured intensities of the CDW peaks are determined [​​ 86] by the lattice CDW order $y$. Fundamentally, the integrated intensity of a CDW peak is proportional to $y^{2}$. Therefore, at first sight, the direct way to test our modeling is to compare the calculated $y^{2}(t)$ with the time evolution of the integrated CDW peak intensity measured in the UED experiment.

However, the problem here is that the experiment 8 measured diffraction intensities only for a two- dimensional slice $(k_{x},k_{z})$ of a three-dimensional reciprocal space $(k_{x},k_{y},k_{z})$, see Fig. 7. In other words, Ref. 8 does not contain direct experimental information about the three-dimensional integrated intensity of the CDW peaks.

At the same time, since, in the experiment [​​ 8], crystal Bragg peaks are resolution-limited, the two-dimensional integrals over the measured sections of these peaks are proportional to the respective integrals over the full three-dimensional reciprocal space. One can now take into account the sum rule implying that the emergence of the CDW order leads to the intensity transfer from the Bragg peaks to the CDW peaks, see Appendix A. Therefore, in our case, the most direct way to extract the value of $y^{2}(t)$ from experiment is to examine the integrated intensity by the Bragg peaks: the suppression of the CDW order leads to increase of the Bragg peak intensity. This relation is quantified in the next subsection.

IV.2.3 Time evolution of the Bragg peaks for the underlying crystal

lattice from UED experiments

For the integrated Bragg peak intensity $I$ in the presence of the CDW order, we use the expression obtained in Ref. 87:

[TABLE]

where $J_{0}$ is the zeroth order Bessel function, and $p$ is a constant. Parameter $W$ accounts for the Debye-Waller suppression of the intensity due to thermal fluctuations of the phonons.

We expect that $W\propto T_{\rm L1}$, and, therefore, write

[TABLE]

where $P=p^{2}/2$, and $S$ is a constant.

As in subsection IV.2.1, to account for the finite experimental temporal resolution, we convolute the rhs of Eq. (33) with the Gaussian filter (31). The final expression used to mimic the actual Bragg peaks dynamics reads

[TABLE]

Here $t_{0}$ is an adjustable parameter shifting the origin of the time axis. This shift is another consequence of the limited time resolution of the experiment [​88].

Function $I(t)$, numerically evaluated with the help of Eq. (34), together with the UED data points, are plotted in Fig. 8. All eight plots in the latter figure were obtained by fitting the experimental points using four adjustable parameters: $\kappa=0.2$, $P=0.1$, $S=3\times 10^{-3}$, and $t_{0}=0.43$ ps. The above value of parameter $\kappa$ was used for all simulations presented in this paper.

Overall, the agreement between the fits and the experiment in Fig. 8 is rather good. For higher excitation densities $F$, the small discrepancy might be due to the fact that approximation (24) for the electronic heat capacity is less accurate at higher temperatures.

IV.2.4 Time evolution of the CDW peak from UED experiments

We now turn to the discussion of the CDW peak. As anticipated in subsection IV.2.2, a straightforward attempt to approximate the experimentally measured $(k_{x},k_{z})$-integrated intensity of a CDW peak by a theoretically computed $y^{2}(t)$ [with the appropriate convolution and time shift, as in Eq. (34)] reveals large discrepancy. Indeed, as one can see from Fig. 9, the calculated curves (red dashed lines) lie significantly higher than the experimental UED data. Below, we test the proposition formulated in Ref. 8 that the deviations between the data and the simulations originate from the fact that the UED measurements access only a two-dimensional $(k_{x},k_{z})$-slice of the three-dimensional $(k_{x},k_{y},k_{z})$-space.

We assume that the intensity of the CDW peak in the reciprocal space can be reasonably approximated by a factorized function

[TABLE]

where $s_{\sigma}(k_{\sigma})$ are the peak shape functions for the respective $k$-space directions, with index $\sigma$ taking values $x$, $y$, or $z$. Following the convention of Ref. 8, axes $x$ and $z$ denote the directions parallel to Te2 planes, while axis $y$ is perpendicular to these planes. (The notations $x$ and $y$ for the spatial axes appear only as subscripts, and are not to be confused with the variables $x$ and $y$ defined by Eq. (9) that represent the electronic and the lattice order parameters.) Functions $s_{\sigma}(k_{\sigma})$ are normalized by the condition

[TABLE]

We also assume that these functions are non-negative and bell-shaped. (In experiment, these functions are fitted by Lorentzians due to the intrinsic profile of the electronic beam.) The two-dimensional integral of the CDW peak reported in Ref. 8 can be written as

[TABLE]

Now we observe that, for common bell-shaped functions, such as Lorentzian or Gaussian, $s_{\sigma}(0)\propto\xi_{\sigma}$, where $\xi_{\sigma}$ is the real-space correlation length in the respective direction. Using this observation and adding explicit time dependencies of the parameters involved, we arrive at the expression

[TABLE]

In this work, we assume that the order parameter is homogeneous is space. Thus, we cannot obtain theoretically $\xi_{y}(t)$. However, we can estimate it on the basis of the assumption that all correlation lengths are determined by the same mechanism and hence are proportional to each other, i.e.

[TABLE]

where $\xi_{\rm exp}(t)$ is the experimentally measured correlation length in the $x$- and $z$-directions, obtained as the inverse of the FWHM of the CDW peak after instrumental resolution is taken into consideration, see Eq. (S4) of the Supplementary Information to Ref. 8.

Finally, taking into account the finite experimental time resolution, as in Eq. (34), the measured integrated intensity can be approximated as

$I^{\rm CDW}_{\rm 2D,exp}(t)=(y^{2}\!*g)(t-t_{0})\,\,\xi_{\rm exp}(t)$. To facilitate the comparison with experiments, we re-express this relation in the following manner

[TABLE]

where the argument $(0^{{}^{-}})$ implies the pre-pulse values of the respective parameters.

In Fig. 9, we test the relation (40) by substituting there the theoretically calculated $y^{2}(t)$. The agreement between the direct measurement and the prediction of Eq. (40) is rather encouraging. This is another consistency check of our modeling.

V Overview of time scales

Here we would like to bring together various aspects of our simulations by attaching concrete time scales to the general non-equilibrium scenario described in Section II.1. In order to be specific, we choose the case of the photoexcitation density $F=3$, which is above the transient-melting threshold $F_{c}=2$ determined in Sec. IV.2.1. (Here and below, the units for $F$ are $10^{20}\text{cm}^{-3}$.) Different stages of the non-equilibrium evolution together with the relevant characteristic times are summarized in Fig. 10.

The fastest time appearing in our description is $\tau_{0}=20$ fs. It enters Eq. (12), and characterizes the relaxation time of the electronic density [​89]. Among time scales of the lattice dynamics, the shortest one is the period of the amplitude mode $2\pi/\omega_{\rm AM}\sim 0.5$ ps. Next is the relaxation time of the CDW phonon mode $1/\gamma_{y}\omega_{0}\sim 1$ ps. Another relevant time having the value of approximately 1 ps is $\tau_{\text{e}}$, which describes the convergence of the electronic temperature $T_{\rm e}$ and the temperature of the hot phonons subsystem $T_{\rm L2}$ (see Fig. 4).

For the photoexcitation density $F=3$, the electronic temperature $T_{\rm e}$ initially jumps to a value significantly higher than $T_{\text{c}}$ and then remains above $T_{\rm c}$ for a time longer than 1 ps. This leaves enough time for the lattice CDW order to relax to zero, which implies complete melting. By the time of about 2 ps, two changes occur: (i) Electronic temperature $T_{\rm e}$ drops below $T_{\rm c}$, which leads to the reappearance of both the electronic and the lattice CDW orders. (ii) Simultaneously, the character of the temperature relaxation changes — the electronic temperature $T_{\rm e}$ and the hot phonons temperature $T_{\rm L2}$, after having approached each other, start decreasing towards the temperature of the rest of the lattice $T_{\rm L1}$ with characteristic time $\tau_{\rm DW}\approx 2$ ps. The recovery of the electronic and the lattice CDW amplitudes becomes eventually completed by the time of about 6 ps. (Due to small overall heating, the CDW amplitude recovers to a value, which is slightly smaller than the pre-pulse one.)

Finally, as indicated in Fig. 10, the longest time scale is supposed to be the one associated with the relaxation of the phase of the CDW order parameter. Since our formalism deals exclusively with the dynamics of the order parameter amplitude, and not the phase, these slow relaxation processes are not covered by our simulations.

VI Discussion

In this paper, we have shown that the combination of the three temperature model and the TDGL equations with two order parameters constitutes a quantitatively adequate description of the amplitude response of CDW materials to a strong femtosecond laser pulse. The successful description of the time evolution of the crystal Bragg peaks intensities shown in Fig. 8 is the testament of the model’s predictive power. Indeed, we were able to simulate a family of eight different time-dependent intensities corresponding to different photoexcitation densities using only one physically important parameter $\kappa$ and three “technical” parameters $P$, $S$, and $t_{0}$. The values of other parameters required by the model were obtained from independent experiments on CDWs in LaTe3 and in other rare-earth tritellurides. We thus demonstrated how to use the available experimental knowledge to systematically extract the parameters for far-from-equilibrium CDW simulations.

We further note that our simulations predict that electronic and lattice modulations are supposed to exhibit an oscillatory dynamics in the vicinity of the CDW melting transition with a frequency which is two times larger than that of the weakly perturbed CDW (see Fig. 5 and the discussion in Sec. IV.2.1). This dynamics would not be directly visible in experiments of Ref. 8 due to finite time resolution but, otherwise, is consistent with the experiments [​​ 85; 57] with the blue bronze K0.3MoO3 and a perovskite-type manganite Pr0.5Ca0.5MnO3. In a related development, the application of the present simulation framework in Ref. 23 pointed at the existence of a dynamical slowing down regime near the CDW melting transition, where the order parameters become “caught” near the metastable maximum of the free energy (3).

In terms of advancing the general knowledge about the far-from-equilibrium CDW dynamics, our simulations shed light on the transient CDW response that is difficult to access experimentally. In particular, the necessity to use two different lattice temperatures — one for hot phonons and the other for the rest of the lattice — confirms previous conjectures [​​79; 80; 61; 81; 82; 63] that the energy transfer from the photo-excited electrons to the phonon bath occurs unevenly among different phonon modes. Overall, the developed theoretical framework should be applicable to other CDW materials and other experimental settings, such as three-pulse experiments of Ref. 3.

Although our treatment only deals with the amplitudes of the order parameters, an important outcome of this work is that, for experimental quantities affected by both the amplitude and the phase relaxation of the CDW order, the lack of the theoretical information about the phase dynamics may be compensated, at least partially, using phenomenology-based approach. An example here is the comparison between the measured 2D-integrated UED intensity of the CDW peak presented in Fig. 9 and the computed intensity. We interpret the difference between the two as being caused by the phase relaxation of the CDW order, which is slower than the amplitude relaxation due to the possible presence of topological defects — in agreement with the analysis of Ref. 8. The consistency of this interpretation is further demonstrated in Fig. 9 by correcting the amplitude-dependent intensity with a factor determined from the experimental knowledge of the peak width, which is, in turn, determined by the phase fluctuations of the order parameter.

As far as the weak points of our modeling are concerned, one of them is its mean-field character. In reality, many CDW systems, including LaTe3, demonstrate pronounced non-mean-field properties near the phase transition [​​ 32; 33; 34]. Fluctuations relative to the mean-field state can be split into two groups: those of the order parameters amplitudes, and those of the phase. The fluctuations of the amplitudes are less worrisome: we expect them to be effectively included into renormalization of the model’s coefficients by replacing some (unknown) “bare” values with the “effective” (observable) values. At the same time, phase fluctuations of the CDW order parameter are of greater concern. They remain largely unaccounted in the present work, which focuses on the CDW amplitudes.

Another limitation of proposed approach is related to the fact that the three-temperature model greatly oversimplifies the kinetic processes in the studied system. However, our success in reproducing the experiments in general and the transient reflectivity experiments in particular suggests that this approximation captures essential physics of the system.

VII Conclusions

We developed a theoretical framework to describe the dynamics of the CDW amplitude after an intense laser pulse. The framework consists of (i) the time-dependent Ginzburg-Landau equations for the electron and lattice CDW amplitudes and (ii) the three-temperature model. We tested the resulting description by comparing the simulations with the available experimental data. The agreement is good, suggesting that the proposed framework can be applied to a broader class of non-equilibrium settings.

VIII Acknowledgements

We thank T. Rohwer, C. Lee, E. J. Sie, E. Baldini, B. Freelon, and H. Zhou for the help in building and acquiring data from the ultra-fast electron diffraction and time-resolved ARPES setups. We acknowledge high-quality samples prepared by J. Straquadine, P. Walmsley, I. R. Fisher, Y.-Q. Bie, and P. Jarillo-Herrero. All authors acknowledge the support of the Skoltech NGP Program (Skoltech-MIT joint project) (theory). N.G., A.Z., and A.K. also acknowledge the support from the Gordon and Betty Moore Foundations EPiQS Initiative grant GBMF4540 (data analysis).

Appendix A Bragg and CDW diffraction peaks

The purpose of this Appendix is to illustrate how the presence of the CDW modifies electron diffraction peaks. We then discuss the role of fluctuations of the phase of the order parameter. The discussion here is a simplified version of a more general treatment of Ref. 87.

We express modulation of the lattice site positions as

[TABLE]

where $\bf u$ and $\bf Q$ are the amplitude and the wave vector of the modulation; ${\bf r}_{n}$ are high-symmetry lattice points. Assuming that the amplitude $|{\bf u}|$ is much smaller than the lattice spacing, we write the density as

[TABLE]

where $\rho_{0}({\bf r})=\sum_{{\bf r}_{n}}\delta({\bf r}-{\bf r}_{n})$ corresponds to the density of unmodulated lattice. By performing the Fourier transformation we obtain:

[TABLE]

where $\rho_{0\bf k}=F_{\bf k}\sum_{\bf b}\delta_{{\bf k},{\bf b}}$ is a sum of sharp peaks located at reciprocal wave vectors $\bf b$ of the underlying crystal lattice. Here $F_{\bf k}$ is the lattice form-factor. Integrating by parts, we obtain

[TABLE]

The terms in Eq. (44) can be combined as follows:

[TABLE]

where

[TABLE]

These terms describe appearance of the CDW peaks with wave-vectors $n{\bf Q},n=1,2,\dots$

Of particular interest to us is the first term in Eq. (45). We note that the presence of the CDW suppresses the amplitudes of the Bragg peaks by an amount

[TABLE]

We use this relation in Sec. IV.2.

It is important for our analysis that Eq. (48) remains valid also when CDW correlations are only short-ranged, while the true long-range CDW order is absent. To show this, we consider more general expression for the ionic density:

[TABLE]

Here $\phi({\bf r}_{n})$ is the phase of the order parameter. We assume that $\phi$ is a slowly varying function of ${\bf r}_{n}$. These variations are often referred to as “phasons”. When $\phi({\bf r}_{n})$ varies as a function of ${\bf r}_{n}$, the CDW order weakens, or disappears completely, and becomes replaced by short-range correlations. Generalizing Eq. (44) to account for the hot phonons, we derive

[TABLE]

Thus, in the presence of the phase variation $\phi({\bf r}_{n})$ the amplitudes of the Bragg peaks remain unchanged, c.f. Eq. (45), while the CDW amplitudes $\rho^{\bf Q}_{\bf k}$ and $\rho^{2\bf Q}_{\bf k}$ become equal to

[TABLE]

We, therefore, conclude that the CDW-induced changes in the Bragg peaks intensities carry information about the short-range CDW correlations. In particular, by means of Eq. (48) one can extract the amplitude $u$ experimentally.

Appendix B Small oscillations near equilibrium state

The TDGL sector of our formalism, Eqs. (12) and (15), contains several unknown coefficients. An important part of our study is the evaluation of these parameters consistent with the available data. An interesting possibility in this regard is to investigate the regime of small oscillations of $x$ and $y$ near the equilibrium state. The resulting theoretically determined frequency and damping factor can be compared with experimental data for the AM oscillation spectrum, which allows us to recover several parameters of our model. Since the calculations for $T<T_{\rm c}$ and $T>T_{\rm c}$ differ, they will be presented separately.

B.1 Oscillations for

$T<T_{\rm c}$

When $T<T_{\rm c}$, both order parameters $x$ and $y$ have non-zero values at equilibrium. In this regime, we parameterize small oscillations as

[TABLE]

where both $\delta x$ and $\delta y$ are complex variables. Writing $\delta x$ and $\delta y$ as sums of real and imaginary parts $\delta x=\delta x^{\prime}+i\delta x^{\prime\prime}$ and $\delta y=\delta y^{\prime}+i\delta y^{\prime\prime}$, we derive the following system of linearized equations

[TABLE]

In these equations, the dynamics of the real and imaginary components are decoupled from each other.

We analyze first the frequencies of the real components in Eqs. (57) and (57). The insertion of an ansatz $\delta x^{\prime}(t)=Xe^{\lambda t}$ and $\delta y^{\prime}(t)=Ye^{\lambda t}$ leads to an equation for $\lambda$:

[TABLE]

where $P(\lambda)$ is a cubic polynomial defined as

[TABLE]

Among three roots of $P(\lambda)$, one is always real and negative. Dependent on parameters, two other roots are either (i) both complex and conjugated to each other, or (ii) both real negative. In case (i), the pair of complex roots represents AM. We identify ${\rm Im}\lambda$ with the frequency $\omega_{\rm AM}$, while $-{\rm Re}\lambda$ is the AM damping parameter $\gamma_{\rm AM}$. The calculated values of $\omega_{\rm AM}$ and $\gamma_{\rm AM}$, as functions of temperature, are plotted in Fig. 11.

Examining panels (a) and (b) of Fig. 11 we notice that the temperature dependence of the AM exhibits two different qualitative regimes determined by the model parameters. Frequency $\omega_{\rm AM}$ plotted in panel (a) remains zero in some finite vicinity of $T_{\rm c}$. As for panel (b), $\omega_{\rm AM}$ never vanishes. Following Refs. 54; 55, where this dichotomy was previously analyzed, we refer to the behavior shown in panel (a) of Fig. 11 as “adiabatic”, while the one shown is in panel (b) is to be called “non-adiabatic”.

To determine the border between the adiabatic and non-adiabatic regimes, we need to analyze $\omega_{\rm AM}$ at $T=T_{\rm c}$. This condition corresponds to $\Theta=0$. As a result, Eq. (58) becomes easily solvable:

[TABLE]

where

[TABLE]

Quantity ${\cal D}$ is important for our analysis. Specifically, Eq. (60) implies that, if ${\cal D}>0$, then $\omega_{\rm AM}=0$, otherwise, it is finite: $\omega_{\rm AM}=|{\rm Im}\lambda_{2,3}|=\frac{1}{2\tau_{0}}\sqrt{-\cal D}$. Therefore, the condition ${\cal D}=0$ separates the adiabatic and non-adiabatic regimes. (In terms of experiment, it might be difficult to detect the difference between a formally adiabatic case $\omega_{\rm AM}(T=T_{\rm c})=0$, and a non-adiabatic case characterized by inequality $\gamma_{\rm AM}(T=T_{\rm c})\gg\omega_{\rm AM}(T=T_{\rm c})$.)

Since the values of $\omega_{\rm AM}$ and $\gamma_{\rm AM}$ at $T=300$ K are known from experiment, see Eq. (20), we can use them to derive constraints on the TDGL parameters. To obtain the constraints, we re-write Eq. (58) as two real-valued equations

[TABLE]

These equations reduce the number of free TDGL coefficients from four ($\omega_{0}$, $\tau_{0}$, $\zeta$, and $\gamma_{y}$) to two. In the main text, we treat $\tau_{0}$ and $\omega_{0}$ as free parameters. Within such a convention, Eqs. (62) and (63) can be used to define two implicit functions $\gamma_{y}=\gamma_{y}(\omega_{0},\tau_{0})$ and $\zeta=\zeta(\omega_{0},\tau_{0})$.

Damping parameter $\gamma_{y}$ must always be non-negative, i.e.

[TABLE]

which further limits the available space for $\omega_{0}$ and $\tau_{0}$ as discussed in the main text, see also Fig. 1.

Now we can analyze Eqs. (57)–(57), which describe the oscillations of the imaginary components $\delta x^{\prime\prime}$ and $\delta y^{\prime\prime}$. One can check that, in this case, we also have three modes whose eigenfrequencies are given by Eq. (60). The zero eigenfrequency represents a Goldstone mode. Within our model small oscillations of $\delta x^{\prime\prime}$ and $\delta y^{\prime\prime}$ have temperature-independent frequencies and damping parameters. This is a consequence of our assumption that quantity $a$ in the Landau functional (3) is the only one dependent on temperature.

As with the real components, the dynamics of the imaginary components $\delta x^{\prime\prime}$ and $\delta y^{\prime\prime}$ is sensitive to the sign of ${\cal D}$. Specifically, in the adiabatic regime, the roots $\lambda_{2,3}$ are both real negative, i.e. the time evolution is overdamped. In the non-adiabatic regime, the roots form a complex conjugate pair, which corresponds to underdamped oscillations.

B.2 Oscillations for

$T>T_{\rm c}$

When $T\geq T_{\rm c}$, the equilibrium values of $x$ and $y$ are zero. Thus, in the regime of linear oscillations in the disordered phase one writes $x(t)=\delta x(t)$ and $y(t)=\delta y(t)$. The resulting linearized equations for the real components coincide with the equations for the imaginary components:

[TABLE]

The eigenfrequencies then satisfy the equation

[TABLE]

Naturally, at the transition (${\Theta}=0$), Eq. (B.2) and Eq. (58) are identical. This ensures that all eigenfrequencies smoothly cross $T_{\rm c}$. Figure 11 shows the numerically calculated eigenfrequencies for both adiabatic and non-adiabatic regimes.

Appendix C Electronic heat capacity

C.1 Temperature dependence of the electronic heat capacity

Here we further motivate Eq. (24) for the temperature dependence of the electronic heat capacity. We already mentioned that LaTe3, despite the presence of the CDW order, is not an insulator, but rather is a metal, with ungapped fragments of the Fermi surface and finite density of states at the Fermi energy $\epsilon_{F}=0$, cf. Fig. 14. Consequently, the low-temperature heat capacity demonstrates [​​ 46] metallic behavior expressed by Eq. (23). However, unlike a “classical” metal for which deviations from linear relation $C_{\rm e}\propto T_{\rm e}$ for $T_{\rm e}\lesssim 2000$ K are generally weak [​​ 90], we expect that the CDW order in LaTe3 affects the validity of Eq. (23) in the above temperature range.

Available thermodynamic and ab initio data support [​​ 46] this expectation: the CDW order suppresses coefficient $c_{0}$ in Eq. (23) almost twofold relative to its value in the hypothetical situation without the CDW order. We assume that this suppression is due to the “pseudo-gapped” single-electron density of state (DOS) $\nu(\epsilon)$ sketched in Fig. 12. In this sketch, the uniform metallic DOS $\nu_{0}$ is modified by the presence of the CDW order. At the Fermi energy $\epsilon_{F}$, we choose it for concreteness to be two times smaller than the bare value $\nu_{0}$. This suppression is caused by the expulsion of the electronic states from the gapped parts of the Fermi surface to higher energies. Since coefficient $c_{0}\propto\nu(\epsilon_{F})$, the value of $c_{0}$ decreases together with $\nu(\epsilon_{F})$. As for the states excluded from the vicinity of $\epsilon_{F}$, they accumulate at [​​ 51] $|\epsilon-\epsilon_{F}|\sim\Delta\sim 0.35$ eV. When $|\epsilon-\epsilon_{F}|\gtrsim 0.6$ eV, the DOS returns to its bare value $\nu_{0}$. The resulting function $\nu(\epsilon)$ exhibits pronounced variations on the scale of hundreds of meV, which leads to a non-linear temperature-dependence of the heat capacity shown in Fig. 13. This temperature dependence is calculated using the expression:

[TABLE]

The plot in Fig. 13 indicates that $C_{\rm e}(T_{\rm e})$ departs from the low-temperature linear dependence, Eq. (23), for $T\gtrsim 300$ K. Figure 13 also shows the plot of $C_{\rm e}(T_{\rm e})$ for a simple piecewise linear function given by Eq. (24). It can be seen in this figure that Eq. (24) adequately approximates $C_{\rm e}(T_{\rm e})$ in the temperature range of interest.

We note that expression (68) is formulated under the assumption that $\nu(\varepsilon)$ is independent of $T_{\rm e}$. This assumption is, likely, violated in LaTe3, because the actual DOS is sensitive to the values of the order parameters $x$ and $y$, both of which are temperature- and time-dependent quantities. Thus, application of Eq. (24) to the non-equilibrium situations should be taken with caution. However, we do expect that the piecewise linear function Eq. (24) would still constitute a reasonable approximation to the actual temperature dependence of the electronic heat capacity. Let us also emphasize that accurate knowledge of the electronic heat capacity is important only during the first rapid stage of the electronic temperature relaxation. At the second slow stage, only lattice contributions to the heat capacity are relevant, see Eq. (27).

C.2 Estimate of electronic temperature from tr-ARPES

The generation of hot carriers after strong photoexcitation is followed by thermalization within the electronic subsystem on a time scale $\lesssim 100$ fs [​​ 92]. Using time- and angle-resolved photoemission spectroscopy (tr-ARPES), one can estimate the electronic temperature after the initial thermalization by fitting the energy distribution of quasiparticles to a Fermi-Dirac distribution [​​ 93].

Figure 14a shows the Fermi surface in LaTe3 at $T=15$ K$\ll T_{\rm c}$ before the arrival of the pump laser pulse (see Ref. 8 for measurement details). The tr-ARPES intensity is absent for the most parts of the Fermi surface due to the opening of the CDW gap. The remaining Fermi surface is consistent with the previous reports [​​ 91; 94]. In order to minimize complications arising from the transient suppression of the CDW gap, we focus on the ungapped part of the Fermi surface at equilibrium. In Fig. 14b, we present an energy-momentum cut through the metallic part of the Fermi surface where the Te $5p_{x}/p_{z}$ bands cross the Fermi level $\epsilon_{F}$. The same cut is shown after photoexcitation at a pump-probe time delay of $t=250$ fs, where states above $\epsilon_{F}$ are transiently populated.

To quantitatively analyze the carrier redistribution after photoexcitation, we plot the energy distribution curves (EDCs) at $k_{z}=0.08$ Å*-1* indicated by the arrow in Fig. 14b. At $t<t_{0}$, there is a sharp cutoff of EDC around $\epsilon_{F}$; this feature is replaced by a long tail at $\epsilon>\epsilon_{F}$ at 250 fs. The temporal evolution of the EDC across $\epsilon_{F}$ can be captured by the following model [​​ 95; 93]:

[TABLE]

where $\nu(\epsilon)$ is the density of states, $f[\cdot]$ is the Fermi-Dirac distribution that depends on the chemical potential $\mu$ and the electronic temperature $T_{\rm e}$. The terms in $\{\cdot\}$ are energy convoluted with a Gaussian kernel $\tilde{g}[\cdot]$, cf. Eq. (31), whose time-dependent [​​ 96; 93] width parameter $\tilde{w}(t)$ arises from the finite energy resolution of the instrument and from spectral broadening due to increased scattering rate after photoexcitation. The density of states $\nu(\epsilon)$ is assumed to remain unchanged over time; it is determined by the EDCs before photoexcitation. This assumption is largely justified, because that particular part of the Fermi surface is minimally affected by the transient suppression and recovery of the CDW gap. To limit the number of free parameters, we adopt a linear approximation $\nu(\epsilon)\approx\nu(\epsilon_{\rm F})+(\epsilon-\epsilon_{\rm F})\alpha_{\nu}$, where $\alpha_{\nu}$ is an adjustable parameter. This, in turn, limits the fitting range as indicated in Fig. 14b, since the strong intensity of a separate band at high binding energy cannot be captured with a linear density of states. In summary, the time-dependent fitting parameters include $\mu(t)$, $T_{\rm e}(t)$, and $\tilde{w}(t)$; the value of $T_{\rm e}$ before the arrival of the laser pulse is fixed to be 15 K, which is the base temperature of the sample during the measurement.

Figure 14c shows the extracted $T_{\rm e}$ throughout the photoexcitation event for three different excitation densities given in terms of the number of absorbed photons per unit volume. Though electronic temperature obtained within a few $\tau_{0}$ after photoexcitation is less reliable due to the non-thermal nature of the carrier distribution [​​ 92], values at longer time delays are indicative of the quasi-thermal state of the electronic subsystem with effective electronic temperature $T_{\rm e}$ and a Fermi-Dirac distribution (see the fits in Fig. 14b). As one expects, higher transient $T_{\rm e}$ is reached at higher excitation density.

Now, using Eq. (19) with the temperature dependence $C_{\rm e}(T_{\rm e})$ given by Eq. (24) and with the values of $\hbar\omega_{\gamma}$, ${\cal V}$, $F$, $T_{\rm env}$, and $T_{\rm e}$ corresponding to the experiment-based plots in Fig. 14c, we obtain the possible range of values $3-5$ mJ$/$mol K2 for the parameter $c$ entering Eq. (24). In the actual simulations, we use $c=4$ mJ$/$mol K2.