text:
# Very large magnetoresistance in lateral ferromagnetic (Ga,Mn)As wires with nanoconstrictions

## Abstract

We have fabricated (Ga,Mn)As nanostructures in which domain walls can be pinned by sub-10 nm constrictions. Controlled by shape anisotropy, we can switch the regions on either side of the constriction to either parallel or antiparallel magnetization. All samples exhibit a positive magnetoresistance, consistent with domain-wall trapping. For metallic samples we find a magnetoresistance up to 8%, which can be understood from spin accumulation. In samples where, due to depletion at the constriction, a tunnel barrier is formed, we observe a magnetoresistance of up to 2000 %.

Spin-valve effects involving magnetic semiconductors are of considerable importance for applications in spintronics. Because the carriers in these materials are holes (with strong spin-orbit coupling), it has been difficult to observe spin dependent scattering in the diffusive regime (GMR, or giant magnetoresistance)akiba . Tunnel magnetoresistance (TMR) experiments have been more successful, yielding spin-valve effects of several tens of percenttmr .

It was recently pointed outflatte  that large MR effects can be expected from domain walls (DW) in magnetic semiconductors due to the large spin polarization in these materials. Domain wall resistance (DWR) in ferromagnetic metals, where small effects are typical, has been intensively investigatedkentreview . Large MR’s, however, have been observed in mechanically manipulated nanojunctionsgarcia , though magnetostriction can play a roleansermet .

Here we address the suggestion of Refs. flatte  in an experimental study of DWR in the ferromagnetic semiconductor (Ga,Mn)As. We use lateral nano-fabricated constrictions in single domain wires to pin the DW and reduce their lengthbruno . This approach facilitates ballistic hole transport through the DW, while the lateral fabrication technology excludes any influence of magnetostriction effectsansermet . We find that DWR leads to very large spin-valve effects in both the GMR-(up to 8 %) and TMR-(2000 %) regimes.

Our samples were made from a thin (19nm) epitaxial layer of Ga0.976Mn0.024As, grown on a semi-insulating GaAs (001) substrate by low temperature molecular beam epitaxy. The carrier density from etch capacitance-voltage calibrations is about $`3\times 10^{20}`$ cm-3, and the sheet resistivity at 4.2 K is about 4.5 k$`\mathrm{\Omega }/\mathrm{}`$. Assuming an effective hole mass as in GaAs, $`m^{}0.5m_o`$, where $`m_o`$ is the free electron mass, these values imply a Fermi energy $`E_F`$ of 150 meV, a Fermi wavelength $`\lambda _F`$ of 6 nm, and a transport mean free path $`l_t`$ of ca. 1 nm. The Curie temperature of the material is 65 KOhnoBook .

We have fabricated transport structures consisting of a central island of 100 nm width and 500 nm length connected to two 400 nm wide and 10 $`\mu `$m long wires by constrictions with widths down to 10 nm or less (Fig.1). The constrictions act as pinning centers for DW. The 400 nm wide wires are contacted by voltage- and current-leads, allowing four-probe transport measurements. The contact pads were defined on the (Ga,Mn)As layer by e-beam lithography, evaporation of W and Au, and lift-off. Subsequently, the wires and constrictions were defined by negative e-beam lithography. Cl2-based dry-etching was used to etch through the (Ga,Mn)As, leaving (Ga,Mn)As underneath the resist and metallization. The long axis of the island is oriented along the - (or equivalent) direction of the (Ga,Mn)As, which is near the magnetic easy axis of this layer, as determined by SQUID and consistent with Ref. roukes . By leaving the resist on the sample we can further etch it and narrow the constrictions.

We study symmetric double constriction samples in order to avoid complications associated with thermoelectric voltages. Moreover, the sample layout is such that the shape anisotropy causes the (magnetically isolated) inner island to switch at different fields than the outer wiresnucl . We verified by SQUID magnetometry that while the wide leads switch magnetization at around 15-20 mT (depending on lithographic parameters), the island exhibits switching fields of order 60-90 mT. The coercive field of the unpatterned epilayer is $``$ 8 mT.

MR measurements are carried out at 4.2 K in a He bath cryostat with a superconducting magnet. The magnetic field is applied parallel to the current direction. Four-probe DC measurements of the magnetoresistance are performed at constant voltage using zero-offset voltage- and current-preamplifiers. During the measurement the magnetic field is varied from full negative saturation of the material to full positive saturation and back.

The inset in Fig. 1 schematically describes the expected MR of our device. Sweeping the field from large negative to positive values causes the outer wires to switch first, inducing antiparallel alignment of the island and the wires. In this state, DW are present at the constrictions and the resistance is increased. At larger positive fields the magnetization of the island is also reversed, leaving all areas aligned in parallel and the resistance returns to its original value. The magnetization reversal is hysteretic, leading to the depicted spin-valve-like MR.

Experimentally, all samples exhibit the expected MR. The amplitude of the effect depends strongly on the resistance of the constrictions. Fig. 2(a) shows the MR of a representative sample with a four-terminal resistance of $``$48 k$`\mathrm{\Omega }`$. The MR is spin-valve-like with the correct hysteresis, and the critical fields agree with the SQUID results. We regard this as evidence that our design successfully incorporates both shape-anisotropy-controlled switching and strong pinning of DW in the constrictions. For comparison, the inset of Fig. 2(a) shows the MR of a 400 nm wire without constrictions which shows only 0.3% bulk-like anisotropic MR, with switching at $``$20 mT. In Fig. 2(a), the maximum MR is $``$1.5%.

Because the resist remains on the sample after fabrication, we can apply a further dry etching step producing a sample with narrower constrictions, but with approximately the same width of island and leads. Applying this procedure to the sample of Fig. 2(a), the additional etching brought the device resistance up to $``$78 k$`\mathrm{\Omega }`$ while the MR effect increased to about 8% (Fig. 2(b)).

The fine structure on the peaks of the as-fabricated sample was completely reproducible, whereas for the re-etched sample it varied from sweep to sweep as evident in the figure. This could suggest the presence of multiple pinning sites near the constrictions. The different pinning sites would yield different geometrical confinement of the DW, thus altering their resistance. Because of the extremely small dimensions realized here, impurities and side wall roughness by etching are likely causes.

A very strong increase in MR is obtained when the constrictions are etched still further. In Fig. 3 we plot the MR of a sample with a zero field resistance of 4 M$`\mathrm{\Omega }`$ and a positive MR of nearly 2000%. From the bulk resistivity of the material, we estimate that the leads and wires in our devices contribute only $``$40 k$`\mathrm{\Omega }`$ to the resistance, implying that the constructions now act as tunnel barriers. In addition, the I-V characteristics of the sample are strongly non-linear with a quadratic dependence of the conductance on bias that is characteristic of tunneling transportbrinkman . This suggests that the observed very large MR could be due to TMR. We also note that, in contrast with the results in Fig. 2, we now observe a hysteretic signal around zero field (although a major jump in resistance still occurs at the 20 mT expected from the wide leads). This observation is also consistent with the presence of tunnel barriers, in that these cause a magnetic decoupling of the island and the wires. We suggest that in Fig. 3 the magnetization of the island no longer switches by the introduction of a DW through a constriction but rather by magnetic rotation, which explains the MR at zero field: If the structure is not perfectly aligned along the easy axis, the relatively wide wires will be magnetized slightly off-axis at zero field while the narrow island is fully dominated by shape anisotropy, so that the relative alignment is not fully parallel. In wider constrictions the magnetic coupling prevents this effect and MR is observed only at finite fields.

We can understand the above observations in a unified manner by assuming that etching causes a gradual depletion of the carrier density at the constrictions. Dry etching of semiconductors damages the walls of the epilayer, and the resulting charged impurities induce sidewall depletion of the interior of the semiconductor, a mechanism that clearly will be most effective at the narrowest parts of the structure, i.e. at the constrictions. In the numerical estimates below, we assume both constrictions have equal resistance. Note that deviations from this assumption have only minor effects on the drawn conclusions.

As a first order approximation, we assume that the DWR is given by the expression of Valet and Fertvalet  for the spin-accumulation-induced resistance at an abrupt junction between two regions of opposite magnetization:

$$\mathrm{\Delta }R2\beta ^2\rho ^{}\mathrm{},$$
(1)
where $`\rho ^{}`$ is the spin-symmetric bulk resistivity in the magnetic material, $`\mathrm{}`$ is the length of the constriction and $`\beta `$ is the spin polarization in the constriction given by

$$\beta =\frac{N_{}(E_F)v_FN_{}(E_F)v_F}{N_{}(E_F)v_F+N_{}(E_F)v_F}=\frac{E_FE_F}{E_F+E_F}.$$
(2)
Here we have used that in a parabolic band model both the densities of states at the Fermi energy $`N_,(E_F)`$ and the Fermi velocitiesvf2  $`v_{F,}`$ are proportional to $`k_F`$ (the arrows refer to the spin subbands). From Eq. 2 one can see the evolution of the resistance and MR with etching. We assume that etching depletes the (Ga,Mn)As, so that $`E_F`$ and $`E_F`$ are reduced, but the exchange splitting $`\mathrm{\Delta }E=E_FE_F`$ remains roughly the same. Hence the numerator of Eq. 2 does not change, but the denominator gets smaller and the polarization increases. We now insert this expression for the polarization into the Valet-Fert expression (Eq. 1) for the MR. Calculating the resistance $`R_c`$ of the constrictions by substracting $``$40 k$`\mathrm{\Omega }`$ from the device resistance, and taking a reasonable valuemcd  of $`\mathrm{\Delta }E30`$ meV for the exchange splitting, we reproduce the observed MR in Fig. 2a using the Fermi energy of 150 meV found for the unetched sample, while the data in Fig. 2b imply a reduction of the Fermi energy to about 90 meV. These values seem quite reasonable, but should, given the many uncertainties and approximations involved, only serve as a rough indication of what may be going on in the sample. Note that the spin polarization of some 20% obtained from these numbers is only a lower limit estimate of the bulk value: we know that the transport mean free path $`l_t`$ of the holes is shorter than the dimensions of the constriction, so we can assume that substantial spin relaxation is taking place.

We now turn to the data in Fig. 3b, where the constrictions clearly are in the tunneling regime. It is tempting to try to model the observed MR in terms of Julliere’s TMR model.julliere  In order to explain a 2000 % MR signal, the Julliere model requires a spin polarization of the contacts of ca. 95 %, much larger than the values found above. This large discrepancy suggests that the model in Ref. julliere  is not applicable here and we have therefore adopted a different approach to modelling the tunneling regime. We also note that a model of the MR effect due to the joint effects of both constrictions petukhov  can be ruled out because of the long distances between the barriers and the very short mean free path. We can therefore safely consider a model based on two independent barriers.

We assume that the etching process creates a shallow barrier of parabolic shape between the two regions of (Ga,Mn)As, as shown in Fig. 3a. We define the barrier height for the majority-spin holes above the chemical potential $`\mu `$ as $`E_B`$. If the barrier is very thin, such that the hole wave functions can penetrate into the barrier region and continue to couple the Mn spins, then it is reasonable to assume that $`\mathrm{\Delta }E`$ in the barrier region is the same as in the bulk. This results in a barrier for minority-spin holes that is higher ($`E_B+\mathrm{\Delta }E`$) than for majority-spin holes ($`E_B`$). As a consequence of the non-abrupt barrier, the thickness of the barrier for minority-spin holes, $`L_P`$, will be greater than that for majority spin holes, $`L_P`$.

In the antiparallel situation depicted on the right, however, the barriers for the two spin channels are the same at approximately $`E_B+\mathrm{\Delta }E/2`$, and their thicknesses are also the same, $`L_A(L_P+L_P)/2`$. The transmission probability $`T`$ through a parabolic barrier issze

$$T=\mathrm{exp}\left(\frac{\pi m^{(1/2)}E_H^{(1/2)}L}{2^{3/2},\mathrm{}}\right)$$
(3)
where $`E_H`$ and $`L`$ are the height and thickness of the barriers respectively, and $`\mathrm{}`$ is the reduced Planck constant. With $`m^{}0.5m_o`$, we have $`T=\mathrm{exp}\left(3.0E_H^{(1/2)}L\right),`$ where $`E_H`$ is in eV and $`L`$ in nm.

We now estimate the values of the parameters required to match the experiment. Assuming that the parabolic shape of the barrier is the same for all situations, there is a uniform relationship between $`L`$ and $`E_H`$ of the form $`L=(\alpha E_H)^{(1/2)}`$, where $`\alpha `$ is constant. This implies that $`T=\mathrm{exp}\left(3.0\alpha ^{(1/2)}E_H\right)`$. From the experiment we have $`T_P/T_A=20`$, and so we choose $`\alpha `$ to satisfy this. With $`\mathrm{\Delta }E=30`$ meV, this yields $`\alpha =4400`$ eV-2, independently of the barrier height for majority spin holes $`E_B`$. We can estimate $`E_B31`$ meV from the resistance of the constrictions. Just as a gauge, the thicknesses of the barriers are then $`11.7`$ nm for the parallel majority case, $`14.3`$ nm for the antiparallel case, and $`16.4`$ nm for the parallel minority case. These numbers all seem reasonable.

A key element of this analysis is that the minority and majority carriers deplete at different positions in the constriction. Depletion at the edge of a (Ga,Mn)As film differs considerably from depletion in the bulk, since Mn spins at the edge remain coupled through the remaining holes to Mn spins which lie effectively within the bulk. The presence of these nearby bulk-like oriented Mn spins produces, through the mediating holes, a large exchange field on the Mn spins at the edge. This in turn induces them to order at local hole concentrations which, in the bulk, would otherwise not lead to ferromagnetism. Hence we argue that magnetically ordered Mn spins, producing an exchange splitting for the holes similar to that in bulk, are present at the edges of the sample where the local hole concentrations is much lower than the bulk.

Finally, we turn to the observation of a voltage dependence shown in the inset of Fig. 2b. The relative amplitude of the MR peak decreases exponentially with increasing bias voltage Vbias. A qualitatively similar behavior is observed for all samples with constrictions that are not in the tunnelling regime. A bias dependence of the current across a ferromagnetic DW was discussed theoretically in Refs. flatte , but not observed previously. By analogy to a semiconducting p-n junction, an exponential increase ($`\mathrm{sinh}(eV_{bias}/kT)`$) of the (electrical) current through the DW is expected. This should lead to an exponential suppression of the MR signal as observed in the inset of Fig. 2. However, instead of the expected slope of e/kT, we find an activation energy of around 11 meV, or 35 kT. Part of the discrepancy stems from the fact that not all the applied bias drops across the constrictions (a factor of 2 in this case), and another factor of 2 comes from the fact that our devices have two constrictions. Nonetheless, the deviation from Refs. flatte  is so large that we assume additional physics is at work here. We suggest that the discrepancy may be caused by high electric fields at the constrictions at elevated bias. In semiconductors, such fields may cause drift effects to dominate the transportyu , causing a strong reduction of the up-stream spin diffusion length. This would substantially modify the exponential increase of the current. It would be of interest to investigate such effects in detail, both experimentally and theoretically.

We thank A. Brataas and C. Timm for useful discussions. This work was supported by the BMBF, the European Commission (the FENIKS consortium), and the DARPA SPINS program.


format: text

text:
# Magnetization of a two-dimensional electron gas with a second filled subband

## Abstract

We have measured the magnetization of a dual-subband two-dimensional electron gas, confined in a GaAs/AlGaAs heterojunction. In contrast to two-dimensional electron gases with a single subband, we observe non-$`1/B`$-periodic, triangularly shaped oscillations of the magnetization with an amplitude significantly less than $`1\mu _\mathrm{B}^{}`$ per electron. All three effects are explained by a field dependent self-consistent model, demonstrating the shape of the magnetization is dominated by oscillations in the confining potential. Additionally, at 1 K, we observe small oscillations at magnetic fields where Landau-levels of the two different subbands cross.

When an extra degree of freedom is added to a two-dimensional electron gas (2DEG), many-body interactions can lead to the formation of novel electronic grounds-states at the crossings of the different energy-levels in the system. Das Sarma and Pinczuk (1997) Two-dimensional electron gases with crossing energy levels can be realized in a variety of systems with different relative sizes of orbital and spin effects, Coulomb energy, and different coupling between the components. Their study has lead to the discovery of many correlated Quantum Hall states, Piazza et al. (1999); De Poortere et al. (2000); Jaroszyński et al. (2002) and much effort is put into unravelling the energy-level structure of these systems.

One way of realizing such a 2D system, is to increase the electron density of a III-V 2DEG such that a second subband becomes occupied. Dual-subband systems realized in a quantum well have recently been studied within this context. Pan et al. (2001); Muraki et al. (2001) A similar system is the dual-subband 2DEG in a GaAs/AlGaAs heterojunction. In transport studies, the multi-subband 2DEG is generally assumed to be a superposition of single 2DEGs: the Landau-level structure is a superposition of Landau fans, separated by the inter-subband spacing calculated self-consistently at zero magnetic field. Ando et al. (1982); Stern and Das Sarma (1984)

In this paper we study the magnetization ($`M`$) of a dual-subband 2DEG. This is a way to directly probe a thermodynamic property, the chemical potential, $`\mu `$. For two-dimensional systems, the Maxwell relation between $`M`$ and $`\mu `$ is reduced to a proportionality. Since the Fermi energy ($`E_\mathrm{F}`$) is equal to $`\mu `$ at low temperatures, the magnetization directly reveals changes in the size as well as the shape of the Fermi energy: $`\mathrm{\Delta }M=\frac{N}{B}\mathrm{\Delta }E_\mathrm{F}`$, where $`N`$ is the total number of electrons.

The magnetization of multi-subband 2DEGs has already attracted some attention, both theoretically Alexandrov and Bratkovsky (2001, 1996) and experimentally, Shepherd et al. (1999) however, these studies have focussed on very high density systems with three or more filled subbands. In this regime changes in the energy-gap between the subbands can be ignored. In this paper we focus on the effect of the filling of a only a second electronic subband on the Fermi energy.

Quantum oscillations in the magnetization of a single 2DEG are well known to be characterized by strictly $`1/B`$-periodic saw-tooth-like oscillations with an amplitude of one effective Bohr magneton ($`\mu _\mathrm{B}^{}`$) per electron. Schoenberg (1984); Eisenstein et al. (1985); Wiegers et al. (1997) Here we will show that this is no longer the case in a multi-component system. Due to a self-consistent, magnetic field dependent redistribution of electrons between the subbands inside the heterojunction, the amplitude of the oscillations becomes considerably reduced, the saw-tooth-like steps are broadened into triangles, and the $`1/B`$-periodicity is lost. Additionally we find that extra magnetization minima appear at low temperature at the Landau-level crossings of the two subbands.

We study the magnetization of two samples with different electron densities, realized in GaAs/AlGaAs heterojunctions grown by molecular beam epitaxy: a single subband 2DEG (sample 1) and a high density 2DEG (sample 2) where two electronic subbands are occupied. Sample 1 has a density of $`4.810^{11}`$ cm-2, and a mobility of $`2.210^6`$ cm2/Vs. Our high electron density sample 2 has a carrier concentration of $`8.010^{11}`$ cm-2 and a mobility of $`1.410^6`$ cm2/Vs. Most of the electrons ($`7.410^{11}`$ cm2, deduced from transport measurements on a reference sample) remain in the lowest subband; the small remaining fraction occupies the second subband. The magnetization experiments were performed using a torsional magnetometer with optical angular detection. Schaapman et al. (2002)

Figure 1 shows the magnetization of sample 1, the single subband 2DEG, as a function of filling factor ($`\nu =hn/eB`$). It displays oscillations periodic in $`\nu `$, i.e. $`1/B`$-periodic oscillations. The steps at the lowest filling factors, where Landau-level broadening has the least influence, are saw-tooth shaped, and, with decreasing temperature, the amplitude saturates to $`1\mu _\mathrm{B}^{}`$ per electron. Although the steps are rather sharp, even at 1.2 K they still have a small, finite width indicating a finite density of states (DOS) in between Landau-levels. Wiegers et al. (1997); Schwarz et al. (2002)

Apart from the clear steps assigned to the Landau gap at even integer filling factors, at 1.2 K (figure 1(b)) additional features appear at odd integer filling factors. They are attributed to the opening of a spin-gap, significantly enhanced due to exchange interaction. Nicholas et al. (1988); Grundler et al. (1998); Meinel et al. (1999)

The second sample, with two filled subbands, displays the magnetization plotted by the solid line in figure 2. The dashed line is a theoretical calculation and will be discussed later. Also for this sample the magnetization oscillates as a function of inverse magnetic field. Closer inspection of the data, however, reveals three distinct differences compared to the single 2DEG. First, the oscillations are no longer sawtooth-like, but instead they are triangularly shaped. Secondly, we find that the oscillations are no longer strictly periodic in $`1/B`$. This becomes clear when we see that while $`\nu =4`$ coincides with an oscillation minimum, there is an increasing discrepancy, and $`\nu =14`$ actually coincides with an oscillation maximum. Finally, the amplitude of the oscillation is about $`0.5\mu _\mathrm{B}^{}`$ per electron, even for the lowest filling factors. This value is significantly less than the $`1\mu _\mathrm{B}^{}`$ per electron observed in figure 1b for the single 2DEG and it remains at this level even for the lowest temperatures (see fig. 3).

In order to understand the behavior of the magnetization, it is important to realize that in a heterojunction the confining potential of the 2DEG is formed by the electrons themselves. In a dual-subband 2DEG redistribution of charge over the two subbands can occur when a magnetic field is applied, resulting in a potential that is not fixed as a function of magnetic field. The wavefunction of the second subband is much more extended than that of the first one, therefore even a small change in its occupation can have profound effects. Since the occupation of the two subbands depends on the magnetic field that quantizes the DOS into Landau levels, and since the shape of the confining potential, the inter-subband spacing and the spacial charge distribution are interdependent, the Schrödinger and Poisson equations have to be solved self-consistently for each value of the magnetic field. Sanchez-Dehesa et al. (1987); Trott et al. (1989)

In our modelling we keep the electron density fixed and assume a Gaussian broadened DOS with a width that increases with the square root of the magnetic field. As the (bare) spin-splitting is too small to have an effect, it is neglected in the calculations. The Landau-level broadening is our, albeit only, fit-parameter. We note that the calculated amplitude decreases with increasing broadening, however, even in the limit of non broadened Landau levels, the calculated amplitude of the oscillations is only $`0.7\mu _\mathrm{B}^{}`$ per electron.

Using the self-consistent model, we have calculated the Fermi energy as a function of magnetic field, from which the magnetization follows directly through the Maxwell proportionality. The dashed line in figure 2 shows the resulting magnetization for a Gaussian Landau-level broadening with $`\mathrm{\Gamma }=0.2\sqrt{B}`$ meV. It is in very good agreement with the experimental data as it reproduces all three observed effects: triangular shape, non-periodicity, and the reduced amplitude.

Inspection of the self-consistent field-dependent modelling of the high-density 2DEG in detail reveals two important points. First, although the number of electrons in the highest subband is small, it remains populated up to high magnetic fields. Secondly, the shape of the Fermi energy (and thus magnetization) is determined by the oscillations in the inter-subband spacing, caused by self-consistent redistribution of electrons over the two subbands. A Landau-level scheme for the dual-subband 2DEG, resulting from the self-consistent model, is depicted in figure 4. While levels originating in the lower subband (solid lines) are linear functions of the magnetic field, the Landau-levels of the higher subband (dashed lines) oscillate according to the inter-subband spacing. Above 1.5 T only the lowest Landau-level of the second subband is populated and the Fermi energy lies continuously within it. It can be clearly seen that the oscillations of the confining potential due to the redistribution of the electrons have a large effect on the energy-level structure.

At this stage it is interesting to remark that consequently the wideness of the magnetization step is mainly caused by this redistribution and not by a finite DOS between two Landau-levels as suggested for a single subband 2DEG. Wiegers et al. (1997); Schwarz et al. (2002) Only in the region nearing $`\nu =4`$ there is, scarcely visible in figure 2, a kink followed by a sharp step, whose finite width is related to this small, extra DOS. These features are not at all visible on the other downward slopes, where the width is determined by the electron redistribution, and including the extra DOS does not influence the shape of the calculated Fermi-energy.

When reducing the temperature to 1 K, additional minima appear in the 2DEG magnetization (figure 3) around filling factor $`\nu =9.6`$ and $`\nu =12.0`$, and $`\nu =14.2`$. Interestingly enough these filling factors coincide with positions where two Landau levels originating from the two subbands cross, indicated by down-arrows in fig. 3.

On the flanks of the triangular oscillations a series of crossings occurs between the lowest Landau-level of the higher subband and Landau-levels with decreasing index of the lower subband as the magnetic field increases (see figure 4). The addition of spin-splitting to this (single-electron) picture results in energy-level schemes as depicted in the inset of fig. 4 for one of the crossings. When spin-splitting is taken into account, there is not a single level-crossing, but a small region where the levels with different spin consecutively cross each other. In this region spin up and spin down do not alternate for increasing energy: the two spin down levels are lowest in energy, the spin up levels the highest.

Although the Landau levels in figure 4 are represented by discrete lines, they are in fact considerably broadened, creating an overlap and giving the electrons some freedom to distribute themselves over the available energy-levels. We suggest this enables electrons to form a novel electronic ground-state that is spin-polarized in the crossing region. Creation of this polarized state would be favored by the system, because exchange interaction significantly reduces the ground-state energy. When the energy gain exceeds the broadening of the energy-levels, the enhanced gap shows up as a minimum in the magnetization.

Although it is evident that at lower magnetic fields (filling factors higher than 14, up-arrows indicate the positions of the level crossings), extra structure cannot be seen due to the broadness of the Landau-levels, extra structure is also too small to be observed at the Landau-level crossings of the lowest filling factors (up-arrows in fig. 3). Clearly the picture of the crossing region sketched above is not yet complete and further experimental and theoretical investigation of this many-body effect is required.

In summary, we have measured the magnetization of the coupled 2DEGs in a dual-subband 2DEG in a GaAs/AlGaAs heterojunction. We find that the de Haas-van Alphen oscillations are changed in three ways compared to those of the single 2DEG. The shape is triangular, the oscillation amplitude is reduced to $`0.5\mu _\mathrm{B}^{}`$ and the oscillations are no longer periodic in $`1/B`$. This behavior is well described by a self-consistent model, taking into account changes of the confining potential with magnetic field. It shows the shape of the Fermi energy and consequently the magnetization is entirely dominated by the oscillations in this potential due to redistribution of electrons over the two subbands. We observe additional magnetization minima at 1 K, which occur at magnetic fields corresponding to the positions where Landau-levels originating in the two different subbands cross. These minima possibly originate from a reduction of the total energy by the formation of a novel, exchange enhanced electronic state at the level-crossing.

###### Acknowledgements.

We would like to thank M. Elliot for fruitful discussion.


text:
# Crossover from Poisson to Wigner-Dyson Level Statistics in Spin Chains with Integrability Breaking1footnote 11footnote 1Submitted to Physical Review B

## I Introduction

The conjecture that the statistical properties of energy levels of chaotic quantum systems may be described in terms of the theory of random matrices is widely accepted in various fields of physics. meh  This however is not a universal property of all complex interacting systems. One example to the contrary is provided by the class of the so-called integrable models, bax  where the behavior of the system is completely described by a large (infinite in the thermodynamic limit) set of conserved quantities. One consequence is that the level-spacing distribution $`P_\mathrm{\Delta }(E)`$ in the case of integrable models is the Poisson distribution ($`\mathrm{\Delta }`$ denotes mean level spacing),

$$P_\mathrm{\Delta }(E)=\frac{1}{\mathrm{\Delta }}e^{E/\mathrm{\Delta }},$$
(1)
whereas in random-matrix theory the distribution takes the Wigner-Dyson form,

$$P_\mathrm{\Delta }(E)=b_\beta \left(\frac{E}{\mathrm{\Delta }}\right)^\beta e^{a_\beta E^2/\mathrm{\Delta }^2},$$
(2)
where $`\beta =1,2,4`$ correspond to orthogonal, unitary and symplectic ensembles respectively and guh  $`b_1=\pi /2`$, $`a_1=\pi /4`$; $`b_2=32/\pi ^23.24`$, $`a_2=4/\pi `$; $`b_4=262144/729\pi ^311.6`$, $`a_4=64/9\pi 2.26`$.

Other statistical properties (for example, the evolution of levels under an external perturbation faa ) also differ for the two cases. One important class of external perturbations is the application of a voltage. The difference in response in this case leads to spectacular differences in transport properties of integrable and non-integrable models. Integrable models have been argued to have an infinite conductivity even at high temperatures, essentially because a typical level has a large response to a voltage, whereas non-integrable models have a finite conductivity because a typical level has a small response. zot ; zo2

While these basic properties have been established for the two generic cases of integrable and non-integrable models, the crossover between these two limits as an integrability-breaking interaction is turned on has not to our knowledge been carefully studied, nor have the implications of the crossover for the finite-size conductivity of nearly integrable systems been determined. Two of us, with N. Andrei, presented a few numerical results in a paper mainly concerned with the charge transport of integrable systems. us1  However, the significance and interpretation of these results was not clear. Song and Shepelyansky son  studied the effects of a random potential on level statistics of 2D interacting Fermions and found evidence for a localization-delocalization transition. However, in their case, the physics of the transition is due to the disorder and thus is different from the situation in integrable models. Earlier work by DiStasio and Zotos distasio95  noted a crossover between Poisson to Wigner-Dyson only for the low-energy part of the spectrum and did not address scaling with system size.

In this paper we will fill these gaps by providing numerical results for finite-size chains with Hamiltonian given by the (integrable) XXZ model plus an integrability-breaking perturbation $`\delta H`$. Our principal results are computations, for finite-length chains, of the crossover from behavior characteristic of Poisson to behavior characteristic of Wigner-Dyson statistics in various statistical measures. These crossovers do not display an obvious universality in the sense that different measures show different behavior depending on the XXZ asymmetry parameter and system size.

Our computations are performed for finite-size systems. An important issue is the behavior in the limit of thermodynamically large system. Extrapolation to the thermodynamic limit proves to be ambiguous for most of the measures we employ (namely, we cannot rule out a saturation of the crossover scales as functions of the system size for chains much longer than those considered in this study), but the data suggest that all the crossover scales vanish at infinite system size.

The rest of the paper is organized as follows. First we discuss the model used in numerical calculations and in particular define numerically the value of the integrability-breaking parameter at which a gap appears in the spectrum. All further considerations will be devoted to the gapless regime. Then we discuss the level-spacing distribution and the correlator of level velocities. The latter is related to the parametric statistics of the system and also to its transport properties. For disordered systems, the correlator of level velocities was shown to correspond to the dimensionless conductance of the system, while if one restricts the analysis to periodic boundary conditions only (see below) it coincides with the Drude weight. Discussion of the Drude weight concludes the paper.

## II The model

We study the effect of integrability breaking on the physical properties of a spin chain. The integrable model we consider is the XXZ chain defined on a $`N`$-site ring with periodic boundary conditions in the presence of external flux $`\varphi `$ threading the ring:

$`H_{XXZ}={\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left(e^{i\phi /N}S_i^+S_{i+1}^{}+e^{i\phi /N}S_i^{}S_{i+1}^+\right)`$
$`+{\displaystyle \underset{i=1}{\overset{N}{}}}J_1S_i^zS_{i+1}^z.`$ (3)
(Alternatively, the flux can be gauged out to the boundary, resulting in twisted boundary conditions.) As is well-known, statistical properties of integrable models are governed by the Poisson distribution, Eq.(1). Transport properties of the model can also be inferred from studying the energy levels of the model, namely by their response to the flux $`\varphi `$. At zero temperature the behavior of the ground-state energy of the system under slow variations of the flux determines the Drude weight or the stiffness kon  $`𝒟_s`$ as

$$𝒟_s=\frac{N}{2}\frac{^2E_0}{\varphi ^2}|_{\varphi 0}.$$

Non-vanishing $`𝒟_s`$ signals ballistic transport in the system. For the XXZ model at $`T=0`$ this is the case sha  for $`1<J_1<1`$, where excitations of the system are gapless. If $`|J_1|>1`$, then the excitation spectrum of the model is gapped,descloiseaux  and $`𝒟_s=0`$. At finite temperatures the above expression for the stiffness can be generalized zot ; zo2 ; us1  to $`𝒟_s=D_1+D_2`$, where

$$D_1=\frac{N}{2\beta }\frac{1}{𝒵}\frac{^2𝒵}{\varphi ^2}|_{\varphi 0}$$

vanishes in the thermodynamic limit, gin  and the rest is positive:

$$D_2=\frac{\beta N}{2}\frac{1}{𝒵}\underset{n}{}\left(\frac{E_n}{\varphi }\right)^2|_{\varphi 0}e^{\beta E_n}.$$
(4)
In the gapless phase of the XXZ model it has been shown zot ; zo2 ; us1  that ballistic transport persists to finite temperatures in the sense that $`D(N)=lim_T\mathrm{}D_2(N)T>0`$. The infinite-temperature limit of this result implies that for a typical level $`dE_n/d\varphi 1/\sqrt{N}`$. At the antiferromagnetic Heisenberg point $`J_1=1`$ the model still has gapless excitations, but results of Fabricius and McCoy mac  suggest that $`D_2`$ vanishes (slowly) as the system size increases. Numerical results of Narozhny et al.us1  were consistent with this suggestion, but the limited range of system sizes attainable precluded a definite statement.

Integrability breaking is introduced by adding the term with next-neighbor coupling

$$\delta H=\underset{i=1}{\overset{N}{}}J_2S_i^zS_{i+2}^z.$$
(5)
This term should be contrasted to that considered by Eggert egg  insofar as it is explicitly not SU(2) invariant. However, away from the Heisenberg point the effect of the interaction Eq. (5) is similar to that of its SU(2)-invariant counterpart: (i) it breaks the integrability of the system, and (ii) for large enough values of $`J_2`$ it causes the system to dimerize, so that the spectrum acquires a gap. The critical value of $`J_2`$ at which the gap opens is of course different from the $`0.24`$ found in Ref. egg, . Our numerical estimates est  suggest a value $`J_2^{(c)}1.1`$ that is a weakly increasing function of both system size and $`J_1`$: for $`N=18`$ and $`J_1=0.2`$, for example, the gap appears to open at $`J_2=1.13\pm 0.01`$. The gap opening limits the range of values of $`J_2`$ under consideration, as we are interested only in properties of the gapless phase; indeed, the measures we consider presently begin showing different behavior for $`J_2>J_2^{(c)}`$. Similarly, the parameter region considered by Faas et al.faa  belongs to the gapped regime, which accounts for certain differences in the behavior of the level statistics reported in Ref. faa,  and in the present paper.

In this paper we study the eigenvalues of $`H=H_{XXZ}+\delta H`$ and their evolution under the change of $`\varphi `$ for the above models with varying $`J_{1,2}`$ and system size. For the non-integrable Hamiltonian $`H_{XXZ}+\delta H`$ we use exact numerical diagonalization to construct the level-spacing distribution and level auto-correlation functions and to evaluate the stiffness $`D_2`$, Eq. (4). The use of exact numerical methods is motivated by the need to obtain the whole spectrum of the model in order (i) to study the statistical properties of the spectrum and (ii) to study the stiffness Eq. (4) at infinite temperature. The drawback of the method is the limitation to small system sizes (we present results for chains of up to $`20`$ sites). For finite system sizes we obtain a detailed characterization of the crossover.

## III Level-spacing distribution

We begin with a brief discussion of the integrable case. The level-spacing distribution for $`J_2=0`$ is the Poisson distribution (shown in Fig.1 by the left solid curve). This illustrates the fact that the integrable system has so many conservation laws that levels essentially do not repel each other. To characterize transport properties of the system we show in Fig. 2 the quantity $`D(N)=D_2(N)T`$ at $`T\mathrm{}`$ for different system sizes and different values of the integrable interaction $`J_1`$ (dashed lines in Fig. 2). $`D`$ is seen to be almost size-independent for the cases $`J_1<1`$, in agreement with previous work, us1  while a weak size dependence is evident in the Heisenberg case $`J_1=1`$. Although this dependence appears to have a positive $`y`$-intercept, we believe that the system size in this study is still too small to make a definite statement regarding the behavior of the Heisenberg model in the thermodynamic limit.twothirds

We turn now to the case of broken integrability. As the integrability-breaking term Eq. (5) is added to the Hamiltonian, energy levels immediately start to repel,rep  and as a consequence immediately $`P_\mathrm{\Delta }(0)=0`$ so that the distribution acquires a peak. As illustrated in Fig. 1, increasing $`J_2`$ shifts the peak to the right until the distribution starts to look like the Wigner-Dyson distribution (shown in Fig. 1 by the right solid curve). At the same time the tail of the distribution changes from the exponential in Eq. (1) to the (asymptotically) Gaussian tail of the Wigner-Dyson distribution.

To quantify this crossover we show the evolution of the peak position and the characteristics of the tail with the change in $`J_2`$ in Figs. 3 and 4. Both exhibit similar features, although the estimates for the crossover scales extracted from the two are numerically different (see table 1 and insets in Figs. 3 and 5).

As shown in Fig. 3, the peak of the distribution grows from zero to the value characteristic of the Wigner-Dyson distribution and then saturates. To estimate the crossover scale $`J_2^{}`$, we fit the data by the hyperbolic tangent fit  of the form $`a\mathrm{tanh}(x/x_0)`$ with $`x_0`$ approximating $`J_2^{}`$. The inset shows the resulting values for $`J_2^{}`$ as a function of the system size. As we noted before, rep  we are restricting our attention to fixed values of the total spin $`S^z`$. However, for the purposes of the finite-size scaling, it makes more sense to compare data with the fixed ratio $`S^z/N`$. One way to see this is to recall that by means of the Jordan-Wigner transformation the spin chain can be mapped onto a system of spinless Fermions.us1  In the Fermion language, $`1/2S^z/N`$ corresponds to the filling fraction. Since it is not possible to keep the ratio $`S^z/N`$ exactly the same for all values of $`N`$ used in this paper, we choose to present the data for two sectors of fixed $`S^z`$ that are closest to the chosen value of $`S^z/N`$. Therefore the inset in Fig. 3 shows two data points for the $`N`$ other than $`N=18`$ (we chose $`S^z/N=1/6`$). The straight lines are just guides to the eye.

To analyze the evolution of the tail, we approximate the intermediate distributions (see Fig. 1) by

$$P_\mathrm{\Delta }(E)\mathrm{exp}\left[a\frac{E}{\mathrm{\Delta }}b\left(\frac{E}{\mathrm{\Delta }}\right)^2\right].$$

Clearly, for Eq. (1) $`a=1`$ and $`b=0`$, while for the orthogonal ensemble, Eq. (2) corresponds to $`a=0`$ and $`b=\pi /4`$. In Fig. 4 we show the evolution of $`b`$ (the fact that plotted values never reach $`\pi /4`$ is an artifact of the calculation). Fitting the curve to a hyperbolic tangent, we can extract an estimate for the crossover value $`J_2^{}(N=20)=0.27`$. This value differs somewhat from the one extracted from the peak position (for the same values of $`N`$, $`S^z`$, and $`J_1`$); see the inset in Fig. 3. The behavior of the tail characteristics with respect to changing system size exhibits the same trend as shown in the inset in Fig. 3 for the peak position: the characteristic scales tend to decay with increasing system size. The naive extrapolation of such a trend is consistent with a statement of vanishing $`J_2^{}`$ as $`N\mathrm{}`$; however, the data are insufficient to prove it.

Another way to quantify the evolution of the level-spacing distribution shown in Fig. 1 is to consider cumulants. cum  Their unbiased estimators (the Fisher statistics fis  $`k_n`$) are easily computed. For our normalized level spacings, the first cumulant (which is equal to the mean) is unity.

The cumulants of a distribution characterize its width (second cumulant, or variance) and shape;gau  beyond perhaps the fifth, numerical cumulants become too sensitive to outliers to be of much use. A study of the cumulants of a distribution is qualitatively similar to our foregoing study of the tails, but it turns out to be simpler numerically. In Fig. 5 we show the unbiased variance estimate, $`k_2`$, as a function of $`J_2`$ for system size $`N=20`$, $`J_1=0.2`$, $`S^z=3`$, momentum $`k=0`$. The theoretical limits should be $`1`$ for the Poisson distribution and $`\pi /41=0.2732`$ (the bottom of the scale) for Wigner-Dyson. The fact that the data points deviate from these ideal values is an artifact of finite sampling. Fitting the curve to a hyperbolic tangent, we estimate the crossover scale, shown in the inset as a function of size (for the same sequence of quantum numbers as in Fig. 3). The inset also shows the crossover scale for the fourth cumulant. As before, the scaling suggests (but does not establish) that the crossover scale $`J_2^{}`$ associated with either cumulant should vanish in the limit of infinite size. It is not possible to determine whether the different measures, $`k_2`$ and $`k_4`$ (we also looked at $`k_3`$), scale in the same way or differently with system size.

The second cumulant $`k_2`$ plays a role similar to that of the parameter $`\eta `$ used in Ref. son,  to estimate the overall “proximity” of the observed distribution to either the Poisson or the Wigner-Dyson limit. In that sense, Fig. 5 shows behavior similar to that found in Ref. son, , although the physics of the evolution of levels is quite different in our case.

## IV Elements of the parametric statistics

More information about the crossover to the chaotic behavior described by the Wigner-Dyson statistics can be extracted from the study of the autocorrelation functions. Here we will discuss the autocorrelation of level velocities

$$C(\varphi )=\frac{1}{\mathrm{\Delta }^2}\frac{E_i(\theta )}{\theta }\frac{E_i(\theta +\varphi )}{\theta }_{\theta ,i},$$
(6)
where the angular brackets indicate averaging over a set of levels and fluxes.faa

A typical form of the autocorrelation function Eq. (6) is shown in Fig. 6. For large values of $`J_2`$, this form resembles the universal correlator characteristic of chaotic systems. faa  However, in the crossover region, $`C(\varphi )`$ deviates from the universal form in a rather complex fashion, which makes a quantitative analysis of the crossover difficult. Therefore we focus on two particular features of the curve, the turning point and $`C(0)`$. We note a feature of the autocorrelation curve shown in Fig. 6 that appears only as the integrability breaking is introduced: all curves for $`J_2>0`$ have a non-zero point of inflection as the autocorrelation decreases with the increase of $`J_2`$, but the autocorrelation function of the integrable system does not have such an inflection point. Consequently, the behavior near zero flux changes from linear ($`C(\varphi )C(0)\varphi `$) for the integrable case to quadratic for $`J_2>0`$.

The autocorrelation function of level velocities at zero flux difference $`C(0)`$ is somewhat similar to the stiffness Eq. (4), the differences being that $`C(0)`$ is also averaged over a set of fluxes, does not contain the extra factors of temperature and system size, and corresponds to a single sector of fixed $`S^z`$. However, in chaotic systems it is $`C(0)`$ that can be related to transport. faa  There it was argued to correspond to the dimensionless conductance.

In Fig. 7 we show the behavior of $`C(0)`$ as a function of $`J_2`$. Clearly, for finite systems, $`C(0)`$ exhibits a well-defined crossover. For the dataset presented in Fig. 7 ($`N=20`$, $`J_1=0.2`$, $`S^z=3`$), the correlator $`C(0)`$ decays as a $`5/2`$-power law after $`J_2`$ exceeds the value $`J_2^{C_0}=0.097\pm 0.004`$ (defined in Fig. 7 as a crossing point of the above power law — the straight line in the log-log scale — with the value at $`J_2=0`$).

The inset shows the crossover scale as a function of the size (in the same manner as the crossover scale extracted from the peak of the level-spacing distribution). The behavior is very similar to that in the inset in Fig. 3 (although numerical values of the crossover scales differ in the two cases). Both would be consistent with the statement that $`J_2^{}0`$ as $`N\mathrm{}`$; however such a conclusion cannot be ascertained on the basis of the data available.

## V Spin stiffness

Now we discuss the effect of the integrability breaking Eq. (5) on the stiffness $`D_2`$. In Fig. 2 we show three sets of data corresponding to three different values of the XXZ anisotropy parameter $`J_1`$. The data illustrate the following tendencies:

(i) For $`0<J_1<1`$ (represented by $`J_1=0.2`$; similar behavior is observed for other values) the data clearly show that very small integrability breaking (characterized by $`J_2=0.01`$) has little effect on the stiffness of the finite chains (which is to be expected). Moreover, the extrapolation to infinite size seems to result in a finite value for the stiffness in a manner similar to that of the integrable model.

(ii) For the Heisenberg model, small integrability breaking again does not have a pronounced effect; however, in this case (even though the extrapolation indicates a finite value for the thermodynamic limit), one cannot be certain of the behavior of the infinite chain.

(iii) For $`J_1=0`$ the situation is different: even a very small amount of integrability breaking leads to a sharp reduction in the stiffness for finite chains. The extrapolation to infinite size is also uncertain. It should be noted, however, that if one compares the behavior of two integrable cases, $`J_1=0`$ and $`J_1>0`$, then a similar picture arises (compare for example the two top dashed lines in Fig. 2). This has to do with the fact that the spin chain at $`J_1=0`$ can be mapped (by means of the Jordan-Wigner transformation) onto a system of free spinless Fermions and as such possesses more symmetries than even the integrable (but interacting) XXZ model.

(iv) When the integrability-breaking parameter is not small, the stiffness decays sharply (in fact, if we were to show a log-log plot, faster than any power law) with system size and clearly extrapolates to zero in the thermodynamic limit. This behavior is qualitatively independent of $`J_1`$.

Thus the data suggest that for finite chains there exists a “critical” value of $`J_2`$ (smaller than the point of the gap opening) beyond which the stiffness tends to vanish. This value is not universal, in the sense that it depends on $`J_1`$. This situation is illustrated in Fig. 8, where we show the dependence of the stiffness on $`J_2`$ for three values of $`J_1`$ and the fixed system size $`N=18`$. For large values of $`J_2`$, all three curves saturate to zero (although the one with $`J_1=1`$ does so faster). For $`J_1=0.2`$ the effect of small $`J_2`$ is rather weak, and the curve exhibits a clear crossover. For $`J_1=0`$ the crossover also appears; however, the value $`0.08`$ to which the curve tends as $`J_20`$ is much smaller than the value at $`J_2=0`$ exactly (which is $`0.25`$ and is thus outside the frame of the plot). The crossover is illustrated in the left inset. The right inset shows the change of the crossover scale with system size for $`J_1=0`$, $`0.2`$, and $`0.4`$. The behavior at $`J_1=1`$ is quite different and in particular does not show an obvious crossover (and is thus not represented in the right inset). This behavior might be related to the conjecture mac  that at the Heisenberg point the stiffness vanishes in the thermodynamic limit even without the integrability breaking. Alternatively, this can reflect the fact that the Heisenberg model is characterized by logarithmic correlations,gin  and thus the small systems considered in this paper are not representative.

## VI Discussion

Prior to performing the calculations one could have had two conflicting expectations for the behavior of the non-integrable system: (i) as soon as the integrability is broken the system becomes chaotic and as a result shows diffusive transport; (ii) there exists a “critical” magnitude of the integrability breaking that separates the chaotic regime from the one that retains some features of the integrable model, in particular, ballistic transport.

The latter picture has an analogy in the localization problem in disordered conductors. alt  The states of an integrable model can be visualized as well-defined localized points in the multidimensional space of the integrals of motion characteristic of the model. These points are well separated due to the quantization of the values of the integrals of motion. Consider now the effect of an infinitesimally small integrability breaking. One can certainly expect the points to spread out into fuzzy spots, but at the same time one might argue that unless the integrability breaking is strong enough, these spots do not overlap. In this regime the system retains some memory of the fact that it was indeed integrable before the extra interaction was turned on. When the integrability breaking is so strong that the spots overlap into a continuum, the system becomes fully chaotic.

The numerical analysis presented in this paper seems roughly consistent with the second possibility for finite chains: a small integrability-breaking term leads to behavior that is close to that of the integrable system. Quantities related to transport, the stiffness $`D_2`$ and the “conductance”, $`C(0)`$, exhibit a reasonably rapid crossover as functions of the strength of the integrability-breaking interaction. The crossover behaviors seem to be different for different quantities. Table 1, for example, illustrates some of this variability for the example of $`J_1=0.2`$, intermediate between the non-interacting model and the Heisenberg point. Similarly, the quantitative characteristics of the level-spacing distribution (namely the peak position and the tail parameter, see Figs. 3 and 4, or cumulants, see Fig. 5) exhibit similar crossovers (the corresponding scales are also included in table 1).

The one exception to this picture is the number of degenerate levels in the system rep  represented by $`P_\mathrm{\Delta }(0)`$. This measure exhibits a jump as infinitesimally small (numerically meaning of the order of the computer precision) $`J_2`$ is introduced (namely $`P_\mathrm{\Delta }(0)`$ vanishes, as illustrated in Fig. 1).

Conclusions for the thermodynamic limit are harder to draw from our data. The variation of crossover scale with system size indicated in the insets to the different figures suggests that the crossover scale vanishes in the limit of infinite size system rather than saturating at non-zero values for $`J_2^{}`$, but the limited range of sizes available to us, along with the absence of a theoretically justified extrapolation to the thermodynamic limit, precludes a definite statement. Constructing a theory of the approach to the infinite-size limit of chains with weak integrability breaking remains an important open problem.

## Acknowledgments

Instructive discussions with B.L. Altshuler, A.A. Nersesyan, and J.K. Looper are gratefully acknowledged. Numerical work was performed at the San Diego Supercomputing Center and the University of Michigan supercomputing facility through NPACI grant CSD268 and at the Research-Oriented Computing Center of the University of South Florida. DAR is a Cottrell Scholar of Research Corporation and also wishes to thank the Abdus Salam International Centre for Theoretical Physics for its hospitality.


text:
# Rate of photon production from hot hadronic matter

## Abstract

Thermal photon emission rates from hot hadronic matter are studied to order $`e^2g^4`$, where $`g`$ indicates a strong-interaction coupling constant. Radiative decay of mesons, Compton and annihilation processes for hadrons, and bremsstrahlung reactions are all considered. Compared to the standard rates from the literature, one finds two orders of magnitude increase for low photon energies stemming mainly from bremsstrahlung and then a modest increase (factor of 2) for intermediate and high energy photons owing to radiative decays for a variety of mesons and from other reactions involving strangeness. These results could have important consequences for electromagnetic radiation studies at RHIC.

Photons carry undistorted detail on strong-interaction dynamics since their mean free paths in hot hadronic matter are much larger than typical system sizes. For this reason, photon measurements in high energy nuclear collisions are deemed the best thermometers of the highest energy densities reached thereef ; es . While pions and protons scatter several times before leaving the collision zones, photons do not, and therefore originate from the entire spacetime hadronic source, and in their asymptotic momentum configurations. Quantitative studies of photon mean free paths in hadronic matter have indeed revealed $`\lambda `$ 2$`\times `$103 fm for a temperature 200 MeVkh\_mfp .

Since production rates are increasing functions of temperature, the inverse slope on a transverse momentum spectrum, for instance, is a good measure of the kinetic properties of the early (thermal) hadronic source. And yet, photons are produced during the entire reaction. Pre-equilibrium dynamics contribute, parton dynamics contribute, and then, hadron dynamics might be considered as a background to the more sought-after quark gluon plasma contributionsmtreview ; cgkhreview . For quantitative interpretation of photon measurements at the Relativistic Heavy Ion Collider (RHIC), a complete determination of all components, including the hadronic background, is clearly a prerequisite.

The early comparison of photon production from quark gluon plasma versus hadronic matter by Kapusta et al.kls91  revealed that the two phases have equal production rates at fixed temperature. Technically, the advancements here were made possible through an application of hard-thermal-loop resummation methodsrp ; bp  and effective hadronic field theory. This physically appealing result is however, somewhat discouraging in terms of allowing distinction between the two phases using photons. The equal-luminous property of the two phases has survived even after next-to-leading-order (NLO) QCD processes have been consideredpa . NNLO parton processes have even been analyzedpa00 ; pa00b ; amy01 , and the new rate is still consistent with the hadron rates. The inevitable conclusion is that photon production as a QGP diagnostic is not as definitive as was originally hoped. However, with the tremendous recent advancements in the determination of QCD processes—to three-loop order, including annihilation with scattering and also bremsstrahlung with multiple scattering interferenceslp ; migdal , it behooves theory to put the hadronic formalism on a more equal footing. The purpose of this paper is to report on photon production from hadronic matter computed at NNLO in the strong-interaction coupling constant, namely, to order $`e^2g^4`$.

The photon production rate to lowest order in the electromagnetic coupling and to all orders in the strong coupling can be compactly written asaw83 ; mt85 ; gk91

$`E_\gamma {\displaystyle \frac{dR}{d^3p_\gamma }}`$ $`=`$ $`{\displaystyle \frac{2g^{\mu \nu }}{\left(2\pi \right)^3}}\mathrm{Im}\mathrm{\Pi }_{\mu \nu }^R(p_\gamma ){\displaystyle \frac{1}{e^{\beta E_\gamma }1}},`$ (1)
where the retarded photon self-energy includes all possible hadron topologies—that is, to arbitrary loop order. The imaginary part corresponds to emission and absorption, of which the former is relevant here. Absorption was studied in Ref. kh\_mfp  and found to be negligible. At the one-loop level, the imaginary part corresponds to such lowest-order reactions as $`\omega \pi \gamma `$ and $`\rho \pi \gamma `$. There is of course a long list of other resonance decays to be considered, as discussed below. Typically, in order to quantify production rates, a model for the strong interaction (e.g. an effective Lagrangian) is necessary. However, the one-loop reactions have the desirable feature that the rate can be computed in a model independent way, requiring only the measured radiative decay rates. If such measurements are not available, then vector meson dominance plus a chiral Lagrangian is an appropriate approach.

For the process $`a1+\gamma `$, where $`a`$ and 1 are mesons with allowed quantum numbers and masses, the thermal rate can be written as

$`E_\gamma {\displaystyle \frac{dR}{d^3p_\gamma }}`$ $`=`$ $`{\displaystyle \frac{𝒩m_a^2\mathrm{\Gamma }(a1\gamma )}{16\pi ^3E_\gamma E_0}}{\displaystyle \underset{E_{\mathrm{min}}}{\overset{\mathrm{}}{}}}𝑑E_af_{BE}(E_a)`$ (2)
$`\times \left[1+f_{BE}(E_aE_\gamma )\right],`$
where $`𝒩`$ is the spin and isospin degeneracy for the parent, $`E_0`$ is the photon energy in the rest frame of the parent, $`E_{\mathrm{min}}=m_a(E_\gamma ^2+E_0^2)/2E_\gamma E_0`$, and $`f_{BE}`$ is the Bose-Einstein distribution to account for quantum statistics and medium-enhancement for species 1. If baryons were studied, e.g. $`\mathrm{\Delta }N\gamma `$, then Fermi-Dirac distributions would be appropriate and a minus sign in the second factor of the integrand would suppress phase-space occupation for species 1. The above integral is typically reported in the literature, and then a numerical exercise ensues. However, there is no need to leave it here since the integral has a closed-form expression with the result feynman\_integral

$`E_\gamma {\displaystyle \frac{dR}{d^3p_\gamma }}`$ $`=`$ $`{\displaystyle \frac{𝒩m_a^2\mathrm{\Gamma }(a1\gamma )}{16\pi ^3E_\gamma E_0}}f_{BE}(E_\gamma )`$ (3)
$`\times T\left[\mathrm{ln}\left({\displaystyle \frac{f_{BE}(E_{\mathrm{min}}E_\gamma )}{f_{BE}(E_{\mathrm{min}})}}\right)E_\gamma /T\right].`$
Empirical radiative decay widths for several low-lying resonances are availablepdg  and included in Table 1. The resulting thermal decay rates are displayed in Fig. 1.

Next-to-leading-order contributions enter at order $`e^2g^2`$ and come from two-loop self-energy structures. Coupling strength and phase space arguments suggest that $`\pi \rho \pi \gamma `$, $`\pi \pi \rho \gamma `$ and $`\rho \pi \pi \gamma `$ will be important contributors. The first channel, which is the Compton type process, will be considered here while the second and third will be discussed later in the higher orders, since they are essentially the on-shell components of the full, higher order contributions. For the Compton process, a particular interaction Lagrangian was assumed in Ref. kls91  describing $`\pi `$-$`\rho `$ dynamics, and then calibrated to $`\rho \pi \pi `$ decay. In terms of practical kinetic theory formalism, the rate for the general reaction $`a+b1+\gamma `$ is

$`E_\gamma {\displaystyle \frac{dR}{d^3p_\gamma }}`$ $`=`$ $`𝒩{\displaystyle \frac{d^3p_a}{(2\pi )^32E_a}f_a\frac{d^3p_b}{(2\pi )^32E_b}f_b\frac{d^3p_1}{(2\pi )^32E_1}\stackrel{~}{f}_1}`$ (4)
$`{\displaystyle \frac{1}{(2\pi )^32}}|\overline{}|^2(2\pi )^4\delta ^4\left(p_a+p_bp_1p_\gamma \right).`$
A numerical task remains: That is, to integrate Eq. (4) over the appropriate phase space. The result for $`\pi \rho \pi \gamma `$ is shown in Fig. 2.

There are many other contributions at the same order with roughly the same phase space properties. Certain strangeness reactions turn out to be strong contributors, as can be seen from the results plotted in Fig. 2. The coupling constants are fitted to empirical hadronic decays where available, and estimated using an SU(3) symmetry in a chiral Lagrangian in the absence of measurements.

The strongest contributor involving strangeness seems to be $`\pi K^{}K\gamma `$, although there are certainly others of some significance. Strangeness reactions of this type go beyond what has been done beforekh94 . Also, these two-loop reactions and others are currently being studied in a parallel investigation, broadened in scope to explore the issue of form factorsmcgillgroup .

A remark is in order here regarding the role of the $`a_1`$ meson. In the reaction $`\pi \rho \pi \gamma `$, the $`a_1`$ appears as an intermediate state, both in the $`s`$ channel and the $`u`$ channel. Early estimates reported that the $`s`$ channel dominated all others, even dominated $`\pi `$ exchangexsb . Later studies using a hidden local symmetry in a chiral Lagrangian were less clear regarding the role of the $`a_1`$ since no unique rate could be establishedsong . However, recent improvements describing $`a_1`$ dynamics with particular concern for the $`D/S`$ ratio in the $`a_1\pi \rho `$ scattering amplitude show a diminished importance of the role of the $`a_1`$ for photon productionsimonandcharles . In Fig. 1, the $`a_1`$ radiative decay has been included while in Fig. 2 only the $`\pi `$-exchange result for the Compton process $`\pi \rho \pi \gamma `$ is shown. This has the slight incomplete feature of ignoring the $`u`$ channel $`a_1`$ exchange and coherence effects, but at the same time has the advantage of removing any possibility for double counting.

Motivated partly by the realization that the QCD processes and the resulting photon rates have been computed with bremsstrahlung included (and the hadron rates up to now have not), one moves next to three-loop topologies. Three-loop self-energy graphs involve four on-shell hadron lines with a variety of possible kinematics, plus the photon. Some of these include processes of the type $`a+b1+2+\gamma `$: hadronic elastic scattering with bremsstrahlung. Based on general principles alone, these contributions are expected to be strong for low photon energies and quickly become less important for increasing photon energy. Indeed, they will even diverge as $`E_\gamma 0`$. When bremsstrahlung processes are included, there is immediately a danger of double counting. For example, $`\pi \pi \rho \gamma `$ and $`\rho \pi \pi \gamma `$ are already part of the full $`\pi \pi \pi \pi \gamma `$ coherent calculation with an intermediate $`\rho `$ meson. Even though the two separate reactions were generated from two-loop self-energy topologies (restricted, however, to on-shell vector mesons), they must not be added on top of the bremsstrahlung contributions. That would be double counting. The strategy taken here is the following. When two-loop contributions are subsets of three-loop contributions, the separate two-loop contributions will not be included.

A treatment of the bremsstrahlung mechanisms in terms of complete phase space is a bit tedious, and so to first assess their relative importance, a soft-photon approximation is adopted here. The rate is thenhge93

$`E_\gamma {\displaystyle \frac{dR}{d^3p_\gamma }}`$ $`=`$ $`𝒩{\displaystyle \frac{T^2}{16\pi ^4}}{\displaystyle \underset{z_{\mathrm{min}}}{\overset{\mathrm{}}{}}}𝑑z\lambda (z^2T^2,m_a^2,m_b^2)K_1(z)`$ (5)
$`\times \left[E_\gamma {\displaystyle \frac{d\sigma }{d^3p_\gamma }}\right]{\displaystyle \frac{R_2(s_2,m_a^2,m_b^2)}{R_2(s,m_a^2,m_b^2)}},`$
where $`z_{\mathrm{min}}T=E_\gamma +\sqrt{M^2+E_\gamma ^2}`$, with $`M=\mathrm{max}(m_a+m_b,m_1+m_2)`$, where the full squared amplitude is approximated as $`||^2`$ $``$ $`|_0|^2`$($`J^2`$), where $`_0`$ is the elastic scattering amplitude, $`J^\mu `$ is the hadronic current including the fundamental charge $`e`$, $`K_1`$ is the modified Bessel function of order 1, and finally, where the ratio of the two-body phase space for $`s_2=s2\sqrt{s}E_\gamma `$ and $`s`$ respectively, is a phase-space correction factor to minimize error at finite $`E_\gamma `$.

Bremsstrahlung reactions are shown in Fig. 3 for pions, rho mesons, $`K`$ and $`K^{}`$ mesons. Interactions are once again modeled with an SU(3) chiral Lagrangian with coupling constants fitted by respecting hadronic phenomenology. As expected, the partial rates are very strong at low photon energy. But the striking feature is the strength surviving even for intermediate photon energies.

Two cautionary remarks must however be made. First, Eq. (5) does not include multiple scattering interference effects of Landau-Pomeranchuk-Migdal, which will reduce the rates when the photon energy is less than the inverse of the mean time between strong interactions and second, the soft photon approximation becomes less reliable with increasing photon energy. Still, one does not expect more than a factor of 2–4 suppression for $`E_\gamma `$ $``$ 0.5 GeV.

Finally, the most useful comparison one can make regarding the relative importance of the higher-order contributions is to add up all the partial rates from radiative decays (order $`e^2`$), Compton processes (order $`e^2g^2`$), and from bremsstrahlung (order $`e^2g^4`$) and then compare to the conventional hadronic rate from the literature (which includes $`\omega \pi ^0\gamma `$, $`\pi \pi \rho \gamma `$, $`\rho \pi \pi \gamma `$ and $`\pi \rho \pi \gamma `$). This comparison is reported in Fig. 4.

A factor of 2 increase is seen over most of the intermediate and high energy region up to 3 GeV. Below 1 GeV photon energy, the increase is rather large, approaching two orders of magnitude around 200 MeV photon energy.

One might summarize this study as follows. Photon production rates in hot hadronic matter have been investigated up to next-to-next-to-leading order in the strong coupling. This means contributions to order $`e^2g^4`$ have been included, most notably, elastic scattering of mesons with accompanying bremsstrahlung. In addition, strangeness reactions were included and found to be somewhat important. Quantitatively, the new photon rate with strangeness and with NNLO hadronic processes is two orders of magnitude greater at low photon energies and a factor of approximately 2 greater at intermediate and high energies as compared with the standard photon production rate from the literaturenad .

Since multiple scattering interferences have not yet been included, the precise increase at low photon energies reported here in the new rate as compared with the old rate is not quantitatively determined. However, the general conclusions are expected to be robust. It will therefore become extremely important to assess such effects in future studies, as well, to incorporate these new rates into dynamical or statistical models for the purposes of comparison with experimental results at RHIC and interpretation regarding QGP diagnostics.

Acknowledgments. I am grateful to Charles Gale and Simon Turbide for useful discussions since they are carrying out a parallel investigation, emphasizing form factor effects, strangeness contributions and axial vector mesons.

This work has been supported in part by the National Science Foundation under grant number PHY-0098760.