text:
# Gate errors in solid state quantum computer architectures

## Abstract

We theoretically consider possible errors in solid state quantum computation due to the interplay of the complex solid state environment and gate imperfections. In particular, we study two examples of gate operations in the opposite ends of the gate speed spectrum, an adiabatic gate operation in electron-spin-based quantum dot quantum computation and a sudden gate operation in Cooper pair box superconducting quantum computation. We evaluate quantitatively the non-adiabatic operation of a two-qubit gate in a two-electron double quantum dot. We also analyze the non-sudden pulse gate in a Cooper-pair-box-based quantum computer model. In both cases our numerical results show strong influences of the higher excited states of the system on the gate operation, clearly demonstrating the importance of a detailed understanding of the relevant Hilbert space structure on the quantum computer operations.

The rapid progress of quantum information theory has prompted intense interest in developing appropriate hardwares for quantum computation and communication . Since quantum coherence needs to be maintained for a quantum computer (QC) to maximize its capability, it has to satisfy stringent requirements . One requirement is that a QC should have a well-delineated Hilbert space, while another requirement is that the quantum bits (qubit) should have long coherence time (i.e. slow decoherence). Indeed, decoherence can be loosely defined as irreversibly losing information from the operational Hilbert space into the environment, which is the rest of the universe from the perspective of a QC system. Furthermore, the requirement of slow decoherence needs to be satisfied not only when the qubits are standing alone, but also when the qubits are interacting with each other or are manipulated externally for single or multi-qubit operations.

Many physical systems with suitable two-level quantum dynamics have been suggested as possible candidates for qubits in a variety of QC architectures. The proposals based on solid state systems, such as the ones based on electron spins , nuclear spins , and superconductors , are particularly attractive because of their potential for scalability. However, solid state structures also present complex environments and fast decoherence rates . For example, for electron spins trapped in quantum dots (QD), their environment includes the underlying crystal lattice in terms of phonons, surrounding electron and nuclear spins, impurities and interfaces, unintentionally trapped charges, and many other defects invariably present in a solid state environment. Even the orbital degrees of freedom of the electron spin qubits themselves are part of the environment. In the case of QC based on Cooper pair boxes (CPB), the environment includes defects in the Josephson junctions, quasiparticles, electromagnetic couplings, etc. From the perspective of qubit Hilbert space, the complex environment in a solid means considerable difficulty in identifying the boundary. The fuzziness in the qubit environment also implies that a gate operation that is not optimized (a gate error) can cause not only incomplete rotations in the qubit Hilbert space, but also leakage into the environment. Thus, to ensure successful operation of a QC, we have to clearly understand how each of the environmental elements affects the QC coherent evolution, and clarify the effects of imperfections in gate operations. In this paper we will investigate specifically the latter effect.

Some of the QC proposals satisfy the requirement for precise coherent manipulation by employing resonant means such as resonant Rabi oscillation and resonant Raman scattering . However, in most solid state QC schemes, non-resonant operations  are crucial components of the QC operation. In the spin-based schemes, adiabatic switching of exchange interaction is a cornerstone of the two-qubit operations, while in the CPB-based scheme, sudden switching of a voltage bias has been proposed and used to perform single-qubit operations. Here we study whether QC system integrity can be faithfully maintained if the exchange interaction switching is not adiabatic in the QD case, or the voltage pulse is not switched on suddenly in the CPB case. More specifically, we would like to explore the gate operation threshold when leakage due to non-resonant operations is sufficiently small so that quantum error correction codes can be effectively applied. A quantitative and qualitative theoretical investigation of this issue is clearly important in understanding the feasibility of proposed QC architecture.

Let us first examine the spin-based quantum dot quantum computer (QDQC) in GaAs. Here effects of non-adiabaticity in the single-qubit operations should be relatively weak because bulk spin-orbit coupling at the bottom of GaAs conduction band is small. If the Rashba spin-orbit coupling can also be minimized for the GaAs heterostructure, the overall spin-orbit coupling should be very small so that applied magnetic field pulses will not be able to mix electron spin and orbital states. For two-qubit operations, the adiabatic condition can be satisfied more easily with larger single electron excitation energy and larger on-site Coulomb interaction, which are both consequences of small QDs. However, the size of a gated quantum dot is limited from below by gate and device dimensions. Thus we would like to quantitatively assess the adiabatic condition for two-qubit operations in a double dot of realistic dimensions. In particular, we would like to determine the minimum duration of two-qubit operations constrained by adiabaticity and the quantum error correction threshold.

To deal with a system described by a slowly-time-dependent Hamiltonian, the adiabatic approximation can be employed, in which the state of a quantum system at time $`t`$ can be expanded on the instantaneous eigenstates at that time : $`\psi (t)=_ic_i(t)u_i(t)`$ and $`H(t)u_i(t)=E_i(t)u_i(t)`$. Substituting these two equations into the time-dependent Schrödinger equation, the time evolution of the coefficients $`c_i(t)`$ can be obtained:

$$\frac{c_k(t)}{t}=\underset{ik}{\overset{N}{}}\frac{c_i(t)}{E_k(t)E_i(t)}k|\frac{H(t)}{t}|i\mathrm{exp}\left\{\frac{1}{i\mathrm{}}_{\mathrm{}}^t(E_i(\tau )E_k(\tau ))𝑑\tau \right\}.$$
(1)
In our calculation for the exchange gate of a QDQC, we consider a two-electron double quantum dot . The basis states $`u_i`$ are the two-electron eigenstates of the double dot. The explicit time-dependence is in the inter-dot barrier height. In other words, our complete confinement potential is given by $`V(𝐫,t)=V_L(𝐫)+V_R(𝐫)+V_b(𝐫,t)`$, where $`V_L(𝐫)`$ and $`V_R(𝐫)`$ are the fixed single quantum dot confinement potentials, while $`V_b(𝐫,t)`$ is a Gaussian barrier potential that takes the form $`V_b(𝐫,t)=\left[V_{max}(V_{max}V_{min})e^{\left(\frac{t}{\tau }\right)^2}\right]e^{\left(x^2/l_x^2+y^2/l_y^2\right)}`$. Here $`V_{max}`$ and $`V_{min}`$ are the highest and lowest barriers respectively, $`l_x`$ and $`l_y`$ defines the width and thickness of the potential barrier between the two QDs located at $`x=\pm a`$ on the $`y`$ axis, and $`\tau `$ is the half width of the Gaussian temporal variation of the barrier height. We have chosen to use Gaussian lineshapes for the confinements, barrier, and time dependence without loss of generality. Indeed, as we point out below, non-Gaussian-shaped pulses might produce more favorable conditions for the operation of two-qubit gates. Now the matrix elements of the time derivative of the Hamiltonian can be expressed completely in terms of the potential barrier matrix elements: $`i|\frac{H(t)}{t}|j=i|\frac{V_b(t)}{t}|j=\left(i|V_{max}|ji|V_b(t)|j\right)\frac{2t}{\tau ^2}.`$ The wavefunction coefficients can then be calculated by plugging these matrix elements into Eq. (1).

We perform a numerical calculation for the evolution of a two-electron double dot using Eq. (1). Initially the system is entirely in the ground singlet state. In other words, $`c_i=\delta _{i0}`$ for all $`i`$. The variation of the potential barrier will not lead to singlet-triplet coupling. Both types of states have about the same excitation gap (smaller of the single particle excitation energy and the on-site Coulomb charging energy) and similar energy spectra, and the exchange coupling does not lead to transitions between singlet and triplet states, so that a pure singlet initial state should give a reliable description of adiabaticity for a general case. By varying the half width $`\tau `$ of the barrier variation time, we quantitatively evaluate the change in the ground singlet state population, thus obtaining a lower limit to the gate operating time using the criterion of quantum error correction code threshold. In the direct integration of Eq. (1), spectrum information such as the barrier matrix elements and eigen-energies are needed. We do not attempt to calculate the state spectrum at all barrier heights, as it would have made the whole computation intractable. Instead, we perform the spectrum calculation at a series of potential barrier heights (at 10 points which correspond to 19 points on the Gaussian temporal profile of the barrier variation), then obtain the rest of the information through interpolation. The energy spectra we use are for a double dot with Gaussian confinements of 30 nm radii and 40 nm inter-dot distance . The energy barrier height $`V_b`$ ranges between 14 meV and 35 meV, corresponding to exchange splitting of 280 $`\mu `$eV to 3.3 $`\mu `$eV.

The result of our calculation is plotted in Fig. 1(a). The leakage (y-axis) is defined as $`1|c_0|^2`$ which is zero before the gate is applied, and should be zero if the gate is perfectly adiabatic. The curve has several interesting features. At short gating times the leakage decreases quite fast, with the straight line in Fig. 1 giving an indication of the speed of this decrease. This decrease rate slows down as the gating time increases. At very long time, adiabatic approximation should be asymptotically exact, so that the $`k`$th excited state occupation takes the form

$$|a_k|^2=\frac{4\mathrm{sin}^2\omega _{k0}t/2}{\mathrm{}^2\omega _{k0}^4}\left|k|\frac{V_b}{t}|0\right|^2.$$
(2)
The magnitude of the probability amplitude $`a_k`$ is in the order of the ratio between the change in the Hamiltonian in one Bohr period $`2\pi /\omega _{k0}`$ and the energy difference $`\mathrm{}\omega _{k0}`$. Thus the leakage should be proportional to $`1/T^2`$ at large $`T`$ if we qualitatively represent $`H/t`$ by $`\mathrm{\Delta }H/T`$ where $`T`$ characterizes the total gating time.

The leakage we obtained oscillates quite rapidly with the gating time. This oscillation is due to the phase factors $`\mathrm{exp}\left\{\frac{i}{\mathrm{}}_{\mathrm{}}^t\mathrm{\Delta }E_{ki}(\tau )𝑑\tau \right\}`$ for all the levels involved. As the gating time is varied, the contribution to the leakage from different higher levels changes rapidly, leading to the fast oscillations.

Figure 1(a) demonstrates that for gating time ($`2\tau `$) longer than 30 ps, leakage in our double dot system should be sufficiently small ($`10^6`$) so that the currently available quantum error correction schemes would be effective. On the other hand, an exchange splitting of 0.1 meV corresponds to about 20 ps gating time for a swap gate  (with rectangular pulse) at the shortest. Therefore, adiabatic condition does not place an extra burden on the operation of the two-qubit gates such as a swap. There is in general no need to significantly increase the gating time in order to accommodate the adiabatic requirement.

Figure 1(b) shows the energy splittings between the ground singlet and the 1st and 2nd excited singlet states. A finite gap between the ground and excited states exists for all barrier heights, thus guaranteeing the relatively easy satisfaction of adiabatic condition in this configuration. The anti-crossing in this figure shows a competition between S-P hybrid states (one electron in the single dot S states while the other in the single dot P states ) and double occupied states (both electrons are in the S state of the same dot) in the double dot. Essentially, this competition indicates the relative importance of single particle excitation (S-P transition) and on-site Coulomb interaction. At high barrier height, the S-P hybrid state has lower energy than the double occupied state, while at lower barrier the 1st excited singlet state is a mixture of the double occupied state and the S-P hybrid state. As the inter-dot barrier is lowered, the two electrons can tunnel between the two dots more easily, thus reducing the effect of on-site Coulomb repulsion.

The duration of the gate voltage pulse that changes the barrier height between two neighboring QDs is not the only parameter that determines whether adiabaticity is satisfied. Pulse shape tailoring is extensively used in optical control experiments , and is presumably applicable in the control of interdot barrier. Our choice of a Gaussian pulse shape is by no means optimal. Therefore, the result shown here provides an upper limit to the adiabatic condition for the adopted parameters in the spin-based QDQC architecture.

One condition for the adiabatic approximation to be valid is that the system has no level crossing during the evolution. In the double QD system, there should be no orbital degeneracy as we lower the inter-dot barrier. The situation bears some similarity to the transition from a hydrogen molecule to a helium atom, although in the operation of a QDQC we would never push the two QDs so close and the barrier so low that the two-electron double dot becomes a two-electron single dot. Here the state spectrum essentially maintains the same structure throughout the gate operation. It would, however, be interesting to explore the state spectrum as the double dot is squeezed towards a single QD, and study whether degeneracy might play a role during such evolution. After all, in a QDQC, QDs are used to label the qubits; while in a two-electron QD, labeling is impossible as the two electrons are indistinguishable.

Let us now turn our attention to the Cooper pair box quantum computer (CPBQC) model. The Hamiltonian of a CPB can be written on the basis of charge number states of the box:

$`H`$ $`=`$ $`4E_C(\widehat{n}n_g)^2E_J\mathrm{cos}\widehat{\varphi }`$ (3)
$`=`$ $`{\displaystyle \underset{n}{}}\left[4E_C(nn_g)^2|nn|{\displaystyle \frac{E_J}{2}}(|nn+1|+|n+1n|)\right],`$ (4)
where $`E_C`$ is the charging energy of a Cooper pair box (CPB), $`E_J`$ is the Josephson coupling between the CPB and an external superconducting lead, $`n_g`$ represents the applied voltage on the CPB in terms of an effective charge number, and $`n`$ refers to the number of Cooper pairs in the box. The stationary solutions of a CPB form energy bands. The two levels for a CPB qubit are the lowest two energy levels at $`n_g=1/2`$, where the eigenstates are approximately $`|=(|0+|1)/\sqrt{2}`$ and $`|=(|0|1)/\sqrt{2}`$ with a splitting of about $`E_J`$.

Most theoretical studies of CPBQC employ a two-level description of the system . However, as we have learned  in the case of QDQC above, higher excited states can play an important role in the dynamics of a qubit, especially when it is subjected to non-resonant operations. One particular operation for CPBQC is the sudden pulse gate to shift $`n_g`$, thus bringing a system from a pure ground state ($`|0`$ at, e.g. $`n_g=1/4`$) to a coherent superpositioned state ($`(|+|)/\sqrt{2}`$ at $`n_g=1/2`$). However, such a simple description of the pulse gate is only valid when $`E_J/E_C0`$. Since $`E_J`$ determines the gate speed of a CPBQC, such a condition is not practical for a realistic QC. Below we explore the situation when $`E_J`$ is not vanishingly small, and study how the pulse gate operation will be affected.

Before considering the pulse gate, we first study the state compositions of a CPB in order to find the limit of a two-level model and show some analytical results in Table I (For these expressions, we include only the lowest excited states). As is shown in Table I, no matter which $`n_g`$ value we choose, there are always unwanted states that mix into the composition of the eigenstates with weights proportional to $`r=E_J/E_C`$. Such mixtures will affect the dynamics of the system. For example, if we prepare a CPB in the ground state $`\psi _0`$ at $`n_g=1/4`$ and apply a sudden pulse gate to shift $`n_g`$ up to $`1/2`$, the state right after the shift will not be $`(|+|)/\sqrt{2}`$, but $`(|+|)/\sqrt{2}+\frac{r}{4}(||)/\sqrt{2}`$. Thus, finite $`r`$ will lead to errors in the pulse gate to the first order of $`r`$.

In real experiments, the pulse gate always has finite rise/fall times (non-sudden). In Ref. , the pulse rise time is in the range of 30 to 40 ps. Such gradual rise and fall of the pulse gate inevitably lead to more errors, which have been considered in the context of two-level systems . Here we calculate the fidelity of the pulse gate taking into account the finite rise/fall time, the higher excited states, and the complete composition of all the eigenstates. The eigenstates are numerically obtained (instead of perturbatively ) on the basis of 21 charge number states ($`|10`$, $`\mathrm{}|0`$, $`|1`$, $`\mathrm{}|10`$). We prepare the initial state as the ground state $`\psi _0`$ near $`n_g=1/4`$, then let the CPB evolve under Hamiltonian (4). The pulse gate represented by $`n_g`$ is turned on with a lineshape ranging from linear to quintic power of time during the rise and fall: $`n_g(t)=n_g(0)+\left(\frac{2t}{t_r}\right)^m(n_g(t_r)n_g(0))/2`$ for $`t[0,t_r/2]`$ and $`n_g(t)=n_g(t_r)\left[\frac{2(t_rt)}{t_r}\right]^m(n_g(t_r)n_g(0))/2`$ for $`t[t_r/2,t_r]`$, where $`t_r`$ is the pulse rise/fall time and the pulse is switched on from $`0`$ to $`t_r`$. For system parameters, we choose the values used in Ref. , with $`E_C=117`$ $`\mu `$eV, $`E_J=51.8`$ $`\mu `$eV, $`n_g(0)=0.255`$, and $`n_g(t_r)=0.5`$.

In Fig. 2 we plot the state fidelity as a function of rise time. Here the state fidelity is defined as the maximum probability for the CPB to be in the first excited state after the pulse gate when $`n_g`$ returns to $`0.255`$. Figure 2 illustrates two important points. First of all, the fidelity of the pulse gate decreases oscillatorily instead of monotonously as the pulse rise time increases. The oscillations (with periods of 30 to 40 ps) in the curves represent the coherent evolution of the CPB during the rise/fall of the pulse voltage: at $`n_g=1/4`$, the dominant energy scale is $`2E_C234\mu `$eV with a corresponding period of 17 ps, while at $`n_g=1/2`$ the dominant energy scale is $`E_J=51.8\mu `$eV with a corresponding oscillation period of 76 ps. The observed oscillation period is essentially an average of these two limits. Secondly, including higher excited states is important in correctly evaluating the fidelity dependence on pulse rise time, as a comparison of the two panels of Fig. 2 clearly demonstrates, especially for higher power of time-dependence. The real pulse shape may not be the same as any of the power law time-dependence we consider here, but the disagreement between the two panels in Fig. 2 clearly shows again the necessity of taking excited states into consideration when evaluating non-resonant gates such the pulse gate here.

Our current study is in essence a quantitative delineation of the Hilbert space for a QDQC (CPBQC) during a two-qubit (single-qubit) operation. Such a study is necessary for other proposed QC schemes as well in order to select the optimal operating regime, eliminate possible errors, and enhance signal to noise ratios in experiments. We have restricted ourselves to $`T=0`$ in our calculations, but our results remain valid at finite (low) temperatures as long as substantial thermal occupation of higher excited levels does not occur in the system.

We thank David DiVincenzo and Belita Koiller for useful discussions, and ARDA and LPS for financial support.


format: text

text:
# Overlap Among States at Different Temperatures in the SK Model

## Abstract

We discuss the issue of temperature chaos in the Sherrington–Kirkpatrick spin glass mean field model. We numerically compute probability distributions of the overlap among (equilibrium) configurations at two different values of the temperature, both in the spin glass phase. The situation on our medium size systems is clearly non-chaotic, but a weak form of chaos could be emerging on very large lattices.

Two years after our paper  on the subject (that we will quote as (A) in the following) we are coming back to the problem of temperature chaos in spin glasses and, more widely, in disordered and complex systems. A problem that is, we believe, still a very open one.

In these two years the problem of temperature chaos has been studied under many new lights. One can fairly say that indications are mixed, with some preference for a no-chaos scenario on medium size systems: the detailed discussion of  (where arguments against the possibility of a strong chaos picture are given) is probably a perfect starting point for the reader interested in the details of the subject.

Perturbation theory in an expansion below $`T_c`$ (in mean field theory ) when pushed to the fifth order shows absence of chaos through highly non-trivial cancelations, although one finds no general feature that could imply that these cancelations will be present at all orders in perturbation theory. The naive TAP equations for the mean field Sherrington–Kirkpatrick model, when solved numerically on systems with order of $`10^2`$ spins, also lead to exclude the presence of temperature chaos .

Bouchaud and Sales have shown that there is no chaos in the REM model (but if one sits exactly at $`T_c`$) . On the contrary Sales and Yoshino have discussed in  the case of DPRM (directed polymers in random media), and have shown that there is temperature chaos in this model, and that a temperature perturbation plays a role very similar for example to a perturbation in the potential: i.e. in the case of DPRM temperature perturbation looks generic and creates chaos. The recent work by Sasaki and Martin  describes situations where chaos is present.

These recent studies add new elements to the many former studies of an interesting problem first discussed by Parisi , and then studied in many other works (see among others). Still the presence or the absence of temperature chaos is one of the few open problems remaining in the Sherrington–Kirkpatrick model.

Let us start by reminding which were the main results of our first paper (A). There we discussed the behavior of the two-temperature overlap

$$q_{T_1,T_2}^{(2),(N)}\overline{\left(\frac{1}{N}\underset{i=1}{\overset{N}{}}\sigma _i^{T_1}\sigma _i^{T_2}\right)^2}$$
(1)
for systems with $`N`$ spins. The over-line is for the average over the quenched disorder, the brackets are for the thermal average, $`\{\sigma _i^{T_1}\}`$ is an equilibrium configuration of Ising spins at temperature $`T_1`$, while $`\{\sigma _i^{T_2}\}`$ is at temperature $`T_2`$. We have considered the Sherrington–Kirkpatrick spin glass mean field model, a diluted finite connectivity mean field spin glass model, and the $`3D`$ Edwards–Anderson spin glass. In a $`T`$-chaotic situation we would expect $`q_{T_1,T_2}^{(2),(N)}`$ to go to zero in the infinite volume limit, as soon as $`T_1T_2`$. We found however that for all models we studied (on systems with up to $`4096`$ sites) $`q_{T_1,T_2}^{(2),(N)}`$ was not small for $`T_2T_1`$ finite and reasonably large. To be more precise we found that in our data we had always

$$q_{T_1^{(min)},T_2}^{(2),(N)}q_{T_2,T_2}^{(2),(N)}>0,\text{ for }T_cT_2>T_1^{(min)}0.4T_c,$$
(2)
even if the value of the l.h.s. was decreasing with increasing volumes. The validity of the relation (2) is suggestive of a very non-chaotic situation (it would hold for an usual ferromagnet), even if from the numerical data of (A) it was clear that asymptotically it may well be violated (that is in any case not enough to imply a $`T`$-chaotic behavior).

So (A) was suggesting a clear absence of chaotic behavior on medium size systems, while showing that for increasing lattice sizes signs of a (modestly) more chaotic behavior were possibly starting to appear. We note here again that systems with a not huge number of spins can be relevant to the physics of spin glasses. The experiments of reference  show that the number of spin involved in the collective behavior observed in a typical spin glass experiment is of the order of $`10^5`$, that is not so far away from the number of spins we can handle numerically today. Accordingly, even if it should emerge that $`T`$-chaos is asymptotically present for very large volumes, the present numerical simulations could turn out to describe reasonably well the experimental regime.

In order to get further informations on the sensitivity of the spin glass phase with respect to temperature, we have decided to measure the full probability distribution of the overlap of configurations equilibrated at different values of the temperature: $`P_{T_1,T_2}(q)`$ (In (A) we restricted ourselves to the second moment of this distribution). In terms of $`P_{T_1,T_2}(q)`$, temperature chaos just means that $`lim_N\mathrm{}P_{T_1,T_2}(q)=\delta (q)`$, for any $`T_1T_2`$ . A more complete information like the one contained in the full probability distribution of the order parameter can indeed allow a more detailed analysis of the scaling behavior, making possible to uncover a wider range of phenomena.

Since the three models studied in (A) were showing a very similar behavior, we focus here on one of them, namely the Sherrington–Kirkpatrick model.

The numerics are very similar to the ones of (A), and we refer to  for a discussion of the details of the simulations. We use binary quenched random couplings, $`J=\pm 1`$. Here we have studied systems with $`N=64`$, $`N=256`$, $`N=1024`$ and $`N=4096`$ spins, down to $`T=0.4=0.4T_c`$. We have used a multi spin (different spins are encoded in the same computer word) version of the parallel tempering Monte Carlo algorithm. Two copies of the system at each one of a set of temperature values ($`0.4,0.4+\mathrm{\Delta }T,0.4+2\mathrm{\Delta }T\mathrm{}`$, with $`\mathrm{\Delta }T=0.025`$ but for $`N=4096`$ where $`\mathrm{\Delta }T=0.0125`$) are brought to equilibrium (in the same realization $`J`$ of the couplings) and are used to compute the $`P_{T_i,T_j}^{(J)}(q)`$ (more precisely the subset with temperatures $`0.40,0.45,0.50,\mathrm{}`$) . An average over $`1024`$ different realizations of the quenched disorder ($`256`$ for $`N=4096`$) is taken to compute the average $`P_{T_i,T_j}(q)`$. This is an order of magnitude more disorder samples than in . For each realization we have performed $`400K`$ sweeps for equilibrium plus $`1000K`$ sweeps for measurements ($`520K`$ only for $`N=4096`$).

We show at first the probability distributions themselves (normalized by $`_1^1P(q)𝑑q=1`$). This distributions have been symmetrized, although they are quite symmetric, even sample by sample. In each of figures 1, 2, 3 and 4 we show three $`P_{T_1,T_2}(q)`$: the non-diagonal one that gives the probability of the overlap of configurations at $`T=0.4`$ (the lower temperature we equilibrate) with configurations at $`T=0.6`$, and the two diagonal probability distributions with $`T_1=T_2=0.4`$ and $`T_1=T_2=0.6`$ respectively. In figure 1 we have data from our smallest lattice with $`N=64`$, in figure 2 from $`N=256`$, in figure 3 from $`N=1024`$, and in figure 4 we have data from the largest lattice, with $`N=4096`$.

On the smallest lattice (figure 1) the peak of the non-diagonal $`P(q)`$ is higher than the one of $`P_{0.6,0.6}`$, and its position is basically at half way between the positions of the maxima of $`P_{0.4,0.4}`$ and $`P_{0.6,0.6}`$. This is typically what would happen in a ferromagnet. For increasing volumes the position of this peak is moving (very slowly) to lower values of $`q`$, but does not seem to approach $`q=0`$ (we will discuss in detail this point in the following).

In figure 5 we show in the same plot the four non-diagonal $`P_{0.4,0.6}`$ for $`N=64`$, $`256`$, $`1024`$ and $`4096`$. These $`P_{0.4,0.6}(q)`$ are normalized and are drawn on the same scale, so that a visual comparison of the four curves is meaningful. Here we see that the $`q0`$ peak of $`P_{0.4,0.6}(q)`$ do increase as function of $`N`$ (but for a negative analysis of this statement see later in the text and figure 7), and that its position shifts only very little towards $`q=0`$. The mass carried by the distribution $`P(q)`$ in $`q0`$ seems to be increasing on the largest lattice, but this effect is affected by a large statistical error. Figure 4 shows indeed wiggles in $`P_{0.4,0.4}(q)`$ that are known not to exist in the infinite volume limit, and are accordingly to be attributed to the limited number of disorder samples. The bump popping up around $`q=0`$ in figure 5 is compatible with the statistical error. It is not visible in our data for larger $`T_1`$ values ($`0.45,0.50,\mathrm{}`$).

In figure 6 we plot the values of $`q`$ where $`P(q)`$ is maximum, $`q_{\text{max}}`$, for the diagonal and non-diagonal $`P(q)`$ at different $`N`$ values. The two points at $`N=\mathrm{}`$ for $`P_{0.4,0.4}(q)`$ and $`P_{0.6,0.6}(q)`$ have been obtained  using the method of : they show how reliable is a linear fit in $`N^{\frac{1}{3}}`$ of our data for the diagonal distributions. The points for $`P_{0.4,0.6}(q)`$ show some curvature, but the effect is not dramatic, and a non-zero value in the limit $`N\mathrm{}`$ is surely favored.

As we have already shown, the peaks of the (temperature) non-diagonal $`P(q)`$ increase when $`N`$ increases, but not as fast as the ones of the diagonal $`P(q)`$. In figure 7 we show the value of $`P(q_{\text{max}})`$ for the different $`N`$ values and for the three $`P(q)`$ we are discussing. The height of the peaks of the two diagonal $`P(q)`$ increases exactly like $`N^{\frac{1}{3}}`$ as it should. The two continuous straight lines are the best linear fit to the diagonal data: the fits turn out to be very good.

From figure 7 the situation of $`P_{0.4,0.6}(q)`$ looks quite different. Here the growth slows down at $`N=1024`$ and is really small on our largest lattice size with $`N=4096`$. We do not believe that one can draw precise quantitative conclusions from figure 7: a (weak) non-chaotic picture could survive in the infinite volume limit, or chaos could appear very slowly (only on very large volumes). Let us spell clearly, in any case, that the mechanism that will govern a possible appearance of chaos will be based on the $`P_{T_1,T_2}(q)`$ having two $`q=0`$ and and at $`q=q^{}>0`$, and by having the $`q=0`$ peak growing for increasing volume at the expense of the $`q^{}`$ peak. In any case our numerical results appear to support the idea that temperature chaos will not be effective on the typical sizes that are relevant in spin glass experiments .

In order to be more quantitative we have also looked at the ratio of the mass carried by $`P_{0.4,0.6}(q)`$ close to $`q=0`$ as compared to the mass carried by the large $`q`$ region. We define

$$k\frac{_m^{+m}P(q)𝑑q}{_{q_{max}}^1P(q)𝑑q},$$
(3)
where $`k`$ is a normalization constant, and we take $`m=0.05`$ (a very similar picture would be obtained for any $`m`$ not too large). $``$ decreases with $`N`$ if the mass of the peaks at large $`q`$ dominates, while it increases when the dominating contribution is the one at $`q0`$.

We plot $``$ for $`P_{0.4,0.6}(q)`$ versus $`N^{\frac{1}{3}}`$ in figure 8. $``$ is constant on the smaller lattices, but starts increasing on the largest lattice size. The error due to sample to sample fluctuations is here very large (the simulation at $`N=4096`$ have been very costly in computer time), and the growth of $``$ is not significant at more than two standard deviations, but an effect is very plausibly there. Again, this is probably an indication toward the fact that temperature chaos could eventually emerge on very large systems.

It is worth noticing that if temperature chaos is (slowly) emerging when increasing $`N`$ this is probably not happening with the position of the two peaks at large $`|q|`$ shifting to $`q=0`$, but with a third peak in $`q=0`$ emerging and eventually becoming the only contribution, and accordingly a discontinuity in $`q_{\text{max}}`$ as function of $`N^{1/3}`$ (Figure 6).

The problem of temperature chaos, already at the mean field level, is turning out to be a hard problem: this is true both for analytical and for numerical computations. Here we have provided some further hints about the behavior of the system in the infinite volume limit: the very large scale, state of the art simulations we have been able to run, give some suggestions, probably hinting in favor of a very weak chaos that would emerge only for very large lattice sizes.

We acknowledge many enlightening discussions with Giorgio Parisi and Silvio Franz. We thank Andrea Crisanti and Tommaso Rizzo for providing us with the $`N=\mathrm{}`$ data for $`q_{\text{max}}`$ for the diagonal $`P(q)`$ of figure 6. The numerical simulations have been run on the CEA Compaq Alpha Server at Grenoble.


text:
# I Introduction

## I Introduction

One of the most interesting achievements in the field of boson condensation was the experimental observation of binary mixture of trapped condensates \[1-4\]. The species in the mixture can be atoms with different F spin orientations \[3-4\] or simply different hyperfine states of the same atom \[1-2\]. In the last case we can have inter-conversion between the components through the coupling of the atoms with an external laser field .

In the case of uncoupled components many effects have been predicted theoretically and determined experimentally such as spatial phase separation , stability properties of the equilibrium states \[6-7\] and symmetry breaking effects . These aspects of binary condensate mixtures have been treated in the literature using various methods such as the Thomas-Fermi (TF) approximation \[8-9\] and numerical solutions of the appropriate coupled Gross-Pitaevskii (GP) equations .

In this paper we use the variational method \[11-13\] to study the equilibrium properties of a coupled two-component BEC. In our exploratory calculation we use the one-dimensional model of reference  and take as variational parameters the population fraction, translation and scaling transformation of the equilibrium state densities, all assumed to have a Gaussian shape. We consider the intensity of the interaction between the condensates equal and the detuning is put equal to zero, leading to equal equilibrium population fraction. In our calculation we investigate the behavior of the spatial phase separation as a function of the trap displacement and the strength of the laser field. The coupling with the laser field have a stabilizer effect in the process of phase separation, opposite to the one coming from the trap displacement.Our paper is organized as follows:In section 2 we show briefly how to use the variational method to study the equilibrium and stability properties of a binary mixture of condensates.Specifically,we derive the general form of the equations of motion and develop a scheme to find the normal modes.In section 3 we present our numerical results and discuss its physical significance.A summary of our numerical results is presented in section 4.

## II Variational study of the equilibrium and stability properties of a binary mixture of condensates.

### A  The variational approach

The starting point of our discussion is the action

$$S=𝑑z𝑑t\left[\underset{j}{}i\mathrm{}\psi _j^{}(z,t)\dot{\psi }_j(z,t)(z,t)\right]$$
(1)
where $`\psi _j(z,t)`$, $`j=a,b`$ are the condensate wave function of each component in the mixture and $`(z,t)`$ is the energy density,

$$(z,t)=\underset{j}{}\psi _j^{}\left[\frac{\mathrm{}^2}{2m}\frac{d^2}{dz^2}+V_{trap}^j\right]\psi _j+\frac{1}{2}\underset{k,j}{}\lambda _{kj}|\psi _j|^2|\psi _k|^2+\mathrm{\Omega }\left[\psi _a^{}\psi _b+\psi _b^{}\psi _a\right].$$
(2)
In the above expression $`V_{trap}^j(z)`$, $`j=a,b`$ are the trapping potential of each component

$$V_{trap}^j=\frac{m}{2}\omega _z^2(z+\gamma _jz_0)^2$$
(3)
where $`\gamma _a=1`$ and $`\gamma _b=1`$, $`z_0`$ is the trap displacement, $`\lambda _{aa}=\lambda _{bb}`$ and $`\lambda _{ab}`$ are the strength of the intraspecies and interspecies interactions, respectively, and the last term comes from the coupling with the laser field responsible by the interspecies tunnelling, with $`\mathrm{\Omega }`$ the Rabi frequency.

Imposing that the action (1) is stationary with respect to a variation of $`\psi _j(z,t)`$ subject only to the normalization constrain $`_j|\psi _j(z,t)|^2=N`$, leads to the coupled time-dependent Gross-Pitaevskii equations for the condensate wave functions ,

$`i\mathrm{}{\displaystyle \frac{\psi _a}{t}}`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2}{z^2}}+V_{trap}^{(a)}+\lambda _{aa}|\psi _a|^2+\lambda _{ab}|\psi _b|^2\right)\psi _a+\mathrm{\Omega }\psi _b`$ (4)
$`i\mathrm{}{\displaystyle \frac{\psi _b}{t}}`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2}{z^2}}+V_{trap}^{(b)}+\lambda _{ab}|\psi _a|^2+\lambda _{bb}|\psi _b|^2\right)\psi _b+\mathrm{\Omega }\psi _a`$ (5)
To get the equations for the stationary states, we look for solutions of equations (4) and (5) of the form $`\psi _j(z,t)=e^{i\frac{\mu t}{\mathrm{}}}\psi _j(z)`$ ,which gives rise to the time-independent GP equations:

$`\mu \psi _a`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{dz^2}}+V_{trap}^{(a)}+\lambda _{aa}|\psi _a|^2+\lambda _{ab}|\psi _b|^2\right)\psi _a+\mathrm{\Omega }\psi _b`$ (6)
$`\mu \psi _b`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{dz^2}}+V_{trap}^{(b)}+\lambda _{ab}|\psi _a|^2+\lambda _{bb}|\psi _b|^2\right)\psi _b+\mathrm{\Omega }\psi _a`$ (7)
We can view the two condensate wave functions as components of a spinor of “quasi-spin” equal to $`\frac{1}{2}`$  which leads immediately to the property that the stationary equations are invariant under a transformation which is the product of a space reflection, $`P_z`$, and an atom exchange, $`P_{ex}`$, where $`P_{ex}=i\mathrm{exp}i\frac{\pi }{2}\sigma _x`$. Therefore, as $`(P_zP_{ex})^2=P_zP_{ex}`$, we can classify the stationary states as even (gerade) and odd (umgerade) under this transformation. In the even class $`\psi _a^g(z)=\psi _b^g(z)`$, whereas in the odd class $`\psi _a^u(z)=\psi _b^u(z)`$. The next step is to solve numerically the coupled stationary equations for the condensate wave functions to map the solutions as a function of $`\mathrm{\Omega }`$ and $`z_0`$. However, before embarking on this task, we found it worthwhile to adopt a simpler approach, the variational method. In the variational method the search for stationary states reduces to finding the stationary points in a finite dimensional energy surface which is much simpler than the corresponding search in an infinite dimensional space, when we solve exactly the coupled GP equations. Besides, the variational solutions can be used as an initial guess in the numerical solution of the coupled GP equations .

In the variational approach we parametrize the time dependence of the condensate wave function through a set of $`2d`$ parameters \[12-13\] which we denote by $`𝐰=(\mathrm{w}_1,\mathrm{w}_2,\mathrm{},\mathrm{w}_{2d})`$,

$$\psi _j(z,t)=\psi _j(z,𝐰(t))$$
(8)
When we replace the condensate wave functions parametrized as in equation (8), in equation (1), the action reduces to a “classical” action, whose variation leads to Hamilton-type equations of motion in terms of these parameters

$$\underset{l}{}\mathrm{\Gamma }_{kl}(𝐰)\dot{\mathrm{w}}_l=\frac{E}{\mathrm{w}_k}(𝐰),$$
(9)
where $`E(𝐰)`$ is the spatial integral of the energy density (2), with the condensate wave functions parametrized as in equation (8)

$$E(𝐰)=𝑑z(z,𝐰)$$
(10)
and the antisymmetric matrix $`\mathrm{\Gamma }_{kl}(𝐰)`$ is given by

$$\mathrm{\Gamma }_{kl}(𝐰)=2\mathrm{Im}\underset{j}{}𝑑z\frac{\psi _j^{}}{\mathrm{w}_k}(z,𝐰)\frac{\psi _j}{\mathrm{w}_l}(z,𝐰)$$
(11)
It follows from equation (9) that the equilibrium configurations are determined by the equations

$$\frac{E}{\mathrm{w}_k}(𝐰^0)=0,k=1,2,\mathrm{},2d.$$
(12)
### B  Normal Modes

To investigate the stability of the equilibrium configurations, we calculate the energies of the normal modes. They are stable if the energies are real and positive and unstable if one of the energies is complex. To find the normal modes, we linearize the equations of motion in the neighborhood of an equilibrium configuration, leading to:

$$\underset{l}{}\mathrm{\Gamma }_{kl}\dot{\overline{\mathrm{w}}_l}(t)=\underset{l}{}\mathrm{H}_{kl}\overline{\mathrm{w}}_l.$$
(13)
In this equation $`\overline{𝐰}`$ are the displacements from equilibrium $`𝐰=𝐰^0+\overline{𝐰}`$ and $`\mathrm{H}_{kl}`$, $`\mathrm{\Gamma }_{kl}`$ are, respectively, the Hessian and the matrix $`𝚪(𝐰)`$ evaluated at the equilibrium configuration, that is

$$\mathrm{H}_{kl}\mathrm{H}_{kl}(𝐰^0)=\frac{^2E}{\mathrm{w}_k\mathrm{w}_l}(𝐰^0)$$
(14)
and

$$\mathrm{\Gamma }_{kl}=\mathrm{\Gamma }_{kl}(𝐰^0)$$
(15)
To proceed, we divide the $`2d`$ parameters into two groups, $`𝐰=(q_1,q_2,\mathrm{},q_d,p_1,p_2,\mathrm{},p_d)=(𝐪,𝐩)`$, of coordinates $`𝐪`$ and momenta $`𝐩`$. Schematically, the reasoning behind this splitting is that the amplitude of the condensate wave function, whose square is the condensate density, depends only on the coordinates while its phase, whose gradient is the velocity field, depends on the momenta and the coordinates.

In terms of the coordinates and the momenta the equations of motion (13), read

$$𝚪\left(\begin{array}{c}\dot{\overline{𝐪}}\hfill \\ \dot{\overline{𝐩}}\hfill \end{array}\right)=𝐇\left(\begin{array}{c}\overline{𝐪}\hfill \\ \overline{𝐩}\hfill \end{array}\right).$$
(16)
where $`𝚪`$ and $`𝐇`$ are, respectively, antisymmetric and symmetric matrices which can be written in terms of four $`d\times d`$ blocks

$$𝚪=\left(\begin{array}{cc}𝚪_{qq}\hfill & 𝚪_{qp}\hfill \\ 𝚪_{pq}\hfill & 𝚪_{pp}\hfill \end{array}\right),𝐇=\left(\begin{array}{cc}𝐇_{qq}\hfill & 𝐇_{qp}\hfill \\ 𝐇_{pq}\hfill & 𝐇_{pp}\hfill \end{array}\right)$$
(17)
where, for example, $`𝚪_{qq}`$ and $`𝐇_{qq}`$ are $`d\times d`$ matrices whose elements are given by equations (11) and (14), where the derivatives are with respect to the coordinates, with an analogous definition for the other $`d\times d`$ matrices.

In our method to find the normal modes, we try to stay as close as possible to the one adopted in the standard case . To begin with, we should find a transformation to a set of new coordinates and momenta

$$\left(\begin{array}{c}\overline{𝐪}\hfill \\ \overline{𝐩}\hfill \end{array}\right)=𝐓^1\left(\begin{array}{c}𝐐\hfill \\ 𝐏\hfill \end{array}\right).$$
(18)
such that

$$𝚪^1\mathrm{𝐇𝐓}^1=𝐓^1\left(\begin{array}{cc}𝐎\hfill & 𝚲\hfill \\ 𝚲\hfill & 𝐎\hfill \end{array}\right)$$
(19)
where $`𝚲`$ is a diagonal $`d\times d`$ matrix whose diagonal elements are the normal mode energies.

In terms of the new coordinates and momenta, the equations of motion, (16), reduce to

$$\left(\begin{array}{c}\dot{𝐐}\hfill \\ \dot{𝐏}\hfill \end{array}\right)=\left(\begin{array}{cc}\mathrm{𝟎}\hfill & 𝚲\hfill \\ 𝚲\hfill & \mathrm{𝟎}\hfill \end{array}\right)\left(\begin{array}{c}𝐐\hfill \\ 𝐏\hfill \end{array}\right).$$
(20)
leading to

$$\dot{Q}_k=\mathrm{\Lambda }_kP_k,\dot{P}_k=\mathrm{\Lambda }_kQ_k$$
(21)
The equations (18)-(19) define the transformation to the normal modes and our scheme to find the matrix $`𝐓^1`$ is a straightforward generalization of the standard case . The starting point is to solve the eigenvalue problem

$$𝚵^1𝐇V^{(k)}=\mathrm{\Lambda }_kV^{(k)}$$
(22)
where $`𝚵`$ is the hermitian and antisymmetric matrix $`𝚵=i𝚪`$

As in the standard case, this eigenvalue problem has two properties: (i) if $`V^{(k)}`$ is an eigenvector with eigenvalue $`\mathrm{\Lambda }_k`$, then $`V^{(k)}`$ is also an eigenvector with eigenvalue $`\mathrm{\Lambda }_k^{}`$. (ii) The eigenvectors with different eigenvalues are $`𝚵`$ orthogonal, that is $`V^{(l)+}𝚵V^{(k)}=0`$, if $`\mathrm{\Lambda }_l\mathrm{\Lambda }_k^{}`$.

If one of the eigenvalues of equation (22) is complex, the system is unstable. If all the eigenvalues are real it is stable and we can find the transformation to the normal modes as follows. We define a matrix $`𝐒^1`$, whose first $`d`$ columns are the $`2d`$ components of the $`d`$ eigenvectors with real positive eigenvalues and the next $`d`$ columns, the $`2d`$ components of the corresponding eigenvectors with real negative energies.

In terms of $`𝐒^1`$, the eigenvalue equations (22) reads

$$𝚵^1\mathrm{𝐇𝐒}^1=𝐒^1\left(\begin{array}{cc}𝚲\hfill & \mathrm{𝟎}\hfill \\ \mathrm{𝟎}\hfill & 𝚲\hfill \end{array}\right)$$
(23)
and the orthogonality of the eigenvectors gives that

$$𝐒^1𝚵𝐒^1=\left(\begin{array}{cc}\mathrm{𝟏}\hfill & \mathrm{𝟎}\hfill \\ \mathrm{𝟎}\hfill & \mathrm{𝟏}\hfill \end{array}\right)$$
(24)
The matrix $`𝐓^1`$ is related to $`𝐒^1`$ by

$$𝐓^1=𝐒^1𝐔$$
(25)
where $`𝐔`$ is the unitary matrix

$$𝐔=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}\mathrm{𝟏}\hfill & 𝐢\hfill \\ \mathrm{𝟏}\hfill & 𝐢\hfill \end{array}\right)$$
(26)
To summarize, we state the main steps in our procedure: (i) First we determine the equilibrium configurations by solving the $`2d`$ equations (12). (ii) Next, for each equilibrium configuration, we solve the eigenvalue equation (22). (iii) When the configuration is stable, we construct the matrix $`𝐒^1`$ from the eigenvectors as indicated above and $`𝐓^1`$ as shown in equation (25). When the system oscillates in a normal mode only one pair $`(Q_k,P_k)`$ is different from zero and equation (18) gives the time evolution of the system in this case.

### C  Gaussian ansatz

Now that we have established the general framework of our calculations, we turn to our specific application where the variational parameters are related to the population fraction, translation and breathing shape oscillation of the condensate densities. Thus, the condensate wave functions are written as

$$\psi _j(z,𝐰)=e^{iF_j(z,𝐰)}A_j(z,𝐰)$$
(27)
where the amplitude of the condensate wave function is parametrized as

$$A_j(z,𝐰)=\sqrt{Nn_j}\frac{1}{\sqrt{\pi ^{1/2}q_{2j}}}e^{\frac{(zq_{1j})^2}{2q_{2j}^2}}$$
(28)
and the phase as

$$F_j(z,𝐰)=p_{1j}(zq_{1j})+\frac{p_{2j}}{2q_{2j}}(zq_{1j})^2+\theta _j$$
(29)
The parameters introduced in the above expression have the following physical interpretation: $`n_j`$ is the population fraction of each condensate,

$$Nn_j=|\psi _j(z,𝐰)|^2𝑑z,$$
(30)
$`q_{1j}`$ is the center of mass of the spatial distribution of each condensate,

$$Nn_jq_{1j}=z|\psi _j(z,𝐰)|^2𝑑z,$$
(31)
and $`q_{2j}/\sqrt{2}`$ is the width of the spatial distribution of each condensate

$$Nn_j\frac{q_{2j}^2}{2}=(zq_{1j})^2|\psi _j(z,𝐰)|^2𝑑z,$$
(32)
The momenta $`p_{1j},p_{2j}`$ and $`\theta _j`$ appear in the phase of the condensate wave functions. The $`p_{1j}`$ is the center of mass moment of each condensate

$$Nn_jp_{1j}=i\psi _j^{}(z,𝐰)\frac{\psi _j}{z}(z,𝐰)𝑑z$$
(33)
and $`p_{2j}`$ is connected to the expectation value of the dilatation operator in the center of mass frame of each condensate

$$\frac{1}{2}Nn_jp_{2j}q_{2j}=\frac{i}{2\mathrm{}}\left(\psi _j^{}(z,𝐰)\frac{\psi _j}{z}(z,𝐰)\frac{\psi _j^{}}{z}(z,𝐰)\psi _j(z,𝐰)\right)(zq_{1j})𝑑z,$$
(34)
showing that $`q_{1j},p_{1j}`$ is related to the translational degrees of freedom and $`q_{2j},p_{2j}`$ to the breathing shape oscillation.

Since $`E(𝐰)`$ depends only on the phase difference between the two condensates, the equations of motion (9) give that the total number of particles is conserved, $`n_a(t)+n_b(t)=1`$. This reduces the number of degrees of freedom to ten, which we take as the relative phase $`\theta =\theta _b\theta _a`$ and the relative population fraction $`n=\frac{n_bn_a}{2}`$, in addition to the $`q_{1j},q_{2j},p_{1j},p_{2j},j=a,b`$.

## III  Numerical Results

To perform the calculations we specify the model parameters . Following reference , we consider our system to be 87Rb and take $`N=2.3\times 10^4`$, $`\nu _z=60`$Hz which gives $`z_{sho}=\sqrt{\mathrm{}/m\omega _z}=1.4\mu `$m and we put all the interaction strength equal to $`\lambda _{aa}=\lambda _{bb}=\lambda _{ab}=17.5`$ $`\mu `$m sec$`{}_{}{}^{1}\mathrm{}`$, in order to reproduce qualitatively the spatial distribution of the condensates shown in FIG.2 of reference .

Our first task is to solve the ten equilibrium equations (12). Six of the equations lead immediately to $`p_{1j}^0=p_{2j}^0=0,j=a,b`$, equal equilibrium population fraction $`n^0=0`$, and $`\theta ^0=0`$ or $`\pi `$. The solution with $`\theta ^0=0`$ belongs to the even class under $`P_zP_{ex}`$ and the ones with $`\theta ^0=\pi `$, to the odd class. We restrict the calculations to the odd class since, as $`\mathrm{\Omega }>0`$, the lowest energy configuration necessarily belongs to this class. Besides, to be an eigenstate of $`P_zP_{ex}`$ the equilibrium parameters should obey the conditions $`q_{1a}^0=q_{1b}^0`$ and $`q_{2a}^0=q_{2b}^0`$. Therefore, to find the equilibrium configurations we calculate the zeros of four functions $`\frac{E}{q_{kj}}(\mathrm{w}^0)=0`$ with $`k=1,2`$ and $`j=a,b`$ and the parameters restricted as indicated in the above discussion.

We characterize the equilibrium configurations (eqc) by the localization of the center of mass of each component and in Fig.1 we have a graph of the relative distance between the centers of mass as a function of $`\mathrm{\Omega }`$, for appropriately chosen values of $`z_0`$. As shown in Fig.1, for $`z_0=0.23z_{sho}`$ we have only one branch of eqc. For small values of $`\mathrm{\Omega }/\mathrm{}\omega _z`$ the condensates are well separated and when the intensity of the laser field increases the overlap between the condensates increases very slowly up to a value of $`\mathrm{\Omega }`$ when there is a sharp transition to a mixed phase.

When we diminish the value of $`z_0`$, Fig.1, the behavior of eqc changes qualitatively. When $`\mathrm{\Omega }/\mathrm{}\omega _z`$ is small the system behaves as in the previous case. However, when $`\mathrm{\Omega }/\mathrm{}\omega _z`$ increases there is a critical value of $`\mathrm{\Omega }`$, $`\mathrm{\Omega }_{c<}`$, where two branches of eqc appear, one stable, the other unstable. The two stable eqc distinguished by the separation of the centers of mass are called, respectively, distant and near stable eqc. When we further increase the value of $`\mathrm{\Omega }`$, the relative distance of the centers of mass in the unstable eqc increases and merges with the stable distant eqc at a critical value of $`\mathrm{\Omega }`$, $`\mathrm{\Omega }_{c>}`$, such that at $`\mathrm{\Omega }>\mathrm{\Omega }_{c>}`$ one is left with only the near stable eqc, where the condensates are mixed. For smaller values of $`z_0`$ we have the same pattern, the values of $`\mathrm{\Omega }_{c<}`$ where we have three branches of eqc and of $`\mathrm{\Omega }_{c>}`$ where we have the merger of the unstable and distant stable eqc diminishing, this effect being less pronounced for the latter.

Also shown in Fig.1 is the graph of the branches of eqc for $`z_0=0`$. We see that already at $`\mathrm{\Omega }`$ near zero we have three branches of eqc, where now in the near stable eqc there is complete overlap between the condensates ($`\psi _{oa}(z)=\psi _{ob}(z)`$). When $`\mathrm{\Omega }`$ increases, the totally mixed eqc remains, with a density profile independent of $`\mathrm{\Omega }`$, whereas the distant stable and the unstable eqc approach each other and merge at $`\mathrm{\Omega }_{c>}`$, such that, at $`\mathrm{\Omega }>\mathrm{\Omega }_{c>}`$ one is left with only the totally mixed eqc. Our results also show that we have spontaneous symmetry breaking effects at $`z_0=0`$ . Indeed, at $`z_0=0`$, our equations are separately invariant by space reflection, $`P_z`$, and atom exchange $`P_{ex}`$. The totally mixed eqc obey these symmetries separately, whereas the distant stable eqc do not, being invariant only by the product of these transformations.

To find the signature of the onset of criticality, we calculate the normal modes along the branches of eqc. We can group the five normal modes into two sets. In one set there are two normal modes which are a linear combination of an out-of-phase oscillation of the centers of mass of each condensate and an in-phase breathing oscillation of the condensate densities with the center of mass of the mixture and the population fraction at its equilibrium values. In the second set we have three normal modes which are a linear combination of an in-phase oscillation of the centers of mass, out-of-phase breathing oscillation of the condensate densities and particle exchange between the condensates. The splitting of the normal modes into these two groups is a general result since it follows from the invariance of the equations (4-5) under $`P_zP_{ex}`$, the normal modes of the first group being even under this transformation and the one of the second, odd.

We found that the signature of criticality in all cases is the collapse of a normal mode which is an out-of-phase oscillation of center of mass (dipole oscillation) of the condensates and an in-phase breathing oscillation of the condensate densities. In fig.2 we have a graph of the energies of the two normal modes of the first group, which involves the dipole oscillation of the condensates, for $`z_0=0.08z_{sho}`$, which is a case where we found the existence of critical points (see Fig.1).

In Fig.2a and Fig.2b we present a graph of the energies of the two normal modes along the distant and near stable eqc. In the distant eqc we see that there is one normal mode with an almost constant energy, Fig.2a, and other whose energy increases beginning from $`\mathrm{\Omega }=0`$, reaches a maximum value and starts to decrease,Fig.2b. When $`\mathrm{\Omega }`$ approaches the critical value $`\mathrm{\Omega }_{c>}`$, where the distant stable eqc disappears, the energies of the two normal modes approaches each other and at $`\mathrm{\Omega }=\mathrm{\Omega }_{c>}`$ one of the energies goes very abruptly to zero. We have a similar behavior in the graph of the normal mode energies along the near stable eqc, Figs.2a-2b. Approaching from above the point where the near stable eqc disappears, the two energies approaches each other and again one the of them goes abruptly to zero at $`\mathrm{\Omega }=\mathrm{\Omega }_{c<}`$. This behavior is completely general near a critical point.

A “scar” of this critical behavior is also present for a value of $`z_0`$ at the interface between values where we have and we do not have critical points, such as $`z_0=0.23z_{sho}`$ (see Fig.1). We see in Figs.2c-2d that, corresponding to the very narrow range of values where the eqc change from separated to mixed, we have also an abrupt change in the values of the two normal modes energies. In Fig.3a we detached the region of the sharp change and we see that it occurs in a very narrow range of values of $`\mathrm{\Omega }`$ and it is a consequence of a strong level repulsion between the two approaching even normal mode energies.

In Fig.3b we illustrate a generic phenomenon that occurs when $`\mathrm{\Omega }0`$, the appearance of a Goldstone zero energy mode. Indeed, when $`\mathrm{\Omega }0`$, the particle number of each component of the mixture is a conserved quantity and since our theory conserves only the total number of atoms, this violation is translated into the appearance of a zero energy mode. Fig.3b shows how the energy of one of odd normal mode goes to zero for $`z_0=0.23z_{sho}`$.

One question left untouched up to now is the identification of the lowest energy configuration when we have many branches. In our model the answer is that, for small $`\mathrm{\Omega }`$, the distant eqc is always the lowest energy configuration, changing to the near eqc at a higher value of $`\mathrm{\Omega }`$, smaller than $`\mathrm{\Omega }_{c>}`$. However, the energy differences are very small, the equilibrium configurations are almost degenerate.

Our conclusions are based in calculations where we took $`\lambda _{aa}=\lambda _{bb}=\lambda _{ab}`$. It is well known that for homogeneous condensate mixtures, the parameters that controls the mechanism of phase separation is $`\lambda _{aa}\lambda _{bb}\lambda _{ab}^2`$ . In coupled mixtures, we add two additional factors which have opposite effects in the mechanism of spatial separation, the trap displacement and the laser coupling field and a point that deserves investigation is how robust are our conclusions when we relax the equal interaction strength condition.

## IV Summary

To summarize, we have studied the equilibrium and stability properties of a coupled two-component BEC, as function of the laser field strength and trap displacement, using the variational method and the one-dimensional model of reference . The laser field has a stabilizer effect in the mechanism of spatial separation of components in the mixture, opposite to the effect of the trap displacement. We found many branches of eqc, with a host of critical points. In all cases the signature of the onset of criticality is the collapse of a normal mode, which is a linear combination of an out-of-phase translation and an in-phase breathing oscillation of the condensate densities.

When the traps are not displaced, we found eqc which exhibits symmetry breaking effects. In principle these eqc with broken symmetry can be reached by, starting at a sufficiently high value of $`\mathrm{\Omega }`$ and $`z_0`$, adiabatically take the limit $`\mathrm{\Omega }0`$ and $`z_00`$. Taking the limit in the opposite order, we end up in the symmetric eqc (see Fig.1).

Undoubtedly our calculations are simple. However, it unveils a very rich structure in systems of coupled condensates, which should be explored experimentally and theoretically by more complete calculations.

ACKNOWLEDGMENT: This work was partially supported by Fundacção de Amparo $`\stackrel{`}{a}`$ Pesquisa do Estado de S$`\stackrel{~}{a}`$o Paulo (FAPESP) under contract number 00/06649-9. EJVP and MSH are supported in part by CNPq. The work of Chi-Yong Lin and Da-Shin Lee was supported in part by the National Science Council, ROC under the Grant NSC-90-2112-M-259-010-Y.EJVP,MSH and AFRTP are members of the INFOQUANTICA contract 62.0057/01-7-PADCTIII/MIL.


text:
# Dielectric behavior of oblate spheroidal particles: Application to erythrocytes suspensions

## I Introduction

Dielectrophoresis and electrorotation (ER) offer a unique capability of monitoring the dielectric properties of dispersions of colloids and biological cells. Under the action of external fields, these particles exhibit rich fluid-dynamic behaviors as well as various dielectric responses. It is of importance to investigate their frequency-dependent responses to ac electric fields, which yields valuable information on the structural (Maxwell-Wagner) polarization effects . The polarization is characterized by a variety of characteristic frequency-dependent changes known as the dielectric dispersion. In particular, for the $`\beta `$-dispersion (also called the Maxwell-Wagner dispersion, ranging from KHz to MHz), it may suffice to focus on the induced dipole moments of the particles. In this work, we will analyze the $`\beta `$-dispersion based on the spectral representation theory.

In the last two decades, various experimental tools have been developed to analyze the polarization of biological cells - dielectric spectroscopy , dielectrophoresis  and ER  techniques. Among these techniques, conventional dielectrophoresis and ER are usually applied to analyze the frequency dependence of translations and rotations of individual cells in an inhomogeneous and rotating electric field, respectively . Moreover, one is able to monitor the cell movements by using automated video analysis  as well as light scattering methods . As a matter of fact, the general cause of ER is the existence of a phase difference between the field-induced dipole moment and the external rotating field, resulting in a desired torque which causes the cells to rotate.

In the dilute limit, the ER of individual cell can be predicted by ignoring the mutual interaction between the cells. However, the cells may aggregate under the influence of the external field. In this case the Brownian motion can be neglected, and the cell system becomes non-dilute even though it is initially dilute. As an initial model, we have recently studied the ER of two approaching spherical particles in the presence of a rotating electric field . We showed that when the two particles approach and finally touch, the mutual polarization interaction between the particles leads to a change in the dipole moment of individual particles and hence the ER spectrum, as compared to that of isolated particles, via the multiple image method .

In this work, we consider further the effects of cell shape on the ER spectrum. In fact, a cell may readily deviate from a perfect spherical shape due to many reasons, e.g., under a hydrostatic pressure . In addition to a prolate spheroidal shape , there exist cells of oblate spheroidal shape, such as human erythrocytes.

Regarding the spectral representation approach , it is a rigorous mathematical formalism of the effective dielectric constant of a two-phase composite material. Also, the approach was extended to deal with a three-phase material . The spectral approach offers the advantage of the separation of material parameters (namely the dielectric constant and conductivity) from the cell structure information, thus simplifying the study. From this approach, one can readily derive the dielectric dispersion spectrum, with the dispersion strength as well as the characteristic frequency being explicitly expressed in terms of the structure parameters and the materials parameters of the cell suspension . In the present work, the analytic expression for the characteristic frequency is derived with the aid of the spectral representation approach, thus simplifying the fitting of experimental data. We show good agreement with experimental data on human erythrocytes .

## II FORMALISM

Consider a suspension of spheroidal particles of complex dielectric constant $`\stackrel{~}{ϵ}_1`$ coated with a shell of $`\stackrel{~}{ϵ}_s`$ dispersed in a host medium of $`\stackrel{~}{ϵ}_2`$, where

$$\stackrel{~}{ϵ}=ϵ+\sigma /i2\pi f$$
(1)
with $`f`$ being the frequency of the applied field $`𝐄_0=E_0\widehat{𝐳}`$. The depolarization factors of the spheroidal particles is described by a sum rule

$$L_z+2L_{xy}=1$$
(2)
where $`L_z(1/3<L_z<1)`$ and $`L_{xy}`$($`=L_x=L_y`$) are the depolarization factors along the $`z`$\- and $`x`$\- (or $`y`$-) axes of the oblate spheroidal particle, respectively. $`L_z=L_{xy}=1/3`$ just indicates sphere, while $`L_z=1`$ disk.

The phenomenon of ER is based on the interaction between a rotating electric field $`𝐄`$ and the induced dipole moment $`𝐌`$. The dipole moment of the particle arises from the induced charges that accumulate at the interface of the particle. In circularly polarized fields, the axis of lowest polarizability, namely $`z`$axis of the present particle, should be oriented perpendicular to the plane of field rotation .

The angle between $`𝐌`$ and $`𝐄`$ is denoted by $`\theta `$, where $`\theta =\omega \times time`$ and $`\omega =2\pi f`$ is angular velocity of the rotating electric field. The torque acting on the particle is given by the vector cross product between the electric field and the dipole moment, so that only the imaginary part of the dipole moment contributes to the ER response. In the steady state, the frequency-dependent rotation speed $`\mathrm{\Omega }(f)`$, which results from the balance between the torque and the viscous drag, is given by

$`\mathrm{\Omega }(f)`$ $`=`$ $`F(ϵ_2,E,\eta )Im[\stackrel{~}{b}_x\mathrm{cos}^2\theta +\stackrel{~}{b}_y\mathrm{sin}^2\theta ],`$ (3)
$`=`$ $`F(ϵ_2,E,\eta )Im[\stackrel{~}{b}_x]=F(ϵ_2,E,\eta )Im[\stackrel{~}{b}_y]F(ϵ_2,E,\eta )Im[\stackrel{~}{b}_{xy}],`$ (4)
where $`Im[\mathrm{}]`$ indicates taking the imaginary parts of $`[\mathrm{}]`$, the angular brackets denote a time average, and $`F`$ is a coefficient which is proportional to the square of magnitude of field $`E^2`$, but being inversely proportional to the dynamic viscosity $`\eta `$ of the host medium. For spherical cells, $`F=ϵ_2E^2/2\eta `$. Regarding $`\eta `$, spin friction expression suffices for spherical particle. However, for spheroidal particles or many particles interacting in a suspension, we must consider the more complicated suspension hydrodynamics. Since the angular velocity of the rotating field is much greater than the electrorotation angular velocity, i.e., $`\omega \mathrm{\Omega }`$, the time averages are just equal to 1/2. For a single coated spheroidal particle, the dipole factor $`\stackrel{~}{b}_z`$ is given by

$$\stackrel{~}{b}_{xy}=\frac{1}{3}\frac{(\stackrel{~}{ϵ}_s\stackrel{~}{ϵ}_2)[\stackrel{~}{ϵ}_s+L_{xy}(\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_s)]+(\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_s)y[\stackrel{~}{ϵ}_s+L_{xy}(\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_s)]}{(\stackrel{~}{ϵ}_s\stackrel{~}{ϵ}_1)(\stackrel{~}{ϵ}_2\stackrel{~}{ϵ}_s)yL_{xy}(1L_{xy})+[\stackrel{~}{ϵ}_s+(\stackrel{~}{ϵ}_1\stackrel{~}{ϵ}_s)L_{xy}][\stackrel{~}{ϵ}_2+(\stackrel{~}{ϵ}_s\stackrel{~}{ϵ}_2)L_{xy}]}$$
(5)
where $`y`$ is the volume ratio of the core to the whole coated spheroid.

We are now in a position to represent $`b_{xy}`$ in the spectral representation. In what follows, we will show that from the spectral representation, we can obtain the analytic expression for the characteristic frequency at which the maximum ER velocity occurs. Let $`\stackrel{~}{ϵ}_1=\stackrel{~}{ϵ}_2(11/\stackrel{~}{s})`$, and assume $`x=\stackrel{~}{ϵ}_s/\stackrel{~}{ϵ}_2`$, we obtain

$$\stackrel{~}{b}_{xy}=NP+\frac{F_1}{\stackrel{~}{s}s_1}$$
(6)
where $`NP`$ denotes the nonresonant part  which vanishes in the limit of unshelled spheroidal cells. In Eq.(6), the various quantities are given by

$`s_1`$ $`=`$ $`{\displaystyle \frac{L_{xy}\{1(1x)[1(1L_{xy})(1y)]\}}{x+(1L_{xy})L_{xy}(1x)^2(1y)}},`$ (7)
$`F_1`$ $`=`$ $`{\displaystyle \frac{x^2y}{3[x+(1L_{xy})L_{xy}(1x)^2(1y)]^2}},`$ (8)
$`NP`$ $`=`$ $`{\displaystyle \frac{[L_{xy}(1x)+x](1x)(1y)}{3[x(1L_{xy})L_{xy}(1x)^2(1y)]}}.`$ (9)
Note that we have assumed $`x`$ to be a real number, that is, $`x\sigma _s/\sigma _2`$, which will be justified below. After substituting $`\stackrel{~}{ϵ}=ϵ+\sigma /i2\pi f`$ into Eq.(6), we rewrite $`\stackrel{~}{b}_{xy}`$ after simple manipulations

$$\stackrel{~}{b}_{xy}=(NP+\frac{F_1}{ss_1})+\frac{\delta ϵ}{1+if/f_c}$$
(10)
with $`s=(1ϵ_1/ϵ_2)^1`$ and $`t=(1\sigma _1/\sigma _2)^1`$, where the dispersion magnitude and characteristic frequency admit respectively

$$\delta ϵ=F_1\frac{st}{(ts_1)(ss_1)},$$
(11)
$$f_c=\frac{1}{2\pi }\frac{\sigma _2}{ϵ_2}\frac{s(ts_1)}{t(ss_1)}.$$
(12)
## III Numerical result

In Fig.1(a)-(c), the dependence of the pole, characteristic frequency and dispersion magnitude on the parameter $`x`$ are investigated for different particle shape. For clarity, we set $`z=1/y`$ from now on. It is evident that both $`s_1`$ and $`f_c`$ decreases monotonically as $`x`$ increases; larger depolarization factor leads to smaller $`s_1`$ and $`f_c`$. However, for increasing $`x`$, $`\delta ϵ`$ increases concomitantly; also, $`\delta ϵ`$ increases while $`L_z`$ increases.

In Fig.1(d)-(f), we investigate the dependence of the pole, characteristic frequency and dispersion magnitude on $`x`$ for different shell thickness. For $`z=1`$, namely without shell, $`s_1=0.3(=L_{xy})`$ always. Increasing $`x`$ leads to decreasing $`s_1`$ or $`f_c`$. It is evident that thicker shell yields larger pole or $`f_c`$ as $`x<1`$, but smaller as $`x>1`$. In the dispersion magnitude plot, $`\delta ϵ`$ is constant for the uncoated particles. In the case of coated particles, increasing $`x`$ yields increasing $`\delta ϵ`$; moreover, larger $`z`$, smaller $`\delta ϵ`$.

In Fig.2, the dependence of material parameters (namely, $`s`$, $`t`$, $`ϵ_2`$, and $`\sigma _2`$) on $`f_c`$ and $`\delta ϵ`$ is discussed. $`f_c`$ is weakly dependent on $`s`$, but strongly dependent on $`t`$, $`ϵ_2`$ and $`\sigma _2`$. Smaller $`|t|`$ (or $`ϵ_2`$), larger $`f_c`$. However, larger $`\sigma _2`$, larger $`f_c`$. On the other hand, $`\delta ϵ`$ is weakly dependent on $`s`$ and $`t`$. Also, it is evident that $`ϵ_2`$ and $`\sigma _2`$ play no role in $`\delta ϵ`$, as Eq.(11) actually shows. Note all curves in Fig.2(g) and (h) are respectively overlapped. In other words, this predicts that $`ϵ_2`$ and $`\sigma _2`$ plays no role in the peak value of the angular velocity (see Fig.4).

In Fig.3, the dependence of $`L_z`$, $`z`$ and $`x`$ on $`Im[b_{xy}]`$ are investigated, respectively. There is always a peak, the location of which is just the characteristic frequency $`f_c`$. As $`L_z`$ (or $`x`$) increases, $`f_c`$ becomes red-shift, while the corresponding peak value increases as well. However, a larger $`z`$ yields smaller peak value, while $`f_c`$ is changed weakly.

In Fig.4, we discuss the dependence of $`s`$, $`t`$, $`ϵ_2`$ and $`\sigma _2`$ on $`Im[b_{xy}]`$. The $`s`$ effect is small so that it can be neglected. Smaller $`|t|`$ leads to larger $`f_c`$; at the same time, increasing peak value appears. Increasing $`ϵ_2`$ yields decreasing $`f_c`$, and the related peak value is unchanged. For increasing $`\sigma _2`$, $`f_c`$ increases concomitantly, but the peak value remains unchanged as well.

In Fig.5, the experimental data is extracted from an experiment on human erythrocytes . Obviously, there is a good agreement between theory and experiment. During this fitting, we model the cell as oblate spheroid ($`L_z=0.65`$).

## IV Discussion and conclusion

In the present work, we have discussed the effect of oblate spheroidal particle shape on the ER spectrum. In reality, there exist many cells in the form of oblate spheroid, such as human erythrocytes . However, there exist only a few theories to discuss the ER spectrum of such cells (e.g.,  and references therein). Our theory is advantageous in that the characteristic frequency at which the maximum rotational velocity occurs, is derived analytically, which simplifies the fitting of experimental data.

We have assumed $`x`$, namely the ratio of the shell to host dielectric constant, to be a real and positive number. In fact, its imaginary part is indeed small. If we had retained the imaginary part of $`x`$, two peaks would appear, and the conductivity-dominated peak would have occurred at substantially lower frequency. Therefore, the neglect of the imaginary part of $`x`$ is to drop the lower frequency peak. Moreover, according to our calculations, there are two (degenerate) sub-dominant poles associated with $`b_{xy}`$ in the spectral representation. Thus, in the present work, only one capacitance-dominated peak has been shown.

In a recent work , we developed simple equations to describe the ER of particles in a suspension from the spectral representation, and obtained a good fitting on a coated-bead ER assay by using the coated spherical model. In this work, an extension to oblate spheroidal case has been made. Our theory serves as a basis which describes the parameter dependence of the polarization and thereby enhances the applicability of various cell models for the analysis of the polarization mechanisms. In this connection, the shell-spheroidal cell model may readily be extended to multi-shell cell model. However, we have shown that the multi-shell nature of the cell may have a minor effect on the ER spectrum .

We have considered the isolated cell case, which is a valid assumption for low concentration of cells. However, for a higher concentration of cells, we should consider the mutual interaction between cells. When the volume fraction of the suspension becomes large, the particles may aggregate in the plane of the rotating applied field, and the mutual interactions between the suspended particles can modify the spin friction, which is a key to determine the angular velocity of ER. For ER of two particles, we have successfully applied the spectral representation to deal with the dispersion frequency . However, the determination of the spin friction is still lacking. For two particles, the basic tool is the reflection method , being analogous to the multiple image method in electrostatics, but being only valid for two particles. For more than two particles, we need a first-principles method, e.g., the Green’s function (Oseen tensor) formulation. For a dilute suspension, however, one may adopt the effective medium theories  to capture the effective viscosity of a suspension, hence modifying further the ER spectrum.

One may argue that normal erythrocytes may deviate from an oblate spheroidal shape. For cells of non-conventional shapes, one may use a first-principles approach , recently developed to deal with particles of rod shape. As shown in Ref., the derivation from spheroidal model can be small when the cells are rotating with their long axes along the applied field.

In summary, we have presented a theoretical study of electrorotation based on the spectral representation theory. From the theoretical analysis, we find that cell shape as well as the coating can change the characteristic frequency. By adjusting the cell shape, dielectric properties and the thickness of the coating, it is possible to obtain good agreement between our theoretical predictions and the experimental data on human erythrocytes suspensions.

## Acknowledgments

This work was supported by the Research Grants Council of the Hong Kong SAR Government under grant CUHK 4245/01P. K.W.Y. acknowledges useful conversation with Professor G. Q. Gu.

FIG.1

FIG.2

FIG.3

FIG.4

FIG.5


text:
# Neutron star properties with in-medium vector mesons

## 1 Introduction

The exploration of the properties of in-medium hadrons is a topic of very great current interest. This is reflected by the number of experiments either planned or recently performed which target at the determination of medium modification of hadrons, especially of $`K^\pm `$ and $`\rho `$ mesons. The first results of the CERES and HELIOS collaborations show a significant enhancement of the production of low-mass lepton pairs  in relativistic nucleus-nucleus collisions, compared to the yield that one expects from proton-proton collisions. Since these lepton pairs are direct signals from the decay of vector mesons within the hot reaction zone, this enhancement points toward modifications of the masses and/or widths of vector ($`\rho `$ and $`\omega `$) mesons in a dense hadronic environment .

Medium modifications of mesons, as described just above, are generally ignored in standard treatments of hot and dense hadronic matter . In the framework of such treatments, the nuclear interaction is described through the exchange of mesons whose properties are entirely medium independent. The most important of these mesons are the $`\pi `$, $`\sigma `$, $`\rho `$ and $`\omega `$ mesons. It seems to be a general feature that standard treatments fail to reproduce the observed enhancement in the production of low-mass lepton pairs . This appears to be different for theoretical treatments which account for in-medium modifications of vector mesons.

In this paper we introduce two such treatments, one based on the Brown-Rho scaling scheme while the other constitutes a novel approach based on the techniques of effective nuclear field theory, which are used to computed corrections to the nuclear equation of state (EOS) that originate from the in-medium modifications of $`\rho `$ and $`\omega `$ vector meson properties. These models for the equation of state are then used to determine several key properties of neutron stars, such a the limiting stellar mass, the mass-radius relationship, and the Keplerian velocity. These quantities are known to be very sensitive to variations in the nuclear equation of state  and, therefore, serve to test as to whether the in-medium vector-meson modifications computed for the two theoretical treatments are compatible with solid astrophysical data on neutron stars.

## 2 Theoretical treatment of vector mesons in matter

### 2.1 Choice of Lagrangian

In order to study the medium modifications on vector mesons, described in section 1, we assume that the interactions producing these modifications are independent of the interactions used in the Walecka lagrangian, so we may take them into account by only modifying the meson properties themselves. For the nuclear interaction we adopt the non-linear Walecka model. In this model the basic degrees of freedom are the nucleon $`N`$, the scalar self-interacting $`\sigma `$ meson, and the $`\omega `$ and $`\rho `$ vector mesons. The lagrangian is given by

$``$ $`=`$ $`\overline{\mathrm{\Psi }}\gamma ^\mu (i_\mu g_\omega \omega _\mu )\mathrm{\Psi }\overline{\mathrm{\Psi }}(m_Ng_\sigma \sigma )\mathrm{\Psi }`$ (1)
$`{\displaystyle \frac{1}{4}}(^\nu \omega ^\mu ^\mu \omega ^\nu )(_\nu \omega _\mu _\mu \omega _\nu )+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _\mu \omega ^\mu `$
$`+{\displaystyle \frac{1}{2}}(_\nu \sigma ^\nu \sigma m_\sigma ^2\sigma ^2){\displaystyle \frac{1}{3}}bm_N(g_\sigma \sigma )^3{\displaystyle \frac{1}{4}}c(g_\sigma \sigma )^4`$
$`{\displaystyle \frac{1}{4}}(^\nu \rho ^\mu ^\mu \rho ^\nu )(_\nu \rho _\mu _\mu \rho _\nu )+{\displaystyle \frac{1}{2}}m_\rho ^2\rho _\mu \rho ^\mu `$
$`g_\rho \rho _\nu ({\displaystyle \frac{1}{2}}\overline{\mathrm{\Psi }}\gamma ^\nu \tau \mathrm{\Psi }+\rho _\mu \times (^\nu \rho ^\mu ^\mu \rho ^\nu )+2g_\rho (\rho ^\nu \times \rho ^\mu )\times \rho _\mu ),`$
where $`\mathrm{\Psi }`$ denotes the nucleon fields, and $`g_\sigma `$, $`g_\omega `$, $`g_\rho `$, $`b`$, and $`c`$ are the meson-nucleon coupling constants . The equations of motion that follow from (1) are solved in the mean-field approximation, where the meson fields are substituted by their classical expectation values. The parameters of the model ($`g_\sigma `$, $`g_\omega `$, $`g_\rho `$, $`b`$, and $`c`$) are fitted to the binding energy ($`16`$ MeV), nuclear incompressibility ($`K=240`$ MeV) and effective nucleon mass ($`m_N^{}=0.8m_N`$) at the saturation density ($`\rho _0=0.168\mathrm{fm}^3`$) of infinite nuclear matter. (A soft equation of state was predicted recently by a novel many-body approach that is based on the principles of chiral perturbation theory .) In passing we mention that we ignore possible medium modifications of the properties of the scalar $`\sigma `$ meson. Being a Goldstone boson, it is not clear how the properties of this meson may change in dense hadronic matter. Application of the Brown-Rho scaling scheme to the $`\sigma `$-meson mass, for instance, leads to an over-binding of nuclear matter . Since we are mainly interested in the high-density part of the equation of state of hadronic matter, which is dominated by the exchange of vector mesons among the nucleons, it is rather save for our purposes to ignore medium effects on the $`\sigma `$ meson.

### 2.2 In-medium modifications of vector meson properties

Next we turn to the two alternative models employed in this paper to explore the in-medium properties of $`\rho `$ and $`\omega `$ mesons in dense hadronic matter. The first model is based on the familiar Brown-Rho scaling  according to which the masses of vector mesons and nucleons in a hadronic medium scale as

$$\frac{m_\omega ^{}(\rho )}{m_\omega }=\frac{m_\rho ^{}(\rho )}{m_\rho }=\frac{m_N^{}(\rho )}{m_N},$$
(2)
where the asterisks denotes the particle masses in the medium. For the motivation of this scaling law we refer to the original literature . The effective mass of the nucleon at a given density, $`m_N^{}(\rho )`$, is taken from the Walecka model,

$$m_N^{}(\rho )=m_Ng_\sigma <\sigma >,$$
(3)
with $`<\sigma >`$ the mean-field ground-state expectation value of the $`\sigma `$ meson. The second model which we study here  is a novel approach based on the techniques of effective nuclear field theory. In the framework of this model, a low-density expansion is performed in order to establish a connection between the meson-nucleon scattering amplitudes and the properties of mesons in nuclear matter . This scheme leads to a systematic expansion for the in-medium self energy of a meson (as well as any other particle) in terms of the corresponding vacuum scattering amplitudes. The vector-meson–nucleon scattering amplitudes are not directly related to experimental observables. Instead they are determined indirectly in a coupled-channel scheme where the $`\rho N`$ and $`\omega N`$ channels enter in intermediate and final states of measured processes. The parameters of this theory are fixed by fitting the available meson-nucleon scattering data in the relevant energy range. We employ these scattering amplitudes to construct the vector-meson self energies in nuclear matter to leading order in density.

The effective field theory provides us with the scattering lengths and spectral functions of vector mesons as a function of the density. For $`\rho `$ and $`\omega `$ vector mesons, the $`\rho `$$`N`$ and $`\omega `$$`N`$ scattering lengths are defined as ($`V=\rho ,\omega `$)

$$a_{VN}=f_{VN}(\sqrt{s}=m_N+m_V),$$
(4)
which leads to

$$a_{\rho N}=(0.1+0.6\mathrm{i})\mathrm{fm}\mathrm{and}a_{\omega N}=(0.5+0.2\mathrm{i})\mathrm{fm}.$$
(5)
To lowest order in density, this corresponds to the following in-medium modifications of the vector meson masses and their widths at nuclear matter density:

$$\mathrm{\Delta }m_\rho 10\mathrm{MeV},\mathrm{\Delta }m_\omega 50\mathrm{MeV},\mathrm{\Delta }\mathrm{\Gamma }_\rho 120\mathrm{MeV},\mathrm{\Delta }\mathrm{\Gamma }_\omega 40\mathrm{MeV}.$$
(6)
We stress that the coupling of the vector mesons to baryon resonances below threshold, which is reflected in the strong energy dependence of the scattering amplitudes, cannot be neglected. One should therefore use the whole vector-meson spectral functions as done in this paper in order to achieve a most reliable description. As it turns out, however, application of the simple mass shift scheme for the in-medium effects on vector mesons leads to results that are rather similar to those obtained by using the full treatment, based on medium dependent spectral functions.

As a side-remark we mention that the validity of these two schemes adopted to include in-medium modifications in the equation of state is naturally limited to moderately compressed hadronic matter. The extent to which they are applicable to very highly compressed hadronic matter as existing in the centers of very heavy neutron stars is an open issue. The situation looks quite promising, though, for neutron stars of average mass, $`M1.4M_{}`$, whose central densities may be just a few times higher than the density of ordinary nuclear matter, depending of the stiffness of the equation of state (cf. section 4).

## 3 Numerical outcome

We begin our discussion of the impact of in-medium modifications of vector mesons on the nuclear equation of state and neutron star properties with fig. 1 which shows the spectral functions of $`\rho `$ and $`\omega `$ mesons for three selected sample densities.

These results are complementary to the recent calculations of ref. . Taking into account the photon induced reactions via a generalized vector meson dominance the latter paper finds weaker medium effects for the $`\rho `$-meson than established here. In the next step we apply the two schemes introduced in section 2.2, i.e. Brown-Rho scaling and modified spectral function method, to determine the

medium induced modifications of vector mesons and calculate the effective nucleon mass in matter, $`M^{}`$, as well as the nuclear equation of state over a considerable range of densities. The outcome is shown in figs. 2 and 3, respectively. It is evident from fig. 2 that in-medium effects do not alter the overall density dependence of $`M^{}`$ in dense hadronic matter: the effective nucleon mass drops monotonically with increasing density, as known from standard calculations . At the density of ordinary nuclear matter, $`\rho _0`$, $`M^{}`$ is basically unchanged from the value where in-medium effects on the vector mesons are ignored. Deviations however begin to show up at densities that are three times $`\rho _0`$ or higher and grow monotonically with density. This feature, together with the pronounced stiffening of the equation of state caused by the in-medium effects, shown in fig. 3, renders such modifications very important for neutron star properties. We will come back to this issue in greater detail in section 4. The stiffening of the equation of state is a consequence of the fact that both energy density and pressure depend on the vector meson masses as $`m^2`$. This proportionality characterizes the high-density regime of the equation of state, which is dominated by the exchange of vector mesons. Quantitatively, the in-medium vector-meson masses $`<m_V^2>`$ were computed self-consistently, with the weights entering in the expectation values given by the spectral functions.

## 4 Neutron star structure

In this section we explore the implications of in-medium modifications of vector mesons for the structure of neutron stars. The properties of such objects provide important constraints on the behavior of dense matter, for any model for the nuclear equation of state which fails to accommodate neutron star masses of at least the mass of the Hulse-Taylor radio pulsar PSR 1913+16, $`(1.444\pm 0.003)M_{}`$, and the rotational periods of the two fastest presently known neutron stars, 1937+21 and 1957+20, each one rotating at 1.6 ms (620 Hz), can definitely be ruled out as a viable candidate for the true equation of state .

From model calculations, it is known that neutron stars possess rather complex interior structures . Only in the most primitive conception, a neutron star is constituted from neutrons. At a more accurate representation, a neutron star contains neutrons ($`n`$) and a certain number of protons ($`p`$) in chemical equilibrium,

$$np+e^{}+\overline{\nu },$$
(7)
which implies for the corresponding chemical potentials that

$$\mu _n=\mu _p+\mu _e+\mu _{\overline{\nu }}.$$
(8)
Since neutrinos do not accumulate in neutron stars older than a few seconds, we set $`\mu _{\overline{\nu }}=0`$. The electric charges of the protons in the star must be balanced to very high precision by electrons ($`e^{}`$), which leads to the additional constraint for the densities of protons and electrons

$$\rho _p=\rho _e.$$
(9)
If the particles in the center of a neutron star would arrange themselves in a way other than dictated by eq. (9), gravity would not be able to bind a neutron star.

Figure 4 shows the equation of state of neutron star matter computed from the lagrangian of eq. (1) subject to the two neutron star matter constraints (8) and (9). As for the nuclear matter case shown in fig. 3, the free-particle case clearly provides the least pressure of all three models. Yet the equation of state is stiff enough to accommodate neutron stars as heavy as about $`2M_{}`$, as will be discussed in more detail below. Another issue concerns the relative proton fraction in neutron star matter, shown in the upper-right panel of fig. 4, which is of central importance for the cooling of neutron stars. As pointed out in , if $`\rho _p/\rho `$ becomes larger than the critical fraction of about 11%, the so-called direct Urca process,

$$np+e+\nu _e,$$
(10)
where neutrons can transform to protons and electrons without the need of a bystander particle, becomes possible in neutron stars. Otherwise neutron stars were to cool via the less efficient standard Urca process,

$$n+nn+p+e+\nu _e.$$
(11)
The direct Urca process is crucial for the temperature evolution of neutron stars as it speeds up the cooling of such objects considerably in comparison with the standard Urca process. Our models for the equation of state predict that the process in (10) can take place at densities around $`2.3\rho _0`$ if in-medium effects on vector mesons are ignored. Accounting for them, on the other hand, brings this threshold down to densities as low as $`1.5\rho _0`$. Because of the rather flat density profiles in the inner parts of neutron stars, this implies that the cores in neutron stars which exhibit the direct Urca process are several kilometers bigger for an equation of state accounting for in-medium effects. As an example of a neutron star that may require very fast cooling as driven by the direct Urca process (10) we quote pulsar PSR 1929+10 .

The stiffness/softness of the equation of state above nuclear matter density intimately manifests itself in the interplay between radius, mass, and central density of a neutron star, as discussed next. We begin with the mass-radius relationship of non-rotating neutron stars computed by solving the Tolman-Oppenheimer-Volkoff equations, which describe the properties of spherically symmetric stars that are in general relativistic, hydrostatic equilibrium . The outcome is shown in fig. 5, which reveals that our collection of equation of states predicts limiting gravitational neutron star masses between about $`2M_{}`$ and $`3M_{}`$, with the high-mass end obtained for the equations of state that include in-medium effects. This range accommodates even the heaviest neutron stars that may exist in stellar binary systems. Examples of which are the X-ray pulsar Vela X-1, whose mass is $`M=1.87_{0.17}^{+0.23}M_{}`$ , and the burster Cygnus X-2 whose mass is $`M=(1.8\pm 0.4)M_{}`$ . Indications for the possible existence of very heavy neutron stars, with masses around $`2M_{}`$, may also come from the observation of quasi-periodic oscillations in luminosity in low-mass X-ray binaries . These mass values are larger than the typical $`(1.36\pm 0.08)M_{}`$ masses found in neutron star binaries, presumably due to accreted matter. It is evident that neutron stars that heavy, provided the mass determinations turn out to be robust, would require extremely stiff equations of state. The inclusion of in-medium effects of vector mesons would alter the equation of state in the right direction.

Another very important conclusion follows from the general feature that stiff equations of state make neutron stars less dense . This effect can be quite dramatic, depending of the stiffness of the equation of state. As shown in fig. 6, the models computed in this paper predict central densities on neutron stars of canonical mass of around $`3\rho _0`$ if no in-medium effects are taken into account. This value drops down to $`1.5`$ to $`2\rho _0`$ if in-medium effects are taken into account. Therefore, if the medium effects as computed in this paper have their correspondence in the full treatment of dense hadronic matter, then the consequences for the composition of neutron stars of masses $`M1.5M_{}`$ would be quite dramatic, for it leaves basically no room for the existence of novel degrees of freedom, such as quark matter or boson condensates, in the centers of such objects. Only neutron stars with masses close to their limiting masses (solid dots in fig. 6) may possibly hide these novel phases of matter.

Rapid neutron star rotation constitutes another important constraint on the stiffness/softness of the nuclear equation of state. The two most rapidly rotating neutron stars presently known, PSR 1937+21 and PSR 1957+20, rotate at 1.58 ms, which corresponds to about 620 rotations per second. Such high rotation rates constitute a considerable fraction of the Kepler frequency, $`\mathrm{\Omega }_\mathrm{K}`$, at which mass shedding from the star’s equator sets it. Evidently, mass shedding sets an absolute limit on rapid rotation which cannot be overcome by any stably rotating object. In oder to determine $`\mathrm{\Omega }_\mathrm{K}`$ for a given model for the equation of state one has to go beyond the Tolman-Oppenheimer-Volkoff treatment and solve Einstein’s field equations

$$R^{\mu \nu }\frac{1}{2}g^{\mu \nu }R=8\pi T^{\mu \nu }(\rho ,P(\rho )),$$
(12)
self-consistently in combination with the general relativistic expression describing the onset of mass-shedding at the star’s equator ,

$$\mathrm{\Omega }_\mathrm{K}=\omega +\frac{\omega ^{}}{2\psi ^{}}+e^{\nu \psi }\sqrt{\frac{\nu ^{}}{\psi ^{}}+\left(\frac{\omega ^{}}{2\psi ^{}}e^{\psi \nu }\right)^2}.$$
(13)
The quantities $`R^{\mu \nu }`$, $`g^{\mu \nu }`$, and $`R`$ denote the Ricci tensor, metric tensor, and Ricci scalar (scalar curvature), respectively. The dependence of the energy-momentum tensor $`T^{\mu \nu }`$ on the equation of state is indicated in eq. (12). The quantities $`\omega `$, $`\nu `$, and $`\psi `$ in eq. (13) denote the frame dragging frequency of the local inertial frames, and the time and space-like metric functions, respectively. The primes denote derivatives with respect to Schwarzschild radial coordinate, and all functions on the right are evaluated at the star’s equator. All the quantities on the right hand side of eq. (13) depend also on $`\mathrm{\Omega }_\mathrm{K}`$, so that it is not an equation for $`\mathrm{\Omega }_\mathrm{K}`$, but a transcendental relationship which the solution of the equations of stellar structure, resulting from eq. (12), must satisfy if the star is rotating at its Kepler frequency.

The numerical outcome for $`\mathrm{\Omega }_\mathrm{K}`$ is shown in fig. 7, where we plot $`\nu _\mathrm{K}`$ ($`2\pi /\mathrm{\Omega }_\mathrm{K}`$) as a function of rotational stellar mass, $`M`$. For neutron stars of canonical mass, $`M1.4M_{}`$, one reads off from this figure that such stars can perform up to $`1000`$ rotations per second before breaking up into pieces. This comfortably accommodates the spin frequencies of the two fastest rotating neutron stars presently known. Rapid rotation in combination with stellar masses provides a double constraint on the equation of state, since the successful model for the nuclear equation of state must account for the smallest observed rotational neutron star periods as well as the heaviest neutron star masses. The constraints on the equation of state become the more stringent the smaller a neutron star’s rotational period and the larger its mass. Searches for such neutron stars are being performed at several radio observatories. It is evident from fig. 7 that only very compact neutron stars, compressed to small radius values, will be able to withstand rapid rotation rates that are significantly higher than 600 Hz. Neutron stars at or close to the limiting mass of a given stellar sequence are such candidates. Depending on stiffness/softness of the equation of state, these may rotate at frequencies as large as 1500 to 1700 Hz, which corresponds to rotational periods between 0.7 and 0.4 milliseconds, respectively. The only other class of stars understood to withstand such rapid rotation are the strange stars, which should replace neutron stars if 3-flavor strange quark matter is the true ground-state of the strong interaction , which, despite of many years of research, is still an open issue. As a final point, we point out that rapid rotation stabilizes a neutron star against gravitational collapse. Figure 6 reveals that the mass that can be supported due to rapid rotation may be up to $`20`$% higher than non-rotating mass value. It is striking that the rotation-induced mass increase is accompanied by a reduction in central star density, contrary to what one would expect from simple mass loading onto a star. The explanation is provided by the centrifugal pressure that builds up inside a rotating star. For rapidly spinning neutron stars, this contribution is quite significant so that considerably less interior pressure (typically 30 to 40%, as can read off from fig. 6 in combination with fig. 4) must be provided by the equation of state itself. As already mentioned in connection with the discussion of the properties of non-rotating neutron stars, lower densities at the center of a neutron star rival with exotic particle processes predicted for high-density matter, like the formation of boson condensation and the transition into quark matter.

## 5 Conclusions

This paper presents an investigation of the effect of in-medium modification of vector mesons on the nuclear equation of state and the properties of neutron stars. The latter serve to test the compatibility of in-medium effects with observed data. We find that in-medium modifications reduce effectively the masses of vector mesons which in turn stiffens the equation of state considerably. This has several intruigung implications for the structure and composition of neutron stars. Firstly, because of the stiffening of the equation of state, even the most massive neutron stars, like Vela X-1 ($`M=1.87_{0.17}^{+0.23}M_{}`$) and the burst source Cygnus X-2 ($`M=(1.8\pm 0.4)M_{}`$), can be easily supported by equations of state accounting for medium-modified vector meson properties. Secondly, knowledge of the limiting neutron star mass is key in order to identify black hole candidates. That is, if the mass of a compact companion of an optical star is determined to exceed the limiting mass of a neutron star, it must be a black hole. The limiting mass of stable neutron stars in our theory is between 2.6 and $`2.8M_{}`$. Hence any binary heavier than these values would be low-mass black holes. The third striking point is that the computed neutron star models have rather low central densities, just a few times the density of ordinary nuclei for canonical neutron star masses around $`M1.5M_{}`$. If this should withstand more elaborate future treatments, it would leave basically no room for the existence of novel degrees of freedom such as quark matter or boson condensates in the centers of neutron stars. Such objects would then merely consist of chemically equilibrated neutrons and protons.

This work were partly done while one of the author (Gy. Wolf) stayed in JAERI as a STA fellow, and supported by Hungarian Research Foundation (OTKA) grants T 26543 and T 30855 grant.

## References