Dissipation of electron-beam-driven plasma wakes

Metre-scale plasma wakefield accelerators have imparted energy gain approaching 10 gigaelectronvolts to single nano-Coulomb electron bunches. To reach useful average currents, however, the enormous energy density that the driver deposits into the wake must be removed efficiently between shots. Yet mechanisms by which wakes dissipate their energy into surrounding plasma remain poorly understood. Here, we report picosecond-time-resolved, grazing-angle optical shadowgraphic measurements and large-scale particle-in-cell simulations of ion channels emerging from broken wakes that electron bunches from the SLAC linac generate in tenuous lithium plasma. Measurements show the channel boundary expands radially at 1 million metres-per-second for over a nanosecond. Simulations show that ions and electrons that the original wake propels outward, carrying 90 percent of its energy, drive this expansion by impact-ionizing surrounding neutral lithium. The results provide a basis for understanding global thermodynamics of multi-GeV plasma accelerators, which underlie their viability for applications demanding high average beam current.


Introduction
Localized deposition of concentrated energy into plasmas by particle bunches or laser pulses underlies fast ignition of laser fusion 1 , formation of plasma waveguides 2 , and wakefield acceleration of electrons and positrons 3 . Subrelativistic particles or laser pulses deposit energy into dense uniform plasmas primarily via collisions 2 , analogous to ohmically heating a resistor, but collisions become inefficient for tenuous plasma and/or relativistic excitation. Relativistic particle bunches (or laser pulses), on the other hand, excite tenuous plasma by creating electron density structures, such as channels 4 or Langmuir waves 3 , via Coulomb (or ponderomotive) forces, analogous to charging a capacitor. Previous experiments in plasmas of millimeter (mm) length and atmospheric electron density (n e ∼ 10 19 cm −3 ) have documented the initial (first few ps) stage of transferring energy stored initially in such structures into long-term ion motion 4,5 .
The emergence of quasi-monoenergetic multi-GeV plasma-wakefield accelerators, in which relativistic particle bunches 6,7 (or laser pulses 8 ) deposit energy into strongly nonlinear wakes in plasma of density n e ∼ 10 17 cm −3 , however, raises the question of energy dissipation with renewed urgency. Multi-GeV plasma accelerators require such low n e in order to avoid depleting the driver's energy before it propagates ∼ 1 m, the length of plasma wake needed to accelerate electrons or positrons to multi-GeV energy 3 . Heat removal by flowing the medium supersonically transverse to the driver propagation direction, which is routine for mm-long, 10 µm-wide, high-n e , gas-jet-based MeV plasma accelerators 3,5 , is not feasible for metre-long, 100 µm-wide, low-n e , multi-GeV accelerators, which require stationary vessels to confine their long, tenuous plasmas [6][7][8] . Such accelerators will therefore require new strategies for managing deposited heat, based on quantitative understanding of energy dissipation. Such understanding is fundamental to achieving luminosities large enough to observe rare processes at the energy frontier. Accelerators achieve high luminosity by delivering focused high-charge, highenergy beams at high repetition rate. Planned conventional machines such as the International Linear Collider 9 or the Compact Linear Collider 10 , for example, are designed to deliver ∼ 10 megawatt (MW) to the interaction point. Current conceptual designs for plasma-based accelerators aim to achieve comparable luminosity by accelerating ∼ nC bunches at ∼ 10 kHz repetition rate (i.e. ∼ 100 µs interbunch spacing) 11,12 . This need has spurred world-wide efforts to develop petawatt-peak-power drive lasers that can operate at multi-kHz repetition rates 13 . However, the problem of dissipating, or otherwise managing, deposited power of order tens of kW per meter over ∼ 100 m of plasma has received comparatively little attention. The present study aims to understand the fundamental processes by which plasma accelerator structures at n e ∼ 10 17 cm −3 dissipate their energy into surrounding plasma, and to evaluate their global energy budget over a nanosecond (ns) time scale. It thus builds a thermodynamic foundation on which future engineering solutions of the heat management problem can be based.
To produce high-quality bunches, plasma accelerators of any n e usually operate in a strongly nonlinear regime in which the driver blows out a "bubble"-like accelerating cavity devoid of plasma electrons in its immediate wake 3 . In such structures the energy density |E| 2 /2 0 of internal wake fields E approaches the rest energy density n e mc 2 of plasma electrons 14 . Simulations predict that such electron wakes can spawn ion wakes of unique structure and dynamics [14][15][16] , early stages of which were observed 5 at high n e . However, no experiments at any n e have yet explored how the enormous energy density stored in a nonlinear wake redistributes over nanoseconds amongst accelerated electrons, undirected hot electrons, freely streaming ions, radiation, electrostatic fields, and ionization of surrounding gas, as well as collective ion motion. Understanding this complex process at n e ∼ 10 17 cm −3 demands experiments with precisely-characterized multi-GeV drive bunches (or petawatt laser pulses), probes that track particle and energy flow over millimeters, and simulation of multifarious plasma processes over nanoseconds. The scale and complexity of the problem rival those of other energy transport problems involving tenuous plasma, such as heating of the solar corona 17 and acceleration of energetic cosmic rays 18 .
Here we present ps-time-resolved optical shadowgraphic measurements of meter-length ion channels that emerge from broken, highly nonlinear plasma wakes. Quantitative observation of the plasma column expansion over nanoseconds to microseconds serves as a calorimeter that determines the fraction of the initial wake energy that the plasma column retains after the wake breaks. Modeling of the energy transfer from initial wake into outward ion motion, and benchmarking of simulation results against measurements, quantifies energy retention in the plasma column and elucidates the physical mechanisms that drive its expansion.

Results
Generation of nonlinear wakes. Experiments were carried out at the SLAC Facility for Advanced aCcelerator Experimental Tests (FACET) 19 . The first 2 km of the SLAC linac delivered drive electron bunches (ebunches) of energy 20 GeV, charge 2.0 nC, rms radius σ r = 30 µm, length σ z = 55 µm to the entrance of the interaction region. Unlike small-scale plasma wakefield accelerator (PWFA) experiments, in which parameters of drive e-bunches delivered from a tabletop laser wakefield accelerator (LWFA) could only be estimated 5 , here e-bunches from SLAC were characterized with high precision (see Methods). This is important for simulating subsequent plasma dynamics accurately over a ns time scale. Drive bunches entered a 150 cm-long column of Li vapor, in which a 120 cm-long region of uniform atomic density n a = 8 × 10 16 cm −3 was centered between 15 cmlong entrance and exit density ramps. This vapor was generated and contained within a cylindrical heat-pipe oven 6,20 of 2.5 cm radius (see Methods), and provided the medium for plasma formation and wake excitation.
Particle-in-cell (PIC) simulations discussed below showed that the SLAC e-bunches singly field-ionized the Li vapor over its entire length out to initial radius ∼ 40 µm and electron density n e = n a , and drove a strongly nonlinear electron density wake consisting of a train of nearly fully blown-out cavities of radius ∼ 20 µm, propagating at ∼ c along the axis of the resulting plasma. The cavities enclosed accelerating fields of magnitude E ∼ mω p c/e ∼ 1 GV/cm, where ω p = [n e e 2 / 0 m] 1/2 is the electron plasma frequency. Although we have the capability at FACET to pre-ionize the lithium vapor along the entire drive beam path out to r ∼ 500 µm, 21 here we retained the Li vapor blanket as an in-situ medium for recording energy transport out of the directly-excited wake. As discussed below, outwardly streaming ions and electrons carry away most of the wake's energy, and expend < 1% of it ionizing Li within the first ∼ 1 ns. Nevertheless, ionization induces a large change in the vapor's refractive index η, enabling us to detect the ionization front with ps time and ∼ 20 µm space resolution, without perturbing overall energy transport dynamics significantly.
When SLAC delivered 2 nC, 20 GeV (i.e. 40 J total energy) drive bunches to the interaction region at repetition rate 1 Hz or less, no heat-pipe temperature rise, nor other time-dependent behavior, was observed upon turning on the beam. Downstream measurements of spent driver energy (see Supplementary Material, Fig. S2) show that ∼ 60% of electrons in each bunch contribute to forming a plasma wake, each losing 11% of its energy in so doingi.e. each bunch deposits 7% of its total incident energy, or 2.6 J, into the 1.2 m-long plasma wake (2.2 J/m). Thus at 1 Hz, the plasma acquires energy at 2.6 W, which is only 0.2% of typical oven heater power (1300 W). When, on the other hand, we increased repetition rate to ∼ 10 Hz, oven temperature typically rose tens of degrees within minutes, even though this repetition rate is ∼ 1000× lower than current design projections 11 . Here we report data acquired at 1 Hz, to avoid drifts in oven temperature during extended data acquisition runs. No witness e-bunch was injected with the driver. Thus plasma wave energy was dissipated entirely within the medium.
Measurement of plasma expansion. To diagnose the evolving radial profile of the plasma column, a collimated probe laser pulse (wavelength λ = 0.8 µm, duration τ = 0.1 ps, transverse width 0.5 cm FWHM) that was electronically synchronized with the drive e-bunch with ∼ 0.1 ps jitter entered one end of the heat pipe 0.8 cm offcenter. It then propagated through a 100-cm-long central section of the Li column at angle θ = 8 mrad to the e-bunch propagation direction [see Fig. 1(a),(b)] at time delay −1 ps < ∆t < 10 µs after the e-bunch, before exiting the other end on the opposite side of center. By using this grazing θ, the largest that the long narrow heat-pipe oven allowed without clipping the probe pulse, we probed the tenuous plasma profile ∼ 100× more sensitively than with a transverse (θ = π/2 rad) probe, because of the long interaction length, at the mild cost of averaging longitudinal (z) density variations n e (z). Moreover, as discussed below, with this geometry we maximized sensitivity to the region of interest -the expanding profile's advancing outer edge -at the cost of insensitivity to its internal structures. Use of an e-bunch driver eliminated strong forward-directed supercontinuum that a laser-driven nonlinear wake generates 22 , which would saturate probe detectors in this near co-propagating pumpprobe geometry.
A lens imaged the portion of the transmitted probe pulse that it collected within its f /40 cone from a vacuum object plane at longitudinal position z = 75 cm near the center of the heat-pipe [labeled "O" in Fig. 1(a)] to a charge-coupled device (CCD) camera. Here z = 0 refers to the foot of the entrance density ramp. Images had ∼ 1 cm depth of field and unity magnification. Unavoidable obstructions blocked ∼ 1/4 of the circular lens aperture [see grey area, Fig. 1(b)], which influenced some nonessential details of images, as discussed below. When the probe arrived before the e-bunch (∆t < 0), only the incident probe pulse profile was observed. At ∆t > 0, an approximately parabola-shaped shadow outlined with alternating bright and dark fringes, resulting from probe refraction and diffraction from the e-bunchexcited plasma column, was observed. Fig. 1(b) shows schematically how the shadow/fringe pattern evolves for fixed ∆t as the probe propagates through the near field of the plasma column. Fig. 1(c) shows how it evolves in the far field as ∆t varies from 100 to 1200 ps. The parabolic shadow expands at ∼ 10 6 m/s, eventually exceeding the cameras field of view for ∆t > 1200 ps. Figs. 1(d)-(f) show calculated images for 3 simulations of plasma evolution, discussed below and in Methods. One horizontal position in these patterns (highlighted by vertical white dashed line in Fig. 1(c)-(f)) corresponds to "O". Fringed regions straddling this plane are out-of-focus images of upstream (z < 75 cm) and downstream (z > 75 cm) slices of the column. They represent images of near-field diffraction patterns of this obliquely illuminated column at each z 23 . At small ∆t, the vertex of the parabola corresponds to z ≈ 75 cm [see 100 ps images in Fig. 1(c)-(f)]. As the plasma column widens, the vertex shifts downstream of z = 75 cm by distance ∆z, and a portion of the parabola of radius r O shifts to z = 75 cm [see 900 ps images in Figs. 1(c)]. This occurs because plasma lensing of the probe at z > 75 cm contributes increasingly to image formation as the column widens.  Fig. 2(a) shows a half cross-section of the plasma column's refractive index profile η(r) at a representative ∆t, with contours (black circles) in ∆η = 10 −5 increments superposed. The bottom half of Fig. 2(a) shows the corresponding intensity profile of the probe pulse at the same z-plane. In this picture, the probe emerges along the page normal, while the plasma column axis is tilted 8 mrad from left (behind page) to right (in front of page). A parabolic shadow (left) develops as probe light refracts away from columns axis, after penetrating to the η − 1 = 10 −5 contour in the plane of incidence. Fringes (right) approximately parallel to the shadow boundary develop as refracted probe light interferes with un-deflected incident probe light. This probe profile is, however, not directly observed. Rather, a lens relays it to a CCD. For an ideal unobstructed lens and cylindrically symmetric plasma column, the detector would record the profile shown in Fig. 2(b). Because the probe passes through the remainder of the plasma column, the image is distorted in two ways from the probe's in-situ shape. First, a nephroid-shaped cusp singularity forms in the image of the parabola's vertex, similar to bright optical caustics that form within the shadow of an obliquely illuminated drinking glass 24 . Second, a second set of interference fringes approximately orthogonal to the shadow boundary develops outside the shadow.
The first of these features proved sensitive to the above-mentioned partial blockage of the lens aperture. Fig. 2(c) shows the reshaped shadow vertex that results when probe propagation calculations take this blockage into account. The second of these features proved sensitive to slight deviations of the simulated column from cylindrical symmetry, and to small longitudinal non-uniformities. Thus these two features were seldom observed in actual images [see Fig. 1(c)]. Nevertheless, the vertex shift ∆z, along with shadow radii r O at "O" and r B at another selected longitudinal location "B", as marked in Fig. 1(c), were consistently observed. Orange data points and grey uncertainty ranges in Fig. 2 show evolving values r O (∆t), r B (∆t) and ∆z(∆t), respectively, for 0 < ∆t < 1200 ps, where "B" here corresponds to z = 55 cm. Black and orange curves in Fig. 2(g), (h) show single-shot and 30-shot averaged lineouts, respectively, of shadows at "O" (top) and "B" (bottom) for ∆t = 100 ps (g) and 1200 ps (h). Fringes outside of, and parallel to, the shadow boundary were observed frequently in both single-shot [ Fig. 2(g), bottom panel, black curve] and multi-shot averaged [ Fig. 1(c), top four panels] data. Probe propagation simulations showed that fringe spacing was sensitive to θ, but agreed well with observed fringe spacing for θ = 8.0 ± 0.1 mrad, thus corroborating the independently measured probe angle. However, observed fringe contrast is invariably lower than calculated, due to imperfect symmetry of the plasma column. In addition, fringes sometimes wash out upon multi-shot averaging, because of shot-to-shot fluctuations in θ and in drive bunch intensity. Thus, r O (∆t), r B (∆t) and ∆z(∆t) constituted the most robust observables for quantitative comparison with simulations.
Qualitative plasma expansion mechanisms. Results in Fig. 2(d)-(e) show that empirical radii r O (∆t) and r B (∆t) grow on average at 1.4 × 10 6 m/s over 1.3 ns, and even accelerate slightly during this interval. Probe propagation simulations show that these radii are approximately twice the plasma column radius r p , here defined as the radius at which n e (r p ) ≈ 0.2n a , implying that r p expands at near-constant velocity v p ≈ 0.7 × 10 6 m/s. Such kinetics rule out expansion driven by electron heat or radiative transport 25 , since heat front velocity would decrease rapidly with time 4 . They also cannot be explained by a radial shock wave at the ion acoustic velocity, which is initially only ∼ 10 4 m/s for our conditions, and would also decrease within ∼ 1 ns as electron temperature cools 2 . Rather, the observed near-constant v p must be attributed to an ionization front driven by charged particles, particularly high-momentum Li ions, streaming freely into surrounding Li vapor 4 . Li ions propagating at the observed v p , for example, would have energy E i ∼ 20 keV. These could be produced if average outward radial electrostatic fields of order 0.01 < E r < 0.05 GV/cm acted on the ions over radial distance of order 200 > r > 40 µm (or time interval 50 > τ > 10 ps) in the collapsed electron wake following its excitation. Such ions have mean free paths of several mm in Li vapor of n a = 0.8 × 10 17 cm −3 , experience only small angle elastic scattering, and lose energy to neutral atoms primarily via impact excitation and ionization, which entails loss per impact only up to the first ionization energy (∼ 5.4 eV) of Li. These characteristics are consistent with near-constant-velocity expansion over a ns time scale. In addition, the ions escort electrons, which maintain charge quasi-neutrality and assist in exciting and ionizing Li atoms.
Simulations of plasma expansion. To understand plasma expansion quantitatively, we carried out PIC simulations of Li plasma dynamics out to ∆t ∼ 1.3 ns. We modeled ionization, self-focusing of the drive bunch, electron wake excitation, and early (∆t ≤ 40 ps) electron wake and ion dynamics in a fully self-consistent manner using OSIRIS 26 in cylindrical geometry (see Methods for details). These simulations modeled radial acceleration of ions and electrons by electrostatic fields of the collapsing electron wake with high space and time resolution, before these particles began to interact significantly with surrounding Li gas. OSIRIS simulations tracked self-consistent driver and plasma evolution for 15 cm of propagation through the gas density up-ramp and an additional 16.9 cm into the 120-cm-long density plateau (i.e. up to z = 31.9 cm).
To simulate long-term plasma evolution, we input the compressed e-bunch and plasma profiles from the selfconsistent OSIRIS simulation output at z = 31.9 cm into the quasi-static, axisymmetric LCODE 27 as initial conditions. The quasi-static approximation reduces dimensionality by one, making LCODE more time-efficient for simulating larger ∆t, at the cost of neglecting selfconsistent evolution of drive bunch and plasma for z > 31.9 cm. Moreover, impact ionization of neutral gas by outwardly streaming electrons and ions comes into play on this time scale, and in fact dominates plasma expansion for ∆t > 100 ps. Accordingly, we included the most important elementary collisional processes -elastic binary collisions 28,29 and impact ionization of neutral lithium atoms by electrons 30 and fast lithium ions 31into LCODE. For electron impact ionization, we included single-step ionization from the neutral lithium ground state and two-step ionization 32,33 via 2P, 3S and 3P excited states (see Methods). Fig. 3 shows OSIRIS simulation results. The electron bunch, initially of bi-Gaussian r, z profile [see Fig. 3(a)], ionized Li according to the Ammosov-Delone-Krainov (ADK) model 34 . Fig. 3(b) shows the driver (blue) and plasma (orange) after the former propagated to z = 31.9 cm. The bunch's leading edge was not dense enough to ionize Li, and thus propagated as in vacuum. Ionization started in the denser center of the bunch, which then drove a nonlinear plasma wake, which in turn compressed the beam waist radially. In Fig. 3(b), the bunch's trailing edge has compressed to ∼ 100 times its initial density. This enabled it to singly-ionize Li gas completely out to r ≈ 40 µm, and partially out to r ≈ 65 µm, and to drive fully blown-out plasma bubbles of radius λ p = 20 µm. Fig. 3(d) shows the radial electric field E r (r, z) remaining at ∆t = 1 ps. Typical of the aftermath of a nonlinear wake, this field is still non-zero 16 . Because of their large mass, plasma ions initially respond to E r (r) , i.e. E r (r, z) averaged longitudinally over ∼ 6λ p , shown in Fig. 3(c). Since E r (r) switches sign at r ≈ 30 µm, it attracts ions initially (i.e. before wave-breaking occurs) at r < 30 µm toward the axis, while pushing ions initially at r > 30 µm outward. The driver also expels a fraction of plasma electrons into surrounding neutral gas, leaving net positive charge in the plasma that further propels outward ion motion. The resulting ion density structure n Li + (r, z), seen in Fig. 3(e,f) at ∆t = 40 ps before (f) and after (e) longitudinal averaging, has a peak on axis. In addition, the outermost ions diffuse outward, moving from ∼ 65 µm [dashed lines in Figs. 3(e,f)] to ∼ 100 µm. These features contrast with observations of Ref. 5 , where a region almost void of plasma formed around the onaxis peak, surrounded by a thin cylindrically symmetric plasma sheath. This difference arises because in our case the initial plasma radius is less than λ p ≈ 120 µm, whereas in Ref. 5 the 75× denser plasma was pre-ionized beyond λ p ≈ 14 µm. Short-term (∆t < 40 ps) evolution of the plasma thus differed from that investigated here. Fig. 3(g) shows the calculated probe diffraction pattern from a longitudinally uniform plasma column of the shape of Fig. 3(e) for our experimental geometry. It is very close to ∆t = 100 ps data [ Fig. 1(c)] and simulation [ Fig. 1(d)]. Indeed, our simulations predict negligible change in this pattern during the interval 40 < ∆t < 100 ps, since radially accelerated ions and electrons have not yet begun to impact-ionize surrounding Li neutrals substantially. Fig. 4 presents LCODE simulation results for ∆t > 100 ps. Fig. 4(a) shows how the plasma density profile n e (r, ∆t) evolves over the interval 100 < ∆t < 1200 ps in the absence of impact ionization processes (hereafter "simulation 1"). At ∆t = 100 ps, the dominant feature is a sharp axial electron density maximum n e (r < 0.01 mm). This is the electron counterpart of the ion density maximum n Li + (r < 0.01 mm) seen in Fig. 3(e)-(f) at ∆t = 40 ps, and maintains its quasi-neutrality. Over the ensuing 1100 ps, this axial peak drops in amplitude and broadens, driven by electrostatic forces. Nevertheless, n e (r > 40 µm) never exceeds 10 16 cm −3 . Fig. 4(b) shows corresponding refractive index profiles η(r, ∆t) (see Methods). The index profile does not change noticeably for r > 40 µm throughout the simulated interval. Since the probe pulse turning radius [η − 1 = 10 −5 , blue dashed line in Fig. 4(b)] occurs at r ≈ 50 µm in our geometry, the probe does not sense index changes occurring at r < 40 µm [shown in Fig. 4(b)]. Hence, simulation 1 predicts no change in probe signatures over the interval 100 < ∆t < 1200 ps, as shown in Fig. 1(d To capture this continuing long-term expansion, we included impact ionization, induced by energetic electrons and ions streaming radially outward from the directly e-beam-ionized plasma [ Fig. 3(b)]. Fig. 4(c) shows evolving n e (r, ∆t) profiles when these processes are restricted to single-step ionization from the neutral Li ground state (hereafter "simulation 2"). The earliest n e (r, ∆t = 100 ps) profile shown, and its corresponding η(r, ∆t = 100 ps) profile in Fig. 4(d), differ only slightly from their counterparts in simulation 1. By ∆t = 400 ps, however, impact ionization has begun to create substantial new plasma in the region 50 < r < 150 µm, which reaches Li vapor density (i.e. n e = n a = 8 × 10 16 cm −3 ) by ∆t = 1200 ps [see Fig. 4(c)]. At the microscopic level, the simulation shows that energetic electrons, which exit the original plasma first, ionize some neutral Li in this region directly, but inefficiently, since their density and ionization cross-section are small. A moving front of fast ions creates most of the new plasma. Once ions appear at a given location, more electrons come, including lower-energy plasma electrons with high impactionization cross-sections. Growth of this electron population triggers near-exponential plasma density growth.

(f)].
To capture these remaining features we included twostep electron impact ionization into LCODE (hereafter "simulation 3"). In these processes, an initial electron impact excited the 2S electron of a ground state Li atom to a 2P, 3S or 3P state. A subsequent impact within the natural lifetime of these states then ionized the atom. Including such processes increased new plasma production significantly, as n e (r, ∆t) plots in Fig. 4(e) show. Simultaneously they hastened growth of probe turning radius, as η(r, ∆t) plots in Fig. 4(f) show. The partial ion density distributions N i (r) in Fig. 5(a) (shades of red) show that e-impact of excited Li* atoms is, in fact, the dominant source of Li + ions at ∆t = 1200 ps. Fig. 5(b) shows, however, that two-step processes become dominant only at ∆t > ∼ 400 ps (see orange, red, green curves). Simulated probe images [ Fig. 1(f)] widen more than twice as rapidly as for simulation 2 [ Fig. 1(e)]. Moreover, a substantial vertex shift ∆z(1200 ps) ≈ 500 µm develops by the end of simulation 3. Simulated average growth of r O (∆t) and r B (∆t) agrees with observed average growth over the interval 100 < ∆t < 1300 ps [see Fig. 2(d), (e)]. Finally, simulated and observed probe image lineouts near the beginning [ Fig. 2(g)] and end [ Fig. 2(h)] of simulation 3 agree well in width and depth, despite discrepancies in fringe amplitude. Thus simulation 3 captures all qualitative features of the data, including the vertex shift ∆z(∆t) that simulation 2 missed, as well as some key quantitative features.
There are several plausible reasons for these discrepancies. First, although incident drive bunches were thoroughly characterized, properties of the bunch after its trailing part focused inside the plasma [ Fig. 3(b)] govern plasma expansion dynamics. Because plasma lensing is nonlinear, small errors in incident bunch properties can lead to large errors in focused bunch properties. For example, simulations assumed axi-symmetric drive bunches, whereas ∼ 10% asymmetries between σ x and σ y , and 10-fold differences between focusing functions β x and β y , were typically present at the plasma entrance, and could have led to asymmetric downstream focusing and plasma expansion. A second probe in an orthogonal plane would help to diagnose such expansion asymme-tries, if present. Similarly, deviations in the longitudinal bunch shape from Gaussian, which were not well characterized, sensitively influence the intra-bunch position at which ionization and self-focusing begin, and in turn the fraction of incident bunch charge that drives a nonlinear wake. This can also lead to significant discrepancies between observed and simulated expansion rates.
In the later part of simulation 3 (600 < ∆t < 1200 ps), the radial slope ∂η/∂r of advancing refractive index profiles at the turning point radius shrinks rapidly [see Fig. 4(f)]. As a result, simulated radii r B (∆t > ∼ 900 ps) in Fig. 2(e) become sensitive to small perturbations in η(r) profiles, and by extension to small deviations in the radial profile of the focused drive bunch. Departures of the incident drive bunch radial profile from its assumed Gaussian shape prior to plasma focusing, and depletion or reshaping of focused drive bunches for z > 30 cm, which are neglected in quasistatic LCODE simulations, are possible sources of such deviations. In addition, neglected drive bunch evolution within the ∼ 100 cm probed region imprints left-right asymmetry onto probe images beyond that currently simulated.
These residual discrepancies indicate that detailed quantitative comparison of simulated and measured nsscale plasma dynamics will require selected improvements to both experiment and simulation, as noted above. Nevertheless, the broad agreement obtained in the spatial and temporal scale of post-wakefield expansion validates the basic plasma/atomic physics on which simulation 3 is based. Its output can thus elucidate additional internal properties of the expanding plasma beyond those that were directly observed.
An example is the plasma's energy budget. According to simulation 3, the fully focused drive bunch [ Fig. 3 Fig. S2). The latter rate is expected to be lower than the former, since energy deposition grows as the bunch approaches full focus at z = 31.9 cm, and decreases thereafter as the drive bunch weakens, evolution that a quasistatic simulation neglects. Fig. 5(c) shows additionally how deposited energy partitions and evolves over 1.3 ns, based on simulation 3 results. Initially (∆t ≤ 1 ps), beyond the horizontal resolution of Fig. 5, most of the deposited energy is stored in electromagnetic fields of the electron wake, with a minor fraction stored in kinetic energy of coherently moving plasma electrons at the bubble boundary, as is typical for wide bubbles 35 . After the wave breaks [∆t ∼ 1 ps, Fig. 3(c)], about 10% of deposited energy transforms to kinetic energy of fast electrons [see Fig. 5(c), light blue] that diverge radially, some forming a tail wave 36 visible in Fig. 3(b). As the fastest electrons escape from the plasma, a radial charge separation field appears [see Fig. 3(c) and Fig. 5(c), green] that holds most remaining electrons within the plasma column 14 . In the first tens of ps (i.e. the range of OSIRIS simulations and the work of Ref. 5 ), electron kinetic energy (blue) and electric field energy (green) comprise most of the plasma columns energy. During the interval 30 < ∆t < 300 ps, the charge separation field accelerates ions, which acquire most of the energy by ∆t ∼ 300 ps [ Fig. 5(c), red]. The fastest ions accelerate beyond 400 keV (i.e. v = 3.6 × 10 6 m/s). Overall 90% of the initially deposited energy remains in or near the plasma column. This contrasts with radially unbounded plasmas, in which fast electrons escape the heated region freely, cooling it to sub-keV temperatures within a few hundred wake periods 15 . Here, plasma energy and its partitioning among plasma species reach steady state for ∆t > 300 ps, and change negligibly through the end of the simulation 3 run at ∆t = 1300 ps. This validates the experiment's premise that the Li blanket records energy transport via ionization without noticeably depleting the energy of the ionizing radiation. It also indicates that the ∼ 10 6 m/s expansion continues unabated well beyond ∆t ≈ 1.3 ns.
In summary, results of this study have identified the principal physical mechanisms, and quantified the dominant dynamical pathways, by which highly nonlinear e-beam-driven wakes in finite-radius n e ∼ 10 17 cm −3 plasmas release their stored electrostatic energy into the surrounding medium. Time-resolved optical diffractometry measurements of the expanding plasma column, in particular, prompted recognition of the critical, previously unrecognized role of ion-mediated impact ionization in driving the plasma radius outward during the first nanosecond. The results also make clear that the plasma columns internal electrostatic fields, even after the original wake breaks, remain not only the principal propellant of outward ion motion, but the principal force responsible for retaining 90% of the wake's initially deposited energy within the plasma column for over a nanosecond. The framework hereby established and validated provides a basis for modeling the global thermodynamics of multi-GeV plasma wakefield accelerators, and for evaluating limits on their repetition rate. Relevant extensions of the experiments and simulations include investigations of laser-driven wakefield accelerators, and the use of varied probe geometries (e.g. larger θ) that will enable space-and time-resolved observation of the plasma columns evolving internal structure.

Methods
Electron drive bunch characterization. The first 2 km of the SLAC linac delivered electron drive bunches to the FACET interaction region located in sector 20. A series of 8-turn toroidal current transformers along the FACET beamline measured the charge of each bunch with 2% accuracy. A synchrotron-x-ray-based spectrometer measured the energy spectrum of each incident bunch non-invasively with ∼ 0.1% resolution 37 . An integrated transition radiation monitor and a transverse deflecting cavity, both located just upstream of the FACET interaction region, measured transverse (σx,y) and longitudinal (σz) dimensions, respectively, of bunches entering the plasma with ∼ 10 µm resolution. We measured σx,y for every shot, σz for selected shots. By modeling the beam focusing optics, beam dimensions inside the interaction region were determined with similar accuracy. Beams incident on FACET typically had asymmetric emittance x ≈ 10 y , so in order to have round beams with σr ≈ σx ≈ σy ≈ 30 µm at the entrance of the plasma, we focused the beam into the plasma with beta-functions βy ≈ 10βx. See Ref. 38 for a detailed overview of FACET beam diagnostics.
Lithium source. Lithium was chosen for the accelerator medium because its low first ionization potential (5.4 eV) allows the drive bunch to singly field-ionize it easily over a 1.5 m path. A heat-pipe oven -consisting of a stainless steel cylindrical tube (length 2 m, inner radius 2.5 cm) heated along its center, and cooled at both ends -generated and confined the lithium gas. The vapor pressure of a melted lithium ingot loaded onto a stainless steel mesh ("wick") lining the inner wall of the hot center generated gas of temperature-controlled density na. Helium buffer gas concentrated at the cold ends confined it longitudinally. The longitudinal density profile na(z) was deduced from the temperature profile T (z) along the length of the oven, measured by inserting a thermocouple probe into the heat-pipe 6,20 . Thermocouple scans with our normal operating heating power of 1340 W yielded a 1.2 mlong central plateau of density na = 8.0 ± 0.2 × 10 16 cm −3 , with 0.15 m long density ramps at each end (see Supplementary Material Fig. S1). With the heat-pipe in steady-state operating mode, we controlled overall Li density primarily by adjusting buffer gas pressure, and the length of the central density plateau by adjusting heater power.
Probe laser and imaging system. Probe pulses (∼ 1 mJ energy, polarized in plane of incidence) were split from the 500 mJ, 50 fs output pulse train of a 10-TW Ti:S laser system 21 . Transverse probe intensity profiles contained hot spots, which were superposed on single-shot images (see Supplementary  Material Fig. S3). Since hot-spots varied from shot to shot, while e-beam and plasma structures were comparatively stable, 30-shot averaging removed probe artifacts from the data, while sharpening details of shadow/fringe patterns [ Fig. 1(c)] for quantitative analysis. Probe angle θ = 8 mrad was chosen to highlight the plasma column's expanding edge, while satisfying space constraints at the ends of the heat-pipe oven, the only optical access to the plasma column. Fig. S4 in Supplementary Material compares paths of 0.8 µm wavelength probe rays through an idealized plasma column, similar to the actual one, for θ = 0.008, 0.02 and π/2 rad. For θ = 0.008 rad, probe rays sample only the plasma columns outer edge, as discussed in connection with main text Fig. 2(a), (b). For θ = 0.02 rad, they penetrate to the plasma axis, enabling probing of interior structures. For θ = π/2 rad, probe rays are nearly un-deflected, rendering the plasma invisible. Transverse probe pulse width (w0 ≈ 0.4 cm), chosen to ensure nearly constant beam size through the heat-pipe, enforced lower limit θmin ≈ 0.007 rad to avoid perturbing the e-beam driver with a probe injection mirror of radius πw0/2, and upper limit θmax ≈ 0.012 rad to avoid clipping the probe beam on the heat pipe's end flange windows of radius 2 cm. A single-element lens (2-inch diameter, 1 m focal length) im-aged probe pulses to 12-bit CCD camera (Allied Vision Technologies Model Manta G-095B, 1292 × 734 pixels, pixel size 4.08 × 4.08 mm).
Simulations. OSIRIS simulations of ionization and wake formation [ Fig. 3(a), (b)] used a moving simulation box of dimensions Lr × Lz = 564 × 940 µm, divided into 0.94 × 0.94 µm cells with 12 × 12 particles per cell. To model ion motion [ Fig. 3(e)] driven by electrostatic energy stored in the nonlinear electron wake [ Fig. 3(c), (d)], we switched to a static Lr × Lz = 0.940 × 3.025 mm simulation box in the Li density plateau divided into 0.47 × 1.21 µm cells with 10 × 10 particles per cell. Impact ionization of neutrals is small on a tens-of-ps time scale, and was not included in OSIRIS simulations.
For LCODE simulations, our simulation window extended laterally out to r = 9.4 mm, with grid size 0.19 µm. The initial plasma consisted of 5 × 10 4 equal-charge macro-particles of each type, which increased as impact ionization proceeded. We implemented elastic collisions via the Takizuka-Abe model 28 modified to include relativistic particles 29 . Approximate cross-sections for electron-and ion-impact ionization came from Ref. 30 and Ref. 31 , respectively. Modeling of two-step ionization of lithium neutrals was based on excitation and ionization cross-sections from Refs. 32,33 .
We simulated probe images [e.g. Fig. 1(d)-(f)] by numerically propagating a probe pulse (λ = 800 nm) at angle θ from the axis through a cylindrically symmetric column consisting of a mixture of atomic Li in 2S ground, and 2P, 3S, 3P excited states, plus singly-ionized Li plasma. We calculated the index of refraction ηi = ηi(ω, r, ∆t) of atomic Li in each of the four states (i = 1, 2, 3, 4) at probe frequency ω = 2πc/λ, radius r, and time delay ∆t from the Lorentz-Lorenz relation where oscillator strengths f k , damping factors γ k and resonant frequencies ω 0k , taken from spectroscopic data in Ref. 39 , are known with ±0.3% uncertainty, yielding uncertainty of order ±2% in ηi − 1 for given density Ni. The refractive index of singly-ionized Li plasma was given a Drude form  (c) corresponding probe profile after non-ideal imaging to detector, showing re-shaping of vertex region due to partial blockage of imaging lens shown in Fig. 1(b).