The Guderley model of a self-similar imploding shock based on the group invariance of the flow equations is a powerful tool in understanding the behavior of converging shock waves. Two modifications described here improve the predictions of observable quantities in spherical-shock wave experiments. First, a noninfinite boundary condition is established by the isentropic release of the outer pressure. Second, a two-temperature system of ions and electrons allows description of higher temperatures while conserving energy and without perturbing the overall hydrodynamics of the solution. These modifications of the Guderley model improve the prediction of the observables in laser driven spherical shock experiments in reference to a one dimensional (1-D) hydrodynamics code.

## I. INTRODUCTION

Converging shock waves appear in diverse locations throughout nature such as in inertial confinement fusion experiments,^{1–3} supernovae explosions,^{4,5} and cavity collapses.^{6} Guderley^{7} found a self-similar solution that can be used to describe all these systems across the different size and temporal scales found in nature. Despite significant work^{8–10} exploring the numerical properties of the Guderley solution and using it to benchmark hydrodynamic codes,^{11–13} little work exists in using the solution to understand experimental results. Work has been done to treat the boundary of self-similar solutions, such as the Guderley, in order to better relate them to numerical hydrodynamic codes.^{12–14} This work describes a new boundary condition and scheme for the treatment of temperature, motivated by experimental comparison.

Shock wave experiments in spherical geometry are challenging to perform and diagnose, and thus direct isolated measurements of extreme physical states or processes accessed with such convergent experiments are usually beyond state of the art techniques, so design and analysis of such convergent experiments often employ radiation hydrodynamic codes.^{2} Semianalytic solutions to hydrodynamic systems offer a valuable baseline for comparison to integrated models.^{12} These solutions are used to inform experimental observations in diverse regimes of spherical shocks, for example, in rarefied gas.^{15} A direct analog to these solutions is challenging to create experimentally,^{16} so modifications of the solution are sometimes necessary^{16,17} to understand experiments outside the realm of the solution such as the ones containing magnetic fields.^{18,19} The scheme presented in this work when combined with previous work,^{15,20–22} makes for a powerful tool in quantitatively determining the states and processes implicit to convergent experiments and is an important step toward constraining physical models from a collection of experimental observables.

### A. Review of the Guderley problem

The Guderley problem of a converging shock wave has been described extensively elsewhere.^{1,8,23} The solution starts with the nonviscous ideal conservation equations of fluid mechanics as follows:

where $u\u2192$ is the velocity field, *P*, *ρ*, and *S* are the gas pressure, density, and entropy, respectively, and the three equations correspond to the conservation of momentum, mass, and energy (entropy), respectively. Additionally, for this solution an ideal adiabatic equation of state is assumed, taking the form

where *V* is the volume, *N* is the number of particles, *k*_{B} is the Boltzmann constant, and *γ* is the ratio of specific heats or the adiabatic index. Finally, invoking spherical or cylindrical symmetry, the problem is reduced to a 1-D radial set of equations of the form

where *n* is a parameter that sets the geometry of the problem taking on the values 1, 2, and 3 for planar, cylindrical, and spherical geometries, respectively. These equations evolve in a self-similar fashion with the similarity coordinate *ξ* defined as

The choice of this coordinate is derived from the scaling symmetries that exist in the equations.^{1}^{,}*ξ* introduces the parameters *ξ*_{0}, *r*_{0}, *t*_{0}, and *α*. Both *ξ*_{0} and *α* help define the shock trajectory and strength, and there is a unique *α* for the choice of *γ*, the ratio of specific heats or the adiabatic index, and the geometric factor, n. It is important to note that *α* is independent of boundary conditions, but *ξ*_{0} depends on the specifics of the boundary, in particular the strength of the pressure wave at infinity. Scaling parameters *r*_{0} and *t*_{0} set the spatial and temporal scales of the system being modeled. This allows the fluid properties to be given in terms of the similarity coordinate. In order to solve the system, only three of the fluid properties (particle velocity, sound speed, density, temperature, pressure, entropy, etc.) need to be defined and the rest can be derived. Different references choose to define different properties. As an example,^{10,24} the flow velocity, density, and sound speed are given as

with two additional parameters: *ρ*_{0} being the reference density and *κ* being the density exponent defining the initial density profile. Inserting the scaled quantities into the flow equations allows one to derive a set of autonomous ordinary differential equations (ODE) of the form

These equations can be solved numerically and then transformed back to the radius and time space.^{1,13,24} The parameters introduced in this solution are given in Table I.

r_{0}
. | Radial scaling factor . |
---|---|

t_{0} | Temporal scaling factor |

ξ_{0} | Shock strength parameter |

γ | Adiabatic index |

n | Geometric factor |

κ | Density profile parameter |

ρ_{0} | Density scaling factor |

r_{0}
. | Radial scaling factor . |
---|---|

t_{0} | Temporal scaling factor |

ξ_{0} | Shock strength parameter |

γ | Adiabatic index |

n | Geometric factor |

κ | Density profile parameter |

ρ_{0} | Density scaling factor |

### B. The Guderley solution

The particle and shock trajectories that result from solving the Guderley problem provide insight into the dynamics of a converging shock wave. These trajectories are shown in Fig. 1 where the temporal axis is shifted by *t*_{c}, the time that the shock reaches the center, *r* = 0. The slope of the profiles is determined by the equation of state, namely, the *γ* used, but the general features of the solution are independent of *γ*. The particle trajectory originating from *r*/*r*_{0} = 1 is the dashed-dotted line. This trajectory has no unique significance in the Guderley solution because the flow extends out to infinity and is self-similar, but in experimental configurations there is usually a reference radius that carries significance such as an outer boundary of the material. This boundary will be discussed in Sec. II.

The density and temperature profiles are critical when describing shock waves in nature. The right-side top and bottom of Fig. 1 show temporal and spatial lineouts of the temperature and density, respectively, denoted by arrows in the particle trajectory plot. The temperature at a fixed radius is relatively steady in time between the two shock events. Conversely, the temperature is extremely steep in its radial profile with the temperature peak at *r* = 0 for all times after the shock collapse. The density is shocked as expected according to the strong shock condition *ρ*/*ρ*_{0} = (*γ* + 1)/(*γ* − 1). After that point, the density continues to rise, a consequence of the material converging into a smaller volume. Finally, after reshock the material relaxes and the density falls off. A single fluid temperature is shown but a scheme for separate treatment of electron and ion temperatures will be addressed in Sec. III.

## II. FINITE BOUNDARY CONDITIONS

The self-similarity of the Guderley solution means the flow extends out to infinity. This implies that the shock wave is supported for all time and that the problem describes a system of infinite mass. In real physical systems, the finite boundary conditions often play a key role in the dynamics, necessitating a treatment of the boundary condition. The simplest condition that can be used is to track the flow from a reference point, for example, *r*_{0}, which may correspond to the outer boundary of an experimental target (dashed-dotted line in Fig. 1) and truncate the solution outside of that boundary. This scheme is valid while the pressure at the boundary is held constant, supporting the wave dynamics for the duration of the shock. The support of the pressure wave in the Guderley solution is a problem when the pressure is not held at the outer boundary indefinitely, such as in laser-driven shock experiments.

When the pressure at the boundary is released, a decompression wave is launched, changing the fluid property profiles. To model this, a released density profile of the form

is assumed, based on previous solutions to isentropic release waves in imploding media using the conservation of mass.^{25} In this work, *ε* will be left as a free parameter that will be system specific. Here *R*_{in} and *R*_{out} are the inner edge and outer edge of the release wave, respectively, and *ρ*_{R} is the value of the density at the inner edge of the release wave, given by the density of the original Guderley solution. *R*_{out} is the boundary of the system beyond which there is no material. *R*_{in} is calculated using the sound speed in the compressed and moving material. Since the sound speed and density are given by the Guderley solution, the only unknown in the density profile is *R*_{out}. To find *R*_{out}, conservation of mass is used. The mass contained in the shell of the material bounded by *R*_{in} and *R*_{c}, the location of the outermost fluid element, must equal the mass within the newly released profile:

where *ρ*_{G}(*r*, *t*_{i}) is the undisturbed radial Guderley density profile at time *t*_{i}. Solving the integral on the right hand side results in a third order polynomial for *R*_{out} given as

Solving this polynomial at each time step gives the trajectory of the outer edge of the release wave. This process is truncated when the inner edge of the release begins to interact with the rebounding shock wave. When this occurs, the trajectory of the outer edge is then extrapolated and the density at the inner edge of the release is solved by rearranging the polynomial to get

The density profiles for the standard Guderley solution and the release modification are compared in Fig. 2.

As the shock moves out reshocking material, it will eventually interact with the released material, modifying the density, pressure, and temperature. This is treated by scaling the density behind the second shock by the ratio of the released density value and the Guderley density value at the location of the reshock:

where +*R*_{in} is referencing the density value on the singly shocked side of the discontinuity (the white region in Figs. 1 and 3). This is done in order to preserve the shock condition which is determined by the equation of state, through *γ*. The different trajectories and density regions are shown in Fig. 3.

Finally, following from the adiabatic relationship *PV*^{γ} = constant and the ideal equation of state, the temperature and pressure profiles in the released material are given by

The modified temperature profile can now be used as the single fluid temperature inherent in the Guderley solution. A two-fluid temperature calculation is discussed in Sec. III.

## III. TWO-FLUID TEMPERATURE TREATMENT

The Guderley problem is usually solved in terms of entropy (as developed in Sec. I A) or sound speed [as given in Eq. (5)] both of which can be related to the temperature of the fluid. Working in the framework already established for the problem, the entropy or sound speed can easily be converted into the temperature through the equation of state, for example, $c=\gamma kBT/m$, where *m* is the mass of the particles being considered.

One primary result from the Guderley solution is that the temperature, while being a single-fluid temperature, reaches extremely high values at the collapse point because of the strong shocks in the system. The combination of high temperatures and short time scales that are associated with shock waves necessitates the treatment of ions and electrons as two fluids with separate temperatures. The presence of these large temperature gradients can be addressed through other methodologies, for example, by including heat conduction.^{26} This work includes a two fluid treatment, specifically with the goal of reproducing experimental observables that have separate dependencies on the electron and ion temperatures as shown in Sec. IV.

The separation of the temperature into two different populations begins with how the energy is partitioned into those populations. Throughout this work, charge neutrality is assumed such that $ne=Z\xafni$; the number density of electrons in a fluid element is equal to the average ionization state times the number density of ions. In shocked material the particles will move with the same velocity so the energy imparted to the particle is proportional to the mass of the particle. In other words, the ratio of the energy goes as the ratio of the masses

where *m*_{e} is the electron mass and $A\xaf$*m*_{u} is the ion mass, with $A\xaf$, the average atomic mass number, and *m*_{u}, the atomic mass unit. The temperature is proportional to the energy for an ideal gas population of particles given by *E* = *Nk*_{B}*T*/(*γ* − 1) (assuming that *γ* is the same for electrons and ions), so the temperature of the population is also proportional to the mass of the particles comprising the population. The shock leaves the electrons and ions with a temperature ratio of

following from the energy of each population, determined by their masses. Here $Z\xaf$ is the average ionization state of the ions, which determines how many electrons there are per ion.

The strong shock wave in the problem leaves the electrons and ions far out of equilibrium. Considering *τ*_{ei} as the characteristic time scale of equilibration between electrons and ions, the electron-electron and ion-ion time scales, *τ*_{ee} and *τ*_{ii}, are much smaller; i.e., they equilibrate faster as a result of the greater efficiency of the energy transfer between collisions of particles with the same mass.^{27} Since this is the case, the electron and ion distributions are assumed to be Maxwellian at all times but relax to different temperatures, *T*_{e} and *T*_{i}.

The equilibration is treated in a way that leaves the single-fluid temperature of the Guderley solution, *T*_{G}, unperturbed, in order not to change the hydrodynamic evolution of the system. The electron and ion temperatures are constrained to obey the following:

i.e., the energy in any given fluid element remains the same and is just partitioned between the electrons and ions. Assuming that *γ*_{i} = *γ*_{e}, $ne=Z\xafni$, and *n* = *n*_{i} + *n*_{e}, the relationship between the temperatures is given by

Using this constraint, the electrons and ions are equilibrated according to^{27}

Here *τ*_{ei} is usually given as a function of electron temperature as well. The Guderley solution gives *T*_{G}(*r*, *t*), so using this and the constraint on the energy partitioning gives an ODE for the electron temperature as a function of time for a given radius *r*_{i},

This can then be numerically solved, and the ion temperature as a function of time is given as $Ti=1+Z\xafTG(ri,t)\u2212Z\xafTe(t)$. The temperature dependence of $Z\xaf$ can also be entered into this equation, but for this work the calculation of $Z\xaf$ was made using the initial ion temperature postshock before equilibration. Figure 4 shows a comparison of the temporal profiles of the Guderley, electron, and ion temperatures.

## IV. EXPERIMENTAL USAGE

The temperature equilibration and release modification of the Guderley solution are used to predict the observables of a high-energy-density experiment. The predictions are compared to the 1-D hydrodynamics code *Lilac*.^{28} The experimental configuration includes a strong laser-driven spherical shock in a solid deuterated polystyrene ball that is 0.890 mm in outer diameter. This is a type of experiment conducted on the *Omega* Laser system^{29} at the Laboratory for Laser Energetics. The key observables of interest are x-rays and deuterium-deuterium fusion neutrons produced around the time of shock collapse. To properly model these in the modified Guderley solution, the addition of an ionization model and x-ray opacity model is necessary; neutron scattering is neglected here due to its negligible effect on calculating the total neutron yield. Neutron scattering is important for higher order moments of the neutron spectrum which are not calculated in the current framework but can be added in future studies. The ionization model from Hu *et al.*^{30} is used, providing a Saha-like framework for calculating the ionization levels in CH (the isotopic difference between CH and CD is neglected). Tabulated opacity values^{31} are used and the emission is assumed to be emitted radially. The choice of ionization model and opacity values are used in order to make comparison to *Lilac*, with the same models being used in both. A benefit of this methodology is that different models can be easily implemented and the sensitivity to these models can be explored. The x-ray emission is calculated assuming only Bremsstrahlung emission,^{1,23,32} due to the lack of spectral lines in the region of experimental observation in the few keV x-ray energy range, according to

and the neutron emission is calculated according to

where ⟨*σv*⟩_{DDn} is the deuterium-deuterium thermal reactivity as given by, for example, Bosch and Hale.^{33} All quantities are in SI units and the resulting units are given in brackets.

### A. Prediction of observables

The parameters of the Guderley solution including the modifications are chosen based on the experimental setup or comparison to the hydrodynamics code. The parameters are given in Table II. Many of the values used are given by nominal target specifications (*n*, *ρ*_{0}, *r*_{0}, *κ*) or laser configurations (*t*_{r}, *t*_{0}); *ε* is chosen based on previous literature on shock releases^{25} and *τ*_{ei} is given by Spitzer^{27} as

where all quantities are in SI units, the resulting units are given in brackets, and ln(Λ) is the Coloumb Logarithm given^{27} as

r_{0} | Radial scaling factor | 443 μm |

ξ_{0} | Shock-strength parameter | 151.6 μm/ns^{α} |

t_{0} | Temporal scaling factor | $(r0/\xi 0)1/\alpha =4.58$ ns |

γ | Adiabatic index | 1.5 |

α | Derived from γ | 0.7044 |

n | Geometric factor | 3 |

κ | Density profile parameter | 0 |

ρ_{0} | Density scaling factor | 1.04 g/cm^{3} |

ε | Release density exponent | 7/2 |

t_{r} | Time release is launched | 0.45t_{0} |

τ_{ei} | Electron-ion equipartition time | Equation (21) |

r_{0} | Radial scaling factor | 443 μm |

ξ_{0} | Shock-strength parameter | 151.6 μm/ns^{α} |

t_{0} | Temporal scaling factor | $(r0/\xi 0)1/\alpha =4.58$ ns |

γ | Adiabatic index | 1.5 |

α | Derived from γ | 0.7044 |

n | Geometric factor | 3 |

κ | Density profile parameter | 0 |

ρ_{0} | Density scaling factor | 1.04 g/cm^{3} |

ε | Release density exponent | 7/2 |

t_{r} | Time release is launched | 0.45t_{0} |

τ_{ei} | Electron-ion equipartition time | Equation (21) |

The choice of *γ* and *ξ*_{0} follows from matching the initial strength (pressure) of the shock and the time of shock collapse between the *Lilac* and the Guderley solution. Once the full solution is established, a set of common experimental observables is calculated in order to compare. These quantities include the total neutron yield, given as the integral of Eq. (20) over volume and time; the total x-ray yield, given as the integral of Eq. (19) over volume, time, and spectral frequency; and the x-ray energy where the peak of the spectrum is located. The x-ray quantities are dependent on the absorption, which is a strong function of the density, while the neutron emission is primarily a function of temperature. The release modification, which primarily changes the density profile, has the most leverage on the x-ray emission quantities, while the temperature equilibration scheme determines the temperature affecting both the x-ray and neutron emission quantities.

### B. Comparison to *Lilac*

The comparison simulation was run in *Lilac* using nonlocal heat transport,^{34} the same opacity^{31} and ionization^{30} values used in the Guderley calculations, and the *sesame* equation-of-state tables.^{35} As was previously stated, the only tuning required for the modified Guderley model to match the *Lilac* results was based on the initial (ablation) pressure and the time of shock collapse. The *Lilac* simulation has an ablation plasma resulting from the direct laser illumination on the outside of the target, which can contribute to self-emission quantities; *Lilac* quantities were calculated using only the region inside of the ablation plasma, i.e., only contributions from the shock-generated states. The resulting trajectories are compared in Fig. 5. The trajectories show good agreement, deviating less than 5% on the ingoing trajectories, with the largest deviations occurring with the release around the time of shock collapse and with the shock postcollapse. The inner release trajectory is calculated according to the sound speed so the difference in trajectory is due to a difference in local sound-speed between *Lilac* and the Guderley model. Since the sound speed is dictated by the equation of state, deviations of the used equation of state from an ideal gas account for the differences. The postcollapse shock trajectory is unmodified by the release in this model, which accounts for the Guderley shock diverging more slowly than the *Lilac* shock postcollapse. This is not, however, a critical period of these experiments because by this time most observable quantities have already been produced. The neutron and x-ray emission quantities are compared in Table III. The quantities are calculated for the *Lilac* simulation and the Guderley solution with and without the release model both for the fully equilibrated case and the case with dynamic equilibration. The neutron yield is primarily a function of the ion temperature, which is notably under-predicted in both *T*_{e} = *T*_{i} cases. The x-ray quantities are functions of electron temperature, but they primarily depend on the mass density due to the attenuation of the x-rays. The models without the boundary treatment underpredict the x-ray yield and overpredict the peak spectral energy because only the high-energy x-rays are able to escape the system. The model with the equilibration and boundary treatment addresses both of these issues and is within 3%, 30%, and 13% for the neutron yield, x-ray yield, and spectral peak, respectively.

Simulation . | Y_{N}
. | Y_{X} (mJ)
. | Spectral peak (keV) . |
---|---|---|---|

Lilac | 7.320 × 10^{6} | 1.151 | 7.1 |

Gud (T_{e} = T_{i}) | 6.625 × 10^{6} | 0.078 | 10.0 |

Gud (τ_{ei}) | 7.742 × 10^{6} | 0.016 | 15.0 |

Gud (T_{e} = T_{i} + Rel) | 5.874 × 10^{6} | 1.410 | 7.0 |

Gud (τ_{ei} + Rel) | 7.116 × 10^{6} | 0.895 | 8.0 |

Simulation . | Y_{N}
. | Y_{X} (mJ)
. | Spectral peak (keV) . |
---|---|---|---|

Lilac | 7.320 × 10^{6} | 1.151 | 7.1 |

Gud (T_{e} = T_{i}) | 6.625 × 10^{6} | 0.078 | 10.0 |

Gud (τ_{ei}) | 7.742 × 10^{6} | 0.016 | 15.0 |

Gud (T_{e} = T_{i} + Rel) | 5.874 × 10^{6} | 1.410 | 7.0 |

Gud (τ_{ei} + Rel) | 7.116 × 10^{6} | 0.895 | 8.0 |

The x-ray emission is dominated by the absorption that is determined by the opacity and areal density (*ρR*) of the material. Figure 6 compares the zeroth and first radial moments of the density distribution as a function of time. The zeroth moment gives the average areal density, ⟨*ρR*⟩, the ratio of the first and zeroth moments gives an average radius, and the second moment is directly proportional to the mass. The first three moments are given by

Since mass is conserved throughout all of the models, the second moments are trivially identical throughout, so the comparison is limited to the zeroth and first moments.

Figure 6 shows that the standard Guderley model greatly overpredicts both of the moments. The zeroth moment of the modified Guderley solution is larger than *Lilac* by about 20% around the time of shock collapse. This is caused by the different timing of the release wave late in time, which is seen in Fig. 5. The difference in the zeroth moment around the time of peak x-ray emission (≈5 ns) corresponds to a greater areal density and is the reason for the lower x-ray yield (Table III) and harder emission predicted by the Guderley solution with release resulting from greater absorption. The first moment is also larger in the Guderley solution with release by about 20% before the time of shock collapse, but after collapse, the value is much closer to *Lilac* below about 5%.

Figure 7 shows how the x-ray yield scales as a function of the initial density for the different model shown in Table III. At lower densities, the attenuation is minimized and the electrons and ions are far out of equilibrium, so the solutions with and without the dynamic equilibration converge to the same values, respectively. At higher densities, the electrons and ions fully equilibrate very quickly but the attenuation is magnified, so the solutions with and without the release model converge to the same values. Only the fully modified Guderley model is able to reproduce the scaling of the x-ray yield over the entire density space considered.

## V. CONCLUSIONS

The Guderley solution of an imploding shock wave with the addition of a boundary condition and electron-ion energy partitioning is effective at reproducing the simulated behavior of a high-energy-density experiment of interest. The particle trajectories of the Guderley solution with the release model boundary condition and *Lilac* simulations are consistent. Self-emission quantities compared between the Guderley solution and the *Lilac* simulation are within about 10% or less for most quantities when the complete Guderley system with boundary condition is considered. The favorable comparison between the hydrodynamics code and the Guderley solution with temperature equilibration and the release model demonstrates that the physics contained in the Guderley model is sufficient to explain most of the system. Additionally, similarity in the experimental quantities predicted by the two models demonstrates that these measurements are not sensitive to the additional physics contained in the hydrodynamics code, such as thermal conduction and radiation transport.

This model and similar models of this type are a critical step toward understanding which physics dominates a particular experimental setup and the physical mechanisms that can be probed using current observations. Additionally, this model has fewer parameters than a typical hydrodynamics code and can be better constrained by a limited set of experimental data when doing a model-fitting through, for example, a Bayesian framework, presenting a path forward in understanding which aspects of fundamental physics are dominant in integrated high-energy-density experiments.

## ACKNOWLEDGMENTS

This material is based upon work supported by the Department of Energy, National Nuclear Security Administration, under Award No. DE-NA0003856, the U.S. Department of Energy, Office of Science, Office of Acquisition and Assistance, under Award No. DE-SC0019269, the University of Rochester, and the New York State Energy Research and Development Authority.

This report was prepared as an account of work sponsored by an agency of the U.S. Government. Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.