2. QUANTUM HYDRODYNAMIC MODEL AND SIMULATION METHOD

A 2D hydrodynamic model is performed to investigate the time-dependent dynamic polarization and heating of two-component (ion and electron) aluminum solid target when irradiated by a proton beam pulse. In a cylindrical coordinate system with {r,z,φ}, φ is neglected in this work, since the system is symmetric with respect to z axis (r = 0). The 2D QHD model is employed to describe the target electron dynamics, with electron velocity field u_e (r,z,t), electron number density n _e (r,z,t), and electron energy density W _e (r,z,t). The dynamics of target ions and beam ions are described by a 2D classical hydrodynamic model, with the target ion velocity field u_i (r,z,t) and number density n _i (r,z,t), the beam ion number density n _b (r,z,t), and velocity field u_b (r,z,t). The hydrodynamic methods involve calculation of the conservation equations of particle number, momentum and energy, where a self-consistent determination of these conservative quantities is required with the help of Poisson's equation. A proton beam pulse with a particle energy of E _b (1.5–10 MeV) and an intensity of N (7.9 × 10⁸) ions that is the total particle number delivered in a single pulse time 10 ps, propagates along the z direction and irradiates perpendicularly on the target. The transverse Gaussian distribution is considered with a focal radius r _b (23 μm), which is defined as a multiple of the root mean square radius. The beam parameters here are selected to be relevant with the experiment carried by Patel et al. (Reference Patel, Mackinnon, Key, Cowan, Foord, Allen, Price and Ruhl2003).

In what follows, the simulation starts at t = 0 at a room temperature. At t = 0, the target ions and electrons are at equilibrium states and the velocities of them are all zero (u_e (r,z,t = 0) = 0,u_i (r,z,t = 0) = 0). The beam is projected from the axis z = 0 at t = 0 with a constant velocity and a initial density n _b0. Our simulation results are obtained for a semi-infinite long cylindrical solid aluminum target with initially equilibrium electron density n _e (r,z,t = 0) = n₀ = 1.18 × 10²³ and ion density n _i satisfying Z _in _i (r,z,t = 0) = n₀. Z _i = 3 is the effective charge state of aluminum target ion, which is equal to the number of valence electrons provided that the bound electrons are not ionized. At the axis z = 0 fixed boundary condition is applied for protons and open boundary condition is used for target ions and electrons. Equilibrium boundary conditions are implemented for z → ∞ or r → ∞ where there is no beam propagation and the target ions and electrons are not disturbed. At the axis r = 0 symmetric boundary condition is presented.

The target electron dynamics are described as:

(1)$$\displaystyle{{\partial n_e } \over {\partial t}}+\nabla \cdot \lpar {n_e {\bf u}_e }\rpar =0\comma \; \eqno \lpar 1\rpar$$

(2)$$\eqalign{ \displaystyle{{\partial {\bf u}_e } \over {\partial t}}+\lpar {{\bf u}_e \cdot \nabla }\rpar {\bf u}_e&=\displaystyle{e \over {m_e }}\nabla \Phi - \displaystyle{{\nabla P_e } \over {m_e n_e }} \cr & \quad+\displaystyle{{\hbar ^2 } \over {2m_e^2 }}\nabla \lpar \displaystyle{1 \over {\sqrt {n_e } }}\nabla ^2 \sqrt {n_e } \rpar - {\rm \gamma} _e {\bf u}_e\comma \; } \eqno \lpar 2\rpar$$

and the energy-balance equation (Chen et al., Reference Chen, Cockburn, Gardner and Jerome1995)

(3)$$\eqalign{ \displaystyle{{\partial W_e } \over {\partial t}}+\nabla \cdot \lpar W_e {\bf u}_e \rpar &=en_e {\bf u}_e \cdot \nabla \Phi - \nabla \cdot \lpar {\bf u}_e P_e \rpar \cr & \quad- \nabla \cdot {\bf q}_e - \displaystyle{{W_e - {\textstyle{3 \over 2}}n_e T_e } \over {{\rm \tau} _w }}\comma \; }\eqno \lpar 3\rpar$$

with $W_e={\textstyle{3 \over 2}}n_e T_e+{\textstyle{1 \over 2}}m_e n_e u_e^2 - {\textstyle{{\hbar ^2 n_e } \over {24m_e }}}\nabla ^2\, log\lpar n_e \rpar +O\lpar \hbar ^4 \rpar $. Thus, the electron temperature is $T_e={\textstyle{2 \over {3n_e }}}\left[{W_e - {\textstyle{1 \over 2}}m_e n_e u_e^2+ {\textstyle{{\hbar ^2 n_e } \over {24m_e }}}\nabla ^2\, log\lpar n_e \rpar - O\lpar \hbar ^4 \rpar } \right]$. In particular, the terms proportional to $\hbar ^2 $ in Eqs. (2) and (3) are regarded as the additional quantum terms. In Eq. (2), the quantum term ${\textstyle{{\hbar ^2 } \over {2m_e^2 }}}\nabla \lpar {\textstyle{1 \over {\sqrt {n_e } }}}\nabla ^2 \sqrt {n_e } \rpar $ is induced by a quantum pressure $P^Q= {\textstyle{{\hbar ^2 } \over {2m_e }}}\lpar \lpar \nabla \sqrt {n_e } \rpar ^2 - \sqrt {n_e } \nabla ^2 \sqrt {n_e } \rpar $, which is related to the derivative of the root density$\sqrt {n_e } $. The target ion dynamics are described as:

(4)$$\displaystyle{{\partial n_i } \over {\partial t}}+\nabla \cdot \lpar {n_i {\bf u}_i }\rpar =0\comma \; \eqno \lpar 4\rpar$$

(5)$$\displaystyle{{\partial {\bf u}_i } \over {\partial t}}+\lpar {{\bf u}_i \cdot \nabla }\rpar {\bf u}_i=- \displaystyle{{Z_i e} \over {m_i }}\nabla \Phi - \displaystyle{{\nabla P_i } \over {m_i n_i }} - {\rm \gamma} _i {\bf u}_i . \eqno \lpar 5\rpar$$

The non-relativistic beam fluid equations are:

(6)$$\displaystyle{{\partial n_b } \over {\partial t}}+\nabla \cdot \lpar {n_b {\bf u}_b }\rpar =0\comma \; \eqno \lpar 6\rpar$$

(7)$$\displaystyle{{\partial {\bf u}_b } \over {\partial t}}+\lpar {{\bf u}_b \cdot \nabla }\rpar {\bf u}_b=- \displaystyle{{Z_b e} \over {m_b }}\nabla \Phi - {\rm \gamma} _b {\bf u}_b . \eqno \lpar 7\rpar$$

The above equations are closed by adding Poisson's equation

(8)$$\nabla ^2 \Phi=4{\rm \pi} e\lpar n_e - Z_i n_i - Z_b n_b \rpar \comma \; \eqno \lpar 8\rpar$$

where $\nabla={\textstyle{\partial \over {\partial r}}}{\bf e}_r+{\textstyle{\partial \over {\partial z}}}{\bf e}_z $. Z _b = 1 is the charge state of beam ion as a proton considered here. Φ is the total electric potential generated by the breaking of local neutrality inside the target by the irradiated beam ions and the collective response of the target ions and electrons to the propagating beam ions, which is usually generated in and around the beam ion position.

Here, m _e is the electron mass, m _i is the target ion mass, m _b is the beam proton mass, e is the elementary charge, q_e, T _e, P _e and P _i are the electron heat flux, temperature, pressure and ion pressure, respectively. γ_e, γ_i, γ_b and, τ_w obtained from our previous work (Zhang et al., Reference Zhang, Song and Wang2011a), are the frictional frequencies with τ_w = 2γ_e, where the subscription of e, i, and b express the electron, ion and beam ion, respectively. The quotidian equation of state (QEOS) model (More et al., Reference More, Warren, Young and Zimmerman1988) is supplied to calculate the electron pressure P _e by finite temperature Thomas-Fermi theory and the ion pressure by modified Cowan model. In particular, the heat flux is specified by Fourier's law ${\bf q}_e=- {\rm \kappa} \nabla T_e $, where the thermal conductivity κ is defined by (Chen et al., Reference Chen, Cockburn, Gardner and Jerome1995). In the right-hand side, the second term of Eq. (2) is semiclassical statistics pressure term, and the last terms of Eqs. (2), (3), (5), and (7) are the frictional forces.

The above Eqs. (1)–(8) involve a self-consistent determination of the target ion and electron density, electron temperature, beam ion density, and the total electric field ${\bf E}=- \nabla \Phi $ with the components E _r and E _z. The Poisson's equation is solved by the successive over relaxation (SOR) method at every time step. The 2D FCT method (Boris et al., Reference Boris, Landsberg, Oran and Gardner1993) is adopted to numerically solve the Eqs. (1)–(7) by split-time integration from the initial time t = 0 when the values of all quantities are known, whose general idea is as follows. For solving Eqs. (1)–(7) a minimum amount of numerical diffusion must be added to assure stability, since these equations are subject to a numerical instability. An antidiffusion term that can ensure stability, which should not generate new maxima or minima in the solution, nor accentuate already existing extrema, is considered to remove this additional diffusion and correct the flux in FCT.

For convenience we introduce the dimensionless variables u/u ₀ → u, l/l ₀ → l, t/t ₀ → t, where u, l, and t stand for any quantity of velocity, length and time, respectively. The units u ₀ = 10⁹, l ₀ = 10⁻³ and t ₀ = l ₀/u ₀ are chosen to bestly meet numerical stability.

3. NUMERICAL RESULTS

Here, we will take the beam particle energy E _b = 1.5 MeV, intensity N = 7.9 × 10⁸, the focal radius r _b = 23 µm, and a fixed pulse time 10 ps in all figures except in specific definition. Note that during the pulse time 10 ps, the beam is consecutive.

Figure 1 shows the distribution of the electron density and temperature along the axis z (r = 0) with and without quantum effect at three different times 0.05, 0.5, and 1 ps. Here the beam parameters are E _b = 1.5 MeV, n _b0 = n ₀/10³. It is interesting to see that at 0.05 ps and 1 ps, there is no quantum effect to influence the target electron density and temperature, while at 0.5 ps the quantum effect is evident. The essence of quantum effect is to make the electrons back to equilibrium states and when the electrons are at a quasi-equilibrium state the quantum effect is very weak, therefore, as seen from Figure 1b, at 0.05 ps the quantum effect is so small (attributed to its physical nature) that can not affect the electron density and temperature. At 1 ps, the numerical characteristics of the quantum terms ${\textstyle{{\hbar }^2 \over {2m_e^2 }}}\nabla \lpar {\textstyle{1 \over {\sqrt {n_e } }}}\nabla ^2 \sqrt {n_e } \rpar $ and $ - \textstyle{{{\hbar }^2 n_e } \over {24m_e }}\nabla ^2\, log\lpar n_e \rpar $ lead the quantum effect to be failing to influence the interacting process. Overall, the quantum effect is obvious only at the very beginning the beam-plasma interaction, however, the afterward plasma evolution will definitely be affected by this process. Furthermore, we also found that including quantum effect can stabilize the numerical calculation here. Therefore, the quantum effect is considered in all following results and figures.

Fig. 1. (Color online) Electron density (a) and temperature (b) versus z at r = 0, without quantum effect (solid line) and with quantum effect (dashed line) at three different times 0.05 ps, 0.5 ps and 1 ps. Initial beam proton energy is E _b = 1.5 MeV and density is _nb0 = n ₀/10³.

Figure 2 represent the z directional (at r = 0) evolution of beam proton density Figure 2a, electron density Figure 2b, and temperature Figure 2c for three different times 10 ps (solid line), 6 ps (dashed line) and 4 ps (dotted line). E _b = 1.5 MeV. At t = 10 ps all beam ions are into the target and when t < 10 ps the only part of the beam ions are into the target. In order to express more clearly, beam ion density and electron density are normalized by initial equilibrium densities. The normalized beam proton densities are similar, while the electron density and temperature are variable at different times with a constant proton energy E _b = 1.5 MeV. As shown in Figure 2c that the pulse duration 10 ps is short enough to allow the central (at r = 0) temperature to reach its maximum value of about 5 eV without undergoing hydrodynamic expansion that would lower the peak temperature from 5 eV. Therefore, the target electrons already reach the WDM states. The WDM temperature range is roughly T _e = 1.0−20 eV for solid aluminum. One can observe that the electron density depends sensitively on the distribution of beam ions, and the electron temperature depends closely on the electron density. The electron density at the proton distributed position forms a peak, therefore, the electrons at the left and right sides of the beam proton positions move to the protons and two valleys are formed. The fluctuation of electron density increases and electron temperature increases too, as time elapses. It is interesting to find that beam proton for the first time to reach its range (z = 20 μm) at 8.2 ps, however, more protons are stopped very near the boundary z = 0 attributed to the very high electron density. Correspondingly, there is a strong interaction process between high density protons and electrons in several μm region near the boundary z = 0, which produces a maxima in electron temperature.

Fig. 2. (Color online) Beam proton density (a), electron density (b) and temperature (c) versus z at r = 0, with 10 ps (solid line), 6 ps (dashed line) and 4 ps (dotted line). Beam proton energy is E _b = 1.5 MeV.

The 1D figures described above can clearly express the spatial variation of the density and temperature in quantities, while 2D spatial density and temperature spectrums are displayed more visually in Figure 3. Here the time is fixed to 10 ps and E _b = 10 MeV. Beam protons impinging on the target let the electrons move to them, which leads to changes in the electron density as well as in the temperature. The beam proton density has its maximum value at the point (r = 0, z = 0), that is, the center of the focus of the incident ion beam, and electron density and temperature have their maximum values near the point. The temperature takes values of T _e = 0.5−4 eV in all heated region where WDM is generated, which is in favor of important experiments on HEDP.

Fig. 3. (Color online) Spacial distribution of beam proton density (a), electron density (b), and temperature (c) at 10 ps. E _b = 10 MeV.

To understand the time and spatial evolutions of the beam-target heating, in Figure 4, we illustrate beam proton density Figure 4a, electron density Figure 4b, and temperature Figure 4c at r = 0. They demonstrate a clear time and space variations and a larger volume of target has been disturbed as time increases. A highest n _b produces a maximum n _e, as well as a maximum T _e, since the beam protons in the target pull the electrons moving to them and also the hydrodynamic expansion is neglected in a picosecond time scale. The number of protons entered in the target is gradually increased during the pulse time 10 ps, as a result, the changes of electron density and temperature increase and the changes can reach maximum values at 10 ps.

Fig. 4. (Color online) Space (z) and time (t) resolved beam proton density (a), electron density (b), and temperature (c) at r = 0. E _b = 1.5 MeV.

To analyze the average influence of the beam and the collective response of the target, an average electron temperature is calculated as a space average in 50 μm × 50 μm (r,z) plane. The influences of N 5(a), r _b 5(b), and E _b 5(c) on the temporal evolution of this average T _e are shown in Figures 5. In Figure 5a one observes a significant increase in the average T _e as N increases. As r _b decreases in Figure 5b, the distribution volume of the total beam particles becomes small leading to a rise in n _b, so that the electron temperature as well as the average T _e increases distinctly. The effect of increasing E _b, seen in Figure 5c, introduces a reduction in T _e, since the beam-target interaction time becomes shorter for higher E _b. Note that the electron temperature of the 50 μm × 50 μm heated target region can reach 0.2 eV for E _b = 1.5 MeV, N = 7.9 × 10⁸, and r _b = 23 μm.

Fig. 5. (Color online) Electron temperature is averaged by a 2D volume of 50 μm × 50 μm at (r,z) plane, defined as Average T _e. Time evolution of average T _e with r _b = 23 μm, E _b = 1.5 MeV (a), E _b = 1.5 MeV, N = 7.9 × 10⁸ (b), and N = 7.9 ×10⁸, r _b = 23 μm (c).

4. SUMMARY

We study a cold solid target heated by a proton beam pulse of particle energy E _b, intensity N and focal radius r _b of transverse Gaussian distribution, with a pulse time 10 ps. The dynamics of target electrons are presented by the 2D QHD theory and the target ions and beam ions dynamics are examined by a 2D classical hydrodynamic model. Numerical simulations for this beam-target heating are carried out by employing the FCT method. As the beam ions penetrate into the target, the target electrons are collectively polarized and the target is isochorically heated in a picosecond time scale. It is seen that the target is heated to several eV, as a result, WDM is generated in the target on a picosecond time scale. Therefore, one can access the very interesting WDM states in this work.