### Nuclear resonant scattering

Table of Contents

The forward transmission of an arbitrary X-ray pulse *E*_{in}(*t*) in an extended resonant sample is^{5,37}

$${E}_{{rm{out}}}(t)={E}_{{rm{in}}}(t),* T(t),$$

(1)

where *T* is a characteristic transmission function and the asterisk denotes a convolution. Neglecting electronic absorption, one can write

$$T(t)=delta (t)+R(t),$$

(2)

where *δ*(*t*) is the Dirac delta function and *R*(*t*) denotes the response function of the nuclear target, that is, the scattered X-rays. This response function *R*(*t*) is directly related to the spatially averaged nuclear magnetic transition dipole moment (langle hat{d}(t)rangle ) induced by the X-ray light in the target, which forms the primary quantity of interest in this work. Realistic transmission functions *T*(*t*), which we use to model the experimental data, can be computed with software packages such as conuss^{38}. The dispersive and absorptive properties of the electronic background are spectrally broad and are included as a constant factor. For a material featuring a single-line resonance, *T*(*t*) can be expressed analytically. Omitting the free phase evolution exp(i*ω*_{0}*t*), we have^{37,39}

$$T(t)=delta (t)-theta (t),frac{gamma {D}_{{rm{M}}}}{2}{{rm{e}}}^{-frac{gamma }{2}t}frac{{J}_{1}(sqrt{gamma {D}_{{rm{M}}}t})}{sqrt{gamma {D}_{{rm{M}}}t}}.$$

(3)

Here, *θ*(*t*) is the Heaviside step function, *J*_{1} is the Bessel function of the first kind, *D*_{M} = *σ*_{0}*fnd* is the Mössbauer optical thickness of the resonant target, *n* is the volume density of the resonant nuclei, *d* the target thickness, *σ*_{0} the cross-section, *f* the Lamb−Mössbauer factor, and *γ* the resonance width. For the data shown in Extended Data Fig. 1, we used Mössbauer thicknesses *D*_{SCU} = 20 and *D*_{T} = 9.2, which optimally mimic the experimentally realized setting. For the analytical calculations, we generalize equation (3) to write the response functions of our SCU and target samples in the coherent control setting as

$$begin{array}{cc}{T}_{{rm{S}}{rm{C}}{rm{U}}}(t) & ,=delta (t)+{{rm{e}}}^{{rm{i}}varphi (t)}{R}_{{rm{S}}{rm{C}}{rm{U}}}(t)\ & ,approx delta (t)-{{rm{e}}}^{{rm{i}}varphi (t)}theta (t)frac{gamma {D}_{{rm{S}}{rm{C}}{rm{U}}}}{2}frac{{J}_{1}(sqrt{gamma {D}_{{rm{S}}{rm{C}}{rm{U}}}t})}{sqrt{gamma {D}_{{rm{S}}{rm{C}}{rm{U}}}t}}{{rm{e}}}^{-frac{gamma }{2}t}[1+{{rm{e}}}^{-{rm{i}}St}],end{array}$$

(4)

and

$${T}_{{rm{T}}}(t,delta )=delta (t)+{R}_{{rm{T}}}(t,delta )=delta (t)-theta (t)frac{gamma {D}_{{rm{T}}}}{2}frac{{J}_{1}(sqrt{gamma {D}_{{rm{T}}}t})}{sqrt{gamma {D}_{{rm{T}}}t}}{{rm{e}}}^{-frac{gamma }{2}t}{{rm{e}}}^{-{rm{i}}delta t}.$$

(5)

Here, the two resonances in the SCU target are accounted for by the part in brackets containing the frequency difference *S* of the two resonances (see Fig. 1c). We note that the ‘approximately equal to’ in equation (4) indicates that this formulation assumes that *S* is sufficiently large to treat the response of the two resonances separately, which is well justified in our SCU sample. The target detuning relative to one of the resonances of the SCU is *δ*.

### SCU operation

Excited by a short *δ*(*t*)-like X-ray pulse and in the case of no motion, the field behind the SCU given in equation (1) reduces to equation (2).

To tune the relative phase between the *δ*(*t*) component and the scattered part *R*_{SCU}(*t*) in equation (4), a motion *x*(*t*) is applied to the SCU. This results in the combined field^{7,10,40}

$${E}_{{rm{SCU}}}(t)={E}_{{rm{exc}}}(t)+{E}_{{rm{control}}}(t)=delta (t)+{{rm{e}}}^{{rm{i}}varphi (t)}{R}_{{rm{SCU}}}(t),$$

(6)

$$varphi (t)=k[x(t)-x(0)],$$

(7)

where *k* = *ω*_{0}/*c* is the wavenumber. In our experiment, we use this double pulse to drive a nuclear target. Again, the downstream X-ray intensity can be computed using equation (1), where *E*_{SCU}(*t*) now takes the role of the input field *E*_{in}(*t*) and *T*(*t*) corresponds to the transmission function of the actual target.

### Target response

To describe the dynamics of the target nuclei without having to impose a particular model, we write the output field behind that target in terms of the input field delivered by the SCU and a scattering component,

$${E}_{{rm{out}}}(t,delta )={E}_{{rm{SCU}}}(t)+alpha langle hat{d}(t,delta )rangle ,$$

(8)

where *α* is a constant. To calculate (langle hat{d}(t,delta )rangle ), using equation (5), we find in Fourier space

$${tilde{E}}_{{rm{out}}}(omega ,delta )={tilde{E}}_{{rm{SCU}}}(omega )+alpha langle hat{d}(omega ,delta )rangle ={tilde{E}}_{{rm{SCU}}}(omega )+{tilde{R}}_{T}(omega ,delta ){tilde{E}}_{{rm{SCU}}}(omega ),$$

(9)

such that

$$langle hat{d}(omega ,delta )rangle =frac{1}{alpha }{tilde{R}}_{T}(omega ,delta ){tilde{E}}_{{rm{SCU}}}(omega ).$$

(10)

Next, we consider the position-resolved dynamics. At a depth *x* inside the target of length *L*, we can write the Fourier transforms of *T*_{T}(*t*, *δ*) and *R*_{T}(*t*, *δ*) as ({tilde{T}}_{{rm{T}}}(omega ,delta ,x)={{rm{e}}}^{a(omega ,delta )x}) and ({tilde{R}}_{{rm{T}}}(omega ,delta ,x)={tilde{T}}_{{rm{T}}}(omega ,delta ,x)-1). (a(omega ,delta )) is the response function for a thin slice of the target, (mathop{R}limits^{ sim }{}_{{rm{T}}}^{{rm{t}}{rm{h}}{rm{i}}{rm{n}}}(omega ,delta )=a(omega ,delta )). For the case of a single target resonance, *a*(*ω*, *δ*) = – i[*D*_{M}*γ* /(4*L*)]/ ((omega -delta +{rm{i}}gamma /2)). Using these definitions, the total field at position *x* in the target is ({tilde{E}}_{{rm{SCU}}}(omega ){tilde{T}}_{{rm{T}}}(omega ,delta ,x)), and the position-dependent transition dipole moment in a thin sample slice at *x* becomes

$$langle hat{d}(omega ,delta ,x)rangle =frac{1}{alpha }{tilde{E}}_{{rm{SCU}}}(omega ){tilde{T}}_{{rm{T}}}(omega ,delta ,x){tilde{R}}_{{rm{T}}}^{{rm{thin}}}(omega ,delta ).$$

(11)

A spatial average of this position-dependent dipole moment over the entire target length is straightforward using

$$overline{{tilde{T}}_{{rm{T}}}(omega ,delta ,x)}=frac{1}{L}{int }_{0}^{L}{tilde{T}}_{{rm{T}}}(omega ,delta ,x){rm{d}}x=frac{{tilde{T}}_{{rm{T}}}(omega ,delta ,L)-1}{a(omega ,delta )}=frac{{tilde{R}}_{{rm{T}}}(omega ,delta ,L)}{{tilde{R}}_{{rm{T}}}^{{rm{thin}}}(omega ,delta )}cdot $$

(12)

Inserting this relation into equation (11), we obtain

$$overline{langle hat{d}(omega ,delta ,x)rangle }=frac{1}{alpha }{tilde{E}}_{{rm{SCU}}}(omega ){tilde{R}}_{{rm{T}}}(omega ,delta ,L)=langle hat{d}(omega ,delta )rangle .$$

(13)

Thus, the quantity (langle hat{d}(t,delta )rangle ) defined in equation (8) is equal to the spatial average over the position-dependent nuclear magnetic transition dipole moment induced by the X-ray light as obtained from a full propagation analysis.

The spatially averaged dipole moment (langle hat{d}(t,delta )rangle ) has the crucial advantage that it can be evaluated without requiring knowledge about the position-dependent dynamics inside the target, which is not accessible in our experiment. From equation (10), we find that measuring the complex field amplitude of the double pulse delivered by the SCU, and determining the target response function *R*_{T}(*t*, *δ*) using fits to its individual response at rest, already allow us to evaluate (langle hat{d}(t,delta )rangle ).

### Intensity crossover

When comparing the two coherent control SCU operations, differences are found in the temporal structure of the X-ray field behind the target (Fig. 2). In particular, the most prominent qualitative feature for the cases considered here is a crossover of the dominating intensity after a certain time. This behaviour can directly be linked to the target dynamics induced by the SCU pulse. Behind both targets, the amplitude at the detector follows from equations (4) and (5) with resonant target (*δ* = 0) as

$${E}_{{rm{out}}}(t)=mathop{underbrace{delta (t)}}limits_{{rm{SCU}},{rm{pulse}},1}+mathop{underbrace{{{rm{e}}}^{{rm{i}}varphi (t)},{R}_{{rm{SCU}}}(t)}}limits_{{rm{SCU}},{rm{pulse}},2}+mathop{underbrace{{R}_{{rm{T}}}(t)}}limits_{{rm{target}},{rm{response}},({rm{SCU}},{rm{pulse}},1)}+mathop{underbrace{{{rm{e}}}^{{rm{i}}varphi (t)},{R}_{{rm{SCU}}}* {R}_{{rm{T}}}(t)}}limits_{{rm{target}},{rm{response}},({rm{SCU}},{rm{pulse}},2)},$$

(14)

where the interpretation of each part is denoted by the underbrace text. The first two contributions correspond to the double-pulse equation (6) delivered by the SCU unit onto the target. The last two contributions are the target response induced by the two parts of the SCU double pulse, where the asterisk denotes the convolution. For definiteness, in the following, we consider the case *D*_{SCU} > *D*_{T}, as in our experiment.

At short times 0 < *γt* ≪ ({D}_{{rm{SCU}}}^{-1}) immediately after the excitation at *t* = 0, the field at the detector for the cases of enhanced emission (that is, *ϕ*(*t*) = 0) and enhanced excitation (that is, (varphi (t > 0)-varphi (0)={rm{pi }})) evaluate to

$${E}_{{rm{o}}{rm{u}}{rm{t}}}^{{rm{e}}{rm{m}}{rm{i}}{rm{s}}{rm{s}}{rm{i}}{rm{o}}{rm{n}}}(t)mathop{longrightarrow }limits^{0 < gamma tll {D}_{{rm{S}}{rm{C}}{rm{U}}}^{-1}}-frac{gamma }{4}({D}_{{rm{T}}}+2{D}_{{rm{S}}{rm{C}}{rm{U}}})+{mathcal{O}}(t),$$

(15)

$${E}_{{rm{o}}{rm{u}}{rm{t}}}^{{rm{e}}{rm{x}}{rm{c}}{rm{i}}{rm{t}}{rm{a}}{rm{t}}{rm{i}}{rm{o}}{rm{n}}}(t)mathop{longrightarrow }limits^{0 < gamma tll {D}_{{rm{S}}{rm{C}}{rm{U}}}^{-1}}-frac{gamma }{4}({D}_{{rm{T}}}-2{D}_{{rm{S}}{rm{C}}{rm{U}}})+{mathcal{O}}(t).$$

(16)

Thus, (|{E}_{{rm{out}}}^{{rm{emission}}}(t){|}^{2} > |{E}_{{rm{out}}}^{{rm{excitation}}}(t){|}^{2}) at early times, that is, the detected intensity is initially higher in the enhanced emission case than in the enhanced excitation case.

To identify the subsequent intensity crossover, we next consider the time evolution of the different contributions to *E*_{out}(*t*). Apart from the oscillation due to the two resonances in the SCU response, *R*_{T}(*t*) and *R*_{SCU}(*t*) are both negative immediately after the excitation at *t* = 0, and then decay in magnitude as time progresses until they vanish at their respective first zeros of the Bessel functions *J*_{1}. Since *D*_{SCU} > *D*_{T}, the zero of the Bessel function in *R*_{SCU}(*t*) is reached first, and we denote this time as *t*_{min}. Up to this point in time, *R*_{T}(*t*) remains negative, and *R*_{SCU} * *R*_{T}(*t*) is positive, as follows immediately from the definition of the convolution. Thus, at *t*_{min},

$${E}_{{rm{out}}}^{{rm{emission}}}({t}_{{rm{min }}})=mathop{underbrace{{R}_{{rm{T}}}({t}_{{rm{min }}})}}limits_{ < 0}+mathop{underbrace{{R}_{{rm{SCU}}}ast {R}_{{rm{T}}}({t}_{{rm{min }}})}}limits_{ > 0},$$

(17)

$${E}_{{rm{o}}{rm{u}}{rm{t}}}^{{rm{e}}{rm{x}}{rm{c}}{rm{i}}{rm{t}}{rm{a}}{rm{t}}{rm{i}}{rm{o}}{rm{n}}}({t}_{min})=mathop{underbrace{{R}_{{rm{T}}}({t}_{min})}}limits_{ < 0}-mathop{underbrace{{R}_{{rm{S}}{rm{C}}{rm{U}}}ast {R}_{{rm{T}}}({t}_{min})}}limits_{ > 0}.$$

(18)

Therefore, (|{E}_{{rm{out}}}^{{rm{excitation}}}({t}_{{rm{min }}}){|}^{2} > |{E}_{{rm{out}}}^{{rm{emission}}}({t}_{{rm{min }}}){|}^{2}), that is, the detected intensity as a function of time is now higher in the enhanced excitation case, thus proving the intensity crossover. Interestingly, at time *t*_{min}, the output field is equal to the scattered response of the target, because the SCU field contribution vanishes.

The intensity crossover observed in the experimental data shown in Fig. 2 and in the full theory calculations shown in Fig. 2 and Extended Data Fig. 1 is therefore directly linked to the coherent control. At early times, the responses of the target and the SCU are in phase for the enhanced emission case and thus add up to a higher initial intensity, whereas in the enhanced excitation case, the SCU and the target contributions have opposite phase owing to the piezo displacement, and therefore destructively interfere to give a lower intensity. Because of this relative phase, the target excitation is increased in one case (excitation), and decreased in the other case (emission). At later times around *t*_{min}, the output field coincides with the response of the target, and the intensity in the enhanced excitation case is higher than in the enhanced emission case. The higher intensity in the enhanced excitation case can thus be directly attributed to a higher absolute value of the induced average target dipole moment compared to the enhanced emission case.

### Propagation effects in the target

In any target of finite thickness, the dynamics of the induced magnetic transition dipole moments will vary as a function of position in the target, since they are driven not only by the externally applied field, but also by the field scattered by the upstream dipoles. To determine these propagation effects, we treat the target as a medium of two-level atoms, and calculate the propagation of the SCU pulse through the target using the Maxwell−Bloch equations in the slowly varying envelope approximation^{40,41}

$$left(frac{partial }{partial x}+frac{1}{c}frac{partial }{partial t}right)varOmega (x,t)=-frac{{rm{i}}gamma {D}_{{rm{M}}}}{2L}{rho }_{eg}(x,t),$$

(19)

where (varOmega (x,t)=-2d {mathcal E} (x,t)/hbar ) is the Rabi frequency, with ( {mathcal E} ) the slowly varying amplitude of the propagating field, *d* the magnetic dipole moment, *L* the target length, and *ρ*_{eg}(*x*, *t*) the density matrix element corresponding to the coherence between the ground state ({|g}rangle ) and the excited state ({|e}rangle ) induced by the propagating field. It follows from the nuclear dynamics governed by the equations of motion for the density operator (hat{rho })

$$frac{partial }{partial t}hat{rho }=frac{1}{{rm{i}}hbar }[V,hat{rho }]-frac{gamma }{2}(|erangle langle e|hat{rho }+hat{rho }|erangle langle e|-2|grangle langle e|hat{rho }|erangle langle g|),$$

(20)

$$V=-hbar Delta |erangle langle e|+frac{hbar }{2}(varOmega (x,t)|erangle langle g|+varOmega {(x,t)}^{ast }|grangle langle e|).$$

(21)

Results of this analysis are shown in Extended Data Figs. 2 and 3 for the parameters relevant to our experiment. Extended Data Fig. 2 shows that there are indeed propagation effects, that is, the dipole dynamics depends on the position in the target because of the light scattered by the upstream nuclei. Nevertheless, the coherent control acts similarly everywhere inside the target: In the enhanced emission case, the excitation due to the first pulse is always rapidly driven back to the ground state by the second pulse. In the enhanced excitation case, the excitation due to the first pulse is always increased by the second pulse. To illustrate this feature in more detail, Extended Data Fig. 3 compares the dipole dynamics at the target entry (*x* = 0), in the middle of the target (*x* = *L*/2), and at the end of the target (*x* = *L*). At all positions in the target, the two coherent control cases are clearly visible. Finally, the results shown as dashed lines in Fig. 3 are obtained by averaging the spatially resolved dipole dynamics in Extended Data Fig. 2 over the sample length.

### Quantum optical two-level model

In the limit of a thin target, the dynamics in the target can be modelled from first principles, using an approach based on a two-level-system description for the resonant target. Even though we do not use this limit in our data analysis, the calculation provides a clear interpretation of the spatially averaged magnetic-dipole moment defined in equation (8) in terms of the microscopic nuclear transition dipole moments, and illustrates how a nuclear two-level quantum system coherently controllable via the double pulses from the SCU can be implemented. In the thin-sample limit and at weak excitation, the two-level description is known to agree with the nuclear resonant scattering approach described above. In the following, we exploit this equivalence to establish an expression for the target dipole moment. The two-level system is formed by one collective ground state (|grangle ) and one collective excited state (|erangle ). The driving with an X-ray field *E*_{in}(*t*) is described by the Hamiltonian

$$H=frac{hbar varOmega (t)}{2}|erangle langle g|+frac{{hbar varOmega }^{ast }(t)}{2}|grangle langle e|,$$

(22)

where (varOmega (t)=-2d{E}_{{rm{i}}{rm{n}}}(t)/hbar ), with *d* being the magnetic dipole moment. Additionally, we include spontaneous decay with rate (tilde{gamma }) in terms of a density operator

$$frac{{rm{d}}}{{rm{d}}t}hat{rho }=frac{1}{{rm{i}}hbar }[H,hat{rho }]-frac{mathop{gamma }limits^{ sim }}{2}(|erangle langle e|hat{rho }+hat{rho }|erangle langle e|-2|grangle langle e|hat{rho }|erangle langle g|).$$

(23)

For weak excitation it is sufficient to consider the coherence (langle {sigma }_{ge}rangle =langle e|hat{rho }|grangle ) only. In the limit (langle g|hat{rho }|grangle =1), (langle e|hat{rho }|erangle =0), we have the equation of motion

$$frac{{rm{d}}}{{rm{d}}t}langle {sigma }_{ge}(t)rangle =-{rm{i}}frac{varOmega (t)}{2}-frac{tilde{gamma }}{2}langle {sigma }_{ge}(t)rangle ,$$

(24)

which is solved by (with initial conditions (langle {sigma }_{ge}rangle =0=varOmega ) at *t* = −∞)

$$begin{array}{cc}langle {sigma }_{ge}(t)rangle & ,=-frac{{rm{i}}}{2}{int }_{-{rm{infty }}}^{t}varOmega ({t}^{{prime} }){{rm{e}}}^{-frac{mathop{gamma }limits^{ sim }}{2}(t-{t}^{{prime} })}{rm{d}}{t}^{{prime} }\ & ,=-frac{{rm{i}}}{2}{int }_{-{rm{infty }}}^{{rm{infty }}}varOmega ({t}^{{prime} })theta (t-{t}^{{prime} }){{rm{e}}}^{-frac{mathop{gamma }limits^{ sim }}{2}(t-{t}^{{prime} })}{rm{d}}{t}^{{prime} }\ & ,=-frac{{rm{i}}}{2}varOmega (t)ast [theta (t){{rm{e}}}^{-frac{mathop{gamma }limits^{ sim }}{2}t}].end{array}$$

(25)

The field behind the two-level system is composed of the initial field and a scattered contribution^{42}

$${E}_{{rm{out}}}^{{rm{TLS}}}(t)={E}_{{rm{in}}}(t)+alpha langle hat{d}(t)rangle ,$$

(26)

where (langle hat{d}rangle =dlangle {sigma }_{ge}rangle ) is the dipole response of the two-level system, and *α* is a constant, also taking into account the extended sample geometry^{43}. In particular, for (alpha =-2{rm{i}}bhbar /{d}^{2}) and (tilde{gamma }=gamma (1+{D}_{M}/4)) we have

$${E}_{{rm{o}}{rm{u}}{rm{t}}}^{{rm{T}}{rm{L}}{rm{S}}}(t)={E}_{{rm{i}}{rm{n}}}(t)ast left[delta (t)-theta (t)frac{gamma {D}_{M}}{4}{{rm{e}}}^{-gamma frac{1+{D}_{M}/4}{2}t}right].$$

(27)

This result is also obtained within the thin-target limit (gamma tll {D}_{M}^{-1}) of the nuclear resonant scattering theory equations (1) and (2). The analytical agreement between the two calculations demonstrates the validity of the two-level-system approach. Comparing equation (26) with equations (1) and (2), we find

$${E}_{{rm{in}}}(t)ast R(t)=alpha langle hat{d}(t)rangle ,$$

(28)

which, together with equation (10), illustrates the relation between the spatially averaged target dipole moment and the microscopic dipole moments, and highlights the correspondence of the response function in the nuclear resonant scattering approach with the time-dependent nuclear dipole moment in the quantum optical model.

### Event-based detection

In our experiment, we make use of an event-based detection system that records, among other quantities, the absolute detection time within the experimental run, the relative detection time after the excitation, and energy information for each photon separately. It thus provides access to two-dimensional time- and energy-resolved spectra, which contain the full holographic (amplitude and phase) information in its interference structures, which furthermore can be split into variable measurement intervals throughout the data analysis. This feature is crucial in two respects. First, the Allan deviation analysis requires an a posteriori splitting of the data into time bins of variable duration *τ*. This splitting is only possible if the arrival time of each photon is stored. Second, we will show below that the time-dependent intensity, which was used in previous experiments, does not provide access to the key observables studied here, namely the complex spatially averaged nuclear magnetic transition dipole moment and the stability of the coherent control scheme. To better appreciate the difference between our event-based detection and the standard time-dependent intensity measurement, it is important to note that in order to determine the nuclear dynamics, we must solve an inverse problem of extracting the nuclear dipole moment from the scattered light. The time-dependent intensity measured in previous works does not provide sufficient information to solve this inverse problem unambiguously, which is a fundamental obstacle in accessing the matter (nuclear) part of the system. We note that a phase determination in nuclear resonant scattering has previously been suggested using a velocity drive as an interferometer and phase shifter^{44}, but this reference neglects the radiative coupling between the analyser and the target. The latter leads to the coherent control reported here.

To illustrate the necessity of our event-based spectroscopy method, we consider the setup used in our experiment, with the three motions shown in Extended Data Fig. 4a. Motion 1 corresponds to a rapid jump shortly after the arrival of the X-ray pulse by half the resonant wavelength *λ*_{0}/2, which leads to the coherent enhanced excitation case. Motion 2 is a similar displacement, but in the opposite direction. Motion 3 modifies motion 1 by an additional linear drift on top of the step-like motion. As discussed in the Methods section ‘Stability and Allan deviation’, such linear drifts are the dominant source of noise expected in our setup, and the drift shown in Extended Data Fig. 4 corresponds to a temporal deviation *ξ* = 25 zs. Our stability analysis is based on the ability to reliably detect drifts of this and smaller magnitude. As shown in Extended Data Fig. 4b, c, the three motions induce different dynamics in the target nuclei, and our experiment aims at detecting these differences. We note that, somewhat counter-intuitively, motions 1 and 2 induce dynamics that do not only differ in phase, but also in the time-dependent magnitude of the induced dipole moments. The reason for this feature is that the two motions include opposite velocities in the approximately step-like part of the motion, leading to opposite transient Doppler shifts, and thus in turn to different spectra of the outgoing double pulses. Thus, the target nuclei experience different driving fields. Motion 3 differs from motion 1 by an additional drift, which translates into a corresponding additional phase dynamics of the induced dipole moments.

Extended Data Fig. 5a shows the theoretical predictions for the time-dependent intensity on resonance, which was used as an observable in previous experiments. The corresponding intensity differences obtained by subtracting the experimentally accessible intensities from each other are shown in Extended Data Fig. 5b. The results for motions 1 and 2 essentially coincide. Motion 3 only differs slightly, in the depth of the beat minima, and is essentially indistinguishable from the other motions, in particular if practical limitations on data acquisition are taken into account. Thus, we conclude that the time-dependent intensity alone is not capable of distinguishing key motions of relevance to our analysis from each other as a matter of principle, and therefore cannot distinguish the different nuclear dynamics induced in the target nuclei.

The event-based detection technique used in our experiment provides time- and energy-resolved spectra as shown in Fig. 2a, b. To illustrate the advantage of this approach, we show relative intensity differences (*I*_{2} – *I*_{1})/(*I*_{1} + *I*_{2}) of the two-dimensional (2D) spectra obtained for motions 1 and 2 in Extended Data Fig. 5c. It can be seen that the two motions lead to rich systematic structure with full visibility. Therefore, through the 2D spectra we can easily distinguish the two motions, whereas the time-dependent intensities on resonance in Extended Data Fig. 5a for the two motions provide insufficient information to distinguish them. Finally, Extended Data Fig. 5d shows the intensity differences of the three motions for sections through the measured 2D spectra at particular Mössbauer drive detunings *δ*. It can be seen that all three motions give rise to substantial intensity differences, which furthermore exhibit characteristic time-dependencies for each detuning separately. In our data analysis, we compute two-dimensional theory spectra and compare them to the entire recorded two-dimensional spectrum at once, thereby including all Mössbauer detunings in a single fit. The rich interference structures encode full tomographic (amplitude and phase) information on the light scattered by the first absorber, and lead to a strong sensitivity of the fit to the slightest deviations in the piezo motion and the nuclear dynamics. These examples clearly show that the time-dependent intensity measured in previous experiments is incapable of distinguishing motions that are crucial to our results, in contrast to the 2D time- and energy-resolved spectra recorded in our experiment.

### Reconstruction of the SCU motion and field

The reconstruction of the SCU motion was performed based on the method in ref. ^{10}. In the experiment, the duration of the periodic motional pattern *x*_{0}(*t*) of the SCU was chosen as a multiple of the synchrotron bunch clock period, and locked to the bunch clock. In this way, stable temporal shifts between the X-ray pulses and the motional pattern could be adjusted. The target was mounted on a Doppler drive, such that the relative detuning *δ* between the target resonance energy and that of the nuclei in the SCU could be tuned via the velocity *v* of the drive. Using our event-based detection system, we recorded two-dimensional time- and velocity-resolved intensities *I*(*t*, *v*) for different temporal shifts of the motional pattern. The set of shifts was chosen in such a way that the recorded time-dependent intensities span the entire motional sequence. Each measurement covers times from 18 ns to 170 ns after the excitation with the initial X-ray pulse, and the velocity was recorded in the range −0.0228 m s^{–1} to 0.0228 m s^{–1}. Using an evolutionary algorithm, we fitted the applied motional sequence to the measured data without imposing a particular model for the motion. In this step, the experimentally measured and the theoretically expected data are compared using a Bayesian log-likelihood method. For this method, we maximized the Bayesian likelihood^{45} under the assumption that the photon counts for each data point in *I*(*t*, *v*) are Poisson distributed^{46}. For a given ideal datum *n*_{theo,i} with index *i*, the probability of obtaining the experimental count number *n*_{exp,i} is then

$$P({n}_{exp ,i}|{n}_{{rm{theo}},i})propto {({n}_{{rm{theo}},i})}^{{n}_{exp ,i}}frac{{{rm{e}}}^{-{n}_{{rm{theo}},i}}}{{n}_{exp ,i},!}.$$

(29)

The likelihood for the whole experimental dataset including all data points *i* is

$$P(exp |{rm{theo}})=prod _{i}P({n}_{exp ,i}|{n}_{{rm{theo}},i}).$$

(30)

Assuming uniform priors^{45}, (P(exp |{rm{theo}})propto P({rm{theo}}|exp )), which allows for the determination of the most likely theoretical prediction given the experimental data. Thus, we calculate all *n*_{theo,i} for each motion obtained during the evolutionary algorithm, and maximize *P*(theo|exp) to choose the most likely motion. As a result of this evolutionary algorithm, we obtain the full periodic motion *x*_{0}(*t*).

### Stability and Allan deviation

The stability of our control scheme is given by the stability of the relative phase between the excitation and the control pulses experienced by the target nuclei. Since the first excitation pulse interacts with the target at *t* ≈ 0, this phase depends on the relative motion of SCU and target during the subsequent 176 ns of each experimental run. In contrast, drifts or perturbations between different runs do not affect the stability. As a result of this relative dependence, in our modelling we can equivalently attribute imperfections in the stability of our setup either to noise or drifts in the relative phase, or to corresponding perturbations in the SCU motion.

To quantify the stability of our coherent control scheme, we use the Allan deviation measure^{28}, which is obtained by the analysis illustrated in Extended Data Fig. 6. The respective recorded datasets are split into non-overlapping samples with equal sampling times *τ*. For example, for *τ* = 10 s the first sample comprises the data taken in the time range 0–10 s, the second sample is formed by the data recorded in the time range 10–20 s, and so forth. For all *N* samples obtained for a given sampling time *τ*, we determine a quantity *ϕ*_{i} characterizing the double-pulse sequence in the interval *i* in terms of a phase deviation as explained below, as well as the corresponding temporal deviation *ξ*_{i}. From the *ϕ*_{i}, the Allan deviation *σ*_{ϕ}(*τ*) can be computed according to

$${sigma }_{varphi }(tau )={left(frac{1}{2(N-1)}mathop{sum }limits_{i=1}^{N-1}{({varphi }_{i+1}-{varphi }_{i})}^{2}right)}^{frac{1}{2}}.$$

(31)

The corresponding Allan deviation *σ*_{ξ}(*τ*) in terms of the temporal deviations *ξ*_{i} is defined analogously. It remains to determine *ϕ*_{i} and *ξ*_{i} from the experimental data as a function of *τ*. However, for short measurement intervals *τ*, the experimental statistics is not sufficient for a full independent recovery of the applied double-pulse sequence. Therefore, we make use of the direct correspondence of the relative double-pulse phase and the SCU motion, and base our analysis on the SCU motion *x*_{0}(*t*) obtained as the best fit for the entire experimental dataset. In the first step, we modify *x*_{0}(*t*) using an error model, which depends on a model parameter specified below. In the second step, we fit the modified motion to the experimental data in each interval *i* of duration *τ* separately, using the model parameter for the fit. In this fit, we use the same Bayesian log-likelihood method as for the recovery of *x*_{0}(*t*). In the third step, we translate the best fit for the model parameter into the desired deviations *ϕ*_{i} and *ξ*_{i} according to the error model.

To derive an error model, we decompose the perturbation δ*x*(*t*) to the motion into frequency components as ({rm{delta }}x(t)={sum }_{omega }{x}_{omega }(0)+{a}_{omega }sin (omega t+{varphi }_{omega })), taking into account offsets *x*_{ω}(0) and relative phases *ϕ*_{ω} for each frequency component *ω* separately. For *ωt* < 1, a series expansion yields ({rm{delta }}x(t)approx {rm{delta }}x(0)+At), where ({rm{delta }}x(0)={sum }_{omega }{x}_{omega }(0)+{a}_{omega }sin ({varphi }_{omega })) and (A={sum }_{omega }{a}_{omega }omega cos ({varphi }_{omega })). Therefore, during each experimental run of 176 ns, perturbations at least for all frequencies well below about 2π/(176 ns) ≈ 10 MHz can together be summarized into a constant offset δ*x*(0) not affecting the relative phase between the two pulses, and a linear drift motion *At* randomly varying from run to run. Therefore, we use *x*_{i}(*t*) = *x*_{0}(*t*) + *A*_{i}(*t*) as our main error model, with the free parameter *A*_{i} characterizing the magnitude of the drift in each interval *i*. The parameter *A*_{i} then translates into the desired deviations as *ϕ*_{i} = *kA*_{i}*t*_{2} and *ξ*_{i} = *A*_{i}*t*_{2}/*c*, where *t*_{2} = 170 ns is the maximum time of our data acquisition, *k* is the X-ray wavenumber, and *c* is the speed of light. With this choice, *ϕ*_{i} and *ξ*_{i} quantify upper bounds for the error acquired due to the drift with parameter *A*_{i} in terms of phase and temporal deviations.

Next to the linear drift motion, we also employed two other noise models to analyse the stability of our data. First, a scaling of the expected motion by a constant factor, (x(t)=(1+s){x}_{0}(t)). For example, in the case of a π phase jump in *x*_{0}(*t*), a scaling by *s* corresponds to a phase deviation of *s*π, or alternatively a temporal shift *sT*/2. This model, for instance, takes into account fluctuations in the voltage applied to the piezo, which to a very good approximation translates into a scaling of the displacement. Second, we superimposed the base motion with a small step-like displacement, *x*(*t*) = *x*_{0}(*t*) + *dθ*(*t* − 0^{+}). The displacement *d* translates into a phase deviation of *kd* or a temporal deviation of *d*/*c*. 0^{+} indicates a time close to zero immediately after the excitation pulse has left the target. This model tests for the presence of potential phase offsets between the excitation and control pulses.

In our analysis we found that the linear model constitutes the dominant type of error. The Allan deviations for the different noise models in the case of coherent enhanced excitation are shown in Extended Data Fig. 7. While the linear noise model predicts an optimum temporal deviation of *σ*_{ξ}(*τ*) ≈ 1 zs for the given data, the uncertainties obtained from the other two models reach well below the zeptosecond scale.

### Detector dead time

In all curves shown in Extended Data Fig. 7 as well in the curves in Fig. 4 we observe unexpected fluctuations in the Allan deviations at sampling times between *τ* ≈ 10 s and *τ* ≈ 60 s. The cause for this phenomenon is a limitation of the employed data acquisition system, which occasionally suffered from dead times of a few tens of seconds, owing to overload resulting from a too-high signal rate. As a result, some data samples with respective sampling times contain only a few or even no counts, which spoils the determination of *y*_{i} and in turn leads to large Allan deviations. This effect can be removed in the data analysis by choosing the samples not according to equal measurement times, but according to equal counts. In other words, instead of the fluctuating count rate in the experiment with its dead times, a constant averaged count rate is assumed. As shown in Extended Data Fig. 8, evaluating the Allan deviation with this method indeed suppresses the fluctuations at intermediate times, which shows that they originate from the detector dead time.

### Systematic deviations throughout the initialization phase

In the Allan deviation shown in Fig. 4, it is not fully clear whether the experimentally achieved stability has already reached its limit, and only an upper bound for possible systematic effects can be given. To interpret this result and to verify our analysis, we artificially introduced systematic deviations, by recording spectra already during the initial time after starting the piezo motion, before the piezo reached stable thermal and mechanical conditions. In this initial time, systematic drifts in the deviations *ϕ*_{i} and *ξ*_{i} as a function of the measurement time may occur. The corresponding results for samples with sampling time 200 s are shown in Extended Data Fig. 9 over the full measurement period, including the initialization phase. We note that in this plot, temporally overlapping samples were analysed, in order to trace the time evolution of the deviation with a high temporal resolution. For example, the first deviation is calculated from data in the time range 0–200 s, the next deviation for the range 1–201 s, and so forth. We find that the deviations systematically drift for an initial period of about 400 s. Afterwards, only small residual fluctuations are observed over the remaining measurement time. The orange and green curves in Fig. 4 compare the Allan deviations with and without this initial phase. It can be seen that the initialization leads to a clear systematic trend of the Allan deviation as compared to the case without the initial phase: the Allan deviation begins to increase again for sampling times exceeding approximately 100 s, which is the expected behaviour in the case of systematic drifts.

### Samples

As resonant nuclear sample we used a single-line stainless-steel foil (Fe_{55}Cr_{25}Ni_{20}), with iron enriched to about 95% in ^{57}Fe and with thickness 1 μm. The X-ray double-pulse sequence was created using an α-iron foil with thickness of about 2 μm, also enriched in ^{57}Fe. An external magnet was used to align its magnetization and the setup was arranged such that only the two Δ*m* = 0 hyperfine transitions of the 14.4-keV resonance in ^{57}Fe were driven. To displace the α-iron foil we employed a piezoelectric transducer consisting of a polyvinylidene fluoride (PVDF) film (thickness 28 μm, model DT1-028K, Measurement Specialties, Inc.). The piezo was glued on an acrylic glass backing and was driven by an arbitrary function generator (model Keysight 81160A-002).