A high temporal resolution numerical algorithm for shock wave velocity diagnosis

We propose a high temporal resolution numerical algorithm for shock wave velocity diagnosis. By analysing variations in the optical path and phasor of a light field, we determine a high temporal resolution shock wave velocity equation for a velocity interferometer system for any reflector (VISAR). The equation can be transformed into matrix form for numerical solution. To solve noise problems, a ‘filtering velocity spectrum’ method is proposed. Analysis of a VISAR data example shows that the resolution precision of shock wave velocity obtained from the numerical algorithm is the same as the temporal resolution of a streak camera. Moreover, it can observe the shock wave in greater detail. This algorithm can be used to observe detailed images and determine the mechanism and evolution of extreme shock waves, as well as provide data for research into hydrodynamic behaviour in inertial confinement fusion.


Derivation of the phase equation
The principle of VISAR can be described by the simple schematic diagram shown in Fig. 1. The x-axis coincides with the slit of the streak camera and the y-axis coincides with the optical axis of VISAR. The etalon is set in branch 1. The solid angle for probe beam collection in VISAR discussed here is less than 7°; thus, the influence of frequency correction can be ignored 18 . First, the probe beam is incident on the shock wave interface through a beam splitter. Because of the optical Doppler effect, the reflected beam carries shock wave velocity information. The difference between the wavelength of the information beam and the probe beam is λ λ Δ = − u c 2 0 , where u is the shock wave velocity, c is the speed of light, and λ 0 is the wavelength of the probe beam. Then, the information beam enters the interferometer through the imaging system and is split into two beams by another splitter. One , where h is the thickness of the etalon, n is the dynamic refractive index of the etalon, and n 0 is the static refractive index of the etalon. Finally, the two beams interfere near the slit of the streak camera. By analysing the recorded phase shift of the interference fringe, the shock wave velocity can be obtained.
In the conventional algorithm, the number of dynamic waves in the etalon is considered to be = τ λ N c , where λ is the dynamic wavelength 22,30 . Thus, the phase shift of the interference fringe is is the static delay time. The wavelength of the VISAR probe beam source at the SG-III prototype is 532 nm, and the etalon is fused silica. Therefore, n 0 = 1.45847 and the correction factor is δ = − = .  (1) can be rewritten as follows: 0 0 0 Equation (2) expresses the relationship between the shock wave velocity and the fringe phase calculated by the conventional algorithm. Although τ and λ are considered variables, the conventional algorithm is established under the condition that the wavelengths of the probe beam in the etalon are equal everywhere. In addition, the optical path change between the target and the camera is not considered, and the optical path change of the target surface to the probe beam source is also neglected. Thus, the influence of the initial phase value of the two branches of the interference is not evaluated. These approximations reduce the temporal resolution of VISAR but can be considered by solving the phase vector equation. The following is a derivation of the relationship between the shock wave velocity and the fringe phase obtained using the phase vector equation.
Assuming the light field direction is perpendicular to the image plane, the phasor of the light field for branch j around the streak camera slit at T can be determined by the optical path: , ϕ x y t ( , , ) 0 0 0 is the initial phase of the light source, Δφ j is the phase shift of the information beam from the detection surface to the streak camera slit of branch j, T j is the generation time for the information beam of branch j, and S j is the optical path between the target and the streak camera at T j .
The following equations are the boundary conditions of equation (3): where T 0 is the static moment before compression, S 0 is the optical path between the target and the streak camera at T 0 , L j is the optical path between the light source and the target at T j , L 0 is the optical path between the light source and the target at T 0 , k is the longitudinal magnification of the VISAR imaging system, θ j is the angle between the light field propagation direction and the y-axis at branch j, (x j , y j ) is the coordinate of the dynamic image, and (x 0j , y 0j ) is the coordinate of the static image. In equation (4), the first term represents the optical path between the shock wave surface and the probe beam source; the second and third terms indicate that the probe beams reflected at different times are simultaneously recorded by the camera due to their passing through different branches; the fourth term represents the optical path between the streak camera slit and the shock wave; and the fifth and sixth terms represent the relationship between the coordinates of the static and dynamic images.
The light field at the slit of the streak camera is superimposed by two branches and the total light field intensity recorded by the camera is The experimental fringe pattern is generally filtered by the Fourier transform method; thus, the processed fringe pattern distribution is mainly determined by the phase change 33 . From equations (4) and (5), and the rotation transformation formulas The purpose of the rotational transformation is to establish the relationship between the image coordinates in branches 1 and 2. Here, (x 0 , y 0 ) is the static coordinate of the image plane when the light field direction coincides with the y-axis, and (0, m) is the intersection of the static images in branches 1 and 2 on the y-axis. Equation (7) has three time-related moments: T, T 1 , and T 2 , whose relationship is described in equation (4). Thus, The change speed of the three time-related moments can then be compared, Here, the change speed of T, T 1 , and T 2 can be regarded as synchronous (non-uniformity is less than one thousandth). Then, considering the above conclusions and differentiating equation (7), we can obtain Integrating equation (10) and transforming it by Taylor expansion, we obtain where Θ ∈ (0, 1). Equation (11) reveals that the conventional equation between shock wave velocity and interfer- is approximately a description of the average velocity in interval τ + t t ( , ) 0 (reflected in the first item of Taylor's expansion), which explains why the temporal resolution of the conventional calculation cannot be less than τ 0 24 . Moreover, the condition for establishing the conventional equation is that the change of u in interval τ + t t ( , ) 0 must not be overly large (reflected in the second item of Taylor's expansion) or overly fast (reflected in the third item of Taylor's expansion). The above conclusion is consistent with previous discussions of the application of the conventional equation to the numerical calculation of shock wave velocity 22,23,31,34 . When the shock wave velocity does not change rapidly, the conventional equation can easily calculate the velocity change caused by the phase change. However, when the shock wave velocity is generated, chased, or even oscillated, the instantaneous solution shown in equation (10) with a high temporal resolution provides a better description. Therefore, accurately solving equation (10) is crucial for shock wave velocity high resolution diagnosis.

Numerical Calculation
Equation (10) describes the relationship between shock wave velocity and the VISAR fringe phase in the form of a differential equation. The VISAR recording based on the streak camera is an averaging process. Here, it is considered that each pixel records the average phase in the camera temporal resolution interval. Then, the solution of equation (10) cannot be less than the temporal resolution of the streak camera; thus, it is necessary to convert equation (10) into a form that can be used for numerical calculations. When ϕ and u respectively represent the arithmetic average of the phase and velocity in the temporal resolution interval, equation (10) can be further changed to where ι is the temporal resolution of the streak camera. Equation (12) shows the average value of the shock wave velocity at each pixel of the camera and that the relationship between the fringe phase and velocity is nonlinear. There are always slight data fluctuations superimposed with the true phase in the phase recording pattern that arise from the recording response, noise, phase extraction, etc. 31 . These fluctuations are random and difficult to completely eliminate. In particular, specific frequencies in the fluctuations may be amplified during the calculation. When the fluctuation is considered, equation (12) can be written as: www.nature.com/scientificreports www.nature.com/scientificreports/ where ϕ x t ( , ) exp and ω x t ( , ) are the recorded phase and the arithmetic average of the phase random fluctuation in interval ι at t respectively. If the velocity and phase values of each time interval are expressed in the above formula, it can be written in a matrix form, Moreover, the shock wave velocity can be obtained from Here, u can be regarded as a one-dimensional matrix or as a velocity curve with respect to time, and the interval time of each element is the temporal resolution. Then, the shock wave velocity equation is transformed into a matrix describing the pixel points. It is then necessary to solve the shock wave velocity with a high temporal resolution (equal to the temporal resolution of the streak camera) through the matrix operation. Moreover, the key to solving the velocity is to remove the influence of noise from the operation process.

Noise Filtering
We determine that A is similar to a Toeplitz matrix 35 ; thus, A −1 is also similar to a Toeplitz matrix and the spectrum of discrete Fourier transforms for most rows of A −1 is approximately equal. Here, we term the discrete Fourier transform result of the velocity curve the 'velocity spectrum' , thus each element of the velocity spectrum of u satisfies the following: . The second term on the right side of equation (16) is the source of amplified noise. It is difficult to calibrate noise in the experiment, but we can design a window such as ∑ ( 1)( 1) to filter it. Then, by performing inverse Fourier transform on the result, the real velocity curve can be obtained. Here, we present an example of the noise filtering.
As shown in Fig. 2(a), we employ an experimental VISAR data set of 20150522037 shooting on the SG-III prototype facility. The target consists of three aluminium steps of different thicknesses. The probe beam reflected from the free surface is recorded by the VISAR. Because the driven laser pulse is modulated, the velocity rises with a certain gradient and gradually reduces owing to impact unloading. A part of the data set is selected for numerical calculation. It should be noted that the detection object of the VISAR can be a free surface or a shock wave surface, which is related to where the probe beam is reflected. The fringe pattern associated with the free surface velocity is chosen here because the data are representative and can be analogized to the shock wave velocity. An image of the phase function can be obtained through image processing (Fig. 2(b)). The temporal resolution of the streak camera is 20 ps. As the thickness of etalon is 17 mm, so the static delay time is τ ι = .
4 378 0 . The discrete Fourier transform spectra of most A −1 rows are equal. Figure 3(a) shows the value of the 200 th , 400 th , and 600 th rows of A −1 and their discrete Fourier transform spectra. The value of Fig. 3(b).
The coordinates of the peaks of the two curves coincide, which can be explained by equation (16). The second term on the right side of equation (16)  There are two aspects of the filtering process. One is that the Fourier transform spectra for the first few rows of A −1 are quite different from ∑ is taken as 1, a phase function is used as the input for both the conventional and proposed algorithms, and the velocity patterns of the entire space can be obtained as shown in Fig. 5.
A comparison of Fig. 5(a,b) shows that the velocity pattern obtained by the two algorithms are almost the same, except for the values in the rapid acceleration process. The proposed algorithm determines the uneven velocity distribution in impact loading, which indicates the possibility of uneven stress on the shock wave interface during acceleration. This uneven area only appears at approximately 400 ps then disappears. The temporal resolution of the conventional algorithm is approximately 100 ps; hence, this phenomenon is not clear. However, as the temporal resolution of the proposed algorithm is 20 ps, it can clearly identify this phenomenon.
It should be noted that the temporal resolution of the shock wave velocity obtained from proposed algorithm is the same as the temporal resolution of the recording device. In practice, the resolution of the streak camera is alterable owing to the parameter adjustment, and the etalon thickness is also variable because of experimental needs. When the etalon delay time is less or not significantly different from the temporal resolution of the streak camera, the difference between the two algorithm calculation results may be difficult to distinguish. However, because the spatial resolution of the camera is fixed, the speed resolution of the system decreases with the thickness of the etalon. Thus, the proposed algorithm has an advantage in terms of temporal resolution when a longer etalon is used and can provide a way to solve the shock wave velocity distribution when shock waves are generating, catching, or loading. On the other hand, the proposed algorithm is difficult to apply to the case where the fringe pattern has a jump. A few large non-random phase errors are amplified in the calculation and cannot be eliminated by existing filtering windows.  www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusion
To improve the temporal resolution of shock wave velocity in inertial confinement fusion, we proposed a numerical algorithm. Firstly, we reviewed the principle of VISAR as a standard piece of equipment for shock wave detection and noted that several conditions of the conventional algorithm reduce the temporal resolution. Then, by solving the phase vector equation of the light field near the detection surface, we obtained a precise differential equation for the shock wave velocity. To solve shock wave velocity, the equation was written in matrix form for numerical operations. We also determined that the noise fluctuation in the phase record is amplified during the calculation, so a filter window was designed for filtering the velocity spectrum. The shock wave velocity calculated by the algorithm has the same temporal resolution as the streak camera and can reflect velocity distribution information that is neglected by the conventional algorithm.
The proposed numerical algorithm provides a way to calculate shock wave generation, catching, and evolution at high temporal resolution, enabling precise observations of extreme fluid dynamic behaviour. In inertial confinement fusion, this algorithm could provide abundant data for researching laser-plasma interactions, hydromechanics instability, and driven laser tuning, as well as contribute to the creation of preparatory conditions for final ignition. Our numerical algorithm can even be extended as a calculation reference for other optical instruments containing an etalon, especially if there is a need for high temporal resolution.