text:
# Families of Matter-Waves for Two-Component Bose-Einstein Condensates

## Abstract

We produce several families of solutions for two-component nonlinear Schrödinger/Gross-Pitaevskii equations. These include domain walls and the first example of an antidark or gray soliton in the one component, bound to a bright or dark soliton in the other. Most of these solutions are linearly stable in their entire domain of existence. Some of them are relevant to nonlinear optics, and all to Bose-Einstein condensates (BECs). In the latter context, we demonstrate robustness of the structures in the presence of parabolic and periodic potentials (corresponding, respectively, to the magnetic trap and optical lattices in BECs).

The recent progress in experimental and theoretical studies of Bose-Einstein condensates (BECs)  has made matter-wave solitons physically relevant objects. One dimensional (1D) dark  and bright  solitons have been observed in experiments, and possibilities for the observation of their multidimensional counterparts were predicted . Further study of matter-wave solitons is a subject of profound interest, not only from a theoretical perspective, but also for applications, as there are possibilities to coherently manipulate such robust structures in matter-wave devices, e.g., atom chips, which are analogs of the existing optical ones . On the other hand, many results obtained for optical solitons as fundamental nonlinear excitations in optical fibers and waveguides (see, e.g., recent reviews ) suggest the possibility to search for similar effects in BECs.

A class of physically important generalizations of the nonlinear Schrödinger (NLS) equation for optical media, or its BEC counterpart, the Gross-Pitaevskii equation (GP), is based on their multi-component versions. In particular, the theoretical work has already gone into studying ground-state solutions  and small-amplitude excitations  of the order parameters in multi-component BECs. Additionally, the structure of binary BECs , including the formation of domain walls in the case of immiscible species, has also been studied ; 1D bound dark-dark  and dark-bright  soliton complexes, as well as spatially periodic states , were predicted too. Experimental results have been reported for mixtures of different spin states of 87Rb  and mixed condensates . Efforts were made to create two-component BECs with different atomic species, such as 41K–87Rb  and 7Li–133Cs .

In this work, we report novel solitons in the context of coupled two-component GP equations. These solutions correspond to new families of solitons even for the NLS equations per se, hence they are also interesting as nonlinear waves in their own right. We start by demonstrating their existence in the context of two coupled NLS equations. Some of them are relevant as new solitons in nonlinear-optical models as well. We also demonstrate that all the solutions proposed herein persist in the presence of the magnetic trap and optical lattice (OL), i.e., parabolic and sinusoidal potentials , which are important ingredients of experimental BEC setups.

Assuming that the nonlinear interactions are weak relative to the confinement in the transverse dimensions, the transverse size of the condensate is much smaller than its lengths. In this case, the BEC is a “cigar-shaped” one, and the GP equations take an effectively 1D form :

$$i\frac{u_j}{t}=\frac{m_1}{m_j}\frac{^2u_j}{x^2}+\underset{k=1}{\overset{2}{}}a_{jk}|u_k|^2u_j+V_j(x)u_j,(j=1,2)$$
(1)
where $`u_j(x,t)`$ are the mean-field wave functions of the two species, $`t`$ and $`x`$ are, respectively, measured in units of $`2/\omega _1`$ and the transverse harmonic-oscillator length $`a_1\sqrt{\mathrm{}/(m_1\omega _1)}`$ ($`m_j`$ and $`\omega _j`$ are the mass and transverse confining frequency of each species), while $`\mathrm{}\omega _1/2`$ is the energy unit. The coefficients $`a_{jk}`$ in Eq. (1), related to the three scattering lengths $`\alpha _{jk}`$ (note that $`\alpha _{12}=\alpha _{21}`$) through $`a_{jk}=4\pi m_1(\alpha _{jk}/a_1)(m_j+m_k)/(m_jm_k)`$, account for collisions between the atoms belonging to the same ($`a_{jj}`$) and different ($`a_{jk}`$, $`jk`$) species; they are counterparts of the, respectively, self-phase and cross-phase modulation in nonlinear optics. While in optics only specific ratios of the nonlinear coefficients are relevant (such as $`a_{ij}/a_{ii}=2`$ or $`a_{ij}/a_{ii}=2/3`$ ), in the BEC context the interactions are tunable , especially because they can be modified by means of the Feshbach resonance (i.e., by magnetic field affecting the sign and magnitude of the scattering length of the interatomic collisions) . The Feshbach resonance allows one to switch between attractive and repulsive interaction , and even to switch it periodically in time, by means of an ac magnetic field, which allows to create a self-confined 2D BEC without the magnetic trap .

In this work, we consider the case with $`m_1=m_2m`$ and $`a_{11}=a_{22}`$, which corresponds to the most experimentally feasible mixture of two different hyperfine states of the same atom species, or, approximately, to different isotopes of the same alkali metal, trapped in the potential including the magnetic trapping and OL components:

$$V_1(x)=V_2(x)V(x)=\left(\mathrm{\Omega }^2/2\right)x^2+V_0\mathrm{sin}^2(kx+\varphi ).$$
(2)
In Eq. (2), $`\mathrm{\Omega }^22\omega _x^2/\omega _{}^2`$ ($`\omega _{1x}=\omega _{2x}\omega _x`$ are the confining frequencies in the axial direction) and $`V_0`$ (measured in units of the recoil energy ) set the respective potential strengths, $`k`$ is the wavenumber of the interference pattern of the laser beams forming the OL, and $`\varphi `$ is an adjustable phase parameter ($`\varphi \{0,\pi /2\}`$). To estimate physical parameters, we resort to a mixture of two different spin states of 87Rb, confined in a trap with the transverse frequency $`\omega _1=183`$ rad/s, which implies that the length and time units are $`2\mu `$m and $`5.46`$ ms, respectively.

We consider a rather general case, in which the inter-atomic interactions in the first species are repulsive \[therefore, we will use the normalization $`a_{11}+1`$ in Eqs. (1)\], while in the other species they may be either attractive or repulsive. As concerns the interactions between the different species, they are, typically, repulsive. Nevertheless, in the case of two different spin states of the same atom species, the Feshbach resonance between such states is possible too (experimental studies of the Feshbach resonance in this case are currently in progress ), therefore attractive inter-species interactions may be relevant, and this case is also considered below. The solutions reported herein, and their existence and stability regimes are summarized in Table 1.

*Types of solitons* *Existence* *Stability*     $`a_{22}`$ $`a_{12}`$ $`a_{22}`$ $`a_{12}`$  Domain wall $`+`$ $`+`$ $`+1`$ $`>1`$  Dark-antidark $`+`$ $`+`$ $`+1`$ $`(0,0.7]`$  Dark-gray $`+`$ $``$ $`+1`$ $`[0.83,0)`$  Bright-antidark $``$ $``$ $`1`$ $`(1,0)`$  Bright-gray $``$ $`+`$ $`1`$ $`>0`$

Table 1: Existence and stability of structures in the binary BEC. In the “existence” column, $`+/`$ indicates the repulsive/attractive character of the respective inter-atomic interaction which is necessary for the solution to exist. The “stability” column indicates the sign of the coefficient $`a_{22}`$ (we normalize $`a_{11}+1`$, and set $`a_{22}=\pm 1`$) and an interval of the values of $`a_{12}`$ for which the solution is stable.

In most cases the existence and stability of the solution families is investigated numerically. The numerical method was implemented as follows: we first seek stationary solutions by means of Newton iterations which are applied to the steady-state equations $`\mu u_j=u_{j,xx}+_{k=1}^2a_{jk}|u_k|^2u_j+V(x)u_j`$ ($`\mu `$ is the chemical potential). Subsequently, we perform the linear-stability analysis of the obtained soliton solutions $`u_j^{(0)}(x)`$, setting the perturbed solution to be $`u_j=e^{i\mu t}\left[u_j^{(0)}(x)+ϵ\left(b_je^{i\omega t}+c_je^{i\omega ^{}t}\right)\right]`$, where $`\omega \omega _r+i\omega _i`$ is a (generally, complex) perturbation eigenfrequency. Then, the ensuing linear stability problem  is solved for the eigenfrequencies and eigenfunctions $`\{b_j,c_j^{}\}`$. Whenever the solution is unstable, we also examine its evolution in direct simulations of the full equations (1), using a fourth-order Runge-Kutta time integrator with the time step $`dt=0.005`$ ($`27.3\mu `$s in physical units). To initiate the instability development, a uniformly distributed random perturbation of amplitude $`10^4`$ was typically added to an unstable solution.

We now examine in detail the solutions shown in Table 1. First, in the absence of the external potential, a family of domain walls can be found in an exact form for the special case, $`a_{12}=3a_{11}=3a_{22}`$:

$$u_j(x,t)=Ae^{i\mu t}\left[1+(1)^j\mathrm{tanh}(\eta x)\right],$$
(3)
where the chemical potential is $`\mu =4a_{12}A^2`$, and $`\eta ^2=2a_{12}A^2`$ (they follow the pattern of domain-wall solutions found long ago in the context of coupled Ginzburg-Landau equations ). These solutions exist only if $`a_{12}>0`$ and $`\mu >0`$. Similar patterns were found in Ref.  and other related structures were also predicted to occur in higher dimensions .

We have confirmed the existence and stability of the domain walls by direct numerical simulations (not shown here), using numerical continuation to extend them to the case $`a_{12}3a_{11}`$, where the analytical solution is not available. We have thus found that the domain walls exist and are stable for values of $`a_{12}`$ down to $`a_{12}=1`$. The case $`a_{12}=2`$ is relevant to nonlinear optics; stability of the domain wall family for this case was suggested by recent numerical results obtained for a similar discrete coupled-NLS model . Here, we find that these solutions are robust as well for other values of $`a_{12}`$, and, as will be shown below, also in the presence of the external potential in Eq. (2).

Proceeding to the other solutions in Table 1, we first note that, upon fixing $`a_{11}+1`$, we set $`a_{22}=\pm 1`$ for the repulsive and attractive interactions in the second species. Although setting $`\left|a_{22}\right|`$ to be $`1`$ formally limits the generality of the results, it has been checked that taking $`\left|a_{22}\right|1`$ produces similar results to those displayed here.

Given that we have assumed the first component as being always self-repulsive, we look for solutions starting from the uncoupled limit ($`a_{12}=0`$) by taking an initially uniform distribution in this component, $`u_11`$. If the second component is also repulsive, we take a dark soliton as its initial configuration. If it is self-attractive, a bright soliton is initially set in it.

The case of the self-repulsion and dark soliton in the second component, with the inter-species coefficient being *repulsive* too, $`a_{12}>0`$ gives rise to a stationary antidark soliton (i.e., a hump on top of a nonzero flat background) in the first component, see the left panel in Fig. 1. It is easy to understand this structure, as the atoms in the first component, being repelled by the matter in the second one, concentrate in an effective potential well generated by the dark soliton (void) in it. Antidark solitons are well-known to occur when higher-order effects (such as third-order dispersion) or a saturable nonlinearity are present in the single-component NLS equations, which is possible in optics . In that case, the antidark soliton is usually described by a KdV-type asymptotic equation (for the elevation on top of the flat background), and is not stationary, running at the respective velocity of sound . The two-species soliton with the antidark component, presented here, is the first example of a stationary antidark soliton, that we are aware of, in a model without higher-order nonlinearities and dispersions. It is also the *first prediction* of antidark solitons in BECs, which suggests that an experimental verification would be of particular interest.

If, on the contrary, the interaction between the two self-repelling species is *attractive* ($`a_{12}<0`$; right panel in Fig. 1), then the void (dark soliton) in the second component effectively repels the matter in the first one, and thus generates a dip, i.e. a gray soliton in it (for a detailed description, see Ref. ). Such solitons exist in the regular NLS equation, but there (as well as in other instances of their presence that we are aware of) they travel at a nonzero speed (the faster the shallower the dip is), while here the gray solitons are stationary.

The bound dark-antidark two-component states persist for $`0<a_{12}0.7`$, while the dark-gray bound-state branch continues down to $`a_{12}=1`$, getting appreciably unstable in the region $`1<a_{12}<0.83`$. In the latter case, the time evolution (not shown here) leads to formation of *moving* dark-gray soliton complexes.

Similar results were obtained for the last two branches of the new solutions: if the second component is now self-attractive, then the attractive interaction between the BEC species ($`a_{22}=1,a_{12}<0`$; see left panel in Fig. 2) gives rise to bright-antidark solitons (whose existence is explained by the fact that the bright soliton attracts matter in the other component). Such a type of two-component solitons was predicted in optics , but there it cannot exist without the third-order dispersion. In our case, this solution branch terminates at $`a_{12}=1`$, since both states become completely flat.

Lastly, if the bright soliton repels matter in the other component, it naturally induces a dip in it, thus generating a bright-gray soliton in the case of $`a_{22}=1,a_{12}>0`$ (right panel Fig. 2). It is noteworthy that this bright-gray branch is extremely robust; we were able to follow such stable solutions down to $`a_{12}=3`$.

We have also examined these novel solutions in the presence of the OL \[i.e., for $`V_00`$, $`\mathrm{\Omega }=0`$ in Eq. (2)\]. We have found them to be typically stable when $`\varphi =0`$, and unstable when $`\varphi =\pi /2`$ (as might be expected, since the solution is then placed, respectively, at a minimum and maximum of the potential). However, there are exceptions to this rule. For instance, dark-antidark bound solitons for $`a_{12}=0.5`$ are unstable at $`\varphi =0`$, in intervals $`0.04V_00.055`$ and $`0.09V_00.29`$. This oscillatory instability, involving a quartet of eigenvalues, will be examined in detail elsewhere.

We have also identified all the solution branches in the presence of the magnetic trap, as well as in the case when both the magnetic trap and the OL are present. The branches are extremely robust in the presence of the magnetic trap. In particular, Fig. 3 shows the branches for the bright-antidark bound state (left panel) and domain walls (right panel), which are always stable in the presence of the combined potential with $`\mathrm{\Omega }^2=0.02`$ and $`V_0=0.5`$. In particular, for a mixture of two different spin states of 87Rb, the confining frequencies corresponding to this value of $`\mathrm{\Omega }`$ are $`\omega _x=18.3`$ rad/s and $`\omega _1=183`$ rad/s in the transverse and axial directions, respectively. Then, the four cases shown in Fig. 3 correspond to the mixture containing, in the first and second species, $`2\times 10^4`$ and $`6\times 10^3`$ atoms ($`a_{12}=0.5`$), $`2\times 10^4`$ and $`3\times 10^3`$ atoms ($`a_{12}=1`$), $`7.4\times 10^3`$ in each component ($`a_{12}=3`$), and, finally, $`7.6\times 10^3`$ atoms in each component ($`a_{12}=2`$), respectively.

Finally, we notice that the increasing interaction between the components can play a stabilizing role for solutions that are unstable in the OL with $`\varphi =\pi /2`$ (Fig. 4). The dark-antidark coupled state serves as an example, being unstable for $`0<a_{12}<0.4`$ and stable for $`a_{12}>0.4`$ (left panel). The instability evolution is shown in the right panel for $`a_{12}=0.2`$: the dark soliton becomes mobile and starts to oscillate in the combined potential. For the physical example mentioned above, the numbers of atoms in the two components corresponding to the latter case are $`10^4`$ and $`1.6\times 10^4`$, respectively.

In conclusion, we have presented a number of novel families of composite solutions in the generic two-species BEC model. In particular, the first possibility to create an antidark soliton in BECs is predicted. In most cases, the compound solitons and domain walls are very robust, keeping the stability in the presence of the parabolic and periodic potentials. In some cases, the new solutions are relevant to optical models as well. It would be of interest to look for such states experimentally.

Acknowledgments: We appreciate a valuable discussion with D.S. Hall. This work was supported by a UMass FRG, NSF-DMS-0204585 and the Eppley Foundation (PGK), the Special Research Account of the University of Athens (HEN, DJF), the Binational (US-Israel) Science Foundation under grant No. 1999459 (BAM) and the San Diego State University Foundation (RCG).


format: text

text:
# Collective resonance modes of Josephson vortices in sandwiched stack of Bi2Sr2CaCu2O8+x intrinsic Josephson junctions

## Abstract

We observed splitting of the low-bias vortex-flow branch in a dense-Josephson-vortex state into multiple sub-branches in current-voltage characteristics of intrinsic Josephson junctions (IJJs) of Bi2Sr2CaCu2O8+x single crystals in the long-junction limit. Each sub-branch corresponds to a plasma mode in serially coupled Josephson junctions. Splitting into low-bias linear sub-branches with a spread in the slopes and the inter-sub-branch mode-switching character are in good quantitative agreement with the prediction of the weak but finite inter-junction capacitive-coupling model incorporated with the inductive coupling. This suggests the importance of the role of the capacitive coupling in accurately describing the vortex dynamics in serially stacked IJJs.

Ever since the Josephson coupling was known to be established in Bi2Sr2CaCu2O8+x (Bi-2212) single crystals along the adjacent superconducting CuO2 bilayers much research interest has been focused on the dynamics of Josephson vortices in these intrinsic Josephson junctions (IJJs) in the long junction limit Int ; Inductive ; Inductive2 . The interest stems mainly from that the Josephson vortex motion in stacked IJJs provides an ideal model to study the coupled nonlinear dynamic phenomena and also provides unique possibilities of applying stacked IJJs to high-frequency active devices such as sub-millimeter oscillators and mixers.

In Bi-2212 IJJs, superconducting CuO2 bilayers are much thinner than the $`c`$-axis London penetration depth, $`\lambda _{ab}`$ ($``$200 nm) Waldmann . Then, in a transverse magnetic field, screening supercurrents are induced in the stacked superconducting layers, which give rise to the inductive coupling between adjacent IJJs Clem . This coupling is strong enough that it leads to mutual phase locking of Josephson vortex motion along the $`c`$ axis of stacked junctions. It has been predicted that inductively coupled and thus collectively moving vortices may resonate with the plasma oscillation modes and induce sub-branch splitting of current-voltage ($`IV`$) curves near the bias voltages corresponding to plasma mode velocities, which was confirmed by previous numerical analysis and observations Sub . This phase-locked coherence of the plasma oscillation may be exploited to develop sub-millimeter oscillator devices.

On the other hand, the thickness of CuO2 bilayers is comparable to the Debye charge screening length $`D`$ so that the local charge neutrality in the bilayers breaks down Machida2 . This nonequilibrium charge variation in superconducting bilayers yields adjacent IJJs are capacitively coupled as well. The purely capacitive coupling effect was treated theoretically before Koyama . Recently, a model of the capacitive coupling incorporated with the strong inductive coupling has been proposed, where the capacitive coupling is taken into account on an equal footing as the inductive one Kim . In spite of the relative weakness, the capacitive coupling is predicted to lead to significant changes in the collective vortex dynamics. The model predicts an $`IV`$ curve splits into linear multiple sub-branches with varied slopes, which corresponds to different plasma modes, in the low-bias region of the collective Josephson vortex flow. The spread of the sub-branch slopes is governed by the strength of the capacitive coupling, represented by a parameter $`\alpha `$ \[=$`ϵD^2/(st)`$; $`ϵ`$ is the dielectric constant of the insulating layer, $`s`$ (=0.3 nm) and $`t`$ (=1.2 nm) are the thickness of the superconducting and the insulating layers, respectively\], which provides a convenient means to determine the strength of the capacitive coupling in IJJs.

The dynamics of Josephson vortices in stacked IJJs is known to be sensitive to the vortex density in the junctions Koshelev0 . The collective resonance is usually revealed in a high field range of $`HH_d`$ \[$`=\varphi _0/2(s+t)\lambda _J`$; $`\varphi _0`$ is the flux quantum and $`\lambda _J`$ is the Josephson penetration depth\]. In this dense-vortex state the non-Josephson-like emission Hechtfischer  and the Shapiro resonance steps Latyshev2  have been observed, confirming the coherence of the Josephson-vortex-lattice motion over the whole stacked IJJs. Especially, for a higher field range of $`HH_{hd}`$ (=$`\varphi _0/(s+t)\lambda _J`$), the transverse inter-vortex spacing becomes close to $`\lambda _J`$ Bulaevskii . In this range the Josephson current along the length of a junction distributes almost sinusoidally so that the collective resonance is expected to be stronger. Little experimental studies, however, have been done for the collective motion in this highly dense-vortex region.

In this study, we investigated details of multiple plasma modes manifested in the Josephson vortex-flow branch (JVFB) of stacks of IJJs in a high magnetic field above $`H_d`$. In all the three samples examined sub-branch splitting was seen to start occurring for a field exceeding $`H_d`$. Especially, in the JVFB for a field above $`H_{hd}`$, we obtained clear hysteretic multiple sub-branches and mode switching between them. The linear sub-branches for different collective vortex modes exhibited a spread in the slopes. The values of the capacitive coupling parameter $`\alpha `$, estimated from the spread, for the three samples under study were within the range of the theoretical prediction of the hybrid-coupling model Kim . In addition, mode velocities estimated from the inter-branch switching voltages were in reasonable quantitative agreement with the prediction of the model.

Slightly overdoped Bi-2212 single crystals used in this study were grown by the solid-state-reaction method Nam Kim . Samples were micropatterned on the crystals into the geometry of a stack of junctions sandwiched between top (400 nm thick) and bottom (100 nm thick) Au electrodes \[see the inset of Fig. 1(a)\]. The sizes of the three stacks were 17$`\times `$1.5 $`\mu `$m2 (SS1), 16$`\times `$1.5 $`\mu `$m2 (SS2), and 16$`\times `$1.4 $`\mu `$m2 (SS3). Since each superconducting bilayer in a stack of IJJs is much thinner than $`\lambda _{ab}`$, Josephson vortices in a usual mesa structure Doh  are supposed to be strongly coupled to those in the basal stack underneath the mesa, which may significantly distort the dynamics of the Josephson vortices in the mesa itself. In this study the basal stack was removed by using the double-side-cleaving technique Wang  to overcome this problem. The detailed sample fabrication process is described elsewhere Bae . Eliminating the basal part turned out to be essential to obtain ideal vortex-flow characteristics of stacked IJJs. This sandwiched-stack geometry allowed to examine the collective vortex dynamics over the whole thickness of a stack. The magnetic field was aligned in parallel with the plane of junctions within the resolution of 0.01 degree to minimize the pinning of Josephson vortices by the formation of pancake vortices in CuO2 bilayers Hechtfischer2 . All the measurements were done at 4.2 K in a two-terminal configuration \[Fig. 1(a)\]. The combined contact resistance for bottom and top interfaces ($``$2 k$`\mathrm{\Omega }`$ for SS1 and SS2, and $``$600 $`\mathrm{\Omega }`$ for SS3) was subtracted numerically.

The sample parameters are listed in Table I. The number of IJJs in a sandwiched stack, $`N`$, was estimated from the number of zero-field quasiparticle curves (not shown). The McCumber parameter $`\beta _c`$ was obtained from the zero-field return current $`I_r`$ \[as denoted by arrows in Fig. 1(a) for decreasing bias\]. Fig. 1 shows the JVFB of sample SS1 with a single up-down bias sweep. The values of $`H_d`$ and $`H_{hd}`$ for SS1 estimated with $`\lambda _J`$ were 2.2 and 4.4 T, respectively. Thus, the data in Figs. 1(a) and 1(b) correspond approximately to the two dense-vortex ranges, $`HH_d`$ and $`HH_{hd}`$, respectively. The average slope, $`\mathrm{\Delta }V/\mathrm{\Delta }I`$, of the JVFB keeps increasing in proportion with the magnetic field intensity. A few voltage jumps are seen in the JVFB for fields beyond $``$2.6 T in the bias range above $``$40 mV. This unusual feature, taking place in the vortex-flow state for a bias below the quasiparticle-state return current $`I_r`$ in each field, is caused by switching of coherent vortex-lattice motion between different plasma modes. Coherently moving vortex lattice remains in a certain transverse plasma mode as long as the vortex-lattice velocity is smaller than the propagation mode velocity of the plasma oscillation. Exceeding the mode velocity for a higher bias current, however, the resonating dynamic state of the Josephson vortex lattice suddenly switches to the adjacent plasma mode with a higher propagation velocity, giving rise to a voltage jump observed. The onset field of the observed voltage jumps almost coincides with the value of $`H_d`$.

Clearer voltage jumps in the JVFB are seen for a higher vortex density in the field range of $`HH_{hd}`$ as in Fig. 1(b), as the resonant coupling between the vortex-lattice motion and the plasma oscillation is strengthened. In this field range the hysteresis in the quasi-particle branches for $`V>V_c^{vf}`$ (denoted by the arrows) in Fig. 1(b) is almost completely suppressed while the voltage jumps become more evident in JVFB. The inset of Fig. 1(b) shows the distinct voltage jumps in the lower bias region for 5.2 T. One notes that, in Fig. 1(b), all the voltage jumps take place in the low-bias flux-flow region bounded by $`V_c^{vf}`$. The voltage interval between the neighboring jumps tends to increase for higher biases, which is consistent with the prediction of both the purely inductive and the inductive-capacitive hybrid coupling models.

Fig. 2(a) displays the JVFB of the sample SS2 (with $`N`$=25) taken for $`H`$=5.8 T. For the measurement the bias current was repeatedly swept up and down so that detailed hysteretic feature of the sub-branches were traced out. These multiple-sweep data clearly indicate that the voltage jumps observed in the single-sweep data as seen in Fig. 1(b) were indeed caused by the inter-sub-branch mode switching. Since the $`H_{hd}`$ for this sample is 5.4 $`T`$ the data in Fig. 2 correspond to the highly dense vortex state of $`HH_{hd}`$, expectedly with many resonant vortex states as in the sample SS1 in the similar field range. About 22 multiple sub-branches were obtained in SS2, which was close to the number of the IJJs $`N`$=25 in the stack. Details of the vortex-flow region of 5.2 T curve of SS1, as shown in Fig. 1(b), are also illustrated in Fig. 2(b) for a quantitative analysis.

In Figs. 2(a) and 2(b), all the discernible sub-branches exhibit linear $`IV`$ characteristics in the low-bias vortex-flow region for $`V<V_c^{vf}`$, the occurrence of which is in agreement with the prediction of the inductive-and-capacitive hybrid coupling model. By contrast, the purely inductive coupling yields sub-branch splitting near the voltages corresponding to the plasma mode velocities. As denoted by the grey region in Fig. 2 the linear sub-branches, when extrapolated to the low-bias region, converges to a single intercept point with a positive value $`I_p`$ on the current axis. It confirms the vortex-flow nature of the sub-branches with finite vortex pinning, that may have been caused by any defects in the IJJs or by any pancake vortices formed in CuO2 layers by a misaligned field component. The pinning current $`I_p`$ was 2.07 $`\mu `$A and 0.64 $`\mu `$A in SS1 and SS2, respectively. For an analysis of the finite spread in the slopes of the linear sub-branches we adopt fitting to the hybrid coupling model Kim , while incorporating the effect of the finite vortex pinning (represented by $`I_p`$) and the change of the vortex density with applied magnetic fields ($`H_{ext}`$) as

$$V=N(II_p)[\frac{R_n}{A_n}+\frac{kH_{ext}}{I_c}].$$
(1)
Here, $`A_n`$=$`1+2\alpha [1\text{cos}(n\pi /(N+1))]`$, $`I_c`$ is the critical Josephson current at a given applied field, and $`k`$ is a constant. Since samples are in the quasiparticle state for $`V>V_c^{vf}`$ we define the current corresponding to $`V_c^{vf}`$ in a given field as the critical current $`I_c`$ in Eq. (1). Also, the vortex velocity for $`I>I_p`$ is assumed to increase linearly with the bias current.

Without the capacitive-coupling effect ($`\alpha `$=0) the low-bias sub-branches with spreaded slopes, as exhibited by the grey regions in Figs. 2(a) and 2(b) for SS1 and SS2, respectively, are supposed to reduce to single lines (pointed by arrows), where the value of $`k`$ determines the slope of the lines. Increasing $`\alpha `$ the spread of sub-branch slopes increases from the $`\alpha `$=0 lines to the left. In the very low bias region the JVFB in both samples becomes chaotic.

The parameter $`\alpha `$ was determined from the best fit of the spread to Eq. (1), where the normal-state resistance $`R_n`$ was fixed separately following the fitting procedure in Ref. Won . The best-fit values of $`\alpha `$ for SS1 at 5.2 T, SS2 at 5.8 T, and SS3 at 4.15 T were 0.35, 0.17, and 0.24, respectively. These values are in the range of theoretical estimate of the hybrid coupling model Kim , $`0.10.4`$, and in reasonable agreement with the recently observed values of $`0.360.44`$ in SmLa1-xSrxCuO4-δ Helm . The corresponding values of the $`k`$ for SS1, SS2, and SS3 were 0.32, 0.14, and 0.38 mV/T, respectively.

Fig. 3 displays the mode-switching voltages (filled circles) and boundary voltage $`V_c^{vf}`$ (filled square) for SS2 as a function of the mode index $`n`$. Out of 25 expected modes corresponding to $`N`$ the three lowest-index modes were not observed, however, presumably because they were located beyond $`V_c^{vf}`$. The solid line is the trace calculated with $`\alpha `$=0.17 and the best fit value for the Swihart velocity $`c_0=`$1.05$`\times `$105 m/s. The resonance voltage of the $`n`$th mode is obtained from the relation Hechtfischer  $`V_n=Nc_nH_{ext}(t+s)`$. The mode velocity of the transverse plasma oscillation ($`c_n`$) is expected to be enhanced by a factor of $`\sqrt{A_n}`$ from that ($`c_n^0`$) of the purely inductive coupling mdoel ($`A_n=1`$) by the nonvanishing capacitive-coupling effect Kim  as in the relation

$$c_n=c_0\sqrt{\frac{A_n}{1\text{cos}(\pi n/(N+1))}}=c_n^0\sqrt{A_n},$$
(2)
where the mode index $`n`$ runs from 1 to $`N`$. The dotted line corresponds to the case of purely inductive coupling ($`\alpha `$=0) with the same Swihart velocity as used for the solid line. One sees that the enhancement of the mode velocities due to nonvanishing capacitive coupling becomes more evident for high mode indices, giving a maximum enhancement for $`c_N`$ as $`c_N=c_N^0[(1+4\alpha )]^{1/2}`$. The inset of Fig. 3 illustrates the best fit of the mode voltages to the purely inductive-coupling model, with the best-fit value of $`c_0`$ to be 1.35$`\times `$105 m/s. In this case, mode voltages for low-index modes are seen to fall significantly below the model expectation. This trend was also observed previously in the motion of microwave-induced Josephson vortices, although the vortex density was much lower than in this study Doh .

The Swihart velocity for SS2, estimated using the resistively shunted junction model, was $`c_0=2\pi f\lambda _J`$=8.7$`\times `$ 104 m/s. In the calculation we used the Josephson plasma frequency $`f_{pl}`$ (=$`\sqrt{I_c/2\pi \mathrm{\Phi }_0C_J}`$) of 55 GHz, estimated from the critical current $`I_c`$=0.358 mA and the junction capacitance $`C_J`$=9 pF determined from $`\beta _c`$ and $`R_n`$. The value of $`c_0`$ is close to the one obtained for the hybrid-coupling model but tends to deviate more from the predication of the inductive-coupling model.

In summary, the observed characteristics of the sub-branch splitting in the vortex-flow region were in reasonable quantitative agreement with the inductive coupling incorporated with the capacitive coupling. The multiple sub-branches in JVFB were shown to be caused by the inter-junction capacitive coupling and the value of the coupling parameter $`\alpha `$ estimated from the spread of multiple branches was in good agreement with the theoretical prediction of the hybrid coupling as well as the previous observation in optical measurement Helm . The quantitative identification of the multiple sub-branches in this study confirmed the significance of the capacitive coupling in addition to the inductive coupling in accurately describing the dynamics of Josephson vortices in serially stacked IJJs. For further detailed analysis of the Josephson vortex dynamics consideration of the charge-hole imbalance coupling Ryndyk1  and the dissipation effect of quasiparticles in the $`ab`$ planeKoshelev  may be required.

We thank Ju H. Kim for valuable discussions. This work was supported by the National Research Laboratory project administrated by KISTEP.


text:
# Analysis of kaon spectra at SIS energies - what remains from the KN potential? 11footnote 1Supported by IN2P3/CNRS and GSI.

## 1 Introduction

One key question in the analysis of sub-threshold kaon production is how to obtain information on the properties of strange mesons in dense nuclear matter. Especially the relation of the optical potential of $`K^+`$ and $`K^{}`$ in nuclear medium to experimental observables like $`K^{}/K^+`$ ratios, mesonic in-plane flow and azimuthal distribution of kaons is subject of vivid discussions, who have triggered a lot of activities on the experimental  and theoretical side .

In this article we study the production of $`K^+`$ in the reaction Au+Au at 1.48 AGeV which corresponds to recent experiments performed by the KaoS collaboration  and by the FOPI collaboration . For this purpose we use the IQMD model  where we have supplemented our standard simulation program by all relevant cross sections for kaon production and annihilation and a (density and momentum dependent) $`KN`$ optical potential. For the latter we use a parametrization resulting from relativistic mean field calculations of Schaffner-Bielich . Detailed description of our simulations can be found in .

## 2 Time evolution of kaon production

Let us first take a short view on the reaction dynamics of a central Au+Au collision at 1.48 AGeV. Figure 1 shows on top of the left hand side the time evolution of the normalized baryonic density (full line) taken at the center of the reaction in a b=0 fm collision. We see a maximum at about 8 fm/c after first contact (0 fm/c) yielding about 2.5 times ground state densities. In the same figure the dotted line shows the time evolution of the transverse energy at the center of the reaction which can be taken as an indicator of temperature. The transverse energies are divided by 0.1 GeV, thus a value of 1 corresponds to 100 MeV. We see a strong rise at very early times (up to about 4 fm/c) followed by a plateau and a moderate decrease. We see that $`K^+`$ are produced quite early (dashed line) and that the maximum of their production co-incidents with the maximum of baryonic densities.

The bottom part of the left hand side of Figure 1 shows the time evolution of the $`K^+`$ production and its decomposition into different channels. the different production channels. We see that $`K^+`$ are produced quite early (full black line ) and that the maximum of their production co-incidents with the maximum of baryonic densities. The different channels show a rather similar shape. The production of $`K^+`$ in direct $`NN`$ collisions (dash-dotted line) acts earliest. It shows it maximum at about 4-6 fm/c, a time when the first violent reactions already took part and the transverse energy (left hand side, dotted line) reached its plateau. At that time first equilibration effects in the collision dynamics can be assumed. The production of $`K^+`$ in $`N\mathrm{\Delta }`$ (dashed line) and $`\mathrm{\Delta }\mathrm{\Delta }`$ collisions shows it maximum a little bit later. These collisions may dominantly take place at high densities in the stopped matter when some equilibrium is already reached. The $`\pi B`$ channel (full gr. line) finally peaks at still later time due to the delayed freeze-out of pions.

## 3 Analysis of kaon spectra

If the kaons are produced at highest densities we may assume that they should carry a signature of the KN potentials. Indeed the penalty of the KN-potentials can clearly be seen in the absolute yield of the produced $`K^+`$. This penalty shifts the threshold for the production toward higher values and therefore reduces the production yields. That effect can be easily seen on the right hand side of figure 1 where we compare laboratory momentum spectra taken at $`40^{}`$ laboratory angle taken from calculations with (open circles) and without (open triangles) KN optical potential to recent data taken by the KaoS collaboration (bullets) . For our calculations we assumed a total cross section of 5.3 barn for minimum bias reactions. We see that a calculation without KN potentials yields higher yields than a calculation with potential. The data favor the calculation with KN potentials. Nevertheless it should be noted that incertitudes on unknown production cross sections like $`N\mathrm{\Delta }NYK`$ may alter this statement .

In order to allow the analysis of high energetic kaons let us now turn to center-of-mass momentum spectra in full $`4\pi `$ as they are shown in figure 2. As we already stated, the penalty of the KN potentials shifts the threshold for the production toward higher values. Furthermore the produced kaons get less energy from the producing collision. However, this energy is regained in the expansion when the kaon leaves the nuclear medium. Therefore it is not very surprising to see on the left hand side of figure 2 that the spectrum of the $`K^+`$ in a calculation with potential (full line with circles) does not very much differ to that in a calculation without KN-potentials (dashed line with triangles), especially at high kaon momenta. Remember that the energy range shown in fig. 2 is larger than that shown on the right hand side of fig. 1. Only at low momenta we see a difference in the absolute values. This difference corresponds to kaons which are produced slightly above the threshold. In calculations without potentials they can be produced but they get little momentum for kinematic reasons. In calculations with potential their production is forbidden due to the enhanced threshold. For high kaon momenta re-scattering plays an important role. If we disable kaon re-scattering the values at high momenta drop down visibly.

The right hand side of figure 2 decomposes the spectrum obtained in a calculation with KN potentials into different time windows. The time-window $`t=`$0-4 fm/c (full line with circles) corresponds to the violent non-equilibrated phase where density and transverse energy is rising (a temperature cannot be defined at that time). We see a strong contribution at high cm momenta. When we go to later time windows the spectra are becoming steeper. The window $`t=48`$fm/c (dotted line with diamonds) corresponds to the time when equilibrium is starting but compression is still increasing. It shows the dominant contribution to the overall spectrum. In the expansion state (t=8-16fm/c, dash-dotted line with square) a strong contribution can be found at lower momenta. The very late kaons ($`t>16`$ fm/c, dashed line with triangles) finally do not contribute significantly. This effect is coherent with the observation that re-scattering is important for the high momentum part of the spectrum. The later the kaons are produced, the less time and less density they find for doing re-scattering. The behavior “the later the production the steeper the spectra” is already known from the analysis of pion spectra  and from perturbative studies of kaon production . It indicates that the detailed analysis of the spectra gives information on the time structure of the collision dynamics rather than on the KN potentials.

## 4 Azimuthal distributions at mid-rapidity

Let us now turn to the azimuthal distributions of kaons at mid-rapidity. For this observable an influence of medium effects has been reported by . Fig 3 shows the azimuthal distributions for calculations with KN potentials (full line with circles) and without KN potentials (dashed line with squares) for central ($`b<5.9`$fm, left hand side) and peripheral ($`b>7.8`$fm, right hand side) collisions. We see a strong dependence of the centrality. While central collisions only show a slight maximum at $`90^{}`$, i.e. out of the reaction plane, this maximum becomes a significant peak for peripheral collisions. This centrality dependence indicates the importance of the spectator matter for understanding this behavior. We have to discuss whether it is the potential or the re-scattering in matter which is important.

Concerning the potential we already see a large influence on the absolute yield of the distribution. In order to compare the shapes of the curves figure 4 shows normalized distributions for peripheral collisions. In central collisions the maximum of the distribution is less significant.

The left hand side of figure 4 shows the influence of different cuts which have been used by experiment. The dashed line with squares shows the distribution integrated over all rapidities, the dotted line with triangles shows the distribution for $`|y(cm)/y_{Proj}(cm)|<0.4`$ and the full line with circles shows the effect of an additional application of a cut in transverse momentum $`0.2`$GeV/c$`p_T0.8`$GeV/c. We see that each cut yields for slightly stronger peak of the squeeze. In the calculations of fig. 3 both cuts have been applied in order to allow a better comparison to experimental data of the KaoS collaboration .

The right hand side of figure 4 compares the azimuthal distribution at mid-rapidity (and with the given $`p_T`$-cut of $`0.2`$GeV/c$`p_T0.8`$GeV/c) obtained in a calculation with (full line with circles) and without (dashed line with squares) KN potentials, both using full re-scattering of the kaons. We see a visible increase of the peak when KN potentials are active. If we use a calculation with KN-potentials but disable re-scattering (dotted line with triangles) we see a slight decrease of the squeeze signal if compared to the calculation with full re-scattering (full line with circles). However, it should be kept in mind, that the influence of the cuts (left hand side) is of about the same size than the discussed effects.

If we now look for the sources, i.e. for the initial momenta of the kaons of the production the curve is flat or even slightly convex. A calculation without potential and re-scattering yields a similarly flat distribution. We conclude that the influence of the nuclear medium - as well potential as re-scattering - is important for explaining the squeeze, both effects contribute. Calculations of  see a stronger contribution of the potentials which might perhaps be due to implementations of Lorentz-type forces.

## 5 Conclusion

We have analyzed azimuthal distributions and energy spectra of kaons in Au+Au collisions at 1.48 AGeV. We find that the KN potential interactions only show a significant influence on the absolute yield of the kaons, which is due to the penalty to be paid when producing a kaon in a dense medium. Concerning the spectra only slight effects in the low momentum range can be reported. However, the dynamics of the collision reaction affects the spectra. Early reactions contribute stronger at high momenta while late reactions contribute more to the low momentum part. Thus, when decomposing the spectra into different time windows we find the spectra become steeper the later the time window is placed. Concerning the azimuthal distributions at mid-rapidity we see a visible influence of the KN potential and of the re-scattering. The influence of re-scattering is in coherence with the centrality dependence of this observable showing a maximum signal at peripheral collisions. In conclusion the shape of spectra and azimuthal distributions give us information about the collision dynamics rather than on a KN optical potential.

## References


text:
# Signatures of Ultra-High Energy Cosmic Ray Composition from Propagation of Nuclei in Intergalactic Photon Fields

## 1 Introduction

After a few decades of observation by several projects, the nature of the UHECR is gradually being revealed. Based on the observational results, there is no doubt of the existence of UHECR up to $`10^{20}`$ eV . This has led to consideration of a large number of production models. Most of the models assume that UHECRs are protons or photons.

Several studies of anisotropy have been done to reveal the origin of the UHECR . Around $`10^{18}`$ eV, a large scale anisotropy which is correlated to the Galactic center was reported by AGASA group . This result, if confirmed, supports an interpretation that cosmic rays below $`10^{18}`$ eV are dominated by a Galactic origin. Above this energy, however, no significant large scale anisotropy has been detected by AGASA .

On the other hand, AGASA reported clusters in the distribution of arrival directions including 4 doublets and 2 triplets with energy above $`4\times 10^{19}`$ eV. This results indicate that the UHECR come from extra-galactic point sources. The cluster events reported by AGASA are composed of particles with different energies. A possible explanation would be neutral particles (like neutrinos) propagating in straight lines from the source. Another possible explanation would be charged particles with large rigidity , i.e. with small charge like protons. If there is only the effect of the Galactic magnetic field then protons of different energy can arrive at same direction . However, these two possibilities cannot explained the three bumps in the energy spectrum of the cluster events between $`10^{18.8}`$ eV and $`10^{19.8}`$ eV  reported by AGASA. This structure is not statistically significant. However, it suggests to us the interesting possibility that UHECR might be nuclei as it will be discussed in this paper. In this scenario, cluster events would be formed by particles with different charge and energy but with the same rigidity. The atomic number of the primary cosmic ray does not affect the energy determination of AGASA which measures the number of charged particles on the ground. However, the capability of AGASA to measure the composition of the primary cosmic rays is not very strong . Nevertheless, the structure of the energy distribution of the cluster events may be used to measure of the composition indirectly as discussed in this paper.

If UHECR come from point sources and their rigidity is large enough or their source close enough so the direction is not lost, a fraction of the particles emitted by the source will arrive within a small angle from the source position. These particles are what we will call in this paper “small deflection angle” particles. Clusters will be formed when the number of particles expected within this angle in our experiment is larger than one. The small deflection angle particles should represent the cluster events observed by AGASA.

UHE protons should interact with Cosmic Microwave Background Radiation (CMBR) and lose energy by photo-pion production. For this reason, $`10^{20}`$ eV protons cannot propagate in intergalactic space beyond 50 Mpc, leading to a cut off in the spectrum below $`10^{20}`$ eV, the so called GZK cut off. On the other hand heavy nuclei will interact during its propagation with the Intergalactic Infrared Background Radiation(IIBR). Unfortunately the IIBR cannot be measured directly so far because of the strong background of infrared photons of Galactic origin. The main component of the IIBR is thought to be red-shifted stellar emission and reemission from interstellar dust. The UHE nucleus should be disintegrated by the IIBR and lose energy by the photo-disintegration interaction. For this reason, it was believed that the attenuation length of UHE nuclei is shorter than that of protons in intergalactic space . However, recent calculations which are based on empirical data show the IIBR density is much lower than originally thought . The lack of absorption in the energy spectrum of very high energy gamma rays from AGNs strongly supports this calculation. Using this result, the attenuation length of the UHE nuclei in the intergalactic photon field was re-estimated  and found to be longer than that of a proton (a few 100 Mpc for $`10^{20}`$ eV Fe ).

We simulate the propagation of UHE nuclei which are emitted from extra-galactic point sources based on a Monte Carlo approach. Using this simulation, we found structures in the energy distribution of the UHE nuclei which are observed as cluster events on the earth. In section 2 we briefly describe the interaction of nuclei with background photons. In section 3 we describe the scattering of the nuclei by the intergalactic magnetic field. In section 4 we show the result of this simulation. We conclude in section 5 and discuss the results briefly.

## 2 Interaction of Nuclei with background photons

There are three processes which affect cosmic rays during intergalactic propagation: energy loss by interaction with cosmic background radiation; deflection by magnetic fields; and energy loss by adiabatic expansion of the universe. The effect of the adiabatic expansion is negligibly small if propagation distances are less than 1 Gpc.

There are four interaction processes of cosmic rays with intergalactic background radiation. The first process, Compton interaction, is negligibly small for our purpose. The second process is pair-production of e+e-. This interaction can occur if the energy of the background photon is greater than 1 MeV in the rest system of the nucleus. If the Lorentz factor of the nucleus is larger than $`10^9`$, it will efficiently interact with the CMBR (Figure 1). In this case, the energy loss rate ($`dE/dt`$) is proportional to the square of the charge of the nucleus. The third process is photo production which mainly affects protons with energy greater than $`10^{20}`$ eV. The interaction of heavier nuclei with the CMBR by this process is negligible for the energy range considered in this paper.

The fourth process is photo-disintegration. The nucleus interacts with the background photon and disintegrates to a lighter nucleus. In this paper, we calculate this process by a Monte-Carlo method according to . These authors parametrized the total cross section $`\sigma (ϵ)`$ as a function of photon energy in the rest frame of the nucleus. Then the probability of photo-disintegration per unit length $`R`$ is calculated by following equation:

$$R=_0^{\mathrm{}}n(ϵ)\left[\frac{_0^{2\gamma ϵ}\sigma (ϵ^{})\frac{ϵ^{}}{ϵ}𝑑ϵ^{}}{2\gamma ϵ}\right]𝑑ϵ$$
(1)
where $`ϵ`$ and $`ϵ^{}`$ are the energies of the photon in the lab frame and in the nucleus rest frame respectively, $`n(ϵ)`$ is the differential number density of photons including the CMBR and the IIBR , and $`\gamma =ϵ^{}/ϵ`$ is the Lorentz factor of the nucleus. The term inside the bracket corresponds to the angle-averaged cross-section for a photon of energy $`ϵ`$. Figure 1 shows characteristic time $`\tau `$ for the energy loss due to photo-disintegration for a single-nucleon emission process as a function of energy. $`\tau `$ is defined by following equation:

$$\frac{1}{\tau }=\frac{1}{E}\frac{dE}{dt}$$
(2)
The effect of the CMBR is maximized when the Lorentz factor $`\gamma `$ is around $`10^{10}`$. In case of $`Fe`$, this value corresponds to an energy of $`10^{21}`$ eV, yielding a value of $`\tau =2\times 10^{14}`$ s (2 Mpc). Because the IIBR density is much lower than the CMBR density, $`\tau `$ increases rapidly at lower energy. $`\tau `$ of $`Fe`$ with $`10^{20}`$ eV is about $`2.5\times 10^{17}`$s (2.5Gpc). In general, one or a few nucleons and a single lighter nucleus are emitted by this interaction. The Lorentz factor is conserved at each photo-disintegration interaction, though it is reduced by the pair-production. Finally, UHE particles will pile up at a Lorentz factor around $`10^9`$. When $`{}_{}{}^{9}Be`$ is disintegrated, one proton and two $`He`$ are emitted. Therefore no nucleus is created with a mass number between 5 and 8.

The interaction probability is calculated for every possible photo-disintegration process involving the emission of one or more nucleons for all nuclei lighter than $`Fe`$. Based on these probabilities, the propagation of the nucleus is simulated by the Monte-Carlo method. The effects of pair production, the adiabatic expansion, and photo-pion production for protons are taken into account analytically in this simulation.

Figure 2 shows the life of an $`Fe`$ nucleus whose initial energy is $`4.9\times 10^{20}`$ eV. The $`Fe`$ is disintegrated during its propagation and turns into a lighter nucleus with the emission of nucleons. When the nucleus is $`{}_{}{}^{42}Ca`$, four protons and four neutrons are emitted. Then a $`Sulfur`$ is created. In comparison with energy, the Lorentz factor changes smoothly because photo-disintegration process does not change the Lorentz factor. This process makes a lighter nucleus and nucleons with same Lorentz factor.

## 3 Scattering by the intergalactic magnetic field

The effect of magnetic fields is described by the rigidity defined as the ratio of energy to charge ($`=E/Z`$). Particles with small rigidity are deflected by the magnetic field and cannot be observed as a cluster. The deflection of magnetic field increases the propagation time, and therefore particles emitted from distant sources may not reach to the observer. We calculate the effect of the intergalactic magnetic field based on a Monte Carlo method.

To simulate the scattering by the intergalactic magnetic field, we assume a Kolmogorov spectrum for the random magnetic field according to the reference of . In this reference, authors divide space in a lattice of 250 kpc cubes. The lattice is filled with a random magnetic field which follows the Kolmogorov spectrum with three logarithmic scales. Three field vectors of random orientation are sampled at scales $`l`$ = 1000, 500, and 250 kpc with amplitudes proportional to $`l^{1/3}`$. The final magnetic field in each 250 kpc cube is vectorial sum of these three vectors. The average magnitude of the magnetic field is assumed to be 1 nG. Particles propagate in spiral trajectories until they leave the lattice. Figure 3 shows examples of the trajectories. $`Fe`$ with energy of $`10^{18}`$ eV does not rapidly lose energy by the interaction with photons, and is trapped by the magnetic field inside 1 Mpc cubes. In case of $`2\times 10^{20}`$ eV, $`Fe`$ is disintegrated by the photons rapidly and the products propagates close to a straight line in the initial direction of the $`Fe`$.

## 4 Energy Distribution

In our simulation, the source energy spectrum $`dN(E)/dE`$ is assumed to have a power law dependence $`E^2`$ and have a cut-off at the energy of $`Z\times 2\times 10^{19}`$ eV, where $`Z`$ is charge of each primary nucleus. Below 300 MeV/nucleon, the cosmic nuclear composition can be divided to 4 types: $`He`$, $`CNO`$, $`NeSi`$, and $`Fe`$. Nuclear abundance between $`He`$ and $`C`$ is relatively low as well as between $`Si`$ and $`Fe`$. We treat the composition as 4 components, $`{}_{4}{}^{2}He`$, $`{}_{14}{}^{7}N`$, $`{}_{24}{}^{12}Mg`$, and $`{}_{56}{}^{26}Fe`$. The fraction of these nuclear abundances in the UHE region are unknown. Therefore we assume equal fractions relative to $`He`$ at the source. 20000 particles are created for each nucleus.

Figure 4 shows the variation of the energy distribution after propagation from a single source initiated by $`Fe`$ through the photon and magnetic field. Nuclei with Lorentz factor above $`10^9`$ interact with the CMBR rapidly. After 1 Mpc propagation a small pile-up is induced below the corresponding energy and the maximum energy of the particle is reduced. Since we assume a cut-off energy in the source spectrum, protons are emitted with energies below $`Z/A\times 2\times 10^{19}(10^{19})`$ eV where $`A`$ is mass number of the primary nucleus. The number of low energy nuclei ($`<10^{19}`$eV) increases up to 10 Mpc from the source since these nuclei are trapped by the magnetic field. These low energy nuclei cannot propagate farther than few 10 Mpc. Therefore, the energy distribution becomes steeper for distances below 10 Mpc and flatter for larger distances. No particles remain after propagation of a few 100 Mpc.

Figure 5 is same result as lower left panel of Figure 4. The energy distributions at a distance of 50 Mpc from the source are shown for different deflection angle (angle between the arrival direction and the source position). Particles with the deflection angle smaller than 2.5, 5, 10, 180 degrees are counted (180 degrees corresponds to the entire sample of particles included in Figure 4). Three bumps appear clearly in the distribution of particles with a small deflection angle. The bump below $`10^{19}`$ eV is composed of protons which are emitted from higher energy nuclei in the vicinity of the observer who is located at 50 Mpc distance from the source. The bump around $`10^{19.4}`$ eV is composed of $`He`$, and the bump above $`10^{19.6}`$ eV is composed of nuclei heavier than $`Be`$.

These bumps can be explained as follows: The rigidity is proportional to $`Z^1`$ and $`E`$, and $`Z`$ is approximately proportional to its mass number $`A`$ while the Lorentz factor is proportional to $`A^1`$ and $`E`$. This means that the Lorentz factor of a particle is approximately proportional to the rigidity. Therefore, if the Lorentz factor of the nucleus is large, the energy is reduced by the photo-disintegration. If the Lorentz factor is small, the particle is deflected by the magnetic field due to its small rigidity. However particles with Lorentz factor of about $`10^9`$ remain after applying the deflection angle cut. As a result, the bumps appear in the energy distribution of the small deflection angle particles at positions dependent on the mass numbers. The combination of the photodisintegration and the magnetic field acts as a “filter” for $`\gamma =10^9`$ particles. The energy distribution of the small deflection angle particles shows clearly the spectral features created by interactions with the intergalactic photon background and magnetic fields during their propagation.

If the sources are distributed uniformly in the Universe, the expected energy distribution measured at the earth can be estimated by adding up linearly the energy distributions corresponding to different distances from the source. In Figure 4 the number of particles around $`10^{20}`$ eV coming from sources at 100 Mpc is 1.5 orders of magnitude smaller than from sources at 100 kpc. When adding linearly the different distances, the number of particles from sources at 100 Mpc is enhanced 3 orders of magnitude compared to sources at 100 kpc. Therefore, the contribution to the energy spectrum at earth is increasing with the distance to the source till a couple of hundred Mpc. On the other hand, most of the low energy particles emitted by distant sources will be trapped by the magnetic field and they will not reach the earth. However, the low energy secondaries resulting from the photodisintegration of high energy heavy nuclei close to earth will regenerate the low energy spectrum coming from distant sources. Therefore, even at low energies we will still be dominated by the contribution of distant sources. The spectrum of small deflection angle particles will follow the same behavior, i.e. the contribution of sources at large distances to the spectrum exceed largely the one from small distant sources, even if the fraction of small deflected particles for closer sources is slightly larger.

It should be noted that the relation between cluster events and small deflection particles is qualitatively clear but difficult to estimate quantitatively. Current work is in progress for a better comparison with experimental data.

Figure 6 shows the result of the expected energy distribution under the assumption that the sources are distributed uniformly in the Universe. The solid line shows the expected energy distribution of all the particles that reach the earth. The expected energy spectrum is steeper than the source spectrum. Above $`10^{20.2}`$ eV, the energy spectrum becomes even steeper and the maximum energy is about $`10^{20.4}`$ eV. The exact details of the upper end of the spectrum will depend on the source spectrum cut off.

The dashed lines in Figure 6 show the expected energy distribution of the small deflection angle particles. This spectrum is much flatter than the source spectrum. Four bumps appear in this figure. The bump just below $`10^{19}`$ eV corresponds to protons and strongly depends on the cut off energy at the source. The bump around $`10^{19.3}`$ eV is composed of $`Deuteron`$ and $`He`$. Every nucleus with a Lorentz factor around $`10^9`$ contributes to the $`He`$ bump after propagation. The bump between $`10^{19.6}`$ eV and $`10^{20}`$ eV correspond to nuclei with mass between $`Be`$ and $`Mg`$. The bump above $`10^{20.1}`$ eV is composed of $`Fe`$ and other nuclei of similar atomic number.

This total spectral structure bears a striking resemblance to the result of AGASA .

## 5 Results and Discussion

We have simulated propagation of UHE nuclei which are emitted from extra-galactic point sources using a Monte Carlo approach. Using this simulation, structures in the energy distribution of the small deflection angle particles have been predicted. In this simulation, we assume that UHECR are nuclei which are emitted from extra-galactic point sources. The nuclei interact with cosmic background radiation and lose energy during their propagation. Based on recent estimations of the IIBR, it has been shown that UHE nucleus can propagate longer distances than a proton with same energy. It is also expected that heavier nuclei are accelerated more efficiently to high energy than lighter nuclei or protons.

We have considered the effect of the intergalactic magnetic fields. We use a Kolmogorov spectrum for random magnetic fields that is 1 nG in average intensity with a 1 Mpc correlation length. The energy distribution of the particles with deflection angles smaller than 10, 5, and 2.5 degrees are investigated. Consequently, for this model no cosmic rays can be observed after propagation of a few hundred Mpc. Lower rigidity particles are deflected by the magnetic field relatively rapidly. Particles with larger Lorentz factor lose energy by photodisintegration. Since the rigidity and the Lorentz factor are proportional to the particle energy and inversely proportional to the mass number, particles which have common a Lorentz factor, or similar rigidity, remain after the deflection angle cut. The $`BeMg`$ bump appears at an energy more than 4 times larger than the $`He`$ bump, and the $`Fe`$ bump appear around 6 times larger energy. The rigidity gives a lower bound and energy loss gives an upper bound on the bumps. The energy distribution of the clustered events observed by AGASA can be explained by this model.

In the simulation, we assume a specific intergalactic magnetic field distribution, cut-off energy and source composition of primary nuclei. The Galactic magnetic field is ignored in this simulation. The events with certain deflection angle are selected and shows the effects due to the intergalactic photons and magnetic fields. This effect depends on the strength and correlation of the intergalactic magnetic field. This dependence gives an ambiguity to the analysis since the deflection angle cut is adjustable. Therefore, this angle should be treated carefully.

This result suggests that if we detect few hundred events above $`10^{20}`$ eV in future experiments, we should be able to extract the structure of the bumps in the energy distribution of cluster events. The Pierre Auger Observatory is quite sufficient for this purpose and may confirm the signature of nuclei in the spectrum in the near future.

## Acknowledgments

We thank Angela Olinto, Alan Watson, James W. Cronin, Maximo Ave, Aaron Chou, and Medina Tanco for useful discussions. This research was supported by grant NSF PHY-0114422 to the Center for Cosmological Physics, University of Chicago, and grant NSF PHY0103717.


text:
# Systematic Study of Electron Localization in an Amorphous Semiconductor

## I INTRODUCTION

Amorphous semiconductors represent a large and important area in materials science, with interest both from the technological and fundamental point of view. Coordination defects in *a*-Si are of key importance to bulk and transport properties. Understanding the nature of defects in *a*-Si can help unlock the mystery behind phenomena like the Staebler-Wronski effect Elliot  and help establish the link between localization of defect states and large electron-phonon coupling. Early theoretical work on defect states in *a*-Si and *a*-Si:H was based on tight-binding methods Biswas ; Allan ; Fedders ; Holender ; Kneif ; Min . For example, Biswas *et al.* Biswas  and Fedders and Carlsson Fedders  investigated the electronic structure of dangling and floating bonds in *a*-Si. They showed that gap defect states associated with dangling bonds are strongly localized on the central atom of the dangling bond Biswas . More recently, density functional calculations of dangling bond states using the local density approximation (LDA) Kohn1  have been performed by Fedders and Drabold Drabold1 . They reported a wave function localization of 10-15% on the central atom in supercell models with one defect and far less on supercell models with many defects due to defect band formation. This finding was at variance with ESR experiments, which showed that over 50% of spin density of the gap state is located on the central atom of the dangling bond Biegelsen ; Umeda . However, recent calculations by Fedders *et al.*Drabold2  using the local spin density approximation (LSDA) have shown that a large spin localization of a defect state does not necessarily imply the existence of a corresponding wave function localization.

Density functional theory (DFT) Kohn2  has enjoyed enormous success in describing the ground state properties and defects for a wide range of materials, and in particular, *a*-Si. Nonetheless, the one-particle Kohn-Sham energies in the theory have no formal justification as quasiparticle energies. However from an empirical point of view, Hybertsen and Louie have shown, using GW calculations, Louie1  that for states close to the fundamental band edges of bulk semiconductors and in particular Si, there is a 99.9% overlap of the quasiparticle wave function with the corresponding Kohn-Sham orbital (GW calculations provide post Hartree self-energy corrections to DFT/LDA). This provides some rationale for interpreting the Kohn-Sham orbitals as quasiparticle states.

Within the density functional framework, there is a general problem for the accurate representation of localized midgap and band tail states in amorphous semiconductors. The reliability of density functionals to correctly reveal the localized behavior of electronic states with respect to its wave function and spin is very important. In particular, the generalized gradient approximation (GGA) and LSDA results for electron localization are sometimes quite different. For example, recent density functional and GW studies of the metal-insulator transition of bcc hydrogen showed that eigenfunctions obtained GGA are more localized Louie2 ; Louie3  and closer to quasiparticle energies and states than LSDA. Also, it was observed that GGA bandgap was systematically larger than LSDA gapLouie2 .

In this paper, we systematically examine the localization of band tail and gap states and its dependence on basis sets, density functionals and spin polarization for three defected models: two 216-atom supercell models of amorphous Si and a 218-atom supercell model of crystalline Si:H with a vacancy. The crystalline model will serve as a benchmark to compare the nature of a dangling bond defect state in an amorphous environment with that of the crystalline phase. We compare the localization of gap and tail states within the LDA, LSDA and GGA for frozen lattices (unrelaxed samples) as well as for samples fully relaxed for a given Hamiltonian. Our motivation for performing the frozen lattice calculation is to systematically investigate the sole effect of different basis sets and density functionals on the localization of states for a fixed configuration. We study the relaxation effects to see the dependence of the local geometry of the defect sites on the different basis orbitals and density functionals and the behavior of localized defect states in a relaxed environment compared to the frozen one. We computed spin and wave function localization for defect states, to determine the correlation between spin density and charge density.

## II METHODOLOGY

### II.1 Models and calculations

Two of the models used here were generated by Barkema and Mousseau Moose  using an improved version of the Wooten, Winer and Weaire (WWW) algorithm Wooten . The details of the construction was reported in Ref. 19 and we just give a brief recap here. Two independent 216-atom models of *a*-Si, each with two dangling bonds, were generated. Atoms are first packed randomly in a box, at crystalline density, with the single constraint that no two atoms be closer than 2.3 Å. A connectivity table was then set up by constructing a loop that passes exactly twice through each atom. To create a dangling bond, a ghost bond (bond with zero energy) is placed between two atoms. The two atoms are chosen in such a way that they are quite close in one model and reasonably far in another model. The network is then relaxed using a fixed list of neighbors and a Keating potential. These steps ensure that the initial state has no trace of crystallinity. The models are then relaxed using a series of WWW moves, using the accelerated algorithm discussed in Ref. 5. The conformations are finally relaxed at zero pressure with a Keating potential. In this work we refer these two 216-atom supercell models of *a*-Si as CLOSE and FAR. In model CLOSE the atoms with the dangling bonds are separated by a distance of 4.6 Å, whereas in the model FAR they are 7.7 Å apart. The third model, which we will refer to as *c*-Si:H, is generated as follows: starting with a 216-atom cell of silicon in the diamond structure, an atom is removed resulting in the creation of a vacancy with four dangling bonds. Three of them were terminated by placing an H atom at about 1.5 Å from the corresponding threefold-coordinated Si atom. This results in one isolated dangling bond. Throughout this work, we will use a cut-off radius $`\mathrm{R}_{\mathrm{Si}\mathrm{Si}}`$ = 2.6 Å to define Si-Si coordination.

Our DFT calculations were performed within the LDA (with and without spin polarization) and the GGA using the code SiestaSiesta1 ; Siesta2 ; Siesta3 . We used the parameterization of Perdew and Zunger Perdew1  for the exchange-correlation functional in all LDA calculations and that of Perdew, Burke and Ernzerhof Perdew2  for the exchange-correlation functional in all GGA calculations. Norm conserving Troullier-Martins pseudopotentials Troullier  factorized in the Kleinman-Bylander form Kleinman  were used to remove core electrons. To describe the valence electrons, we used atomic orbitals basis set consisting of finite-range numerical pseudoatomic wave functions of the form proposed by Sankey and Niklewski Sankey . We employ single-$`\zeta `$ (SZ) and double-$`\zeta `$ with polarization functions (DZP) basis sets on all the atomsbasis . We solved the self consistent Kohn-Sham equations by direct diagonalization of the Hamiltonian. The $`\mathrm{\Gamma }`$ point was used to sample the Brillouin zone in all calculations.

### II.2 Quantifying the degree of localization

To characterize the spatial extent of an electronic state we use three measures of localization. Two of them depend on the Mulliken (point) charge Szabo  $`\text{q}_i(\mathrm{E})`$ residing at an atomic site $`i`$ for an eigenstate with energy eigenvalue $`\mathrm{E}`$. The charges are normalized ($`_{i=1}^\mathrm{N}\text{q}_i(\mathrm{E})=1`$, where $`\mathrm{N}`$ is the total number of atoms in the supercell). The third is a measure of the spread of the wave function in real space, as the the second moment of the radial displacement about its center (variance). We briefly describe below each measure and its interpretation.

The first quantity, the conventional inverse participation ratio (IPR)

$$\mathrm{I}(\mathrm{E})=\underset{\mathrm{i}=1}{\overset{\mathrm{N}}{}}(\text{q}_\mathrm{i}(\mathrm{E}))^2,$$
(1)
is a measure of the inverse of the number of sites involved in the state with energy $`\mathrm{E}`$. For a uniformly extended state, the Mulliken charge contribution per site is uniform and equals $`\frac{1}{\mathrm{N}}`$ and so $`\mathrm{I}(\mathrm{E})=\frac{1}{\mathrm{N}}`$. For an ideally localized state, only one atomic site contributes all the charge and so $`\mathrm{I}(\mathrm{E})=1`$. This implies that high values of $`\mathrm{I}(\mathrm{E})`$ correspond to localized states and low values of correspond to delocalized states. Despite its ubiquity, the IPR is an *ad hoc* measure of localization. The second quantity has been used by Lewis *et. al* Sankey2 ; Sankey3  to measure the localization of DNA electron states. It utilizes the information entropy Shannon ; Jaynes  $`S(\mathrm{E})`$ of an eigenstate with energy $`\mathrm{E}`$ and is given by

$$\mathrm{W}(\mathrm{E})=\mathrm{exp}(\mathrm{S}(\mathrm{E}))$$
(2)
where $`S(\mathrm{E})=_{\mathrm{i}=1}^\mathrm{N}\text{q}_\mathrm{i}(\mathrm{E})\mathrm{ln}\text{q}_\mathrm{i}(\mathrm{E})`$. In information theory, the entropy is a measure of disorder or randomness in the {$`\text{q}_i`$}. As such, it has more *a priori* justification for use as a measure of localization. This is especially true if we interpret localization as a departure from randomnessed . We interpret $`\mathrm{W}(\mathrm{E})`$ as an estimate of the *number of accessible atoms* for a given state or equivalently, the number of atoms a particular electronic state can reach. It follows that $`\mathrm{W}(\mathrm{E})`$ ranges from the total number of atoms for a uniformly extended state (that is, all the atoms are accessible) to unity for an ideally localized state (that is, only one atom is accessible). For purposes of discussions and comparisons to IPR, we use $`\frac{1}{\mathrm{W}(\mathrm{E})}`$ for our work.

The third quantity we use to determine the localization of the wave functions is a measure of the spread in real space, as defined in the context of localized wave functions in molecules boys  and maximally localized Wannier functions in solids. mlwf ; wannier  This definition assumes that the wave function has a well defined center, and that the average of the position operator can be defined. This is not the case for our eigenstates, which correspond to the $`\mathrm{\Gamma }`$ point and therefore are periodic in the supercell. However, since we are interested in states which are localized (exponentially) within the supercell, we can still use the definition of the spread, by integrating over the cubic supercell volume centered at the defect site:

$$\sigma ^2=[𝐫^2𝐫^2].$$
(3)
We also compute the amount of the norm of the wave function that is around the defect center, by integrating over a sphere $`\mathrm{\Omega }`$ centered at the defect and which has the same diameter as the lattice constant (so it is the largest sphere that fits into the supercell volume). Clearly, the defect states will not be completely localized within the sphere, and the norm of the wave function inside the integration volume, $`q_\mathrm{\Omega }`$, will be smaller than one. This is also a measure of the amount of wave function which is outside the integration sphere: for an extended state, the amount of charge will be close to one half, since the sphere encloses approximately half the supercell volume. For localized states, we typically capture more than 80% of the wave function norm. We will characterize the localization by using both the spread $`\sigma ^2`$ (larger values correspond to less localization) and by the norm $`q_\mathrm{\Omega }`$.

## III Results

### III.1 Frozen lattice calculations

#### III.1.1 Wave function localization

In order to understand the effect of basis set and density functional on the localization of defects states, we first compute the electronic structure of the three models with fixed geometry, using the different basis sets and density functionals. We begin our discussion with simplest of the three models, c-Si:H. The position and localization for individual states are reported in Figs. 1 and 2, where the IPR and $`\frac{1}{\mathrm{W}(\mathrm{E})}`$, respectively, are used to measure localization. Each spike indicates an energy eigenvalue. With both measures of localization, and for all basis sets and density functionals considered, we observe a highly localized state near at Fermi level, with all the other states being extended. This state is centered on the dangling bond atom, which contributes with a Mulliken charge of $`0.54e`$ and $`0.57e`$ in the SZ-LDA and SZ-GGA cases respectively, and $`0.27e`$ and $`0.31e`$ in the DZP-LDA and DZP-GGA cases respectively. The rest of the state is mainly localized in the neighbor atoms to the dangling bond. Both the IPR and $`\frac{1}{\mathrm{W}(\mathrm{E})}`$ provide the same qualitative picture of the localized state, and its evolution with basis set and density functional. In particular, localization decreases strongly when we go from the SZ basis to a more complete DZP basis. It also increases but only by a small amount when using GGA instead of LDA. These are general trends which we will also observe with the other models, as we will see below. We note that the decrease in the measure of localization is not proportional to the decrease in the Mulliken charges at the defect state, as both definitions of localization are nonlinear. While the charge at the dangling bond site is reduced by a factor of two when moving from SZ to DZP, the localization measure decreases by roughly a factor of three, both for IPR and $`\frac{1}{\mathrm{W}(\mathrm{E})}`$.

Unlike $`\mathrm{I}(\mathrm{E})`$ and $`\mathrm{W}(\mathrm{E})`$ that are point estimators of the localization (in the sense that they use only the Mulliken charge at each atomic site), $`\sigma ^2`$ is a more physical representation. The results obtained for $`\sigma ^2`$ for *c*-Si:H with the lattice frozen are presented in Table 1. The spread for the localized state $`\mathrm{M}`$ is simple to compute since it is unimodal (peaked at only one site and therefore having a well defined center). From $`\sigma ^2`$, we also see monotonic increase in the spread and a decrease in the total charge in the localization volume as the basis sets are increased (see Table 1). The GGA states show a slightly smaller spread than the corresponding LDA states. In Fig. 3 we show snapshots of the isosurface of the wave function for the localized midgap state of *c*-Si:H within two approximations. We see a dangling bond orbital confined to a small region in space in the SZ case implying that the state is well localized. In the DZP case, we see pieces of the surface in the vacancy and other neighboring atoms besides the dangling bond orbital making it less localized compared to the SZ case.

Next, we analyze the localized nature of the states for the amorphous model CLOSE. In this system, we expect to see two localized states in the gap, corresponding to the two dangling bonds present in the structure. Indeed, we see two highly localized gap states H and L in the IPR shown in Fig. 4, with an energy splitting which is just over a tenth of an eV. The state H is the highest occupied molecular orbital (HOMO) and L is the lowest unoccupied molecular orbital (LUMO). Due to the small distance between the two dangling bonds, the localized H and L states are bonding and antibonding combinations of the dangling bond states, and therefore both H and L have nearly equal weights on the two defect sites. In the SZ case, the total Mulliken charge contributions from the two dangling bond sites for each of the two localized wave functions range between $`0.52e`$ and $`0.64e`$. The charge concentrations drop to the range $`0.29e`$-$`0.40e`$ in the DZP case. In Fig. 5, we show the results for $`\frac{1}{\mathrm{W}(\mathrm{E})}`$. We see a sharp drop in the number of atoms a given eigenstate can reach as the energy changes from the edge of either the valence or conduction band into the gap. We again see that all the features of $`\mathrm{I}(\mathrm{E})`$ are reproduced in $`\frac{1}{\mathrm{W}(\mathrm{E})}`$. For each localized state, both measures decrease by approximately a factor of 2 when the basis sets are increased from SZ to DZP.

In Fig. 6, we show the IPR for the model FAR. As in the CLOSE case, both the HOMO and LUMO states are localized. Now, however, since the distance between the two dangling bond sites is larger, the HOMO-LUMO splitting is much smaller, only $`6`$ meV. The HOMO has now most of its weight on one of the dangling bonds, whereas the LUMO is mostly localized in the other. The trend in localization of the gap states is similar to the the trend observed in the other two two models, decreasing strongly with more complete basis sets.

Since the wave functions associated with the gap states in the unrelaxed CLOSE and FAR models do not have a single center, but are peaked at the two dangling bonds, we will not analyze the localization by means of the spread in these cases. However, as we show in next section, relaxation leads in some cases to localization of the wave functions around one of the dangling bonds, and this will allow us to use the spread in such cases to quantify localization.

Our frozen lattice calculations consistently show that increasing the basis set decreases the localization. Although this is not unexpected, the huge decrease in the localization of the wave function as the size of the basis functions increase from SZ to DZP is quite surprising. The fact that the results are consistent in both the amorphous and crystalline system makes it even more interesting and general. A plausible explanation for this effect is that the energy gap is clearly reduced as the basis set is more complete. The localized states are then closer to the band edges, and therefore are more able to mix with the extended bulk states, becoming more delocalized. Obviously, the delocalization will not proceed indefinitely upon improvement of the basis set, but will converge as the basis approaches completeness.

#### III.1.2 Spin localization versus wave function localization

Fedders and co-workers Drabold2  have shown that, in order to correlate the degree of localization from dangling bond states with ESR experiments, it is not enough to look at the wave functions, but to the net spin polarization near the dangling bond. The reason is that the spin density has also contributions from electronic states other than the localized defect wave function, which contribute to make the spin polarization more localized than the specific localized state wave function. In order to confirm this result (obtained by Fedders et al. on cells of a-Si:H) in our structural models, we performed calculations allowing for spin polarization in our frozen lattice models, using the DZP basis set. Except for the $`c`$-Si:H case, where there is one unpair electron that yields a net spin polarization, we were not able to find a spin polarized solution for any of the amorphous cells. The reason is the existence of two interacting dangling bonds, which favors the formation of a spin singlet with two electrons paired. In order to force the appearance of a spin moment in our models, we introduce an unpaired spin by removing a single electron from the system.

In the model *c*-Si:H with, we find a contribution of almost 50% to net spin by the central dangling bond and its neighbors (the central atom alone contributing 38%). However, the Mulliken charge contribution to the wave function of the corresponding localized state from the defect site is only $`0.29e`$. The hydrogen-terminated dangling bond sites also contribute about 10% of the net spin. The remainder is somewhat distributed uniformly at the other sites. In CLOSE, about 57% of the net spin polarization was due the one dangling bond and its three neighbors. The other dangling bond contributed only 6% to the net spin with essentially zero contribution coming from the neighbors. In FAR, about 54% of the net spin localization sits on the isolated dangling bond and its nearest neighbors. Our results are in good agreement with the previous LSDA calculations by Fedders et al.Drabold2 , and in reasonable agreement with the experiment Biegelsen ; Umeda .

Our results confirm that, for a dangling bond defect state, there is a rather large difference between spin localization and wave function localization. In particular, the degree of spin localization is greater than that of the wave function localization at the dangling bond site. To our knowledge, no experimental methods exist for measuring the extent wave function localization on the dangling bond orbital as opposed to spin

### III.2 Relaxation effects

#### III.2.1 Geometry of defect sites, density functionals and basis sets

In this section we discuss the geometry around the defect sites before and after relaxation. The details of the local geometry are very important in determining the local electronic structure and the strain around the defect site, and here we study the dependence with varying basis and density functional. We relax all the models using a conjugate gradient optimization until the largest atomic force is smaller than 0.04 $`\mathrm{eV}`$/Å.

Relaxation effects for the simple dangling bond defect in *c*-Si:H are relatively small. There is no major rearrangements in bonding, but only a relaxation of the surroundings of the vacancy site.

In the unrelaxed CLOSE model, the dangling bonds were originally separated by a distance of 4.6 Å. After SZ basis relaxations, both with LDA and GGA, the defect sites come closer, to a distance of about 2.6 Å, to form a highly strained bond. For the more complete DZP basis set, the two defect sites also approach each other, but they continue being undercoordinated, so the two distinct dangling bonds remain present.

In the case of the FAR model, one of the two well separated dangling bonds forms a strained bond with a neighboring atom after a SZ-LDA relaxation. The dangling bond is therefore terminated, and a floating bond is introduced. The other dangling bond remains present. For the SZ-GGA relaxation, the two dangling bonds still continue to exist, but the one that was terminated in the SZ-LDA case also approaches a neighbor and tries to form a bond. After the DZP relaxations, the dangling bonds still exist both in LDA and GGA.

Our results indicate that SZ basis tends to favor tetrahedral bond formation whereas DZP qualitatively preserves the original structure with the dangling bonds present. Also, the SZ tends to favor the transformation of dangling bonds into floating bonds. The results therefore suggest that the richer DZP basis set is necessary for an accurate description of the geometry of both isolated and clustered dangling bonds in amorphous silicon. The SZ basis is not flexible enough to provide sufficient freedom to describe the different shape of the wave function at the dangling bonds compared to covalent $`sp^3`$ bonds (for which the basis is ideally suited). Therefore, it tends to favor the disappearance of the dangling bonds through annihilation with other dangling bonds or formation of floating bonds with other already fourfold coordinated atoms.

#### III.2.2 Localization, density functionals and basis sets

We now consider in detailed the trends in the localization behavior of the electronic states in the gap after full structural relaxations. We also examine the role the defect site in a relaxed environment plays in the localization of gap states. We study the localization using the IPR for each of the three fully relaxed models within the different approximations.

We first studied the simple dangling bond defect in *c*-Si:H. Fig. 7 shows that, for the relaxed structures, the localization behavior of the midgap state is density functional dependent but rather basis set independent, contrary to the results obtained in the frozen lattice calculations (Fig. 1). In other words, within the same density functional approximation, SZ yields a similar wave function localization as DZP for the simple defect in the relaxed crystalline environment. We also see that the GGA defect state is more localized than the LDA defect state for a given basis set. The IPR values obtained with DZP are, nevertheless, almost unchanged upon relaxation, the difference between the unrelaxed and relaxed geometries occurring mainly for the SZ basis. The analysis from the real space spread, shown in Table 2, confirms the results obtained via the IPR. We see again the similarity between SZ and DZP for a given functional, with the GGA states having a smaller spread than LDA.

The IPR for the fully relaxed CLOSE is shown in Fig. 8. As discussed in the previous subsection, structural relaxation for this model is basis dependent. We first see that the splitting between the HOMO and LUMO levels is now much larger than in the unrelaxes case. This is attributed to the fact that, in order to minimize the total energy, occupied defect states move down in energy and closer into the valence band edge, whereas the unoccupied states do not affect the energy and thus can move towards the conduction band edge. The second observation is that now the localization of both HOMO and LUMO has decreased considerably compared to the unrelaxed case. Again, this is a consequence of the levels being closer to the band edges, mixing more strongly with the delocalized bulk states and therefore becoming less localized. The effect is larger for the HOMO, which is the one that adjusts its shape to optimize the total energy. A third observation is the appearance of increasingly localized states in the band edges, corresponding to bulk states which start becoming localized and form the precursor of band tails. This effect is originated from the strain field imposed by the relaxation of the sites around the defects. Therefore, the presence of defects like dangling bonds in amorphous silicon also brings the appearance of band tails of weakly localized states due to the introduction of stress in the surroundings of the defect. This supports results from previous work Drabold3  that there is not a one to one correspondence between spectral and geometric defects. The localization of the tail states decays as one moves deeper into the conduction and valence regions, as was previously observed by Dong and Drabold Dong  using a simple orthogonal tight-binding Hamiltonian on a much larger 4096-atom model of *a*-Si.

We now focus on the evolution of localization with basis set and density functional for the relaxed CLOSE model. As we observed with the $`c`$-Si:H case, the HOMO level becomes less localized upon relaxation, specially in the case of the SZ basis, for which the IPR is reduced by more than a factor of two. The degree of HOMO localization predicted by the SZ and DZP bases for the relaxed structure is therefore very similar. For the LUMO, the difference between SZ and DZP is still large, as in the unrelaxed case. The high IPR values associated with the LUMO in the SZ cases are primarily due to strain (as a result of the bond between the two neighboring dangling bonds atoms). The spread of the LUMO is reported in Table 3. We again observe the common trends in the spread: GGA states show slightly less spread compared to LDA, and SZ basis set yield more localized states than DZP basis set.

In order to get a pictorial representation of the localized states in this relaxed model and the evolution with basis set and density functional, we assign different colors to each site according to its Mulliken charge contribution to the given eigenstate. We depict this spatial feature by showing only 65% of the total charge for the LUMO in Fig. 9. We observe a small network connection of atoms for the localized states in Figs. 9(a) and 9(b) but the connectivity spreads out in a rather 1D fashion, mimicking a chain in Figs. 9(c) and 9(d). The small size of our cell does not allow us to immediately visualize a localized region containing completely the cluster of atoms, which can be done in larger supercell containing thousands of atoms Dong .

We now turn to the relaxed FAR model. In Fig. 10, we have plotted the IPR for the this model, for different density functionals and basis sets. Again, we see the same general features that we pointed out in the relaxed CLOSE case. First, the splitting between HOMO and LUMO is much larger than in the unrelaxed case. Both states, and specially the HOMO, become more delocalized upon relaxation, with the notable exception of the SZ-GGA case, which yields a strongly localized HOMO state. We also see the formation of band tail states, and even the presence of strongly localized states in the gap above the LUMO, due to strong relaxation induced strain fields. Again, in this model, we observe that localization is stronger for GGA than with LDA, and that the difference in localization between SZ and DZP bases is much reduced upon relaxation. The difference in localization for the SZ-GGA and SZ-LDA cases can be explained by the fact that, as mentioned previously, SZ-LDA relaxation results in the disappearance of a dangling bond and the formation of a floating bond, which are known to be less localized than dangling bond defects.

Finally, we visualize a chosen state (HOMO in this case) for the relaxed FAR model using a color coding in Fig. 11 The Mulliken charge concentrations on the atoms changes from confined cluster-like character (or equivalently short 1D strings) in Fig. 11(a) to long string-like character of atoms in Fig. 11(d) as one tunes the basis and functional from SZ-GGA through to DZP-LDA . The ”tinker-toy” character can be attributed to less localized states and it is mainly due to weak quantum mechanical mixing. This behavior has been observed earlier by Drabold et. alDong ; Nak

## IV CONCLUSION

We have performed a first principles electronic structure calculation on three Si supercells: two 216-atom supercells of amorphous silicon with two dangling bonds and one 218-atom supercell of hydrogenated crystalline silicon with a void. Depending on the initial distance between the dangling bonds, the two a-Si models have been classified as CLOSE and FAR. We examined the nature of localized band tail and gaps states within the LDA and GGA using both minimal SZ basis and more complete DZP basis with particular attention paid to relaxation effects. Spin localization and wave function localization for dangling bond defect states states has also been studied. We computed the wave function localization as the spread of the wave function in real space and via other measures that utilize the Mulliken charges.

For the frozen lattice calculations, we find that the localization of wave functions associated with defect states decrease with larger basis sets and has enhancement of localization using GGA compared to LDA for all the models. The reduction in charges at the atomic sites for a larger basis set can be attributed to the hybridization between the atomic orbitals, providing the electronic charges more degrees of freedom to redistribute themselves. This is reflected in a smaller distance between the defect states and the band gap edges, that also favors delocalization. Unpaired spin LSDA calculations performed on frozen lattices showed that the degrees of spin and wave function localization are different. In particular, degree of spin localization at a dangling bond site is far greater than the degree of wave function localization. The difference between the localization of a defect state in a fully relaxed and frozen systems is non-trivial, especially in the minimal basis calculations. In particular, there was a considerable reduction in localization (as measured using the Mulliken charge) for the relaxed systems compared to the frozen lattices. We also conclude that a large basis set (DZP in this case) is necessary for an accurate description of both the geometry and localized states associated with defect sites.

###### Acknowledgements.

We are grateful to Professor Normand Mousseau for sending us his models of amorphous silicon. We also acknowledge the support of the National Science Foundation under Grant Nos. DMR-0074624, DMR-0310933 and DMR-0205858. P.O. acknowledges support for his research visit to Ohio University from the Programa de Movilidad de Investigadores of Ministerio de Educación y Cultura of Spain.