text:
# The magnetisation of protoplanetary disks

## 1 Introduction

Meteors are commonly considered to be relics of the protoplanetary nebula in our early solar system. Of particular importance in the context of the early evolution of the solar system is the remanent magnetisation of meteorite material that indicates the presence of relatively strong magnetic fields ($`B=110`$ G) in protoplanetary nebula (Levy & Sonett, 1978; Sugiura et al., 1979). The formation of meteorite chondrules may also demand for magnetic fields as the relevant energy source (Sonett 1979; Levy & Araki 1989). Thus, one faces the question of the origin of such strong magnetic fields in protoplanetary disks. Recent numerical calculations of collapsing magnetised molecular cloud cores have proven that ambipolar diffusion can solve the magnetic flux loss problem of the formation of protostellar objects (, 1993; Desch & Mouschovias, 2001). The magnetic flux can diffuse to the region around the protostellar object and thus result in a magnetic field of the protoplanetary disk of some 0.1 G (Desch & Mouschovias, 2001). The fate of this field after the stellar object has fully evolved and accretion flow in the disk sets in is, however, not clear. In this contribution we consider disk magnetisation in a later stage of the stellar evolution, namely after a stellar object with a dipole magnetic field has formed.

A characteristic feature of partially ionised plasmas is the possibility that magnetic fields can be generated without any seed fields. As shown by (Huba & Fedder 1993) a relative plasma-neutral gas flow with non-vanishing vorticity results inevitably in the generation of magnetic fields. This effect is described by the generalised Ohm’s law on the grounds of a multi-fluid description. The magnetisation mechanism works significantly more efficiently in a dusty plasma (Birk et al. 1996; Birk et al. 2001) where the inertia effects and collisional momentum transfer are dominated by the heavy dust grains.

In this contribution we show by means of multi-fluid simulations that relative shear flows lead to a very fast three-dimensional magnetic field self-generation in dusty protoplanetary disks. Magnetic field strengths up to one Gauss are generated on the time scale of 1yr on spatial scales of some astronomical units (AU).

## 2 Self-magnetisation of protoplanetary accretion disks

We consider the following balance equations that govern the low-frequency macroscopic dynamics of a partially ionised dusty plasma (see Birk et al. 1996; Birk et al. 2001)

$$\frac{\rho _d}{t}=(\rho _d𝐯_𝐝)$$
(1)
$$\frac{\rho _i}{t}=\left(\rho _i𝐯_𝐢\right)$$
(2)
$$\frac{\rho _n}{t}=(\rho _n𝐯_𝐧)$$
(3)
$`{\displaystyle \frac{(\rho _d𝐯_𝐝)}{t}}=`$ $``$ $`(\rho _d𝐯_𝐝𝐯_𝐝)(p_e+p_i+p_d)+(\rho _i+\rho _d)𝐠`$ (4)
$`+`$ $`(\times 𝐁)\times 𝐁\nu _{dn}\rho _d(𝐯_𝐝𝐯_𝐧)`$
$``$ $`\nu _{in}\rho _i\left(𝐯_𝐢𝐯_𝐧\right)\nu _{en}\rho _e\left(𝐯_𝐞𝐯_𝐧\right)`$
$`{\displaystyle \frac{(\rho _n𝐯_𝐧)}{t}}=`$ $``$ $`(\rho _n𝐯_𝐧𝐯_𝐧)p_n+\rho _n𝐠`$ (5)
$`+`$ $`\nu _{dn}\rho _d(𝐯_𝐝𝐯_𝐧)+\nu _{in}\rho _i\left(𝐯_𝐢𝐯_𝐧\right)`$
$`+`$ $`\nu _{en}\rho _e\left(𝐯_𝐞𝐯_𝐧\right)`$
$`{\displaystyle \frac{1}{\gamma _e1}}{\displaystyle \frac{p_e}{t}}=`$ $``$ $`{\displaystyle \frac{1}{\gamma _e1}}(p_e𝐯_𝐞)p_e𝐯_𝐞+\rho _n\nu _{ne}(𝐯_𝐞𝐯_𝐧)^2`$ (6)
$``$ $`2{\displaystyle \frac{\rho _d\nu _{de}}{m_d}}\left({\displaystyle \frac{k_BT_e}{\gamma _e1}}{\displaystyle \frac{k_BT_d}{\gamma _d1}}\right)2{\displaystyle \frac{\rho _i\nu _{ie}}{m_i}}\left({\displaystyle \frac{k_BT_e}{\gamma _e1}}{\displaystyle \frac{k_BT_i}{\gamma _i1}}\right)`$
$``$ $`2{\displaystyle \frac{\rho _n\nu _{ne}}{m_n}}\left({\displaystyle \frac{k_BT_e}{\gamma _e1}}{\displaystyle \frac{k_BT_n}{\gamma _n1}}\right)+\rho _d\nu _{de}\left(𝐯_𝐞𝐯_𝐝\right)^2`$
$`+`$ $`\rho _i\nu _{ie}(𝐯_𝐞𝐯_𝐢)^2`$
$`{\displaystyle \frac{1}{\gamma _i1}}{\displaystyle \frac{p_i}{t}}=`$ $``$ $`{\displaystyle \frac{1}{\gamma _i1}}(p_i𝐯_𝐢)p_i𝐯_𝐢+{\displaystyle \frac{m_n}{m_i+m_n}}\rho _i\nu _{in}(𝐯_𝐢𝐯_𝐧)^2`$ (7)
$`+`$ $`{\displaystyle \frac{m_d}{m_d+m_i}}\rho _i\nu _{id}(𝐯_𝐢𝐯_𝐝)^22{\displaystyle \frac{\rho _i\nu _{in}}{m_i+m_n}}\left({\displaystyle \frac{k_BT_i}{\gamma _i1}}{\displaystyle \frac{k_BT_n}{\gamma _n1}}\right)`$
$``$ $`2{\displaystyle \frac{\rho _i\nu _{id}}{m_i+m_d}}\left({\displaystyle \frac{k_BT_i}{\gamma _i1}}{\displaystyle \frac{k_BT_d}{\gamma _d1}}\right)2{\displaystyle \frac{\rho _i\nu _{ie}}{m_i}}\left({\displaystyle \frac{k_BT_i}{\gamma _i1}}{\displaystyle \frac{k_BT_e}{\gamma _e1}}\right)`$
$`+`$ $`{\displaystyle \frac{m_e}{m_i}}\rho _i\nu _{ie}(𝐯_𝐢𝐯_𝐞)^2`$
$`{\displaystyle \frac{1}{\gamma _d1}}{\displaystyle \frac{p_d}{t}}=`$ $``$ $`{\displaystyle \frac{1}{\gamma _d1}}(p_d𝐯_𝐝)p_d𝐯_𝐝+{\displaystyle \frac{m_e}{m_d}}\rho _d\nu _{de}(𝐯_𝐝𝐯_𝐞)^2`$ (8)
$`+`$ $`{\displaystyle \frac{m_i}{m_i+m_d}}\rho _d\nu _{di}(𝐯_𝐝𝐯_𝐢)^2+{\displaystyle \frac{m_n}{m_n+m_d}}\rho _d\nu _{dn}(𝐯_𝐝𝐯_𝐧)^2`$
$``$ $`2{\displaystyle \frac{\rho _d\nu _{di}}{m_i+m_d}}\left({\displaystyle \frac{k_BT_d}{\gamma _d1}}{\displaystyle \frac{k_BT_i}{\gamma _i1}}\right)2{\displaystyle \frac{\rho _d\nu _{dn}}{m_n+m_d}}\left({\displaystyle \frac{k_BT_d}{\gamma _d1}}{\displaystyle \frac{k_BT_n}{\gamma _n1}}\right)`$
$``$ $`2{\displaystyle \frac{\rho _d\nu _{de}}{m_d}}\left({\displaystyle \frac{k_BT_d}{\gamma _d1}}{\displaystyle \frac{k_BT_e}{\gamma _e1}}\right)`$
$`{\displaystyle \frac{1}{\gamma _n1}}{\displaystyle \frac{p_n}{t}}=`$ $``$ $`{\displaystyle \frac{1}{\gamma _n1}}(p_n𝐯_𝐧)p_n𝐯_𝐧`$ (9)
$`+`$ $`{\displaystyle \frac{m_i}{m_i+m_n}}\rho _n\nu _{ni}(𝐯_𝐢\mathrm{-}𝐯_𝐧)^2+{\displaystyle \frac{m_e}{m_n}}\rho _n\nu _{ne}\left(𝐯_𝐞\mathrm{-}𝐯_𝐧\right)^2`$
$`+`$ $`{\displaystyle \frac{m_d}{m_n+m_d}}\rho _n\nu _{nd}(𝐯_𝐧𝐯_𝐝)^22{\displaystyle \frac{\rho _n\nu _{ni}}{m_i+m_n}}\left({\displaystyle \frac{k_BT_n}{\gamma _n1}}{\displaystyle \frac{k_BT_i}{\gamma _i1}}\right)`$
$``$ $`2{\displaystyle \frac{\rho _n\nu _{nd}}{m_n+m_d}}\left({\displaystyle \frac{k_BT_n}{\gamma _n1}}{\displaystyle \frac{k_BT_d}{\gamma _d1}}\right)`$
$``$ $`2{\displaystyle \frac{\rho _n\nu _{ne}}{m_n}}\left({\displaystyle \frac{k_BT_n}{\gamma _n1}}{\displaystyle \frac{k_BT_e}{\gamma _e1}}\right)`$
$`{\displaystyle \frac{𝐁}{t}}=`$ $``$ $`{\displaystyle \frac{m_ic}{e}}\times \left({\displaystyle \frac{p_i}{\rho _i}}\right)+{\displaystyle \frac{m_i}{m_d}}z_d\times \left({\displaystyle \frac{\rho _d}{\rho _i}}𝐯_𝐝\times 𝐁\right)`$ (10)
$`+`$ $`{\displaystyle \frac{m_ic}{4\pi e}}\times \left(\left({\displaystyle \frac{\times 𝐁}{\rho _i}}\right)\times 𝐁\right)\eta \mathrm{\Delta }𝐁`$
$``$ $`{\displaystyle \frac{m_ic}{e}}\times \left[{\displaystyle \frac{n_d}{n_i}}\left(\left(z_d{\displaystyle \frac{n_i}{n_d}}\right)\nu _{id}+z_d\nu _{in}\right)𝐯_d\nu _{in}𝐯_n\right]`$
where $`\rho _\alpha `$, $`𝐯_\alpha `$, $`m_\alpha `$, $`p_\alpha `$ and $`𝐁`$ are the respective ($`\alpha =d,i,e,n`$) mass density, bulk velocity, particle mass (the ions are assumed to be protons), thermal pressure ($`p_\alpha =\rho _\alpha k_BT_\alpha /m_\alpha `$, where $`k_B`$ is the Boltzmann constant and $`T_\alpha `$ is the temperature) and magnetic field strength. The adiabatic ratios of specific heats are chosen as $`\gamma _\alpha =5/3`$. The collision frequencies for collisions (elastic collisions as well as charge exchange) between species $`\alpha `$ and $`\beta `$ are denoted by $`\nu _{\alpha \beta }`$ and satisfy $`\rho _\alpha \nu _{\alpha \beta }=\rho _\beta \nu _{\beta \alpha }`$. By $`\eta `$ the collisional diffusivity ($`\eta =m_ic^2(\nu _{id}+\nu _{in})/4\pi n_ie^2`$) is denoted. Ionisation and recombination are neglected in the present context. The electron density is calculated from the quasi-neutrality condition ($`\rho _i/m_i\rho _e/m_ez_d\rho _d/m_d=0`$) where $`z_d`$ is the negative dust charge number. The momentum transfer equations that include external gravitational forces with the acceleration $`𝐠`$ are formulated as the sum of the balance equations of all charged components under the assumption of inertialess ions and electrons (as compared to the massive dust grains). Most important for the self-generation of magnetic fields is the induction equation (see also Ciolek & Mouschovias, 1993) which may be derived from the inertialess ion momentum equation. The final term in Eq. (10) is the relevant magnetic field self-generation term. The mass, momentum and pressure balance equations as well as the induction equation are integrated numerically by the DENISIS code (Schröer et al. 1998) under the assumption of ionisation equilibrium. This explicit finite differences code is based on a modified Leapfrog algorithm. Extensive tests and applications of the code can be found in the papers by Schröer et al. (1998), Birk et al. (2001) and Birk and Wiechen (2002).

Although the present study is motivated by the local information of the magnetised meteorite material in the solar system we assume that our solar system is a standard case for normal stellar systems. Thus, in the following we consider a simple protoplanetary accretion disk model consisting of electrons, protons, neutrals and charged dust grains using typical dynamical parameters for circumstellar disks around T Tauri stars (Montmerle 1991; Königl 1994). In particular, we choose the following parameters: $`n_i=1.01n_d=0.1n_n=10^{12}\mathrm{cm}^3`$, $`z_d=1`$, $`m_d=100m_i`$, and for the temperatures $`T_n=T_i=T_d=500`$K. The ion- neutral collision frequencies can be estimated as (e.g. Krall & Trivelpiece 1986) $`\nu _{in}\sigma (T_i)n_n\sqrt{kT_i/m_i}`$ where $`k`$ is the Boltzmann constant and $`\sigma `$ denotes the cross section. The latter depends on the ion temperature in a way that results in a collision frequency that is independent on the ion temperature (for more details on the collision theory see, e.g., Mitchner & Kruger 1973). We use $`\sigma =510^{15}\mathrm{cm}^2`$ as constant approximation of sufficient accuracy in our context. The ion-dust collision frequency is $`\nu _{id}4\sqrt{2\pi }n_dz_d^2e^4\mathrm{ln}(\lambda _D/a)/3m_i^{1/2}(kT_i)^{3/2}`$ (Benkadda et al. 1996) where $`\lambda _D=\sqrt{kT_d/4\pi n_dz_d^2e^2}`$ denotes the dust Debye screening length. For the dust grain size we choose $`a=10^6`$cm. The actual choice of the grain size has no significant influence on the induction equation, since $`a`$ only enters in the logarithm of the formula for the ion-dust frequency and does not play an important role for the magnetic field generation due to the strong electron depletion ($`z_dn_i/n_d0`$ in Eq. 10). On the other hand, the momentum transfer between the dust and the neutral fluid components depends on the grain radius. However, the continuous acceleration (accretion flow) guarantees a finite relative dust-neutral shear flow. Without this energy source larger dust grains would result in faster saturation and weaker generated magnetic fields.

Within a region of 10 AU we consider an initially homogenous protoplanetary disk in which the dust grains as well as the ions are chosen to be in corotation within a region of about 1 AU (Fig.1, upper plot). The corotation of the charged dust grains is caused by a stellar dipole field. The neutrals do not couple directly to the magnetic field but via momentum transfer with the charged particles. Outside this corotation radius the dust is in Keplerian motion (Li 1996). The neutral gas is assumed to be in Keplerian motion with the exception of a very inner region of 0.2 AU where the neutral gas is assumed to be in corotation to avoid singularities (Fig.1, lower plot). The typical velocities near the transition regions between corotation and Keplerian motion are about $`110^6`$ cm s-1 for the dust and $`1.510^7`$ cm s-1 for the neutral gas, respectively.

The orders of magnitude of these velocities are in a rough agreement with a central mass object of about 2 solar masses, so that we include an external gravitational force corresponding to such a central mass. The initial configuration is characterised by the immediate onset of idealised accretion flow rather than to be in a stationary state. The gravitational field is the ultimate energy source for the generation of the magnetic field.

The simulation box is given by $`x[25,25]`$, $`y[25,25]`$ and $`z[0,10]`$ in units of 0.2 AU. The equidistant grid is resolved by $`77\times 77\times 31`$ grid points. The boundary conditions are symmetric at all boundaries, i.e. the values of the scalar quantities and vector components at the numerical boundaries correspond to the values at the respective last inner grid points.

As the upper plot of Fig. 2 indicates, a swirl-like dust flow in the central disk in the $`x`$-$`y`$-plane has built up after $`0.7`$yr. The curl of the relative dust-neutral gas flow (Fig. 2, lower plot) is sustained by the gravitational pull and results in an ongoing magnetisation of the disk.

The magnetic field after $`t=0.7`$yr and $`t=1.4`$yr is shown in Figs. 3 and 4, respectively. The strengths of the different field components are roughly comparable. After $`t=0.7`$yr, the amplitude of $`B_z=0.8`$G is 5 times the amplitudes of $`B_x`$ and $`B_y`$ and after $`t=1.4`$yr, the factor is about 2 ($`B_z=1.2`$G). The field generation is caused not only by the poloidal flows that develop from the initial configuration but also by vertical relative shear flows arising in the dusty accretion disk.

The magnetic energy density $`W_B=B^2/8\pi `$ at different heights of the disk are shown in Fig. 5. Obviously, the inner part of the disk is completely magnetised in all three dimensions.

Since the ultimate source of the self-magnetisation is the stellar gravitation the process lasts as long as the accretion flow in the circumstellar disk is present. The numerical simulation, however, fails when the central density gradient becomes too steep (after about 200 dynamical times). This is a consequence of our simulation approach that does not properly deal with the accretion of matter (loss of plasma to the stellar object) and the formation of bipolar jets. A global multi-fluid model including all relevant physics is clearly beyond the scope of the present state of the art.

## 3 Discussion

We studied the magnetisation of protoplanetary disks that consist of dusty plasmas. Different from previous studies on the magnetisation of the environment of collapsing molecular cloud cores (, 1993; Desch & Mouschovias, 2001) we dwell on the magnetic field generation in fully developed accretion disks. The process is driven by accretion, i.e. the ultimate source of energy for the magnetic field is the stellar gravitational field, it results in an ordered poloidal magnetic field of some Gauss. It is a very fast and efficient mechanism of field generation. Also it is a process of permanent magnetisation as long as some accretion flow is present.

The numerical simulations show that magnetic fields are generated in circumstellar dusty disks by macroscopic relative shear motions between the different species. This self-magnetisation mechanism that works in addition to the formation of turbulent magnetic fields due to the interaction of the stellar field with a turbulent accretion disk (Bardou & Heyvaerts 1996).

We note that in our studies we did not include a stellar magnetic field explicitly. Such a field should give rise to a (sheared) wind component close to the star and thus, possibly to the further magnetisation of a significant fraction of the magnetosphere. In this contribution we concentrated on the generation of magnetic fields in the accretion disks. Further studies that include stellar fields are a promising task for the future.

In our simulations we did not include a gap between the central stellar object and the disk for technical reasons. However, since we were interested in the in situ magnetisation of the disk material this idealisation does not influence our findings.

A further idealisation of our study is the initial density homogeneity of the disk. However, in the course of the dynamical evolution radial profiles form in all fluid components due to the external gravitation. More realistic density profiles perpendicular to the disk would demand for the inclusion of self-gravitation. The inclusion of this further energy source that also may contribute to the self- magnetisation of the disk is beyond the scope of the present investigations but could be of interest for future numerical modelling. Parameter studies, however, have shown that the efficiency of the process of magnetisation is not significantly influenced by the variation of the initial number densities over a wide parameter range.

Given the uncertainties of the actual plasma parameters the high efficiency of the discussed magnetisation guarantees that protoplanetary disks are expected to generate their own magnetic fields of the order of 0.1 – 1 Gauss. In the presence of a small radial magnetic field component the Balbus-Hawley shear instability (Balbus & Hawley 1991) can, in principle, result in a faster than linear growth of the poloidal magnetic field component provided the appropriate onset criterion is fulfilled. This criterion, however, is not known for multifluid systems. On the other hand, localized reconnection may result in a redistribution of magnetic flux. Investigations on this process are beyond the scope of the present contribution and will be considered in further more detailed work.

If circumstellar disks carry their own significant magnetic fields, dynamic phenomena in the coupled star-disk system are rather complex. The disk field should influence characteristic phenomena as, e.g., accretion, winds and jets and magnetic activity, significantly.


format: text

text:
# ITEP-LAT/2002–29 KANAZAWA 02–41 Finite temperature phase transition in lattice QCD with 𝑁_𝑓=2 nonperturbatively improved Wilson fermions at 𝑁_𝑡=8

## 1 Introduction

Determination of the critical temperature of the chiral phase transition is one of the important nonperturbative problems in QCD to be addressed in lattice simulations. Bielefeld group and CP-PACS collaborations using two different types of improved lattice action for fermions were able to estimate $`T_c`$ in the chiral limit and their values are in good agreement . Still there are many sources of systematic uncertainties and new computations of $`T_c`$ with different actions are useful as an additional check. We made first large scale simulations of the nonperturbatively $`O(a)`$ improved Wilson fermion action. Moreover we performed simulations with the lattice spacing $`a`$ substantially smaller than in studies by Bielefeld and CP-PACS groups.

The following fermionic action is employed in our study:

$$S_F=S_F^{(0)}\frac{\mathrm{i}}{2}\kappa gc_{sw}a^5\underset{x}{}\overline{\psi }(x)\sigma _{\mu \nu }F_{\mu \nu }\psi (x),$$
(1)
where $`S_F^{(0)}`$ is the original Wilson action, $`c_{sw}`$ is the clover coefficient determined nonperturbatively . We use Wilson gauge field action.

The action (1) has been used in $`T=0`$ studies of lattice QCD by UKQCD and QCDSF collaborations . These studies confirmed that $`O(a)`$ lattice artifacts are suppressed as expected. To fix the physical scale and $`m_\pi /m_\rho `$ ratio we used results obtained by these collaborations . So far only $`N_t=4`$ and $`6`$ finite temperature results obtained with action (1) are available . These results were obtained at a rather large quark mass ($`m_\pi /m_\rho >0.85`$). In this work me make simulations with $`N_t=8`$ what allows us to decrease $`m_\pi /m_\rho `$ down to $`0.8`$. We choose the spatial extension of the lattice $`N_s=16`$ as a compromise between computational burden and need to reduce the finite size effects.

It is known that the order parameter of the finite temperature phase transition in quenched QCD is the Polyakov loop and the corresponding symmetry is global $`Z(3)`$ symmetry. In the chiral QCD the order parameter of the chiral symmetry breaking transition is the chiral condensate $`<\overline{\psi }\psi >`$. As numerical results show  both order parameters can be used to locate the transition point at intermediate values of the quark mass. We use Polyakov loop susceptibility to determine the transition temperature.

## 2 Simulation details

Hybrid Monte Carlo algorithm with parameters $`\delta \tau =0.0125`$, $`n_\tau =20`$, providing acceptance rates of about 70% is used in our simulations. The simulations were done on SR8000 Hitachi at KEK, Tsukuba and MVS 1000M at Joint Supercomputer Center, Moscow. For analysis SX5 NEC at RCNP and PC-cluster at ITP, Kanazawa were employed. Our code performs at the speed of 2.4 GFlops per node on Hitachi computer. We needed from 1000 ($`\tau N_{traj}n_\tau \delta \tau =250`$) to 3000 ($`\tau =750`$) trajectories for thermalization, depending on $`\kappa `$ and $`\beta `$. For runs started from configurations generated at the adjacent $`\kappa `$ this value was much lower. We determined the transition temperature at two values of $`\beta `$, 5.2 and 5.25, varying $`\kappa `$ as shown in Table 1.

We fixed the maximally abelian (MA) gauge on generated configurations. The simulated annealing algorithm has been used to fix the gauge. The advantages of this algorithm in comparison with the usual iterative algorithm has been demonstrated in $`SU(2)`$ pure gauge theory for MA gauge  and Maximal Center gauge .

Abelian and monopole Polyakov loops and their susceptibilities were measured on gauge fixed configurations. We found that for abelian and monopole observables the signal/noise ratio was better than that for gauge invariant nonabelian observables. This observation is in agreement with the results from quenched QCD at $`T>0`$ and both quenched and unquenched QCD at $`T=0`$. In particular the maximum of the Polyakov loop susceptibility is essentially better separated from the rest of the data for the monopole Polyakov loop than for nonabelian Polyakov loop. It is important that the transition temperature values determined by both susceptibilities are the same.

## 3 Results and Conclusions

In Fig.1 we show results for average of various kinds of Polyakov loops. One can see that $`P`$ is a smooth function of $`\kappa `$ and that abelian and monopole Polyakov loops behavior is qualitatively the same as behavior of the nonabelian Polyakov loop while the photon Polyakov loop is almost constant across the transition.

Transition point $`\kappa `$, denoted as $`\kappa _t`$, has been determined from the maximum of the Polyakov loop susceptibility, see Fig 2. As can be seen from Fig.2 the susceptibilities for abelian and monopole Polyakov loops have maxima at the same value of $`\kappa `$. We found that $`\kappa _t=0.1344(1)`$ at $`\beta =5.2`$ and $`\kappa _t=0.1341(1)`$ at $`\beta =5.25`$, see Fig.2. This was transformed into transition temperature with the help of the interpolation formula for $`r_0/a`$ . The results are: $`T_cr_0=0.54(2)`$ and $`0.56(2)`$. In physical units, taking $`1/r_0=394`$ MeV, we obtained $`T_c=213(10)`$ and $`222(10)`$MeV, respectively. Using again data from  we estimate $`m_\pi /m_\rho =0.78,0.82`$ at the transition points.

In Fig.3 we show our results for transition temperature in comparison with results of refs.  and . The ratio $`T_c/\sqrt{\sigma }`$ was computed for $`\sqrt{\sigma }=425`$ MeV used in . We conclude that our results are in good qualitative agreement with previous results . This agreement implies that the dependence of the transition temperature on the lattice spacing is rather weak.

For our lattice with $`N_s/N_t=2`$ the question of the finite volume effects is very important. To check this effect we made simulations at $`\beta =5.2,\kappa =0.1343`$ on $`N_s=24`$ lattice. The $`\kappa `$ value was chosen close to the transition point where finite volume effects should be more pronounced. We found that both average of the Polyakov loop and its susceptibility are the same within error bars as in our main simulations on $`N_s=16`$ lattice. This implies that the finite volume effects should not spoil our conclusions. We are planning to perform simulations on $`24^310`$ lattice at $`m_\pi /m_\rho <0.7`$ to determine $`T_c`$ closer to the chiral limit.

Acknowledgements
This work is partially supported by grants INTAS-00-00111, RFBR 02-02-17308, RFBR 01-02-117456, RFBR 00-15-96-786 and CRDF award RPI-2364-MO-02. We are very obliged to the staff of the Joint Supercomputer Center at Moscow and especially to A.V. Zabrodin for the help in computations. We thank Ph. De Forcrand, V. Mitrjushkin, M. Müller-Preussker for useful discussions. M.Ch. is supported by JSPS Fellowship grant No. P01023. V.B. acknowledges support by JSPS during period of work on this project.


text:
# Quantum state tomography of dissociating molecules

## Abstract

Using tomographic reconstruction we determine the complete internuclear quantum state, represented by the Wigner function, of a dissociating I2 molecule based on femtosecond time resolved position and momentum distributions of the atomic fragments. The experimental data are recorded by timed ionization of the photofragments with an intense 20 fs laser pulse. Our reconstruction method, which relies on Jaynes’ maximum entropy principle, will also be applicable to time resolved position or momentum data obtained with other experimental techniques.

Quantum state tomography derives its name from the tomographic technique in medical diagnostics by which three-dimensional images of the inner parts of an object or a person can be derived from two-dimensional NMR or X-ray pictures obtained from different directions, a technique for which Cormack and Hounsfield were awarded the Nobel prize in physiology and medicine in 1979. In quantum physics, we often deal with a phase space, characterizing the position and velocity distribution of particles, and the aim of quantum state tomography is to determine this distribution from only position and momentum measurements on the particles. As an example, the oscillatory motion of a particle in a quadratic potential is described as a simple rotation in phase space, and hence the mathematical Radon-transform technique Natterer , used in physiology and medicine can be used to extract the phase space distribution of such a particle if only position measurements are made sufficiently often over a time corresponding to the oscillation period. This technique has been succesfully applied to analyze bound molecular states, where the internuclear distance oscillates harmonically Whamsley , and it has been used for non-classical states of light Raymer ; Mlynek , which are formally described as harmonically trapped particles. A number of extensions have been made to incorporate anharmonicities in the trapping potential, and to use mathematical reconstruction techniques which match specific detection schemes Leonhardt ; Vogel ; Wineland ; Haroche . In the present work we determine the internuclear quantum state of dissociating $`I_2`$ molecules. The atomic fragments are created in a state with positive energy, and they fly apart with no confining force between them. The situation thus differs from the one of trapped motion, and it calls for a new theoretical approach.

Experimentally, $`I_2`$ molecules in a molecular beam are photodissociated by irradiation with a 100-fs-long pump pulse centered at 488 nm (see Fig. 1). To determine the internuclear separation (position) and velocitiy (momentum) of the dissociating molecules at time $`\tau `$ after excitation the atomic $`I`$ fragments are ionized by a 20-fs-long probe pulse, and the velocities of the $`I^+`$ ions are measured by a 2-dimensional ion detector skovsen . If only one of the iodine atoms of a dissociating molecule is ionized the $`I^+`$ velocity is equal to half of the internuclear velocity, approximately 4.0 Å/ps. On the 2-D image, shown in Fig. 1, these ions constitute the innermost pair of half rings ($`I^+`$-$`I`$). If, instead, both iodine atoms are ionized the velocity of the $`I^+`$ ions is increased due to their internal Coulomb repulsion. On the 2-D image these ions constitute the outermost pair of half rings ($`I^+`$-$`I^+`$). Using Coulomb’s law the $`I^+`$ velocity distribution from this channel gives directly the internuclear distribution skovsen .

The atomic density is sampled on a discrete position grid with $`N`$ grid points, where $`N=51`$ in the example shown in Fig. 2. For the data presented here, the reconstruction is based on position distributions recorded at $`\tau `$ = 2, 3, 4, 5 picoseconds (ps) and the momentum distribution (which is identical at the four different delays). In quantum tomography on harmonic oscillator states the fast and slow velocity components are slushing back and forth, and the time dependent position distribution maps the phase space distribution from all angles. In free expansion, slow particles will not be able to catch up on fast ones, and we have access to only sheared pictures of the distribution where the high velocity components are shifted towards larger displacements than the low velocity components. This implies that a range of angular viewing angles is missing, and we cannot apply the Radon transform to reconstruct the phase space distribution.

Furthermore, the data-sets are noisy due to the finite counts in different position intervals, and our task is to establish the most reliable estimate of $`\rho `$ in best possible (rather than perfect) agreement with the observed data. Following Drobny and Buzek Drobny , we apply Jaynes’ maximum entropy principle which at the same time provides the smallest deviation from the measured data and the largest possible von Neumann entropy Drobny ; Jaynes : $`S(\rho )=Tr(\rho \mathrm{ln}\rho )`$. Drobny and Buzek have used the maximum entropy principle for quantum state tomography on harmonically trapped atoms detected at too few instants of time to allow a normal reconstruction Drobnyatom , and we have reformulated their method to deal with free atoms.

To calculate the von Neumann entropy, one may use the familiar expression in terms of the eigenvalues $`p_j`$ of the density matrix, $`S=_jp_j\mathrm{ln}p_j`$, but as it was shown by Jaynes Jaynes , it is not necessary to compute the entropy of $`\rho `$ in order to use the maximum entropy principle: The density matrix with largest entropy which conforms best with measured values $`a_i`$ for a set of observables $`A_i`$, can be written on the form

$$\rho =\frac{1}{Z}\mathrm{exp}(\underset{i}{}\lambda _iA_i),$$
(1)
where $`\lambda _i`$ are variables that have to be adjusted to satisfy the agreement with the measured data, and $`Z`$ normalizes the density matrix to unit trace.

A position distribution is the set of expectation values of projection operators on position eigenstates $`|xx|`$. We determine the position distribution at different times, and we need a formal representation of the projection operator $`|x,tx,t|`$ where $`|x,t`$ is a position eigenstate at time $`t`$. For free particles the momentum $`p=\mathrm{}k`$ is a conserved quantity, and it is convenient to represent the projection operators in a momentum representation, since momentum eigenstates accumulate only a trivial phase factor with time, $`|k,t=\mathrm{exp}(i\mathrm{}k^2t/2M)|k,t=0`$:

$`|xx|={\displaystyle 𝑑k_1𝑑k_2|k_1k_1|xx|k_2k_2|}`$
$`={\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑k_1𝑑k_2\mathrm{exp}(i(k_2k_1)x)|k_1k_2|}.`$ (2)
Both the spatial distribution at various times and the momentum distribution, which is independent of time, are measured, and the observables $`A_i`$ in Eq.(1) are precisely the operators $`|x,tx,t|`$ and $`|kk|`$, which are now given explicity as $`N`$ by $`N`$ matrices in a discrete basis. The exponential in Eq.(1) is understood as a matrix power series, and it has to be evaluated numerically because the momentum projection operators and the position projection operators at different times do not commute.

In our numerical reconstruction we used four position distributions and one momentum distributions, and we thus identified $`551`$ variables $`\lambda _i`$ so that the expectation values $`A_i=`$Tr$`(\rho A_i)`$ are as close as possible to the measured distributions. We quantify the agreement with the measured data $`a_i`$ by the sum of the squares of all deviations $`\mathrm{\Delta }=_i(a_iA_i)^2`$. A numerical routine identifies the global minimum of $`\mathrm{\Delta }`$ as a function of the $`\lambda _i`$ variables, and we thus obtain the explicit expression for the density matrix in Eq.(1).

To illustrate our reconstruction of the internuclear quantum state of the dissociated molecules we use the Wigner function $`W(x,p)`$ Hillery  rather than the N by N complex density matrix, discussed above. The Wigner function is obtained by a Fourier transform with respect to the difference between the two momentum arguments in the density matrix. While still fully characterizing the quantum state the Wigner function is a real function, and it has the advantage of presenting at one glance the joint position and momentum distribution. Figure 2 shows the Wigner function, $`W(x,p)`$, at the earliest (2 ps) and the latest time (5 ps) for which position distributions were recorded. Note that we present the momentum dependence in units of the velocity $`v=p/M`$ of the atomic fragments to facilitate comparison of the two panels. The two Wigner functions in Fig. 2 a) and b) have the same velocity distributions, but they are sheared differently in phase space due to the free motion of the two iodine photofragments. On the side panels the marginal position and velocity distribution, obtained by integrating $`W(x,p)`$ over velocity and position, respectively, are compared with the experimental results (open circles). The agreement is excellent, and similar good agreement is found between the measured data and the reconstructed state at the intermediate times at 3 and 4 ps.

We stress that the Wigner function contains more information than what can be obtained directly from a single measured position and velocity distribution. For instance, both the measured position and velocity distribution at 2 ps consists of two peaks (side panels of Fig. 2 (A)), but without further information it is not possible to determine if the molecules at the larger internuclear separation (the peak around 13.6 Å ) are also the ones moving with the largest internuclear velocity (the peak around 4.6 Å/ps). The tomographic reconstruction of the Wigner function makes use of our knowledge of the position distributions measured at 3, 4 and 5 ps, and, indeed, the two separated peaks in $`W(x,p)`$ at 2 ps show that the iodine atoms roughly 13.6 Å apart have a velocity distribution centered around 4.6 Å/ps whereas the iodine atoms apparoximately 12.4 Å apart have a velocity distribution centered at 4.0 Å/ps. The physical origin of the two internuclear velocity components is the incoherent population of v = 0 and v = 1 vibrational levels of the electronic ground state of the I2 molecules prior to dissociation by the pump pulse skovsen .

It has been emphasized in the literature Hillery  that the Wigner function is not a real probability distribution because, due to complimentarity and Heisenbergs’s uncertainty relation, it is not meaningful to assign precise values for both the position and velocity of a particle. This has as its most striking consequence that states exist for which $`W(x,p)`$ attains negative values, but these negative values always occur in phase space regions with area less than Planck’s constant $`\mathrm{}=h/2\pi `$. Dissociation of molecules by two phase-locked femtosecond laser pulses will produce photofragments in a coherent superposition of two localized wave packets skovsen , and we are currently working on improvement of the experiments to enable reconstruction of the Wigner function with a sufficient resolution to display the negativities expected in this case. Such an experiment constitutes a fs time resolved analogue to double-slit atom interferometer studies Mlyneknature .

Quantum state tomography has been established in quantum optics as a remarkable diagnostics tool that has provided clear illustrations of the experimental ability to control and manipulate selected states of fundamental quantum systems such as a single field mode and single ions and atoms Mlynek ; Vogel ; Drobnyatom ; Wineland ; Haroche . Our experimentally reconstructed density matrix or Wigner function provides the complete information about the non-trivial quantum state of the internuclear motion of the dissociating $`I_2`$ molecules. Although our experimental method is limited to small molecules, we believe that quantum state tomography will be applicable also to study chemical reactions in large molecules through time resolved position data obtained, e.g., by the emerging femtosecond electron and X-ray diffraction techniques Ihee ; Reutze ; Corkum ; Service . This will open for comparisons with theoretical calculations at a much more detailed level than probability distributions of single observables, which have formed the basis of comparison so far in femtosecond time resolved chemical reaction dynamics. Finally, we note that modern femtosecond laser technology enables controlled shaping of both the electronic and the atomic structure of molecules through irradiation with sequences of tailored laser pulses Rabitz ; Levis ; Gerber . To fully exploit the capabilities of such quantum manipulation it is crucial to completely characterize the molecular quantum states formed.

We acknowledge the support from the Carlsberg Foundation and The Danish Natural Science Council.


text:
# Conditions of Dynamical Stability for the HD 160691 Planetary System

## 1 Introduction

In our previous paper (Kiseleva-Eggleton et al. 2002) we applied the new technique invented by Cincotta & Simó (2000) and called MEGNO (the acronym of Mean Exponential Growth factor of Nearby Orbits) to a wide neighbourhood of orbital parameters determined using standard two-body Keplerian fits for the recently discovered multi-planetary systems HD 12661, HD 38529, HD 37124, and HD 160691 in order to distinguish between regular and chaotic planetary configurations. We showed that the currently announced orbital parameters place these systems in very different situations from the point of view of dynamical stability. While HD 38529 and HD 37124 are located within large stability zones in the phase space around their determined orbits, the orbital parameters of the HD 12661 planets are located in a border region between stable and unstable dynamical regimes, so while its currently determined orbital parameters produce stable orbits, a minor change within the margin of error of just one parameter may result in a chaotic dynamical system. The orbits in HD 160691 (Jones et al. 2002) at first appeared highly unstable, but using MEGNO we were able to identify a few stability zones in a parameter space which included the parameters not determined from observations, such as the relative inclination $`i_r`$ between the two planetary orbits and the longitudes of the ascending nodes $`\mathrm{\Omega }`$. All these stable configurations are associated with the 2:1 mean motion resonance. The present paper is wholly devoted to a detailed and complete dynamical analysis of the HD 160691 planetary system by taking into account the angular orbital parameters not constrained by observational data ($`i_r`$, $`\mathrm{\Omega }`$), as well as $`\mathrm{sin}i_l`$ (the line-of-sight inclination factor) and the resulting different planetary masses. We have also explored the space of the mean anomalies $`M`$ of the two planets taking initially into account the time of periastron passage $`\tau _{per}`$ given by Jones et al. 2002 (see Table 1).

In one of our earlier papers (Goździewski et al. 2002) we have clearly identified (using MEGNO) the exact location of the 2:1 mean motion resonance and its width for the Gliese 876 planetary system. A recent study by Hadjidemetriou (2002) of periodic orbits in this resonance predicts stable and unstable configurations of planetary systems depending on the hierarchy of planetary masses and eccentricities. The Gliese 876 system where $`m_b<m_c`$ and $`e_b>e_c`$ is, according to this study, a stable configuration. In contrast, the hierarchies of the HD 160691 system are inverse, i.e. $`m_b>m_c`$ and $`e_b<e_c`$. Hadjidemetriou (2002) found that a planetary system at the 2:1 resonance where the inner planet is much more massive than the outer planet is unstable for all values of the eccentricities. However, in the HD 160691 system the two planetary masses are comparable and this leaves the question about its stability open.

In this work we have explored the parameter space available for planets in HD 160691 in order to determine the stability conditions for this system. We have to notice that the orbital parameters of both planets are rather speculative due to the insufficient amount of observations, and even the existence of the second planetary companion is not yet fully confirmed (Jones et al. 2002, Butler et al. 2001). Nevertheless, the mechanism we present in this paper is probably the key to the existence of planetary systems like HD 160691. We have clearly identified the exact location of the 2:1 mean motion resonance, its width, and the secular resonance in apsidal longitudes preserving the stability related to the mean motion resonance.

## 2 Method

A classical method that allows us to distinguish between regular and chaotic dynamical states is the method of Lyapunov Characteristic Numbers (LCN). The estimation of LCN usually requires computations over long evolutionary times, sometimes much longer than the lifetime of the system studied. Let us note that chaotic in the Poincaré sense means that the dynamical behavior is not quasi-periodic (the conventional definition useful for conservative dynamical systems), and does not necessarily mean that the system will disintegrate during any limited period of time. Let us state in addition that we use the property of stability in the Poisson sense : stability is related to the preservation of a certain neighbourhood relative to the initial position of the trajectory. In conservative systems, quasi-periodic orbits remain always confined within certain limits; in this sense they are stable.

In the present work, we use two different methods in order to identify the dynamical state of the HD 160691 system : the MEGNO and FLI (the acronym of Fast Lyapunov Indicator) techniques. These two methods converge faster and are more sensitive than the LCN technique.

FLI is the method introduced by Froeschlé et al. (1997), permitting us to distinguish qualitatively between regular and chaotic motion in a dynamical system (see for example Pilat-Lohinger & Dvorak 2002). MEGNO is a new method developed by Cincotta & Simó (2000) that we have already successfully applied to the study of dynamical stability of extrasolar planetary systems (see e.g. Goździewski et al. 2001, 2002, Kiseleva-Eggleton et al. 2002). This method provides relevant information about the global dynamics and the fine structure of the phase space, and it yields simultaneously a good estimate of the LCN with a comparatively small computational effort (Cincotta & Giordano 2000). The MEGNO is an alternative technique that proves to be efficient for investigation of both ordered and stochastic components of phase space (Cincotta, Giordano & Simó 2002). It provides a clear picture of resonance structures, location of stable and unstable periodic orbits, as well as a measure of hyperbolicity in chaotic domains (i.e. the rate of divergence of unstable orbits) which coincides with that given by the Lyapunov characteristic number.

## 3 Stability conditions

In order to identify the different dynamical behaviors in the parameter space, we use MEGNO which provides the exact location of stable and unstable orbits as well as a measure of hyperbolicity. Figures 1 and 2 show the dynamical state of the HD 160691 system as a function of both orbital semi-major axes $`a_b`$ and $`a_c`$ (Fig. 1) and eccentricities $`e_b`$ and $`e_c`$ (Fig. 2). All other orbital parameters were taken from Table 1 with the addition of values of $`\mathrm{\Omega }_{(b,c)}`$ and $`i_{(b,c)}`$ undetermined from observations, and with $`\tau _{per}`$ replaced by the corresponding calculated values of mean anomalies $`M=(2\pi /P)(t\tau _{per})`$ (with an initial time of reference $`t=\tau _{per}^c`$, one obtains $`M_b=156.3^{}`$ and $`M_c=0^{}`$). In order to avoid a dynamical behavior bound to the plane, the two initial orbital inclinations $`i_b`$ and $`i_c`$ are taken slightly different ($`0^{}`$ and $`1^{}`$ respectively); due to gravitational interactions of the 3-body problem, the relative inclination $`i_r`$ is free to evolve in the 3-D space. In Fig. 2 as well as in other Figures in this paper $`a_c`$ was taken to be 2.381 AU according to the Table footnote. In all our Figures the intersection of horizontal and vertical lines indicates the ‘observational’ initial parameters taken from Table 1. By using the MEGNO indicator characteristic value $`<Y>`$ (Cincotta & Simó 2000), filled and open circles in Fig. 1a (and all other Figures plotted in 2-D in this paper) indicate stable orbits ($`<Y>=2\pm 3\%`$, and $`<Y>=2\pm 5\%`$ respectively), while small dots (not surrounded by circles) indicate highly unstable orbits ($`<Y>2`$). One can see that the system is globally highly unstable, except for a clearly marked stability region. Fig. 1a is a cross section of Fig. 1b (plotted in 3-D) in the plane $`<Y>=2`$ with projections of the different points according to their nature. In Fig. 1b, the peaks indicate the magnitude of instability. The stability strip in Fig. 1a corresponds to a stability valley in 3-D (Fig. 1b). However, its fine structure cannot be seen in this 3-D graph because of the scale factor in the $`z`$ axis due to the high magnitude of the instability peaks.

The stability strip is associated with the 2:1 mean motion resonance. Generally, the presence of a resonance structures the phase or parameter space in different regions, regular and irregular (i.e. chaotic). In Fig. 1b, such different regions appear clearly. The regions on both sides of the stability valley present an irregular distribution of the values of $`<Y>`$ meaning very different magnitudes of instability. The valley on the other hand presents a regular picture at the used scale, which is confirmed by the cross section (a) of (b). The strip indeed obtained in Fig. 1a, composed of filled and open circles, is sufficiently homogeneous and dense to be a stability zone. It is rather wide ($`0.1`$ AU in both semi-major axes), and besides one order of magnitude larger than the one obtained for the Gliese 876 planetary system (also in 2:1 mean motion resonance, see Goździewski, Bois, & Maciejewski, 2002). However, for HD 160691 this stability valley may completely disappear for a small change in angular elements. The stability zone is indeed only permitted for a particular geometrical configuration of the orbits (determined by the combination of the elements $`\mathrm{\Omega }`$, $`\omega `$, and $`i_r`$) combined with particular relative positions of the two planets on their orbits determined by the $`M`$ elements. For example, for a shift only on the $`M_b`$ element without losing the 2:1 mean motion resonance, close approaches at the periastron occur in the end, and consequently the dynamical behavior becomes totally and highly unstable as it is shown in Figure 3 (computed with $`M_b=0^{}`$ instead of $`156.3^{}`$). The same situation appears with a comparable shift on the $`\stackrel{~}{\omega }`$ elements. For instance, with $`\stackrel{~}{\omega }_b=0`$, 100, or $`200^{}`$ instead of $`320^{}`$, the stability valley is lost. The whole stabilizing mechanism depends on particular combinations of the elements $`\mathrm{\Omega }`$, $`\omega `$, $`i_r`$, and $`M`$. We will discuss these stability conditions in detail below.

Figure 2 shows two stability regions in the $`[e_b,e_c]`$ parameter space : (i) a small one with few points where both eccentricities are small, and (ii) a much larger region for high values of the outer orbit ($`e_c>0.5`$). Let us notice that in the latter, the eccentricity of the inner orbit $`e_b`$ is relatively unimportant when $`e_c>0.7`$. Let us emphasize that the large region with a sizeable number of stable points decreases significantly when $`a_c`$ is taken to be 2.3 instead of 2.381 AU (see Fig. 3d in our previous paper Kiseleva-Eggleton et al. 2002). In the present Figure 2, due to the resolution of the grid, the dynamical behavior of the‘observational’ point in eccentricities (Table 1) does not appear. However, with a higher resolution (i.e. more integration points), computations confirm the existence of a stable point at the intersection of the horizontal and vertical lines. In addition, few unstable points are also very close. Let us recall that some zones (notably on the sides of stability structures) may present a sensitiveness to initial conditions for stable versus unstable points in a very close vicinity between them. This situation reminds us of the dynamical state of the HD 12661 planetary system (Kiseleva-Eggleton et al. 2002).

Figure 4 shows how the dynamical state of the HD 160691 system depends on the relative inclination $`i_r`$ between the planetary orbits and on the angular elements $`\omega `$ and $`\mathrm{\Omega }`$ for both planets $`b`$ and $`c`$. In the general case of a 2-planet system, the relative inclination $`i_r`$ is related to the longitudes of ascending nodes $`\mathrm{\Omega }_b`$ and $`\mathrm{\Omega }_c`$ as follows :

$$\mathrm{cos}i_r=\mathrm{cos}i_b\mathrm{cos}i_c+\mathrm{sin}i_b\mathrm{sin}i_c\mathrm{cos}(\mathrm{\Omega }_c\mathrm{\Omega }_b)$$

However, in our computations, for a variation of $`i_b`$ or $`i_c`$, the inclination of the other orbit is taken equal to zero. As a consequence, $`i_r=i_b`$ or $`i_c`$, and then the inclinations on the vertical axes of the Figure 4 may be simply read as a scanning of the relative inclination $`i_r`$ between the orbits. It is important to remark that the distributions of structures are relatively similar in the pairs $`[\omega _b,i_b]`$ and $`[\omega _c,i_c]`$, as well as in $`[\mathrm{\Omega }_b,i_b]`$ and $`[\mathrm{\Omega }_c,i_c]`$, with stable orbits appearing in strips of a comparable width $`100^{}`$ around the nominal values of $`\omega _b`$ and $`\omega _c`$ (Table 1) for all values of $`i_r`$. There are also additional stability zones, but very narrow, for small values of $`i_r<20^{}`$ in both maps related to $`\omega `$, as well as in the two others related to $`\mathrm{\Omega }`$. The nominal values for $`\omega _b`$ and $`\omega _c`$, namely 320 and 99 degrees respectively, are located well inside the stable strips.

Figures 5 and 6 present the stability maps of the $`[\stackrel{~}{\omega }_b,\stackrel{~}{\omega }_c]`$ (where $`\stackrel{~}{\omega }=\mathrm{\Omega }+\omega `$) and $`[M_b,M_c]`$ parameter spaces, respectively. Fig. 5 shows a linear relationship between the longitudes of periastron $`\stackrel{~}{\omega }_b`$ and $`\stackrel{~}{\omega }_c`$ for stable configurations. This means that $`\stackrel{~}{\omega }_b`$ and $`\stackrel{~}{\omega }_c`$ precess at the same rate. Fig. 6 shows stable strips in the mean anomaly $`[M_b,M_c]`$ plane that permit us to identify adequate relative orbital positions for stable configurations at the 2:1 resonance in mean motions. Figures 5 and 6 also show that the HD 160691 planetary system is dynamically stable only when the corresponding orbital parameters are simultaneously represented by stable points inside the strips in both $`[\stackrel{~}{\omega }_b,\stackrel{~}{\omega }_c]`$ and $`[M_b,M_c]`$ planes. If only one of the two conditions is realized, the map $`[a_b,a_c]`$ exhibits a totally unstable dynamical state without any particular stability structure (see Fig. 3). We have also computed the $`[\stackrel{~}{\omega }_b,\stackrel{~}{\omega }_c]`$ stability maps for different relative inclinations (from 0 to $`90\mathrm{deg}`$, not shown in this paper), and consequently we may point out that the whole mechanism of stabilization can be successfully applied even to high relative inclinations if $`\mathrm{\Omega },\omega ,M`$ are determined in an appropriate way : $`\stackrel{~}{\omega }_b`$ and $`\stackrel{~}{\omega }_c`$ on average precess at the same rate while the planets b and c, in 2:1 mean motion resonance, are located on their orbits in such a way that there are no close approaches at their periastron. The favourable respective positions are defined by the stable strips in the mean anomaly $`[M_b,M_c]`$ plane (cf. Fig. 6).

The stability requires a particular geometrical configuration of the orbits (defined by elements $`\mathrm{\Omega }`$, $`\omega `$, and $`i_r`$), added to particular positions of the two planets on their orbits defined by the $`M`$ element. In addition to the 2:1 mean motion resonance, the planets b and c have to be initially located in favourable respective positions, for example at their apoapses, in order to avoid close approaches at their periapses. Following analogous properties analyzed by Lee & Peale (2002) for the GJ 876 planetary system, we have found in the HD 160691 system the simultaneous librations of the two mean motion resonance variables $`\theta _1=\lambda _b2\lambda _c+\stackrel{~}{\omega }_b`$ and $`\theta _2=\lambda _b2\lambda _c+\stackrel{~}{\omega }_c`$ (where $`\lambda =M+\stackrel{~}{\omega }`$), while the secular resonance variable

$$\theta _3=\theta _1\theta _2=\stackrel{~}{\omega }_b\stackrel{~}{\omega }_c$$

librates about $`180^{}`$ with a period of about 2800 years and an amplitude of 80 degrees (see Fig. 7). $`\theta _3`$ being only librating (i.e. without secular term), the two orbital planes on average precess at the same rate. This is a mechanism of secular resonance in periastron, which as a consequence ensures the maintenance of the orbital topology as well as the respective motions of the planets, that is to say in the present case, the 2:1 mean motion resonance without close approaches of the planets at their periapses. In other words, this mechanism, when combined with adequate relative positions of the planets on the initial osculating orbits, can be understood as the condition for the preservation of the dynamical stability related to the 2:1 mean motion resonance. Because the two variables $`\theta _1`$ and $`\theta _2`$ librate around $`180^{}`$ and $`0^{}`$, respectively (see Fig. 7), the lines of apsides are anti-aligned. This situation is different from the GJ 876 configuration in which the apsidal lines are aligned with $`\theta _1`$, $`\theta _2`$, and $`\theta _3`$ all librating about $`0^{}`$ (Lee & Peale 2002, Goździewski, Bois & Maciejewski 2002). Let us note that in the familiar Io-Europa 2:1 resonance, the very small eccentricities lead to a geometry where $`\theta _1`$ is librating about $`0^{}`$, while $`\theta _2`$ and $`\theta _3`$ are librating about $`180^{}`$ (Lee & Peale 2002). In this case, conjunctions occur when Io is near periapse and Europa near apoapse.

The whole stability mechanism for the HD 160691 system allows it to avoid close approaches between planets, especially at their periapses (or near them). If it is not the case, the dynamical behavior is wholly unstable since the close approaches near periapses occur regularly (see Fig. 3 where $`M_b=M_c=0^{}`$). The original orbital topology of the HD 160691 planetary system associated to the 2:1 mean motion resonance can be written as $`P_c(ap)SP_b(ap)`$ 111This 3-D topology with anti-aligned apsidal lines is out of the scope of the various 2:1 resonance configurations analyzed by Hadjidemetriou (2002) (2 types of configuration) and Hadjidemetriou & Psychogiou (in preparation) (4 types in the planar case)., which means that the planets c and b may be considered as initially located at their apoastron around the central star S (state 1). Taking into account the 2:1 mean motion resonance, the present topology is equivalent to $`SP_c(per)P_b(ap)`$ : after one revolution of the planet b, the planet c is at its periapse (state 2). This situation corresponds (by taking into account the important amplitude of the secular resonance variable $`\theta _3`$, namely $`80\mathrm{deg}`$) to the nominal anomalies given in Table 1 ($`M_b160\mathrm{deg},M_c=0`$). After one revolution more of the planet b, the planets are again near their apoastron (state 1). Besides, with the 2nd state, we understand that the present topology may drive to a solid stability where the sum of both eccentricities is particularly important (i.e. $`e_b+e_c>0.7`$, high $`e_c`$ helps) (see Fig. 2). Because of the high eccentricities of the orbits, and despite relatively small semi-major axes, the relative distances between the two planets may remain sufficiently large over the whole evolutionary time scale of the system.

We have also tested the robustness of this mechanism when the values of the planetary masses were progressively increased. In Figure 8, computed by the FLI technique in the $`[a_b,a_c]`$ parameter space, the two masses are consecutively multiplied by a factor $`k=`$ 1 (a), 1.3 (b), 1.4 (c), 1.5 (d). The total destruction of the stable structure is reached for $`k=1.5`$, which corresponds to a line-of-sight inclination factor $`\mathrm{sin}i_l=2/3`$ (i.e. $`i_l=41.8^{}`$), as we know only the lower limit of planetary masses $`m_0`$ from the function $`m_0=m_P\mathrm{sin}i_l`$ where $`m_P`$ is the real mass. This result is in agreement with a numerical analysis by Kiseleva-Eggleton & Bois (2001) of multi-planetary systems of $`\upsilon `$ And-type which showed that the dynamical regime, and consequently the lifetime of the system, depends strongly on the mass hierarchy as well as on the absolute values of planetary masses. In other words, assuming that the observed planetary system is dynamically stable, we could probably determine a window for the possible values of line-of-sight inclinations $`i_l`$, and as a consequence give upper limits on the planetary masses.

## 4 Conclusion

Using the MEGNO technique of global dynamics analysis we scanned the most relevant cases of the orbital parameter space for the HD 160691 planetary system. We have found the existence of a stability zone ruled by a mechanism which involves angular elements of the system. This stability zone is indeed due to the 2:1 mean motion resonance coupled with adequate relative positions of the planets on their orbits avoiding close approaches at their periastron, the two apsidal lines being anti-aligned. The mechanism is not lost during the dynamical evolution of the system due to an apsidal secular resonance : the mean motion resonance variables are librating while the longitudes of periapse on average precess at the same rate.

We conclude that in order to be dynamically stable, the HD 160691 planetary system has to satisfy the following conditions : (1) a 2:1 mean motion resonance, (2) combined with an apsidal secular resonance, (3) in a configuration $`P_c(ap)SP_b(ap)`$ (i.e. an apsidal antialignment), (4) and specific conditions on the respective sizes of the eccentricities. High eccentricity for the outer orbit ($`e_c>0.52`$) is the most probable necessary condition, while the eccentricity of the inner orbit $`e_b`$ becomes relatively unimportant when $`e_c>0.7`$. These four conditions, taking also into account various relative inclinations between the two orbits, determine the dynamical behavior of the system in such a way that the planets are never too close to each other. In the end, the HD 160691 system, where the resonance variables $`\theta _1`$ and $`\theta _3`$ librate about $`180^{}`$ while $`\theta _2`$ librates about $`0^{}`$, has revealed resources of the 2:1 orbital resonances that have not been observed nor analyzed before. The present orbital topology \[$`P_c(ap)SP_b(ap)`$\] combined with the 2:1 orbital resonance is indeed more wide-ranging than the particular case of the HD 160691 planetary system. It is a new theoretical possibility suitable for a stable regime despite relatively small semi-major axes with respect to the important masses in interactions.222Whereas we are finishing the present paper, we learn that an apsidal antialignment has been proposed for the dynamical stability of the HD 82943 planetary system related to a 2:1 orbital resonance (Ji et al. 2003).

Combining our MEGNO maps (confirmed with the FLI method) for different pairs of parameters, it is possible to converge towards stability conditions related to different combinations of the angular parameters. As the observational determination of the elements of the HD 160691 system is far from being finalized, we hope that the maps presented here will be useful for testing both future observations and different parameter fitting techniques. In this respect, let us mention the fitting method of Laughlin & Chambers (2001), suitable for resonant interactions between the planets and where the true masses can be determined by eliminating the indeterminacy in $`\mathrm{sin}i_l`$ inherent in fits that assume independent Keplerian motions. Let us note here that we have already tested some of the new fits for HD 160691 obtained by Eugenio Rivera (private communication), and found that the only stable systems with $`k<2`$ that do not meet our stability conditions would be ones with a very massive substellar distant companion ($`m_c50m_J`$) on a very large orbit ($`a_c30`$ AU). This type of system is very different from the systems with Jupiter-mass planets on close orbits that we have discussed in this paper.

LK-E thanks the University Bordeaux 1 for a short-term research fellowship in Bordeaux and John Hadjidemetriou for the text of his yet unpublished paper. EB and LK-E thank Rudolf Dvorak and the Institute of Astronomy, University of Vienna, for hospitality. EP-L wishes to acknowledge the support by the Austrian FWF (Hertha Firnberg Project T122). The authors thank Eugenio Rivera for providing us with his latest unpublished orbital fits, and for useful discussions. We thank Peter Eggleton for his comments and English corrections, and the referee Georges Contopoulos whose questions and comments greatly helped to improve the paper.


text:
# Nonlinear Thomas-Fermi-Poisson theory of screening for a Hall bar under strong magnetic fields

## 1 Introduction

As a consequence of the highly degenerate, Landau-quantised energy levels, a two-dimensional electron gas (2DEG) in a strong perpendicular magnetic field has unusual low-temperature screening properties . In an inhomogeneous 2DEG with long-range density fluctuations, pinning of Landau levels (LLs) leads to quasi metallic regions with high density of states (DOS) at the Fermi energy and therefore nearly perfect screening ability. These so called “compressible” regions coexist with quasi-insulating “incompressible” regions separating compressible regions with different LLs at the Fermi energy. In the incompressible regions the Fermi energy falls into the gap between two adjacent LLs and, in thermodynamic equilibrium at low temperatures, no redistribution of electrons is possible which could contribute to screening. Also, in the incompressible regions the electron density $`n_{el}(𝐫)`$ is constant, corresponding to total filling of an integer number of LLs, while in the compressible regions $`n_{el}(𝐫)`$ adjusts itself so that the self-consistent electrostatic potential energy $`V(𝐫)`$ of an electron differs from the Fermi energy, more precisely from the constant electrochemical potential $`\mu ^{}`$ by a Landau energy $`\mathrm{}\omega _c(n+1/2)`$, where $`\omega _c=eB/m`$ is the cyclotron frequency in the magnetic field $`B`$. As a consequence, $`V(𝐫)`$ differs between different compressible regions by integer multiples of $`\mathrm{}\omega _c`$. Landau level pinning and the interplay of compressible and incompressible regions will strongly affect the screening properties of the 2DEG.

In a bounded 2D geometry with translation invariance in one direction, the compressible and incompressible regions degenerate to strips parallel to the boundary. For half-space and Hall-bar geometries models with planar charge distributions have been proposed that allow closed solutions of Poisson’s equation, i.e., the calculation of the potential for given electron density, and estimates of position and widths of the incompressible strips have been given . By adding the non-linear Thomas-Fermi approximation for the calculation of the electron density from the potential, that work was extended to a self-consistent approach which allows to calculate both electron density and electrostatic potential for arbitrary temperature . This approach shows that the existence and the width of incompressible strips depends sensitively on temperature. It also permits to calculate their position and width for given background charges without additional assumptions. We adopt this approach to investigate systematically the screening of a harmonic external potential in such a confined 2DEG.

## 2 The Thomas-Fermi-Poisson (TFP) approach

In order to employ the approach presented in Ref. , we assume that our 2DEG lies in the $`z=0`$ plane, laterally confined by in-plane gates located at $`x<d`$ and $`x>d`$, which are kept at voltages $`V_L`$ and $`V_R`$, respectively. In this work we assume a symmetric depletion on both sides with the depleted strips having a width of $`|db|`$, and we will consider only symmetric potentials $`V(x)=V(x)`$, notably with $`V(\pm d)=V_L=V_R=0`$, which leave the electron density profile symmetric, $`n_{el}(x)=n_{el}(x)`$. Then, solution of Poisson’s equation yields for the electrostatic potential energy of an electron

$`V(x)={\displaystyle \frac{2e}{\overline{\kappa }}}{\displaystyle _d^d}𝑑tK(x,t)\rho (t),K(x,t)=\mathrm{ln}\left|{\displaystyle \frac{\sqrt{(d^2x^2)(d^2t^2)}+d^2tx}{(xt)d}}\right|,`$ (1)
where $`\rho (x)`$ is the surface charge density in the strip $`|x|<d`$ and $`\overline{\kappa }`$ is an average dielectric constant . If we assume a homogeneous surface density $`n_0`$ of positive background charges, $`\rho (x)=e[n_0n_{el}(x)]`$, we obtain

$`V(x)=V_{bg}(x)+V_H(x),\text{where}V_{bg}(x)=E_0\sqrt{1(x/d)^2}`$ (2)
with $`E_0=2\pi e^2n_0d/\overline{\kappa }`$ is the confining potential due to the background and $`V_H(x)`$ the Hartree contribution due to the electron density, given by Eq. (1) with $`\rho (x)`$ replaced by $`en_{el}(x)`$. Assuming that $`V(x)`$ varies on a characteristic length much larger than typical quantum lengths, notably the magnetic length $`l_m=\sqrt{(\mathrm{}/eB)}`$, one can calculate $`n_{el}(x)`$ for given $`V(x)`$ in the Thomas-Fermi approximation

$`n_{el}(x)={\displaystyle 𝑑ED(E)f\left([E+V(x)\mu ^{}]/k_BT\right)},\text{with}f(ϵ)=[1+e^ϵ]^1`$ (3)
the Fermi function, $`\mu ^{}`$ the electrochemical potential and $`T`$ the temperature. Here $`D(E)`$ is the density of states of the 2DEG, which we take to be the bare spin-degenerate Landau DOS,

$`D(E)=(\pi l_m^2)^1{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\delta \left(E\mathrm{}\omega _c(n+1/2)\right).`$ (4)
This completes the TFP scheme for the self-consistent calculation of $`n_{el}(x)`$ and $`V(x)`$ for given temperature and magnetic field and homogeneous background charge. To study screening effects, we add to $`V_{bg}(x)`$ in Eq. (2) the harmonic external potential

$`V_m(x)=V_0\mathrm{cos}(k_\lambda x),\text{with}k_\lambda =(\lambda +1/2)\pi /d,`$ (5)
where $`\lambda `$ is an integer (to preserve the boundary conditions), calculate the resulting self-consistent potential and electron density as function of the modulation strength $`V_0`$, and compare these results with those for $`V_0=0`$. \[In principle $`V_m(x)`$ can be generated by a redistribution of the background charge density.\]

For technical reasons we always start the iterative solution of the self-consistent TFP scheme with a calculation for vanishing temperature and magnetic field and for homogeneous background, i.e., we take $`n_{el}(x)=D_0[E_FV(x)]\theta \left(E_FV(x)\right)`$ with $`D_0=m/\pi \mathrm{}^2`$ as 2D DOS and define the symmetric density profile by the requirement $`V(b)=V(b)=E_F=\mu ^{}(T=0)`$. Then Eq. (1) reduces to a linear integral equation,

$`V(x)=E_0\sqrt{1\left({\displaystyle \frac{x}{d}}\right)^2}+{\displaystyle \frac{1}{\pi a_0}}{\displaystyle _b^b}𝑑t\left[E_FV(t)\right]K(x,t)`$ (6)
with $`a_0=\overline{\kappa }/2\pi e^2D_0`$ the bare screening length. \[Due to a misprint, in the corresponding Eq. (15) of Ref.  the factor $`\pi `$ in the denominator is missing.\] The choice of $`b`$ defines the average density, which is then kept fixed for the calculations at finite temperature, magnetic field, and modulation.

## 3 Results and discussion

In the numerical calculations we keep some parameters fixed. So we divide the interval $`d<x<d`$ into $`500`$ subintervals and calculate $`V(x)`$ and $`n_{el}(x)`$ on the corresponding equidistant mesh points $`x=x_n`$. The depletion length is defined by taking $`b/d=0.9`$ and the bar width by $`\pi a_0/d=0.01`$. For GaAs values ($`\overline{\kappa }=12.4`$) $`a_05`$nm, this means $`2d3\mu `$m, which is small compared to typical Hall bars used in experiments. Calculation for smaller $`a_0/d`$, i.e., larger $`d`$ requires more subintervals, and thus larger storage and computation time, to achieve the same accuracy near the incompressible strips, but yields no qualitative changes. Therefore, we fix the mentioned parameter values and obtain from the solution of the system of linear equations defined by Eq. (6) an average electron density that yields the average LL filling factor $`\overline{\nu }=2`$ for a cyclotron energy $`\mathrm{}\omega _c=\mathrm{\Omega }_20.2311\times 10^2E_0`$. Starting from this solution we find the solution of the non-linear TFP scheme at finite $`T`$ and $`B`$ implementing the Newton-Raphson method . For the bare modulation potential (5) we take $`\lambda =2`$, so that $`V_m(x)`$ has a maximum at $`x=0`$ and two minima at $`x=\pm 0.4d`$, i.e. well inside the Hall bar, whereas near the maxima at $`x=\pm 0.8d`$ the confinement potential increases strongly and the electron density starts to decrease, so that the screening properties are obscured by boundary effects.

Figure 1 shows the potentials and electron densities of the Hall bar calculated for several modulation amplitudes at strong magnetic fields and low temperatures. On the left panel we consider the case $`V_0=0`$ and show the temperature dependence in the inset, while the right panel shows the evolution of the potential and electron density with increasing modulation amplitude. All the electron densities are expressed in terms of local filling given by $`\nu (x)=2\pi n_{el}(x)l_m^2`$ and potentials are normalised by $`E_0`$.

Figure 1a presents the electron density, with two incompressible strips located symmetrically near $`x/d`$=$`\pm 0.75`$ and the surrounding compressible regions. With decreasing temperature the incompressible strips are stronger pronounced (see inset Fig. 1a) and so are the related potential steps (inset Fig. 1b). As mentioned in the introduction, in compressible regions one of the LLs is pinned at $`\mu ^{}`$, leading to nearly perfect screening and constant potential, while in an incompressible region $`\mu ^{}`$ falls into a gap between the LLs making redistribution of electrons, and hence screening, energetically impossible and $`n_{el}(x)`$ constant. We demonstrate these features in Fig. 1b, where we exhibit $`\mu ^{}`$ (dash-dotted line), the first two LLs (dashed lines with opaque circles), and the total potential (solid line).

In the right panels we show electron density and potential profiles for various modulation amplitudes $`V_0`$ at a low temperature and a high magnetic field, at which the average filling factor in the centre region is around 2.5, i.e., the $`n=1`$ LL is pinned to $`\mu ^{}`$. With increasing $`V_0`$ the external potential (energy) is raised in the centre and lowered near $`x/d`$=$`\pm 0.4`$, which leads to a lowering of $`n_{el}(x)`$ near $`x=0`$ and an increase near $`x/d`$=$`\pm 0.4`$. This redistribution of charges screens the potential at the low temperature of Fig. 1c,d so effectively, that the resulting weak modulation of the self-consistent potential for $`V_0/E_0=0.01`$ and 0.02 is not seen on the scale of Fig. 1d. At larger modulation, however, the filling factor $`\nu (0)`$ in the centre falls below 2, which means that near $`x=0`$ the self-consistent potential must rise so much that the lowest LL $`n=0`$ reaches $`\mu ^{}`$ (as is demonstrated for $`V_0/E_0=0.07`$ in Fig. 1d). Therefore, on both sides of the centre region incompressible strips must develop across which the self-consistent potential changes by $`\mathrm{}\omega _c`$.

We see from these results that the screening is highly non-linear, since the change of the self-consistent potential produced by the harmonic modulation (5) is step-like rather than cosine-like. Moreover, for weak modulation $`V_0/E_00.02`$ the resulting variation $`\mathrm{\Delta }V=V(0)V(0.4d)`$ of the self-consistent potential is so small that it is not detectable on the scale of Fig. 1d, whereas $`\mathrm{\Delta }V`$ for $`V_0/E_00.04`$ assumes the constant value $`\mathrm{}\omega _c`$. Apparently this stepwise increase of $`\mathrm{\Delta }V`$ with $`V_0`$ is a consequence of the pinning of LLs to the electrochemical potential (locally perfect screening), which is stronger pronounced at lower temperatures. We now consider this effect in more detail.

### 3.1 Tuning the modulation strength at fixed magnetic fields

In Fig. 2a we show the “variance” of the screened potential as function of the amplitude $`V_0`$ of the applied modulation for several values of the magnetic field and a fixed low temperature, and in Fig. 2b the same dependence for a fixed value of the cyclotron energy and several temperatures. We define this variance as the difference between the maximum (at $`x=0`$) and the minimum (near $`x=0.4d`$) of the difference $`V(x;V_0)V(x;0)`$ of the self-consistent potentials calculated with and without modulation, respectively.

For a fixed value of the magnetic field (i.e., of $`\mathrm{\Omega }\mathrm{}\omega _c`$), we obtain a stepwise increasing curve, as we expect from the discussion of Fig. 1c,d, which corresponds to $`\mathrm{\Omega }=\mathrm{\Omega }_2`$. The corresponding (thick solid) curve in Fig. 2a shows two broadened steps around $`V_0/E_0=0.025`$ and 0.075. We have already discussed the reason for the step at the lower $`V_0`$ value and for the plateaus below and above this step. With further increasing $`V_0`$ the filling factor near the density maxima becomes larger than 4 (see Fig. 1c) while the density in the centre is still finite. Then new incompressible strips with local filling $`\nu (x)=4`$ develop on both sides of the density maxima, which corresponds to a decrease of the self-consistent potential by the amount $`\mathrm{\Omega }`$ to new local minima near $`x=\pm 0.4d`$. In Fig. 2a this leads to the plateau of height $`2\mathrm{\Omega }`$. As $`V_0`$ increases further, the electron density $`n_{el}(0)`$ in the centre vanishes and no further redistribution of electrons from the centre region to the region of maximum density is possible. Then the high screening ability of the electron system breaks down and the variance increases with further increasing $`V_0`$ much more rapidly (see behaviour at $`V_0/E_0>0.125`$). This breakdown happens for $`\mathrm{\Omega }=\mathrm{\Omega }_2`$ when the filling factor at the density maxima is somewhat larger than 4 ($`5`$).

We can now understand all the curves in Fig. 2a. Since, apart from the fine structure near incompressible strips, the global density profile depends only weakly on the magnetic field, for the lowest curve (long-dashed, $`\mathrm{\Omega }=3\mathrm{\Omega }_2`$) typical filling factors $`\nu (x)`$ are by a factor 3 smaller than for the case $`\mathrm{\Omega }=\mathrm{\Omega }_2`$ we have just discussed. Therefore the breakdown situation is reached when the filling factor near the density maxima is of the order $`5/3<2`$. Thus the breakdown sets in before an incompressible strip and thus a plateau of height $`3\mathrm{\Omega }_2`$ of the variance can develop. For $`\mathrm{\Omega }=2\mathrm{\Omega }_2`$ \[short-dashed line, label (i)\] the local filling factors $`\nu (x)`$ are typically about half as large as in the case $`\mathrm{\Omega }=\mathrm{\Omega }_2`$. Thus, $`\nu (x)<2`$ for weak modulation, and incompressible strips with $`\nu (x)=2`$ can occur only near the density maxima for sufficiently strong modulation, and a single plateau of height $`2\mathrm{\Omega }_2`$ occurs below the breakdown regime. A similar situation is met for $`\mathrm{\Omega }=1.5\mathrm{\Omega }_2`$ \[label (ii)\], but since now $`\nu (x)`$ is about $`2/3`$ times that for $`\mathrm{\Omega }=\mathrm{\Omega }_2`$, the threshold of the plateau is reached at a considerably weaker modulation. Finally, for $`\mathrm{\Omega }=0.75\mathrm{\Omega }_2`$ we obtain three plateaus in the pre-breakdown regime, corresponding to the successive development of incompressible strips with local filling factor 4 (maxima), 2 (minimum), and 6 (maxima).

We should mention that the behaviour of the variance curves is not always that regular. For special values of $`\mathrm{\Omega }`$ it may happen that incompressible strips near the centre and near the density maxima occur at the same modulation strength. Then steps of height $`2\mathrm{\Omega }`$ occur. We met such a situation for $`\mathrm{\Omega }=0.5\mathrm{\Omega }_2`$ (not shown in Fig. 2a), which yields nearly the same trace as $`\mathrm{\Omega }=\mathrm{\Omega }_2`$.

For the smaller $`\mathrm{\Omega }`$ values, the variance curves oscillate around the dash-dotted line in Fig. 2a, which represents the corresponding result $`2V_0/[E_0\epsilon (q)]`$ of linear Thomas-Fermi screening  for zero $`T`$ and $`B`$, with $`\epsilon (q)=1+1/(a_0q)=41`$ for $`q=2.5\pi /d`$. This linear screening approximation breaks down when the amplitude of the screened potential (i.e. half the variance) becomes equal to the Fermi energy of the unmodulated system , which can be estimated from Fig. 1 as $`3\times 10^3E_0`$, corresponding to $`V_00.12E_0`$. Finally we should mention that the small-variance regime in Fig. 2a also corresponds to linear screening, but now with $`\epsilon (q)=1+(D_T/D_0)/(a_0q)`$, where $`D_T`$ is the thermodynamic DOS in a LL, which can be estimated by $`D_T/D_0\mathrm{}\omega _c/(4k_BT)`$ . Thus, the slope of the variance curves is by a factor $`4k_BT/\mathrm{\Omega }1`$ smaller than that of the dash-dotted line.

As shown in Fig. 2b, with increasing temperature the plateaus of the variance-versus-$`V_0`$ curves get a finite slope and the steps are smeared out. This is an immediate consequence of the temperature dependence of the incompressible strips (see Fig. 1a). At the highest shown temperature screening becomes again linear and independent of $`B`$, but now the temperature is so high that the thermodynamic DOS is noticeably smaller than $`D_0`$. A clear indication of steps is seen only for $`k_BT/\mathrm{}\omega _c1/25`$ .

### 3.2 Sweeping the magnetic field at fixed modulation strength

From the above considerations we expect that, for fixed modulation strength $`V_0`$, the variance as a function of the cyclotron energy $`\mathrm{\Omega }`$ will mainly follow integer multiples of $`\mathrm{\Omega }`$ with abrupt transitions between neighbouring integers corresponding to the steps in Fig. 2a. These abrupt transitions will be similar to the familiar magnetic-field-dependent jumps of the chemical potential of a homogeneous 2DEG at fixed density, which have nothing to do with the periodic modulation and should be separated from the modulation-induced structures. To this end, we first discuss the chemical potential.

In a Hall bar with spatially varying electron density and electrostatic fields the electrochemical potential $`\mu ^{}`$ is a thermodynamic quantity being constant in equilibrium. As an analogue of the chemical potential of a homogeneous 2DEG we define $`\mu (x)=\mu ^{}V(x)`$ and denote in the following $`\mu (0)`$ as the chemical potential. In Fig. 3a we show, together with the LLs $`\mathrm{\Omega }(n+1/2)`$, a plot of $`\mu (0)`$ versus $`\mathrm{\Omega }`$, which exhibits a saw-tooth shape as known from the unbounded 2DEG. To understand this, we plot in Fig. 3b and 3c filling factor and self-consistent potential $`V(x)`$ for four magnetic field values marked in Fig. 3a.

For the lowest selected $`B`$ value \[opaque circle in Fig. 3a, label (i)\], in the centre region the LL with $`n=1`$ is pinned to $`\mu ^{}`$, the filling in the centre is $`\nu (0)>2`$, and the incompressible strips with $`\nu (x)=2`$ are close to the edges. As $`\mathrm{\Omega }`$ increases these strips move towards the centre, and the compressible centre region, with $`\nu (x)>2`$ and a flat potential minimum, shrinks \[case (ii), filled square\]. In the transition region at still larger $`\mathrm{\Omega }`$ \[(iii), opaque diamond\] the strips merge to an incompressible centre region, where the potential minimum is no longer flat, and its depth (measured from the adjacent compressible strips) is smaller that $`\mathrm{\Omega }`$. As $`\mathrm{\Omega }`$ sweeps through the transition region, the width of the incompressible centre region and the depth of the potential minimum shrink to zero. When the transition is completed \[(iv), filled triangle\], we have again a broad compressible centre region in which the next lower LL (here $`n=0`$) is pinned to $`\mu ^{}`$. Thus, in the confined 2DEG in a Hall bar, the jump of the chemical potential from a LL to the next lower one is realised by a drastic change of the position dependence of the self-consistent potential, which in our idealised model happens in the centre.

To separate these effects from the modulation-induced screening effects, we define now the variance of the self-consistent potential as $`var=V(0)V(x_{min})`$, where $`\pm x_{min}\pm 0.4d`$ are the positions of the minima of $`V(x)`$ in the presence of modulation, and $`V(0)`$ is the local maximum at the centre. In Fig. 4a we plot $`var`$ versus $`\mathrm{\Omega }`$, together with $`\mu (0)`$. As expected, $`var`$ follows widely integer multiples of $`\mathrm{\Omega }`$, whereas $`\mu (0)`$ follows half-integer ones. But, whereas $`\mu (0)`$ jumps with increasing $`\mathrm{\Omega }`$ always to the next lower LL, $`var`$ can also jump to a higher multiple of $`\mathrm{\Omega }`$. To understand this, we show in Fig. 4b-4e density and potential profiles for the six $`\mathrm{\Omega }`$ values marked in Fig. 4a.

In case (i) (solid lines in Fig. 4b and 4c) the centre region is compressible with $`\nu (0)`$ slightly larger than 2 and with the $`n=1`$ LL pinned to $`\mu ^{}`$, whereas the compressible regions near the density maxima have $`\nu (x)>4`$ ($`\nu (x_{min})5`$) and the $`n=2`$ LL is pinned to $`\mu ^{}`$. This yields $`var\mathrm{\Omega }`$. For slightly larger $`\mathrm{\Omega }`$ \[(ii), dashed lines in 4b and 4c\], an incompressible strip with filling factor 2 and a local maximum of $`V(x)`$ at $`x=0`$ develops in the centre, while the situation near the density maxima is unchanged. This leads to $`\mathrm{\Omega }<var<2\mathrm{\Omega }`$. At still slightly larger $`\mathrm{\Omega }`$ \[(iii), dotted lines in 4b and 4c\], the region near the density maxima is still qualitatively unchanged, whereas in the centre a compressible region with $`\nu (x)<2`$ develops, where the $`n=0`$ LL is pinned to $`\mu ^{}`$, leading to $`var2\mathrm{\Omega }`$. As $`\mathrm{\Omega }`$ increases further \[(iv) – (vi)\], the situation in the centre remains qualitatively the same, whereas the filling factor at the density maxima changes from $`\nu (x_{min})>4`$ to $`\nu (x_{min})<4`$, so that the incompressible strips with $`\nu (x)=4`$ merge \[(v), dashed lines in 4d and 4e\] and finally disappear, accompanied by the disappearance of the potential step across these strips. Thus, we have again $`var\mathrm{\Omega }`$, until for much larger $`\mathrm{\Omega }`$ the maximum filling factor becomes smaller than 2 and $`var`$ goes back to the very small value corresponding to linear screening in the lowest Landau level.

In summary, the low-temperature screening properties of a confined inhomogeneous 2DEG in strong magnetic fields can be understood from a few basic facts: Incompressible regions with integer values of the local filling factor (even integer values for a spin-degenerate 2DEG), and thus constant electron density, are accompanied by steps of the self-consistent potential between the adjacent compressible regions, where adjacent Landau levels are pinned to the constant electrochemical potential, so that the step height is $`\mathrm{}\omega _c`$, the cyclotron energy. The modulation-induced variation of the self-consistent potential (“variance”) is, therefore, usually an integer multiple of $`\mathrm{}\omega _c`$ and changes by $`\pm \mathrm{}\omega _c`$ if incompressible regions appear or disappear (near local maxima or minima of the electron density) due to a variation of modulation strength or magnetic field. Screening of the applied modulation potential breaks down if the electron density near the density minima becomes so small that no further redistribution of electrons from the local potential maxima to the local potential minima is possible.

## Acknowledgements

We gratefully acknowledge helpful discussions with E. Ahlswede and J. Weis, and financial support by the Deutsche Forschungsgemeinschaft, SP “Quanten-Hall-Systeme” GE306/4-1.

## References


text:
# Normed convergence property for hypergroups admitting an invariant measure

## Abstract

We prove that a hypergroup admitting a countable basis and an invariant Haar measure has normed convergence property if and only if it is compact.

AMS Mathematics Subject Classification 2000: 43A62, 60B15.

Keywords: Hypergroups, probability measures, weak topology and invariant measure.

Let $`K`$ be a (locally compact) hypergroup and $`𝒫(K)`$ be the space of all regular Borel probability measures on $`K`$ equipped with the weak topology. For $`xK`$, $`\delta _x`$ denotes the probability measure concentrated at $`x`$ and for $`\mu ^1(K)`$, $`x\mu `$ and $`\mu x`$ denote $`\delta _x\mu `$ and $`\mu \delta _x`$ respectively. Also, for $`\mu _1,\mu _2𝒫(K)`$, $`\mu _1\mu _2`$ denotes $`\mu _1\mu _2`$.

###### Definition

A hypergroup $`K`$ is said to have normed convergence property if for any sequence $`(\mu _n)`$ in $`𝒫(K)`$, there exists a sequence $`(x_n)`$ in $`K`$ such that $`(x_n\nu _n)`$ is relatively compact where $`\nu _n=\mu _n\mathrm{}\mu _2\mu _1`$ for all $`n1`$.

It is well-known that a locally compact group has normed convergence property if and only if it is compact (see Theorem 2.3.2 of ) and a Moore hypergroup has normed convergence property if and only if it is compact (see ). In this short note we extend the afore-mentioned result to hypergroups admitting an invariant Haar measure.

Let $`K`$ be a hypergroup admitting a right-invariant Haar measure $`m`$. Let $`L^p(K)`$ be the space of p-integrable functions on $`K`$ for $`1p<\mathrm{}`$. Let $`^1(K)`$ be the space of all non-negative measures on $`K`$ with total measure less than are equal to $`1`$. Now, for any $`\mu ^1(K)`$ and for $`1p\mathrm{}`$, define $`P_\mu :L^p(K)L^p(K)`$ by

$$P_\mu (f)(x)=f(xy)𝑑\mu (y)=f(y)d(x\mu )(y)$$

where $`f(xy)=f(z)d(\delta _x\delta _y)(z)`$, for all $`xK`$. It is easy to see that $`P_\mu `$ is a contraction on $`L^p(K)`$ for any $`\mu ^1(K)`$ and any $`1p\mathrm{}`$. The map $`\mu P_\mu `$ is an injection, so we may view $`^1(K)`$ as a subspace of bounded linear operators on $`L^2(K)`$.

We now recall the following topologies on $`^1(K)`$:

1. weak topology on $`^1(K)`$ is the weakest topology for which the function $`\mu \mu (f)`$ is continuous for all bounded continuous functions $`f`$ on $`K`$;
2. vague topology on $`^1(K)`$ is the weakest topology for which the function $`\mu \mu (f)`$ is continuous for all continuous functions $`f`$ with compact support;
3. weak operator topology (respectively, strong operator topology) on $`^1(K)`$ is the subspace topology inherited from the space of bounded linear operators on $`L^2(K)`$ with weak operator topology (respectively, strong operator topology).

For a locally compact group $`G`$, it is well-known that on $`^1(G)`$ the vague and weak operator topologies coincide and on $`𝒫(G)`$, the weak, vague and weak operator and strong operator topologies coincide (see Lemma 6.1.23 of ). The following extension of this result may be proved arguing as in Lemma 6.1.23 of  and we omit the details.

###### Proposition 1

Let $`K`$ be a hypergroup admitting a countable basis and a right-invariant measure. Then we have the following:

1. on $`^1(K)`$, the vague and weak operator topologies coincide;
2. on $`𝒫(K)`$, the weak, vague and weak operator and strong operator topologies coincide.

In particular, $`^1(K)`$ is compact in the weak operator topology.

We next prove the following:

###### Lemma 1

Let $`K`$ be as in Proposition 1. Suppose $`(\mu _n)`$ is a sequence in $`𝒫(K)`$ such that $`(x_n\mu _n)`$ is relatively compact. Then there exists a $`fL^2(K)`$ such that $`P_{\mu _n}(f)_2\to ̸0`$.

###### Proof

Since $`(x_n\mu _n)`$ is relatively compact, there exists a $`\delta >0`$ and a compact set $`C`$ such that $`x_n\mu _n(C)>\delta `$ for all $`n1`$.

Let $`f`$ be a continuous function on $`K`$ with compact support such that $`0f1`$ and $`f(x)=1`$ for all $`xC`$. Then

$$P_{\mu _n}(f)(x_n)=f(y)d(x_n\mu _n)(y)>\delta $$

for $`n1`$. Let $`(h_j)`$ be a compactly supported left approximate identity in the space of continuous functions vanishing at infinity. Then as in Proposition 3.1 of , we have

$$P_{\mu _n}(f)_{\mathrm{}}h_jff_{\mathrm{}}+P_{\mu _n}(f)_2h_j_2$$

where $`h_jf(y)=h_j(x)f(xy)𝑑m(x)`$. Since $`P_{\mu _n}(f)_{\mathrm{}}>\delta >0`$ for all $`n1`$, we get that $`P_{\mu _n}(f)_2\to ̸0`$.

###### Theorem 1

Let $`K`$ be a hypergroup admitting a right-invariant Haar measure and a countable basis. Then $`K`$ has normed convergence property if and only if $`K`$ is compact.

###### Proof

It is well-known that compact hypergroups have normed convergence property. We now prove the ’only if’ part. Suppose $`K`$ has normed convergence property. Let $`\mu 𝒫(K)`$ be a symmetric measure such that the support of $`\mu `$ is $`K`$. Then by assumption, there exists a sequence $`(x_n)`$ such that $`(x_n\mu ^n)`$ is relatively compact. By Lemma 1, we get that

$$P_\mu ^n(f)_2\to ̸0$$

as $`n\mathrm{}`$ for some continuous function $`f`$ with compact support. By convergence of alternating sequences (see ), in the strong operator topology

$$P_\mu ^{2n}P$$

where $`P`$ is a bounded linear operator on $`L^2(K)`$. This implies that $`P`$ is a self-adjoint projection. By Proposition 1, there exists a $`\rho ^1(K)`$ such that $`P=P_\rho `$. Since $`P_\mu ^n(f)\to ̸0`$, $`\rho 0`$ and $`\rho ^2=\rho `$. This means that $`\rho `$ is the normalized Haar measure on a compact subhypergroup $`H`$ of $`K`$ (see ). Since $`P_\mu ^2P=P`$, we get that the support of $`\mu ^2`$ is contained in $`H`$ (see ). Since the support of $`\mu `$ is $`K`$, we get that $`K=H`$ and hence $`K`$ is compact.

C. Robinson Edward Raja,
Stat-Math Unit,
Indian Statistical Institute,
8th Mile Mysore Road,
R. V. College Post,
Bangalore - 560 059.
India. e-mail: creraja@isibang.ac.in
or
e-mail: raja$`\mathrm{\_}`$robinson@hotmail.com