Numerical Study of Proton Impact Rates on a Solid Target in Pyroelectric Fusion: Effects of Electrode Geometry and Crystal Heating Temperature

Abstract

Numerical simulations were conducted for a pyroelectric fusion device with three electrode geometries to investigate proton generation, acceleration, and impact with a solid ¹¹B target. Evaluation of the electric field distributions showed that the maximum potential scales with the crystal temperature difference as ϕ_max ≃ 238 ΔT [kV]. Analysis of electron-impact ionization clarified the spatial regions where protons are generated and begin to accelerate, revealing that a considerable fraction of protons is laterally deflected by the electric field structure and fail to reach the target. The proton impact rate was systematically examined by varying the electrode geometry and heating conditions. The results indicate that the impact rate is maximized at moderate temperature differences, whereas excessively strong electric fields shift the ionization region outward and reduce proton delivery to the target. Further analysis suggests that rapid electron acceleration near the electrode shortens the residence time within the energy range favorable for ionization, thereby suppressing proton generation under high heating conditions. These results demonstrate that proton impact efficiency in pyroelectric fusion devices is strongly governed by the interplay between electrode geometry, ionization location, and operating temperature. Appropriate control of these parameters is therefore essential for improving target-directed proton irradiation.

1.  Introduction

Pyroelectric crystals are known to generate high electric potentials due to temperature-induced variations in their polarization state. This unique property has garnered significant interest in their application in compact particle accelerator systems. Commercialized pyroelectric X-ray sources are already in use in analytical techniques such as X-ray absorption spectroscopy and X-ray fluorescence analysis [1]. Furthermore, these crystals are capable of accelerating ions upon heating. Naranjo et al. first demonstrated deuterium-deuterium (D-D) fusion reactions by accelerating deuterium ions using a lithium tantalate (LiTaO₃) crystal and directing the resulting beam onto an erbium deuteride (ErD₂) target [2]. Following this breakthrough, numerous studies have investigated the feasibility of pyroelectric fusion, particularly as a means to develop compact neutron sources [3–6].

Among the various nuclear fusion reactions under investigation, the proton-boron (p-¹¹B) reaction has attracted considerable attention due to its advantageous fuel-handling characteristics and aneutronic nature. However, the fusion cross-section for p-¹¹B is relatively low—approximately 1.4 barns at a proton energy of 675 keV [7]—compared to the 5–6 barns observed for deuterium-tritium (D-T) fusion at 60–70 keV, making ignition substantially more difficult. To address the engineering challenges associated with this reaction, compact experimental platforms are essential for evaluating conditions such as the reaction environment and energy conversion processes. Pyroelectric crystal-based beam accelerators represent a promising candidate for such platforms.

This study conducts numerical simulations of proton acceleration using pyroelectric crystals, with the following three objectives: (i) to propose and evaluate a mechanism for proton generation within a hydrogen-filled device, (ii) to analyze the incidence rate of protons on a solid ¹¹B target via charged-particle trajectory calculations, and (iii) to develop a method for tracking the stopping points of individual charged particles. These simulations aim to identify key issues relevant to the future development of compact proton–boron fusion accelerator systems.

2.  Numerical Procedure

The simulation domain consists of a cylindrical metallic chamber filled with hydrogen gas at 300 K and 1 Pa. All metallic surfaces, including the chamber wall, the lower face of the pyroelectric crystal, and the electrode surfaces, are treated as ideal conductors and assigned the boundary condition $\varphi$ = 0, where $\varphi$ denotes the electric potential. In contrast, the upper surface of the crystal carries a temperature-dependent polarization charge density $\sigma =\Gamma \Delta T$ [8], where $\Gamma$ is the pyroelectric coefficient and $\Delta T$ is crystal heating temperature. The surface of a pyroelectric crystal carries bound charges originating from its spontaneous polarization, even in the absence of any temperature variation. When a metallic electrode is attached to the crystal surface, the electrode enforces an equipotential condition, and a compensating negative charge distribution appears on the electrode to balance the polarization charge. Upon heating the crystal, the magnitude of the spontaneous polarization decreases, whereas the compensating negative charge on the electrode does not immediately dissipate. As a result, the electrode temporarily retains an apparent excess of negative charge. To neutralize this imbalance, positive free charges in the surrounding environment are attracted to and accumulate on the electrode surface. This sequence of charge redistribution constitutes the mechanism by which a high electric potential associated with the pyroelectric polarization effectively manifests on the metallic electrode. A schematic of this process, including the Gaussian surface used in later analysis, is shown in Fig. 1.

Fig. 1.  Surface charge formation induced by pyroelectric polarization change.

In the evaluation of the electric field using a Gaussian surface, the surface area of the metal electrode generally differs from that of the crystal’s upper face. It is therefore necessary to treat these areas separately. In the present study, this difference in surface area is explicitly considered when computing the normal component of the electric field, leading to the following expression:

  

$$\begin{array}{ccc}\hfill {\int }_S^{}{\mathit{E}}_{\mathbf{n}}\mathit{d}\mathit{S}& \hfill =\hfill & \frac{\Gamma \Delta T{S}_{\mathrm{c}}}{{\epsilon }_0}\hfill \\ \hfill E_{\mathrm{n}}{S}_{\mathrm{e}}& \hfill =\hfill & \frac{\Gamma \Delta T{S}_{\mathrm{c}}}{{\epsilon }_0}\hfill \\ \hfill E_{\mathrm{n}}& \hfill =\hfill & \frac{S_{\mathrm{c}}}{S_{\mathrm{e}}}\frac{\Gamma \Delta T}{{\epsilon }_0}.\hfill \end{array}$$(1)

Here, $S_{\mathrm{e}}$ and $S_{\mathrm{c}}$ denote the surface areas of the electrode and the crystal’s upper face, respectively; $E_{\mathrm{n}}$ represents the normal component of the electric field at the crystal surface, and ${\epsilon }_0$ is the vacuum permittivity. Because the electrode surface area varies depending on the electrode geometry, a correction must be applied to the evaluation of the normal electric field component. As illustrated in Fig. 1, the Gaussian surface in this study is constructed so as to enclose only the upper side of the electrode surface. The value obtained under this configuration is adopted as the normal component of the electric field at the electrode surface. In addition, the normal electric field in the z-direction is assumed to become approximately uniform once a certain distance away from the electrode surface. This assumption was introduced because, in the present work, we position the study as a fundamental investigation aimed at capturing the overall characteristics of the potential distribution and particle behavior. The validity of this treatment will be examined in greater detail as the model is further developed and refined in future work.

Because the electric field is obtained through an approximate evaluation, this value is subsequently used to compute the maximum potential, which is required for the normalization of later parameters. In this study, the normal electric field $E_{\mathrm{n}}$ is further examined from the standpoint of nondimensionalization, and its value is converted into a dimensionless form through unit unification. Since the electric field has the dimension of volts per meter [V m⁻¹], the normalization is performed following the procedure described below. Notably, the normalization formula inherently contains the maximum potential, allowing its value to be determined by inverting the expression:

  

$$\overline{E_z}=\frac{E_{\mathrm{n}}}{{\varphi }_{\text{max}}/z_{\mathrm{d}}}.$$(2)

In Eq. (2), ${\varphi }_{\text{max}}$ denotes the maximum electric potential, and $z_{\mathrm{d}}$ represents the height of the device. Furthermore, the normal component of the electric field is assumed to be aligned with the z-direction of the apparatus. Based on this assumption, the normalization formula for the electric field is defined as follows,

  

$$\overline{E_z}=-\frac{\mathit{\partial }\overline{\varphi }}{\mathit{\partial }\overline{z}}.$$(3)

Using Eqs. (2) and (3) to calculate the maximum potential yields the following expression,

  

$${\varphi }_{\text{max}}=-\frac{S_{\mathrm{c}}}{S_{\mathrm{e}}}\frac{\Gamma \Delta T{z}_{\mathrm{d}}}{{\epsilon }_0}\frac{1}{\frac{\mathit{\partial }\overline{\varphi }}{\mathit{\partial }\overline{z}}}.$$(4)

This value is then used in all subsequent normalization procedures.

To enable consistent computational analysis, all physical parameters in this study were normalized. Normalization transforms physical quantities into dimensionless variables, typically scaled from 0 to 1, thereby eliminating unit dependency and enabling uniform comparison across simulations conducted under varying conditions. The example for normalization is as follows:

  

$$\overline{\varphi }=\frac{\varphi }{{\varphi }_{\text{max}}}\left(\begin{array}{ccccc}\hfill 0\hfill & \hfill \le \hfill & \hfill \varphi \hfill & \hfill \le \hfill & \hfill {\varphi }_{max}\hfill \\ \hfill 0\hfill & \hfill \le \hfill & \hfill \overline{\varphi }\hfill & \hfill \le \hfill & \hfill 1\hfill \end{array}\right).$$(5)

The key parameter required in Eq. (5) is the maximum potential, ${\varphi }_{\text{max}}$. In this study, the electric potential $\varphi$, defined over the range 0 to ${\varphi }_{\text{max}}$, is normalized to the dimensionless quantity $\overline{\varphi }$ within the interval 0 to 1 and subsequently used in later computations. To obtain the potential distribution, Laplace’s equation is employed. Under the assumption that the potential is time-independent and that charge is present at a specified location, the electric potential is defined through the following equation:

  

$$\Delta \overline{\varphi }=\frac{1}{\overline{r}}\frac{\mathit{\partial }}{\mathit{\partial }\overline{r}}\left(\overline{r}\frac{\mathit{\partial }\overline{\varphi }}{\mathit{\partial }\overline{r}}\right)+\frac{{\mathit{\partial }}^2\overline{\varphi }}{\mathit{\partial }{\overline{z}}^2}=0.$$(6)

The normalized coordinates $\bar{r}$ and $\bar{z}$ are defined by scaling the radial and axial coordinates with respect to the device radius $r_{\mathrm{d}}$ and height $z_{\mathrm{d}}$, respectively. In the numerical implementation, the mesh resolution in the r–z plane is specified by the number of mesh, which is set to 256 in the present study. The detailed computational configuration and boundary conditions are summarized in Table 1 and illustrated in Fig. 2.

Table 1. List of numerical parameters: physical constants, parameters, and variables.

Constants	Symbol	Value
Proton mass (kg)	$m_{\mathrm{p}}$	1.67 × 10−27
Electron mass (kg)	$m_{\mathrm{e}}$	9.11 × 10−31
Elementary charge (C)	$e$	1.60 × 10−19
Speed of light (m/s)	$c$	3.00 × 108
Pyroelectric coefficient (C/m2·K)	$\Gamma$	1.90 × 104
Boltzmann constant (J/K)	$k_{\mathrm{B}}$	1.38 × 10−23
Vacuum permittivity (F/m)	${\epsilon }_0$	8.85 × 10−12
Parameters
Device radius (m)	$r_{\mathrm{d}}$	0.0765
Device height (m)	$z_{\mathrm{d}}$	0.460
Aspect ratio of the device (-)	$e_0$	0.166
Crystal Base Area (m2)	$S_{\mathrm{c}}$	1.00 × 10−4 π
11B Target Base Area (m2)	$S_{\mathrm{t}}$	9.00 × 10−4 π
Hydrogen gas temperature (K)	$T_{{\mathrm{H}}_2}$	300
Hydrogen gas pressure (Pa)	$P_{{\mathrm{H}}_2}$	1
Variables
Crystal heating temperature (K)	$\Delta T$	0.01–100
Surface charge on the crystal (C)	$Q$	—
Maximum Electric Potential (V)	${\varphi }_{\text{max}}$	—

Fig. 2.  Schematic diagram of the cylindrical device in the r–z plane, including charge placement coordinates and boundary conditions.

In Fig. 2, the schematic of the device cross-section in the r–z plane is shown along with the assumed boundary conditions and charge locations. The initial mesh is defined on the left-hand side of the domain, which corresponds to a critical region for field formation. The potential is computed in the r–z plane, assuming axial symmetry and a null electric field in the azimuthal direction (θ-direction). The resulting potential distribution is then extended uniformly in the θ-direction over 360°, establishing the three-dimensional electric field inside the device. The red line in the figure indicates the upper surface of the crystal, where positive surface charge accumulates, while the blue lines denote the grounded surfaces—namely, the lower crystal surface and the cylindrical device wall. In the present simulation, the device radius and height were set to 0.0765 and 0.46 m, respectively.

The temperature difference range considered in this study was determined based on practical and physical constraints associated with the pyroelectric crystal and the heating method. The lower limit of $\Delta T$ = 0.01 K was selected in consideration of the achievable temperature controllability of the crystal, as maintaining and resolving temperature variations smaller than this value becomes increasingly difficult in practice. The upper limit of $\Delta T$ = 100 K was chosen due to the difficulty of achieving precise and stable temperature control using a heater-based heating method, as well as to ensure that the crystal temperature remains sufficiently below the Curie temperature of the pyroelectric material (approximately 900 K), beyond which the pyroelectric effect vanishes.

In this study, the pyroelectric crystal is modeled as a cylinder with a radius of 10 mm and a height of 30 mm, depicted in yellow in Fig. 3.

Fig. 3.  Schematic cross-sectional views of the experimental setup in the x–z plane. The left, center, and right panels show configurations with a disk electrode, a disk + needle electrode, and a disk + cylinder electrode, respectively. The yellow region indicates the pyroelectric crystal.

To investigate how electrode configuration affects the proton impact rate on a solid ¹¹B target, three types of electrodes were considered. The first configuration, referred to as the “Disk electrode,” serves as the baseline geometry and represents a flat planar conductor placed on top of the crystal surface. In previous studies, a metallic needle was mounted vertically at the center of this disk to locally intensify the electric field and promote field ionization near the tip; this second configuration is defined here as the “Needle electrode.” Furthermore, to direct the electric field toward the target, a third configuration was introduced: the “Cylinder electrode,” which incorporates a hollow cylindrical conductor surrounding the disk periphery. All three configurations were implemented under otherwise identical boundary and material conditions, and their effects on the resulting electric field structures were comparatively analyzed. The geometric differences among these configurations, particularly in terms of field localization and directionality, are critical in determining the efficiency of field-driven ionization and particle acceleration mechanisms. The target position was kept identical across all three electrode configurations to ensure consistent target conditions. Under this unified setup, the fraction of protons impacting the target was calculated and compared for each configuration.

Returning to Fig. 2, the Laplace’s equation is solved over a discretized mesh structure. To facilitate this, Eq. (6) is reformulated into a form suitable for finite-difference operations. The indices i and j correspond to the mesh points in the radial and axial directions, respectively. The governing equation is discretized using a second-order finite difference method, which enables the numerical computation of the potential distribution.

Building on this discretized potential field, the next stage of the computation involves evaluating charged-particle trajectories. For this purpose, the electric field must be interpolated at particle positions, and the equations of motion must be integrated in time. Because the particles may attain high kinetic energies during acceleration, relativistic effects must be considered in the trajectory analysis. However, an important consideration is the possibility that charged particles may reach velocities approaching the speed of light. In particular, lightweight particles such as electrons can reach relativistic speeds—approximately 94.1% of the speed of light—at an energy level of around 1 MeV. If classical mechanics is applied under such conditions, particles may unrealistically exceed the speed of light, leading to physically incorrect results. To address this issue, relativistic effects were incorporated into the calculations. Specifically, the Lorentz factor γ was introduced, and the relativistically corrected velocity components in the x, y, and z directions were defined as follows.

  

$$\begin{array}{ccc}\hfill \gamma & \hfill =\hfill & \frac{1}{\sqrt{1-\frac{v^2}{c^2}}}\hfill \\ & \hfill =\hfill & \frac{1}{\sqrt{1-\frac{|v_x^2+v_y^2+v_z^2|}{c^2}}}\hfill \\ & \hfill =\hfill & \frac{1}{\sqrt{1-\frac{({\xi }^2+{\eta }^2+{\zeta }^2)}{(c^2+{\xi }^2+{\eta }^2+{\zeta }^2)}}}\hfill \end{array}\begin{array}{ccccc}\hfill \xi & \hfill =\hfill & \hfill \frac{v_x}{\sqrt{1-\frac{v^2}{c^2}}}\hfill & \hfill =\hfill & \gamma v_x\hfill \\ \hfill \eta & \hfill =\hfill & \hfill \frac{v_y}{\sqrt{1-\frac{v^2}{c^2}}}\hfill & \hfill =\hfill & \gamma v_y\hfill \\ \hfill \zeta & \hfill =\hfill & \hfill \frac{v_z}{\sqrt{1-\frac{v^2}{c^2}}}\hfill & \hfill =\hfill & \gamma v_z\hfill & .\end{array}$$(7)

It should be noted that Eq. (7) represents the non-normalized form of the equation of motion. To ensure consistency within the overall computational framework and to enable a unified analysis of particle dynamics, the equation must be normalized. The normalized form of the equation of motion, together with the corresponding normalization scheme, is presented as follows:

  

$$\begin{array}{c}\begin{array}{ccc}\hfill \frac{d\overline{x}}{d\overline{t}}& \hfill =\hfill & \frac{1}{e_0}\frac{\overline{\xi }}{\gamma }\hfill \\ \hfill \frac{d\overline{y}}{d\overline{t}}& \hfill =\hfill & \frac{1}{e_0}\frac{\overline{\eta }}{\gamma }\hfill \\ \hfill \frac{d\overline{z}}{d\overline{t}}& \hfill =\hfill & \frac{\overline{\zeta }}{\gamma }\hfill \end{array}\begin{array}{ccc}\hfill \frac{d\overline{\xi }}{d\overline{t}}& \hfill =\hfill & -\frac{1}{e_0}\frac{\mathit{\partial }\overline{\varphi }}{\mathit{\partial }\overline{r}}\cos \theta \hfill \\ \hfill \frac{d\overline{\eta }}{d\overline{t}}& \hfill =\hfill & -\frac{1}{e_0}\frac{\mathit{\partial }\overline{\varphi }}{\mathit{\partial }\overline{r}}\sin \theta \hfill \\ \hfill \frac{d\overline{\zeta }}{d\overline{t}}& \hfill =\hfill & -\frac{\mathit{\partial }\overline{\varphi }}{\mathit{\partial }\overline{z}}.\hfill \end{array}\\ \begin{array}{c}\left(\begin{array}{c}\hfill \overline{\varphi }=\frac{\varphi }{{\varphi }_{\text{max}}},v_0=\sqrt{\frac{q{\varphi }_{\text{max}}}{m}},\tau =\frac{z_{\mathrm{d}}}{v_0},\overline{t}=\frac{t}{\tau },e_0=\frac{r_{\mathrm{d}}}{z_{\mathrm{d}}}\hfill \\ \hfill \overline{r}=\frac{r}{r_{\mathrm{d}}},\overline{x}=\frac{x}{r_{\mathrm{d}}},\overline{y}=\frac{y}{r_{\mathrm{d}}},\overline{z}=\frac{z}{z_{\mathrm{d}}},\overline{c}=\frac{c}{v_0},\overline{v_x}=\frac{v_x}{v_0}\hfill \\ \hfill \overline{v_y}=\frac{v_y}{v_0},\overline{v_z}=\frac{v_z}{v_0},\overline{\xi }=\frac{\xi }{v_0},\overline{\eta }=\frac{\eta }{v_0},\overline{\zeta }=\frac{\zeta }{v_0}\hfill \end{array}\right)\end{array}\end{array}$$(8)

Using the normalized equation of motion given in Eq. (8), the trajectory of the charged particle was numerically calculated by applying the fourth-order Runge-Kutta method. The azimuthal angle θ was calculated from the x and y coordinates using the arctangent function. In this study, the origin of the three-dimensional Cartesian coordinate system is defined at the center of the bottom surface of the device.

In the trajectory calculation of charged particles, the initial velocity and initial position are critical parameters. The initial velocity was assumed to follow the Maxwell-Boltzmann distribution. Since this is a probability density function, the most probable speed is zero. To randomly sample the velocity components from this distribution, the Box-Muller method was employed.

The initial positions of electrons and protons were determined using different approaches. For electrons, hydrogen atoms were assumed to be uniformly distributed within the gas-filled region, leading to a uniform spatial distribution of electrons throughout the device volume. An equal number of electrons were initially placed along the x-, y-, and z-directions, and particles whose positions overlapped with the outer region of the cylindrical chamber or internal structures such as electrodes and the crystal were excluded. In this study, 31 electrons were uniformly placed along each spatial direction, after which those violating the spatial constraints were removed.

The initial positions of protons were determined based on ionization events induced by electron-hydrogen collisions. The electron-impact ionization cross-section of hydrogen depends on the kinetic energy of incident electrons and reaches a maximum at approximately 70–80 eV, while decreasing rapidly at higher energies. In the present model, the ionization probability was evaluated using the energy of each electron, and protons were generated probabilistically according to the corresponding cross-section values. The energy dependence of the ionization cross-section is shown in Fig. 4, and the associated reaction process is given below.

  

$$\mathrm{H}+{\mathrm{e}}^{\mathrm{-}}\to \mathrm{p}+2{\mathrm{e}}^{\mathrm{-}}.$$(9)

Fig. 4.  Electron-impact ionization cross-section of hydrogen atoms.

It is important to note that the electron kinetic energy relevant to ionization is defined in terms of the relative motion between the electron and the hydrogen atom. Therefore, the kinetic energy should be evaluated based on their relative velocity, rather than the absolute velocity of the electron alone. At high electron energies, the non-relativistic approximation may become invalid. Therefore, relativistic corrections were introduced in this study. Specifically, the relativistic kinetic energy was calculated based on the relative velocity between the electron and the hydrogen atom, and this energy value was employed for ionization assessment. Using the aforementioned procedures, the initial positions of electrons and protons were determined. A schematic summary of these initial conditions is presented in the following Fig. 5.

Fig. 5.  Schematic illustration of the initialization procedure for charged particle parameters.

3.  Results & Discussion

A color contour of the calculated potential distribution is presented in Fig. 6. The black lines in Fig. 6 represent the crystal and the electrode structures. In the potential distribution map, the red regions indicate high-potential areas, while the blue regions correspond to low-potential zones. In all three cases, the maximum potential is approximately ${\varphi }_{\text{max}}$ ∼ 238 $\Delta T$ [kV], and the potential distribution in Fig. 6 is normalized by ${\varphi }_{\text{max}}$. The electric field is concentrated near the electrode in all cases, with the cylinder electrode exhibiting a high-potential region extending into the cylinder. The high-potential regions are primarily concentrated near the electrodes, whereas the majority of the device interior exhibits a significantly low potential.

Fig. 6.  Color contour plots of the electric potential in the x–z plane. The arrangement of the panels is the same as in Fig. 3.

The electrostatic field shown in Fig. 6 accelerates electrons within the simulation domain. At the initial time $t$ = 0, electrons are assumed to be uniformly distributed throughout the computational region. As the simulation proceeds, the accelerated electrons collide with hydrogen atoms in the background gas, leading to ionization and subsequent proton generation. The electron-impact ionization of hydrogen is modeled probabilistically using pseudorandom numbers based on the ionization cross section, following the method described in [9].

Next, the ionization positions of hydrogen atoms due to electron impact—that is, the initial generation locations of protons—are presented.

A comparison of Figs. 7–9 reveals that the overall behavior is generally consistent across the cases. A clear temperature dependence is observed, indicating that the ionization region shifts progressively farther away from the electrode as the heating temperature increases. This trend can be attributed to the fact that a higher crystal temperature leads to an increased potential near the electrode, which enhances the potential difference over a given distance and consequently strengthens the electric field. The resulting stronger field accelerates electrons more rapidly, driving them into a higher-energy regime in which the ionization cross-section of hydrogen becomes significantly smaller, thereby suppressing ionization near the electrode.

Fig. 7.  Proton generation locations near the disk electrode. The left, center, and right panels correspond to crystal heating temperatures of ΔT = 0.01, 5, and 100 K.

Fig. 8.  Proton generation locations near the needle electrode. The left, center, and right panels correspond to crystal heating temperatures of ΔT = 0.01, 5, and 100 K.

Fig. 9.  Proton generation locations near the cylinder electrode. The left, center, and right panels correspond to crystal heating temperatures of ΔT = 0.01, 5, and 100 K.

Thus far, the calculated initial proton generation coordinates have been presented, and their spatial characteristics have been examined. Next, to illustrate how the proton impact rate varies with electrode geometry under different crystal temperatures, Fig. 10 summarizes the impact rates obtained for each electrode and temperature condition.

Fig. 10.  Dependence of proton impact rate on the target as a function of crystal heating temperature for the three electrode configurations.

The proton impact rate is defined as the fraction of generated protons that reach and collide with the solid boron target. For both the disk and needle electrodes, the highest impact rate of approximately 37% is observed at a heating temperature of 5 K, whereas the cylindrical electrode exhibits its maximum impact rate at 1 K. In contrast, when the crystal heating temperature approaches the lower or upper bounds of the temperature range considered in this study (e.g., $\Delta T$ = 0.01 or 100 K), the proton impact rate decreases markedly for all electrode configurations. As indicated in Figs. 7–9, this behavior is associated with temperature-dependent variations in the ionization position. At the lower end of the temperature range, the pyroelectrically induced electric field is relatively weak, causing protons to disperse and resulting in reduced directionality toward the target. Conversely, near the upper end of the temperature range, strong electric fields lead to spatially broadened ionization regions, which similarly degrade into proton directionality. Consequently, temperature conditions near the bounds of the investigated range are unfavorable for efficient proton irradiation of the solid ¹¹B target. These results indicate that selecting an appropriate heating temperature within the considered range is essential for enhancing proton directionality and improving irradiation efficiency.

4.  Comparison with Previous Study

In this section, the effects of electrode geometry and crystal heating conditions on the proton impact rate are discussed based on the present simulation results, with reference to previous study.

This prior study represents a significant experimental investigation focused on D-D pyroelectric fusion. In the experiment, the researchers conducted comparative tests by altering the geometry of the pyroelectric crystal itself, as shown in Fig. 11. A deuterated material was used as the target, and the interior of the chamber was filled with deuterium gas. The crystal temperature was set to increase by 100 K over a period of five minutes. A.M. Kovanen et al., concluded that neutron emission was observed, indicating that D-D fusion reactions had occurred [10]. However, the reaction yield was extremely low, rendering the system unsuitable for energy extraction. Instead, they suggested potential applications as a neutron source.

Fig. 11.  Schematic diagram of the experimental setup in a previous study [10].

Nevertheless, when these findings are compared with the results of our current simulations, several questions and inconsistencies emerge. In the following sections, we examine these issues in greater detail and provide a critical discussion based on our analysis.

As shown in Fig. 11, the pyroelectric crystal used in the previous experimental study had a radius and thickness of 20 mm. The copper disk electrode had a radius that was 2 mm smaller than that of the crystal. However, detailed parameters such as the overall dimensions of the apparatus, the length of the needle electrode, and the temperature variation of the crystal with the target attached were not explicitly specified in the original paper. Therefore, in the present discussion, the target potential is assumed to be 0 V. The following analysis is based on assumptions derived from the schematic description of the experimental setup.

In the previous study, a needle electrode was employed on the pyroelectric crystal, and the crystal heating temperature was set to 100 K. When compared with the results of the present numerical simulations, these experimental conditions appear to be unfavorable for achieving efficient positive-ion irradiation of a solid target. In this paper, three potential issues that may have limited the performance of the previous experiment are discussed.

The first issue concerns needle electrode geometry, which may have prevented most of the positive ions generated from reaching the target. As shown in Fig. 5, charged particles are accelerated perpendicular to the equipotential lines; consequently, only ions located directly above the needle tip along the central axis can propagate toward the target, while the remainder are deflected outward. In the present simulations, a needle radius of 1.5 mm is assumed, whereas the previous experiment employed an extremely sharp needle with a radius of approximately 70 nm. Under such conditions, it is highly probable that most positive ions generated in the experiment were deflected away from the target.

The second issue relates to the location of proton generation due to ionization. In the previous experiment, the ionization region may have been positioned too far from the electrode, such that the majority of generated ions did not experience sufficient acceleration by the pyroelectric electric field. As indicated in Fig. 8, the ionization region for the needle electrode configuration is located at positions considerably distant from the pyroelectric crystal, where the influence of the crystal-derived electric field is expected to be relatively weak. Consequently, under the needle electrode configuration with a crystal heating temperature of 100 K, the directionality of positive ions toward the target may have been degraded, resulting in a reduced impact rate.

The third potential issue arises from an unfavorable combination of the ionization position, system scale, and crystal heating temperature. As shown in Fig. 8, the ionization region in the needle electrode configuration is located approximately 0.1–0.15 m away from the electrode, whereas the crystal-to-target distance in the experiment was only about 0.07 m. This suggests that electrons may have reached the electrode surface before inducing sufficient ionization. Electron-impact ionization of hydrogen is most efficient at electron energies of approximately 50–70 eV; however, under high heating conditions, a steep potential gradient forms near the electrode, causing electrons to be rapidly accelerated over a short distance. As a result, electrons traverse the energy range favorable for ionization with minimal residence time and are quickly transferred to higher-energy regimes, leading to a reduced probability of ionization. This tendency is not expected to differ significantly for deuterium.

Considering the three limitations discussed above—namely, positive-ion deflection, inappropriate ionization locations, and the potential failure of sufficient ion generation—it is expected that improvements in the experimental conditions could increase the number of positive ions impinging on the target and consequently enhance the beam output. If the acceleration conditions for positive ions can be appropriately optimized, pyroelectric fusion may offer a viable pathway toward compact and low-cost fusion experimental devices. For instance, if positive ions can be reliably generated and positioned directly above the pyroelectric crystal, they can be efficiently accelerated, leading to a substantial increase in the number of positive ions colliding with the solid target. These considerations indicate that significant room for improvement remains in the current approach, and future studies should systematically explore such possibilities.

5.  Conclusion

In this study, numerical simulations were carried out for a pyroelectric fusion system with three different electrode geometries in order to investigate proton generation positions, particle trajectories, and impact rates with a solid ¹¹B target. The electric field distributions were evaluated for each configuration, and the maximum potential was found to scale with the crystal temperature difference as ${\varphi }_{\text{max}}$ ≃ 238 $\Delta T$ [kV]. Analysis of the ionization positions clarified the regions where protons begin to be accelerated, indicating that a considerable fraction of generated protons may deviate laterally and fail to reach the target. The proton impact rate with the target was systematically examined as a function of electrode geometry and heating temperature difference ΔT. The results show that the impact rate strongly depends on these parameters, with the highest values obtained under moderate heating conditions rather than at large temperature differences. Excessively large potentials tend to shift the ionization region radially outward, thereby reducing the probability that protons are directed toward the target. This finding suggests that reducing the heating temperature can, counterintuitively, improve the proton impact rate. Overall, the present results demonstrate that the efficiency of pyroelectric fusion devices is highly sensitive to device configuration and operational parameters. By appropriately controlling the ionization region, optimizing the heating conditions, and refining target placement, an increase in the number of positive ions colliding with the target may be achievable. Although the present study is fundamental in nature, it highlights the importance of systematic design optimization and suggests that pyroelectric fusion systems retain significant potential for further technical improvement.

 Acknowledgements

This work was performed with the support and under the auspices of the NIFS Collaboration Research program (NIFS24KISC009 and NIFS24KIPC003).

References

• [1]   J.  Kawai et al., X-Ray Spectrom. 41, 216 (2012).
• [2]   B.  Naranjo et al., Nature 434, 1115 (2005).
• [3]   J.  Geuther et al., Phys. Rev. Lett. 96, 054803 (2006).
• [4]   D.  Gillich et al., Nucl. Instrum. Methods Phys. Res. A 602, 306 (2009).
• [5]   M.M.  Nasseri, Nucl. Instrum. Methods Phys. Res. B 362, 45 (2015).
• [6]   S.  Mohtashami et al., J. Appl. Phys. 135, 200701 (2024).
• [7]   S.H.  Sikora et al., J. Fusion Energy 35, 538 (2016).
• [8]   K.  Hanamoto et al., Nucl. Instrum. Methods Phys. Res. A 669, 66 (2012).
• [9]   T.  Takahashi et al., Phys. Plasmas 11, 3131 (2004).
• [10]   A.M.  Kovanen et al., D-D nuclear fusion using different size pyroelectric crystals, (IEEE, 2009) pp. 1989–1992.