text:
# References

MERGING THE CEM2K AND LAQGSM CODES WITH GEM2 TO DESCRIBE FISSION AND LIGHT-FRAGMENT PRODUCTION

S. G. Mashnik1, K. K. Gudima2, and A. J. Sierk1

1Los Alamos National Laboratory, Los Alamos, NM 87545, USA

2 Institute of Applied Physics, Academy of Science of Moldova, Kishinev, MD-2028, Moldova

Abstract

We present the current status of the improved Cascade-Exciton Model (CEM) code CEM2k and of the Los Alamos version of the Quark-Gluon String Model code LAQGSM. To describe fission and light-fragment (heavier than He4) production, both CEM2k and LAQGSM have been merged with the GEM2 code of Furihata. We present some results on proton- and deuteron-induced spallation, fission, and fragmentation reactions predicted by these extended versions of CEM2k and LAQGSM. We show that merging CEM2k and LAQGSM with GEM2 allows us to describe many fission and fragmentation reactions in addition to the spallation reactions which are already relatively well described. Nevertheless, the standard version of GEM2 does not provide a completely satisfactory description of complex particle spectra, heavy-fragment emission, and spallation yields, and is not yet a reliable tool for applications. We conclude that we may choose to use a model similar to the GEM2 approach in our codes, but it must be significantly extended and further improved. We observe that it is not sufficient to analyze only A and Z distributions of the product yields when evaluating this type of model, as is often done in the literature; instead it is important to study all the separate isotopic yields as well as the spectra of light particles and fragments.

Introduction

During recent years, for a number of applications like Accelerator Transmutation of nuclear Wastes (ATW), Accelerator Production of Tritium (APT), Rare Isotope Accelerator (RIA), Proton Radiography (PRAD), and others projects, we have developed at the Los Alamos National Laboratory an improved version of the Cascade-Exciton Model (CEM), contained in the code CEM2k, to describe nucleon-induced reactions at incident energies up to 5 GeV  and the Los Alamos version of the Quark-Gluon String Model, realized in the high-energy code LAQGSM , able to describe both particle- and nucleus-induced reactions at energies up to about 1 TeV/nucleon.

In our original motivation, different versions of the CEM and LAQGSM codes were developed to reliably describe the yields of spallation products and spectra of secondary particles, without a special emphasis on complex-particle and light-fragment emission or on fission fragments in reactions with heavy targets. In fact, the initial versions of the CEM2k and LAQGSM codes simulate spallation only and do not calculate the process of fission, and do not provide fission fragments and a further possible evaporation of particles from them. Thus, in simulating the compound stage of a reaction, when these codes encounter a fission, they simply tabulate this event (that permits calculation of fission cross sections and fissility) and finish the calculation of this event without a subsequent treatment of fission fragments. To be able to describe nuclide production in the fission region, these codes have to be extended by incorporating a model of high energy fission (e.g., in the transport code MCNPX , where CEM2k and, initially, its precursor, CEM97 , are used, they are supplemented by Atchison’s fission model .

Since many nuclear and astrophysical applications require reliable data also on complex particles (gas production) and light and/or fission fragment production, we addressed these questions by further development of CEM2k and LAQGSM codes. We tried different ways of solving these problems and as a first attempt to describe with our codes both emission of intermediate-mass fragments heavier than 4He and production of heavy fragments from fission, we merged CEM2k and LAQGSM with the Generalized Evaporation Model (GEM) code by Furihata . We have benchmarked our codes on all proton-nucleus and nucleus-nucleus reactions measured recently at GSI (Darmstadt, Germany) and on many other different reactions at lower and higher energies measured earlier at other laboratories. The size of the present paper allows us to present here only a few results, which we choose to be for the GSI measurements on interaction of 208Pb beams with p  and d  targets. Results for other reactions may be found in .

CEM2k and LAQGSM Codes

A detailed description of the initial version of the CEM may be found in Ref. , therefore we outline here only its basic assumptions. The CEM assumes that reactions occur in three stages. The first stage is the IntraNuclear Cascade (INC) in which primary particles can be re-scattered and produce secondary particles several times prior to absorption by or escape from the nucleus. The excited residual nucleus remaining after the cascade determines the particle-hole configuration that is the starting point for the preequilibrium stage of the reaction. The subsequent relaxation of the nuclear excitation is treated in terms of an improved Modified Exciton Model (MEM) of preequilibrium decay followed by the equilibrium evaporative final stage of the reaction. Generally, all three stages contribute to experimentally measured outcomes.

The improved cascade-exciton model in the code CEM2k differs from the older CEM95 version (which is available free from the NEA/OECD, Paris)  by incorporating new approximations for the elementary cross sections used in the cascade, using more precise values for nuclear masses and pairing energies, employing a corrected systematics for the level-density parameters, adjusting the cross sections for pion absorption on quasi-deuteron pairs inside a nucleus, including the Pauli principle in the preequilibrium calculation, and improving the calculation of the fission widths. Implementation of significant refinements and improvements in the algorithms of many subroutines led to a decrease of the computing time by up to a factor of 6 for heavy nuclei, which is very important when performing simulations with transport codes. Essentially, CEM2k  has a longer cascade stage, less preequilibrium emission, and a longer evaporation stage with a higher excitation energy, as compared to its precursors CEM97  and CEM95 . Besides the changes to CEM97 and CEM95 mentioned above, we also made a number of other improvements and refinements, such as: (i) imposing momentum-energy conservation for each simulated event (the Monte Carlo algorithm previously used in CEM provides momentum-energy conservation only statistically, on the average, but not exactly for the cascade stage of each event), (ii) using real binding energies for nucleons at the cascade stage instead of the approximation of a constant separation energy of 7 MeV used in previous versions of the CEM, (iii) using reduced masses of particles in the calculation of their emission widths instead of using the approximation of no recoil used previously, and (iv) a better approximation of the total reaction cross sections. On the whole, this set of improvements leads to a much better description of particle spectra and yields of residual nuclei and a better agreement with available data for a variety of reactions. Details, examples, and further references may be found in Refs. .

The Los Alamos version of the Quark-Gluon String Model (LAQGSM)  is the next generation of the Quark-Gluon String Model (QGSM) by Amelin et al. (see  and references therein) and is intended to describe both particle- and nucleus-induced reactions at energies up to about 1 TeV/nucleon. The core of the QGSM is built on a time-dependent version of the intranuclear cascade model developed at Dubna, often referred in the literature simply as the Dubna intranuclear Cascade Model (DCM) (see  and references therein). The DCM models interactions of fast cascade particles (“participants”) with nucleon spectators of both the target and projectile nuclei and includes interactions of two participants (cascade particles) as well. It uses experimental cross sections (or those calculated by the Quark-Gluon String Model for energies above 4.5 GeV/nucleon) for these elementary interactions to simulate angular and energy distributions of cascade particles, also considering the Pauli exclusion principle. When the cascade stage of a reaction is completed, QGSM uses the coalescence model described in  to “create” high-energy d, t, 3He, and 4He by final state interactions among emitted cascade nucleons, already outside of the colliding nuclei. After calculating the coalescence stage of a reaction, QGSM moves to the description of the last slow stages of the interaction, namely to preequilibrium decay and evaporation, with a possible competition of fission using the standard version of the CEM . But if the residual nuclei have atomic numbers with $`A13`$, QGSM uses the Fermi break-up model to calculate their further disintegration instead of using the preequilibrium and evaporation models. LAQGSM differs from QGSM by replacing the preequilibrium and evaporation parts of QGSM described according to the standard CEM  with the new physics from CEM2k  and has a number of improvements and refinements in the cascade and Fermi break-up models (in the current version of LAQGSM, we use the Fermi break-up model only for $`A12`$). A detailed description of LAQGSM and further references may be found in .

We have benchmarked CEM2k and LAQGSM against most available experimental data and have compared our results with predictions of other current models used by the nuclear community. Figure 1 shows examples of calculated neutron spectra from the interaction of protons with 208Pb at 0.8 and 1.5 GeV compared with experimental data , while Figure 2 gives examples of neutron spectra measured by Nakamura’s group (see  and references therein) from 560 MeV/nucleon Ar beams on C, Cu, and Pb targets compared with our LAQGSM results and predictions by QMD  and HIC  from Iwata et al. . We see that our codes describe well neutron spectra both for proton- and nucleus-nucleus reactions and agree with the data no worse than other models do. Similar results are obtained for other reactions for which we found measured data.

Merging CEM2k and LAQGSM with GEM2

The Generalized Evaporation Model (GEM)  is an extension by Furihata of the Dostrovsky et al.  evaporation model as implemented in LAHET  to include up to 66 types of particles and fragments that can be evaporated from an excited compound nucleus plus a modification of the version of Atchison’s fission model  used in LAHET. Many of the parameters were adjusted for a better description of fission reactions when using it in conjunction with the extended evaporation model. We merged GEM2 (the last update of the GEM code) with CEM2k and LAQGSM as follows: we calculate the cascade and preequilibrium stages of a reaction with our CEM2k or LAQGSM, then we describe the subsequent evaporation of particles and fragments and fission from the remaining excited compound nuclei using GEM2. To understand the role of preequilibrium particle emission, we performed calculations of all the reactions we tested both with emission of preequilibrium particles and without them, i.e., going directly to GEM2 after the intranuclear cascade stage of a reaction described by CEM2k or LAQGSM.

A very detailed description of the GEM, together with a large amount of results obtained for many reactions using GEM coupled either with the Bertini or ISABEL INC models in LAHET are published by Furihata ; many useful details of GEM2 may be found also in our paper . Therefore, we list briefly below only the main features of GEM and forward readers interested in details to Refs. . Furihata did not change in GEM the general algorithms used in LAHET to simulate evaporation and fission. The decay widths of evaporated particles and fragments are estimated using the classical Weisskopf-Ewing statistical theory .

The new ingredient in GEM in comparison with LAHET which considers evaporation of only 6 particles (n, p, d, t, 3He, and 4He) is that Furihata included the possibility of evaporation of up to 66 types of particles and fragments (both in the ground and excited states) and incorporated into GEM several sets of parameters used to calculate inverse cross sections and Coulomb barriers for each ejectile (we use here only default parameters of GEM2). The 66 ejectiles considered by GEM are: n, p, d, t, 3,4,6,8He, 6-9Li, 7,9-12Be, 8,10-13B, 10-16C, 12-17N, 14-20O, 17-21F, 18-24Ne, 21-25Na, and 22-28Mg.

The fission model used in GEM is based on Atchison’s model  as implemented in LAHET , often referred in the literature as the Rutherford Appleton Laboratory (RAL) model, which is where Atchison developed it. There are two choices of parameters for the fission model: one of them is the original parameter set by Atchison  as implemented in LAHET , and the other is a parameter set evaluated by Furihata , used here as a default of GEM2.

The Atchison fission model is designed to only describe fission of nuclei with $`Z70`$. It assumes that fission competes only with neutron emission, i.e., from the widths $`\mathrm{\Gamma }_j`$ of n, p, d, t, 3He, and 4He, the RAL code calculates the probability of evaporation of any particle. When a charged particle is selected to be evaporated, no fission competition is taken into account. When a neutron is selected to be evaporated, the code does not actually simulate its evaporation, instead it considers that fission may compete, and chooses either fission or evaporation of a neutron according to the fission probability $`P_f`$. This quantity is treated by the RAL code differently for the elements above and below $`Z=89`$. The reasons Atchison split the calculation of the fission probability $`P_f`$ are: (1) there is very little experimental information on fission in the region $`Z=85`$ to 88, (2) the marked rise in the fission barrier for nuclei with $`Z^2/A`$ below about 34 (see Fig. 2 in ) together with the disappearance of asymmetric mass splitting, indicates that a change in the character of the fission process occurs. If experimental information were available, a split between regions about $`Z^2/A34`$ would more sensible .

1) $`70Z_j88`$. For fissioning nuclei with $`70Z_j88`$, GEM uses the original Atchison calculation of the neutron emission width $`\mathrm{\Gamma }_n`$ and fission width $`\mathrm{\Gamma }_f`$ to estimate the fission probability as

$$P_f=\frac{\mathrm{\Gamma }_f}{\mathrm{\Gamma }_f+\mathrm{\Gamma }_n}=\frac{1}{1+\mathrm{\Gamma }_n/\mathrm{\Gamma }_f}.$$
(1)
Atchison uses  the Weisskopf and Ewing statistical model  with an energy-independent pre-exponential factor for the level density (see Eq. (6) in ) and Dostrovsky’s  inverse cross section for neutrons and estimates the neutron width $`\mathrm{\Gamma }_n`$ as

$`\mathrm{\Gamma }_n=0.352\left(1.68J_0+1.93A_i^{1/3}J_1+A_i^{2/3}(0.76J_10.05J_0)\right),`$ (2)
where $`J_0`$ and $`J_1`$ are functions of the level density parameter $`a_n`$ and $`s_n\left(=2\sqrt{a_n(EQ_n\delta )}\right)`$ as

$$J_0=\frac{(s_n1)e^{s_n}+1}{2a_n}\text{ and }J_1=\frac{(2s_n^26s_n+6)e^{s_n}+s_n^26}{8a_n^2}.$$

Note that the RAL model uses a fixed value for the level density parameter $`a_n`$, namely

$$a_n=(A_i1)/8,$$
(3)
and this approximation is kept in GEM when calculating the fission probability according to Eq. (1), though it differs from the Gilbert-Cameron-Cook-Ignatyuk (GCCI) parameterization (see Eq. (7) in ) used in GEM to calculate particle evaporation widths. The fission width for nuclei with $`70Z_j88`$ is calculated in the RAL model and in GEM as

$$\mathrm{\Gamma }_f=\frac{(s_f1)e^{s_f}+1}{a_f},$$
(4)
where $`s_f=2\sqrt{a_f(EB_f\delta )}`$ and the level density parameter in the fission mode $`a_f`$ is fitted by Atchison  to describe the measured $`\mathrm{\Gamma }_f/\mathrm{\Gamma }_n`$ as:

$$a_f=a_n\left(1.08926+0.01098(\chi 31.08551)^2\right),$$
(5)
and $`\chi =Z^2/A`$. The fission barriers $`B_f`$ \[MeV\] are estimated as

$$B_f=Q_n+321.216.7\frac{Z_i^2}{A}+0.218\left(\frac{Z_i^2}{A_i}\right)^2.$$
(6)
Note that neither the angular momentum nor the excitation energy of the nucleus are taken into account in the estimate of the fission barriers.

2) $`Z_j89`$. For heavy fissioning nuclei with $`Z_j89`$ GEM follows the RAL model  and does not calculate at all the fission width $`\mathrm{\Gamma }_f`$ and does not use Eq. (1) to estimate the fission probability $`P_f`$. Instead, the following semi-empirical expression obtained by Atchison  by approximating the experimental values of $`\mathrm{\Gamma }_n/\mathrm{\Gamma }_f`$ published by Vandenbosch and Huizenga  is used to calculate the fission probability:

$$\mathrm{log}(\mathrm{\Gamma }_n/\mathrm{\Gamma }_f)=C(Z_i)(A_iA_0(Z_i)),$$
(7)
where $`C(Z)`$ and $`A_0(Z)`$ are constants dependent on the nuclear charge $`Z`$ only. The values of these constants are those used in the current version of LAHET  and are tabulated in Table 5 of Ref.  (note that some adjustments of these values have been done since Atchison’s papers  were published).

The selection of the mass of the fission fragments depends on whether the fission is symmetric or asymmetric. For a pre-fission nucleus with $`Z_i^2/A_i35`$, only symmetric fission is allowed. For $`Z_i^2/A_i>35`$, both symmetric and asymmetric fission are allowed, depending on the excitation energy $`E`$ of the fissioning nucleus. No new parameters were determined for asymmetric fission in GEM.

For nuclei with $`Z_i^2/A_i>35`$, whether the fission is symmetric or not is determined by the asymmetric fission probability $`P_{asy}`$

$$P_{asy}=\frac{4870e^{0.36E}}{1+4870e^{0.36E}}.$$
(8)
For asymmetric fission, the mass of one of the post-fission fragments $`A_1`$ is selected from a Gaussian distribution of mean $`A_f=140`$ and width $`\sigma _M=6.5`$. The mass of the second fragment is $`A_2=A_iA_1`$.

For symmetric fission, $`A_1`$ is selected from the Gaussian distribution of mean $`A_f=A_i/2`$ and two options for the width $`\sigma _M`$ as described in .

The charge distribution of fission fragments is assumed to be a Gaussian distribution of mean $`Z_f`$ and width $`\sigma _Z`$. $`Z_f`$ is expressed as

$$Z_f=\frac{Z_i+Z_1^{}Z_2^{}}{2}\text{ , where }Z_l^{}=\frac{65.5A_l}{131+A_l^{2/3}}\text{ , and }l=1\text{ or 2}.$$

The original Atchison model uses $`\sigma _Z=2.0`$. An investigation by Furihata  suggests that $`\sigma _Z=0.75`$ provides a better agreement with data; therefore $`\sigma _Z=0.75`$ is used in GEM2 and in all our calculations.

The kinetic energy of fission fragments \[MeV\] is determined by a Gaussian distribution with mean $`ϵ_f`$ and width $`\sigma _{ϵ_f}`$. The original parameters in the Atchison model are:

$$ϵ_f=0.133Z_i^2/A_i^{1/3}11.4\text{, and }\sigma _{ϵ_f}=0.084ϵ_f.$$

Furihata’s parameters in GEM2, which we also use, are:

$$ϵ_f=\{\begin{array}{cc}0.131Z_i^2/A_i^{1/3}\text{, for }Z_i^2/A_i^{1/3}900,\hfill & \\ 0.104Z_i^2/A_i^{1/3}+24.3\text{, for }900<Z_i^2/A_i^{1/3}1800,\hfill & \end{array}$$

and

$$\sigma _{ϵ_f}=\{\begin{array}{cc}C_1(Z_i^2/A_i^{1/3}1000)+C_2\text{, for }Z_i^2/A_i^{1/3}>1000,\hfill & \\ C_2\text{, for }Z_i^2/A_i^{1/3}1000,\hfill & \end{array}$$

where $`C_1=5.70\times 10^4`$ and $`C_2=86.5`$. More details may be found in .

We note that Atchison also has modified his original version using recent data and published  an improved (and more complicated) parameterization for many quantities and distributions in his model, but these modifications  are not yet included either in LAHET or in GEM2.

We have merged the GEM2 code with our CEM2k and LAQGSM, initially keeping all the default options in GEM2. We began by concentrating on an analysis of the recent GSI measurements in inverse kinematics as the richest and best data set for testing this kind of model. As mentioned above, to understand the role of preequilibrium particle emission, we performed calculations of all the reactions we tested both taking into account preequilibrium particle emission and ignoring it, i.e., going directly to GEM2 after the intranuclear cascade stage of a reaction described by CEM2k or LAQGSM.

If we merge GEM2 with CEM2k without any modifications, the new code does not describe correctly the fission cross section (and the yields of fission fragments) whether we take into account preequilibrium emission (see the short-dashed line on Fig. 3) or not (see the long-dashed line on Fig. 3). Such results were anticipated, as Atchison fitted the parameters of his RAL fission model when it was coupled with the Bertini INC  which differs from our INC. In addition, he did not model preequilibrium emission. Therefore, the distributions of fissioning nuclei in $`A`$, $`Z`$, and excitation energy $`E^{}`$ simulated by Atchison differ significantly from the distributions we get; as a consequence, all the fission characteristics are also different.

Furihata used GEM2 coupled either with the Bertini INC  or with the ISABEL  INC code, which also differs from our INC, and did not include preequilibrium particle emission. Therefore the real fissioning nuclei simulated by Furihata differ from the ones in our simulations, and the parameters adjusted by Furihata to work the best with her INC should not be the best for us. To get a good description of the fission cross section (and fission-fragment yields) we need to modify at least one parameter in GEM2, namely to adjust the level density parameter $`a_f`$ to get the correct fission cross section (see Eq. (5)), in the case of fissioning nuclei with $`Z88`$ (pre-actinides), and the parameter $`C(Z)`$ (see Eq. (7)) for fissioning nuclei with $`Z>88`$ (actinides). From the dashed lines on Fig. 5 we see that we need to enlarge $`a_f`$ in our code to get a proper fission cross section when we include preequilibrium emission (the excitation energy of our fissioning nuclei and their $`A`$ and $`Z`$ are smaller than provided by the Bertini or ISABEL INC without preequilibrium), and we need to decrease $`a_f`$ in the case without preequilibrium. By increasing $`a_f`$ by 1.06 compared with the original RAL and GEM2 value ($`a_f^{CEM}/a_f^{RAL}=1.06`$), we are able to reproduce correctly with CEM2k+GEM2 the fission cross section for this reaction when we take into account preequilibrium emission (below, we label such results as “with Prec”). In the case with no preequilibrium emission, a proper fission cross section is obtained for $`a_f^{CEM}/a_f^{RAL}=0.947`$ (we label such results as “no Prec”). We choose these values for $`a_f`$ for all our further CEM2k+GEM2 calculations of this reaction and do not change any other parameters.

The solid lines in Fig. 3 show results with these values of $`a_f`$. One can see that the “no Prec” version provide a good description of both the mass and charge distributions and agrees better with the data for these characteristics than the “with Prec” version (that is not true for isotopic distributions of individual elements, as we show below). The “with Prec” version reproduces correctly the position of the maximum in both $`A`$ and $`Z`$ distributions and the yields of fission fragments not too far from these maximums, but the calculated distributions are narrower than the experimental ones. This is again because both Atchison and Furihata fitted their $`A`$ and $`Z`$ distributions using models without preequilibrium emission, which provide higher values for the excitation energy, $`A`$, and $`Z`$ of fissioning nuclei. This means that to get a good description of $`A`$ and $`Z`$ distributions for fission fragments using GEM2 in CEM2k “with Prec”, we would need to modify the $`A`$ and $`Z`$ distributions of fission fragments in GEM2, making them wider. This would take us beyond the scope of the present work and here we do not vary any more parameters than we have already discussed.

Fig. 4 shows the GSI measurements  of the $`A`$ and $`Z`$ distributions of the kinetic energy of products from the same reaction compared with our CEM2k+GEM2 calculations both with and without preequilibrium emission. Both versions of our calculations are in reasonable agreement with the data.

Mass and charge distributions of the yields or kinetic energies of the nuclides produced show only general trends and are not sensitive enough to the details of a reaction. It is much more informative to study the characteristics of individual nuclides and particles produced in a reaction. Fig. 5 shows a comparison of the experimental data on production yields of thirteen separate isotopes with $`Z`$ lying from 22 to 82 from the same reaction measured at GSI  with our calculations using both the “with Prec” (upper plot) and “no Prec” (middle plot) versions.

The agreement (or disagreement) of our calculations with these data is different from what we have for the integral $`A`$ or $`Z`$ distributions in Figs. 3 and 4: We see that for the isotopes produced in the spallation region (not too far from the target) and for fission fragments in the region with the maximum yield, the version “with Prec” agree much better with the data than the version “no Prec”. Only for production of isotopes at the border between spallation and fission and between fission and fragmentation does the version “with Prec” underestimates the data, due to too narrow $`A`$ and $`Z`$ distributions in the simulation of fission fragments, as we discussed previously. The “no Prec” version agrees better with the data in these transition regions but are in worse agreement for isotopes both in the spallation region and in the middle of the fission region. We conclude that if a model agrees well with some $`A`$ or $`Z`$ distributions it does not necessarily mean that it also describes well production of separate isotopes. In other words, integral $`A`$ and $`Z`$ distributions are not sensitive enough to develop and test such models, a practice which is often used in the literature.

The lower plot in Fig. 5 shows results of calculations with a version of CEM2k+GEM2 with reduced preequilibrium emission. We prefer to discus this version together with results by LAQGSM+GEM2, and will return to this plot later.

It is more difficult for any model to describe correctly the energy dependence for the production cross sections of different isotopes, i.e., excitation functions. We calculated using both the “with Prec” and “no Prec” versions of CEM2k+GEM2 all the excitation functions for the same reaction, p + Pb, for proton energies from 10 MeV to 3 GeV and compared our results with all available data from our compilation referred to here as T-16 Library (“T16 Lib”) . Only several typical examples from our comparison are shown below.

Fig. 6 shows two examples of excitation functions for the production of several isotopes in the spallation (first two columns of plots) and fission (the last two columns of

plots regions). One can see a not too good but still reasonable agreement of both calculations with many data (note that most of the data were measured for natPb targets, while our calculations were done for 208Pb). We see that merging CEM2k with GEM2 allows us to reasonably describe yields of fission fragments, while in the old standard CEM2k we do not have any fission fragments and are unable to describe such reactions at all. We see that as shown in Fig. 5 for a single proton energy of 1 GeV, the “with Prec” version agrees better with the data in the whole energy region both for spallation products and for the production of most of the fission fragments. Only on the border between fission and fragmentation regions (46Sc and 60Co) does the “no Prec” version agree much better with the data than the “with Prec” version; the reason for this we have already discussed. Similar results were obtained for excitation functions of many other isotopes in the spallation and fission regions. On the whole, the version “with Prec” reproduces most of the experimental excitation functions better that the version “no Prec”.

In Fig. 7 we show examples of excitation functions for the production of light fragments (the first two columns of plots), in the fragmentation region, that are produced in CEM2k+GEM2 only via evaporation (the contribution to the yield of these isotopes from fission or deep spallation is negligible), and of nucleons and complex particles up to $`\alpha `$ (the last two columns of plots). We see that with the “no Prec” version, GEM2 reproduces correctly the yields of light fragments 6He, 9Li, 7Be, 13N, and 18F, and not so well the excitation functions for heavier fragments like 22Na. With increasing mass of the fragment, the calculations progressively underestimate their yields. Note, that in , we got very similar results for excitation functions for the p+Au reaction. The version

“with Prec” strongly underestimates the yields of all these fragments, and this is again not surprising, as Furihata developed her model and fitted all parameters without taking into account preequilibrium processes.

Undeniably, the parameters determining the yields of evaporated fragments in GEM (inverse cross sections and Coulomb barriers) could be adjusted to get a good agreement with the data for the yields of light fragments with the version “with Prec” (see, e. g., how Furihata and Nakamura addressed this problem in  for their version of a code without preequilibrium emission). This is not the aim of our present work and we will not do this here. Even if we were to do this, we expect in advance to get similar results to those we got for the “no Prec” version: It would be possible to describe correctly the yields and spectra of light fragments but not of heavy fragments like 24Na and 28Mg. To describe such heavy fragments (not only their yields, but also their spectra) the model would need to be improved further, by considering other mechanisms for heavy fragment production in addition to the evaporation process taken into account by GEM2.

Finally, the two last columns of plots in Fig. 7 show excitation functions for emission of nucleons and complex particles up to $`\alpha `$ for this reaction. Note that the data for these excitation functions are not so extensive and precise as we have for heavier products: many data points were obtained by integration (plus extrapolation) of the spectra of particles measured only at several angles and only for a limited range of energy. But even from a comparison with these sparce and imprecise data we see that the “with Prec” version describes these excitation function better than the “no Prec” version, just as we found in  for the p+Au reaction. This is an expected result as the high Coulomb barriers for heavy nuclear targets oppose evaporation of low energy charged particles and the main contribution to their yields comes from preequilibrium emission from highly excited pre-compound nuclei.

For completeness sake, we show here also an example of results from a calculation with the merged CEM2k+GEM2 code of a reaction on an actinide, p(100 MeV) + 238U. Generally, to get for actinides a proper fission cross section, we need to adjust in GEM2 the parameters $`C(Z)`$ (or, also $`A_0(Z)`$) in Eq. (7), as they were fitted by Atchison to work the best with Bertini’s INC and we have in CEM2k our own INC. As mentioned above, for actinides, Eq. (1) is not used in GEM2 and $`a_f`$ is not used in any calculations, therefore we do not need to adjust $`a_f/a_n`$, for fissioning nuclei with $`Z>88`$. We found that for this particular reaction, p(100 MeV) + 238U, we get with CEM2k+GEM2 a fission cross section in agreement with the data without any adjustments of the parameter $`C(Z)`$ in GEM2, i.e., we can use just the default parameters of GEM2. Nevertheless, our results for other reactions show that for higher energies of the incident protons or for other target-nuclei, the parameter $`C(Z)`$ has to be fitted to get a correct fission cross section when GEM2 is coupled either with CEM2k or with LAQGSM. In addition, we should mention that for reactions on actinides at intermediate or high energies, the parameter $`a_f^{CEM}/a_f^{RAL}`$ should also be fitted along with $`C(Z)`$. In some simulated events several protons can be emitted at the cascade and preequilibrium stages of the reaction, as well as at the evaporation stage, before the compound nucleus actually fissions (also complex particles can be emitted before fission), and the charge of the fissioning nucleus can have $`Z88`$, even when the initial charge of the target has $`Z>88`$. At the same time, for $`Z88`$, due to charge exchange reactions, the charge of the fissioning nucleus may exceed 88, so that we would need to fit as well $`C(Z)^{CEM}/C(Z)^{RAL}`$. This is a peculiarity of treating the fission probability $`P_f`$ differently for the elements above and below $`Z=89`$ in the Atchison model.

Fig. 8 shows mass distributions of products from p(100 MeV) + 238U calculated with both versions of CEM2k+GEM2 compared to the available experimental data  and with results by the phenomenological code CYF of Wahl  (short dashed lines). We need to mention that these data are not as good for testing and developing models as are the GSI data measured in inverse kinematics for the p + Pb reaction discussed above: All the data shown in Fig. 8 were obtained by the $`\gamma `$-spectrometry method. Only some of the produced isotopes were measured, and most of the data were measured for the cumulative yields. To get the “experimental” A-distribution, we summed for each $`A`$ the available data taking care to not sum the individual cross sections already included in some cumulative yields; but the resulting A-distribution is still not complete, as many isotopes were not measured. This means that some theoretical values can be above the experimental data (where some isotopes were not measured) without necessarily implying disagreement between calculations and measurements.

One can see that both the “with Prec” and “no Prec” versions of CEM2k+GEM2 describe equally well the mass distribution of the products from this reaction; therefore it is not possible to choose between either of the versions from a comparison with these not very informative data. We see that results by the phenomenological code CYF of Wahl  also agree well with the A-distribution of measured fission products. It is more

useful to compare calculations with the yields of individually measured isotopes. In Fig. 9, we compare our calculations with all cross sections measured by Titarenko  for each nuclide separately, where we can compare our results with the data (we do not include in our comparison the nuclides measured only either in their isomer or ground states, as our model does not provide such information: CEM2k+GEM2 provides only yields for the sum of isotope production cross sections both in their ground and excited states). We see that on the whole, the “with Prec” version agrees better with most of the individually measured cross sections than the “no Prec” version and for many of the measured isotopes the disagreement is less than a factor of two, with several of our calculated points coinciding with the data within the plotting accuracy. Nevertheless, for several isotopes like 72Ga, 227Thc, 228Pac, and 238Np, we see some big disagreements. For comparison, we also show in Fig. 9 calculations by the phenomenological code YIELDX of Silberberg, Tsao, and Barghouty  and using its 2000 updated version, YIELDX2k , and with the phenomenological code CYF by Wahl65 often used in applications. We see that these phenomenological systematics fail to describe the production of many isotopes from this reaction, indicating that we cannot rely on phenomenological systematics and must develop reliable models to be used in applications.

Very similar results for the studied proton-induced reactions with the ones shown above by CEM2k+GEM2 were obtained also by merging LAQGSM with GEM2. We will not discuss such them here; instead, we will apply the merged LAQGSM+GEM2 code to analyze reactions on the same 208Pb-target (projectile, in the case of the inverse kinematics of the GSI measurements ), but induced by deuterons.

A- and Z-distributions of the nuclides produced in the d(1 GeV/nucleon) + Pb reaction calculated with the LAQGSM+GEM2 code are compared in Fig. 10 with the GSI data . As in case of proton-induced reactions, we get a correct fission cross section only by adjusting the ratio of level density parameters in the fission and evaporation channels $`a_f/a_n`$ in Eq.(5), or, equivalently, by adjusting $`a_f^{CEM}/a_f^{RAL}`$, if we use $`a_f^{RAL}`$ provided by Eq. (5). We get a correct fission cross section for $`a_f^{CEM}/a_f^{RAL}=1.15`$ in the “with Prec” case (short-dashed lines in Fig. 10) and for $`a_f^{CEM}/a_f^{RAL}=0.97`$ in the case of “no Prec” (long-dashed lines in Fig. 10) when calculating with LAQGSM+GEM2 both p- and d-induced reactions on Pb at 1 GeV/nucleon. Note that these values differ slightly from the ones obtained above for the CEM2k+GEM2 code (see Fig. 3). This can be easy understood, as the INC of LAQGSM differs from the INC of CEM2k, therefore the nuclei simulated by LAQGSM+GEM2 and CEM2k+GEM2 to really fission after the cascade and preequilibrium stages of reaction and after evaporation of several particles from compound nuclei before fission have slightly different average values of A, Z, and $`E^{}`$.

The agreement (or disagreement) of the LAQGSM+GEM2 results with these data is similar to what we get for the p+Pb reaction: the “no Prec” version provides a good description of both mass and charge distributions and agrees better with the data for these characteristics than the “with Prec” version. These results, together with the ones discussed above for p+Pb reactions, as well as results for many other p+A and A+A reactions measured at GSI at about 1 GeV/nucleon and analyzed with our codes suggest us that we need to take into account preequilibrium processes, but we need less emission of preequilibrium particles than provided by out standard CEM2k and LAQGSM codes. We address this problem following Veselsky . It is assumed that the ratio of the number of quasiparticles (excitons) $`n`$ at each preequilibrium reaction stage to the number of excitons in equilibrium configuration $`n_{eq}`$, corresponding to the same excitation energy, to be the crucial parameter for determining of the probability of preequilibrium emission $`P_{pre}`$. This probability for a given preequilibrium reaction stage is evaluated using the formula

$$P_{pre}(n/n_{eq})=1\mathrm{exp}\left(\frac{(n/n_{eq}1)}{2\sigma _{pre}^2}\right)$$
(9)
for $`nn_{eq}`$ and equal to zero for $`n>n_{eq}`$. The basic assumption leading to Eq. (9) is that $`P_{pre}`$ depends exclusively on the ration $`n/n_{eq}`$ as can be deduced from the results of Böhning  where the density of particle-hole states is approximately described using

a Gaussian centered at $`n_{eq}`$. The parameter $`\sigma _{pre}`$ is a free parameter and no dependence on excitation energy is assumed . Our calculations of several reactions using different values of $`\sigma _{pre}`$ show that an overall reasonable agreement with available data can be obtained using $`\sigma _{pre}=0.4`$ or 0.5 (see Fig. 11). We choose the fixed value $`\sigma _{pre}=0.4`$ both for CEM2k+GEM2 and LAQGSM+GEM2 for all our further calculations.

Results with reduced preequilibrium emission according to Eq. (9) using $`\sigma _{pre}=0.4`$ are shown in Fig. 10 with solid lines (and in the lower plot of Fig. 5 for p+Pb reactions); they agree better with the data than either “with Prec” or “no Prec” standard results.

Fig. 11 shows a comparison of the experimental data on production yields of ten separate isotopes with $`Z`$ lying from 35 to 82 from the same reaction measured at GSI  with our LAQGSM+GEM2 calculations using both the “with Prec” and “no Prec” (upper plot) versions, and when calculated with the version of reduced preequilibrium emission according to Eq. (9) (lower plot). One can see that as in the case of the p+Pb reaction (see lower plot in Fig. 5), calculations with reduced preequilibrium emission agree much better with experimental data than either standard “with Prec” or “no Prec” versions. Similar results were obtained for all other p+A and A+A reactions measured at GSI and calculated by our codes. Therefore we choose the version with reduced preequilibrium emission with a fixed value of $`\sigma _{pre}`$ of 0.4 as our “standard” version both for CEM2k+GEM2 and LAQGSM+GEM2, and all our further calculations will be done in this approach.

Fig. 12 shows a comparison of the experimental  A-distribution of the mean kinetic energies of the spallation-residue products from the same reaction with different calculations by LAQGSM+GEM2. One can see that for this particular reactions all calculations agree well with the data, with a slightly better agreement for the version with reduced preequilibrium emission. Similar results were obtained for the Z-distribution of products from this reaction.

Finally, for completeness sake, the last four figures (Figs. 13–16) show a comparison of all cross sections measured at GSI for the production of nuclides both from spallation and fission reactions from interaction of 208Pb beams with p and d targets compared with out “standard” LAQGSM+GEM2 results. On can see a very good overall agreement, taking into account that we use here our “standard” version of the code and no further fitting of any parameters is done. Similar results were obtained with standard versions of CEM2k+GEM2 and LAQGSM+GEM2 for many other p+A and A+A reactions measured recently at GSI. Such results will be presented in a separate publication.

Further Work

Merging the Generalized Evaporation Model code GEM2 by Furihata  with our CEM2k and LAQGSM codes allows us to describe reasonably well many fission and fragmentation reactions in addition to the spallation reactions already described well by CEM2k and LAQGSM. But to do so, ones need to provide first a correct description of the fission cross sections by fitting in GEM2 the ration $`a_f/a_n`$ for pre-actinides and/or the parameter C(Z) in Eq. (7) for actinides. This is not a serious problem and we have derived already new approximations for $`a_f/a_n`$ and for parameter C(Z) in Eq. (7) that provide for CEM2k+GEM2 and LAQGSM+GEM2 correct fission cross section calculation both for pre-actinides and actinides for a large range of projectile energies. Such new approximations will be presented together with examples of our results in a separate publication.

Some reactions, like production of fission fragments at the borders between fission and fragmentation or between fission and emission of heavy fragments like Na and Mg are poorly described by CEM2k+GEM2 and LAQGSM+GEM2 in their current versions. This disagreement does not discourage us; the results of the present work suggest that some of the fission and evaporation parameters of GEM2 can be adjusted to get a much better description of all reactions. This lends credibility to such an approach.

There is one more drawback of this approach to be mentioned: considering evaporation of up to 66 particles in GEM becomes extremely time consuming when calculating reactions with heavy targets at high incident energies. But even this disadvantage may be mitigated by the performance of modern computers. We have nevertheless some more serious doubts about the current version of GEM2 related to its lack of self-consistency, e.g.:

1) using different, not physically related parameterizations for inverse cross sections and Coulomb barriers for different particles and fragments;

2) using different level density parameters for the same compound nuclei when calculating evaporation and estimating fission probability from the widths of neutron evaporation and fission;

3) different, and purely phenomenological treatments of fission for pre-actinide and actinide nuclei;

4) not taking into account at all the angular momentum of compound and fissioning nuclei;

5) rough estimations for the fission barriers and level density parameters, etc.

This means that an approach like GEM2 can in principle be used to describe fission and evaporation of particles and fragments heavier that 4He after the INC and preequilibrium parts of CEM2k, LAQGSM, and other models, but it should be considerably improved striving first to progressively incorporate better physics, and only after that adjusting to selected data.

The results of the present work and from  show that on the whole, the merged CEM2k+GEM2 and LAQGSM+GEM2 codes agree better with most of tested experimental data when we take into account the preequilibrium emission of particles, than when we neglect completely preequilibrium processes. But there is still a not so clear question here; we have had some indications for many years that CEM accounts for too many preequilibrium particles, at least at energies above the pion-production threshold. To solve this problem, as a “zero-step” approximation, in our original CEM2k version  we neglected in the exciton model the transitions that decrease or do not change the number of excitons $`\mathrm{\Delta }_n=2`$ and $`\mathrm{\Delta }_n=0`$, shortening in this way the preequilibrium stage of reactions. We do not like this approach as it is somehow arbitrary, “ad hoc”, even though this “never come back” approximation is used in some popular codes like LAHET  and FLUKA . In the present work we removed this arbitrary condition in CEM2k and the “with Prec” version takes into account all the preequilibrium transitions $`\mathrm{\Delta }_n=+2`$, 0, and -2, making the preequilibrium stage of a reaction longer and increasing the number of emitted preequilibrium particles. The results of the present work indicate to us once again that we need to take into account the preequilibrium stage in reactions, but we need less particle emission than we calculate at this stage with the standard version of CEM. In the present paper we address this problem reducing the number of emitted preequilibrium particles following Veselsky , as described above.

Besides GEM2, we have investigated the well known code GEMINI by Charity  as an alternative way to describe production of various fragments by merging GEMINI with both our CEM2k and LAQGSM, and we tested also the thermodynamical fission model by Stepanov  with its own parameterizations for mass and charge widths, level-density parameters, fission barriers, etc., merging it with both CEM2k and LAQGSM to describe fission. In addition, we have started to extend CEM2k and LAQGSM and to develop our own fission model, as briefly noted in . The preliminary results we found for spallation, fission, and fragmentation products from several reactions we tested so far using these approaches are very promising and we will present our results from these studies in several separate papers.

Acknowledgment

We thank Dr. Shiori Furihata for several useful discussions, for sending us her Generalized Evaporation Model code GEM, and providing us further with its updated version, GEM2. We are grateful to Dr. William B. Wilson for providing us with results of his calculation the reaction p(100 MeV) + 238U with the codes CYF and YIELDX2k. We thank Prof. Takashi Nakamura and Drs. Yoshiyuki Iwata and Hiroshi Iwase for sending us numerical values of their measured neutron spectra from heavy-ion reactions and results of their calculations with QMD and HIC.

This study was supported by the U. S. Department of Energy and by the Moldovan-U. S. Bilateral Grants Program, CRDF Project MP2-3025.


format: text

text:
# Breakdown of a Mott insulator — non-adiabatic tunneling mechanism

## Abstract

Time-dependent nonequilibrium properties of a strongly correlated electron system driven by large electric fields is obtained by means of solving the time-dependent Schrödinger equation for the many-body wave function numerically in one dimension. While the insulator-to-metal transition depends on the electric field and the interaction, the metallization is found to be described in terms of an universal Landau-Zener quantum tunneling among the many-body levels. These processes induces current oscillation for small systems, while give rise to finite resistivity through dissipation for larger systems/on longer time scales.

Introduction Nonlinear responses and time-dependent phenomena are an open frontier in the physics of strongly correlated systems. Namely, while the physics of electron correlation has been intensively studied in the last decades in the context of the high-Tc cuprates and related systems, the efforts were concentrated on equilibrium or linear-response properties. Properties beyond the linear response or time-dependent properties have been much discussedNoneq , but remains a challenging problem, which is, theoretically, due to the difficulty in dealing simultaneously with the many-body effect (correlation) and finite current (non-equilibrium), which are both non-perturbative physics. So we definitely require a clear picture for this fundamental problem, especially for low-dimensional materials.

In this Letter we consider the Hubbard model under constant driving forces in one dimension. The ground state of the half-filled Hubbard system is a Mott insulatorLieb  for arbitrary strengths of the electron-electron repulsion $`U>0`$, while the state is metallic when the band filling is shifted away from the half-filling by doping carriers. So the question we pursuit here is: what will happen if we destroy the Mott insulator by applying strong electric field, instead of doping?

There are existing theoretical approaches that employ the Bethe-ansatz method, and the closing of the Mott gap has been discussedFukui ; Deguchi:1998GC . In these studies, however, electric fields are not actually applied, but the left-going and right-going hopping terms are made different instead. It is rather difficult to relate this artificial, and non-Hermitian, model with a system in an electric field. That is why we have here opted for actually applying electric fields for the first time to keep track of the time evolution of the many-body wave functions and levels. For that we make use of a numerical integration of the time-dependent Schrödinger’s equation, where we have also looked at excited states, which should be essential in examining non-equilibrium phenomena.

We have found that the electric field, if strong enough, breaks the Mott-insulator phase. While the critical field strength required for the breakdown of the Mott insulator depends sensitively on the magnitude of the electron-electron interaction, we propose here that the mechanism for the metallization can be viewed as the non-adiabatic tunneling between the many-body levels. We have verified this by confirming numerically that a universal Landau-Zener quantum tunneling governs the nonlinear conduction (e.g., the $`IV`$ characteristics representing the insulator-to-metal transition), where the fit is surprisingly good although the Landau-Zener formalism is originally intended for one-body problems while the problem at hand, being many-body, involves the Hilbert space with huge dimensions.

Just after the metallization, we have a self-induced current oscillation. While this should be realistic for mesoscopic systems, for large systems or on a longer time scale, a novel, “Ohmic” conduction is found to result. This occurs despite the absence of disorder and the heatbath degrees of freedom, but a series of non-adiabatic tunneling among many-body states are responsible for the dissipation effect. Indeed, the breakdown of a one-dimensional (1D) tag  as well as two-dimensional (2D) yama  Mott insulator has been experimentally studied, and, among other interesting phenomena including spontaneous density-pattern formations,yama ; tag ; Kuma  a seemingly Ohmic conduction was found for a rather wide range of external electric field, until the Ohmic conduction is eventually broken.

Formulation If one wants to apply a finite electric field to a system, one problem is how to treat the electrodes. To get rid of ambiguities arising from this, we have opted here for a periodic system (a ring in 1D), where the electric field $`F`$ is applied via a time-dependent AB flux $`\mathrm{\Phi }(t)=eLFt`$ piercing the ring (inset of Fig.1(a); $`L`$: sample length). This will lead to a circular electromotive force due to Faraday’s law. The flux makes the hopping integral in a tight-binding model complex, where the Hamiltonian is

$`H(t)`$ $`=`$ $`{\displaystyle \frac{W}{4}}{\displaystyle \underset{i,\sigma }{}}(e^{2\pi i\mathrm{\Phi }(t)/N}c_{i+1\sigma }^{}c_{i\sigma }+\mathrm{h}.\mathrm{c}.)`$ (1)
$`+`$ $`U{\displaystyle \underset{i}{}}n_in_i.`$
Here $`W`$ is the bandwidth, $`U`$ the electron-electron repulsion, $`N=L/a`$ the total number of sites with $`a`$ being the lattice constant. Hereafter we take the unit in which $`e=a=h=1`$. The electric field in this formalism amounts to a time-dependent deformation of the Hamiltonian.

We have then to solve the time-dependent Schrödinger equation, $`i\frac{d}{dt}|\mathrm{\Psi }(t)=H(t)|\mathrm{\Psi }(t)`$, which governs time evolution of the quantum system at absolute zero temperature, starting from the ground state of the Hamiltonian at $`t=0`$, $`|\mathrm{\Psi }(t=0)|\mathrm{\Psi }_0`$, which is obtained with the Lanczos method. The time-integration of the state vector, which, being many-body, has a huge dimension, requires a reliable algorithm. So we adopt here Cranck-Nicholson’s methodCranckNicholson  that gurantees the unitary time evolution, where the time evolution is put in a form, $`|\mathrm{\Psi }(t+\mathrm{\Delta }t)=e^{i_t^{t+\mathrm{\Delta }t}H(t)}dt|\mathrm{\Psi }(t)[1i\mathrm{\Delta }t/2H(t+\mathrm{\Delta }t/2)]/[1+i\mathrm{\Delta }t/2H(t+\mathrm{\Delta }t/2)]|\mathrm{\Psi }(t),`$ which is unitary by definition. Here the time step is taken to be small enough ($`dt=1.0\times 10^2`$ with the time in units of $`4\mathrm{}/W`$ hereafter) to ensure convergence for $`N10`$ site systems, for which the dimension of the Hamiltonian is $`10^4`$. We have concentrated on the total $`S^z=0`$ subspace with $`N_{}=N_{}=N/2`$.

Result We first plot in Fig.1(a) the result for the expectation value of the current density averaged over the sites,

$$J=\frac{W}{4N}\underset{i,\sigma }{}(ie^{2\pi i\mathrm{\Phi }/N}c_{i+1\sigma }^{}c_{i\sigma }+\mathrm{h}.\mathrm{c}.).$$
(2)
The behavior of $`J(t)`$ for various values of $`U`$ with a fixed value of the electric field $`F`$ is seen to fall upon three regimes when $`U`$ is varied: A perfect metallic behavior ($`J(t)t`$) when the electrons are free ($`U/W=0`$), while when the interaction is strong enough ($`U/W1`$) the current has a zero expectation value. For an intermediate regime of $`U/W`$ we have finite $`J`$’s with some oscillations in the current for finite systems. By contrast, a non-half-filled system has the time evolution distinct from the ground-state behavior. The difference has its root in the spectral property as will be discussed later.

Now, Fig.1(b) plots the time evolution of the current when the electric field $`F`$ is varied with a fixed value of $`U/W`$, again for the half-filled case. The result may be summarized as follows:

Small $`F`$ regime (Mott insulator)
A drastic difference between the half-filled and doped systems appears for small $`F`$. When half-filled, $`J(t)`$ in the limit of $`F0`$ smoothly approaches the $`J(t)=0`$ behavior of the ground state (Mott insulator). Here $`J(t)`$ is the time-averaged current. On top of the $`J(t)=0`$ an oscillatory behavior with the period of $`\mathrm{\Phi }_0(e/h=1`$: flux quantum) is seen, which is nothing but the AB oscillation (a saw-tooth, due to a symmetry about $`\mathrm{\Phi }=\mathrm{\Phi }_0/2`$).

Moderate $`F`$ regime (metal)
In this regime, the current in the half-filled case shows an oscillatory behavior (see typically the $`LF=0.008`$ data in Fig.1(b)).

Large $`F`$ regime (perfect metal)
When the electric field $`F`$ becomes large enough, the system becomes a metal, in which $`J(t)t`$ for $`\mathrm{\Delta }\mathrm{\Phi }<\mathrm{\Phi }_0N/4`$. A further oscillation in $`J(t)`$, with a long period ($`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_0N`$), is seen, which we will discuss later.

The $`F`$-dependence of $`J`$ are displayed in Fig.2 as the $`I`$-$`V`$ characteristics for various values of $`U`$. Here the time average $`J=_0^{\mathrm{\Delta }t}J(t)𝑑t/\mathrm{\Delta }t(LF\mathrm{\Delta }t=N\mathrm{\Phi }_0/4)`$ of the current density $`J(t)`$ is taken over one forth of the extended AB period ($`0\mathrm{\Phi }N\mathrm{\Phi }_0/4`$), since we are interested in the rise in the current which should represent the behavior in the thermodynamic limit. We can see that $`J`$ becomes nonzero rather abruptly at the metallization, where the threshold electric field increases with $`U/W`$, thereby the $`F`$\- dependence becoming weaker. Just after the metallization some oscillation (in the $`F`$\- dependence this time) is seen for finite systems.

Non-adiabatic tunneling In order to understand the physics underlying these time evolutions and their $`F`$ dependence, let us here evoke the notion of non-adiabatic tunneling, originally conceived for one-body problems by Landau, Zener, and Stűckelberg (LZS)Landau ; Zener ; St . The LZS theory considers the quantum tunneling beyond the adiabatic approximation. Namely, when a parameter defining the Hamiltonian is varied infinitely slowly (adiabatically), we can just plot a set of energy levels against the parameter, which in general contain level anticrossings, since the two levels should repel with each other unless they are allowed to cross due to e.g. symmetry reasons. An initial state that start from one of the lines should evolve with time sticking to that line (adiabatic theorem, corresponding to $`p=0`$ in Fig.3(a)). When the parameter is varied with a finite velocity, the state can make a transition across the level anticrossing with a finite probability $`p(0`$ in Fig.3(a)), where the transition is caused by a quantum mechanical tunneling across the gap $`\mathrm{\Delta }E`$. The transition probability $`p`$ depends in the LZS theory on the speed ($`\delta \dot{E}`$) the two energy level approach as

$$p=\mathrm{exp}\left(2\pi \frac{(\mathrm{\Delta }E)^2}{\delta \dot{E}}\right)=\mathrm{exp}\left(2\pi \frac{(\mathrm{\Delta }E)^2}{\frac{d\delta E}{d\mathrm{\Phi }}LF}\right).$$
(3)
Here $`\delta E`$ is the difference between the “unperturbed”, crossing energy levels (dashed lines in Fig.3(a)) and $`\delta \dot{E}d\delta E/dt=(d\delta E/d\mathrm{\Phi })\dot{\mathrm{\Phi }}`$ with $`LF=d\mathrm{\Phi }/dt`$. We can immediately see that the process is non-perturbative in $`F`$, since $`p`$ is singular in $`F`$.

Although the original LZS theory is devised for one-body systems, there is no reason why we cannot apply it to many-body systems, as demonstrated for a spin system by Miyashita et alraedt ; Miya . So here we apply the concept to the Hubbard model, which is, to our knowledge, the first time the LZS theory is applied to interacting electron systems. In order to check the validity of the LZS picture, we have first calculated the transition probability $`p`$. The level anticrossing we focus on is the first one encountered by the ground state at $`\mathrm{\Phi }=\mathrm{\Phi }_0/2`$ in the level flow (marked with a double circle in Fig.3(c)). In Fig.4(a) we plot $`|\mathrm{\Psi }_0|\mathrm{\Psi }(t)|^2`$, the weight of the ground state around the level anticrossing, for various values of $`U`$ and $`F`$. From this we have obtained the transition probability $`p`$ as the asymptotic value of $`|\mathrm{\Psi }_0|\mathrm{\Psi }(t)|^2`$asymptotic . Figure 4(b) plots -Log $`p`$ as a function of the LZS parameter, $`(\mathrm{\Delta }E)^2/(\delta \dot{E})`$, with $`\mathrm{\Delta }E`$ now defined as the interaction-originated one. We can see a remarkably accurate linear dependence on the LZS parameter, which clearly indicates that the LZS theory is applicable to the many-body system we have at hand. So the field $`F`$ and the interaction $`U`$ enter into the problem only as a combination $`(\mathrm{\Delta }E(U))^2/(\frac{d\delta E}{d\mathrm{\Phi }}LF)`$ (where $`\mathrm{\Delta }E(U)`$ is an increasing function of $`U`$; see Fig.3(b)).

For large enough electric field $`F`$ the nonadiabatic tunneling is so effective that the state goes straight across each crossing with probability close to unity, so we end up with a long period ($`\mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }_0N`$), an analogue of the “extended AB period” discussed for electron systems by Kusakabe and two of the present authorsKabe ; Ari , a notion originally proposed for a spin (Heisenberg) system by Sutherland Sut . In the thermodynamic limit ($`N\mathrm{}`$) the period becomes $`\mathrm{}`$.

Encouraged by this, we have then re-plotted the time-averaged current, $`J`$, against the inverse of the LZS parameter $`\frac{d\delta E}{d\mathrm{\Phi }}/(\mathrm{\Delta }E)^2\times LF`$, in Fig.5. Dramatically, all of the curves, which appeared quite different for different values of $`U/W`$ in the raw $`I`$-$`V`$ characteristics (Fig.2), fall on a single, universal curve within a reasonable error when plotted against the LZS parameter. Specifically, a threshold between the insulating behavior and the dissipative metallic one is clearly seen at around $`F0.5\times (\mathrm{\Delta }E)^2/\left(\frac{d\delta E}{d\mathrm{\Phi }}L\right)`$.

After the metallization the current is seen to behave roughly linearly with $`F`$. This seemingly “Ohmic” behavior is far from trivial but rather surprising, since we have a many-body but clean system. This implies that (i) after many level crossings the system reaches a steady state, and (ii) the many-body gaps at these level crossings have similar magnitudes as the first one ($`\mathrm{\Delta }E`$ above). We can in fact recall that the non-adiabatic tunneling is a quantum version of dissipation, where different quantum states become mixed after level anticrossing. To get an insight into the quantum dissipation let us examine the nature of the many-body gap $`\mathrm{\Delta }E(U)`$ in the half-filled Hubbard model, which is a charge gap characteristic to the half-filled Hubbard model and vanishes when doped (Fig.3(c)). The total momenta of the anticrossing states (the ground state and the excited state with one pair of complex charge rapiditiesWoy ) differ by $`2k_F`$ where $`k_F`$ is the Fermi wavenumber (for small $`U`$; replace $`k_F`$ with the quasimomentum for general $`U`$). At half-filling, these Umklapp processes take place and the momentum of the many-body state is dissipated. In other words, the role of heatbath degrees of freedomCaLe  is played by the many-body system itself. Let us add that the current oscillation (typically seen for $`LF=0.008`$ in Fig.1(b)) is due to the kicks from Umklapp processes and should be observed in small systems with strong electron correlation. The threshold electric force, $`F0.5(\mathrm{\Delta }E)^2/[(d\delta E/d\mathrm{\Phi })L]`$, has a similar order of magnitude as the critical field observed in the experimenttag . The threshold electric force, $`F(\mathrm{\Delta }E)^2/[(d\delta E/d\mathrm{\Phi })L]`$, should approach an asymptotic value in the thermodynamic limit, since $`(d\delta E/d\mathrm{\Phi })|_{\mathrm{\Phi }=\mathrm{\Phi }_0/2}((d^2E/d\mathrm{\Phi }^2)|_{\mathrm{\Phi }=0}`$ Drude weight) $`1/L`$ should cancel the $`L`$ in the denominator. Exactly how the thermodynamic limit is reached and how the seemingly Ohmic conduction occurs is an important future problem.

We thank Koichi Kusakabe for discussions, and TO acknowledges Naoki Watanabe for a help on the time-dependent numerical simulation.