Numerical simulation of the charge balance and temperature evolution in an electron beam ion trap

A computer code has been developed to simulate and study the evolution of ion charge states inside the trap region of an electron beam ion trap. In addition to atomic physics phenomena previously included in similar codes such as electron impact ionization, radiative recombination, and charge exchange, several aspects of the relevant physics such as dielectronic recombination, ionization heating, and ion cloud expansion have been included for the first time in the model. The code was developed using object oriented concepts with database support, making it readable, accurate, and well organized. The simulation results show a good agreement with various experiments, and give useful information for selection of operating conditions and experiment design.


I. INTRODUCTION TO ELECTRON BEAM ION TRAPS
Highly charged ions (HCIs) have various applications in physics [1].Accurate knowledge about HCIs and their interactions with electrons gives us useful information for many kinds of plasmas, including those applicable to astrophysics, nuclear fusion, microelectronics, and medicine.
Cold HCIs are essential for precise spectroscopic measurement and in recent years they have been a topic of considerable interest [2].The electron beam ion trap (EBIT) is considered to be an effective way to produce and study HCIs.It is a compact and relatively inexpensive device.As the name implies, it produces and traps very highly ionized ions using an electron beam, and was designed from the beginning for x-ray spectroscopic measurements of the trapped ions although it has subsequently been used for many other studies involving HCIs [2].
In an EBIT, an electron beam passes through a trapped ion cloud and, due to collisions with the ions, various atomic processes occur.The principle was developed from the electron beam ion source (EBIS) as first built by Donets [3,4].
A high quality electron gun is used to produce a highdensity electron beam which is compressed by the magnetic field of a solenoid.Highly charged ions are created by electron impact ionization in the trap region and can be extracted.The EBIS was conceived as a source of highly charged ions.Much subsequent work has focused on the development either of higher voltage machines to produce higher charge states or of smaller, cheaper sources.
The first electron beam ion trap (EBIT) was developed by Levine et al. [5] as a development of the EBIS.The main difference is that the EBIT had a much shorter trap length compared to a typical EBIS.The EBIT was designed so that the physics inside the trap could be studied in greater detail.A pair of coils was used instead of a continuous solenoid.This allows direct spectroscopic access to the trap region perpendicular to the electron beam.One further possible advantage of a shorter trap is the suppression of certain types of plasma instabilities.Levine et al. [5] concluded that the rotational electron-ion two stream instability was suppressed for short trap lengths.
The original design was used to develop several other EBITs around the world including Oxford [6], the National Institute of Standards Technology (NIST) [7], and Berlin [8].A second EBIT was built at Lawrence Livermore National Laboratory (LLNL) and the original device upgraded for higher energy operation and was called ''Super-EBIT'' [9].The main difference between a normal EBIT and a super-EBIT type of machine is the power supply arrangement.In a super-EBIT type of machine, the gun and collector assemblies are floated at a high negative voltage with respect to the laboratory earth.The Tokyo EBIT [10], the Freiburg EBIT [11] (now operating in Heidelberg), and more recently the Shanghai EBIT [12] were designed using the same principle but with some new features.In recent years, other designs have been developed including smaller, more compact EBIT devices with permanent magnets, such as the Belfast EBIT [13] and the Dresden EBIT [14].
An EBIT consists of three parts: electron gun, trap, and collector.The electron gun is used to generate an electron beam, which is subsequently accelerated to high energy and compressed to have a high density.This beam passes through the trap and finally the collector collects the electrons after deceleration.
Figure 1 shows the relationship between the various concepts underlying the machine physics of EBITs.The high current (typically 10 s or 100 s of mA), high energy (typically 3-200 keV) electron beam is responsible for ion creation through electron impact ionization, one of the charge-changing reactions shown of Fig. 1.The beam acts as a heater to increase the ion temperatures through electron beam heating and ionization heating.In the mean time, it also provides the radial potential valley for ion trapping, which gives rise to the ion spatial distributions.
It is a convenient and reasonable assumption that the trap environment is cylindrically symmetric.The ions move about the trap with a characteristic time scale of 100 ns or less, being subject to an electric field due to the drift tubes and the electron beam and the magnetic field applied to compress the electron beam.These ions undergo ion-ion collisions, typically on the ms time scale, resulting in energy exchange and cross-field diffusion.The energies of ions of a given species and charge state are treated as being characterized by a single temperature, regardless of the degree of freedom [15][16][17][18], an assumption which is justified because the ion-ion collisions occur much quicker than the processes whereby ions change charge state.
Typically on a 10 ms to 10 s time scale, depending on the exact conditions, the ions are involved in various reactions that change their charge state.Predominant among these reactions are electron impact ionization, radiative recombination (both occurring inside the electron beam), and charge exchange due to interaction with neutral gas.Again depending on the exact conditions, an individual ion can be trapped for times as long as several hours.It is this hierarchy of time scales that implicitly underpins most of the simulations of EBIT charge-balance dynamics.The cross-field diffusion acts to allow the ions to form a Boltzmann distribution as described by Eq. ( 18) below.The validity of this assumption has been demonstrated elsewhere [16,17].
In this paper, the physics occurring inside an EBIT is discussed and a new simulation code, numerical simulation of charge balance (NSCB), is described.Some new features are included in the code.In Sec.II, we show the basic physics process in the EBIT.In Sec.III, the new features of the code are described.In Sec.IV, some results are compared with experiments.In Sec.V we outline some extensions of the work while in Sec.VI we summarize our conclusions.

II. THE PREVIOUS STUDY OF PHYSICS INSIDE EBITS
The first evolution model to describe charge-balance dynamics inside EBITs was that of Penetrante et al. [19], with subsequent refinements being made by Margolis [6], Fussmann et al. [16], and Liu et al. [20].The evolution in this model is determined by various atomic processes, which can be described as a charge evolution staircase shown in Fig. 2. In the figure, ''EI'' denotes electron impact ionization, ''RR'' denotes radiative recombination, ''DR'' denotes dielectronic recombination, the resonant counterpart to RR, and ''CX'' denotes charge exchange, whereby an ion collides with a neutral atom, and an electron is transferred from the atom to the ion.
There are two directions in the charge evolution model.In the ionization direction, the low charged ions lose electrons and become highly charged ions.In the recombination direction, the ions capture electrons and become lower charged ions or neutral atoms.Considering also the escape of ions, an evolution equation was established by Penetrante et al. [19]: where i are the reaction rates for electron impact ionization, radiative recombination, and charge exchange with neutral atoms, respectively, for an ion in charge state q i .R ESC i is the escape rate.N i is the ion density.The reaction rate denotes how quickly ions of one charge state change to another charge state due to a particular process.The EI reaction rate can be written as where J e is the electron beam current density, e is the electron charge, fðe; iÞ is the overlap factor with the electron beam, and EI i is the EI cross sections which can be calculated, for example, using the relativistically corrected Lotz formula [21][22][23][24].The RR reaction rate can be written in exactly the same manner, replacing EI by RR throughout.RR cross sections can be calculated by the Kim-Pratt formula [25].Similarly, DR rates can be calculated from DR cross sections although there is not a general formula describing all instances of this process.The charge exchange rate is given by where N 0 is the neutral gas density, CX i is the charge exchange cross section which can be calculated by the Mu ¨ller formula [26], and " i is the ions' average speed (averaged over a Maxwellian distribution): where m is the ion's mass.Multiple electron capture [27,28] is not included in the current version.The overlap factor fðe; iÞ is defined as the ratio of the number of ions of a given charge state in the electron beam to the total number of ions of that charge state.As such, it is a factor that acts to scale both heating by the electron beam and charge-changing reactions that involve the electron beam (RR, DR, and EI).These causal relationships are shown alongside the appropriate arrows in Fig. 1.
In the treatment of Penetrante et al. [19], the overlap factor was determined by the characteristic radius r i , given by r i ¼ r e ffiffiffiffiffiffiffiffiffi kT i q i eV e q , when q i eV e > kT i , and where r i ¼ r e exp½ 1 2 ð kT i q i eV e À 1Þ when q i eV e < kT i .V e is the potential difference between beam center and beam edge.The overlap factor is given by fðe; iÞ ¼ r e r i 2 : (5) The escape rate is obtained from the Fokker-Planck equation as where !i ¼ q i eV=kT i , V is the potential trap depth either axially or radially.v i is the total collision rate, which is given by summing up all the collision rates with other ions, where subscript i, j means ions in charged states q i and q j , respectively.The collision rate v ij can be expressed as where lnÃ ij is the ion-ion Coulomb logarithm for these two charge states [16].The total escape rate is the sum of the axial and radial escape rates.The energy evolution is based on the following equation: where the first term is the Spitzer heating, it comes from collision between electrons and ions, Here m is the electron mass, M is the mass of the ion, and E e is the electron beam energy.The collision rate v ei can be written as where e is the electron velocity and lnÃ i is the electronion Coulomb logarithm for the ith ion.Here we assume that e = i ) 1. " 0 is the permittivity of free space.The second term in right side of ( 9) is to account for the total energy transfer between different charge states: The last term of the right side of ( 9) is the escape cooling, comprising effects of axial and radial escape.It can be derived from the Fokker-Plank equation: Equations ( 1) and ( 9) form the set of coupled nonlinear evolution equations which govern the charge-balance and temperature evolution in an EBIT.
Fussmann et al. [16] made further improvements to this evolution model, the overlap factor and escape rate being modified as outlined below.One of the essential insights was that the magnetic field plays no role in determining the ion distribution.Cross-field diffusion, mediated by ion-ion collisions, happens sufficiently quickly that radial trapping occurs due to the electrostatic field contributions (predominantly the electron beam).
To obtain the overlap factor, the radial potential distribution should be calculated first.It can be obtained from the Poisson equation, which can be written in cylindrical coordinates (the angular term is eliminated because of symmetry): where ðrÞ is the distribution of charge density.
For an electron beam alone (i.e.ignoring the spacecharge of any trapped ions) writing the electron beam radius as r e , assuming the radial electron density is well described by a Gaussian distribution, the radial density distribution can be written as where Q e is the linear charge density: I e is the beam current, c is the speed of light, is the relativistic velocity ratio, the relativistic parameter is determined by the electron beam energy: Considering the cylindrical symmetry, the electron beam generates a potential well in the radial direction, which can be used to trap the positive ions.As the total charge of the trapped ions increases, the radial potential well becomes shallower due to the trapped ions neutralizing the space charge of the electron beam.The shape of this potential and the distribution of ions have an effect on each other.It is much more difficult to obtain a self-consistent potential when the trap is loaded with ions.
Assuming diffusion across the magnetic field lines happens sufficiently rapidly [16], the ions obey the Boltzmann distribution, for the charge states q i , the charge density distribution is where N i ð0Þ and kT i are the axial density and temperature for the ions in state q i , respectively.Adding the positive charge term to the right-hand side of the (14) gives The solution of Eq. ( 19) can be calculated numerically; Fig. 3 shows the numerical result of the potential and density distributions in the trap.Line 1 is the potential distribution without positive ions, while line 2 is the potential distribution with ions; lines 3 and 4 show the ions' density distributions subject to potential distributions 1 and 2, respectively.The calculation was performed by the following conditions: electron beam energy, 2.65 keV; electron beam radius, 35 m; electron beam current, 68 mA.The positive charged ions here used were Si 14þ at an axial density of 2:5 Â 10 9 cm À3 ; the temperature is 81 eV.
Assuming a square spatial distribution for the electron beam and a Boltzmann distribution for the ions, for a trapped ion in charge state q i , it can be written as where r dt is the radius of the drift tube.
A new term to better express the escape rate has been developed in [16]; it is given by where !i ¼ q i eV=kT i .Either axial or radial escape rate can be calculated by changing V to be either the axial trap depth or radial trap depth, respectively.The total escape rate is the sum of the axial and radial escape rates.

III. INTRODUCTION TO THE NSCB CODE
A. The overview of the model A new simulation code, numerical simulation of charge balance in an EBIT (NSCB), has been developed.The evolution model is based on those already elucidated by Penetrante et al. [18] and Fussmann, Biedermann, and Radtke [16], with many new features added.The code was organized using object oriented methods and with a database to store some physics parameters so they can be conveniently changed by the user.
In the NSCB model, the number density evolution equation is given by Compared with (1), the new term dielectronic recombination (DR) rate R DR i is present in the equations.This is an important process in many kinds of plasma.Under certain conditions it has a large effect on the charge-balance evolution in an EBIT.DR into very highly charged ions has been observed in high temperature fusion and astrophysical plasmas.It is an energy dependent process, which can only occur when the free electron energy added to the binding energy of the level in which it is captured is equal to the excitation energy of the bound electron.Experimentally determined resonance strengths [29] were stored in a database.
The energy evolution equation is given by The third term of the right side accounts for ionization heating [18], which is included for the first time in the NSCB code.
The expansion rate and the magnetic confinement are also considered in the model.These new features are described in the following paragraphs.

B. Heating capacity evolution
The heat capacity relates the energy input to the resultant increase in the ion temperatures, or the energy released when their temperatures reduce.It is defined in the usual way as the energy required per unit temperature rise.In this paper, we use the electron volt for both energy and temperature; hence the resultant heat capacity is dimensionless.
In [30], a sudden increase in the yield of x rays was observed as the trap depth was lowered in a controlled manner.The sudden transition underlying this phenomenon can be explained through a calculation of the evolution of the heat capacity as has been performed with the NSCB code.If a Boltzmann distribution is assumed for ions with charge q and temperature kT, the average particle potential for the charge state q is It is useful to define a quantity ¼ qV 0 =kT, where V 0 is the difference between the center of the beam and its edge and the critical temperature as kT c ¼ qV 0 .At this temperature in the uncompensated limit (i.e.where space charge due to the trapped ions is neglected), ¼ 1, the average potential will become infinite; this is the mathematical signature of a phase transition.In practice this infinity can be avoided because the integrals of Eq. ( 24) have an upper limit of the drift tube radius.Nevertheless, the heat capacity varies as a function of temperature, going through a maximum at the critical temperature as is shown in Fig. 4 below.It is possible that there are two different modes of EBIT operation, one where the ions are cold and predominantly escape axially (i.e.> 1) and a second where the ions are dynamically streaming out predominantly escaping radially (i.e.< 1).Through inclusion of the heat capacity in our model, we connect these two possible behaviors for the first time in a single simulation code.
When the ionic space charge is included, in the asymptotic limit the average potential takes the same form with V 0 replaced by V 0 ð1 À Þ; here ¼ P qN q =N e is the ratio of positive charge to negative charge.The critical temperature is given by kT c ¼ qV 0 ð1 À Þ.The (dimensionless) heat capacity is therefore calculated by where NUMERICAL SIMULATION OF THE . . .Phys.Rev. ST Accel.Beams 12, 014401 (2009) 014401-5 Figure 4 shows the relationship between C V and the scaled temperature kT=kT c .The calculation was done for argon ions, Ar 18þ .
In the high temperature limit, the heat capacity has the value of 3=2.This is because the free particle has 3 degrees of freedom.If the temperature is low, the trapping potential has the form of a harmonic oscillator in two directions perpendicular to the electron beam, which adds two new degrees of freedom to particles; the heat capacity has the value of 5=2.Near kT=kT c ¼ 1, the heat capacity is very large because a large amount of energy is required to drive the system through the phase transition.In the model, the heat capacities are calculated numerically.The numerical values are then used as parameters affecting the evolution of ion cloud expansion.Although this phenomenon has been discussed previously [31,32], the NSCB code is the first code to incorporate them into a full description of EBIT charge dynamics.

C. Ion cloud expansion
In the number density evolution equations (22), N i is unambiguously defined as the number density on the axis (i.e.r ¼ 0).As its temperature increases, an ion cloud will expand, which causes a reduction of the axial number density so that the total number of trapped ions of a given species is conserved when the temperature changes, ignoring the effect of all other processes.
For computational efficiency and stability, an effective volume is used here to calculate the expansion revision while the full ion distribution is retained for all other purposes.Defining the effective volume as where l dt is the axis trap length, r eff is the effective radius: With the temperature rising, the effective radius becomes larger and results in an increase of the effective volume, which in turn causes a decrease of the axial number density.The revised axial number density can be expressed as where V eff is the effective volume in the previous time step, and V 0 eff is the present effective volume.N is the number density before expansion revision.Implicit in our approach for revising the axial number density as described by Eq. ( 29) is the use of a uniform density distribution since we use an effective volume.However, for this computational step we found results when working with the full distribution of Eq. ( 18) were essentially the same although they were computationally much more expensive and less stable.
The expansion of the volume also causes a decrease of the temperature as work must be done in expanding the ion cloud.Consider that the total number of ions is conserved.The process is an adiabatic process.According to the characteristic equation for the adiabatic process, where V is the effective volume from Eq. ( 24).In the trap, the heat capacity C V for an ion cloud is a function of temperature.According to thermodynamics [33], can be determined by the ratio of the constant pressure heat capacity C P and the constant volume heat capacity C V .As discussed above, both of the heat capacities are dimensionless in our model.Hence, so where V eff is the previous effective volume, and V 0 eff is the present effective volume.T is the temperature without expansion revision.
The heat capacity C V can be calculated numerically as discussed in Sec.III B. The NSCB code contains the first reported inclusion of the effect of ion cloud expansion on the temperature in such simulations.In the NSCB model, the formula ( 21) is used to calculate the axis escape, but for the radial escape, the effect of magnetic confinement has been included.Extending the discussion of Fussmann et al. [16] about the axial escape rate, the radial escape rate is given by the same formula as that for axial escape rate but the magnetic confinement is added to ! i to form a new term !rad i .The differential of the effective radial trapping voltage of the ions is given by where v is the velocity of the ions.Because the magnetic field is axial, B ¼ Be z , the item in brackets can be expressed as v Â B ¼ B e r þ B r e , and dr ¼ dre r , so Assuming a Maxwellian temperature distribution and supposing equipartition, " z ¼ " r ¼ " , the average veloc- q .If the velocity of some ions is large enough, they will escape.An ion must overcome the confinement energy to escape, which requires that the kinetic energy obeys the inequality where where V is the voltage difference between beam center and trap edge, i.e., the inner radius of the central drift tube.Hence, Now, assuming a 1-dimensional Boltzmann distribution for the ions' number density, for radial escape, the radial escape rate is Equation (38) is the new radial escape rate term, which has been derived for the first time in this work.In the new program this radial escape term is used with the factor ! rad i determined from (37).

E. Ionization heating
Ionization heating [18,34] occurs when ionizing reactions occur at different electrostatic potentials within the region where the ions are trapped.For example, consider two ions each of charge number q i , each with the same total energy (sum of kinetic and potential energy), undergoing ionization at different radii and hence at different values of electrostatic potential.Before the ionization events each one has the same energy However, after the ionization the total energies become Since the potentials are not equal, these two energies are also not equal, differing by an amount ÁE ¼ e½ðr 1 Þ À ðr 2 Þ.Subsequent motion resulting in collisions will redistribute the energy, resulting in an increase in temperature.The phenomenon can happen anywhere inside the electron beam and is not limited to the simple case described above.Supposing the ions have a Boltzmann distribution throughout a top-hat electron beam and they are situated in the trap potential shape given in Sec.II, then the average change in energy per ionization event is given by Z r e 0 ðrÞ exp À q i ðrÞ kT i 2rdr: To find the correct evolution rate, consider the following: It is the energy increase of the ions of charge state i due to ionization of charge state i À 1 which contributes to ionization heating.Analogously, recombination of charge state i þ 1 produces an energy increase of the ions of charge state i.The rate of the total energy increasing is given by Using (2), this expression can be written as Formula ( 41) is used to calculate the ionization heating.It will affect the evolution of ion temperatures.The NSCB code contains the first reported inclusion of ionization/ recombination heating in such simulations.

IV. THE RESULTS AND DISCUSSION
An implicit integral method [35] is used to solve the set of evolution equations ( 22) and (23).It is difficult to calculate a charge-balance evolution result which exactly coincides with the experiment without adjustment of some NUMERICAL SIMULATION OF THE . . .Phys.Rev. ST Accel.Beams 12, 014401 (2009) of the input parameters.There are two main reasons for this.First, not all parameters are exactly known; for example, the neutral gas density inside the trap is not known precisely.Second, the quality of available cross sections also affects the calculation.Although we can store some reliable values to the database interrogated by the code, in general for the processes of interest reliable data is still very scarce, and some scaling laws play an important role.Indeed one of the motivations for developing this code is to determine which cross sections EBIT behavior is sensitive to in various circumstances and hence develop new means to evaluate those cross sections.Previous comparisons with EBIT charge-balance simulations have made comparison to the equilibrium charge balance using extracted ions [19,20].It is not possible to infer absolute populations or densities from experiments that infer the equilibrium charge balance through extracted ion measurements due to uncertainties in the extraction efficiency.This problem is further exacerbated for measurements of the charge-balance dynamics since the temperature of the ions and hence the extraction efficiency that charge state changes with time.We prefer to make a more exact comparison with the explicit charge-balance dynamics derived through selected photonic measurements.Because there are cases where one has available reliable recombination cross sections and the detection efficiency can be well characterized it is possible to infer the absolute number of ions of a given charge state inside the electron beam.Assuming a characteristic electron beam radius this absolute number can be transformed into an absolute density.If the x-ray data is collected so that both the x-ray energy and the time between the x-ray detection and the last ion injection are recorded for each x ray detected, then these absolute measurements can be made including the explicit time dependence of the chargebalance dynamics.Clearly this approach provides a much more exacting test of the model and hence it is the one we have chosen to use.
Simulations have been performed concerning the time dependence of the evolution of hydrogenlike and bare molybdenum ions at fixed beam energy and concerning the time dependence of the evolution of lithiumlike to bare titanium ions while the beam energy was being swept across various dielectronic recombination resonances.In the following paragraphs, the results of these simulations will be discussed along with comparison between them and the corresponding measurements.

A. Molybdenum fixed energy simulation
Electron impact ionization of H-like molybdenum ions was measured by Watanabe et al. [36] using the Tokyo EBIT.The technique they used is an extension of the method used by Marrs et al. [37].At equilibrium, the charge balance of H-like ions and bare ions is determined by the EI of H-like ions and RR of bare ions.
To keep the trapped ions longer, they are cooled by evaporative cooling, whereby a low Z inert gas, such as neon, is injected into the trap region.Once ionized, the low Z ions collide with the trapped ions, taking heat from them and then escaping from the trap.As a result, the temperature of the trapped ions reduces and they stay in the trap longer.
During the experiment, the Mo ions were injected into the EBIT using a metal vapor vacuum arc (MEVVA) source, while neural neon was introduced from an injection line, which is designed to deliver gases directly into the trap region for evaporative cooling.The MEVVA injects a short burst of predominantly singly charged ions just after it fires, whereas the gas injector continuously admits neutral gas.Appropriate source terms were used in the simulations to mimic this behavior.
X rays were detected by the Ge detector and the time between the time MEVVA injection and the x-ray detection was recorded for each x ray along with the pulse height of the detection event from which the x-ray energy could be derived.This was achieved using a multiparameter data acquisition system.Since the RR peaks due to n ¼ 1 recombination of hydrogenlike and bare ions are well resolved and the recombination cross sections are available [35], this data was processed to give the number of hydrogenlike and bare Mo ions in the electron beam as a function of time since injection.This processing has been done in order to compare with predictions made by the NSCB code.
The simulation of the evolution of Mo in the EBIT has been performed.The cross sections provided by [36] were stored in the database with the other cross sections used being derived from the scaling laws discussed in the Appendix.The machine conditions used correspond to those of the experiment; namely, axial magnetic field, 4.0 T; electron beam energy, 64.4 keV; electron beam current, 170 mA; trap potential, 50 V.The above conditions are known from the operating parameters of the machine.The other conditions, such as neutral gas density, initial Mo þ density, and electron beam radius, however are not known precisely and must be treated as adjustable parameters.
The number of the ions inside the EBIT is estimated by the x-ray photons detected by the detector.Nð90 Þ is the number of photons detected at 90 in a solid angle d, which is proportional to the number of ions that have reacted in a volume V [15,38]: where p is the detection efficiency, which is estimated to be $0:5-0:6, n i and n e are the ion and electron densities, respectively, " v e is the average electron velocity of the electron beam, and is the observation time.In the Tokyo EBIT with the detector concerned, the solid angle d is approximately 4:0 Â 10 À3 , Conv ð90 Þ corresponds to the 90 observation convoluted cross sections [36].
X. LU AND F. J. CURRELL Phys.Rev. ST Accel.Beams 12, 014401 (2009) 014401-8 The initial Mo þ density was inferred by the H-like and bare ion's density from [36]; it was estimated to be 4:4 Â 10 8 cm À3 .The neutral neon density was estimated by the pressure of the final stage of the injector and its geometry.In the simulation, it was assumed to be 5:0 Â 10 9 cm À3 although the simulations were not sensitive to the adjustment of this parameter over a wide range.
The evolution is affected by the beam radius since this determines the electron density and hence all rates for processes involving interactions with the electron beam.Simulations were performed with all other parameters being held constant as described above, using the following values for the beam radius: 48, 52, 58, 62, and 66 m.As shown in Figs. 5 and 6, for the experimental result, Hlike ions reach the equilibrium within 3 s after trapping has started, while bare ions take about 7 s.The simulations show the evolution of Mo ions in the EBIT.Seen from Figs. 5 and 6, with different beam radii, the evolution rates are different.The smaller beam radius results in high electron density and quicker evolution rates as expected.This trend provides a way to estimate the beam radius.Comparison with the measured onset of bare ions suggests the beam radius is 58 AE 4 m.
The beam radius can also be estimated in a simpler but more approximate manner from the ion evolution.For the bare ions, early on in the evolution ionization is the predominant process in the trap.The density of bare ions must increase at this stage according to The right side of ( 43) is the ionization rate.An approximate value for this rate can be deduced from Fig. 4 at a time of 4 s to have a value of about 2:3 Â 10 4 cm À3 s À1 .As can be seen from Fig. 5, the H-like ions have a constant density after about 4 s, with a value of about 7:1 Â 10 4 cm À3 .The overlap factor can safely be assumed to be close to 1 since shallow trapping conditions were used in the experiment and indeed this is in agreement with the findings of our more detailed simulations used to produce Figs. 5 and 6.Using the cross section and beam current provided before, the beam radius is estimated to be about 57 m, close to the value from our simulation.Note however that this is a crude estimate based only on a few of the measured data points and hence will have a large uncertainty associated with it.
For the H-like ions, the simulation results reach equilibrium earlier than experiment.One reason is that the quality of the cross sections used in the calculation might be questionable as these cross sections were derived from scaling laws.The EI cross section for the H-like ions is provided by [34]; for other charge states, they are provided by the Lotz formula.The RR cross sections for bare and Hlike Mo are provided by [34], which is reliable, but for other charge states, they are calculated by the Kim-Pratt FIG. 5.The simulation results for the evolution of bare Mo ions inside the EBIT with machine parameters as described in the main text.The evolutions were calculated using different values for the beam radius as described in the text.The lines (from slowest to fastest onset) correspond to a beam radius of 66, 62, 58, 52, and 48 m, respectively.The symbols ''j'' are the experimental data with the errors being determined statistically.The slow falloff at times greater than 14 seconds is due to Ba ions displacing those of Mo, an effect unaccounted for in the simulations.formula.Usually the results from the Kim-Pratt formula are smaller than the actual cross sections, particularly for highly charged systems.This means the H-like ions accumulate in the trap more quickly than is predicted by the simulations.After a time, some of the ions are ionized to become bare ions and the number of H-like ions reduces.Note that the bottleneck in production of bare ions is the ionization of hydrogenlike ions, not the preceding ionization steps.Hence, the time evolution of bare Mo is much less sensitive to the quality of these cross sections, being dominated by the H-like EI cross section which we took from [36].Our simulations then suggest that EI cross section measurements for few-electron highly charged ions should be made to help develop better scaling laws.Figure 7 shows the simulation results of the temperature evolution of bare Mo ions, taken from the same simulations as were used to produce Figs. 5 and 6.The curves correspond to different electron density in the beam, with the beam radii listed in the figure caption.As expected, the higher electron densities give rise to a greater degree of Spitzer heating and hence result in faster ion temperature rise and higher equilibrium ion temperatures.During the period 0-3 s, the temperature rises quickly, so the collision rates with ions of lower charge state become large.The escape of those particles takes energy away, causing the ion temperature to decrease.After 9 s, the heating and cooling reach equilibrium, passing through more than one maximum of temperature for the higher electron beam densities.

B. Titanium evolution simulation with varying beam energy
The event mode files were reanalyzed from the titanium DR experiment of O'Rourke et al. [39].Again in the event files x-ray events were time stamped with the time since injection and the beam energy for each x-ray event was also recorded.In the experiment, the low charged Ti ions were injected to the trap from a MEVVA.The potential on the trap was set to 4 kV with respect to earth, while the potential on the electron gun was fixed at À3 kV.The electron beam current was 50 mA.After injection, a sawtoothlike waveform with an amplitude of 8 kV and a frequency of 250 Hz was applied to the trap power supply so the beam energy was quickly swept between 11 and 3 kV with a period of 4 ms.This means the beam was sufficiently energetic to create bare ions for some of the time.
The neutral neon was injected continuously as an evaporative coolant species during experiment.After 1 s the trap was emptied and the Ti was reinjected from the MEVVA to avoid heavier impurity ions, such as barium from the cathode, accumulating in the trap region.
The total number density can be determined by summing the densities of bare, H-like, He-like, and Li-like ions after 0.66 s, to be 4:2 Â 10 8 cm À3 [39].
Based on the experiment, a simulation of the evolution of Ti ions in the EBIT has been performed.The machine conditions used were: magnetic field, 4 T; trap potential, 400 V; electron beam current, 50 mA; and electron beam radius, 30 m.An impulsive source term for low charged Ti ions was used to correspond to them being injected from the MEVVA as in the experiment.A constant neutral neon gas source term corresponding to a density of 5:0 Â 10 9 cm À3 was also used to correspond to its injection continuously to cool down the Ti ions throughout the experiment.Also in the simulation the electron beam energy scanned from 11 to 3 keV with a period of 4 ms.The DR resonances were included with resonant strengths being taken from [39].The resonances were treated as having a width of 50 eV, to account for the energy spread of the electron beam.Since resonant strengths are used, provided the beam energy steps are significantly smaller than the width used for the resonances, the final result is insensitive to the value of this width.
The EI cross sections were calculated by the Lotz formula and the RR cross sections were calculated by the Kim-Pratt formula as given in the Appendix.The DR resonant strengths for He-like and H-like ions given in [39] were also included in the calculation.Figure 8 shows both the theoretical and the experimental evolution of titanium in the EBIT.When calculating the theoretical evolution the effect of DR was included.Figure 9 shows the simulation without DR.
In Fig. 8, the theoretical lines fit to the experimental data reasonably well although there is still some disagreement.Both of the theoretical results and experimental data show that Li-like ions density drops after about 0.1 s; the He-like ions density reaches a maximum at 0.18 s and then drops down slightly; the H-like ions density increases slightly before 0.6 s and then keeps flat; the bare ions density keeps increasing slowly.The simulation shows lower predicted densities of Hlike ions than that observed in the experiments.There are several possible reasons.First, extraction of data for the Lilike ions from the KLL peak (using standard inverse Auger notation) is difficult due to it overlapping with other stronger peaks.Second, not all the parameters are known exactly, such as neutral gas density and the cross sections.The calculation is quite dependent on the quality of the cross sections available.At present, they are calculated from the scaling laws so the quality is unsure.For example, the cross section from the Lotz formula for high energy, highly charged ions is smaller than the real cross section [6].
Figure 9 shows the result without DR.Compared with experiment, the theoretical evolution onsets for Li-like and He-like are quicker so the densities of H-like and bare ions increase quickly.The reason is that without DR, the recombination rates are reduced and the effect of the ionization rates is enhanced.The evolution towards H-like and bare ions becomes quicker.Comparing Figs. 8 and 9, it is found that the simulation with DR agrees with the experiment better.It also shows that the DR processes affect the ions' density evolution significantly and should be included in situations where the beam energy spends time coincident with DR resonance energies for abundant charge states within the charge balance.
Figure 10 shows the simulation results of temperature evolution of Ti ions from Li-like to bare.For times much greater than 0.1 s, the lines cannot be clearly distinguished.The expansion of part of this graph, shown in Fig. 11, reveals that each temperature is oscillating with the period of the beam energy ramp.As in the Mo ions case, the ion temperatures rise quickly due to electron beam heating, and then reduce due to the energy exchange with low charged ions.After 0.2 s, the heating and cooling reach equilibrium.In contrast to the Mo case, the Ti ion temperatures show a wavelike behavior with a period the same as the electron beam energy waveform used to scan through the DR resonances.This can be understood as an effect of the energy dependence of electron beam heating as given in Eq. (10).Although this effect has never been observed directly in EBITs, it has important consequences for experiments that seek to determine cross sections or resonance strengths by comparing measurements made at different energies.The variation of ion cloud temperatures as a function of electron energy (as the ramp waveform is applied) leads to a corresponding variation in electron-ion overlap factors.This effect might then result in systematic errors when making such comparisons.Simulations such as those presented here can then provide a useful guide as to the importance of such systematic errors.Besides the beam energy, the code also has the ability to calculate the evolution with linear changing of other parameters such as trap depth or beam current.

V. ONGOING AND FUTURE WORK
The code is undergoing beta testing prior to our releasing it into the public domain.EBITs are typically run in conditions of ultrahigh vacuum when the recombination rate dominates over the charge exchange rate.Hence, we considered it important that the model correctly accounted for recombination, including dielectronic recombination where necessary.This approach is validated by the results shown in Figs. 8 and 9 where inclusion of DR modified the densities of some species by about a factor of 2. In comparison, Liu et al. [19] have considered the importance of multiple charge exchange and compared an equilibrium charge balance with and without inclusion of this effect.As expected, the effect of multiple charge exchange is to modify the charge balance to lower charge and broaden it somewhat.Again densities change by as much as a factor of 2 in some cases.Clearly then a full model of EBIT charge-balance dynamics needs to include both DR and multiple charge exchange.
It is however difficult to make quantitative comparisons with real EBIT measurements of the importance of charge exchange since the exact neutral density and neutral-ion overlap geometry are subject to significant uncertainties.However, recently Watanabe et al. [40] have published results that contain dips that are artifacts due to multiple charge exchange or escape; the theory does not allow us to determine which.Supposing the dips to be due to multiple charge exchange, their depths are a sensitive probe to the role of multiple charge exchange within the EBIT environment.We are currently incorporating multiple charge exchange into our code in the manner described by Liu et al. [20] and once this is complete we will then compare the results of our model to results of the type shown by Watanabe et al.
Potential users of the code can contact its author, Lu at luxj@uestc.edu.cn for a copy and installation instructions.We anticipate creating a web-based version in the near future.

VI. CONCLUSION
The charge-balance simulations are of great merit to study the physics processes inside the EBIT.They give us useful information to understand the atomic reactions in the trap region.Compared with previous models, the NSCB simulation code includes some processes which were omitted in previous work, such as inclusion of DR, ion cloud expansion, and ionization heating.A database is used to store the cross sections and some other important parameters.Besides the scaling laws, the cross sections can be determined using theoretical methods or from experiment and then transferred into the database.Use of these more precise cross sections improves the accuracy of the simulation.During the development of the code, object-oriented programming methods were used, which is a better way to organize and maintain the code.
Simulations are important for planning future machine building and experiment design.The optimum conditions can be decided from the simulations.There are several parameters which are unknown or difficult to measure such as beam radius and initial particle densities.To estimate the value of these uncertain parameters, the simulations can be made under different conditions and the optimum value determined.In this paper, we have shown how to decide the electron beam radius from the time dependence of bare Mo while the hydrogenlike Mo time dependence could be used to critique the cross sections used.In comparison to an experiment involving DR of Ti, the simulations with and without the inclusion of DR were presented.From the comparison with experiment it is clear that DR should be included in such circumstances.nancial support to do this work.We also gratefully acknowledge Gerd Fussmann and his colleagues for sending us a copy of their EBIT modeling code during the development of the NSCB code.

APPENDIX A: THE CROSS SECTIONS USED IN THE MODEL
In the NSCB code, a database is established which stores the ionization potential for each charge state and which also has the capacity to store the EI, RR, and CX cross sections.The cross sections can be obtained either by experiment or by theoretical methods such as R-matrix or distorted wave approximation and then used for the EBIT modeling.If the ideal data cannot be found in the database, a scaling law is applied to calculate the cross sections.The various scaling laws used are reviewed in this Appendix.

The electron impaction ionization (EI) cross section
If other cross sections are not available for EI then the relativistic Lotz formula is used.In Lotz formula [22,23], the ionization cross section of an ion from charge state i to i þ 1 is given: where E e is the incident electron beam energy.This summation is over subshells of the initial ion, & j and I j being the number of electrons in the jth subshell and the subshell binding energy, respectively, with the ionization potentials for every charge state being obtained from the database.
Based on comparison to experimental results, Lotz gives the parameters a ij , b ij , c ij , for specific ions.For other ions, those parameters are just treated as a ij ¼ 4:5 Â 10 À14 cm 2 eV 2 and b ij ¼ c ij ¼ 0 [6].
The relativistic correction factor R, which was first calculated by Gryzinski [24,25], is also included in the code: 3=2 ; (A2) here, " and are the projectile electron energy and target binding energy, respectively.

Radiative recombination
For radiative recombination, if a cross section is not provided in the database then the Kim and Pratt simple scaling law is used since this law approximates the radiative recombination cross sections well [24]: where is the fine-structure constant and e ¼ @=mc is the Compton wavelength of the electrons.For ions of nuclear charge Z and ionic charge q i , ¼ 2Z 2 eff I H E e ; where Z eff ¼ 1 2 ðZ þ q i Þ, and I H is the ionization potential of the hydrogen atom.The effective quantum number ðn v Þ eff may be written as where n v is the principal quantum number of the valence shell, and W nv is the ratio of the number of unoccupied states to the total number of states in the valence shell.

Charge exchange
Charge exchange (CX) occurs between the positive ions and neutral gas.Inside the trap region, neutral low-Z gas is injected to collide with ions for evaporative cooling and the CX rate is large for highly charged ions.In this process, the neutral atom loses an electron, the ion gains an electron, and its charge state is reduced by one.Only single charge exchange is considered here.A single electron transfers from neutral atom to the ions and remains captured.
A simple scaling rule for one electron transfer cross sections was obtained by Mu ¨ller and Salzborn [27]: where E 0 is the ionization potential (in eV) of the neutral atom.

Dielectronic recombination
At present no scaling law is provided and this process is only considered if the corresponding elements of the database are filled.

FIG. 1 .
FIG.1.An overview of the relationships between various processes occurring inside an EBIT[15].

FIG. 2 .
FIG. 2. Diagram indicating the main couplings considered effecting the change of numbers of trapped species (horizontal or curved arrows) and energy (vertical arrows).The definitions for the abbreviations RR, DR, EI, and CX are given in the main text.Each charge state is also energetically coupled to all the others through ion-ion collisions although this coupling has been omitted for clarity.

FIG. 3 .
FIG. 3. Potential and ion distribution vs radius.Line 1: the potential distribution without including the effect of positively charged ions; line 2: the potential distribution including the effect of positively charged ions; line 3: the ion density distribution with the potential distribution 1; line 4: the ion density distribution with the potential distribution 2.

Figures 5
and 6 show the calculated and measured evolutions of the bare and the H-like Mo ions.The vertical coordinate shows the number of ions of a given species inside the electron beam, and the horizontal coordinate shows the time since injection.The various simulations using different values of the beam radius are compared to the experimental data in these figures.

FIG. 6 .
FIG.6.The simulation results for the evolution of H-like Mo ions inside the EBIT with machine parameters as described in the main text.The evolutions are calculated using different values for the beam radius.The lines (from slowest to fastest onset) correspond to a beam radius of 66, 62, 58, 52, and 48 m, respectively.The symbols ''j'' are the experimental data with the errors being determined statistically.

FIG. 7 .
FIG.7.The simulation results for the temperature evolution of bare Mo ions inside the EBIT with machine parameters as described in the main text.The lines from top to bottom correspond to a beam radius of 48, 52, 58, 62, and 66 m, respectively.

FIG. 10 .
FIG.10.The temperature evolution lines of the lithiumlike to bare titanium ions in the Tokyo EBIT with the effect of DR included.Machine parameters are as described in the main text.The lines from top to down correspond to bare, H-like, He-like, and Li-like, respectively.

FIG. 11 .
FIG.11.The simulation results of Ti ions temperature in a small time interval.The curves from top to down correspond to bare, H-like, He-like, and Li-like ions, respectively.The period of these oscillations is that of the waveform used to drive the electron beam energy due to the change in heating [see Eq. (10)] and is not governed by the faster ion-ion collision time scale (typically 0.1-1 ms).