text:
# Parameter study of the diamagnetic relativistic pulse accelerator (DRPA) in slab geometry I: Dependence on initial frequency ratio and slab width

## I Introduction

In astrophysical, solar and space plasmas, the acceleration of high energy particles is an intriguing and unsolved problem. Magnetic reconnection, wave turbulence and collisionless shocks are often invoked for the acceleration of plasma particles, from cosmic rays, solar flares to gamma-ray burstspira00  and astrophysical jets.

A new mechanism of particle energization by freely expanding, hot, strongly magnetized collisionless plasmas was recently discovered by Liang et al. using two-and-a-half dimensional (2-1/2D) particle-in-cell (PIC) simulations lian03 . This mechanism, called the diamagnetic relativistic pulse accelerator (DRPA), is able to convert most of the initial magnetic energy into the ultra-relativistic directed kinetic energy of a small fraction of the surface particles. Moreover, Liang and Nishimuralian032  have found that the late-time plasma pulse of the DRPA reproduces many unique features of cosmic gamma-ray bursts. Hence the DRPA may shed important new light on the unsolved problem of gamma-ray bursts. When a hot, strongly magnetized ($`\beta `$=thermal pressure/plasma pressure $`1`$) collisionless plasma expands into a vacuum or low-density surrounding, an electromagnetic pulse is formed at the expansion surface, which traps and accelerates a small fraction of surface particles via the pondermotive forcelian03 ; mill59 , and the magnetic energy is efficiently converted into directed particle energy .

In our previous Letterslian03 ; nish03 , the DRPA was studied by performing 2-1/2D PIC simulations for both electron-positron (mass ratio is $`m_i/m_e=1`$) and electron-ion ($`m_i/m_e=100`$) plasmas with a specific set of initial parameters, and only carried out to 1000$`\mathrm{\Omega }_e^1`$ (where $`\mathrm{\Omega }_e=eB_0/m_e`$ is the electron cyclotron frequency defined by the initial magnetic field). In this paper, we will study the parameter dependence and the long-term evolution of the DRPA of electron-positron plasmas by means of a series of large-scale PIC simulations.

The outline of this paper is as follows. Section II describes the simulation model and the simulation parameters. The benchmark results of Liang et al.lian03  are reviewed and summarized in Section III.1 for comparison with other cases. The long-term evolution of the DRPA is described in Sec. III.2. The initial parameter dependence of DRPA is discussed in Secs. III.3 and III.4. Section IV gives a summary of our results.

## II Simulation model

In our particle simulation codenish02 , we use a 2-1/2D explicit simulation scheme based on the PIC method for time advancing of plasma particles. In this method, spatial grids are introduced to calculate the field quantities, and the grid separations are uniform, $`\mathrm{\Delta }x=\mathrm{\Delta }z=\lambda _e`$, where $`\lambda _e`$ is the electron Debye length defined by $`c/\omega _{pe}`$ ($`c`$ is the speed of light and the electron plasma frequency is $`\omega _{pe}=\sqrt{e^2n_0/ϵ_0m_e}`$; $`e`$ is elementary charge, $`n_0`$ is the initial electron density, and $`ϵ_0`$ is the dielectric constant of vacuum). In our simulations, the thermal speed of electrons is almost equal to $`c`$ because the plasma temperature is $`k_BT_e=k_BT_p=5MeV`$, where the subscripts $`e`$ and $`p`$ refer to electrons and positrons. The simulation domain on the $`xz`$ plane is $`L_x/2xL_x/2`$ and $`L_z/2zL_z/2`$.

We treat plasma particles as superparticles, with 20 superparticles per cell (if the superparticles are spread uniformly on the entire simulation domain) for each component (electrons and positrons). Most of the particle densities, however, concentrate on the expanding surface in our study and we are mainly interested in particle behaviors near the surface. Thus, this number of superparticles is enough to accurately resolve the physics of surface particle acceleration. To prevent large violations of Gauss’s law due to numerical noises, the Marder’s method for the electric field correctionmard87  is adopted at every time step in the code.

We use a doubly periodic system in $`x`$ and $`z`$ directions, and the system length is $`L_x=3840\mathrm{\Delta }x`$ (Run A), $`3840\mathrm{\Delta }x`$ (Run B), $`1920\mathrm{\Delta }x`$ (Run C), and $`7680\mathrm{\Delta }x`$ (run D) and $`L_z=10\mathrm{\Delta }z`$ (for all four Runs). Initially, the electron and positron distributions are assumed to be relativistic Maxwellian with spatially uniform temperature, $`k_BT_e=k_BT_p=5MeV`$. The spatial distribution of initial plasmas is assumed to be a uniform slab with the lengths $`12\mathrm{\Delta }x\times 10\mathrm{\Delta }z`$,(runs A,B,C) and $`60\mathrm{\Delta }x\times 10\mathrm{\Delta }z`$ (run D) and the plasmas are located in the center of the system. The background magnetic field $`𝑩_0=(0,B_0,0)`$ initially exists only inside the plasma.

We will first compare three simulation results by changing the frequency ratio $`\omega _{pe}/\mathrm{\Omega }_e`$. The values are $`\omega _{pe}/\mathrm{\Omega }_e=0.104`$ (Run A), 1.0 (Run B), and 0.01 (Run C). At the end we will also discuss the result of changing the initial slab thickness (Run D).

## III Results

### III.1 Diamagnetic relativistic pulse accelerator

First, we illustrate the mechanism of the diamagnetic relativistic pulse accelerator (DRPA) by reviewing the early ($`t\mathrm{\Omega }_e10^3`$) results of the benchmark Run A (Liang et al., 2003lian03 ). In Run A, the magnetic field energy is initially almost equal to the thermal energy of plasma particles, i.e., electron and positron beta is $`\beta _e=\beta _p=0.551`$ ($`\beta =2\mu _0n_0k_BT/B_0^2`$, where $`\mu _0`$ is the magnetic permeability of vacuum).

In Figure 1(a), the system-integrated energies of the magnetic field (line 1), electric field (line 2), electrons (line 3), and positrons (line 4) are shown as functions of time. The initial static magnetic field rapidly transforms into two oppositely propagating electromagnetic pulses, which propagate towards the vacuum region. The plasma simultaneously begins to expand.

We see the rapid decrease of the magnetic energy and the increase of the electric field energy due to the generation of the electromagnetic (EM) pulse during the period shown as Phase 1 in Fig. 1(a). In Phase 2, both the magnetic and electric field energies begin to decrease and they are converted into the electron and positron kinetic energies. Also, the EM pulse with the amplitude $`cB_yE_z`$ is formed at the expansion surface during this period. In Phase 3, the energy conversion from the field energies to the particle kinetic energies still continues although the conversion rate becomes slower compared with that in Phase 2.

Figure 2 illustrates the spatial profiles of the (a) magnetic field $`B_y`$, (b) electric field $`E_z`$, (c) current density $`J_z`$, and (d) electron density $`n_e`$ at $`t\mathrm{\Omega }_e=0`$, 80, 500, and 1000 for Run A. The magnetic field and the electron density for $`x0`$ are identical, on the other hand, the electric field and the current density become antisymmetric about $`x=0`$. The values of $`E_z`$ and $`J_z`$ are zero at $`t\mathrm{\Omega }_e=0`$. Only the profiles near the expansion surface are drawn at $`t\mathrm{\Omega }_e=500`$ and 1000 to avoid making the figure cluttery. At $`t\mathrm{\Omega }_e=80`$, the EM pulse is formed and begin to expand with plasmas into the vacuum. The gradient drift in the z-direction induces a current density. At $`t\mathrm{\Omega }_e=500`$, the EM pulse has been formed at the expansion surface. Most of the plasma particles follow the EM pulse propagation and the plasma density has a peak corresponding to particles trapped by the pondermotive force of the pulse. The amplitude of the EM pulse is gradually damped as EM energy is transformed into particle energy.

The $`xz`$ configuration of electrons is plotted at three different times in Figure 3. The electron-positron slab with the internal magnetic field expands rapidly into the vacuum region. The behavior of positrons is identical with that of electrons. The expansion speed of the plasma surface in the $`x`$ direction is almost the speed of light, $`V_xc`$. Several simulation studies of non-relativistic plasma expansion across an initially uniform magnetic field have been carried outsydo83 ; wins88 ; sgro89 ; gisl89 , and the lower hybrid drift instability has been found to deform the structure of the expanding plasma surface. As shown in Fig. 3, there appears no instability at the plasma surface in our result. There are two main differences between our simulation and the previous studiessydo83 ; wins88 ; sgro89 ; gisl89 : (a) in their studies, an initial magnetic field exists uniformly all over the simulation plane and the plasma expands (or diffuses) across the field, (b) the electron-ion plasma was considered in their studies. Thus, the charge separation due to the large mass ratio causes the different drift speed between electrons and ions, which makes the lower hybrid drift wave unstable. These differences may explain the absence of instability in our simulation. In addition, the effective Larmor radius of relativistic plasmas given by $`\gamma (c/\mathrm{\Omega }_e)`$ is larger than our system length in the $`z`$ direction at late times. Thus, any effect of the Larmor radius play no important role in our simulations, and we cannot see any effect caused by the Larmor motion of relativistic plasmas.

Figure 4 shows phase plots of electrons at four different times in the space (a) $`xP_x`$ and (b) $`xP_z`$ for Run A. It is found that most of the electron population moves to the expansion surface and the momentum $`P_x`$ of the surface electrons increases as time elapses. The averaged momentum reaches $`<P_x>100m_ec`$ and the maximum momentum of energetic particles reaches $`P_{x,max}200m_ec`$ at $`t\mathrm{\Omega }_e=1000`$. On the other hand, the momentum $`P_z`$ indicates no dramatic increase even at late times. This means that the plasma particles are effectively accelerated into the $`x`$ direction rather than the $`z`$ direction.

### III.2 Long-term evolution of the DRPA

Let us focus on the long-term evolution of the DRPA in Run A. Figure 5 shows the spatial profiles of the (a) magnetic field $`B_y`$, (b) electric field $`E_z`$, (c) current density $`J_z`$, (d) electron density $`n_e`$, and (e) the phase plots $`xP_x`$ of electrons for Run A at $`t\mathrm{\Omega }_e=6000`$, 10000, 14000, and 18000.

Only the portion near the expansion front at each time has been displayed in each panel. In Figs. 5(a) and 5(b), the wave packet is dispersively spreading and the EM pulse amplitude gradually decreases as time elapses. The reduction of the pulse amplitude corresponds to the conversion of the EM energy into the kinetic energy of the expanding plasma. In Fig. 5(e), we recognize the electron momentum $`P_x`$ near the first pulse of the wave packet dramatically increases as time goes on. The maximum value of $`P_x`$ reaches more than $`500m_ec`$ at $`t\mathrm{\Omega }_e=18000`$.

Another significant feature of the late time behavior in Run A is the bifurcation of plasma density. In the initial phase of the DRPA, the electron density near the wave packet has only one peak as shown in Fig. 2(d). At late times in Fig. 5(d), the peak has been split into several peaks. The bifurcation of the plasma density profile results from the development of new transverse current and reverse currents behind the front, which forms new traps for particles in the pondermotive potential.

Figure 6(a) shows the momentum distributions of electrons which are located in the vicinity of the EM pulse at $`t\mathrm{\Omega }_e=1000`$, 6000, and 10000. The momentum peak at which the distribution takes a maximum value tends to increase as time elapses, i.e., $`P_x45m_ec`$ at $`t\mathrm{\Omega }_e=1000`$, $`P_x60m_ec`$ at $`t\mathrm{\Omega }_e=6000`$, and $`P_x80m_ec`$ at $`t\mathrm{\Omega }_e=10000`$ in Fig. 6(a). Also in Figure 7(a), the temporal evolution of this momentum is shown until $`t\mathrm{\Omega }_e=20000`$ for Run A. This temporal increment of the momentum verifies that the average electon near the pulse is gradually accelerated. In addition, we can see in Fig. 6(a) that the tail part of the distribution develops into a quasi-power-law. At $`t\mathrm{\Omega }_e=10000`$, a second peak of the distribution appears at a high value $`P_x`$, which is consisted of non-thermal electrons. In short, some of the surface electrons is more accelerated than other surface electrons. This multi-peak structure of the momentum distribution may be related to the bifurcation of the electron density shown in Fig. 5(d).

### III.3 Frequency ratio dependence of the DRPA

In the previous sections, we summarize the results of Run A (with the frequency ratio $`\omega _{pe}/\mathrm{\Omega }_e=0.105`$). Next, we examine the frequency ratio dependence of the DRPA. The frequency ratios discussed here are $`\omega _{pe}/\mathrm{\Omega }_e=1.0`$ (Run B) and $`\omega _{pe}/\mathrm{\Omega }_e=0.01`$ (Run C). Other fundamental simulation parameters are exactly same as those of Run A.

Figure 8 displays the spatial profiles of the same quantities shown in Fig. 2, in the case of Run B, at different times. Figs. 8(a), (b), (c), and (e) mainly show the region near the EM pulse at each time, and the values of $`E_z`$ and $`J_z`$ are zero at $`t\mathrm{\Omega }_e=0`$. The electron density $`n_e`$ at $`t\mathrm{\Omega }_e=0`$ is not shown in Fig. 8(d) since the initial density $`n_e/n_0=1`$ is much larger than the values at the other three times. Compared with the results of Run A, no apparent EM pulse with large amplitudes can be formed at the expansion surface in Figs. 8(a) and 8(b). In Fig. 8(d), we find that the electrons are distributed in a broad region, rather than confined in the EM pulse. We see no sign of bifurcation of the density pulse like in Run A. In Run B, initially the thermal energy of plasma particles is much larger than the magnetic field energy because of $`\beta _e=\beta _p=50`$. Temporal evolution of each energy component is shown in Fig. 1(b). In Run A, the field energy is converted into the kinetic energy of surface particles by the DRPA, thus the total field energy decreases as time elapses. In Run B, however, the field energy continues to increase as shown in Fig. 1(b). This increase of the EM field energy is mainly caused by the field fluctuation in the regions away from the expansion surface.

We focus on the amplitude of the propagating magnetic field formed at the expansion surface. Temporal evolution of the peak value of the magnetic field pulse $`B_y`$ is shown in Figure 9(a). As for Run A, the peak value dramatically decreases during $`0<t\mathrm{\Omega }_e<80`$ (this period corresponds to Phase 1 in Fig. 1) and then continues to decrease at a constant damping rate (this period corresponds to Phases 2 and 3 in Fig. 1). We can see the similar tendency in the temporal evolution of the peak value for Run B during $`0<t\mathrm{\Omega }_e<500`$, but the damping rate of the peak value is much larger than that of Run A. At $`t\mathrm{\Omega }_e500`$, the peak value becomes $`B_y<0.3B_0`$. This means that the magnetic field energy in Run B is dissipated much faster than that in Run A. Although the magnetic field energy is rapidly converted into particle energy, the total kinetic energy of Run B hardly changes since the initial plasma beta is very high. Therefore, the evolution can be described as the free-expansion of the plasma rather than the expansion accompanied by the acceleration of the EM pulse. The electron momentum peak of Run B, at which the momentum distribution takes a maximum, is shown as a function of time in Fig. 7(b). The momentum remains around $`P_x15m_ec`$ even when the amplitude of the magnetic field pulse becomes $`B_y0.3B_0`$ at $`t\mathrm{\Omega }_e=500`$ because the converted magnetic field energy is too small to accelerate all the plasma particles via the DRPA. Thus, the electrons at late times exhibit almost uniform distribution by the free-expansion and the density pulse is not formed in Run B, as shown in Fig. 8(d). We see the evolution of the electron momentum distribution for Run B in Fig. 6(b). The distribution keeps a Maxwellian form even at late times. This means that the plasma is consisted almost only of the thermal component.

Figure 10 shows the results of Run C ($`\omega _{pe}/\mathrm{\Omega }_e=0.01`$ and $`\beta _e=\beta _p=0.005`$). In this run the initial magnetic field energy is much greater than the thermal energy of plasma particles. The converted energy of the EM pulse can be effectively used for the particle acceleration \[see Fig. 1(c)\]. In Fig. 10(d), most of the electrons has been accelerated owing to the EM pulse and the electron density is much more concentrated at the expanding surface. Temporal evolution of the maximum amplitude of the pulse magnetic field in Run C is also plotted in Fig. 9(a). The period during $`0<t\mathrm{\Omega }_e<750`$ represents an initial relaxation phase in Run C, which corresponds to Phase 1 in Fig. 1(a) of Run A. After that, the peak value stays almost constant, $`B_y0.62B_0`$. For comparison, the result of a vacuum EM pulse without plasma is shown in Fig. 9(a). In this run without plasma, we set up the initial magnetic field $`B_y`$ to be $`B_y(x)=B_0`$ for $`|x|<3\mathrm{\Delta }x`$ and $`B_y(x)=0`$ for $`|x|3\mathrm{\Delta }x`$, and the electric field is initially zero everywhere. The EM pulse generated by the spatial gradient of the initial magnetic field propagates according to the wave equation in a vacuum. It is found that the amplitude of the EM wave is almost unchanged after the initial relaxation phase ($`0<t\mathrm{\Omega }_e<80`$) in Fig. 9(a). The EM pulse in Run C is very similar to the pulse without plasma since the initial plasma beta is small, $`\beta =0.005`$, and only a small portion of the total magnetic field energy is converted into the plasma kinetic energy.

The evolution of the electron momentum distribution for Run C is shown in Fig. 6(c). Since most of electrons are accelerated by the DRPA, the tail part of the distribution exhibits a power-law with sharp slope. In Fig. 7(c), the peak momentum of the surface electrons increases rapidly as time elapses, compared to Runs A and B with higher frequency ratios. The average momentum of electrons reaches $`<P_x>1500m_ec`$ at $`t\mathrm{\Omega }_e=30000`$.

### III.4 Dependence on the initial plasma width

In this section, we examine how the width of the initial plasma can affect the DRPA. The widths of the density pulse of expanding plasmas tend to scale with the width of the initial plasma, as we expect from causality. The propagating EM pulses have the same tendency since the initial magnetic field has the same width as the plasma density. We have carried out Run D, in which the frequency ratio is $`\omega _{pe}/\mathrm{\Omega }_e=0.105`$ and the initial width of plasma is five times that of Run A.

The temporal evolution of each energy is shown in Fig. 1(d). The apparent energy conversion from field energies into particle energies due to the DRPA can be seen until $`t\mathrm{\Omega }_e5000`$, which is similar to the evolution in Fig. 1(a). After that time, the total EM energy slightly increases and the total particle energies tends to decrease. These slight energy changes against the DRPA is caused by the local energy changes in the regions away from the expansion surface. In those regions, the particle acceleration due to the DRPA hardly works efficiently. Even though the total particle energies slightly decreases, we will show that the plasmas near the propagating EM pulse continue to be energized after $`t\mathrm{\Omega }_e5000`$.

The electron momentum distributions for Run D are displayed in Fig. 6(d) at three different times. As time elapses, the average momentum of surface electrons increases, and we recognize a lot of peaks in the high momentum part at $`t\mathrm{\Omega }_e=10000`$. As we will show later in Figure 11, these peaks correspond to a repeated bifurcation of electron density at the expansion surface. In Fig. 7(a), the results of Run D are also plotted. There is no big difference between the peak electron momenta of Runs A and D, except in D the acceleration seems more delayed due to the prolific early bifurcations. However the maximum energy achieved by the most energetic particles is higher in D than in A. This is likely caused by the larger coherence length for acceleration in D due to the wider EM pulse.

Fig. 11 shows the spatial profiles of the (a) magnetic field $`B_y`$ and (b) electron density $`n_e`$, and (c) the phase plots $`xP_x`$ of electrons for Run D at $`t\mathrm{\Omega }_e=6000`$, 10000, 14000, and 18000. Only the portion near the expansion surface at each time has been displayed in each panel. Compared with the electron density of Run A in Fig. 5(b), the width of the density of the expanding electrons $`w`$ becomes larger, e.g., $`w2000c/\mathrm{\Omega }_e`$ at $`t\mathrm{\Omega }_e=18000`$ in Fig. 11(b). We also recognize many more bifurcations of the electron density at late times.

## IV Summary

We have used two-and-a-half-dimensional large-scale particle-in-cell plasma simulations to investigate the particle energization of expanding relativistic electron-positron plasmas. When magnetized plasmas with high temperature expand into a vacuum, the particles at the expansion surface are accelerated in the expansion direction. This energization mechanism is called DRPA (diamagnetic relativistic pulse accelerator)lian03 .

In this paper, we mainly discuss the long-term evolution and the initial parameter dependences of DRPA. In the case of the electron-positron plasma with the frequency ratio $`\omega _{pe}/\mathrm{\Omega }_e=0.105`$ and the temperature $`k_BT_e=k_BT_p=5MeV`$, an electromagnetic (EM) pulse with large amplitude is generated and propagates into a vacuum. The plasma particles are trapped and accelerated by the pondermotive force. As time elapses, since the magnetic field energy is converted into the kinetic energy of the particles, the amplitude of the EM pulse is decreased and the particles near the expansion surface are energized. The expanding plasma density is concentrated on the vicinity of the expansion surface at the early phase of the expansion. And at the late phase, we observe the bifurcation of plasma density due to the formation of new transverse currents.

The properties of DRPA strongly depend on the frequency ratio $`\omega _{pe}/\mathrm{\Omega }_e`$. In the case of a high frequency ratio ($`\omega _{pe}/\mathrm{\Omega }_e=1.0`$), the EM pulse is rapidly damped and the plasma is expanding freely with little acceleration since the energy of the plasma is much larger than the magnetic energy. On the other hand, in the case of a low frequency ratio ($`\omega _{pe}/\mathrm{\Omega }_e=0.01`$), the magnetic energy is much higher than the plasma thermal energy. Most of the expanding plasma is trapped by the EM pulse and the plasma is efficiently energized at late times, achieving higher Lorentz factors than in the other cases.

In addition to the frequency ratio, we have studied the effect of varying the initial width of the plasma slab. As expected, the width of the asymptotic EM pulse tends to be proportional to the initial width of the plasma slab. Both the EM pulse and the plasma density are broadened, and the density pulse develops many more peaks at late times due to prolific bifurcations. We also find that the maximum energy of the accelerated particles is higher in this case due to the longer pulse length.

Finally, we note the system length in the $`z`$ direction in our two-dimensional simulation domain. As we mentioned in Sec. III.1, the effective Larmor radius of relativistic electron-positron plasmas $`\gamma (c/\mathrm{\Omega }_e)`$ becomes comparable or larger than $`L_z`$ at late times. Thus, the present system length $`L_z`$ will be insufficient to fully cover the physical phenomena scaled by the Larmor radius if there exist some instabilities which can develop along $`L_z`$ on that scale. However, to our knowledge, previously observed instabilities in non-relativistic simulationssydo83 ; wins88 ; sgro89 , such as the lower hybrid drift instability, do not grow in our electron-positron relativistic plasmas. Future advanced simulations will be able to give the clear answer to this problem through further simulation studies with larger system lengths.

###### Acknowledgements.

The work of K.N. was performed under the auspices of the U. S. DOE and was supported by the DOE Office of Basic Energy Sciences, Division of Engineering and Geosciences, the LDRD Program at Los Alamos, and the Sun-Earth Connections Theory Program of NASA. E.L. was partially supported by NASA grant NAG5-7980 and LLNL contract B510243.


format: text

text:
# Roto-vibrational spectrum and Wigner crystallization in two-electron parabolic quantum dots

## I Introduction

For a low enough electron density, WignerW1  predicted that electrons should localize creating an ordered spatial structure, thenceforth named a Wigner crystal, that breaks the complete translational symmetry of the homogeneous electron gas (also see Ref. Mah, ). Indeed, the formation of the Wigner crystal was observed in two-dimensional (2D) distributions of electrons on the surface of liquid helium.2  A phase transition, induced by the magnetic field, from an electron liquid to a crystalline structure has also been reported for a 2D electron plasma at a GaAs/AlGaAs heterojunction.3

The existence of different phases in quantum dots, where a few electrons are confined into a small space, has become a topical subject in mesoscopic physics (see, for a recent review, Ref. RM, ). In fact, the high controllability of quantum dots suggests that these systems could provide an attractive opportunity to achieve the appropriate conditions for localized states. It is precisely to stress this controllability that the names artificial atoms and quantum dots have been coined.

There is a general persuasion that the Wigner crystallization in quantum dots, whose localized states are referred to as Wigner molecules, should occur at significantly larger densities than in the 2D bulk. It is based on the argument that in quantum dots potential-energy contributions can easily exceed the kinetic terms and, therefore, electronic motion can be effectively quenched by manipulating the external confinement and/or an applied magnetic field. As for the homogeneous gas, one would expect that in crystallized states the kinetic energy is solely that of the vibrational zero-point motion of the electrons around their mean positions, much smaller than the interaction (potential) energy. Various approaches including ab initio calculations within diffusion and path integral Monte Carlo methods, Hartree-Fock and spin-density functional methods etc have been applied to analyze the onset of the crystallization.RM  However, a nonambiguous theoretical result that would justify the above conjecture for a zero magnetic field is lacking. The case with an intense magnetic field is better understood since the magnetic field induces an edge reconstruction, beginning with the appearance of localized vortices on the outer region, that ultimately propagates to all the dot for very high $`B`$’s.RMM ; RM

In the simpler case of a two-electron 2D quantum dot at zero magnetic field, Yannouleas and LandmanYL  pointed out that the excited-state energies of this system closely follow the rotor sequence when the repulsion-to-confinement ratio, as given by the Wigner parameter $`R_W`$, is large enough ($`200`$). This was shown to be a proof of the crystallization of the two electrons on fixed positions in a reference frame which is rotating. Quite remarkably, the hypothesized rotating Wigner molecule fulfills at the same time the strict symmetry conditions of quantum mechanics –circularity in this case– and the obvious preference for opposite positions when repulsion is large enough. This is a major difference from the above mentioned bulk case where a Hamiltonian symmetry (translation) is broken by the crystallized state. For Wigner molecules, symmetries are preserved in the laboratory frame and one must consider an intrinsic (rotating) frame to ’see’ the underlying deformation. A similar situation is found for particular states of two-electron atoms that have been much investigated in physical chemistry (we address the reader to the review paper by BerryBer ). For the two-electron quantum dot, however, the crystallization condition from Ref. YL, , $`R_W200`$, looks disappointing since it seems unrealistic to achieve such a value experimentally.

Although the exact ground-state wave function of the two-electron artificial atom can be obtained, at least numerically, it may seem paradoxical that one also needs the excited states in order to ascertain the existence of a crystallization. In fact, this inability to disentangle the system’s intrinsic structure from its full wave function in a clear way can be taken as a weakness of the ab initio, symmetry preserving, approaches. In general, even in cases when the exact ground- and excited-state wave functions and energies are known, an intrinsic deformation can only be inferred by comparing with the result of simpler models in which either symmetries are relaxed or the intrinsic structure is imposed. A clear example of the former approach is given by the unrestricted Hartree-Fock (HF) method for the ground stateKoo ; Yan  followed by the random-phase approximation (RPA) for excitations.rpa  Conversely, the roto-vibrational model of Wendler et al.Wen  for two electrons in a ring could be included in the latter category.

One should be aware that when symmetries are relaxed, as in the Hartree-Fock approach, artifacts or unphysical properties may appear. In a recent contribution Reusch and GrabertRG  discussed the validity of the latter, drawing special attention to the caution with which one must take Hartree-Fock predictions on symmetry breaking, in agreement with the results presented below. Therefore, a complete physical understanding requires both exact results and model solutions. This way the system’s intrinsic deformations are physically understood while, at the same time, artifacts can be safely discarded. A paradigmatic case where the proposed analysis can be performed is given by the two-electron 2D parabolic quantum dot. The separation of center-of-mass and relative coordinates along with the circular symmetry restriction allows the reduction of the multidimensional Schrödinger equation to just a radial one, easily solvable numerically. On the other hand, the Hartree-Fock and RPA solutions without any symmetry restriction can also be obtained. A most convenient basis for this latter calculations is given by the Fock-Darwin orbitals in terms of which one can analytically develop much of the required algebra. It is our aim in this work to determine the crystallization onset of two-electron parabolic dots by recourse to the three different approaches referred to above; namely, (a) an analytical roto-vibrational model, (b) a numerical solution of the Schrödinger equation, and (c) symmetry unrestricted Hartree-Fock and random-phase approximations.

Hereafter, we refer to the solution of the Schrödinger equation for the two-electron parabolic dot as the exact solution. It should be pointed out that, as shown by Taut,Taut  this Schrödinger equation is analytically solvable only for particular confinement/interaction strengths. For general values of this quantity a numerical treatment is required. As mentioned above, the most straightforward one is an integration of the radial equationYL ; Zhu  but, nevertheless, other methods such as diagonalization in a basisPfan ; Wag  and the Monte Carlo method Ped ; Loz  have also been applied. One of us has used the so-called oscillator representation method, perturbatively treating the residual interaction, to derive analytical expressions for the energy levels.Din

The paper is organized as follows. Section II introduces the magneto-parabolic units that allow one to trace the evolution of ground and excited states of artificial atoms at various conditions. An analytical roto-vibrational model for the two-electron parabolic quantum dot is described in Sec. III. Section IV provides details of our numerical calculation of exact solutions and compares these solutions with those of the roto-vibrational model. Section V analyzes the reliability of the Hartree-Fock and RPA results for the present system. A short summary is finally drawn in Sec. VI.

## II Magneto-parabolic unit system

We consider two electrons with a Coulomb interaction between them. The electrons move in the $`xy`$ plane where a circular parabolic confinement induces the formation of an electron island. The system is also subject to an external magnetic field applied in the vertical direction ($`z`$). The full Hamiltonian thus reads

$``$ $`=`$ $`{\displaystyle \underset{i=1,2}{}}\left[{\displaystyle \frac{1}{2m}}\left(𝐩+{\displaystyle \frac{e}{c}}𝐀\right)^2+{\displaystyle \frac{1}{2}}m\omega _0^2r^2\right]_i+{\displaystyle \frac{e^2}{\kappa r_{12}}}`$ (1)
$`+`$ $`g^{}\mu _BBS_z.`$
In Eq. (1), $`m`$, $`\kappa `$, and $`g^{}`$ are the electron’s effective mass, the dielectric constant, and the effective gyromagnetic factor, respectively, and we have used planar polar coordinates ($`r^2=x^2+y^2`$). The two contributions within the square brackets are, respectively, the generalized kinetic energy in terms of the vector potential $`𝐀`$, and the external confinement. Within the so-called symmetric gauge one has $`𝐀(x,y)=B/2(y,x)`$. The next contribution is the Coulomb repulsion and, finally, the last term is the Zeeman energy involving the total spin operator $`S_z`$ and the universal Bohr’s magneton $`\mu _B=e\mathrm{}/2m_ec`$.

It is well known that in the chosen gauge the magnetic field contributions can be recast into the form of an effective parabolic confinement of frequency $`\mathrm{\Omega }=\sqrt{\omega _0^2+\omega _c^2/4}`$, where $`\omega _c=eB/mc`$ is the cyclotron frequency, and an additional angular-momentum-dependent term $`(\mathrm{}\omega _c/2)\mathrm{}_z`$ (cf. Refs. Hawrylak, ; Chak, ). The magneto-parabolic units (mpu) we shall use consist of taking $`\mathrm{}\mathrm{\Omega }`$ as the energy unit and $`\mathrm{}_\mathrm{\Omega }\sqrt{\mathrm{}/(m\mathrm{\Omega })}`$ as the length unit. In addition, one also imposes $`\mathrm{}`$ as angular momentum unit which, in turn, fixes the time unit $`\tau _\mathrm{\Omega }=1/\mathrm{\Omega }`$. Summarizing in the standard abuse of notation we may write $`\mathrm{}=\mathrm{\Omega }=\mathrm{}_\mathrm{\Omega }=1\mathrm{mpu}`$. This is a natural choice for magneto-parabolic confinements and it allows one to express the spatial part of the Hamiltonian in terms of only two adimensional parameters, namely,

$`R_{\mathrm{𝑚𝑝}}`$ $`=`$ $`{\displaystyle \frac{e^2/(\kappa \mathrm{}_\mathrm{\Omega })}{\mathrm{}\mathrm{\Omega }}},`$ (2)
$`W_{\mathrm{𝑚𝑝}}`$ $`=`$ $`{\displaystyle \frac{\omega _c}{\mathrm{\Omega }}}.`$ (3)
Note that $`R_{\mathrm{𝑚𝑝}}`$ and $`W_{\mathrm{𝑚𝑝}}`$ give, respectively, the ratios of Coulomb interaction strength and cyclotron frequency to effective confinement. In the absence of a magnetic field $`R_{\mathrm{𝑚𝑝}}`$ coincides with the so-called Wigner parameter $`R_W`$ of Ref. YL, . Also note that $`W_{\mathrm{𝑚𝑝}}`$ has a maximal value $`W_{\mathrm{𝑚𝑝}}=2`$ that corresponds to a zero confinement $`\omega _0=0`$. We also mention that Reusch and Grabert used these parameters in their recent Hartree-Fock calculations.RG

Within the mpu system the Hamiltonian reads

$``$ $`=`$ $`{\displaystyle \underset{i=1,2}{}}\left[{\displaystyle \frac{1}{2}}^2+{\displaystyle \frac{1}{2}}r^2+{\displaystyle \frac{W_{\mathrm{𝑚𝑝}}}{2}}\mathrm{}_z\right]_i+{\displaystyle \frac{R_{\mathrm{𝑚𝑝}}}{r_{12}}}`$ (4)
$`+`$ $`{\displaystyle \frac{g^{}m^{}}{2}}W_{\mathrm{𝑚𝑝}}S_z,`$
where $`m^{}`$ is the adimensional ratio of effective to bare electron mass, i.e., $`m^{}=m/m_e`$. The passage from the mpu system, with a given $`R_{\mathrm{𝑚𝑝}}`$ and $`W_{\mathrm{𝑚𝑝}}`$, to physical units requires the knowledge of the effective mass $`m`$ and the dielectric constant $`\kappa `$. More specifically, with fixed $`m`$ and $`\kappa `$ one can invert Eqs. (2) and (3) for the effective confinement $`\mathrm{\Omega }`$ and cyclotron $`\omega _c`$ frequencies or, equivalently, for the external confinement $`\omega _0`$ and the magnetic field $`B`$. In terms of these results the physical values of the energy ($`\mathrm{}\mathrm{\Omega }`$) and length ($`\mathrm{}_\mathrm{\Omega }`$) units are readily found.

Figure 1 shows the equivalence between the adimensional parameters $`(R_{\mathrm{𝑚𝑝}},W_{\mathrm{𝑚𝑝}})`$ of the mpu Hamiltonian and the physical values of $`(\mathrm{}\omega _0,B)`$ for the case of a GaAs host semiconductor having $`m=0.067m_e`$ and $`\kappa =12.4`$. The advantage of working with the mpu system becomes obvious when realizing that, with unit redefinition, the same numerical results can be equally applied to a variety of confinements, magnetic fields, and material parameters (effective mass and dielectric constants). Therefore, the model acquires a certain degree of universality. We also expect that quantum dot properties such as scaling laws or phase diagrams will be better displayed in terms of the adimensional mpu parameters.

## III A roto-vibrational model

Taking Eq. 4 as a starting point and introducing the standard center of mass $`(R,\mathrm{\Theta })`$ and relative $`(r,\theta )`$ coordinates it is well known that the Hamiltonian separates and, therefore, that the wave function factorizes. The center of mass (CM) problem is that of a single particle in a harmonic potential and magnetic field, having an analytical solution in terms of Fock-Darwin orbitals and energies $`\epsilon _{NM}^{(\mathrm{CM})}=2N+|M|+1+MW_{\mathrm{𝑚𝑝}}/2`$, with $`(N,M)`$ the radial and angular momentum center-of-mass quantum numbers. Focusing next on the relative problem, one introduces the wave function $`e^{im\theta }u_{nm}(r)/\sqrt{r}`$ having good $`\mathrm{}_z`$ angular momentum ($`m`$), and an additional quantum number $`n`$ whose meaning will be clarified below. The equation for the unknown $`u_{nm}(r)`$ reads

$$u_{nm}^{\prime \prime }+\left[\stackrel{~}{\epsilon }_{nm}^{(\mathrm{rm})}\left(\frac{1}{4}r^2+\frac{R_{\mathrm{𝑚𝑝}}}{r}+\frac{m^21/4}{r^2}\right)\right]u_{nm}=0,$$
(5)
where we have defined $`\stackrel{~}{\epsilon }_{nm}^{(\mathrm{rm})}=\epsilon _{nm}^{(\mathrm{rm})}mW_{\mathrm{𝑚𝑝}}/2`$ in terms of the relative-motion energy $`\epsilon _{nm}^{(\mathrm{rm})}`$ and the mpu parameter $`W_{\mathrm{𝑚𝑝}}`$.

Equation (5) will be the basis of our roto-vibrational model. Note that it resembles a Schrödinger one-dimensional equation with an effective potential

$$V_{\mathrm{𝑒𝑓𝑓}}(r)=\frac{1}{4}r^2+\frac{R_{\mathrm{𝑚𝑝}}}{r}+\frac{m^21/4}{r^2}$$
(6)
that includes the rotational motion term $`m^2/r^2`$ characterized by the angular momentum quantum number $`m`$. We can expect a rigid-rotor behaviour if $`V_{\mathrm{𝑒𝑓𝑓}}(r)`$ has a deep minimum at a particular value $`r=r_0`$. When this occurs the situation resembles that of diatomic molecules like H2, where the potential well for nuclear motion is described by the Morse potential (see Ref. BJ, ). In the present case the effective potential indeed has a minimum although it is in general rather shallow. This property is responsible for the coupling between rotation and vibration or, equivalently, for the floppiness of the rotating molecule mentioned in Ref. YL, .

The minimum condition on $`V_{\mathrm{𝑒𝑓𝑓}}(r)`$ yields the rotor radius from

$$\frac{r_0}{2}\frac{R_{\mathrm{𝑚𝑝}}}{r_0^2}\frac{2(m^21/4)}{r_0^3}=0.$$
(7)
Neglecting the third contribution on the left-hand-side, an assumption that will be valid for large enough $`r_0`$, one finds the asymptotic law $`r_0(2R_{\mathrm{𝑚𝑝}})^{1/3}`$. Now, expanding to second order around $`r_0`$ we approximate

$`V_{\mathrm{𝑒𝑓𝑓}}(r)`$ $``$ $`V_{\mathrm{𝑒𝑓𝑓}}(r_0)+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{3}{2}}+2{\displaystyle \frac{m^21/4}{r_0^4}}\right)(rr_0)^2`$ (8)
$`=`$ $`\mathrm{𝑐𝑜𝑛𝑠𝑡}.+{\displaystyle \frac{1}{2}}k(rr_0)^2,`$
a result that, when substituted into Eq. (5) for the round parentheses, leads immediately to the analytical predictionana

$`\stackrel{~}{\epsilon }_{nm}^{(\mathrm{rm})}`$ $`=`$ $`{\displaystyle \frac{1}{4}}r_0^2+{\displaystyle \frac{R_{\mathrm{𝑚𝑝}}}{r_0}}+{\displaystyle \frac{m^21/4}{r_0^2}}`$ (9)
$`+`$ $`\left(n+{\displaystyle \frac{1}{2}}\right)\sqrt{3+4{\displaystyle \frac{m^21/4}{r_0^4}}}.`$
The result embodied by Eq. (9) has a clear physical interpretation. It contains a rotor-like contribution, $`m^2/(2𝒥)`$, with a moment of inertia given by $`𝒥=r_0^2/2`$ mpu, and a vibrational one characterized by a quantum number $`n`$. The vibrational frequency $`\omega _{\mathrm{𝑣𝑖𝑏}}=\sqrt{k/\mu }`$ ($`\mu =1/2`$) is given by the last square-root factor. Similar to atomic molecules, there is roto-vibrational coupling, since the vibration frequency depends on $`m`$ and, in addition, centrifugal distortion since $`r_0`$ also depends on $`m`$. For large enough values of $`R_{\mathrm{𝑚𝑝}}`$, implying large $`r_0`$ and therefore small average densities, the centrifugal distortion disappears and one has $`r_0(2R_{\mathrm{𝑚𝑝}})^{1/3}`$ for all $`m`$’s. In this limit rotational terms become negligible, as well as roto-vibrational ones. Thus, Eq. (9) reduces to a simple $`m`$-independent asymptotic expression

$$\stackrel{~}{\epsilon }_n^{(\mathrm{rm})}=\frac{3}{2^{4/3}}R_{\mathrm{𝑚𝑝}}^{2/3}+\sqrt{3}\left(n+\frac{1}{2}\right).$$
(10)
When adding the magnetic field, the roto-vibrational energy becomes

$`\epsilon _{nm}^{(\mathrm{rm})}`$ $`=`$ $`\stackrel{~}{\epsilon }_{nm}^{(\mathrm{rm})}+mW_{\mathrm{𝑚𝑝}}/2`$ (11)
$``$ $`{\displaystyle \frac{(m+W_{\mathrm{𝑚𝑝}}r_0^2/4)^2}{r_0^2}}+\left(n+{\displaystyle \frac{1}{2}}\right)\omega _{\mathrm{𝑣𝑖𝑏}},`$
in agreement with the expectations from Ref. Wen,  for two interacting electrons in a quantum ring. Pet  It is worth stressing that since the $`W_{\mathrm{𝑚𝑝}}`$ dependence only amounts to an energy shift of the Eq. (5) eigenvalue, the radial function $`u_{nm}(r)`$ does not depend on $`W_{\mathrm{𝑚𝑝}}`$. Therefore, one may conclude that the roto-vibrational properties are not affected by magnetic fields, for a fixed $`(n,m)`$ state. Of course, since the energy shift varies for different states, the magnetic field will modify the ordering of energy levels. For instance, the level crossings as a function of $`W_{\mathrm{𝑚𝑝}}`$ will cause the ground state to have a nonvanishing $`m`$-value. This actually explains the buildup of increasing permanent currents in the dot’s ground state.

The results from this section will be validated by comparing with the exact ones below. The roto-vibrational model presented here allows one to determine the crystallization onset from the criterion that rotation and vibration motions decouple when intrinsic-frame electron localization sets in. Conversely, when the coupling is strong the system could be represented by either a vibrating rotor or a rotating vibrator and, therefore, the situation can not be clearly resolved. It is also worth stressing that the roto-vibrational model describes all possible excitations of the relative-motion problem. For this particular system, this amounts to a description of all the excitations since the center-of-mass and spin degrees of freedom can be analytically integrated out.

## IV Onset of Wigner crystallization

### IV.1 Exact solutions

We have solved the radial equation for the relative problem, Eq. (5), numerically. Several standard approaches to eigenvalue problems with boundary conditions can be used for this purpose. Specifically, we have applied the so-called matching method where one integrates with the Numerov or Runge-Kutta algorithms from the origin $`r=0`$ outwards assuming an $`r^{|m|+1/2}`$ behaviour. Additionally, imposing an exponential decay law $`\mathrm{exp}(r^2/4)`$ for large $`r`$, inwards integration is performed and the two solutions are required to match at an intermediate $`r`$-point. To ascertain the numerical result for the eigenvalue, as a control we also used the method of ’node-counting’, where only outwards integration is performed and the eigenvalue is found from the condition that $`u_{nm}(r)`$ increases by one the number of radial nodes from the required value when the energy exceeds the correct eigenvalue by an infinitesimal amount. In principle, the node-counting method assures the correct boundary condition at $`r\mathrm{}`$ automatically. However, in practice, there always remains a small difference between the numerical and exact eigenvalues, responsible for a deviation from the exponential decay from some (large) $`r`$ onwards. In spite of this possible difference in asymptotic behavior, the two methods (matching and nodecounting) provide to a high accuracy the same eigenvalue.

With the above numerical methods an exploration of a part of the $`R_{\mathrm{𝑚𝑝}}W_{\mathrm{𝑚𝑝}}`$ plane has been performed. Figure 2 summarizes the results for ground-state energy and angular momentum. Note that the ground state always has $`(N,M)=(0,0)`$ for the CM quantum numbers and that, because of symmetry, even (odd) $`m`$ states are associated with singlet (triplet) total spin. Clear singlet-triplet oscillations, as studied in the literature (cf. Ref. Wag, ), are seen in the lower panel of Fig. 2. It is also worth mentioning that the ground-state energy contour lines are piece-wise linear as a result of the simple dependence on the $`W_{\mathrm{𝑚𝑝}}`$ parameter. Indeed, the energy is totally $`W_{\mathrm{𝑚𝑝}}`$-independent for $`m=0`$ and a fixed $`R_{\mathrm{𝑚𝑝}}`$.

Figure 3 displays the radial wave functions $`u_{00}(r)/\sqrt{r}`$ for two different values of $`R_{\mathrm{𝑚𝑝}}`$ as well as the corresponding effective potentials $`V_{\mathrm{𝑒𝑓𝑓}}`$ defined by Eq. (6). We note that, in agreement with the discussion of the preceding section, when increasing $`R_{\mathrm{𝑚𝑝}}`$ the effective-potential minimum moves outwards and it effectively binds the lower states to its neighborhood. When this occurs the radial probability is strongly quenched at small $`r`$ (lower panel) and the scenario indeed resembles the familiar one from the physics of diatomic molecules. A more detailed comparison of the roto-vibrational model with the exact results will be done in the next subsection.

### IV.2 Crystallization criterion

The validity of the roto-vibrational model presented in Sec. III is proved by the results of Fig. 4. As shown in the upper panel, the error of the analytical prediction for $`\stackrel{~}{\epsilon }_{nm}`$ is important only when $`R_{mp}`$ is small. As a matter of fact, for $`R_{mp}>2`$ the discrepancy for $`\stackrel{~}{\epsilon }_{00}`$ is always below 2%, even with the asymptotic expression Eq. (10). Note that the relative errors slightly increase with increasing $`n`$ (dashed curves), although the analytical approximation can still be considered quite good. The analogous comparison for $`m>0`$ (not shown) yields smaller relative errors than those for $`m=0`$. The lower panel of Fig. 4 analyzes the excitation energies as measured from $`\stackrel{~}{\epsilon }_{00}`$, again showing an excellent agreement between the exact results and the ones obtained with the analytical model, with small deviations only at small $`R_{\mathrm{𝑚𝑝}}`$. Based on the asymptotic expression, Eq. (10), we have chosen $`1/R_{\mathrm{𝑚𝑝}}^{2/3}`$ as an independent variable in order to better display the linear behavior associated with the rigid rotor at large $`R_{mp}`$.

We shall rely on the high accuracy of the roto-vibrational model to provide a quantitative measure of the crystallization onset. Our criterion will be the following: the two-electron parabolic system is assumed to be crystallized as a rotating Wigner molecule when the roto-vibrational coupling falls below a given percentage (typically chosen as 2 or 3 %). The roto-vibrational coupling is defined as

$$\gamma (m)=100\frac{\omega _{\mathrm{𝑣𝑖𝑏}}(m;R_{\mathrm{𝑚𝑝}})\sqrt{3}}{\sqrt{3}},$$
(12)
where $`\omega _{\mathrm{𝑣𝑖𝑏}}`$ is the square-root factor in Eq. (9) and $`\sqrt{3}`$ is the limit of this quantity for $`R_{\mathrm{𝑚𝑝}}\mathrm{}`$. Using a 3% condition the crystallization onset for each angular momentum is given by the $`R_{\mathrm{𝑚𝑝}}`$ value where each $`m`$ curve of Fig. 5 enters the shaded region. Note that the crystallization onset moves towards higher values as the angular momentum is increased, reflecting the property that roto-vibrational coupling is stronger for high-$`m`$ states. It is also worth mentioning that, since not all $`m`$-states are simultaneously crystallized, in practice the rotational bands will gradually degrade for increasing $`m`$’s as a consequence of the roto-vibrational coupling. In agreement with the bulk gas situation, the crystallized states with the proposed criterion are characterized by having a potential energy that largely exceeds the kinetic one, as can be easily checked from Eq. (9).

The crystallization properties of $`m`$ and $`m`$ states are identical since one can easily check that $`u_{nm}(r;R_{\mathrm{𝑚𝑝}})=u_{nm}(r;R_{\mathrm{𝑚𝑝}})`$. Therefore, taking into account the variations in ground state angular momentum we can draw the boundaries for Wigner crystallization in the $`R_{\mathrm{𝑚𝑝}}W_{\mathrm{𝑚𝑝}}`$ plane, i.e., the crystallization phase diagram (see the lower panel of Fig. 2). Of course, if instead of a 3% threshold for roto-vibrational decoupling one chooses a different value the crystallization onset will vary, although as shown in Fig. 5, for low $`m`$’s the crystallization is not crucially dependent on the precise percentage in the range 2 to 4%. Actually, it should be more appropriate to speak of crystallization onset for a given percentage of roto-vibrational decoupling than of an absolute value.

It should be noted that the onset of Wigner crystallization may be studied by means of the conditional probability distribution (CPD, cf. Refs. Yan, ; Mik, ) for finding one electron at positon $`𝐫_1`$ given that the second electron is at position $`𝐫_2`$. As soon as the interaction is switched on, the CPD, which for a system with only two electrons is just the modulus squared of the wave-function $`|\mathrm{\Psi }(𝐫_1,𝐫_2)|^2`$, exhibits the formation of a molecular-like state (see Fig. 6). It is difficult, however, to use this measure alone as a conclusive evidence of the formation of the Wigner molecule in two-electron dots. The difficulty arises from the fact that even weak interactions, for which we should not expect a crystallization, yield the formation of a hole around the electron at $`𝐫_2`$ (the correlation hole) and a maximum at $`𝐫_1=𝐫_2`$. Indeed, as seen from Fig. 6, the results for $`R_{mp}=1`$ hint at a molecular state even though according to our analysis the roto-vibrational coupling is still strong and this state is not associated with a crystallized phase (see Fig. 5). When $`R_{mp}`$ increases the depletion of the CPD around the fixed electron becomes much stronger, which is in qualitative agreement with the crystallization trend from our analysis.

## V Hartree-Fock and RPA approaches

In this section we discuss the results obtained within the symmetry unrestricted Hartree-Fock method for the ground state and the corresponding RPA for excitations. We recall that the HF and RPA approaches were devised for the analysis of many-body systems. Therefore, their application to a two-electron quantum dot is merely an exploratory approach to the qualitative features of the Wigner crystallization rather than a quantitative description of the above exact results.

We have solved the HF problem in the Fock-Darwin basis which diagonalizes the square bracket in Eq. (1), namely $`\{|a\eta ;a=1,\mathrm{}𝒩;\eta =,\}`$, where $`a`$ labels the orbital part and $`\eta `$ the spin. Our basis has been optimized such a way that we consider the 70 lowest Fock-Darwin states, of the non-interacting energy level scheme, for a chosen value of the magnetic field. An arbitrary single-particle orbital $`|i`$ is then expanded as $`|i=_{a\eta }B_{a\eta }^{(i)}|a\eta `$. In the chosen basis, the HF equations are written as a system of nonlinear eigenvalue equations for the matrix of $`B`$ coefficients (see, for instance, details in Ref. rpa, ).

We have imposed good $`s_z`$ HF orbitals leaving totally unspecified the remaining spatial symmetries. Note that the Slater determinant built with these single-particle orbitals will be an eigenstate of the total $`S_z`$ operator but not, in general, of $`𝐒^2`$. This, as we shall see in the results, is intimately connected with the prediction of broken circular symmetry.

In Fig. 7 we show the HF phase diagram in the $`R_{\mathrm{𝑚𝑝}}W_{\mathrm{𝑚𝑝}}`$ plane (lower panel) and the corresponding total energies (upper panel). The total energies resemble those obtained in the exact treatment, with approximate piecewise linear regions between orbital angular momentum transition lines. As expected the actual values at a given point in the diagram lie slightly above the corresponding exact results. In the lower panel different gray regions reflect a measurenot1  of the deviation from circularity of the ground state density, with the lightest intensity corresponding to a circular (nonbroken symmetry) solution and more intense gray levels to noncircular (broken symmetry) results. The contour lines show the total orbital angular momentum $`L_z`$. Regions I and III are of circular symmetry and for them $`L_z`$ has a good quantum number, taking the values 0 and -1, respectively. For the rest of the diagram (regions II, IV and V) the contour lines only indicate the expectation value of $`L_z`$ but, since the mean field in not circularly symmetric, this is no longer a good quantum number. The dotted curve separates the states having total spin projection $`S_z=0`$ (below) from $`S_z=1`$ (above).

For $`S_z=0`$ configurations, the HF method predicts a broken-symmetry solution when $`R_{\mathrm{𝑚𝑝}}`$ exceeds a value $`1`$, somehow below the onset of crystallization obtained in Sec. IV ($`R_{\mathrm{𝑚𝑝}}2.2`$). In this region (II) the HF solution is indeed a mixture of singlet and triplet states, as can be verified by computing the total spin dispersion $`\mathrm{\Delta }𝐒^2`$. The corresponding spatial density is built from two opposite and localized single electron orbitals. It is instructive to compare the HF solutions in region II with those obtained using a total-spin conserving ansatz

$`\mathrm{\Psi }_{\mathrm{𝑠𝑛𝑔}}(𝐫_1,𝐫_2,\eta _1,\eta _2)`$ $`=`$ $`\varphi (𝐫_1)\varphi (𝐫_2)\chi _{\mathrm{𝑠𝑛𝑔}}(\eta _1,\eta _2)`$ (13)
$`\mathrm{\Psi }_{\mathrm{𝑡𝑟𝑝}}(𝐫_1,𝐫_2,\eta _1,\eta _2)`$ $`=`$ $`𝒜[\varphi _1(𝐫_1),\varphi _2(𝐫_2)]`$ (14)
$`\times \chi _{\mathrm{𝑡𝑟𝑝}}(\eta _1,\eta _2)`$
for singlet ($`\mathrm{\Psi }_{\mathrm{𝑠𝑛𝑔}}`$) and triplet ($`\mathrm{\Psi }_{\mathrm{𝑡𝑟𝑝}}`$) states, where $`𝒜[\varphi _1(𝐫_1),\varphi _2(𝐫_2)]`$ denotes the antisymmetrized product of the two orbitals $`\varphi _1`$ and $`\varphi _2`$ while the $`\chi `$’s are the well known singlet and triplet spin states.

As a sample result the total energies obtained for $`(W_{\mathrm{𝑚𝑝}},R_{\mathrm{𝑚𝑝}})=(0,2)`$ are $`E_{\mathrm{𝑒𝑥𝑎𝑐𝑡}}=3.720`$, $`E_{HF}=4.034`$, $`E_{\mathrm{𝑠𝑛𝑔}}=4.185`$ and $`E_{\mathrm{𝑡𝑟𝑝}}=4.168`$; i.e., by requiring total $`𝐒^2`$ conservation the mean field energy raises considerably. In addition, while the HF solution breaks circular symmetry, both $`\mathrm{\Psi }_{\mathrm{𝑠𝑛𝑔}}`$ and $`\mathrm{\Psi }_{\mathrm{𝑡𝑟𝑝}}`$ yield circular densities because of the spatial dependence of the $`\varphi `$’s in Eqs. (13) and (14). We have also checked that these results are equivalent to those obtained using the Lipkin-Nogami projection method Lip ; Nog  for the effective Hamiltonian $`H_{\mathrm{eff}}=\lambda 𝐒^2`$ in order to restore $`𝐒^2`$ symmetry approximately. The above ansatz for states with good total spin are examples of the use of constraints in mean field approaches, which necessarily raise the energy above the mean field minimum. One could also project the symmetry-unrestricted HF orbitals as discussed in Ref. YanP, . In the latter case, however, the wave function is no longer a single Slater determinant but rather a sum of few determinants of the corresponding symmetry operator. Therefore, the ground state energy with the restored symmetry is no longer bound by the mean-field minimum and thus can be closer to the absolute minimum imposed by the variational principle.

With the above results, we conclude that the lowest HF solution in region II requires a simultaneous breaking of the spin and space symmetries. Taking into account the results of the preceding sections we can say that singlet-triplet mixing in region II is an artifact of the HF solution. However, space symmetry breaking in this region is a true physical mechanism indicating an intrinsic structural change in the exact wave function, as supported by the roto-vibrational model discussed above. We believe this peculiar combination of artifact and physics is due to the smallness of the configuration space for a two-electron system. We must also point out that at high $`R_{\mathrm{𝑚𝑝}}`$ and/or $`W_{\mathrm{𝑚𝑝}}`$ (regions IV and V) the HF prediction fails to match the results of the analysis given in Sec. IV.

In contrast to the exact results of Fig. 2 where a region (although narrow) with $`m=2`$ corresponding to a singlet state appears at large $`W_{\mathrm{𝑚𝑝}}`$, the HF approximation predicts only one $`(S_z=0)(S_z=1)`$ transition as $`W_{\mathrm{𝑚𝑝}}`$ increases. This can be understood as an overestimation of the exchange energy in the HF model which tends to favor spin alignment whenever orbital overlapping occurs, as it does in region III with circular orbitals, as well as in regions IV and V with two-lobed orbitals.

We consider next the results obtained by solving the RPA for excitations. As discussed in Ref. rpa,  the RPA determines the moment of inertia associated with the collective rotation of a deformed HF structure (see Eq. (18) in Ref. rpa, ). Figure 8 presents the evolution with $`R_{\mathrm{𝑚𝑝}}`$, at $`W_{\mathrm{𝑚𝑝}}=0`$, of the RPA moment of inertia $`𝒥_{RPA}`$ (circles). For comparison, the values computed through the solution of Eq. (7) of the roto-vibrational model, $`𝒥=r_0^2/2`$, are also shown (solid lines). Each $`𝒥`$-line starts at the crystallization onset for the corresponding angular momentum using the criterion of Sec. IV.B. Note that $`𝒥_{RPA}`$ remains null until the HF solution breaks rotational symmetry at $`R_{\mathrm{𝑚𝑝}}1`$, from there on it reasonably agrees with the exact values, somehow averaging the exact results for different $`m`$’s. The molecule stretching, yielding larger $`r_0`$’s ($`𝒥`$’s) as $`m`$ increases, is obviously outside the RPA. All $`𝒥`$ values slowly converge to a common result with increasing $`R_{\mathrm{𝑚𝑝}}`$, i.e., to an exact rigid-rotor behavior. For comparison, the dashed line represents the asymptotic value corresponding to $`r_0(2R_{\mathrm{𝑚𝑝}})^{1/3}`$. The nice qualitative agreement between $`𝒥_{RPA}`$ and $`𝒥`$ is giving additional support to the above conclusion that space symmetry breaking in region II of the HF phase diagram indicates a genuine physical effect and, thus, it also supports the overall picture of a rotating Wigner molecule.

Although the RPA restores circularity on the deformed HF mean field,rpa  associated with the $`L_z`$ operator, it is not able to restore the symmetry related to $`𝐒^2`$, since the latter one is a two-body operator which is beyond the RPA treatment of broken symmetries. A side effect of the HF spin artifact in region II is that when spin-flip bosonic pairs are included in the RPA, generalizing our previous calculation,rpa  the correlation energy in this symmetry-broken phase is badly overestimated, i.e., it is between two and three times larger than the exact value. This does not occur, however, in the circularly symmetric regions.

Focusing now on the RPA vibron states, we must also distinguish excitations associated with spin and space (charge) degrees of freedom. While the RPA qualitatively describes all excitations in regions I and III, it fails for spin excitations in phase II. Obviously, this is due to the HF spin artifact in this region. It also fails for space excitations of phases IV and V. Finally, we end this section by pointing out that, in spite of the spin artifact of phase II, in pure triplet phases RPA reproduces the exact spin precession frequencies known from the theory of magnetic resonance.Sli  That is, a pure spin-flip state (precessional mode) is expected at the Larmor energy $`\mathrm{}\omega _{}=g^{}\mu _BB=\frac{1}{2}g^{}m^{}W_{\mathrm{𝑚𝑝}}`$ mpu. Indeed, within the RPA, i.e., in the quasi-boson approximation,BR ; rpa  one finds

$$[,O_{}^+]=\mathrm{}\omega _{}O_{}^+,$$
(15)
with the vibron operator for the Larmor mode

$$O_{}^+=\frac{S_x+iS_y}{\sqrt{2S_z}}.$$
(16)
In Eq. (16) $`S_z`$ is the HF expectation value of the $`S_z`$ operator. In fact, the Larmor mode at $`\mathrm{}\omega _{}`$ appears whenever the ground state has $`S_z0`$ and it is normally the lowest excitation of the system.

## VI Summary

We performed a systematic study of the evolution of ground and excited states of two-electron quantum dots subject to an external magnetic field. The analysis has been done in terms of magneto-parabolic units and the associated parameters $`(R_{\mathrm{𝑚𝑝}},W_{\mathrm{𝑚𝑝}})`$ that give, respectively, the ratios of Coulomb interaction strength and cyclotron frequency to effective confinement. The ground state calculations are summarized in a phase diagram that can be equally applied to a variety of confinements, magnetic fields and material parameters.

We suggested an analytical model for the interpretation of the exact results, including roto-vibrational coupling and centrifugal distortion (molecule stretching). Within this roto-vibrational model we proposed a criterion to determine the onset of Wigner crystallization based on the decoupling of rotational and vibrational motions. For a 3%-decoupling threshold we found that Wigner crystallization appears, for zero-angular-momentum states, when $`R_{\mathrm{𝑚𝑝}}`$ exceeds a value $`2`$. States with larger $`m`$’s crystallize at higher $`R_{mp}`$ values. In agreement with the homogeneous gas situation the potential energy of the crystallized states is much larger than the kinetic energy, the latter one being solely due to the vibrational zero-point motion of the electrons.

The HF calculations predict that crystallization for $`S_z=0`$ occurs when $`R_{\mathrm{𝑚𝑝}}>1`$, the new phase (II) being in an artificial mixture of singlet and triplet spin states. The space symmetry breaking in phase II is a genuine physical effect but the spin mixture is an artifact due to the smallness of the configuration space for a two-electron system. Other HF symmetry-broken phases (IV and V) do not agree with the exact results. The RPA moment of inertia qualitatively agrees with the result from the roto-vibrational model, although the molecule centrifugal distortion is missed. On the other hand, the RPA produces reliable results for space (charge) excitations in regions I, II and III, as well as for spin excitations when the HF solutions possess good $`𝐒^2`$ and $`S_z`$ quantum numbers (regions I and III). We would expect a broader applicability of the many-body theories (HF+RPA) for larger systems. Work along this line is in progress. In conclusion, the combined use of exact and model calculations allowed us to ascertain the existence of a rotating Wigner molecule in a two-electron dot for relatively large electron densities or, equivalently, small $`R_{\mathrm{𝑚𝑝}}`$ parameters.

###### Acknowledgements.

This work was supported by Grant No. BFM2002-03241 from DGI (Spain). R. G. N. gratefully acknowledges support from the Ramón y Cajal programme (Spain).


text:
# A High-resolution Spectrum of the R CrB Star V2552 Ophiuchi

## 1 Introduction

The R Coronae Borealis variables are a select group of distinctive supergiants. A Galactic census presently recognizes no more than about 30 R CrB stars. Their especial distinctions are a severe deficiency of hydrogen in their atmosphere, and a propensity to fade by up to 6 - 8 magnitudes at irregular and unpredictable times. The origin of the hydrogen deficiency (and other composition anomalies) and the mechanism responsible for the characteristic fading of the star remain the subjects of theoretical and observational enquiries. Given the small number of R CrB stars, one seeks to enlarge the sample in the hope that a new discovery may yield a fruitful clue to outstanding questions.

The latest addition to the R CrB census is V2552 Oph (Kazarovets et al. 2003). This star CD -22 12017 was identified as a R CrB star from a light curve showing two R CrB-like fadings (Kato & Katsumi 2003). This identification was confirmed by Hesselbach, Clayton, & Smith (2002) from photometry and a low resolution spectrum. In this paper, we report on the first high-resolution optical spectrum of V2552 Oph, and show that it is remarkably similar to the spectrum of the eponym with a few striking differences.

## 2 Observations

Our spectrum of V2552 Oph was obtained on June 11, 2003 (JD 2452801.809) with the McDonald Observatory’s 2.7m telescope and the 2dcoudé echelle spectrograph (Tull et al. 1995) at a resolving power of 60,000. The spectrum runs from about 3800 Å to 10,000 Å with gaps between the orders beyond 5800 Å. The signal-to-noise (S/N) ratio in the continuum is 95 or greater over much of the observed wavelength interval longward of 5500Å. The S/N ratio decreases to shorter wavelengths. V2552 Oph appeared to be at maximum light. Standard reductions of the CCD frames give the radial velocity of V2552 Oph as 60.5$`\pm 0.9`$ km s-1 from a set of 51 unblended photospheric lines.

## 3 V2552 Oph and R CrB: Spectra and Abundances

By inspection of the spectra, it was apparent that V2552 Oph has a spectrum remarkably similar to that of R CrB. This is illustrated in Figures 1 and 2. The similarity suggests that the stars have similar atmospheres and compositions. Lines which are of very similar strength in V2552 Oph and R CrB belong to the following atoms and ions: H i, C i, Na i, Mg i, Al i, Si i, Si ii, S i, K i, Ca i, Sc ii, Ti ii, Fe i, Fe ii, Ni i, Zn i, Ba ii, and La ii. The similarity for C i is a direct consequence of the fact that the neutral carbon atom is the principal contributor of continuous opacity (Schönberner 1975). The similarity for other atoms and ions with a mix of low to high excitation lines is an indication that the effective temperature ($`T_{\mathrm{eff}}`$), surface gravity ($`g`$), microturbulence ($`\xi _t`$), and composition are also similar for the pair of stars.

Obvious differences in line strength for atoms and ions are noted which imply differences in abundance:

* The Li i 6707 Å resonance doublet is strong in R CrB but absent from V 2552 Oph (Figure 1).
* The N i lines, all of high-excitation, are considerably strengthened in V2552 Oph. Figure 1 shows the N I lines at 6706.1Å and 6708.7Å to be enhanced in V2552 Oph relative to R CrB.
* The O i high excitation and the low excitation 6363Å \[O i\] lines are weaker in V2552 Oph. The fact that the low and high excitation are all weaker suggests that the oxygen abundance is lower for V2552 Oph.
* Lines of Sr ii, Y ii, and Zr ii are much stronger in V2552 Oph (Figure 2).

The lines of the C2 Swan bands are much weaker in V2552 Oph than in R CrB. The rotational excitation of the C2 molecules is representative of the photosphere.

An abundance analysis of V2552 Oph and R CrB was made using the family of model atmospheres employed by Asplund et al. (2000). The R CrB spectrum is one obtained in 1995 August with the 2dcoudé at the same spectral resolution as our spectrum of V2552 Oph. In light of the spectral similarities, we adopted as the initial (and final) model one with the atmospheric parameters chosen by Asplund et al. for R CrB: $`T_{\mathrm{eff}}=6750`$ K, $`\mathrm{log}g`$ = 0.5 cgs, and $`\xi _t=7`$ km s-1, and an assumed ratio C/He = 0.01 by number of atoms. Standard spectroscopic checks for $`T_{\mathrm{eff}}`$, $`g`$, and $`\xi _t`$ confirm the adopted parameters to be appropriate for V2552 Oph. Results of the abundance analysis are given in Table 1 where we list results for V2552 Oph and R CrB from our spectra, and results for R CrB from Asplund et al. Selection of lines and atomic data (i.e., the $`gf`$-value of a line) are taken largely from Asplund et al. For a few Fe i lines, $`gf`$-values were taken from Reddy et al. (2003). In general, the same lines were used for both stars and, hence, the abundance analyses are accurate in a differential sense. There is very good agreement between the two R CrB analyses; the difference seen for oxygen and barium is due to an improved selection of lines in the case of oxygen and to an erroneous equivalent width for one Ba ii line used in Asplund et al.’s analysis.

Table 1 shows that for many elements there is an uncanny match of the abundances for V2552 Oph and R CrB: differences are 0.05 dex or less in many cases. This degree of similarity is surely fortuitous. Abundance differences for N, Y, and Zr average about 0.6 dex in favor of V2552 Oph and the difference for O is 0.6 dex in favor of R CrB.

## 4 V2552 Oph among the R CrBs

If classified by composition, the F-G type R CrB stars belong to one of two classes - a majority and a minority class (Lambert & Rao 1994). The minority R CrB stars are marked by an approximately solar Si abundance and a severe Fe deficiency or, equivalently, an extraordinarily high Si/Fe \- see Lambert & Rao (1994) and Asplund et al. (2000) for a fuller characterization of the minority class. Majority members have a solar or slightly subsolar Si abundance and a mild Fe deficiency resulting in a less unusual Si/Fe ratio (\[Si/Fe\] $``$ +0.5). V2552 Oph with \[Si/Fe\] = $`+`$0.5 and a 1 dex Fe underabundance becomes the newest member of the majority.

Within the majority, there are significant abundance differences for some elements, notably those few elements for which V2552 Oph and R CrB differ: Li, N, O, Sr, Y, and Zr. Except for a coupling of Sr, Y, and Zr, these elements do not vary in a correlated manner. Among majority R CrBs previously analysed, Y Mus and FH Sct closely resemble V2552 Oph in composition across the elements from H to La. In particular, the trio have no detectable Li, similar low abundances of O, and high Sr, Y, and Zr abundances (relative to Fe and to Ba and La). We refer the reader to Asplund et al. (2000) for a discussion of what may be learnt from the compositions about the origins of these variable stars.

## 5 Interstellar and circumstellar gas

Absorption components of likely interstellar origin are present in the Na D lines (Figure 3). Stellar Na D lines at +60 km s-1 are accompanied by strong components extending from -10.5 to 0.0 km s-1. The strongest components are seen in the K i 7699Å line at -7.1 and -0.6 km s-1. Diffuse interstellar bands are present in V2552 Oph’s spectrum at 5780, 6284, 6379, and 6613Å. Adopting the wavelengths given by Herbig (1995), we calculate a heliocentric velocity of about -7 km s-1 from these bands. Herbig’s (1993) relation between the equivalent width of a diffuse interstellar band and reddening gives EB-V = 0.4 for V2552 Oph.

This estimate, the apparent visual magnitude (V = 10.7, Kato & Katsumi 2003) combined with the absolute magnitude (MV = -5) for a F-type R CrB star in the LMC (Alcock et al. 2001) puts V2552 Oph at a distance of about 7.7 kpc. At this distance, interstellar gas following galactic rotation has a heliocentric velocity of about -11 km s-1. This velocity shifts to the red for smaller distances: it is -2.3 km s-1 for a distance of 1.1 kpc. The similarity between this range of predicted velocities and the measured velocities suggests that the Na D and K i sharp components and the diffuse interstellar bands at velocities of -10 to 0 km s-1 are interstellar in origin.

In addition to the interstellar components, there are sharp components in the Na D lines at -38.4, -52.9, and a weak one at -69.0 km s-1. These are possibly of circumstellar origin. If so, they represent clouds expanding away from the star at -98, -113, and -129 km s-1. These velocities exceed the surface escape velocity for a R CrB star. If it is assumed that the two stronger components represent gas ejected at the time of the deep double minimum reported by Kato & Katsumi (2003), the radial distances traveled are 93, and 135 stellar radii, respectively, where a stellar radius is estimated at 80 solar radii. Such high velocity ejections are commonly seen following a deep minimum.

## 6 Circumstellar dust

Declines of the R CrB stars are caused by soot clouds in front of the star. A characteristic of the stars is an infrared excess from circumstellar dust (Feast et al. 1997). Many R CrB stars show an infrared excess at the K band. Most have an infrared excess at longer wavelengths. A few have a strong excess at 60$`\mu `$m and 100$`\mu `$m from very cold dust in an extended shell (Rao & Nandy 1986; Walker 1986).

To assess the infrared excess of V2552 Oph, we searched the 2MASS All-sky Catalog of Point Sources (Cutri et al. 2003) for V2552 Oph’s presence at J, H, and K. A source was found coincident to within 0.4 arc sec with V2552 Oph’s position on an ESO R band frame. The 2MASS magnitudes are J = 8.69, H = 8.39, and K = 8.16. The J-H and H-K color indices, after correction for interstellar reddening corresponding to EB-V = 0.41, lie on the blackbody color-color plot at about the temperature of the photosphere, i.e., there is no infrared excess detectable at J, H, and K. V2552 Oph could not be located in the IRAS Faint Source Catalog.

V2552 Oph is one of few R CrB stars not to show an excess in these bands. In this respect, it is similar to Y Mus (Feast 1997). Kilkenny & Whittet (1984) did find an infrared excess for Y Mus at 5$`\mu `$m and 10$`\mu `$m. The absence for Y Mus of an infrared excess at K and shorter wavelengths has been linked to the fact that the last minimum occurred several decades ago, i.e., recent soot production has ceased or been reduced to low levels. This explanation is not easily transferred to V2552 Oph which has experienced one or two fadings recently. It is presumably nothing more than a coincidence that V2552 Oph resembles Y Mus both in a lack of infrared excess at K and in composition.

## 7 Concluding Remarks

Addition of V2552 Oph to the rare class of R CrB variables is attributable to the Japanese amateur astronomer K. Haseda whose announcement identifiying the star as a variable included the suggestion that it might be a R CrB variable.111 http://www.kusastro.kyoto-u.ac.jp/vsnet/Mail/alert6000/msg00226.html Our high-resolution spectrum fully confirms the identification of V2552 Oph from photometry and low-resolution spectroscopy that the star is a R CrB variable. We show that it is a member of the majority class of these variables, as introduced by Lambert & Rao (1994). Within that class, V2552 Oph closely resembles Y Mus and FH Sct in showing a below average oxygen abundance and an above average abundance of Sr, Y, and Zr. In this respect, V2552 Oph confirms the spectrum of compositions offered by the majority stars but provides no novel clues to the origins of the anomalous compositions of these remarkable variables. Perhaps, the next discovered R CrB variable will be the Rosetta Stone needed to unlock the chain tying the abundances to the origins of these stars.

We thank B.E. Reddy for help with the acquisition and reduction of the spectrum of V2552 Oph, G. Pandey for the model atmosphere, and D. Yong for the preparation of the figures.

TABLE I
Elemental Abundances for V2552 Oph and R CrB

|  | R CrB |  | V2552 Oph | n a | Sun b |
| --- | --- | --- | --- | --- | --- |
| Species | Asplund et al. | present |  |  |  |
| H i | 6.9 | 6.86 | 6.66 | 1, 1 | 12.00 |
| Li i | 2.8 | 2.55 | $`<`$0.96 | 1, 1 | 3.35 |
| C i | 9.2 | 9.11 $`\pm `$ 0.25 | 9.05 $`\pm `$ 0.24 | 16, 17 | 8.46 |
| N i | 8.4 | 8.42 $`\pm `$ 0.44 | 8.96 $`\pm `$ 0.28 | 5, 8 | 7.90 |
| O i | 9.0 | 8.60 $`\pm `$ 0.37 | 7.96 $`\pm `$ 0.21 | 6, 6 | 8.76 |
| Na i | 6.1 | 6.11 $`\pm `$ 0.05 | 6.14 $`\pm `$ 0.12 | 3, 3 | 6.37 |
| Mg i |  | 6.81 $`\pm `$ 0.26 | 6.77 $`\pm `$ 0.09 | 2, 2 | 7.62 |
| Al i | 5.8 | 5.76 $`\pm `$ 0.13 | 5.78 $`\pm `$ 0.11 | 2, 2 | 6.54 |
| Si i | 7.2 | 6.97 $`\pm `$ 0.21 | 6.97 $`\pm `$ 0.22 | 4, 5 | 7.61 |
| Si ii |  | 7.34 $`\pm `$ 0.15 | 7.67 $`\pm `$ 0.28 | 2, 2 |  |
| S i | 6.8 | 6.70 $`\pm `$ 0.13 | 6.77 $`\pm `$ 0.09 | 7, 7 | 7.26 |
| K i |  | 4.77 | 4.48 | 1, 1 | 5.18 |
| Ca i | 5.3 | 5.32 $`\pm `$ 0.33 | 5.34 $`\pm `$ 0.27 | 5, 4 | 6.41 |
| Sc ii |  | 2.89 $`\pm `$ 0.30 | 2.91 $`\pm `$ 0.23 | 3, 3 | 3.15 |
| Ti ii |  | 4.05 $`\pm `$ 0.30 | 4.11 $`\pm `$ 0.37 | 2, 2 | 5.07 |
| Fe i | 6.5 | 6.40 $`\pm `$ 0.24 | 6.37 $`\pm `$ 0.27 | 20, 21 | 7.54 |
| Fe ii |  | 6.40 $`\pm `$ 0.13 | 6.42 $`\pm `$ 0.09 | 10, 10 |  |
| Ni i | 5.5 | 5.49 $`\pm `$ 0.11 | 5.55 $`\pm `$ 0.18 | 4, 5 | 6.29 |
| Cu i |  |  | 3.95 $`\pm `$ 0.16 | 2, 2 | 4.34 |
| Zn i |  | 4.16 | 4.31 | 1, 1 | 4.70 |
| Sr ii |  |  | 4.28 | 1 | 2.99 |
| Y ii | 1.5 | 1.55 $`\pm `$ 0.18 | 2.29 $`\pm `$ 0.09 | 3, 3 | 2.28 |
| Zr ii |  | 1.84 $`\pm `$ 0.19 | 2.25 $`\pm `$ 0.19 | 3, 3 | 2.67 |
| Ba ii | 1.6 | 1.13 $`\pm `$ 0.16 | 1.03 $`\pm `$ 0.09 | 3, 3 | 2.25 |
| La ii |  | 0.64 $`\pm `$ 0.38 | 0.63 $`\pm `$ 0.43 | 2, 2 | 1.25 |

a n = number of lines used in the analyses: (a, b) where a and b are number

of lines used for R CrB and V2552 Oph, respectively.

b Recommended abundances for the solar system from Lodders

(2003, her Table 2).