Neutron sources have been widely applied in numerous fields, ranging from detection technology1–3 to medicine4 and laboratory astrophysics,5–7 owing to the neutron’s unique features of electrical neutrality and deep penetration capability. Currently, photonuclear reactions (γ, n)8,9 and ion-acceleration-induced nuclear reactions10,11 are the two main methods of neutron production. Compared with traditional spallation neutron sources based on large particle accelerators,12 laser-driven neutron sources (LDNSs) have a number of advantages, being able to generate neutron pulses with high brightness, compactness, and short duration at relatively low cost.13–19 Typical LDNSs are based on nuclear fusion20 or ion acceleration, including Coulomb explosions in cluster targets21 and the pitcher–catcher model.11,22 In the case of sources based on ion acceleration, the ions are accelerated via target normal sheath acceleration (TNSA),23 breakout afterburner acceleration (BOA),11,24 or shock acceleration.25 Currently, LDNSs are still limited by low conversion efficiency, large energy spread, and wide divergence angles. Also, they are often produced with a duration beyond dozens of picoseconds.

Bulk acceleration in foam targets is another LDNS scheme.26 It has been verified theoretically and experimentally26,27 that a remarkable enhancement of production efficiency is possible with this approach, since the majority of deuterons are accelerated inside the target rather than at its surface, and they collide with both accelerated and background deuterons. The produced neutrons are mainly isotropic, since the density of a foam target is randomly distributed with a certain degree of isotropy. On the other hand, periodic plasma density distributions or plasma gratings can be formed via optical traps (OTs),28–31 produced using two counter-propagating laser pulses. Ion acceleration inside a plasma grating can be realized in the whole plasma volume, since the accelerating electric fields are generated inside the target through modulation of the density distribution. This may lead to a high flux of ion acceleration and subsequently high-flux neutron generation.

In this work, we propose a new scheme for an LDNS to generate quasi-monoenergetic high-flux neutrons with short temporal duration and high anisotropy. This is achieved via the interaction of intense laser pulses with a local near-critical density (NCD) deuteron plasma. In this scheme, a pair of counter-propagating laser pulses at moderate intensity first interact with the plasma to form plasma gratings and meanwhile cause pre-acceleration of deuteron ions. Then, the target is irradiated by a relativistic intense laser pulse to boost the ion energy and trigger nuclear reactions D(d,n)³He to release neutrons. Both particle-in-cell (PIC) simulations and Monte Carlo post-processor simulation, as well as a theoretical analysis, are carried out to illustrate these processes. The potential of the proposed scheme for LDNS is demonstrated, with its merits including high peak fluxes with compact target size, quasi-monoenergetic spectra, relatively small angular distributions, and sub-picosecond duration. The remainder of the paper is organized as follows. In Sec. II, the model of pre-acceleration of deuteron ions in OTs is described. In Sec. III, ion energy boosting with a relativistic intense laser is shown, and a theoretical model is presented. In Sec. IV, neutron generation is calculated. The paper is concluded in Sec. V.

The production of a neutron source from an OT involves three stages: formation of the OT and plasma density gratings for pre-acceleration of deuteron ions, ion energy boosting by a relativistic intense laser, and neutron generation via DD fusion. An illustrative plot of the scheme is shown in Fig. 1. First, to form an OT and produce plasma density modulation, two counter-propagating laser pulses with the same frequency and amplitude overlap inside the NCD plasma. During this overlapping of the two pulses, a standing wave is formed, leading to a modulated laser intensity distribution with a certain period in the plasma.28 This intensity modulation creates a strong ponderomotive force, which pushes plasma electrons rapidly to form an electron density grating. Once this grating has been formed, strong charge separation fields are created, and these drag plasma ions.29,30 This results in the formation of a plasma grating, which is almost charge-neutral. The thermal pressure resulting from the existence of a thermal potential expels electrons, and the resulting strong self-generated electric field efficiently reflects ions in the opposite direction. During this stage, deuteron ions are pre-accelerated to the tens of keV level. Second, a relativistic laser pulse is injected into the plasma grating. As the grating is an underdense plasma, a large number of ions are further accelerated to the MeV level in both the forward and backward directions. Third, DD fusion produces neutrons, predominantly in the forward and backward directions.

In the following, we describe the pre-acceleration of deuteron ions. We assume that two identical laser pulses counter-propagate along the x axis in a uniform plasma slab with density n₀, with the laser pulses being described by normalized vector potentials a = a₀ cos(kx ∓ ω₀t). Here, the peak laser amplitude a₀ = eE₀/m_eω₀c, $k=k_0{(1-n_0/n_c)}^{1/2}$, k₀ = ω₀/c, ω₀ is the laser frequency, and the corresponding critical density $n_c=m_e{\omega }_0^2/4\pi e^2$. In the overlap region, a periodic high-intensity modulation and a corresponding ponderomotive potential are then created, which can be described by ${\varphi }_p=(m_e{c}^2/2e)a_0^2\cos (2kx)$. This leads to a rapid response by the electrons, which become trapped in the intensity valleys and bunched. Owing to the much greater mass of the ions, a large electrostatic field will form between the bunched electrons and background ions. The electrostatic force then drags ions into the valleys. Finally, a quasi-neutral OT, also called a plasma density grating, is formed. However, such a spatially periodic modulation of the plasma density is not stable, especially when the electron temperature grows.32 With the continuous accumulation of electrons in the valleys, the thermal pressure will become dominant and expel the electrons from the density peaks to the background, leading to localized non-neutrality and the formation of electrostatic dipole fields. Such self-generated fields in the OT will then reflect and pre-accelerate the remaining unbunched deuterons to both sides.

It should be mentioned that the laser intensity for the formation of an OT must be well below unity, a₀ < 1, and the plasma must be below the critical density. In this case, the laser pulses can propagate through the plasma and overlap inside the target without significant reflection. Such a laser–plasma interaction configuration is accessible with multiple tunable high-power lasers.33,34 The NCD plasma can be obtained by ablation of a thin-foil target or with high-density gas nozzles.

To illustrate the formation of an OT and pre-acceleration of deuteron ions, a one-dimensional (1D) PIC simulation using EPOCH35 is conducted. In this simulation, we take the two laser pulses to have the same intensity and wavelength. The laser profile is given by a = a₀ sin²(πt/2τ₀), where τ₀ = 30T₀ is the pulse duration in full width at half maximum (FWHM), T₀ ≈ 3.5 fs is the laser cycle, and a₀ ≈ 0.14, corresponding to an intensity I = 2.45 × 10¹⁶ W/cm² for a laser wavelength λ₀ = 1.057 μm. The target is a fully ionized deuteron plasma with initial density n₀ = 0.3n_c, located within 5λ₀ ≤ x ≤ 35λ₀. Only electrons and deuterons (m_e:M_D = 1:3672) are considered, and it is assumed that their initial temperatures are 1 and 10 eV, respectively. The simulation box size is 40λ₀, with a spatial resolution of 0.002λ₀ and 200 macroparticles per cell for both species. The results are shown in Fig. 2.

When the two laser pulses have fully overlapped, the electrons start to accumulate. In particular, the electrons will first form a density grating, which causes charge separation as shown in Fig. 2(a) at t = 50T₀, since the high mass of the ions means that it is difficult for the laser ponderomotive force to push them in a short time. Note that a closeup of only three periodic layers is shown here. Once the charge separation field has been established and the two laser pulses cross each other, ions will be moved by the charge separation field to form an ion density grating. During this process, the maximum velocity reached by the ions is given by$$\frac{v_{i,\max }}{c}=\frac{4kc{\tau }_0}{{(1+k/k_0)}^2}\frac{m_e}{M_D}{a}_0^2.$$(1)Meanwhile, electrons will follow the ion motion. Finally, a quasi-charge-neutral OT will be formed as shown in Fig. 2(b). On the other hand, as the OT grows, the peak electron density n_e and temperature T_e at the density peaks will become very high. For simplicity, we assume that this process is adiabatic, and we then have $T_e\sim {(n_e/n_0)}^2{T}_{e,0}$, where n₀ and T_e,0 are the initial electron density and temperature, respectively. This implies that the thermal potential ${\varphi }_{th}=3T_{e,0}{n}_e^2/2e{n}_0^2$ grows and that electrons in the density peaks are expelled by thermal pressure. Given that the thermal velocity of the electrons is much larger than that of the ions, the ion density peak will be narrower and higher than the electron density peak, as can be seen from the blue dashed line and red dotted line, respectively, in Fig. 2(c). This leads to the formation of a new bipolar charge separation field around each plasma density peak, corresponding to the E_x field observed after t = 160T₀ in Fig. 2(a). The critical point at which the self-generated field appears can be obtained by equating the kinetic energy gained by the ions from the early E_x field caused by the laser ponderomotive force to the thermal potential barrier:$$\frac{n_{e,\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}{n_0}=\sqrt{\frac{M_D{v}_{i,\max }^2}{3T_{e0}}}.$$(2)On the basis of Eqs. (1) and (2), the critical point for electrons is found to be n_e,crit = 2.42n_c, which matches the simulation results. Such a self-generated field will resist and stop the incoming ions, causing the accumulation of deuterons in each sides, which corresponds to the double spikes in the ion density distribution in Fig. 2(d). As the imbalance between the electron and ion densities increases, the electrostatic field will increase to 10¹¹ V/m and then begin to reflect and accelerate ions in reverse directions to both sides, causing the “X” pattern in ion phase space shown in Fig. 2(e). The existence and properties of the self-generated field illustrated here indicate that the assumption of quasi-neutrality of the OT during its evolution that was made in previous theoretical work is no longer valid. The presence of this field also provides an opportunity for the reflected deuterons to collide with each other.

The description above is of the early bulk acceleration in OT. However, although the self-generated fields can pre-accelerate and reflect most ions, considering that the strength of OT lasers is limited, the final collision velocity of deuterons will be much lower than that required for a high nuclear reaction cross-section. To enhance the ion energy, we introduce additional intense lasers to boost the energy of reflected ions. Since the density is NCD (∼γn_c) and the width of each layer in the OT is comparable to the laser skin depth, with the help of relativistic induced transparency (RIT), the boost lasers can continuously penetrate through each layer and accelerate electrons and ions there, which makes the OT a perfect medium for ion acceleration via mechanisms like breakout afterburner (BOA)36,37 and RIT24 acceleration.

Both 1D and 2D PIC simulations are performed to illustrate the ion energy boost. The boost lasers have wavelength and intensity of 1.057 μm and 10²⁰ W/cm², respectively, and are horizontally injected from both sides. The boost laser duration is 26 fs, with the same profile as the OT lasers. The laser spot size is 10 μm in FWHM, with a Gaussian distribution in the 2D PIC simulations. The results are shown below. In Fig. 3, we verify the existence of OT self-generated fields in the 2D case and compare the maximum velocities with and without energy boosting. With the introduction of the boost lasers, the electrostatic fields between OT layers increase to 3 × 10¹² V/m, as the electrons are further expelled by the intense lasers. The reflected ions are then accelerated from 5 × 10⁻³c to 4 × 10⁻²c. To further illustrate the acceleration features in 2D, we show in Fig. 4 the boosted ion energy angular distributions based on the transverse and longitudinal velocities. It is found that the ions are accelerated mainly along the directions of laser propagation. In addition, the ion energy spectra show that compared with the transverse direction (90°), the majority of the high-energy ions are emitted longitudinally (0°) and these are also the main contributors to neutron generation. Meanwhile, an arched platform is found in the 0° spectrum around 2 MeV, corresponding to the ion acceleration to a velocity of 0.04c by the boost lasers. Although the ions can be accelerated transversely in the 2D geometry, owing to the transverse ponderomotive force of the laser pulses, the ion transverse acceleration is not remarkable in our model. This can be attributed to two causes. First, unlike a foam target, no transverse structures or density distributions are involved in the 1D plasma grating, which reduces the transverse acceleration effects of both electrons and ions. Second, given that the pre-acceleration effects are small, as is the distance between each layer (∼λ₀), the time of ion acceleration along the transverse direction before collision can be ignored, and so the ions remain to be mainly accelerated longitudinally before they collide with each other and release neutrons. On the basis of the above analysis, it is clear that the ion acceleration features of our model in the 2D case will be comparable to those in the 1D configuration, as found in our simulations.

Our 1D simulation results illustrate even more clearly how the accelerated ions evolve as they counter-propagate between each layer. Figure 5(a) corresponds to the first collision of accelerated ions. The maximum relative velocity is about 0.07–0.08c, and at the collision point the ions are highly compressed, with a large accumulation of density around 10n_c. As shown in Fig. 5(b), secondary counter-propagating collisions can be triggered at later times, as the E_x field (red dashed line) will continuously drag ions and form a plateau feature in phase space. These results also illustrate that there will be two contributions to the final neutrons produced in OT via DD collisions: from the counter-propagating ions with high energy level and from the interaction between boosted ions and background ions.

To better estimate this enhanced acceleration in the OT, we construct a toy model that gives a semiquantitative description. Given that the laser can penetrate the OT, the dynamics of electrons will be directly influenced by the laser pulse. Therefore, we adjust the electron temperature equation $T_{\text{hot}}^{\prime }$ by including the effects of both the thermal energy E₀ and the ponderomotive potential Φ_p:$$T_{\text{hot}}^{\prime }=m_e{c}^2\left(\sqrt{1+(1-\mathrm{\Gamma })a^2}+\sqrt{1+\mathrm{\Gamma }{a}^2}-2\right).$$(3)Here, $\mathrm{\Gamma }={(1+{\pi }^2{\xi }^2)}^{-1}$ is a transmission coefficient38 that determines the relative contributions of thermal motion E₀ = m_ec²(γ_t − 1) and the direct influence of the laser ponderomotive potential Φ_p = m_ec²(γ_p − 1), where ${\gamma }_t=\sqrt{1+(1-\mathrm{\Gamma })a^2}$ and ${\gamma }_p=\sqrt{1+\mathrm{\Gamma }{a}^2}$ are the electron kinetic energies associated with the incident and transmitted laser pulses.39,40ξ = (d/λ)n/γn_c is a dimensionless parameter related to the optical thickness and relativistic effects,39 with d being the layer width and n the density. It is obvious that when the target layer is too thick or of too high a density, ξ will tend to infinity, and Eq. (3) will reduce to the expression for the case of TNSA. We then express the corresponding electrostatic field while still in the theoretical framework of TNSA41 for simplicity: $E(x)=2k_B{T}_{\text{hot}}^{\prime }/e(x+\sqrt{2e_N}{\lambda }_D)$. Given that for most cases of bulk acceleration, the layer spacing L is small and the interaction is fast, we assume that E_x at the middle point between each layer can be used to approximately calculate the continuous work that the inhomogeneous electrostatic field does on ions before they collide. The interaction time can then be approximated as that required for the ion acoustic wave to expand within a distance of half of the layer spacing, namely, Δt = Δx/C_s. Finally by conservation of momentum, the velocity that deuterons gain in the OT under the action of the boost lasers before they collide can be expressed as $v_{i,\max }=e{E}_x(T_{\text{hot}}^{\prime })\mathrm{\Delta }t/M_D$.

We then apply the toy model to illustrate the dependence of the boosted ion velocity on the bulk acceleration in the OT, which is of significance for parameter optimization in neutron generation. A series of 1D simulations are performed as shown in Fig. 6. We first change the time difference between the injection of the boost lasers and the end of the OT laser, $\mathrm{\Delta }t=t_{\text{boost}}^S-t_{\text{OT}}^E$, to modulate the interacted OT density n_peak. In Fig. 6(a), when the boost lasers are injected earlier, the OT still grows, with a low density value of 2n_c. As the time difference Δt increases, the density peak approaches its maximum value of about 11n_c, and starts to decrease as the OT decays. The corresponding boosted velocities in both the simulation and the toy model exhibit the same trend, with the maximum velocity reaching 0.04c. The next parameter to be changed is the layer spacing L. The distance between each layer of the OT is given by $L=2\pi /2k={\lambda }_{\text{OT}}/2{(1-n_0/n_c)}^{1/2}$. As shown in Fig. 6(b), by increasing the initial density of the plasma n₀, it is possible to enlarge the layer distance L from 0.5λ_OT to 1.12λ_OT, and the velocity also increases, to 0.039c since the acceleration length grows. Both the toy model and simulation results in Fig. 6(c) show that the boosted velocity increases as the intensity of the boost lasers grows. However, the velocity appears to become saturated around 0.05c.

An interesting phenomenon can be observed in Fig. 6(d) when we fix a_OT, a_boost, and the normalized initial density n₀/n_c(λ_OT) = 0.8. As the OT laser wavelength λ_OT increases, which leads to an increase in the layer spacing L, i.e., the acceleration length, we might expect the final boosted velocity to grow as well. However, both model and simulation results indicate that for bulk acceleration in the OT, the maximum velocity remains the same, at around 0.05c. This is because although the layer distance is increased by λ_OT, the critical density decreases as $n_c\propto 1/{\lambda }_{\text{OT}}^2$, which decreases both of the plasma density values n₀ and n_peak. This “velocity locking” revealed by both the toy model and simulation suggests that it should be possible to use a shorter wavelength in OT lasers so that the allowed maximum initial density can be increased. This leads to a larger collision density for triggering fusion reactions, and meanwhile the OT lasers can still interfere freely inside the plasma, with the final collision velocity and corresponding neutron pulse duration remaining unchanged as long as the boost laser and target scale are fixed.

A PIC code from EPOCH and a Monte Carlo (MC) post-processor code that we have constructed are applied to calculate the generation of neutrons in the OT by bulk acceleration. On the basis of information about D⁺ positions, momenta, and densities given by the PIC simulation, the MC code divides the macroparticles into different cells. All of the particles in the same cell are chosen in pairs randomly to calculate the corresponding volumetric neutron production rate$$R=\frac{n_1{n}_2}{1+{\delta }_{12}}\sigma v,$$(4)where n₁ and n₂ are the number densities of two reacting deuteron macroparticles, $v=|{\bar{\mathbf{v}}}_1-{\bar{\mathbf{v}}}_2|$ is their relative velocity, δ₁₂ is the Kronecker symbol, and σ is the total cross section. The cross section σ is the integral of the differential cross section dσ(θ)/dΩ over solid angle, which can be obtained via interpolation using a data from a nuclear reaction database for different ion collision energies E_r and neutron emission angles θ_n referred to the relative velocity direction of colliding ions. The neutron energy is then calculated as10$$E_n=\frac{1}{8}{E}_r{\left[{\left(2+\frac{19.6\mathrm{M}\mathrm{e}\mathrm{V}}{E_r}+{\cos }^2{\theta }_n\right)}^{1/2}+\cos {\theta }_n\right]}^2.$$(5)Equation (4) indicates that the quality of a neutron source can be improved by increasing the deuteron collision velocity v_c or density n. On the basis of the above analysis and the toy model, the optimized parameters that we choose for simulation are as follows. The two OT lasers have low intensity a_OT = 0.19 and short wavelength λ_OT = 0.5 μm, with duration 0.25 ps. Two boost lasers are then injected from both sides of the plasma for better symmetry in neutron generation, with regular laser wavelength λ_boost = 1.057 μm, high intensity I_boost = 6 × 10²⁰ W/cm², and short duration 70 fs. The time difference between the two laser groups is Δt = −125 fs. The simulation box length is 30 μm, with plasma located within 2.5–27.5 μm, and with high initial density 0.8n_c(λ_OT) and low temperature for larger OT density.28 To estimate the realistic neutron number in a 3D situation, we assume that the produced neutron volume V has cylindrical symmetry in the 1D case, with laser spot radius 10 μm.

In Fig. 7(a), the black line shows the temporal evolution of the neutron production rate P = ∫RdV (black line), which has two spikes with values of 4.3 × 10¹⁸n/s and 3.7 × 10¹⁸n/s, respectively. The first spike marked by the red arrow at 0.5 ps corresponds to collisions of energy-boosted deuterons, which counter-propagate. The spatial distributions of the volumetric production rate R at this moment is shown in the insert, where each peak of R corresponds to a collision point with high compression density and high relative velocity. This is exactly the same situation as shown in Fig. 5(a), where the majority of counter-propagating ions reflected by the OT field have been accelerated to high energy by the boost lasers, and the density increases as the ions begin to accumulate at the collision points. After that, the encountered ions are passing through each other and mainly interact with the stationary ions in the background, causing an instant drop in the production rate as the relative velocity and collision density decrease, until a turning point appears, as marked by the blue arrow in Fig. 7(a). At this point, a secondary intense nuclear reaction starts to be triggered, corresponding to collisions between the neighboring counter-propagating ions, which can be explained by reference to Fig. 5(b). Compared with the first main collisions, the second ones have a relatively low collision density, since the whole OT is strongly disturbed by the intense boost lasers. Meanwhile, the accelerated ions have a velocity plateau in phase space, owing to the electrostatic field distribution, which leads to a lower but wider secondary spike in the production rate. There are no more subsequent neutron peaks, because of the total collapse of the OT. The final neutron yields N_n reaches 1.3 × 10⁶ within 800 fs, as shown by the red dashed line in Fig. 7(a), which indicates that an OT neutron source has a high growth rate.

As well as the two peaks in the production rate, we can see that its time duration is δT_n = 390 fs in FWHM, which is much shorter compared with other LDNS schemes. To further verify the results for the duration, we check the neutron energy spectra shown in Fig. 7(b). Given that our PIC simulation is 1D3V, we can obtain the velocities of particles in the transverse directions to calculate the emission angle. The simulation shows that at θ = 0°, the neutron spectrum has a monoenergetic peak ranging mainly from 1.8 to 4.6 MeV. For θ = 90°, the spectrum has a sharper distribution but a narrower energy range. The insert in Fig. 7(b) shows the neutron spectrum dN/dE integrated over emission angles between 0° and 180° and the energy range from E_nl = 1.64 MeV to E_nh = 8 MeV. The neutron source pulse duration can then be estimated as follows:42$$2\delta T_n(\mathrm{p}\mathrm{s})=73d(\mathrm{m}\mathrm{m})\left(\frac{1}{\sqrt{E_{nl}(\mathrm{M}\mathrm{e}\mathrm{V})}}-\frac{1}{\sqrt{E_{nh}(\mathrm{M}\mathrm{e}\mathrm{V})}}\right).$$(6)For a target size with d = 25 μm, the corresponding neutron full width duration is about 2δT_n = 0.79 ps, which is consistent with the simulation results, as Fig. 7(a) shows. The short duration can be explained as follows. First, given that our scheme is more like a mutual pitcher–catcher process because of the periodic acceleration and collision inside the plasma, the OT neutron source can achieve nuclear reaction more efficiently, even though the plasma scale is quite small. The latter directly determines the neutron pulse duration. Second, the pre-acceleration process in the second step reduces the response time of the ions when they are further accelerated by the boost lasers, and so the reactions are more rapid. Third, the boost lasers are intense enough that the OT structures inside the NCD plasma are disturbed significantly, causing its total collapse and the cessation of nuclear reactions in a short period of time.

Another interesting phenomenon observed in Fig. 7(b) is that neutrons emitted at 90° have a narrower spectrum, and the central energy has a left offset compared with emission at 0°. To illustrate this, we further output neutron spectra for emission angles from 0° to 90°, as shown in Fig. 8. It is found that as the emission angle increases, the central energy (black line) and spectral range (blue lines) decrease. The neutron production (red line), however, grows when θ is beyond 60°. These results can be explained as being due to the two contributions to neutron generation in the OT, as shown in Fig. 5. For the collisions between counter-propagating ions that are boosted synchronously to comparable energy in opposite directions, the corresponding generated neutrons are emitted randomly in the full 4π solid angle, since the collisions are elastic and the central mass velocity of collided ions is almost zero. Our simulation shows that the neutrons generated in radial directions from 60° to 90° maintain a large peak production with narrow energy spectrum. For the collisions between boosted ions and stationary ions in the background, given that their relative velocity $d\bar{v}$ has a certain directionality and the boosted ions are continuously accelerated, the corresponding neutrons are mainly emitted along the x axis with a wider energy distribution. The offset of the central energy in the neutron spectra can be explained directly by Eq. (5), which indicates that θ increases as the neutron energy E_n decreases. These results illustrate unique features of our OT neutron source: by choosing different angles, we can obtain different-width neutron spectra with different central energies, while maintaining the peak neutron flux at a high level. This could be interesting for applications. Besides, it also provides an unique perspective for further experimental verification of the OT neutron source.

Figure 7(c) plots the angular distributions of the total neutron yield. It is found that the OT neutron source has a good anisotropy and high yields, with 2.6 × 10⁵n/sr forward neutrons emitted at 0°, and the total number of generated neutrons is N_tot = ∫Pdt = ∫(dN/dΩ)2π sin θdθ ∼ 10⁶. In addition, as shown in the insert of this figure, the high symmetry in neutron emission promises higher-efficiency production of neutrons in both forward and backward directions. The proportion of neutrons emitted along the axial direction with this scheme is much larger than with previously proposed LDNS schemes such as the pitcher–catcher scheme or photonuclear reactions. This can be explained by the fact that ion collisions occur mainly along the x axis in our scheme, owing to the bulk acceleration mechanism, and the OT structures ensure that the energy converted from the lasers mainly leads to longitudinal acceleration of ions.

In addition, we have examined the main results on neutron generation obtained from 2D PIC simulations. To better compare the results obtained in the 1D situation with a uniform transverse distribution within a volume of diameter 20 μm, in 2D, we consider a laser with a super-Gaussian profile and a diameter of 20 μm. It is found that in 2D, as shown in Figs. 9(a) and 9(b), the spatial distributions of the neutron volumetric production rate R at the ion collision time 0.5 ps has similar features to the 1D case, with periodic multiple peaks along the longitudinal direction at around 1 × 10²⁸ cm⁻³ s⁻¹, while also exhibiting transverse distribution features. We show the angular distribution of the neutron yield at the collision time in Fig. 9(c). In the 2D situation, the angular distribution maintains good anisotropy, with the neutrons being emitted mainly along the laser propagation direction, owing to the dominant longitudinal acceleration of ions, as discussed in Sec. III. The ratio of the neutron yields at 0° and 90° is 3.1, which is close to the 1D result. The maximum neutron yield along 0° at the collision time is dN/dΩ ∼ 900, which is lower than in the 1D case: dN/dΩ ∼ 2000. The maximum peak production rate in the 2D case is about 2.0 × 10¹⁸ s⁻¹, which is half of that in the 1D case. On the basis of this analysis, we believe that the basic features of our model in 2D are similar to those in the 1D case. Although the final neutron production in 2D is smaller than that in 1D, this discrepancy is acceptable, since the values of the volumetric neutron yields are comparable for the two cases.

On the basis of the dependence of the ion collision velocity on bulk acceleration in the OT described in Sec. III, we further compare neutron generation under different laser parameters. Our 1D3V simulation results show that for two boost lasers with lower intensity of I_boost = 1 × 10²⁰ W/cm², the corresponding production rate decreases to about 1 × 10¹⁸n/s. This is due to the reduced laser penetration in the high-density OT as the intensity drops, with the consequence that the ions inside the target fail to be fully boosted. Similar results have been observed when only one boost laser is injected, in which case the production rate shrinks by almost half compared with the situation with two boost lasers, while the duration is extended. In addition, on the basis of the “velocity lock” pattern found in OT bulk acceleration, we further compare the final neutron results for different wavelengths of the OT lasers, with the normalized laser amplitudes for both the OT and boost lasers and the normalized plasma density remaining the same. As shown in Fig. 10, as the OT laser wavelength λ_OT decreases, the maximum production rate P_max increases rapidly, while the neutron source duration δT_n fluctuates slightly around 400 fs. This is because as λ_OT decreases, the corresponding plasma initial density 0.8n_c increases, since the critical density limit is enhanced by $n_c\propto 1/{\lambda }_{\text{OT}}^2$. On this basis, the volumetric neutron production rates should be increased by X⁴ times if λ_OT shrinks X times, according to Eq. (4), which implies that $R\propto n^2\sim n_c^2\sim {\lambda }_{\text{OT}}^{-4}$. This is exactly the same trend found in our simulation. On the other hand, the duration is not dramatically influenced by the OT laser wavelength. This is due to the “velocity lock” feature, as we mentioned above. Once the boosted velocity for collision ions and the whole length of the target have been fixed, the final neutron source duration is also locked. It should be mentioned that the wavelength of the boost laser in the case of λ_OT = 0.4 μm is 0.8 μm, for better penetrability in the high-density plasma grating. Such parameters are experimentally accessible, and indeed the former wavelength can be considered as corresponding to second-harmonic generation (SHG) of the latter.

Finally, we compare our neutron source results with those of previous studies with regard to a few parameters, as shown in Table I. Overall, one of the most remarkable features of the OT neutron source is its extremely short duration, down to the femtosecond scale, while maintaining a high production rate and high anisotropy. Compared with the pitcher–catcher scheme for LDNS, the OT neutron source is based upon bulk acceleration, which is able to more efficiently accelerate ions inside the target rather than at the target surface. In addition, nuclear reactions occur directly within the bulk target, which also reduces the source scale and neutron emission duration.

We have proposed a theoretical scheme for a compact LDNS based on laser interaction with near-critical-density deuteron plasma. It is consists of three steps. First, deuteron ions are pre-accelerated to the 10 keV level in an OT or plasma grating formed by two counter-propagating laser pulses (called OT lasers in this paper) at moderate intensity. Second, a huge number of deuteron ions in the OT are further accelerated to the MeV level via bulk acceleration by another one or two boosting lasers (called boost lasers) at relativistic intensity, where ion acceleration occurs in both the forward and backward directions. Third, nuclear reactions are triggered by accelerated deuteron ions, leading to bursts of short-pulse neutron emission with high flux dominantly along the boost laser directions. Numerical simulations have been carried out to illustrate these processes. By constructing a toy model, we have optimized parameters for bulk acceleration in the OT and neutron emission. In particular, it has been found that the maximum ion energy is not sensitive to the OT laser wavelength, provided the relative plasma density n₀/n_c is fixed. This feature of ion velocity locking suggests that it should be possible to use two-color lasers, i.e., with the OT lasers and boost lasers having different wavelengths. In a typical example, with a plasma size of $\sim 30\mu$m, a plasma density of 0.8n_c, and boost lasers at an intensity of $\sim 10^{20}{\mathrm{W}/\mathrm{c}\mathrm{m}}^2$, it is possible to achieve a high flux ($>10^{18}n$/s), short duration (∼400 fs), and good directivity. Such compact ultrafast neutron sources may find applications in neutron radiography, nuclear research, and diverse scientific fields.

Acknowledgment. This work is supported by the National Natural Science Foundation of China (Grant Nos. 11991074, 11975154, 12135009, 12005287, and 12225505) and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA25050100).