Quantum coherent diffractive imaging

Coherent diffractive imaging (CDI) has enabled the structural analysis of individual free nanoparticles in a single shot and offers the tracking of their light induced dynamics with unprecedented spatial and temporal resolution. The retrieval of structural information from scattering images in CDI experiments is so far mostly based on the assumption of classical scattering in linear response. For strong laser fields and in particular for resonant excitations, both the linear and the classical description may no longer be valid as population depletion and stimulated emission become important. Here we develop a density matrix-based scattering model in order to include such quantum effects in the local medium response and explore the transition from linear to non-linear CDI for the resonant scattering from Helium nanodroplets. We find substantial departures from the linear response case already for experimentally reachable pulse parameters and show that non-linear spatio-temporal excitation dynamics results in rich features in the scattering, leading to the proposal of quantum coherent diffractive imaging as a promising novel branch in strong-field XUV and x-ray physics.


Introduction
Single-shot coherent diffractive imaging (CDI) has enabled the ultrafast characterization of the structure and dynamics of individual isolated nanoparticles in free flight [1][2][3][4]. This unique capability has led to various fundamental scientific findings, including the shape analysis of viruses [5], tracking of laser-driven damage in clusters on femtosecond time scales [6], the discovery of previously unresolved metastable shapes in the growth process of metal clusters [7], and the direct visualization of quantum vortices in superfluid helium nanodroplets [8], to name only a few.
The intense, ultrashort short-wavelength pulses required to record meaningful single-shot scattering patterns are typically only available at XUV and x-ray free electrons lasers (FELs) [9,10]. While CDI with hard x-rays for maximal spatial resolution is likely to remain the domain of XFELs, lab-based sources have become sufficiently brilliant to enable CDI in the XUV domain. Despite lower nominal spatial resolution, XUV-CDI enables the measurement of wide-angle scattering and thus provides three-dimensional structural information [7], which is particularly valuable for the characterization of non-reproducible targets [11].
Only recently, lab-based single-shot CDI of single Helium nanodroplets could be demonstrated with a laser-driven table-top high-harmonic generation (HHG) source [12]. Besides marking a breakthrough regarding accessibility, laser-based sources open unprecedented conceptual opportunities for CDI applications. The most obvious one results from the exquisite temporal control of optical laser systems and promises extreme time resolution that is likely to reach the sub-fs range for CDI with attosecond XUV pulses [13]. The second aspect, which is relevant also for FELs if timing is less critical, is associated with the prospect of flexible multi-color CDI scenarios. Note that, despite limitations regarding pulse duration, seeded FELs can provide widely adjustable multi-color beams, even with accurate phase control [14]. HHG-based multi-color scattering patterns have already been shown to enable the single-shot characterization of optical properties of non-reproducible targets [12]. Most importantly, because of the opportunity to provide intense pulses with durations comparable to the coherence times of XUV excitations, XUV-CDI with HHG sources (or appropriately seeded FELs) is expected to open unprecedented routes to monitor and temporally probe quantum coherent dynamics in excited nanostructures.
So far, the analysis of CDI experiments is mostly based on the assumptions of classical scattering and a linear medium response. In this case the scattering process is described by solving Maxwell's equations using a linear complex refractive index [7] or tabulated atomic scattering factors [13] or by employing appropriate approximations such as the first Born approximation [15]. However, for strong laser fields and especially for excitations near atomic resonances, the linear description may no longer be valid as effects like population transfer and depletion may become non-negligible. Specifically, the importance of quantum coherent dynamics, such as Rabi-cycling, stimulated emission, and induced transparency in the context of CDI are essentially unknown. Non-linear effects and transient changes of the optical response due to ionization have been included, e.g. via an effective damage-induced reduction of atomic form factors [16] or using explicit classical models for the resulting microscopic plasma motion [17]. Recent studies on the prospects of incoherent scattering from fluorescence and Compton scattering [18,19] and superfluorescence in the saturation regime [20] indicate the potential of using quantum effects in CDI scenarios. The so far missing firm understanding of the impact of spatio-temporal coherent bound-state dynamics is of paramount importance to (i) justify applicability of established CDI reconstruction methods employing linear scattering or (ii) to identify scenarios that enable the monitoring of associated non-linear processes, including the prospect of tracking the spatio-temporal excitation transport in extended quantum systems.
Here we present a systematic exploratory theoretical analysis to address the above two central questions. We develop a consistent semi-classical model for the electromagnetic field propagation in the presence of quantum coherent bound state dynamics and apply the approach to the benchmark scenario of XUV scattering from Helium nanodroplets. Specifically, we solve Maxwell's equations in the time domain using a quantum mechanical local polarization response treated with a few-level density matrix model. The parameters of the model are matched to the linear optical properties, which we describe by a Drude-Lorentz reference model, ensuring convergence to the correct linear limit for low intensity. Our systematic analysis of the scattering from a sub-µm Helium droplet for resonant excitation of the 1s 2 to 1s2p transition shows severe departures from the linear response behavior for laser parameters that are feasible with available laser technology. Our theory predicts substantial inelastic wide-angle scattering features that are reminiscent of the Mollow triplet [21] in atomic systems and are shown to reflect inhomogenous Rabi-cycling and excitation waves.
The remainder of this paper is structured as follows. In section 2 we introduce the two employed local polarization models and their self-consistent solution along with the electromagnetic field propagation. Section 3 contains a linear-response benchmark calculation to validate the few-level description and presents the main results for non-linear excitation. Finally, we summarize the main conclusions and provide an outlook for future work in section 4.

Methods
The basis of our semi-classical scattering model is the explicit propagation of the electromagnetic field according to Maxwell's curl equations in the absence of free charges and for non-magnetic materialṡ Here E and B are the spatially averaged electric and magnetic fields (continuum picture), ε 0 and µ 0 are the permittivity and permeability of free space, and P is the polarization. The latter has to be propagated self-consistently with the fields, contains the complete medium response, and is non-zero only inside the target. The equations are solved in 3D for a rectangular arena employing the Finite-Differences-Time-Domain method (FDTD) [22] and using Yee staggering of all discretized fields [23].
Two models are used to describe the time-dependent electronic polarization response, i.e. (i) a standard classical linear response model and (ii) a quantum model that is based on a few-level density matrix description of the electron dynamics.

Drude-Lorentz model
As the standard linear response reference we employ the Drude-Lorentz model (DLM) that mimics the local dispersive medium response via a set of damped harmonic oscillators. For a single oscillator per atom and an isotropic response, the atomic electric dipole p(t) is described by the classical equation of motion where E(t) is the electric field at the location of the oscillator, ω 0 and ν are the oscillator's eigenfrequency and the associated collision frequency, q e and m e are the electron's charge and mass, and f is the oscillator strength. For an atomic number density n 0 the associated polarization P(ω) = n 0 p(ω) in the Fourier domain ( ∂ ∂t → −iω) becomes where ε0me is the plasma frequency. Although we do not consider free charges, the plasma frequency can be used as a measure for the density to preserve the analogy to the Drude formula. The resulting relative permittivity ε r (ω) = 1 + χ(ω) finally establishes the link to the refractive index where δ and β describe its real and imaginary parts. A desired spectral dispersion can now be mimicked by a certain set of oscillator parameters (ω 0 , ν, f ). The generalization to multiple fractional oscillators (multiple resonances) is straightforward by representing the susceptibility as a sum χ(ω) = i χ i (ω), where each term reflects the response of an individual fractional oscillator according to equation 4. Note that this widely used superposition of uncoupled oscillators provides a valid description of the time-domain polarization only in the limit of linear response. For a generalization to a non-linear Lorentz model see references [24,25].
In this work we consider a single oscillator in the DLM and solve the equation of motion (equation 3) for the field-driven local dipole response in each cell self-consistently with the field evolution in equations 1-2, where the local polarization velocityṖ establishes the feedback to the field dynamics. This approach corresponds to the well-known auxiliary differential equation (ADE) method for the description of dispersion [22].

Few-level density matrix model
For the quantum mechanical description of the local response including coherence and non-linearity we use a few-level density matrix (FLDM) model that is sketched in figure 1.
The active Hilbert space is constructed from a ground state |0⟩ with energy E 0 that reflects He 1s 2 and the threefold degenerate first excited 1s2p-states (|1⟩, |2⟩, |3⟩) with energies E 1 = E 2 = E 3 . Within the atomic orbital representation, the dipole selection rules associated with the 1 s-2p transition lead to a particularly simple structure of the matrix elements µ ij = ⟨i|μ|j⟩ of the dipole operator in the Schrödinger picturê where q e is the charge of the electron, µ is the magnitude of the (isotropic) transition dipole moment and e i with (i = x, y, z) are the Cartesian unit vectors. The associated expectation value of the dipole for a given density operatorρ follows as The system Hamiltonian contains the unperturbed atomic part (withĤ (0) |i⟩ = E i |i⟩) and the interaction term in dipole approximation according tô Here |0⟩ is the ground state and |1⟩, |2⟩, and |3⟩ are the three degenerate excited state orbitals (2p). Solid arrows indicate the coupling to the respective field components and dashed arrows indicate relaxation (γr) and decoherence (γ d ).
The time evolution of the density matrix is evaluated in the interaction picture to remove the trivial unperturbed dynamics from the propagation. Matrix elements of a given operatorÂ in the Schrödinger picture are transformed into the interaction picture representation (henceforward indicated with a tilde) viã with the transition frequency ω mn = Em−En ℏ . In the evolution described by the generalized von-Neumann equationρ the last phenomenological dissipation term has been added to describe decoherence and relaxation towards the equilibrium stateρ eq. . The latter we consider to reflect the fully occupied ground state. In the absence of interaction, the dynamics is determined only by relaxation and decoherence. We assume symmetric relaxation and decoherence behavior between ground and excited states with where the rates for decoherence and relaxation are denoted as γ d and γ r , respectively. It should be emphasized that our assumptions of an isotropic response and equivalent excited states requires vanishing decoherence among excited states. Further we assume that relaxation, e.g. due to radiative decay, is negligible on the time scales of the considered fs scattering dynamics and set γ r = 0. We like to note that in principle the inclusion of additional relaxation channels is straightforward. For the considered 1s2p excited state, the only additional possible channel would be the 1s2p to 1s2s relaxation [26]. As our scenario is focused on the coherent dynamics, its impact is considered negligible as long as the relaxation time is large compared to the dephasing time.
The response properties of the FLDM model are determined by the excitation energy E 1 − E 0 = ℏω 10 , the magnitude of the transition dipole µ, and the decoherence rate γ d . In analogy to the DLM description, an individual local few-level system is propagated in each FDTD cell using the corresponding electric field to solve the evolution for the polarization.
The resulting polarization velocityṖ serves as the feedback for the field propagation, in close analogy to the DLM case in equation 6.

Definition and matching of the model parameters
Both the DLM and the FLDM model describe the response properties via respective sets of oscillator parameters, where the DLM parameters are directly linked to the refractive index and the dispersion of the model system. In this work the response near the first excitation of Helium (1s 2 → 1s2p) will be modelled with a single oscillator, which requires to fix the three DLM parameters ω 0 , ν, f. In order to describe liquid Helium we use parameters (ℏω 0 = 21.61 eV, ℏν = 0.43 eV, and f = 0.49) that where fitted to experimental data in references [27,28]. The resulting dispersion described by the DLM for an atomic density of n 0 = 0.022 Å −3 is compared to the experimental data from references [27,28] in figure 2. At this point, all parameters of the DLM description are fixed. In order to determine the FLDM model parameters (γ d , ℏω 01 , µ), we require that the resulting linear optical properties must match the DLM description. Therefore the weak-field limit of the FLDM response has to be determined, which, because of isotropy of the linear response in our model, need to be considered only along one axis. e.g. the x-axis. The evolution of the expectation value of the dipole acceleration in the weak field limit (ρ 00 = 1 and ρ jj ≈ 0 for j = 1, 2, 3) can then be written in a similar way as for the classical oscillator in equation 3⟨ The expectation values of acceleration, velocity and magnitude of the dipole obviously coincide with the classical result if the prefactors are matched as indicated. Solving the matching conditions for the FLDM parameters yields where e is the elementary charge. The FLDM parameters (ℏγ d = 0.215 eV, ℏω 10 = 21.61 eV, µ = 0.2939 eÅ = 1.41 D) are now fully matched to above DLM parameters.

Results and discussion
To investigate the impact of the non-linear local polarization response on the scattering process and the final diffraction pattern for a realistic scenario, we considered a spherical Helium droplet with radius R = 300 nm as a model system. The sphere is centered in the computational domain that is an equidistant, cubic grid of 251 × 251 × 251 cells with a spatial resolution of ∆x = 2.8 nm, see figure 3. The resonant incident pulse with central photon energy of ℏω inc = 21.61 eV and flat phase front is linearly polarized along the x-axis, Figure 3. Sketch of the scattering scenario described in the simulations. A spherical helium droplet with radius R = 300 nm is placed in a cubic arena (left). The incident pulse with flat spatial intensity profile and Gaussian temporal profile propagates along the z-axis and is polarized along the x-axis. The equivalent source method, sampled on a control surface (green cube), is used to determine spectrally resolved far-field scattering images (as indicated). Radial scattering intensity profiles (right) along the qx axis for qy = 0 for low incident field (linear response limit) are calculated using the FLDM method (blue line) and DLM (red line).
propagates along the z-axis and is described with a Gaussian temporal envelope with 8 fs duration (FWHM). Such a short duration is required to study coherent dynamics on the timescale of the dephasing time.

Spectrally selective scattering images and linear response benchmark
The incidence pulse is injected into the FDTD arena using the total-field-scattered-field scheme. The choice of the DLM or FLDM models determines only the engine that is used for the integration of the polarization response. The scattered fields are sampled on a cubic control surface (green) and decomposed into spectral components on-the-fly via a discrete Fourier transform. The spectral amplitudes and phases of the electric and magnetic fields on the control surface are used to construct electric and magnetic surface currents that represent an equivalent source [22]. From the latter the scattering in the far-field is sampled for each spectral component, which yields spectrally-resolved scattering images. Note that this method allows the description of the scattering in full solid angle. Summation over individual spectral contributions in real space results in the final predicted scattering images, see figure 3. To verify consistence of our description for both polarization models, we calculated scattering images with the DLM and FLDM simulations for low pulse intensity. The comparison of the respective (spectrally integrated) scattering profiles in figure 3 documents the agreement in the linear response limit and validates our method. The scattering images show the typical ring structure for spherical droplets [11].

Population and electric field dynamics for non-linear excitaion
Pronounced non-linear features appear in the prediction of the FLDM model as the intensity of the incident pulse is increased. In the following we discuss as a representative case the excitation with pulse intensity The resulting snapshots of the ground state depletion in the propagation-polarization plane (figure 4(a)-(d)) display curved regions of similar depletion level that form excitation (or depletion) slabs whose structure approximately reflect the droplet curvature. The evolution of the depletion maps reveals traveling excitation waves (from left to right) with a focus-like convergence towards the back of the droplet. The corresponding snapshots of the local polarization (figure 4(e)-(h)) show that the phase fronts of the dipole excitation within the individual slabs are flat and do not reflect the curvature of the excitation slabs. However, the slab structure remains visible in the polarization signal in terms of the nodes of the polarization amplitude (white areas) that are associated with maxima and minima of the depletion (individual level population).
The detailed evolution of the level population in the droplet center in figure 4(i) shows distinct signatures of Rabi-flopping between the ground and excited states |0⟩ and |1⟩ due to the coupling to the x-component of the electric field. Several Rabi-cycles are predicted before the coherent state dynamics ceases near the end of the pulse. Note that decoherence leads to a steady decrease of the amplitude of the oscillating population, underlining the necessity of short pulses for the realization of substantial Rabi cycling. The associated induced polarization in the droplet center displayed in figure 4(j) shows a clear fingerprint of the Rabi flopping and the associated dressed states, i.e. a signal modulation with the period of the Rabi cycle, with The detailed time-evolution of the population (i) and polarization (j) corresponds to the droplet center. The gray shaded area indicates the incident free-space intensity envelope at the center plane of the droplet. The spectrum of the polarization signal from panel (j) is shown in panel (k), documenting the spectral splitting and the formation of satellites that are up-/down shifted from the central frequency ω inc by roughly the Rabi frequency (as indicated). minimal amplitude for dominant population of only one state and maximal amplitude for equal population. The underlying field induced level splitting is expressed by the redistribution of the spectral polarization signal in two dominant satellites that are up-/downshifted from the central frequency of the incidence beam by the Rabi frequency Ω = µE 0 where E 0 is the laser electric field amplitude, see figure 4(k). The saturation effects, the decoherence, as well as the Rabi-cycling can be expected to cause substantial changes in the scattering images compared to the linear response result. In any case, at least for situations similar to the considered resonant scenario, the coherent state dynamics must obviously be taken into account already for intensities approaching 10 14 W cm −2 , which are reachable in experiments already.

Nonlinear coherent diffractive imaging
As the next step we investigate the structure and evolution of the scattering images in the regime of non-linear excitation. Figures 5(a)-(c) show three representative examples of spectrally integrated scattering patterns for different excitation intensities. Three main trends can be identified from the comparison. First, the total intensity of the scattering images shows a saturation behavior, as documented by the transition from a constant to a decreasing scattering cross section (∝ 1/I) at high intensity, see figure 5(k). Note that the switching depends on the system size. Second, the spatial structure of the fringes becomes blurred with increasing intensity, exhibits a shift towards larger angles, and shows a modified decay behavior, with reduced angular decay at high excitation intensity, see figure 5(j). Third, the spectra of the scattered fields widen substantially when compared with the linear-response results, see evolution of the full spectra (black line) in figures 5(d)-(f). The high intensity result shows clearly the emergence of an inelastically scattered radiation component in terms of satellites at ω inc ± Ω that is associated with the non-linear scattering process.
The saturation effect, as the most obvious feature, can be explained at least qualitatively from the fact that the polarization amplitude in the few-level system has an upper bound (∝ µ). As a result, also the scattering from a pulse with finite duration must have an upper limit, impeding a linear scaling with the incidence pulse intensity at some point. The detailed modification of the fringe structure is less obvious. The dressed states created by the Rabi-cycling can be considered as an induced transparency, which will most strongly affect the droplet surface such that the effective sphere radius shrinks. This argument is consistent with the slight increase of the fringe spacing in figure 5(j). For the blurring of the fringe features itself we consider two possible contributions. The first is the increased spectral width of the scattered fields. For linear scattering, spectral broadening reduces the fringe contrast in scattering patterns [13] due to the mismatch in the fringe separations (Mie-scattering) associated with the different spectral contributions. However, another effect which we expect to be even more important for our case is connected with the structure of the excitation slabs.
In order to motivate the slab structure effect we consider a slicing of the droplet into curved discs that resemble the typical thickness and curvature of the excitation slabs. The passage of the train of excitation waves through the spatially fixed slices leads to the emission of a Rabi splitted spectrum. As individual slabs have different lateral sizes and thus create different fringe spacing in the diffraction pattern, the fringe structure in the full image is expected to blur if adjacent slices emit incoherently. This argument, however, would imply that the inelastically scattered radiation would show a weaker angular decay than the elastic (volume) contribution. The underlying logic is that a disc-like object shows a weaker angular decay in the scattering profile than a sphere. In the limit of the Born approximation the profiles of a disc and a sphere decay with q −3 and q −4 , respectively. Though these values are modified in the presence of absorption, the trend of weaker angular decay for disc-like scatterers should remain. In fact, the frequency resolved, integrated scattering signal in figures 5(e),(f) documents that the relative contribution of the inelastic signal is substantially higher for large scattering angles. Furthermore, the inelastic contribution shows a spectral asymmetry, which, however, also depends on scattering angle. In particular, in figure 5(f), the downshifted branch is dominant over the up-shifted feature for large scattering angles (blue shaded area), while the asymmetry is reversed for the small angle signal (red shaded area). A detailed analysis of this behavior, which also exhibits an intensity dependence, is beyond the scope of the current work.
Finally, if the inelastic signals of the scattered fields can be associated predominantly with the Rabi-cycling in the excitation-slabs, the spectrally selective analysis of the scattering images may provide rich information about the underlying spatio-temporal excitation dynamics. Figures 5(g)-(i) display a corresponding analysis from the high intensity case in figures 5(c),(f) and show the scattering images associated with only the scattered radiation at the central frequency and for frequencies around the two Rabi satellites, respectively. The deviation of the individual patterns, including the elastic one in figure 5(h), from the linear response result highlights the requirement to include the non-linear response in the analysis of scattering data. Moreover, the fundamental structural differences between the individual patterns show that the non-linear response leaves a pronounced spectral fingerprint in the scattering patterns, opening a promising perspective for time-resolved characterization of the excitation dynamics.

Conclusions
In summary, we developed a density matrix-based description of coherent diffractive imaging including quantum coherent bound-state dynamics and employed it to the scenario of resonance scattering from Helium nanodroplets. Although our model was used here for the describtion of a single resonance, it will be straightforward to describe more complex systems by adding additional active states. Our model is constructed such that it converges to the correct linear-response result for given optical parameters. For the investigated Helium droplets, our prediction at high intensity reveals a substantial impact of the non-linear quantum response. We find spatio-temporal excitation waves, Rabi-cycling, inelastic features in the scattered radiation, and pronounced deviations of the scattering patterns from the linear response case. Even though the presented application is based on strong approximations (single resonance, ionization neglected, correlated decay processes neglected), our findings show that quantum effects cannot be neglected for scenarios that are already experimentally in reach. The specific parameters of a corresponding experiment strongly depend on the particular objective. If coherent dynamics and Rabi satellites are of major interest, the laser pulse should be short enough to compete with the dephasing time but long enough to realize several Rabi cycles. If the goal is the mere observation of non-linear features, the pulse duration is less critical and just the intensity has to be high enough to achieve substantial excited state population. Finally, the rich features in the non-linear scattering reveal a promising potential for the characterization of complex quantum dynamics at the nanoscale via quantum coherent diffractive imaging and justify further in-depth analysis.