text:
# New quantum phase transitions in the two-dimensional 𝐽₁-𝐽₂ model

\[

## Abstract

We analyze the phase diagram of the frustrated Heisenberg antiferromagnet, the $`J_1J_2`$ model, in two dimensions. Two quantum phase transitions in the model are already known: the second order transition from the Néel state to the spin liquid state at $`(J_2/J_1)_{c2}=0.38`$, and the first order transition from the spin liquid state to the collinear state at $`(J_2/J_1)_{c4}=0.60`$. We have found evidence for two new second order phase transitions: the transition from the spin columnar dimerized state to the state with plaquette type modulation at $`(J_2/J_1)_{c3}=0.50\pm 0.02`$, and the transition from the simple Néel state to the Néel state with spin columnar dimerization at $`(J_2/J_1)_{c1}=0.34\pm 0.04`$. We also present an independent calculation of $`(J_2/J_1)_{c2}=0.38`$ using a new approach.

\]

The nature of the quantum disordered phases of low-dimensional quantum antiferromagnets is a topic of fundamental importance for the physics of quantum magnetism . Such phases can result from mobile holes in an antiferromagnetic background as in the $`tJ`$ or Hubbard model at finite doping. Alternatively, competition of purely magnetic interactions can also lead to destruction of long-range order. A typical example of the second kind is the $`J_1J_2`$ model which exhibits a quantum disordered (spin-liquid) phase due to second-neighbor frustrating interactions. Even though it has been intensively studied during the last ten years, the $`J_1J_2`$ model apparently still holds many secrets. This model is also an ideal testing ground for the theory of quantum phase transitions because it has very complex dynamics and contains a variety of transitions. Exact diagonalization studies  have shown that the excitation spectrum of the model is quite complex and that finite-size effects are large . Spin-wave like expansions around the simple Néel state (which occurs for small frustration) naturally cannot give any information about the ground state at stronger frustration, and consequently non-perturbative methods are needed to analyze the latter regime.

An important insight into the disordered regime was achieved by field-theory methods  and dimer series expansions . The above works have established the range of the disordered regime, $`0.38<g<0.60`$ ($`g=J_2/J_1`$), and have also shown that the ground state in this regime is dominated by short-range singlet (dimer) formation in a given pattern (see Fig.1). The stability of such a configuration implies that the lattice symmetry is spontaneously broken and the ground state is four-fold degenerate. This picture is somewhat similar to the one dimensional situation, where the Lieb-Schultz-Mattis theorem guarantees that a gapped phase always breaks the translational symmetry and is doubly degenerate, whereas gapless excitations correspond to a unique ground state .

Two very recent calculations  performed by Green function Monte Carlo methods have raised new questions on the structure of the intermediate phase. The authors of Ref.  claim stability of the “plaquette RVB” state at $`g0.5`$. Reference  comes to a different conclusion: there is a columnar spin dimerized state with plaquette type modulation along the columns. An additional very interesting observation  is that the columnar spin dimerization penetrates into the Néel phase to $`g0.3`$. To conclude the list of observations which do not agree with a simple spin liquid with columnar dimerization we mention the divergence in the plaquette susceptibility found in Ref.  at $`g0.5`$.

In the present paper we elucidate all the above questions and come to the conclusion that two additional quantum critical points exist in the phase diagram of the system. These critical points correspond to a new generic type of second order quantum phase transition considered in Ref. . At each of the critical points points there is condensation of some singlet excitation and the critical dynamics is described by the nonlinear O(1) $`\sigma `$-model.

The Hamiltonian of the $`J_1J_2`$ model reads:

$$H=J_1\underset{nn}{}𝐒_i𝐒_j+J_2\underset{nnn}{}𝐒_i𝐒_j,$$
(1)
where $`J_1`$ is the nearest-neighbor, and $`J_2`$ is the frustrating next-nearest-neighbor Heisenberg exchange on a square lattice (see Fig.1). Both couplings are antiferromagnetic, i.e. $`J_{1,2}>0`$ and the spins $`S_i=1/2`$. We also use the notation $`g=J_2/J_1`$. The spin columnar dimerization at $`g>g_{c2}`$ is well established  and therefore we start our consideration from this state shown schematically in Fig.1. If there is an instability with respect to some kind of additional ordering then the gap in the spectrum of some singlet excitation must vanish at the corresponding critical point . We do not have a reliable technique for direct calculation of the singlet gap, but we do have a well developed series expansion technique for calculation of static susceptibilities. A static susceptibility is proportional to the corresponding Green’s function at zero frequency

$$\chi _𝐪G_𝐪(\omega =0)Z_𝐪/\omega _𝐪^2,$$
(2)
where $`\omega _𝐪`$ is the quasiparticle energy, and $`Z_𝐪`$ is the quasiparticle residue. So at the critical point $`1/\chi `$ must vanish approximately as $`(gg_c)^\gamma `$, with $`\gamma 2(\nu \eta )`$, where $`\nu `$ is the critical index for the spectral gap, $`\mathrm{\Delta }(gg_c)^\nu `$, and $`2\eta `$ is the critical index for the quasiparticle residue, $`Z(gg_c)^{2\eta }`$

To analyze possible plaquette type modulation we calculate the susceptibility of the spin columnar dimerized state with respect to the field

$$F_P=\underset{i,j}{}(1)^j𝐒_{i,j}𝐒_{i,j+1},$$
(3)
which breaks the translational symmetry in the direction perpendicular to the dimers. The series has been computed up to the seventh order in the dimerization parameter. Results for $`1/\chi _P`$ are shown in Fig. 2. The value of $`1/\chi _P`$ vanishes at $`g_{c3}=0.50\pm 0.02`$ and this is the critical point for the second order quantum phase transition from a simple columnar dimerized state to the 8-fold degenerate columnar dimerized state with plaquette type bond modulation in the direction perpendicular to the dimers suggested in Ref. . This phase transition is of the generic type considered in Ref.  and therefore it is described by 2D nonlinear O(1) $`\sigma `$-model. The critical indexes for this model are : $`\nu 0.630`$, $`\eta 0.034`$. Therefore one shall expect $`\gamma =2(\nu \eta )1.19`$. On the other hand the Dlog Padé approximants to the series $`\chi _P`$ give $`\gamma =0.9\pm 0.1`$. This is fair agreement, and we offer an explanation for the small discrepancy. The phase transition is related to the condensation of some singlet excitation which can be considered as a bound state of triplet excitations.

$$|s=a_2|tt+a_3|ttt+a_4|tttt+\mathrm{}$$
(4)
We would like to stress that there is very strong mixing between two-triplet and multi-triplet bound states. This mixing was the reason why vanishing of the singlet gap at $`g=g_{c3}`$ was missed in Reference . In Ref.  analysis of the singlet excitation was based on a two-particle Bethe-Salpeter equation with further account of multi-particle contributions as a small perturbation. This assumption was wrong because of the strong mixing. So at $`g=g_{c3}`$ we have condensation of effectively a multi-particle bound state with relatively small two particle component. The mixing between two-particle and multi-particle components of the singlet excitation varies with $`g`$ and this effect cannot be taken into account in the non-linear $`\sigma `$-model which assumes condensation of an “elementary” (=structureless) field. Ultimately very close to the critical point the variation of the mixing can be neglected and one shall expect restoration of the pure $`\sigma `$-model behavior. However it happens in so narrow vicinity of the critical point that the present numerical data cannot assess it.

Let us consider now the appearance of spin dimer order at $`g=g_{c1}`$ as $`g`$ is increased from small values. A scenario put forward some time ago  and based on the analysis of the $`Sp(N)`$, $`N\mathrm{}`$ field theory suggests that the dimer order appears simultaneously with disappearance of the Néel order, $`g_{c1}=g_{c2}`$. The dynamics in the vicinity of the critical point is described by the nonlinear O(3) $`\sigma `$-model in spite of an additional dimer order parameter. The additional gapless excitation is irrelevant to the critical dynamics because this excitation has extremely large size : $`r1/(gg_c)^M`$, $`M1`$. Another possibility is that $`g_{c1}<g_{c2}`$ and hence there are two separated quantum phase transitions . The transition at $`g_{c2}=0.38`$ is still described by the nonlinear O(3) $`\sigma `$-model, while the transition at $`g_{c1}`$ is of the O(1)$`\times `$O(1)-type. So in the vicinity of the point $`g_{c1}`$ there is an additional effective singlet field which can condense either at momentum $`𝐤=(\pi ,0)`$ or $`(0,\pi )`$. The sign degeneracy of the scalar field together with the momentum degeneracy gives a four-fold degenerate ground state which exactly corresponds to the degeneracy of the spin-dimerized state. A recent work based on the Green function Monte Carlo method  gives a hint in favor of this picture.

Let us give the precise meaning to the terms relevant and irrelevant singlet excitation. We consider a quantum critical point at which the singlet gap $`\mathrm{\Delta }_s`$ vanishes. An external field which is coupled to the singlet excitation, $`s|F|00`$, is applied to the system. If the corresponding susceptibility given by eq. (2) is diverging at the critical point we call this singlet excitation “relevant”. If the susceptibility is not diverging we call the singlet excitation “irrelevant”. It is clear that for an irrelevant excitation the quasiparticle residue $`Z`$ vanishes faster than $`\mathrm{\Delta }_s^2`$.

To analyze the problem of spin dimer order we calculate the susceptibility of the Néel state with respect to the external field which probes spin columnar dimerization.

$$F_D=\underset{i,j}{}(1)^i𝐒_{i,j}𝐒_{i+1,j}.$$
(5)
In this calculation we use the usual Ising series expansion up to seventh order. Note that in spite of the similarity between (3) and (5) these are two quite different situations. The field (3) assumes that the dimers aligned along the i-direction already exist and it probes a possible modulation in the j-direction. The field (5) is applied to the Néel state and therefore it does not assume any dimer order. The values of $`1/\chi _D`$ versus $`g`$ are plotted in Fig. 3. It is clear that $`1/\chi _D`$ vanishes somewhere in the interval $`g[0.3,0.4]`$, but the data is not precise enough to distinguish $`g_{c1}`$ from $`g_{c2}`$. To distinguish between the two scenarios discussed above we have to realize that Fig. 3 clearly indicates the relevant singlet excitation. In the case of an irrelevant singlet  the quasiparticle residue is extremely small, $`Z(gg_c)^M`$, $`M1`$ and hence the susceptibility has no divergence at the critical point. Thus we conclude from Fig. 3 that $`g_{c1}=0.34\pm 0.04`$, and that $`g_{c1}<g_{c2}`$, so there is a region $`g_{c1}<g<g_{c2}`$ where the spin columnar dimer order and the Néel order coexist. The critical dynamics at $`g_{c1}`$ is described by the relevant gapless singlet excitation. There is no doubt that the irrelevant gapless singlet excitation at $`gg_{c2}`$ also exists, but it has an exponentially small residue  and hence its contribution to the susceptibility is negligible.

The final result we report here is a new way of estimating $`g_{c2}0.38`$. The previous best calculation  was based on vanishing of the triplet gap in the spin liquid phase. A previous attempt to estimate $`g_{c2}`$ by Ising expansions for the staggered magnetization in the Néel phase showed the magnetization vanishing around 0.4, but the series were erratic in this region and the precision low. The new estimate is based on Ising expansions in the Néel phase for the 1st, 2nd and 3rd neighbor correlators $`S_i^xS_j^x`$ and $`S_i^zS_j^z`$, where $`z`$ is the direction of staggered magnetization. The series has been computed up to order 9 for 1st and 2nd neighbor correlators and to order 7 for 3rd neighbor correlator. The differences of these correlators are shown in Fig.4. The transition point is identified by the condition $`S_i^xS_j^x=S_i^zS_j^z`$, corresponding to restoration of spin rotational symmetry to the ground state. This gives $`g_{c2}0.38(3)`$, in excellent agreement with previous results.

In conclusion, the zero temperature phase diagram and the excitation spectra of the $`J_1J_2`$ model are shown schematically in Fig. 5. There are four critical points: $`g_{c1}=0.34\pm 0.04`$, $`g_{c2}=0.38`$, $`g_{c3}=0.50\pm 0.02`$, $`g_{c4}=0.60`$. The states are: $`g<g_{c1}`$ \- the simple Néel state, $`g_{c1}<g<g_{c2}`$ \- the columnar dimerized Néel state, $`g_{c2}<g<g_{c3}`$ \- the simple columnar dimerized spin liquid, $`g_{c3}<g<g_{c4}`$ \- the columnar dimerized spin liquid with plaquette type modulation, $`g>g_{c4}`$ \- the collinear state. The transitions at $`g_{c1}`$ and $`g_{c3}`$ are second order phase transitions of the O(1)$`\times `$O(1) and O(1) symmetry classes correspondingly. Energies of the relevant singlet excitations vanish at the critical points. The transition at $`g_{c2}`$ is a second order phase transition described by the nonlinear O(3) $`\sigma `$-model. The energy of the triplet excitation vanishes at $`gg_{c2}`$; at the critical point there is also a singlet excitation with zero gap, but this singlet is irrelevant. The transition at $`g_{c4}`$ is probably of first order, but is very close to second order.

This work forms part of a research project supported by a grant from the Australian Research Council. We thank V. N. Kotov and J. M. J. van Leeuwen for very helpful stimulating discussions. We are also grateful to S. Sachdev for very important comments. The computation has been performed on Silicon Graphics Power Challenge and Convex machines. We thank the New South Wales Centre for Parallel Computing for facilities and assistance with the calculations.


format: text

text:
# New results from the HRX BL Lac sample

## 1. Introduction

BL Lac objects are rare AGN, which are thought to be oriented towards us with their jet, thus showing high polarisation, strong variability at all wavelengths and non-thermal featureless spectra. Still there are several open questions about the nature of the BL Lac objects, e.g. whether there are differences between the BL Lac objects which are selected due to their strong radio-emission (RBL) and those selected because of their X-ray brightness (XBL). The difference between those two classes might be mainly caused by different peak frequencies of the two components in the BL Lac spectrum: XBL have usually a higher peak frequency of the synchrotron branch than RBL and thus we have high-frequency (HBL) and low frequency (LBL) peaked BL Lac objects (see e.g. Fossati et al. 1998). Between those two classes one can find intermediate objects (IBL) whose properties show them to be the link between HBL and LBL. We study the X-ray spectra of IBL and HBL and fit them with a single-power law with absorption by neutral hydrogen to investigate any possible spectral differences between these two samples.

## 2. The Hamburg X-ray bright BL Lac sample

This work is based on the bright sources from ROSAT All Sky Survey (Voges et al. 1996), which have an X-ray flux $`f_X>1110^{13}\mathrm{erg}\mathrm{cm}^2\mathrm{sec}^1`$ in the hard PSPC energy band ($`0.52.0`$ keV). A first complete sample of 39 BL Lac objects (Bade et al. 1998) in an area of $`2800\mathrm{deg}^2`$ is based on the Hamburg/RASS X-ray bright sample (HRX, Cordis et al. in preparation). A larger sample of candidates has been derived by cross-correlating the X-ray sources from RASS in a larger area ($`4500\mathrm{deg}^2`$) with radio catalogues (NVSS, FIRST: radio flux limit $`2.5\mathrm{mJy}`$). These 262 sources, if not already classified, have been included in follow-up spectroscopy with the 3.5m telescope on Calar Alto. Up to now this statistically complete sample is 95 % classified and contains 72 BL Lac objects, both of the HBL and the IBL type. For 61 of them we have already determined the redshift ($`z<0.9`$; Beckmann 1999). We will use in the following this sample of 72 BL Lacs.

## 3. X–ray spectra

Since the only available spectral measures are the two hardness ratios given in the RASS-BSC, we can only assume a simple spectral shape and apply the method of Schartel et al. (1996) to determine its parameters. We therefore assume that the spectrum is in the form of a single-power law absorbed by two different contributions: the galactic absorption and an “intrinsic” absorption that is given by the difference between the best fit absorption and the Galactic one ($`N_{H,intrinsic}=N_{H,freefitted}N_{H,galactic}`$). The galactic values were taken from Dickey and Lockman (1990). We find a strong correlation of the “intrinsic” absorption with the spectral slope (Fig. 1) on a $`>99\%`$ confidence level also when taking into account the large errors in $`N_{H,intrinsic}`$.

The negative values of $`N_{H,intrinsic}`$ seem to be caused by large errors; if we just take into account objects with an error of $`N_{H,intrinsic}<0.3`$, we do not have any “negative absorption” at all. Since we do not see any correlation to other observable parameters (such as flux or luminosity), we make the hypothesis that this is an effect of curvature in the X-ray spectra (see e.g. Urry et al. 1996). Thus the flat X-ray spectra are well described by a single power law with galactic low energy absorption ($`N_{H,intrinsic}0`$), while for steep spectra the curvature is significant and an additional “absorption” is needed to fit the X-ray data to a single-power law. In this model, a large value of $`N_{H,intrinsic}`$ would be explained by a convex spectrum. Nevertheless, true intrinsic absorption can not be ruled out in principle.

## 4. $`\alpha _{OX}`$ \- $`\alpha _X`$ relation for HRX-BL Lac

A relation, which is found for BL Lac objects, is the dependence of the spectral slope $`\alpha _X`$ on $`\alpha _{\mathrm{OX}}`$111we define $`\alpha _{\mathrm{OX}}`$ as the power law index between $`1`$ keV and $`4400`$ Å with $`f_\nu \nu ^{\alpha _{\mathrm{OX}}}`$, that describes the X-ray dominance over the optical brightness. This effect, which shows objects with a higher X-ray dominance having flatter X-ray spectra, was detected in several X-ray bright samples of BL Lac objects (e.g. Wolter et al. 1998, Padovani & Giommi 1995, Comastri et al. 1995).
Based on the RASS data we determined this relation for the HRX-BL sample (Fig.2). The effect of spectral flattening in the X-ray region with decreasing $`\alpha _{\mathrm{OX}}`$ is detectable, even though less significant than in previous works. An analysis, which takes into account the errors in $`\alpha _{\mathrm{OX}}`$ and $`\alpha _\mathrm{X}`$ gives a significance for correlation of both values on a $`>93\%`$ level.

## 5. Dependence of X-ray luminosity on $`\alpha _{OX}`$

Another physical parameter which seems to be correlated to the X-ray dominance is the X-ray luminosity of the BL Lac objects. Because we know the redshift and thus the luminosity for more than $`80\%`$ of our objects, it is possible to study this relation. Figure 3 shows the dependency of $`\mathrm{log}L_X`$ on $`\alpha _{OX}`$. There seem to be no objects, which are X-ray dominant (low $`\alpha _{OX}`$) and have a low X-ray luminosity. This result is solid, because we do not miss any redshift for objects with $`\alpha _{OX}<0.9`$ (Beckmann 1999). This effect could be explained in the view of the unified schemes (e.g. Padovani & Giommi 1995). The HBL are more X-ray dominated (low $`\alpha _{OX}`$) and the X-ray luminosity is therefore high because the X-ray band is in the vicinity of the peak frequency. On the other hand, LBL (high $`\alpha _{OX}`$) have their maximum of the synchrotron emission in the optical region and are X-ray faint.

## REFERENCES

Bade, Beckmann, Douglas, et al., 1998, A&A 334, 459
Beckmann, 1999, in: PASPC Vol. 159, eds. Takalo & Silanpää (astro-ph/9810154)
Comastri, Molendi & Ghisellini, 1995, MNRAS, 277, 297
Dickey and Lockman 1990, ARA&A 28, 215
Fossati, Maraschi, Celotti, et al., 1998, MNRAS 299, 433
Padovani and Giommi, 1995, MNRAS 277, 1477
Schartel, Walter, Fink, and Trümper 1996, A&A 307, 33
Urry, Sambruna, Worral et al. 1996, ApJ 463, 424
Voges, Aschenbach, Boller et al. 1996, IAU Circ. 6420
Wolter, Comastri, Ghisellini et al., 1998, A&A 335, 899


text:
# Order and Chaos in the Local Disc Stellar Kinematics

## 1. Introduction

The observed local stellar velocity distribution in the $`UV`$ plane, where $`U`$ is the radial velocity positive towards the anti-centre and $`V`$ the azimuthal velocity positive towards galactic rotation, betrays a large stream of old disc stars with an asymmetric drift $`S4050`$ km s-1 and $`\overline{U}>0`$ (e.g. Dehnen 1998), which hereafter will be referred to as the “Herculis” stream according to a co-moving Eggen group. One of the most reasonable interpretation of this stream relies on the influence of the Galactic bar, which is believed to have a corotation radius $`R_{\mathrm{CR}}3.55`$ kpc and a major-axis presently leading by $`\phi 15^{}45^{}`$ the Sun’s azimuth. Raboud et al. (1998) have shown that $`N`$-body models of the Milky Way do produce a mean outward motion of disc particles at realistic positions for the Sun relative to the bar and suggest that the Herculis stream involves stars on hot orbits, i.e. orbits with a high enough value of Jacobi’s integral to cross the corotation radius, whereas Dehnen (2000, hereafter D2000) relates this velocity mode and the main low-asymmetric drift mode of the observed $`UV`$ distribution to the coexistence near the outer Lindblad resonance (OLR) of the inner anti-bar and outer bar elongated periodic orbits replacing the circular orbit close to the OLR in a rotating barred potential (i.e. the same idea introduced by Kalnajs 1991 to explain the Hyades and Sirius streams), and attributes the valley at $`S30`$ km s-1 between the two modes (respectively “OLR” and “LSR” mode in his terminology) to stars on unstable OLR orbits. This paper analyses the role of bar induced chaotic orbits in this context and confronts both alternatives using a similar 2D modelling technique as in Dehnen’s work. The details are presented in Fux (2000).

## 2. Resonances and periodic orbits

The working potential is the same as in D2000, consisting of the sum of an axisymmetric logarithmic potential with constant circular velocity $`V_{}`$ and of a rotating quadrupole term with a bar strength $`F=0.15`$, where $`F`$ is defined as the maximum azimuthal force on the circle passing through the ends of the bar normalised by the axisymmetric radial force at same radius ($`F=8.89\alpha `$ in Dehnen’s notation). The position in space will be parametrised by the galactocentric distance relative to the OLR radius $`R_{\mathrm{OLR}}=(1+1/\sqrt{2})R_{\mathrm{CR}}`$ and the angle $`\phi `$ relative to the major axis of bar. A distance of $`R/R_{\mathrm{OLR}}=1.1`$ implies $`R_{\mathrm{CR}}=4.26`$ kpc if $`R=R_{}=8`$ kpc.

Figure 1 (left) shows the characteristic curves in the Hamiltonian$`x`$ plane (at $`y=0`$ and $`\dot{y}>0`$) of the main periodic orbits existing near the OLR and closing after one rotation in the rotating frame, using the labelling of Contopoulos & Grosbol (1989). The curves essentially resemble those of the circular orbit and of the resonant orbits in the axisymmetric limit, except near the resonances of the circular orbit where bifurcations may transform into gaps, as is the case at the OLR (e.g. Sellwood & Wilkinson 1993). Figure 1 also displays some members of the $`x_1(1)`$ and $`x_1(2)`$ families. The thick orbits are examples of the perpendicularly oriented orbits that replace the circular orbit near the OLR and which are responsible for the bimodality in the observed $`UV`$ distribution according to D2000. In both families, the orbits develop loops at high eccentricity, i.e. on the part of the characteristic curves emerging leftwards from the zero velocity curve and where $`\dot{x}0`$, and the $`x_1(1)`$ orbits are stable over a remarkable eccentricity range. Figure 2 shows the location in the $`UV`$ plane of the $`4/1`$, $`2/1`$ (OLR) and $`1/1`$ resonances in the axisymmetric part of the potential for different space positions, and of the periodic orbits passing through these positions in the total axisymmetric+bar potential. At $`\phi =90^{}`$ most resonant periodic orbits are stable, except for instance the $`x_1^{}(2)`$ orbit at $`R/R_{\mathrm{OLR}}>1.1`$, whereas at $`\phi =0`$, most of them are unstable.

## 3. Hot orbits

In a rotating barred potential, the value of the Hamiltonian $`H=\frac{1}{2}\dot{\stackrel{}{x}}^2+\mathrm{\Phi }_{\mathrm{eff}}(\stackrel{}{x})`$, where $`\mathrm{\Phi }_{\mathrm{eff}}(\stackrel{}{x})`$ is the effective potential, is the only classical integral of motion, known as Jacobi’s integral. Two critical values of this integral are related to the Lagrangian points $`L_{1/2}`$ and $`L_{4/5}`$, i.e. respectively the saddle points and maxima of $`\mathrm{\Phi }_{\mathrm{eff}}(\stackrel{}{x})`$. If $`H>H_{12}\mathrm{\Phi }_{\mathrm{eff}}(L_{1/2})`$, an orbit is susceptible to cross the corotation radius, defining the hot orbit category, and if $`H>H_{45}\mathrm{\Phi }_{\mathrm{eff}}(L_{4/5})`$, an orbit is no longer spatially confined by this integral.

If $`U`$ and $`V`$ are measured relative to the galactic centre, the contours of constant $`H`$ in the $`UV`$ plane are circles centred on $`(V,U)=(R\mathrm{\Omega }_\mathrm{P},0)`$, where $`\mathrm{\Omega }_\mathrm{P}`$ is the pattern speed of the bar, and of radius $`\sqrt{2(H\mathrm{\Phi }_{\mathrm{eff}})}`$. Figure 2 gives the $`H_{12}`$ and $`H_{45}`$ contours as a function of space position for the working potential. The hot orbits are those outside (low $`V`$) the former contour. At $`R=R_{\mathrm{OLR}}`$ and $`\phi =30^{}`$, the $`H_{12}`$ contour crosses the $`V`$-axis at $`S=V_{}V0.177V_{}=35`$ km s-1 for $`V_{}=200`$ km s-1, placing the Herculis stream in the hot orbit zone.

## 4. Regular and chaotic regions of phase-space

The next step after the periodic and hot orbits is to determine the extension in phase-space of the regular orbits trapped around the stable periodic orbits and of the chaotic regions. A convenient tool for this purpose are the Liapunov exponents, which describe the mean exponential rate of divergence of two trajectories with infinitesimally different initial conditions. The largest such exponent, $`\lambda _1`$, determines the degree of stochasticity of the orbits and vanishes for regular orbits. Figure 2 gives the divergence timescale $`t_\lambda 1/\lambda _1`$ of the orbits as a function of planar velocity and space position.

The Liapunov divergence timescale can be as low as a few orbital times for chaotic orbits. As expected, the stable periodic orbits always lie within regular regions and the unstable orbits within chaotic regions. Furthermore, the $`H_{12}`$ contour appears as the average transition from regular motion at low $`H`$ to chaotic motion at high $`H`$, and thus most hot orbits are chaotic. Due to the four-fold symmetry of the potential, the diagrams at $`\phi =0`$ and $`\phi =90^{}`$ are symmetric with respect to $`U=0`$, and more generally, diagrams at supplementary angles are anti-symmetric to each other in $`U`$. At $`\phi =90^{}`$, the $`2/1`$ and $`1/1`$ resonance curves are mainly embedded within broad arcs of regular orbits spaced by chaotic regions at large $`H`$, while at $`\phi =0`$, these regular arcs occur between the resonance curves. At intermediate angles, the regular and chaotic regions become offset from the resonance curves and the $`U`$-symmetry brakes. In particular, for $`\phi 30^{}`$ and for the displayed radial range, a prominent region of regular orbits trapped around stable eccentric $`x_1(1)`$ orbits extends down to negative $`U`$, bounded roughly by the OLR curve on the right and penetrating well inside the hot orbit zone, whereas the positive $`U`$ part of the OLR curve is surrounded by a wide chaotic region extending somewhat inside the $`H_{12}`$ contour and coinciding very well with the observed $`UV`$ location of the Herculis stream. Note that since 2D axisymmetric potentials are always integrable, all chaotic regions reported here owe uniquely to the bar.

A detailed investigation reveals that the regular orbits associated with the stable low eccentricity periodic orbits of the $`x_1(2)`$ family exist for $`R/R_{\mathrm{OLR}}<1.1`$ and over an angle range around $`\phi =90^{}`$ increasing as $`R`$ decreases. Elsewhere, unstable $`x_1^{}(2)`$ orbits have completely substituted to the regular $`x_1(2)`$ islands. In particular, there is no $`x_1(2)`$ streaming motion expected at $`R/R_{\mathrm{OLR}}=1.1`$ and $`\phi =30^{}`$ in the present working potential.

## 5. Phase-space crowding

Finally, we want to know how nature populates the available orbits. This is done by test particle simulations starting with exactly the same initial axisymmetric distribution functions $`f_{}`$ and growing the bar the same way (i.e. over $`2`$ bar rotations) as in Dehnen’s default simulations. However, we do not apply the claimed high-resolution and noise-free backward integration technique based on the conservation of the phase-space density along orbits. The reason is because much longer integration times than the few bar rotations adopted in D2000 are necessary to get phase-mixed and quasi-steady evolved distribution functions, and that for integration times larger than about $`10`$ bar rotations, the fine-grained $`UV`$ distributions resulting from this technique become highly fluctuating on small scales and even noisy in the chaotic regions and thus must be smoothed to yield the physical coarse-grained distributions. Therefore, we decided to use the standard forward integration technique, sampling $`f_{}`$ with $`N=10^6`$ particles, and to average the velocity distributions, derived within finite volumes of radius $`R/80`$, over $`10`$ bar rotations to improve the statistics and also to minimise the phase-mixing problem. The results are depicted in Figure 3. Here the distance parameters in $`f_{}`$ always scale as $`R`$ and the truly changing parameters are those of the bar. Hence frames at different $`R/R_{\mathrm{OLR}}`$ come from distinct simulations.

The $`UV`$ distributions seem to present the same symmetry properties in $`U`$ as described in the previous section for the Liapunov divergence timescale. This is obviously not the case in D2000, as a consequence of the too short integration. The high angular momentum peak of the velocity distributions also coincides better with the trace of the non-resonant $`x_1(1)`$ orbits at $`U=0`$. The regular low eccentricity $`x_1(2)`$ regions also produce a peak in the velocity distributions, which always lies inside the $`H_{12}`$ contour.

In the hot orbit zone, the chaotic regions appear to be more heavily crowded than the regular regions. This is because the phase-space density in regular regions remains roughly constant during the simulations, whereas there is a net migration of particles on hot chaotic orbits from the inner to the outer space regions, and because the chaotic and regular regions are decoupled from each other, i.e. chaotic orbits cannot visit the regions occupied by regular orbits and vice versa. Hence at large $`H`$, the velocity distributions are dominated by chaotic orbits which are forced to avoid less populated regular regions. In the parameter range considered here, these forbidden regions are essentially the broad regular orbit regions around the stable eccentric $`x_1(1)`$ orbits. In particular, at $`R/R_{\mathrm{OLR}}=1.1`$ and $`\phi =30^{}`$, this $`x_1(1)`$ region lies at $`U<0`$, as mentioned above, and thus the main overdensity of hot chaotic orbits occurs at positive $`U`$, providing a natural interpretation of the Herculis stream. The depression on the right of this stream may be attributed to the decline of the hot orbit population as $`HH_{12}`$. The high eccentricity $`x_1(1)`$ region is part of Dehnen’s valley, but represents regular OLR orbits, and the unstable $`x_1^{}(2)`$ orbit does not always fall within this valley (e.g. at $`R/R_{\mathrm{OLR}}=1.2`$ and $`\phi =30^{}`$).

Another property of the hot chaotic regions is that the contours of the velocity distribution follow the lines of constant Hamiltonian. This can be understood in terms of the Jeans theorem (see also the answer to Dehnen’s question): in two-dimensions, chaotic orbits have no other integral of motion than Jacobi’s integral and thus the distribution function in their neighbourhood depends only on this integral, i.e. on the value of the Hamiltonian.

## References

Contopoulos G., Grosbol P. 1989, A&AR 1, 261

Dehnen W. 1998, AJ 115, 2384

Dehnen W. 2000, AJ 119, 800 (D2000)

Fux R. 2000, A&A in preparation

Kalnajs A.J. 1991. In: Sundelius B. (ed.) Proc. Varberg Conf., Dynamics of Disc Galaxies. Göteborg, p. 323

Raboud D., Grenon M., Martinet L., Fux R., Udry S. 1998, A&A 335, L61

Sellwood J.A., Wilkinson A. 1993, Rep. Prog. Phys. 56, 173

## Questions

### Ortwin Gerhard:

The spiral arms outside the Milky Way’s bar region must modify the population seen in the $`UV`$ diagram. Could you say how?

### Roger Fux:

The effect of spiral arms is not yet very clear. As evolving transient features, they certainly introduce a time dependency of the $`UV`$ distributions and probably increase the amount of chaos. Since Jacobi’s integral is then no longer conserved, they may also introduce some diffusion from the regular regions to the chaotic regions and vice versa (and the same remark is true for massive compact objects like giant molecular clouds). $`N`$-body simulations displaying spiral arms reveal time averaged $`UV`$ distributions that are qualitatively similar to those of the test particle simulations presented here, but with an azimuthal phase shift apparently related to the twist of the potential well induced by the spiral arms outside the bar.

### Walter Dehnen:

How can stars stay on chaotic orbits, which due to their $`H`$ are not bound to the galaxy?

### Roger Fux:

The timescale for chaotic hot orbits with moderate $`H`$ to escape the galaxy is much larger than the orbital timescale and even the Hubble time (see for instance Kaufmann & Contopoulos 1996, A&A 309, 381). In particular, once phase-mixed, the distribution function in the neighbourhood of such orbits may be considered as almost steady, which is a necessary condition to apply the Jeans theorem.


text: # 1 Surface-Brightness Profile Parameters Of The Control Galaxies


text:
# 1 Function f(x) in the critical case 𝑝=0.899.

Holography in an Early Universe with Asymmetric Inflation

Elcio Abdalla and L. Alejandro Correa-Borbonet

Instituto de Física, Universidade de São Paulo,

C.P.66.318, CEP 05315-970, São Paulo, Brazil

## Abstract

We discuss the validity of the holographic principle in a $`(4+n)`$ dimensional universe in an asymmetric inflationary phase.

PACS number(s): 98.80.Cq 04.50.+h

Recently a new possible solution to the hierarchy problem has been proposed , where the fundamental Planck scale goes down to the level of the Tev region, provided that there are compact extra dimensions where gravitationals fields propagate. In such a case, the size of these extra dimensions should be in the submillimiter scale in order to conform to the phenomenological constraints. Indeed, using Gauss’s law in $`(4+n)`$ dimensions it is possible to relate the Planck scale in the $`(4+n)`$ dimensional theory($`M_{}`$) with that of the $`4`$-dimensional theory ($`M_{Pl}`$) by means of the size of the $`n`$ compactified dimensions ($`b_0`$). It is easy to find

$$M_{Pl}^2=(b_0)^nM_{}^{n+2}$$
(1)
Assuming the usual value $`M_{Pl}10^{19}`$ GeV and $`M_{}1`$ TeV, we find $`b_010^{17+\frac{30}{n}}`$ cm. For $`n=2`$ we obtain the (reasonable) value $`b_01`$ mm. We would thus find deviation from Newtonian gravitation in a range within experimental search in the near future.

There are several papers discussing the theoretical and phenomenological aspects of this idea . From the fundamental point of view, gravity can propagate freely around the whole set of dimensions in the submillimeter region, but the Standard Model (SM) fields get confined in the $`4`$-dimensional world. This is feasible if we consider that the SM fields are localized on a $`3`$-brane immersed in the higher-dimensional space, as in the case of D-branes occurring in type I/II string theory,. Thus a very attractive model of early universe has been proposed , and the underlying cosmology that occurs for small extra dimensions can be of wide interest.

In this model the dynamics of the dimensions as described by the Einstein equations is the central question. They are divided in a (so called wall) four dimensional world, and the $`n`$-internal dimensions. In concrete terms, the picture proposed is that of an early inflationary era followed by a long epoch where the scale factor of the brane-universe undergoes a slow contraction while the internal dimensions continue to expand towards their stabilized value. When the contraction ends the universe goes asymptotically to the radiation dominated behaviour as described by the FRW model. Although details of the dynamics underlying the transition are not known, a reasonable alternative explanation of the dynamics leading to the Cosmological standard model may arises.

On the other hand, we have witnessed a great activity on the holographic principle  which states that all information about physical processes in cosmology can be stored on the surface enclosing the physical universe. This idea has received strong support from another conjeture known as the AdS/CFT correspondence. The conjetured equivalence between Supergravity on D-dimensional AdS space and conformal field theory on the $`(D1)`$-dimensional boundary has been proven very useful to gain information about supergravity in the bulk from knowledge about its boundary. Moreover, the holographic principle has been extended to cosmology and there have been many attempts to apply different formulations of this principle to various cosmological models.

In this paper we will try to check the compatibility of the holographic principle with the cosmological model proposed by Arkani-Hamed et al . Specifically the cosmological version of the holographic principle proposed by Fischler and Susskind will be checked. Here we will restrict to the study of the contraction epoch with two extra dimensions, in agreement with the afore mentioned phenomenological constraints.

The evolution of the scale factors $`a(t)`$ corresponding to the $`3`$-dimensional physical (observable) space, and $`b(t)`$, corresponding to the internal dimensions, is described by the Einstein equations as given by

$$6H_a^2+2H_b^2+12H_aH_b=\frac{V}{M_{}^4b^2},$$
(2)
$$\frac{\ddot{b}}{b}+H_b^2+3H_bH_a=\frac{1}{M_{}^4b^2}(\frac{V}{2}\frac{b}{8}\frac{V}{b}),$$
(3)
$$\frac{\ddot{a}}{a}+2H_a^2+2H_bH_a=\frac{b}{8M_{}^4b^2}\frac{V}{b}.$$
(4)
where we have introduced the Hubble parameters $`H_a\dot{a}/a`$ and $`H_b\dot{b}/b`$ for the two scale factors, the overdot denoting the derivative with respect to $`t`$. Here, $`V`$ is the stabilizing potential, which is a function of the scale factor $`b`$. These quantities are in the socalled string (brane) frame. The reason is that the kinematics in this frame is automatically expressed in terms of the units measured by the observers which live on the wall.

This system has been obtained from the original Einstein equations after neglecting the energy density $`\rho `$ and pressure $`p`$ (see Eq($`57`$) in ). Using the COBE data we can prove that values of $`\rho `$ and $`p`$ are much smaller than the potential $`V`$. This potential, in the semiclassical limit, may be viewed as an expansion in inverse powers of the scale factor $`b(t)`$. Simplifying, it has been taken as being a monomial form, $`V=Wb^p`$, where $`W`$ is a dimensionful parameter.

These equations can be solved exactly. Defining

$$a=a_ie^{\alpha (t)}b=b_ie^{\beta (t)}$$
(5)
and going to the new time variable

$$d\tau =e^{3\alpha 2\beta }dt$$
(6)
the previous system, after some algebra, can be reduced to a functional constraint equation

$$4\alpha +\frac{4p}{4+p}\beta =C_1+C_2\tau $$
(7)
together with a simple first order differential equation,

$$X^2=\frac{\mathrm{\Delta }}{4}\omega e^X+\frac{3(4+p)^2}{32}C_2^2$$
(8)
where

$$X=6\alpha +(2p)\beta $$
(9)
and $`C_1`$, $`C_2`$ are integration constants. The other two parameters that appear in the system are $`\mathrm{\Delta }=88pp^2`$ , $`\omega =\frac{W}{M_{}^4b_i^{2+p}}`$. The prime denote the derivative with respect to $`\tau `$.

The solutions of the system take different form, controlled by whether $`C_2`$ vanishes or not, and by the value of $`\mathrm{\Delta }`$. It is found that the solutions can be classified in four types:

$$(1)\mathrm{\Delta }=0;(2)\mathrm{\Delta }>0,C_2=0;(3)\mathrm{\Delta }>0,C_2<0;(4)\mathrm{\Delta }<0,C_2<0$$
(10)
The critical value $`\mathrm{\Delta }=0`$ corresponds to $`p0.899`$. In this case the solution is

$$a=a_i(\frac{t}{t_i})^k,$$
(11)
$$b=b_i(\frac{t}{t_i})^l,$$
(12)
with $`k=0.1266`$ and $`l=0.69`$. Such values satisfy the two algebraic equations, namely

$$3k+2l=1,$$
(13)
and

$$3k^2+2l^2=1,$$
(14)
implying that in this era the solution describes a Kasner universe.

Applying the F-S criteria for anisotropic universes the following entropy/area relation is obtained, with the bar quantities defined in the Einstein frame,

$$S/A=\frac{\overline{\sigma }(\overline{R}_{H,\overline{a}})^3(\overline{R}_{H,\overline{b}})^2}{[(\overline{a}\overline{R}_{H,\overline{a}})^3(\overline{b}\overline{R}_{H,\overline{b}})^2]^{4/5}},$$
(15)
where the comoving size of the horizon is

$$\overline{R}_{H,\overline{a}}(\overline{t})=_{\overline{t}_i}^{\overline{t}}\frac{d\overline{t}^{}}{\overline{a}(\overline{t}^{})},$$
(16)
and $`\overline{\sigma }`$ is the comoving entropy density.

In order to calculate the S/A relation $`(\text{15})`$ we use the expression that relate the metric in the Einstein frame with the one in the string frame

$$\overline{g}_{\mu \nu }=b^2g_{\mu \nu }.$$
(17)
Using this map, we get

$$\overline{a}=ba\overline{b}=bbd\overline{t}=bdt.$$
(18)
Therefore we have

$$\overline{R}_{H,\overline{a}}=R_{H,a}\overline{R}_{H,\overline{b}}=R_{H,b}.$$
(19)
Now, the S/A relation can be rewritten as

$$S/A=\frac{1}{b^4}\frac{\sigma (R_{H,a})^3(R_{H,b})^2}{[(aR_{H,a})^3(bR_{H,b})^2]^{4/5}}.$$
(20)
The authors of Ref  argued that when the Fischler-Susskind criteria is applied to inflationary cosmology the evolution of the universe should be considered after reheating. For the model of Arkani-Hamed et $`al`$ we will take the initial time $`t_i`$ as corresponding to the end of the de Sitter phase, because the model does not have the usual post-inflationary reheating of standard inflationary cosmology. We also should remember that the end of the de Sitter phase is the beginning of the contraction epoch.

Taking into account the scale factors $`(\text{11})`$ and $`(\text{12})`$ we get, after some manipulations,

$$S/A=\frac{\sigma }{a_i^3b_i^6}c_1t_if(x),$$
(21)
where

$$c_1=\frac{1}{(1k)^{3/5}}\frac{1}{(1l)^{2/5}}=1,4873$$
(22)
and

$$f(x)=\frac{1}{x^{4l}}\frac{(x^{1k}1)^3(x^{1l}1)^2}{(xx^k)^{12/5}(xx^l)^{8/5}},$$
(23)
where $`x=t/t_i`$.

The shape of $`f(x)`$ is illustrated in figure 1. Using $`(\text{13})`$ it can be checked that the asymptotic value of $`f(x)`$ is zero. Therefore if the holographic bound holds when $`f(x)`$ has its maximum value, later on, the bound will be satisfied even better.

The solutions in the case $`\mathrm{\Delta }>0`$($`0p<0.899`$) and $`C_2=0`$ are similar but the corresponding powers are

$$v=\frac{2p}{(2+p)(4+p)},$$
(24)
$$u=\frac{2}{2+p}.$$
(25)
In this case the powers satisfy

$$3v+2u>1.$$
(26)
Repeating the same procedure as the previous case we find for the entropy to area ratio the result

$$S/A=\frac{\sigma }{a_i^3b_i^6}c_2t_ig(x),$$
(27)
where for $`p<0.899`$

$$c_2=\frac{1}{(1v)^{3/5}}\frac{1}{(1u)^{2/5}},$$
(28)
$$g(x)=\frac{1}{x^{4u}}\frac{(x^{1v}1)^3(x^{1u}1)^2}{(xx^v)^{12/5}(xx^u)^{8/5}},$$
(29)
and for $`p=0`$ ,

$$c_2=1$$
(30)
$$g(x)=\frac{1}{x^4}\frac{(x1)^{3/5}(ln(x))^{2/5}}{x^{8/5}}.$$
(31)
In figure 2 the function $`g(x)`$ is shown for diferents values of the parameter $`p`$($`0,0.5,0.7`$). Also in this case the asymptotic value of $`g(x)`$ is zero. In figure 3 a comparation between the functions $`f(x)`$ and $`g(x)`$ is displayed.

In the case ($`\mathrm{\Delta }<0(p>0.899),C_2<0`$)the Kasner solution applies asymptotically (see Ref ). In other words, the scale factors $`a(t)`$ and $`b(t)`$ take the form $`(\text{11}),(\text{12})`$ for large times. This implies that the ratio entropy/area will coincide with the $`\mathrm{\Delta }=0`$ case in the remote future (see figure 4).

The last solution ($`\mathrm{\Delta }>0,C_2<0`$) is similarly reduced to the case ($`\mathrm{\Delta }>0,C_2=0`$). Then, the scale factors will have the powers $`(\text{24}),(\text{25})`$ and the entropy/area relation will be described by equation $`(\text{27})`$ in the remote future.

In conclusion, we have found that the holographic principle can be satisfied in the contraction epoch for the model of Arkani Hamed et $`al`$. The entropy to area ratio has a similar asymptotic behavior for all the solutions. As this model presents some novelties for inflationary cosmology, the results that we have obtained could be useful in order to improve our understanding of the holographic principle in relation to cosmology.

ACKNOWLEDGMENT: This work was partially supported by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPQ) and CAPES (Ministry of Education, Brazil).


text:
# Nonequilibrium molecular dynamics simulation of rapid directional solidification.

## I Introduction

Directional solidification (DS) of binary alloys is a reference experimental method to conduct carefully controlled tests of industrial casting. Besides their practical interest, DS experiments also bring insight into the fundamental study of basic instability morphologies of solid-liquid interfaces, such as cells or dendrites . By the combination of theoretical and numerical methods, much progress has been now achieved toward solving this difficult physical problem. However, most of the theoretical effort has concentrated so far on continuous descriptions of the underlying phenomena . In contrast, attempts to attack the problem at the atomistic level remain very few. To this extent, the present study may be considered as a first step to bridge the gap between both points of view.

Because of the micrometer size of the growth structures, a direct and quantitative atomistic simulation is still not at hand on this scale. Nevertheless, in the limit of large pulling velocities, a microscopic technique like Molecular Dynamics (MD) can still be used to follow atomistic phenomena occuring close to the solid-liquid interface. On one hand, MD was used to simulate laser-pulsed melting, for which the velocity is not controlled, but rather governed by heat diffusion. . On the other hand, a first MD simulation of directed growth has been recently reported by Coura and al. . They considered the growth of a solid from a fluid phase with a density ten times smaller, a case which is probably more relevant to deposition from a vapor.

In this paper, we present non-equilibrium molecular dynamic simulations of directional solidification in two dimensions. We restrict ourselves to the case of rapid solidification with a large temperature gradient, for which the interface remains plane on the atomistic scale, so that microstructures will not be considered here. After the description of the simulation details in section II, we present in section III the results obtained for the segregation profiles, the segregation coefficient and the advection velocity. In the last section, we finally discuss how this nonequilibrium simulation can be extended in the future to study different microscopic mechanisms involved in DS.

## II Simulation details

Two atoms $`i`$ and $`j`$ separated by a distance $`r`$ interact via the well-known Lennard-Jones potential :

$$u_{ij}(r)=4ϵ_{ij}[(\frac{\sigma _{ij}}{r})^{12}(\frac{\sigma _{ij}}{r})^6]$$
(1)
The interaction energies between pairs of solvent-solvent, solute-solute and solvent-solute atoms are respectively denoted by $`u_{11}`$, $`u_{22}`$ and $`u_{12}`$. The solvent potential parameters are chosen in order to describe Argon properties ($`ϵ_{11}=120`$K, $`\sigma _{11}=3.405\mathrm{\AA }`$). For the solute, we take $`ϵ_{22}=0.5ϵ_{11}`$ and $`\sigma _{22}=\sigma _{11}`$. The cross-species parameters are fixed using the Lorentz-Berthelot rules  :

$$ϵ_{12}=\sqrt{ϵ_{11}ϵ_{22}}$$
(2)
and

$$\sigma _{12}=(\sigma _{11}+\sigma _{22})/2.$$
(3)
All the interactions are truncated at a cut-off radius $`r_c=2.5\sigma _{11}`$ and the equations of motion are integrated using the Beeman algorithm  with a time step $`\delta t=0.01`$ps. The particles coordinates are defined in a reference frame moving at the pulling velocity $`V`$ in the $`x`$ direction and periodic boundary conditions (PBC) are applied in the $`x`$ and $`y`$ directions. After integration of the dynamical equations over a time $`\delta t`$, pulling is implemented by adding an increment $`V\delta t`$ to the $`x`$ coordinate of each atom. In the reference frame, the solid-liquid interface is thus immobile when the stationnary state is reached.

To simulate heat transport from the furnace to the system, four regions of fixed temperature are used (Fig. 1). Regions I, II, III and IV are centered at fixed positions, $`x_I`$, $`x_{II}`$, $`x_{III}`$ and $`x_{IV}`$ and have a width of $`20\mathrm{\AA }`$ for regions I and IV and $`10\mathrm{\AA }`$ for regions II and III. In each region the temperature is kept constant by using a classical velocity rescaling. To maintain the solidification front between region II and III, we impose $`T_I=T_{II}<T_m`$ and $`T_{III}=T_{IV}>T_m`$, $`T_m`$ being the melting temperature of the alloy.

In Fig. 1 we represent the temperature gradient obtained after equilibration in the simulation box. The large difference between $`T_{II}`$ and $`T_{III}`$ and hence the large gradient permits to localize the interface easily. With such a high gradient, instabilities cannot develop, which ensures the stability of planar interfaces. For the same reason we also take a small width $`L_y`$ for the simulation box in order to reduce the natural roughness of the interface. More realistic systems, as compared to experiments, would correspond to $`x_I=x_{II}`$ and $`x_{III}=x_{IV}`$, together with a larger value of $`L_y`$. Because of the PBC in the $`x`$ direction, we simultaneously have a solidification and a melting front, as in a melting zone experiment. We concentrate here on the solidification front and we use a wide liquid zone ($`x_{III}<x<x_{IV}`$) to allow solute diffusion. Fixing the temperature in four regions instead of two reduces considerably the rescaling of velocities within the liquid region. This is helpful to suppress pertubations and artefacts during the computation of microscopic quantities. The density difference between the solid and the liquid ($`20\%`$) is sufficient to induce advection of the liquid towards the front. A part of the associated momentum is then transmitted to the solid layer which in turn acquires a translationnal motion in the $`x`$ direction. To avoid this finite-size effect, we rescale to zero the mean velocity within the deeper part of the solid (region I).

Our system contains about $`2000`$ atoms, 10 percent of which are solute atoms. Its size, $`L_x=400\mathrm{\AA }`$ in length and $`L_y=60\mathrm{\AA }`$ in width, is relatively modest as compared to the size of systems currently used in MD. The reason is that the computational effort is here essentially spent in the time length of the simulation. To obtain good statistics, the simulation has to be long enough to allow each atom to perform several solidification-melting cycles. Since ten cycles require a time of $`10L_x/V`$, for the slowest velocity studied here ($`V=10`$cms-1), this represents a simulation time of $`4\times 10^4`$s ($`4\times 10^8`$ MD steps). To increase the performance of our code we then adopt the ’cell lists’ method  in which the box is divided into cells with a size slighty larger than the cutoff radius $`r_c`$.

## III Results

### A Concentration profiles

The partition ratio, or segregation coefficient, is defined as

$$k=c_s^i/c_l^i,$$
(4)
$`c_l^i`$ and $`c_s^i`$ being the concentrations at the interface, respectively in the liquid and the solid. With our choice of the potential parameters ($`ϵ_{12}<ϵ_{11}`$), $`k`$ is expected to be less than one , so that solute accumulates in the liquid near the front. We first illustrate how the simulation reproduces solute rejection. The pulling velocity is fixed to $`V=1`$ms-1 and the solute atoms are initially placed at random in the simulation box. Fig. 2a is a snapshot of the system in its initial configuration. After a time $`\mathrm{\Delta }t=5`$ns, which corresponds to a spatial translation of the furnace of approximatly $`L_x/8`$, the second snapshot (Fig 2b) shows that almost all the solvent atoms initially close to the front are incorporated into the solid. Conversely, because of their poorer solubility in the solid phase, a majority of the solute atoms remain in the liquid phase. As compared to a solvent atom, it takes a longer time for a solute atom to cross the interface and hence the solute concentration is higher at the interface, as expected here. To quantify this segregation, we compute the profile concentration by averaging over $`10^4`$ uncorrelated spatial configurations (Fig. 3). The diffusion equation for the solute in the moving frame reads

$$\frac{c}{t}=D\frac{^2c}{x^2}+V\frac{c}{x}.$$
(5)
In the stationary regime, it gives the well known theoretical profile for the solute concentration in the liquid,

$$c(x)=c_s^i+(c_l^ic_s^i)\mathrm{exp}(x/ł_s)$$
(6)
where $`ł_s=D/V`$ is the diffusion length and the origin of the $`x`$ axis is placed at the front. A good agreement is found between this expression and our simulations, the best fit giving a diffusion coefficient $`D_{fit}=0.3\mathrm{\AA }^2\mathrm{ps}^1`$. To verify this value, we independently compute the diffusion coefficient profile $`D(x)`$, with the help of Einstein’s relation between the mean square displacement of the atoms and $`D`$,

$$<|𝐫(0)𝐫(t)|^2>=4Dt.$$
(7)
In Fig. 4, we see that the diffusion coefficient increases abruptly between zero in the solid phase to roughly $`0.4\mathrm{\AA }^2\mathrm{ps}^1`$ in the liquid near the interface. This variation of $`D`$ takes place over a distance of $`30\mathrm{\AA }`$ that can be considered as a good aproximation of the interface width. In Eq. (5), $`D`$ was assumed constant. We can check that the fitted value $`D_{fit}`$ approximately corresponds to the average of $`D(x)`$ within the liquid region near the interface. The same analysis is repeated for different pulling velocities between $`V=0.1`$ms-1 and $`V=5`$ms-1. For each velocity, we compute the concentration profile and extract $`D_{fit}`$ from the diffusion length. In Fig. 5, $`D_{fit}`$ is plotted as a function of the pulling velocity. As discussed below, the segregation coefficient tends to unity for the largest $`V`$ and the calculation of $`D_{fit}`$ is thus restricted to velocities $`V<3`$ms-1. A reasonable agreement is found in each case between $`D_{fit}`$ and the average value of the $`D(x)`$ profile.

In Fig 6, we plot the two concentration profiles obtained for $`V=1`$ms-1 and $`V=0.5`$ms-1. As expected, the smaller the velocity, the larger $`ł_s`$. If the agreement with the exponential law is good for the largest velocity, we notice that it is less satisfactory for the smaller one. This can be understood by a second look at Fig. 4. We observe that the thermal gradient induces an increase of the diffusion coefficient, as we go deeper in the liquid. As a consequence, we also have an increase of the diffusion length, which explains the disagreement observed for the largest $`x`$ values.

### B Segregation coefficient

Now that we verified the ability of our method to reproduce both the segregation at the interface and the dependence of the diffusion length on the pulling velocity, we examine the behavior of the segregation coefficient $`k`$ as a function of $`V`$. We see in Fig. 6 that the liquid concentration at the interface, $`c_l^i`$, is larger for a lower value of $`V`$. This can be easily explained by the fact that, at large pulling velocities, the solute does not have enough time to fully diffuse over the interface width before it is incorporated into the solid. As a consequence, some solute is trapped in the solid and segregation at the interface is lowered. Several models have been proposed to describe solute-trapping . The usual qualitative criterion for segregation is to compare two different characteristic times. One, $`t_V=\delta _i/V`$, corresponds to the growth of the solid over a distance comparable to the interface width, $`\delta _i`$. The second, $`t_D=\delta _i^2/D`$, is the time needed for a solute atom to diffuse over the same distance $`\delta _i`$. The Aziz model  predicts then

$$k=\frac{k_e+\beta }{1+\beta }$$
(8)
for the segregation coefficient. Here $`k_e`$ is the segregation coefficient at equilibrium (i.e., when $`V=0`$) and

$$\beta =t_D/t_V=\delta _i/ł_s$$
(9)
is the ratio of the two characteristic times (or lengths) defined above. Two remarks arise at this point. The first one concerns the neighborhood of the absolute stability threshold, $`V=V_a`$, where the diffusion length becomes comparable to the capillary length, $`l_sd_0`$. When $`d_0\delta _i`$, one also has $`\beta 1`$, so that $`kk_e`$. In this case, the absolute stability threshold is well described by the usual linear stability analysis (see ). This is not true anymore, when $`d_0`$ becomes comparable to $`\delta _i`$, as it is very likely the case for some reference alloys . The second remark is that the assumption of continuous growth, on which Eq. (8) relies, is actually verified by the rough solid-liquid interfaces produced with our Lennard-Jones potential.

In Fig 7 we plot our estimates for $`k`$ as a function of the pulling velocity. The line is a best fit to the Aziz model. The parameter $`\beta `$ is estimated by using $`\delta _i=30\mathrm{\AA }`$, the characteristic width of the interface extracted from the $`D(x)`$ profile in Fig.4 and the value of the diffusion coefficient, $`D=0.2\mathrm{\AA }^2`$ps-1 is estimated at the center of this interface region. The agreement between the simulated values and the model is rather good for a fit with a single free parameter. We obtain the value $`k_e=0.503`$ which is difficult to confirm because it is not possible to simulate segregation at much lower velocities with current computing power. Simulations at a fixed pulling velocity thus provide an indirect way of determining $`k_e`$. The fact that the interface is moving permits all atoms, and especially the solute ones, to be part of the fluid phase and then to speed up relaxation to the equilibrium situation. For an immobile interface, one would have to wait for the very slow diffusion of solute atoms within the solid phase. The solid-liquid equilibrium for 3D binary Lennard-Jones mixtures has been determined by Monte-Carlo simulations . It would be interesting to have similar 2D results available to test our indirect evaluation of $`k_e`$.

### C Advection velocity

We finally study the advection of the fluid near the solidification front. The Lennard-Jones mixture considered in this work produces a substantial density difference between the solid and the liquid. This is illustrated in Fig. 8, where we plot the density profile along the $`x`$ direction perpendicular to the interface. From these data, we estimate both densities at the interface, $`\rho _s^i`$ and $`\rho _l^i`$. We find that they do not depend on the pulling velocity and that the normalized density difference is also a constant, $`(\rho _s^i\rho _l^i)/\rho _l^i0.24`$. This rather high value for a Lennard-Jones potential is due to the large temperature gradient imposed at the interface. As a consequence, the solid cannot grow unless the liquid is advected toward the front with a velocity $`V_{ad}`$. In the laboratory frame, the solid is at rest, while the liquid has an advection velocity $`V_{ad}`$ all along the $`x`$ axix (Fig. 9). Mass conservation at the interface imposes that $`V_{ad}^i`$ varies with $`V`$ as

$$V_{ad}^i=\frac{\rho _s^i\rho _l^i}{\rho _l^i}V.$$
(10)
A linear variation is effectively recovered in our simulations (Fig. 10), with a slope comparing well to the estimate extracted from the density profile (see Fig. 8). It is interesting to note that there is no deviation from the linear law, even at the highest velocities studied in this work. On one side, one could have expected that the onset of solute trapping would modify the densities. This is not the case here, probably because the solvent and solute atoms have the same atomic radius, $`\sigma _{22}=\sigma _{11}`$. On the other side, a deviation would also appear as $`V`$ is increased, if the solid was created with more and more defects. We would then have a decrease of $`V_{ad}^i`$ because the density of the solid would also decrease. This type of deviation is not observed any more in our simulations. We finally see in Fig. 9 that the advection velocity increases in magnitude with $`x`$. As we go deeper in the liquid, the density progressively decreases in the thermal gradient and an extra advection velocity is thus necessary to compensate the density drop between two adjacent slices of alloy.

## IV Conclusion

With the help of MD simulations, we directly observed and analyzed quantitatively two important phenomena occuring at the atomistic level in rapid DS. The first one is the cross-over between the regime of solute segregation and the regime of solute trapping, which sets in as the pulling velocity, $`V`$, is increased. The good quantitative agreement of our estimates for the segregation coefficient, $`k(V)`$, with existing theories  confirmed that MD is a valuable simulation tool for DS. The second point is the advection of the liquid phase towards the solidification front. The advection velocity at the interface was shown to remain proportional to $`V`$, in the whole range of values considered here. This effect is more likely masked by stronger phenoma like convection in the usual experimental set-ups. It may however become relevant in micro-gravity or for thin sample DS.

The use of Lennard-Jones potentials is rather restrictive to compare the present results with experiments. In the future, potentials based on empirical descriptions, like the embedded atom method , the glue model , or the effective medium theory  will be used to study technical materials like metals. This will allow us to explore important issues like anisotropy and facetting, or the role of defects and constraints. In parallel, it will be necessary to simulate more extended systems, in order to explore the front instabilities from an atomistic point of view. In this last category, we can cite the oscillations of the solidification front just above the absolute stability velocity .

###### Acknowledgements.

It is a pleasure for us to thank B. Billia, J. Jamgotchian, and H. Nguyen Thi for valuable discussions, as well as K. Kassner for a critical reading of the manuscript.


text:
# 1 a

## Acknowledgments

It is a pleasure to thank N. Armesto, J. A. Casado, E. G. Ferreiro, U. Heinz, A. B. Kaidalov, C. Pajares and J. Tran Thanh Van for discussions. A.C. acknowledges partial support from a NATO grant OUTR.LG 971390. C.A.S. thanks Ministerio de Educación y Cultura of Spain for financial support.


text:
# Untitled Document

MAGNETIC FIELD OF RELATIVISTIC NONLINEAR PLASMA WAVE

Arsen G. Khachatryan

Yerevan Physics Institute, Alikhanian Brothers Street 2,

Yerevan 375036, Armenia

e-mail: khachatr@moon.yerphi.am

Longitudinal and transverse behavior of magnetic field of relativistic nonlinear three-dimensional plasma wave is investigated. It is shown that the magnetic field of the wave is different from zero and performs higher frequency oscillations compared to the plasma electron frequency. An increase in the nonlinearity leads to strengthening of magnetic field. The oscillations of magnetic field in transverse direction arise, that caused by the phase front curving of nonlinear plasma wave. The numerical results well conform with predictions of the analytical consideration of weakly-nonlinear case.

PACS number(s): 52.35.Mw, 52.40.Mj, 52.40.Nk

The progress in the technology of ultrahigh intensity lasers and high current relativistic charged bunch sources permits the use of laser pulses1 or charged bunches2 for excitation of strong plasma waves. The excited plasma waves can be used, for example, for acceleration of charged particles and focusing of bunches.2 The amplitude of longitudinal electric field $`E_{\mathrm{max}}`$ of relativistic plasma waves excited in cold plasma is limited by the relativistic wave-breaking field3 $`E_{rel}=[2(\gamma 1)]^{1/2}E_{WB}/\beta `$, where $`\gamma =(1\beta ^2)^{1/2}`$ is a relativistic factor, $`\beta =v_{ph}/c`$ is a dimensionless phase velocity of the wave, $`E_{WB}=m_e\omega _{pe}v_{ph}/e`$ ($`E_{WB}[V/cm]0.96n_p^{1/2}[cm^3]`$) is the conventional nonrelativistic wave-breaking field, $`\omega _{pe}=(4\pi n_pe^2/m_e)^{1/2}`$ is the electron plasma frequency,  $`n_p`$ is the equilibrium density of plasma electrons, $`m_e`$ and $`e`$ are the mass and absolute value of the electron charge.

The linear plasma wave theory is valid when $`E_{\mathrm{max}}E_{WB}`$. It is well known that a three-dimensional linear plasma wave in a cold plasma and in the absence of external fields is potential.4 In this case the magnetic field of plasma wave is zero. The magnetic field is absent also in one-dimensional nonlinear case due to the symmetry of the problem. However, magnetic field of nonlinear three-dimensional plasma wave with the relativistic phase velocity ($`\beta 1`$) was not studied up to now and our aim is to investigate this problem.

We shall study nonlinear plasma waves (wake waves) excited in cold plasma by relativistic electron bunches or intense laser pulses (drivers) and suppose the azimuthal symmetry of the problem. In this case for non-zero components of plasma electrons momentum and electromagnetic field of the wave we have the following set of equations5,6:

$$\beta \frac{P_z}{z}\frac{\gamma _e}{z}\beta ^2E_z=0,$$
(1)
$$\beta \frac{P_r}{z}\frac{\gamma _e}{r}\beta ^2E_r=0,$$
(2)
$$\frac{H_\theta }{z}+\beta \frac{E_r}{z}+\beta _rN_e=0,$$
(3)
$$_{}H_\theta +\beta \frac{E_z}{z}+\beta _zN_e+\beta \alpha =0,$$
(4)
$$\beta \frac{H_\theta }{z}\frac{E_r}{z}+\frac{E_z}{r}=0,$$
(5)
$$N_e=1\alpha _{}E_r\frac{E_z}{z}.$$
(6)
As usual, the Eqs. (1) and (2) were derived taking into account that the curl of the generalized momentum is zero, $`\beta ^2𝐇rot𝐏=0`$, or in our case

$$\beta ^2H_\theta +\frac{P_z}{r}\frac{P_r}{z}=0.$$
(7)
In Eqs. (1) - (6) $`\gamma _e=(1+P_z^2+P_r^2+a^2/2)^{1/2}`$ , $`\beta _{z,\text{ }r}=P_{z,\text{ }r}`$ $`/\gamma _e`$ and $`N_e=n_e/n_p`$ are respectively a relativistic factor, dimensionless components of velocity and dimensionless density of plasma electrons, $`z=(\omega _{pe}/v_{ph})(Zv_{ph}t),`$ $`\alpha =n_b(z,r)/n_p`$, $`n_b`$ is density of bunch electrons, $`a=eE_0(z,r)/m_ec\omega _0`$, where $`E_0`$ and $`\omega _0`$ are the amplitude and frequency of laser pulse, $`_{}=/r+1/r`$. Also the following dimensionless variables have been used: the space variables are normalized on $`\lambda _p/2\pi =v_{ph}/\omega _{pe}`$, where $`\lambda _p`$ is the linear plasma wavelength, the momenta and velocities - respectively on $`m_ec`$ and the velocity of light and the strengths of electric and magnetic fields - on the nonrelativistic wave-breaking field $`E_{WB}`$. In the general case the analytical consideration of the problem seems impossible. First we consider weakly-nonlinear case and then present numerical results.

In the weakly-nonlinear case we suppose $`\beta =1`$ and use the generalization of the well known expansion7 which was used to study one-dimensional nonlinear relativistic-plasma-waves8:

$`u(r,z)`$ $`=`$ $`\epsilon u_1(r,\mathrm{\Psi })+\epsilon ^2u_2(r,\mathrm{\Psi })+\epsilon ^3u_3(r,\mathrm{\Psi })+\mathrm{},`$ (8)
$`\mathrm{\Psi }/z`$ $`=`$ $`1+\epsilon k_1(r)+\epsilon ^2k_2(r)+\mathrm{},`$
where $`u`$ stands for normalized values $`P_{z,r}`$, $`E_{z,r}`$ or $`H_\theta `$, $`\epsilon =E_z^{\mathrm{max}}1`$ is the small parameter, $`P_{z,r}`$, $`E_{z,r}`$, $`H_\theta 1`$, $`\mathrm{\Lambda }_p`$ is the nonlinear plasma wavelength. Substituting expansion (8) into equations (1)-(5) we have:

$$\frac{P_{zi}}{\mathrm{\Psi }}E_{zi}=S_{1,i},$$
(9)
$$\frac{P_{ri}}{\mathrm{\Psi }}E_{ri}=S_{2,i}$$
(10)
$$\frac{H_{\theta i}}{\mathrm{\Psi }}+\frac{E_{ri}}{\mathrm{\Psi }}+P_{ri}=S_{3,i},$$
(11)
$$_{}H_{\theta i}+\frac{E_{zi}}{\mathrm{\Psi }}+P_{zi}=S_{4,i},$$
(12)
$$\frac{H_{\theta i}}{\mathrm{\Psi }}\frac{E_{ri}}{\mathrm{\Psi }}+\frac{E_{zi}}{\mathrm{\Psi }}=S_{5,i},$$
(13)
where subscript $`i`$ denotes the order of approximation and $`S_{15,i}`$ are the nonlinear functions of $`u_{i1}`$, $`u_{i2}`$, $`\mathrm{}`$, and $`u_1`$. In the linear case ($`i=1`$) $`S_{15,1}=0`$ and one can obtain well known solutions (see, e.g., Ref. 9): $`P_{z1}=R\mathrm{sin}\mathrm{\Psi },`$ $`P_{r1}=(dR/dr)\mathrm{cos}\mathrm{\Psi },`$ $`E_{z1}=R\mathrm{cos}\mathrm{\Psi },`$ $`E_{r1}=(dR/dr)\mathrm{sin}\mathrm{\Psi },`$ and $`H_{\theta 1}=0`$, where $`R(r)`$ depends on the radial distribution in the exciting source, i. e., on $`\alpha (r)`$ or $`a^2(r)`$. In the second approximation, from the condition of absence of resonance terms (proportional to $`\mathrm{sin}\mathrm{\Psi }`$ or $`\mathrm{cos}\mathrm{\Psi }`$) one can find that $`k_1=0`$. Then

$`S_{1,2}`$ $`=`$ $`(1/2)(P_{r1}^2+P_{z1}^2)/\mathrm{\Psi },`$
$`S_{2,2}`$ $`=`$ $`(1/2)(P_{r1}^2+P_{z1}^2)/r,`$
$`S_{3,2}`$ $`=`$ $`P_{r1}(_{}E_{r1}+E_{z1}/\mathrm{\Psi }),`$
$`S_{4,2}`$ $`=`$ $`P_{z1}(_{}E_{r1}+E_{z1}/\mathrm{\Psi }),`$
$`S_{5,2}`$ $`=`$ $`0,`$
and for magnetic field strength we obtain the following equation:

$$[(/r)_{}1]H_{\theta 2}=S_{4,2}/rS_{3,2}/\mathrm{\Psi }.$$
(14)
Solution of this equation is

$$H_{\theta 2}=f_1(r)+f_2(r)\mathrm{cos}(2\mathrm{\Psi }),$$
(15)
where $`f_{1,2}(r)`$ satisfies to equations

$`(\mathrm{\Delta }_{}r^21)f_1`$ $`=`$ $`(1/2)(d/dr)[R(\mathrm{\Delta }_{}1)R],`$
$`(\mathrm{\Delta }_{}r^21)f_2`$ $`=`$ $`(1/2)[(dR/dr)\mathrm{\Delta }_{}RR(d/dr)\mathrm{\Delta }_{}R],`$
here $`\mathrm{\Delta }_{}=_{}(d/dr)`$ is the transverse part of Laplacian. In the $`ith`$ approximation $`H_{\theta i}\mathrm{cos}(i\mathrm{\Psi })`$. The dependence $`k_2(r)`$ can be obtained in the third approximation. So, the weakly-nonlinear theory predicts that the magnetic field of 3D nonlinear plasma wave is not zero, in opposite to the 3D linear case and to 1D nonlinear one, and oscillates at the harmonics of the linear plasma frequency; the nonlinear wavelength changes in the radial direction.

We have solved Eqs. (1) - (6) numerically choosing the Gaussian profile of the driver both in longitudinal and transverse directions:

$$A(z,r)=A_0\mathrm{exp}[(zz_0)^2/\sigma _z^2]\mathrm{exp}(r^2/\sigma _r^2),$$
(16)
where $`A(z,r)`$ stands for $`\alpha `$ or $`a^2`$. Fig. 1 shows three-dimensional linear plasma wave excited by relativistic electron bunch. The linear numerical solution obtained well agreed with predictions of the linear theory.2 As is seen in Figure 1(b), in the linear case the magnetic field excited in plasma is localized in the range occupied by the bunch and in the wake $`E_z^{\mathrm{max}}H_\theta ^{\mathrm{max}}0`$, that correspond to the linear theory.

In the nonlinear regime the behavior of plasma waves is qualitatively changed. In Fig. 2 one can see a nonlinear plasma wave that is excited by an intense laser pulse. The main difference here from the linear case is the change of shape and length of the wave with the radial coordinate $`r`$ \[see Fig. 2(a)\]. The magnetic field strength in the nonlinear plasma wave as shown in Figure 2(b) is different from zero and has the magnitude comparable with that of other components of the field. The magnitude of higher frequency oscillations (as compared to the plasma frequency) performed by the magnetic field along $`z`$ grows in proportion to the nonlinearity. Such a behavior of magnetic field is a purely nonlinear effect. Indeed, in the linear case $`H_\theta =0`$, i.e., the contribution of momentum components at the plasma frequency in the expression (7) is compensated. The nonlinearity of the wave implies a rise of higher harmonics in $`P_z`$ and $`P_r`$. According to (7), the rise of magnetic field is due to these harmonics and this accounts for frequent oscillations seen in Figure 2(b), that conforms with the weakly-nonlinear analytical consideration presented above. On the other hand, the non-zero magnetic field means, according to (7), that the motion of plasma electrons in the nonlinear wave is turbulent ($`rot𝐏0`$). The degree of turbulence (the measure of which is $`H_\theta `$) grows in proportion to the nonlinearity.

It is easy to see that due to the dependence of the wavelength on$`r`$, the field in the radial direction grows more chaotic as the distance from the driver increases. In fact, the oscillations of plasma for different $`r`$are ”started” behind the driver with nearly equal phases but different wavelengths \[see Fig. 2(a)\]. As $`\left|z\right|`$ increases, the change of phase in transverse direction (for fixed $`z`$) becomes more and more marked. This leads to a curving of the phase front and to ”oscillations” in the transverse direction.6,10 The radial behavior of the longitudinal electric field and magnetic field of the nonlinear plasma wave is presented in Fig. 3. Qualitatively the radial dependence of the field differs from that of the linear case by the change of sign and ”steepening” of fields along $`r`$ and that is connected with the curvature of wave phase front.

So, we have found that the magnetic field of nonlinear plasma wave is different from zero and performs higher frequency oscillations as compared to the plasma frequency. The latter qualitatively differs from the linear case and one-dimensional nonlinear one. The obtained numerical results well agreed with the predictions of weakly-nonlinear theory.

REFERENCES

1T. Tajima and J. M. Dawson, Phys. Rev. Lett. 43, 267 (1979).

2R. D. Ruth, A. W. Chao, P. L. Morton, and P. B. Wilson, Part. Accel. 17, 171 (1985); P. Chen, Part. Accel. 20, 171 (1987).

3A. I. Akhiezer and R .V. Polovin, Zh. Eksp. Teor. Fiz. 30, 915 (1956) \[Sov. Phys. JETP 3, 696 (1956)\].

4N. A. Krall and A. W. Trivelpice, Principles of Plasma Physics (McGraw-Hill Book Co., NY, 1973).

5B. N. Breizman, T. Tajima, D. L. Fisher, and P. Z. Chebotaev, In: Research Trends in Physics: Coherent Radiation and Particle Acceleration , edited by A. Prokhorov (American Institute of Physics, New York, 1992), pp. 263-287; K. V. Lotov, Phys. Plasmas 5, 785 (1998).

6A. G. Khachatryan and S. S. Elbakian. Two-dimensional nonlinear regime in the Plasma Wakefield Accelerator, Proceedings PAC’99, edited by A.Luccio and W.MacKay (IEEE, Piscataway, NJ, 1999), pp. 3663-3665.

7G. B. Witham, Linear and Nonlinear Waves (Wiley, New York, 1974), Chap. 13.

8A. G. Khachatryan, Phys. Rev. E 58, 7799 (1998).

9R. Keinigs and M. E. Jones, Phys. Fluids 30, 252 (1987); A. G. Khachatryan, A. Ts. Amatuni, S. S. Elbakian, and E. V. Sekhpossian, Plasma Phys. Rep. 22, 576 (1996).

10A. G. Khachatryan, Phys. Rev. E 60, 6210 (1999).

FIGURE CAPTIONS

Figure 1. The linear plasma wave excited by an electron bunch with gaussian profile;$`\alpha _0=0.1`$, $`\sigma _z=2`$,$`\sigma _r=0.5`$, $`\gamma =10`$. (a). The dimensionless strength of longitudinal electric field. 1—$`E_z(z)`$ at the axis, $`r=0`$; 2 — $`r=1`$; 3 — $`r=2`$. (b). The strength of azimuthal magnetic field. 1 — $`r=0.2`$; 2 — $`r=0.75`$; 3 — $`r=2`$.

Figure 2. The nonlinear plasma wave excited by laser pulse. The pulse parameters are : $`a_0^2=3.6`$, $`\sigma _z=2`$,$`\sigma _r=5`$, $`\gamma =10`$. (a). The longitudinal electric field $`E_z`$. $`r=0,2,4`$ and $`5`$ in the order of magnitude reduction. (b). The magnetic field strength. 1 — $`r=2`$; 2 — $`r=4`$.

Figure 3. The radial behavior of field for the case given in Figure 2 for $`z=25.`$ 1 — longitudinal electric field $`E_z(z=25,r)`$; 2 — magnetic field $`H_\theta (z=25,r)`$.


text:
# Anisotropic Disorder in High Mobility 2D Heterostructures and its Correlation to Electron Transport

## Abstract

Surface morphology of high mobility heterostructures is examined and correlated with d.c. transport. All samples examined show evidence of lines in the \[1$`\overline{1}`$0\] direction with roughness ranging from small amplitude features to severe anisotropic ridges. Transport in these samples is consistent with that in samples having artificially induced 1D charge modulations. The native surface properties reflect a prevalent, anisotropic disorder affecting 2D electron conduction. Importantly, the native lines are orthogonal to the stripes theoretically proposed to explain high Landau level transport anisotropies.

Correlated electron effects in 2D heterstructures can be observed because the disorder of these structures has been understood and subsequently controlled to facilitate low energy electronic interactions. The fractional quantum Hall effect(FQHE) was observed in modulation doped heterostructures as this specific device provides spatial separation of ionized impurities and the interacting 2D electron system (2DES). Molecular beam epitaxy (MBE) has produced extraordinarily pure systems with electron mobilities exceeding 107 cm2/V-sec. Such structures exhibit phenomena ranging from composite fermion Fermi surface effects at $`\nu `$=1/2, to transport anisotropies at high Landau levels, proposed to be electron stripe and bubble phases. In the composite fermion system, density variations play a particularly significant role as disorder due to the association of charge and fictitious magnetic field. Establishing sources of lateral charge variation within the high mobility interface is therefore important for understanding disorder influencing the correlated 2DES.

The surface of high mobility heterostructures provides an accessible plane to study possible density related disorder. Surface morphology can vary under different MBE growth conditions, and growth anisotropies can occur along certain crystallographic directions. The length scales of these variations are large, on the order of microns, and so are imminently measureable. Key questions are then, do these large-scale imperfections occur in high mobility heterostructures, and do they influence the electronic properties of the 2D layer?

To answer this we have examined surface morphology using atomic force microscopy and studied transport of a series of high mobility heterostructures. In an extensive pool of high mobility wafers we observed a range of surface properties related to known crystal growth instabilities. Our principal finding is that all surfaces are observed to have some formation of ridges along the \[1$`\overline{1}`$0\] direction. These lines exist within a range of large-scale topographic fluctuations that vary from nearly isotropic mounds to stark anisotropic ridges consistently along the \[1$`\overline{1}`$0\] direction. Transport in samples with severe surface roughness is shown to be similar to that observed in heterostructures with artificially induced 1D charge density modulations. In samples with the less severe native surface roughness, transport anisotropies can be identified in the lowest Landau level which are consistent with the surface morphology.

A focus of this study is the surface morphology of ultra-high mobility samples that demonstrate both extensive FQHE features and high Landau level transport anisotropies. These high quality samples show i) FQHE states manifest properties associated with the anisotropic growth, and ii) the native surface lines, which we associate with charge modulation as shown by transport, are oriented orthogonally to the stripes proposed in the theoretical models of high Landau level anisotropies.

The MBE-produced wafers in this study were all grown on  oriented substrates with a buffer layer of minimally 1 $`\mu `$m thick GaAs or AlGaAs/GaAs superlattice. Transport properties and surface morphology were independent of whether heterostructures are single interfaces or quantum wells. The 2D layer resides less than 300 nm below the surface for all samples studied. The wafers focused upon here (and shown in Figures 1a-c) were produced under essentially the same conditions, indicating that the properties determining surface morphology are not well understood.

A large number of heterstructure wafers with mobilities greater than 8$`\times `$106cm2/V-sec were examined, and shown here are samples which are representative of the full range of surface morphologies. Light microscopy was also used to examine heterostructure surfaces; using phase contrast optics, most of the larger surface features can be detected readily at relatively low magnifications.

Surface morphology was examined using atomic force microscopy operating in tapping mode. Topographic images (Figure 1, right column) depict surface height variations. The change in cantilever oscillation amplitude, proportional to the derivative of topographic changes, is also shown (Fig. 1, left column). Transport was performed using standard lock-in techniques at low frequencies.

Figure 1 shows surface morphology for heterostructures ranging from smoothest (a) to roughest (d). All samples show features, lines, along the \[1$`\overline{1}`$0\] direction. In the lowest overall amplitude surface topography, Figure 1(a), the dominant morphology is a coarsely isotropic “orange peel” effect with a diameter of about 4 $`\mu `$m and height range of less than 40 Å. Within these large features, fine lines corresponding to atomic scale terracing are observed alligned in the \[1$`\overline{1}`$0\] direction. This predilection toward line formation in the \[1$`\overline{1}`$0\] direction becomes more apparent as the overall surface roughness increases. In Figure 1(b) the orange peel texture has evolved to a surface with mounds clearly oriented in the \[1$`\overline{1}`$0\] direction, and with vertical excursion of almost 100 Å. Figure 1(c) shows another representative sample in the range of common surface morphologies. Severe ridges are present with about the same period as that shown in Figure 1(b) and with almost the same vertical excursion. It is important to note that wafers generally display morphology ranging from that shown in (a) to that shown in (c). The final AFM image (Figure 1d) shows an extreme case of surface roughness: this was induced by growing a particulary thick buffer layer of about 10 $`\mu `$m, capped by a standard interface. The vertical excursion is large, up to 300 Å, with line formation distinctly in the \[1$`\overline{1}`$0\] direction. Note that ridges can traverse the full field of the scan.

The common property of the above described range of morphologies is the formation of lines in the \[1$`\overline{1}`$0\] direction. This effect is due to a previously recognized MBE growth instability  with fast growth rate along the direction of the As dangling bonds. Extensive growth or lowered temperatures during growth can accentuate this effect.

A second mechanism contributing to this formation of lines in the \[1$`\overline{1}`$0\] direction is from the AlGaAs layering. A temperature and Al concentration dependent growth preference in the \[1$`\overline{1}`$0\] direction has been documented. It provides a smaller length scale ridge formation than the above described mechanism, and is at work in these samples due to the presence of GaAs/AlGaAs layering. Suppressing this contribution to ridge formation involves growth at the appropriate temperature and Al concentration; suppression of the previously noted growth instability is more difficult.

Transport results show that the surface morphology is associated with charge variation at the 2D electron layer. This precept is supported by the data of Figure 2: it shows typical transport data from the sample of Figure 1d having severe natural lines, with current driven across and along the \[1$`\overline{1}`$0\] direction. The key features in this data are the large minimum in resistance at 1/2 filling factor for current driven across \[1$`\overline{1}`$0\] (across the native lines), accompanied by a peak at 1/2 for current driven along \[1$`\overline{1}`$0\]. Note also that current driven across the native lines produces narrow FQHE states; conversely, current driven along the lines produces relatively widened FQHE states.

These effects are similar to those observed in transport in the lowest Landau level for a 1D charge density modulation artificially induced on a heterostructure: results from such a sample are shown in Figure 3. Artificial 1D charge density modulation is achieved by etching an array of lines into the surface of the heterostructure. The etch depth can be small, typically about 150 Å, to induce transport consequences. The artificially modulated system shows a peak at 1/2 for current driven along the imposed lines and a large minimum for current driven across the lines, just as seen in the sample of Figure 1(d) with native lines. In addition, the FQHE states show relative narrowing when current is driven across the modulation lines and broadening when driven along the modulation lines; again, this is similar to the effect observed for the respective transport directions in the native line sample of Figure 1(d). The central minimum and peak at 1/2 for the two transport directions in the artificially modulated sample of Figure 3 show small features at and symmetrically about 1/2. The lateral minima (peaks) correspond to geometric resonance of the charged carrier motion with the regular period of the modulation, and the central minimum (peak) at 1/2 is due to “snake” states. These occur as a consequence of large radius composite fermion cyclotron orbit motion along an effective B-field gradient centered about zero effective magnetic field, as expected to exist along the lines of a charge density gradient. Note the presence of such snake state transport at 1/2 in the sample with strong native lines; the geometric resonances are absent due to the aperiodic nature of the native charge modulation. Other features are common to the magnetotransport spectra of Figures 2 and 3; for example note the large peak between $`\nu =`$ 3/2 and 1 in the red traces. The presence of these multiple specific features in both the native line sample and the artificially modulated sample clearly show empirically that the native surface lines in Figure 1d are indeed associated with charge modulations in the 2D electron layer.

Less coarse surface effects also produce transport consequences in the lowest Landau level. Shown in Figure 4(a) is the magnetotransport around 3/2 for the sample pictured in Figure 1(c). The native modulations of this sample induce transport properties similar to those noted in the artificially charge modulated system: the FQHE states are wider in the orientation of current along the \[1$`\overline{1}`$0\] direction (along the lines). This broadening of about $`\delta B/B`$ 0.5% reflects a spread in the density of $`\delta n/n`$ 0.5%. This property was found in all of the large number of samples of similar surface morphologies that were examined. Additionally, current driven along the lines results in a less well developed minimum at 3/2 and poorly formed high order fractions when compared to transport in the orthogonal direction. These specific effects are also observed in samples with less severe surface morphology. The sample shown in Figure 1b demonstrates similar broadening of the FQHE for current along the \[1$`\overline{1}`$0\] direction, particularly at low temperatures ($`<`$ 100 mK). Samples with morphology approaching that of Figure 1a can likewise demonstrate these effects at very low temperatures.

A striking finding of our study is that samples with relatively severe surface mophology demonstrate beautifully the transport anisotropies in the high Landau levels that have recently been documented and explained using a new striped phase ground state. Figure 4(b) shows transport in two directions from the sample of Figure 1(c), which has substantial native line formation in the \[1$`\overline{1}`$0\] direction. Figure 4(b) shows the characteristic deep minimum at 9/2 for current driven across the \[1$`\overline{1}`$0\] direction (across the native lines) and also shows a peak at 9/2 for current along the native lines. Empirically this is the same as the transport properties observed in the lowest Landau level for the artificially modulated 2DES. Most importantly, the native lines observed in our morphology studies, which apparently act as charge lines as supported by the sample’s lowest Landau level transport, are orthogonal to the theoretically proposed stripes used to explain the high Landau level transport. The sample from Figure 1b, which has less severe but still significant anisotropy in the surface morphology shows similar high Landau level anisotropies.

In summary, surface morphology of high mobility heterostructures has been examined and found to show line formation in the \[1$`\overline{1}`$0\] direction with varying degrees of severity over the full complement of wafers tested. These surface properties are correlated to transport results, with even mild surface roughness associated with distinct transport features. These transport features are like those seen specifically in samples with artificial 1D charge density modulations, suggesting similar transport mechanisms are at work in both systems. These pervasive native lines represent an important form of systematic, anisotropic disorder that is of particular significance in ultra-high mobility samples, as it is here that their consequences can be readily observed.

This finding of native density anisotropy in ultra-high mobility samples has particular pertinence in considering transport in high Landau levels, since the native lines are orthogonal to the direction of the stripes proposed theoretically to explain the transport anisotropy. It is tempting to consider that the same mechanism responsible for effects observed around 1/2 in the artificial 1D charge modulated samples is at play at 9/2 in the ultra-high mobility samples, with the 1D charge modulation provided by the native line formation. In spite of numerous empirical similarities, however, this is unlikely as Fermi surface formation would be required at the high Landau level half-filling. The native lines may have some influence in breaking the symmetry of the system, and so may induce orientation of the proposed ground state. This remains to be tested.

We show that surface morphology and its commensurate charge distribution have an impact on correlations in the 2DES. Further advancement in heterostructures designed to support electron correlation effects must overcome this systematic disorder.