text:
# Bound State Calculations for Three Atoms Without Explicit Partial Wave Decomposition

## I Introduction

In recent years the 4He trimer has been the center of several theoretical investigations (see, for example, Refs. Fedorov ; Roudnev ; Motoold  and references therein). From all methods employed in these studies, the Faddeev equations method is perhaps the most attractive since it reduces the Schrödinger equation for three particles into a system of integral or differential equations which can be used to study bound states and scattering processes in a rigorous way.

The differential form of the Faddeev equations has been proposed long ago by Noyes and Fiedeldey Noyes . Since then, these equations have been used in bound states calculations in Nuclear and Coulomb systems GignLav ; MerkGignLav ; Friar ; KviHu ; YkFl  – just to mention a few references. The derivation and discussion by Merkuriev of the asymptotic boundary condition of the differential Faddeev equations (DFE) Merkuriev , paved the way to use the DFE not only in bound-state calculations but in investigations of scattering processes as well. For atomic systems, however, some peculiarities of the inter-atomic interactions resulted in a limited use of DFE in both bound and scattering calculations.

One of the peculiarities is that the interaction among atoms often contains very strong repulsion at short distances which is difficult to handle numerically. In addition, it generates strong short range correlations implying that in atomic systems one should take into account many partial waves in order to achieve convergence. These problems can hardly be solved by using the computational power of modern computers and requires instead the development of appropriate numerical techniques.

One such technique is the so-called tensor-trick proposed by the Groningen group Groning ; Groning1 . Application of the method in Nuclear and Coulomb systems showed that high accuracy calculations can be made using low computer power. Another method is that of Cartesian coordinates applied for the first time in three body Faddeev calculations in Carbonell . The latter method is suitable in describing the long-range behavior of weakly three-body bound states correctly. This feature is crucial in calculating excited states close to the two-body threshold. Yet another proposed method of solving the three-body problem is that of Ref. Kvits  in which no explicit partial wave decomposition is required and only the total angular momentum of the system is used. Within this total angular momentum representation method the Faddeev operator has a very simple form that allows the construction of effective numerical schemes.

In the present work we have combined the aforementioned methods and, moreover, we propose a new approach to overcome the difficulties arising from the short range repulsion. We applied the overall procedure in bound state calculations of the He3 system using various potentials and compare our results with other results previously obtained by various groups.

In Sec. II we describe the equations to be solved and formulate the appropriate boundary conditions. In Sec. III the numerical method is presented. Our results are given in Sec. IV while our conclusion are summarized in Sec. V. Some technical details for the field practitioners are shifted in the appendix.

## II Formalism

When considering a three-body system it is convenient to denote the two-body subsystems by $`\alpha `$, $`\beta `$ and $`\gamma `$ and introduce Jacobi coordinates which in the configuration space are defined by

$$\begin{array}{ccc}& & \\ \hfill 𝐱_\alpha & =& \left(\frac{2m_\beta m_\gamma }{m_\beta +m_\gamma }\right)^{1/2}(𝐫_\beta 𝐫_\gamma )\hfill \\ \hfill 𝐲_\alpha & =& \left(\frac{2m_\alpha (m_\beta +m_\gamma )}{m_\alpha +m_\beta +m_\gamma }\right)\left(𝐫_\alpha \frac{m_\beta 𝐫_\beta +m_\gamma 𝐫_\gamma }{m_\beta +m_\gamma }\right)\hfill \end{array}$$
(1)
The coordinates corresponding to other pairs can be obtained by cyclic permutations of the subscripts $`\alpha `$, $`\beta `$, and $`\gamma `$, their relation being

$$\left(\begin{array}{c}𝐱_\beta \\ 𝐲_\beta \end{array}\right)=\left(\begin{array}{cc}S_{\beta \alpha }^{11}& S_{\beta \alpha }^{12}\\ S_{\beta \alpha }^{21}& S_{\beta \alpha }^{22}\end{array}\right)\left(\begin{array}{c}𝐱_\alpha \\ 𝐲_\alpha \end{array}\right)$$
(2)
where $`S_{\beta \alpha }^{ij}`$ are coefficients which depend on the masses of the particles MerkFadd . For identical particles these coefficients are

$$S_{\beta \alpha }^{11}=S_{\beta \alpha }^{22}=\frac{1}{2},S_{\beta \alpha }^{12}=S_{\beta \alpha }^{21}=\frac{\sqrt{3}}{2}$$

We assume that the Hamiltonian of the system involves only two-body interactions

$$H=H_0+\underset{i}{}V_i(𝐱_i),$$
(3)
where $`H_0=\mathrm{\Delta }_x\mathrm{\Delta }_y`$ is the Hamiltonian of three free particles and $`V_i(𝐱_i)`$ are the two-body potentials.

The wave function $`\mathrm{\Psi }`$ of the system can be expressed in terms of the three Faddeev components $`\mathrm{\Phi }_i`$

$$\mathrm{\Psi }(𝐱,𝐲)=\underset{i}{}\mathrm{\Phi }_i(𝐱_i,𝐲_i)$$
(4)
satisfying the Faddeev equations

$$(\mathrm{\Delta }_x\mathrm{\Delta }_y+V_i(𝐱_i)E)\mathrm{\Phi }_i(𝐱_i,𝐲_i)=V_i(𝐱_i)\underset{ji}{}\mathrm{\Phi }_j(𝐱_j,𝐲_j)$$
(5)
where $`E`$ is the energy of the system. In what follows we shall restrict ourselves, without loss of generality, to three identical bosons which is the case under consideration. In such a case the Faddeev components have the same functional form in their own coordinates and thus the system (5) is reduced to one equation only

$$(H_0+VE)\mathrm{\Phi }=V(C^++C^{})\mathrm{\Phi }.$$
(6)
In this equation $`H_0=\mathrm{\Delta }_x\mathrm{\Delta }_y`$, $`C^+`$ and $`C^{}`$ are the cyclic and anticyclic permutation operators respectively, $`\mathrm{\Phi }`$ is one of the Faddeev components $`\mathrm{\Phi }_i`$ written in its own coordinates.

The potential energy of the system is invariant with respect to rotations. This makes it possible to separate out the degrees of freedom corresponding to rotations of the system. This can be achieved by expanding the Faddeev component $`\mathrm{\Phi }`$ in terms of eigenfunctions of the total angular momentum, i.e. Wigner functions $`D_{mn}^L(g)`$ Kvits ,

$$\mathrm{\Phi }(𝐱,𝐲)=\underset{L,m,n}{}\frac{\varphi ^{Lmn}(x,y,z)}{xy}D_{mn}^L(g).$$
(7)
Here $`gSO(3)`$ are the coordinates describing collective angular motion of the system and $`\varphi ^{Lmn}(x,y,z)`$ are the projections of the Faddeev component in subspaces with fixed angular momentum. The components of the projections depend on the intrinsic coordinates

$$x=|𝐱|,y=|𝐲|,z=\frac{(𝐱,𝐲)}{xy},x,y[0,\mathrm{}),z(1,1)$$
(8)
which describe the intrinsic state of the cluster. Since we can fix the total angular momentum $`L`$ of the system, the corresponding projection of the free Hamiltonian can be written as

$$H_0^L=D^L(g^1)xyH_0\frac{1}{xy}D^L(g)$$
(9)
where $`D^L(g)`$ stands for a matrix constructed from the Wigner functions Kvits . In the case of zero total angular momentum the explicit expression for $`H_0^0`$ reads

$$H_0^0=\frac{^2}{x^2}\frac{^2}{y^2}(\frac{1}{x^2}+\frac{1}{y^2})\frac{}{z}(1z^2)\frac{}{z}.$$
(10)
The corresponding projection of Eq. (6) takes the form

$$(H_0^0+V(x)E)\varphi ^0(x,y,z)=V(x)P\varphi ^0(x,y,z),$$
(11)
where

$$P\varphi ^0(x,y,z)xy(\frac{\varphi ^0(x^+,y^+,z^+)}{x^+y^+}+\frac{\varphi ^0(x^{},y^{},z^{}))}{x^{}y^{}})$$

and $`x^\pm (x,y,z)`$, $`y^\pm (x,y,z)`$, and $`z^\pm (x,y,z)`$ are

$$\begin{array}{ccc}\hfill x^\pm (x,y,z)& =& (\frac{1}{4}x^2+\frac{3}{4}y^2\frac{\sqrt{3}}{2}xyz)^{1/2},\hfill \\ \hfill y^\pm (x,y,z)& =& (\frac{3}{4}x^2+\frac{1}{4}y^2\pm \frac{\sqrt{3}}{2}xyz)^{1/2},\hfill \\ \hfill z^\pm (x,y,z)& =& \frac{\pm {\displaystyle \frac{\sqrt{3}}{4}}x^2{\displaystyle \frac{\sqrt{3}}{4}}y^2{\displaystyle \frac{1}{2}}xyz}{x^\pm (x,y,z)y^\pm (x,y,z)}.\hfill \end{array}$$
(12)
Assuming that in each two-body subsystem only one bound state exists, we can write the asymptotic boundary conditions for the Faddeev component $`\varphi ^0`$ as follows

$$\varphi ^0(x,y,z)\phi _2(x)\mathrm{e}^{k_yy}+A(x/y,z)\frac{\mathrm{e}^{k_3(x^2+y^2)^{1/2}}}{(x^2+y^2)^{1/4}},$$
(13)
where $`\phi _2(x)`$ denotes the wave function of the two-body bound state in the two-body subsystem, $`k_y=\sqrt{E_2E_3}`$, $`k_3=\sqrt{E_3}`$, $`E_2`$ is the two–body bound state energy, and $`E_3`$ the energy of the three-body system. The first term corresponds to virtual decay into a particle and a two-body bound system, usually denoted as 2+1, while the second term corresponds to a virtual decay with an amplitude $`A(x/y,z)`$ into three single particles denoted as 1+1+1. The term corresponding to the latter configuration decreases much faster than in the 2+1.

In the present work the term corresponding to 1+1+1 virtual decay is neglected and thus the asymptotic boundary conditions for the Faddeev component at sufficiently large distances $`R_x`$ and $`R_y`$ read

$$\frac{}{x}\mathrm{ln}\varphi ^0(x,y,z)|_{x=R_x}=k_xi\sqrt{E_2},\frac{}{y}\mathrm{ln}\varphi ^0(x,y,z)|_{y=R_y}=k_y.$$
(14)
## III Numerical Procedure

The first numerical calculation with the DFE were made by Laverne and Gignoux in the early seventies GignLav . Since then, many numerical methods were proposed. In the present work we shall employ: i) Quintique splines KviHu  together with an orthogonal collocation deBoorSwartz  procedure which allows us to construct a linear system of equations corresponding to Eq. (11). In this way the bound state problem can be transformed to a generalized eigenvalue problem. ii) The tensor-trick Groning  to solve the eigenvalue problem using the restarting Arnoldi algorithm Saad .

Any regular solution of Eq. (11) fulfilling the boundary conditions (14) can be approximated by an expansion in terms of the basic functions (see Appendix A)

$$\varphi ^0(x,y,z)=\underset{i}{}f_iB_i(x,y,z)$$
(15)
where $`i`$ stands for a multi–index $`\{i_x,i_y,i_z\}`$ and $`B_i(x,y,z)=B_{i_x}(x)B_{i_y}(y)B_{i_z}(z)`$. Following the procedure described in the Appendix A we obtain an equation for the coefficients $`f_i`$:

$$(\widehat{H}_FE\widehat{I})𝐟=0,$$
(16)
where $`\widehat{I}`$ is the unit matrix and $`\widehat{H}_F`$ is the discrete analog of the Faddeev operator (see Appendix A)

$$\widehat{H}_F=\widehat{H}_0^0+\widehat{V}(\widehat{I}+\widehat{P}).$$
(17)
When investigating nuclear systems with short-range potentials, such as the Malfliet-Tjon V (MT-V) potential mt , the spectrum of the system (16) can be calculated by direct application of the restarting Arnoldi or biorthogonal Lanczos algorithm. However, an additional regularization is usually needed to accelerate the convergence. For this, one may split the operator $`\widehat{H}_F`$ in Eq. (16) into two parts, $`\widehat{H}_F=\widehat{H}_1+\widehat{H}_2`$, where $`\widehat{H}_1`$, is such that the operator $`(\widehat{H}_1E\widehat{I})^1`$ can be explicitly constructed while $`\widehat{H}_2`$ can be considered as a remainder. The technique of explicit inversion of an operator $`\widehat{H}_1E\widehat{I}`$ is known as tensor-trick Groning  and is briefly described in Appendix B for Cartesian coordinates.

Following the standard procedure we write for the eigenvalue problem

$$(\widehat{H}_1E\widehat{I})^1\widehat{H}_2𝐟=\lambda (E)𝐟.$$
(18)
where for physical solutions we have $`\lambda (E)=1`$. Eq. (18) can be solved using the Arnoldi algorithm Saad .

Although the above scheme can be applied effectively for a wide range of potentials, it is not satisfactory for interatomic potentials with a strong repulsive core. A typical example is the helium dimer potential which has a core with enormously strong repulsion and an extremely weak attractive well just enough to hold a two-body bound state. In what follows we give a brief analysis of the difficulties that arise and propose a recipe to overcome these difficulties.

Since to reproduce the correct long-range behavior of the Faddeev components Cartesian coordinates are used, the natural choice of the operator $`\widehat{H}_1`$ in (18) is

$$\widehat{H}_1=\widehat{H}_0^0+\widehat{V}.$$

In this case the problem of calculating bound states is reduced to calculation of the largest eigenvalues of the operator

$$\widehat{L}(z)(\widehat{H}_0+\widehat{V}z)^1\widehat{V}(\widehat{C}^++\widehat{C}^{})$$
(19)
To demonstrate how the presence of a short-range repulsion manifests itself in the spectral properties of the operators $`\widehat{L}(z)`$, one must compare the spectra of potentials with and without a strong repulsive core. Two such potentials are the interatomic He–He TTYPT potential TTY  and the nucleon-nucleon (MT-V) potential mt  the spectra of which are shown in Fig. 1. Two main features of the spectrum of the TTYPT potential should be pointed out, namely, that the spectrum contains numerous large negative eigenvalues and that for all values of the energy $`z`$, $`\widehat{L}(z)`$ has a number of eigenvalues close to one. According to the estimations of the Arnoldi algorithm convergence rate (see Appendix C), these are the features that make the convergence unacceptably slow.

In order to improve convergence we first look for the origin of the negative part of the spectrum. For this we consider the values $`z_i`$ such as $`\mathrm{\Lambda }_i(z_i)1`$ of the operator $`\widehat{L}(z)`$. Evidently, they are simultaneously eigenvalues of the operator

$$\widehat{H}_{\mathrm{negative}}\widehat{H}_0+\widehat{V}\widehat{V}\widehat{P}.$$

Keeping in mind that the operator $`\widehat{V}`$ stands for a potential with extremely strong short-range repulsion, we conclude that the term $`\widehat{V}\widehat{P}`$ produces extremely strong short-range attraction. Therefore, the operator $`\widehat{H}_{\mathrm{negative}}`$ contains a large number of eigenvalues in the range $`[\mathrm{},E_2]`$ and their presence suppresses the convergence of the Arnoldi algorithm for values close to one.

To eliminate the negative eigenvalues from the spectrum, we introduce, instead of $`\widehat{L}(z)`$, the operator

$$\widehat{M}(z)(\widehat{H}_0+\widehat{V}+\widehat{V}_mz)^1\left(\widehat{V}_m\widehat{V}\widehat{P}\right),$$
(20)
where $`\widehat{V}_m`$ corresponds to some strong short-range repulsive potential. It will be referred to as modifying potential. It can be shown that the eigenvalues $`\mu _i(z_i)=1`$ of the operator $`\widehat{M}(z)`$ correspond to the eigenvalues $`z_i`$ of the original Faddeev operator and this does not depend upon the modifying potential. In contrast, the eigenvalues $`\mu _i(z)=1`$ correspond to the eigenvalues $`z_i`$ of the equation

$$\left[\widehat{H}_0+2\widehat{V}_m+\widehat{V}(1\widehat{P})z_i\right]𝐟=0.$$
(21)
Obviously, these eigenvalues depend on the choice of the modifying potential $`\widehat{V}_m`$ a proper choice of which will eliminate all eigenvalues of the equation (21) in the range $`z_i(\mathrm{},E_2)`$. Thus, according to Arnoldi algorithm convergence rate estimations, this feature of $`\widehat{M}(z)`$ should dramatically improve the convergence to the maximal eigenvalue close to one.

The convergence rate can be further improved if instead of the operator $`\widehat{M}`$ a properly selected function of $`\widehat{M}`$ is used in the calculations. One such function can be an even power of $`\widehat{M}(z)`$ or the normalized Chebyshev polynomial of the second kind $`U_k(\widehat{M}(z))/U_k(1)`$. It allows an upward shifting of the lower bound of the spectrum and the increase of the separation between the physically interesting eigenvalues close to one and the rest part of the spectrum. According to the convergence rate estimations these features allow the reduction of the computation time and memory requirements considerably.

## IV Results

We applied the aforementioned numerical procedure to study the spectrum of the He3 trimer system. For this purpose, we employed the most recent He–He potentials, namely, the SAPT1, SAPT2, and HFD-B3-FC11 potentials Az97 , and compare the results with those previously obtained Roudnev  using the LM2M2 Az91  and TTYPT TTY  interactions, as well as with those obtained via other methods.

In Table 1 we present the results for the trimer ground and excited states. We also present, in the same table, the results obtained for the dimer binding energy $`E_2`$ calculated on the same grid used in three-body calculations as well as the exact results for the dimer $`E_2^{\mathrm{ex}}`$. The difference between these values can be regarded as the lower bound for the error of our approach.

In Tables 2 and 3 we demonstrate the convergence, for the LM2M2 potential, of the calculated energies with respect to the number of grid points used. It is seen that convergence is achieved with a relatively low number of the $`z`$ points ($`12`$ points) while for the other coordinates the number is much larger ($`100`$ points).

Our results for the LM2M2 potential are given in Table 4 together with those obtained by Motovilov et. al. Moto1 ; Moto2  via Faddeev-type equations (boundary condition model (BCM)) and by Nielsen et al. Fedorov  via the hyperspherical adiabatic approach. It is seen that an overall good agreement is achieved with these approaches. However, for the ground state a slightly better agreement is achieved with the results of Ref. Moto2  while for the excited state is with those of Ref. Fedorov . Comparing our results with the results of Moto2  we should emphasize, that the angular basis used by Motovilov et. al. Moto2  coincide with our spline basis for the simplest grid that can be used in the present method. Performing calculations with this simplified grid, we recovered all the digits of the ground state energy reported in Moto2 . This confirms the accuracy of their result and the suitability of the present method for bound state calculations. The better agreement for the excited state with the one obtained in Ref. Fedorov  can be attributed to the 2+1 contribution to the Faddeev component for the excited states. Taking into account this term in the BCM is very difficult, whereas the hyperspherical adiabatic approach of Nielsen et. al. Fedorov  is more suitable for this purpose. However, the ground state energy of Fedorov  is about 1% less than our result. Consideration of the geometric properties of helium trimer can clarify the possible nature of the latter difference.

The characteristic size of the bound states of the trimer can be estimated either by calculating $`r`$ or $`r^2^{1/2}`$. The results obtained for these radii are presented in Tables 5 and 6. It is seen that they are approximately 10 times less than the radii of the dimer molecule. However, the size of the excited state has the same order of magnitude as that of the dimer.

The wave function of the system can be easily obtained from the Faddeev component,

$$\psi (x,y,z)=\varphi (x,y,z)+xy\left(\frac{\varphi (x^+,y^+,z^+)}{x^+y^+}+\frac{\varphi (x^{},y^{},z^{})}{x^{}y^{}}\right),$$

where $`x^\pm `$, $`y^\pm `$ and $`z^\pm `$ are defined by (12). The most intuitive way to visualize the results of the calculations is to draw the one-particle density function defined as

$$\rho (𝐫_1)=d𝐫_2d𝐫_3|\mathrm{\Psi }(𝐫_1,𝐫_2,𝐫_3)|^2,$$

where

$$\mathrm{\Psi }(𝐫_1,𝐫_2,𝐫_3)=\frac{\psi (x(𝐫_1,𝐫_2,𝐫_3),y(𝐫_1,𝐫_2,𝐫_3),z(𝐫_1,𝐫_2,𝐫_3))}{2\pi x(𝐫_1,𝐫_2,𝐫_3)y(𝐫_1,𝐫_2,𝐫_3)}.$$

The functions $`x(𝐫_1,𝐫_2,𝐫_3)`$, $`y(𝐫_1,𝐫_2,𝐫_3)`$, and $`z(𝐫_1,𝐫_2,𝐫_3)`$ are obtained from the Jacobi coordinates (1) and (8) and the wave function $`\psi (x,y,z)`$ is normalized to unity. Due to the symmetry properties of the system, the one-particle density function depends only on the coordinate $`r=|𝐫_1|`$. Taking into account the relation $`𝐲_1=\sqrt{3}𝐫_1`$ we get

$$\rho (r)=\frac{\sqrt{3}}{4\pi ^2r^2}dxdz|\psi (x,r\sqrt{3},z)|^2.$$

Omitting the integration over $`z`$, we obtain the conditional density function $`\rho (r,z)`$ describing a spatial distribution for particle one when the other two particles are located along a fixed axis. It is useful to plot this function in coordinates $`(r_l,r_a)`$ such that $`r_l=rz`$ is a projection of the position of particle 1 onto the axis connecting the other particles and

$$r_a=\frac{z}{|z|}r(1z^2)^{\frac{1}{2}}$$

is a projection to the orthogonal axis. Three-dimensional plots of the function $`\rho ((r_l^2+r_a^2)^{1/2},\mathrm{cos}\mathrm{arctan}r_l/r_a)`$ corresponding to the ground and excited states of the trimer calculated with the LM2M2 potential are presented in Figs. 2 and 3 respectively. The conditional density function of the ground state decreases in all directions in a similar way. The density function of the excited state has two distinguishable maxima and exhibits the linear structure of the cluster. This structure has a simple physical explanation: The most probable positions of a particle in the excited state lie around the other two particles and when the latter particles are well separated the third one forms a dimer-like bound state with each of them. This interpretation agrees with the clusterisation coefficients presented in the Table 7. These coefficients are calculated as a norm of the function $`f_c`$ defined by

$$f_c(y)=dxdz\varphi (x,y,z)\phi _2(x),$$
(22)
where $`\phi _2(x)`$ is the dimer wave function. The values of $`|f_c(y)|^2`$, shown in the Table 7, demonstrate the dominating role of a two-body contribution to the trimer excited state whereas in the ground state this contribution is rather small. We could suppose that this dominating contribution of the cluster wave in the excited state ensured the fast convergence of the hyperspherical adiabatic expansion Fedorov  to the correct value, but in order to get the same order of accuracy for the ground state, possibly more basic functions should be taken into account.

The advantage of using Faddeev equations over the Schrödinger one in bound-state calculations, is deduced from the results shown in Tables 8 and 9. In the latter table we present the contribution of different angular states to the Faddeev component and to the wave function calculated as

$`C_n`$ $`=`$ $`|f_n(x,y)|^2`$
$`f_n(x,y)`$ $`=`$ $`{\displaystyle _1^1}dzF(x,y,z)P_n(z),`$ (23)
where $`P_n(z)`$, $`n=0,2,4,\mathrm{}`$, are the Legendre polynomials and $`F(x,y,z)`$ is the Faddeev component or the wave function. The angular coefficients for the Faddeev component decrease much faster than the wave function coefficients. This could explain the difference between the estimations of the trimer mean square radius reported in Roudnev  and the one presented here. In the previous paper the same angular basis for the Faddeev component and the wave function has been used. However, it turned out that the contribution of the higher partial waves to the wave function is not negligible. This leads, for instance, to the change of the mean square radius of the excited state for the LM2M2 potential from 60.85 Å to 59.3 Å. This observation agrees with the results of Nielsen et al. Fedorov , who reported the value 60.86 Å using essentially the same angular basis both for the wave function and for the Faddeev component. The Table 8 also demonstrates that more angular functions should be taken into account in the ground state calculations.

## V Conclusions

We present a method which can be used to perform accurate three-body bound states calculations. It is specially suited for systems where the interparticle forces have a strong repulsive core. The method is based on Faddeev equations in the total angular momentum representation and without any further partial wave decomposition. The equations are expressed in terms of Cartesian coordinates which are best suited to describe the long-range behavior of weakly three-body bound states. Combined with the tensor-trick preconditioning for large linear systems and the Arnoldi’s algorithm for eigenanalysis it provided us accurate results for the the He3 trimer bound and excited states.

Results obtained with the most recent realistic intermolecular forces (SAPT1, SAPT2, HFD-B3-FC11), as well as with earlier forces (LM2M2, TTYPT) indicate that only two bound states exist. The properties of these states are very different with the ground state being strongly bound while the binding energy for the excited state is comparable to that of the dimer. The latter implies that a dimmer cluster within the molecule is well formed.

The sizes of these two states also differs much. The characteristic size of the ground state either estimated by $`r`$ or $`r^2^{1/2}`$ is $`10`$ times less than the size of dimer molecule, but the size of the excited state has the same order of magnitude with that of the dimer. This estimation make it necessary to check for the absence of trimers in the experimental media during the measurement of dimer properties and vice versa.

## Appendix A Reduction to a Linear System

Suppose $`\tau _t=\{t_0,t_1,\mathrm{},t_{N_t}\}`$ is a partition in the range $`[t_0,t_{N_t}]`$ in coordinate $`t`$ and let $`S_{5,2}(\tau _t)`$ be the space of quintique Hermite splines associated with this partition i.e. piecewise fifth order polynomial functions having two continuous derivatives in this range. The partition $`\tau _t`$ will be refereed to as the base grid. The set of basic functions in the space $`S_{5,2}(\tau _t)`$ can be defined by the following conditions

$`S_i^1(t_i)=1,_tS_i^1(t_i)=0,_t^2S_i^1(t_i)=0`$
$`S_i^2(t_i)=0,_tS_i^2(t_i)=1,_t^2S_i^2(t_i)=0`$ (24)
$`S_i^3(t_i)=0,_tS_i^3(t_i)=0,_t^2S_i^3(t_i)=1`$
$`t[t_{\mathrm{max}(0,i1),},t_{\mathrm{min}(i+1,N_t)}]S_i^j(t)=0`$.

For simplicity we omit the subscript $`\tau _{N_t}`$ corresponding to any particular choice of a grid in each coordinate. Let the basic functions be arranged as follows:

$$\widehat{B}_i(t)=S_k^j(t),k=0\mathrm{}N_t,j=1\mathrm{}3,i=3k+j1$$
(25)
The solution of the Eq. (11) must satisfy Eq. (14) and vanish on the planes $`x=0`$ and $`y=0`$. To achieve this we modify the set of basic functions (24) by constructing linear combinations $`B_i(t)`$ which for $`t=x,y`$ read

$`B_i(t)\widehat{B}_i(t),i=1\mathrm{}3N_t1`$
(26)
$`B_{3N_t}(t)\widehat{B}_{3N_t}(t)k_t\widehat{B}_{3N_t+1}(t)+k_t^2\widehat{B}_{3N_t+2}(t)`$
and for the coordinate $`z`$ they coincide with the initial basic functions (25) $`B_i(z)=\widehat{B}_i(z)`$. We look for an approximate solution of Eq. (11) in the form of expansion in terms of the basic functions (26)

$$\varphi ^0(x,y,z)=\underset{i}{}f_iB_i(x,y,z)$$
(27)
where $`i`$ stands for the multi–index $`\{i_x,i_y,i_z\}`$ and $`B_i(x,y,z)=B_{i_x}(x)B_{i_y}(y)B_{i_z}(z)`$. Since the basic functions (26) fulfill the boundary conditions (14) these conditions are also satisfied by the approximate component $`\varphi ^0`$.

To obtain a system of equations for the coefficients $`f_i`$ the expansion (27) should be substituted into Eq. (11) with some subsequent projection procedure. A reasonable choice of such a procedure is the method of orthogonal collocations which do not require any integration. The highest order of approximation is achieved by a careful choice of collocation points for a given basic grid $`\tau `$ deBoorSwartz . For the grids $`\tau _{N_x}`$, $`\tau _{N_y}`$ and $`\tau _{N_z}`$ one should construct the grids of collocation points $`\tau _{N_x^c}=\{x_1^c,x_2^c,\mathrm{},x_{N_x^c}^c\}`$, $`\tau _{N_y^c}=\{y_1^c,y_2^c,\mathrm{},y_{N_x^c}^c\}`$, and $`\tau _{N_z^c}=\{z_1^c,z_2^c,\mathrm{},z_{N_x^c}^c\}`$ containing $`N_x^c=3N_x`$, $`N_y^c=3N_y^c`$, and $`N_z^c=3(N_z+1)`$ points respectively. Given these collocation grids the matrix elements of the operators involved in Eq. (11) can be easily constructed. For example, the matrix elements of the identity operator are $`[I]_{ij}=B_i(x_{j_x}^c,y_{j_y}^c,z_{j_z}^c)`$ with $`i=\{i_x,i_y,i_z\}`$ and $`j=\{j_x,j_y,j_z\}`$.

Substituting the expansion (27) into Eq. (11) followed by the orthogonal collocation procedure one gets a system of linear algebraic equations

$$(\widehat{H}_0^0+\widehat{V}(\widehat{I}+\widehat{P})E\widehat{I})𝐟=0.$$
(28)
where $`\widehat{H}_0^0`$, $`\widehat{V}`$, $`\widehat{I}`$, and $`\widehat{P}`$ are the matrices corresponding to the operators of Eq. (11) while $`𝐟`$ is the coefficient vector of the expansion (27).

The dimension of this system is rather large and equals $`N=N_x^cN_y^cN_z^c`$. Therefore, we use this system of equations to find only its lowest generalized eigenvalues and eigenvectors.

## Appendix B The Tensor-trick for Cartesian Coordinates

The inversion procedure of the operator $`(\widehat{H}_1E\widehat{I})`$ is crucial for the performance of the algorithm and therefore a brief description of a tensor-trick algorithm for the DFE in Cartesian coordinates will be given. Furthermore, as the matrices involved are large, certain comments will be passed concerning the efficiency in computation.

Consider the Hamiltonian $`\widehat{H}_1=\widehat{H}_0^0+\widehat{V}\widehat{I}`$. The structure of this operator is

$`\widehat{I}`$ $`=`$ $`\widehat{I}_x\widehat{I}_y\widehat{I}_z`$
$`\widehat{H}_1`$ $`=`$ $`\widehat{D}_x\widehat{I}_y\widehat{I}_z+\widehat{I}_x\widehat{D}_y\widehat{I}_z`$
$`+`$ $`\widehat{I}_x\widehat{I}_y\widehat{1}_z(\widehat{L}_x\widehat{1}_y+\widehat{1}_x\widehat{L}_y)\widehat{D}_z+\widehat{V}\widehat{I}_x\widehat{I}_y\widehat{I}_z.`$
The $`\widehat{I}_t`$ are the matrices corresponding to the identity operators acting in coordinate $`t`$, $`\widehat{1}_t`$ is the unit matrix, and $`\widehat{D}_t`$ contains the elements of the corresponding differential operators. The tensor structure of the operator $`\widehat{H}`$ allows a diagonalisation of it without solving the spectral problem in the $`N`$–dimensional space. To perform this diagonalisation one introduces the following matrices $`\widehat{B}_z`$, $`\widehat{B}_x^i`$, $`\widehat{B}_y^i`$ constructed from the generalized eigenvectors of the operators acting in different coordinates,

$`\widehat{B}_z`$ $`:`$ $`\widehat{D}_z\widehat{B}_z=\widehat{L}_z\widehat{I}_z\widehat{B}_z`$
$`\widehat{B}_x^i`$ $`:`$ $`(\widehat{D}_x+\widehat{V}\widehat{I}_x+l_i\widehat{L}_x\widehat{I}_x)\widehat{B}_x^i=\widehat{X}^i\widehat{I}_x\widehat{B}_x^i`$
$`\widehat{B}_y^i`$ $`:`$ $`(\widehat{D}_y+l_i\widehat{L}_y\widehat{I}_y)\widehat{B}_y^i=\widehat{Y}^i\widehat{I}_y\widehat{B}_y^i`$
where $`\widehat{L}_z=\mathrm{diag}\{l_1,l_2,\mathrm{},l_{Nz}\}`$, $`\widehat{X}^i=\mathrm{diag}\{E_{x1}^i,E_{x2}^i,\mathrm{},E_{xN_x}^i\}`$, and $`\widehat{Y}^i=\mathrm{diag}\{E_{y1}^i,E_{y2}^i,\mathrm{},E_{yN_y}^i\}`$, and $`E_{tN_t}^i`$ are the eigenvalues of the corresponding operators. One also defines the matrices $`\widehat{B}_z^{}`$, $`\widehat{B}_{x^i}^{}`$, $`\widehat{B}_{i_y}^{}`$ constructed of the eigenvectors of the adjoined equations. Expanding the expression $`\widehat{B}^{}(\widehat{H}_1E\widehat{I})\widehat{B}`$, with $`\widehat{B}`$ and $`\widehat{B}^{}`$ being given by

$`\widehat{B}`$ $``$ $`{\displaystyle \underset{i=1}{\overset{N_z}{}}}(\widehat{B}_x^i\widehat{B}_y^i)\widehat{[B}_z]_i`$
$`\widehat{B}^{}`$ $``$ $`{\displaystyle \underset{i=1}{\overset{N_z}{}}}(\widehat{B^i}_x^{}\widehat{B^i}_y^{})\widehat{[B}_z^{}]_i,`$
and taking into account the relation $`\widehat{B}^{}\widehat{B}=\widehat{1}_x\widehat{1}_y\widehat{1}_z`$, one can prove that the matrices $`\widehat{B}`$ and $`\widehat{B}^{}`$ are the desired matrices diagonalising the operator $`(\widehat{H}_1E\widehat{I})`$. Once this operator is diagonalised it can then be easily inverted

$`\widehat{B}^{}(\widehat{H}_1E\widehat{I})\widehat{B}=\widehat{G},`$
$`(\widehat{H}_1E\widehat{I})^1=\widehat{B}\widehat{G}^1\widehat{B}^{}`$
$`\widehat{G}=\mathrm{diag}\{g_{111},g_{112},\mathrm{},g_{i_xi_yi_z},\mathrm{},g_{N_xN_yN_z}\}`$
$`g_{i_xi_yi_z}=E_{xi_x}^{i_z}+E_{xi_y}^{i_z}E`$
Thus, in principle, the task of construction of the operator $`(\widehat{H}_1E\widehat{I})^1`$ is over. However, to use computer power effectively one must care about optimal implementation of all the operations involved in the computation of the product of the operator $`\widehat{h}^1(E)(\widehat{H}_1E\widehat{I})^1`$ and a vector $`𝐮`$, the number of arithmetic operations should be minimized. To implement the operation effectively the structure of the operator $`\widehat{h}^1(E)`$ can be used

$$\widehat{h}^1(E)\widehat{B}\widehat{G}^1(E)\widehat{B}^{}.$$

The matrix $`\widehat{G}^1(E)`$ has explicit diagonal structure, and the calculation of the product $`\widehat{G}^1𝐮`$ requires only $`O(N_xN_yN_z)`$ operations. Since however, the matrices $`\widehat{B}`$ and $`\widehat{B}^{}`$ are dense, and in principle the calculation of the products $`\widehat{B}^{()}𝐮`$ requires $`O((N_xN_yN_z)^2)`$ arithmetic operations, that could essentially decrease the effectiveness of the method. However, taking into account the structure of the matrices $`\widehat{B}`$ and $`\widehat{B}^{}`$ one can improve the performance. Consider the procedure of tensor multiplication of two matrices by a vector $`𝐮`$:

$$𝐯=\widehat{B_x}\widehat{B_y}𝐮,$$

and use the identity

$$\widehat{B_x}\widehat{B_y}=\widehat{I}\widehat{B_y}\widehat{B_x}\widehat{I}.$$

Obviously, even the calculation of the product $`\widehat{B_x}\widehat{B_y}𝐮`$ requires $`N_x^2N_y^2`$ operations, the calculation of the products $`\widehat{I}\widehat{B_y}𝐮`$ and $`\widehat{B_x}\widehat{I}𝐮`$ require correspondingly only $`N_x^2N_y`$ and $`N_xN_y^2`$ operations. Therefore a multiplication of a tensor product by a vector would cost $`N_xN_y(N_x+N_y)`$ operations. Similar relations can be written for a direct sum. Thus the correct implementation of a direct sum multiplication by a vector can reduce the computation costs $`N_xN_yN_z/(N_x+N_y+N_z)`$ times.

The computation costs can be reduced even further. One of the possibilities is to abandon full diagonalisation of the operator $`\widehat{h}(z)`$ and to modify the computation scheme. For this, let us construct the matrices diagonalising only the terms corresponding to the coordinates $`y`$ and $`z`$

$`\stackrel{~}{B}{\displaystyle \underset{i=1}{\overset{N_z}{}}}(\widehat{I}\widehat{B^i}_y)\widehat{[B_z]}_i,`$
$`\stackrel{~}{B}^{}{\displaystyle \underset{i=1}{\overset{N_z}{}}}(\widehat{I}\widehat{B^i}_y^{})\widehat{[B_z^{}]}_i.`$
In this case the matrix $`\stackrel{~}{G}(z)\stackrel{~}{B}^{}\widehat{h}(z)\stackrel{~}{B}`$ is band because of the fact that the operator is local and the basic functions have well localized support (26). The multiplication $`\widehat{L}(z)𝐮`$ can be performed in two stages:

* Solve a system of equations with a band matrix $`\stackrel{~}{G}(E)`$

$$\stackrel{~}{G}(z)\stackrel{~}{𝐮}=\stackrel{~}{B}^{}(\widehat{V}\widehat{P}+\widehat{V}_c\widehat{S})𝐮;$$
* calculate $`𝐯\widehat{L}(z)𝐮`$ as

$$𝐯=\stackrel{~}{B}\stackrel{~}{𝐮}.$$

Taking into account the tensor structure of the matrices $`\stackrel{~}{B}^{}`$ and $`\stackrel{~}{B}`$ this computation scheme requires only $`O(N_xN_yN_z(N_y+N_z))`$ operations.

## Appendix C  Convergence rates for the Arnoldi algorithm

The key point of this algorithm is the repeated action of the operator under consideration on a sequence of vectors generating a Krylov subspace where the desired eigenvector lies. To use this algorithm to solve the eigenvalue problem (18) one must be able to calculate the action of the operators $`(\widehat{H}_1E\widehat{I})^1`$ and $`\widehat{H}_2`$ on an arbitrary vector $`𝐮`$. Thus the problems of inversion of the operator $`(\widehat{H}_1E\widehat{I})`$ and storage of the matrix $`\widehat{H}_2`$ should be solved first. In Ref. Groning  hyperspherical coordinates were used which lead to an effective storage scheme for the matrices $`\widehat{H}_2=\widehat{V}(\widehat{I}+\widehat{P})`$. However, one can not include pair interaction into the invertible part $`(\widehat{H}_1E\widehat{I})`$ using hyperspherical coordinates. On the other hand, if Cartesian coordinates are used then one has difficulties with the storage of the matrix $`\widehat{P}`$ since these coordinates do not have the same symmetry properties as the operator $`P`$ and this results to an irregular structure of $`\widehat{P}`$. In our case the situation is saved by the total angular momentum representation. In contrast to bispherical expansion used in Groning , the operator $`P`$ is local and does not contain any terms requiring integration. As a result the matrix $`\widehat{P}`$ is sparse and can be either effectively stored using any standard storage scheme for sparse matrices sparce  or not stored at all. In the last case the matrix elements are calculated “on fly” when required.

The estimations of the convergence rate of the Arnoldi algorithm Saad  are essential in understanding the problems that arise when the tensor-trick algorithm is applied. Suppose $`\widehat{L}`$ is a linear operator in a finite-dimensional space with

$$\widehat{L}𝐮_i=\lambda _i𝐮_i.$$

Let $`𝒦_m(𝐯)`$ be a Krylov subspace of dimension $`m`$ constructed from a starting vector $`𝐯`$ of the operator $`\widehat{L}`$ and $`𝒫_m`$ be an orthogonal projector to $`𝒦_m`$ space. The convergence rate of the algorithm can be expressed in terms of the norm

$$\epsilon _i^{(m)}(1𝒫_m)𝐮_i$$

of the eigenvector $`u_i`$ projection onto the orthogonal completion to $`𝒦_m`$. Suppose $`\alpha _k`$ are the expansion coefficients of the starting vector $`𝐯`$ in terms of the eigenvectors $`𝐮_i`$,

$$𝐯=\underset{i=1}{\overset{n}{}}\alpha _i𝐮_i.$$

Then the following relation is fulfilled

$$\epsilon _i^{(m)}\xi _iϵ_i^{(m)},$$

where $`\xi _i=1+_{k=1}^n|\alpha _k|/|\alpha _i|`$ and $`ϵ_i^{(m)}`$ depends on the spectral properties of the operator $`\widehat{L}`$. For the eigenvalues lying in the external part of the spectrum some practically important estimations are known. Suppose that all the eigenvalues of $`\widehat{L}`$ but $`\lambda _1`$ lye inside of an ellipse with the center $`c`$, largest semi-axis $`a`$, and focuses $`c\pm e`$. It can then be shown that $`ϵ_1^{(m)}`$ satisfy the following estimation Saad

$$|ϵ_1^{(m)}|\frac{C_{m1}({\displaystyle \frac{a}{e}})}{C_{m1}({\displaystyle \frac{\lambda _1c}{e}})},$$

where $`C_m(t)`$ are Chebyshev polynomials of the first kind.

Consider now the convergence rate of Arnoldi algorithm for the operator $`\widehat{L}(z)`$ (19) for the value of z corresponding to the ground state of a three-body system. In this case the maximal eigenvalue $`\lambda _1=1`$. Having defined the separation distance between the largest eigenvalue and the rest part of the spectrum $`s|1\lambda _2|`$, the distance between the maximal eigenvalues in the rest part $`d|\lambda _2\lambda _n|`$ and neglecting the imaginary part of the spectrum, we obtain the estimation

$$|ϵ_1^{(m)}||C_{m1}(1+\frac{2s}{d})|^1.$$

Therefore the parameter managing the convergence rate is the ratio $`s/d`$ of the separation to the size of the rest part of the spectrum. The convergence rate could be improved by enlarging the separation and by enhancing the remainder part of the spectrum localization.


format: text

text:
# Partition function zeros of the one-dimensional Potts model: The recursive method

## 1 Introduction

It is well established that the thermodynamic properties of a physical system can be derived from a knowledge of the partition function. Since the discovery of statistical mechanics it has been a central theme to understand the mechanism how the analytic partition function for a finite-size system acquires a singularity in the thermodynamic limit if the system undergoes a phase transition. The answer to this quest was given in 1952 by Lee and Yang in their famous papers . They considered the partition function of the Ising model as a polynomial in activity ($`\mathrm{exp}(2H/kT)`$, where H is a magnetic field) and studied the distribution of zeros of the partition function in the complex activity plane. It was shown that phase transitions occur in the systems where continuous distribution of zeros of the partition function cuts the real axis in the thermodynamic limit. Also proved, was the circle theorem, which states that zeros of the partition function of the ferromagnetic Ising model lie on the unit circle in the complex activity plane (Yang-Lee zeros). Later, Fisher  initiated a study of zeros of the partition function in the complex temperature plane (Fisher zeros). Fisher showed that complex temperature zeros of the partition function of the Ising model in zero magnetic field on a square lattice lie on two circles $`|v\pm 1|=\sqrt{2}`$, where $`v=\mathrm{tanh}(J/2kT)`$. Since that time there has been a considerable body of work studying the partition function zeros of the Ising and Potts  models on various regular lattices, most notably in recent years . Zeros of the partition function were also studied for spin models defined on hierarchical  and recursive  lattices, random graphs  and aperiodic systems , spin glasses , percolation and Self-Organized Criticality models .

Nowadays, the investigation of zeros of the partition function becomes a powerful tool for studying phase transitions and critical phenomena. Recently, much attention has been attracted to the study of zeros of the partition function in one-dimensional systems in connection with helix-coil transitions in biological macromolecules . Hence, it is very important to study the general properties of zeros of the partition functions for different one-dimensional systems.

In 1994, Glumac and Uzelac  using the transfer matrix method studied the Yang-Lee zeros of the one-dimensional ferromagnetic Potts model for non-integer values of $`Q0`$. They showed that for $`0<Q<1`$ the Yang-Lee zeros are located on a real interval and for low temperatures these are located partially on the real axis and in complex conjugate pairs on the activity plane. Later on, Monroe  numerically studied the Yang-Lee zeros of the Potts model for some particular values of $`Q`$. Then, Kim and Creswick  showed that for $`Q>1`$ the Yang-Lee zeros lie on a circle with radius $`R`$, where $`R<1`$ for $`1<Q<2`$, $`R>1`$ for $`Q>2`$ and $`R=2`$ for $`Q=2`$. Only recently the full picture of Yang-Lee zeros of the ferromagnetic Potts model for $`0<Q<1`$ has been found .

In this paper, the dynamical systems theory is used to study the Yang-Lee, Fisher and Potts zeros111Zeros of the partition function considered as a function of complex $`Q`$. of the partition function for the one-dimensional $`Q`$-state Potts model. In Section 2, a recurrence relation for the partition function is derived. It is shown that zeros of the partition function may be associated with neutral fixed points of the recurrence relation. A general equation for zeros of the partition function is derived. Formulae for the free energy and the density of zeros of the partition function are found. In Sections 3 and 4 the method developed in Section 2 used to study the Yang-Lee and Potts zeros of ferromagnetic and antiferromagnetic Potts models. A classification of Yang-Lee and Potts zeros is given. In Section 5, the Fisher zeros in a nonzero magnetic field are investigated. It is shown that the Fisher zeros in a nonzero magnetic field located on several lines in the complex temperature plane and the density of Fisher zeros is singular at the edge singularity points with the same critical exponent as that of both the Yang-Lee and Potts zeros, $`\sigma =\frac{1}{2}`$.

## 2 Zeros of the partition function of the Potts model

The Hamiltonian of the one-dimensional $`Q`$-state Potts model in a magnetic field is defined as follows

$$=\text{}\underset{<ij>}{}\delta (\sigma _i,\sigma _j)\text{}\underset{i}{}\delta (\sigma _i,0).$$
(1)
where $`\delta `$ is the Kroneker delta function, $`\sigma _i`$ denotes the Potts variable at site $`i`$ and takes the values $`0,1,2,\mathrm{},Q1`$. The first sum in the r.h.s. of (1) goes over all edges and the second one over all sites on the lattice. For $`\text{}>0`$ the model is ferromagnetic and for $`\text{}<0`$ is antiferromagnetic. Note that due to the symmetry, the Hamiltonian (1) is the same if the external field  is applied to any spin state, namely, if $`\delta (\sigma _i,0)`$ in (1) is replaced by $`\delta (\sigma _i,\alpha )`$ for any $`\alpha =1,\mathrm{\hspace{0.17em}2},\mathrm{},Q1`$. For $`Q=2`$ the Potts model corresponds to the Ising model and in order to keep the analogy with the Ising model we designate  as a magnetic field. We may assume the cyclic boundary condition $`\sigma _n=\sigma _n`$ and that the number of sites is $`2n+1`$ without losing a generality. Cutting the lattice at the central site $`\sigma _0`$ will separate it into two branches I and II with equal statistical weights $`g_n(\sigma _0)`$ (Figure 1).

Cutting the branch I (II) at the site $`\sigma _1`$ ($`\sigma _1`$), the recurrence relation for $`g_n(\sigma )`$ may be found as

$$g_n(\sigma _0)=\underset{\sigma _1=0}{\overset{Q1}{}}\mathrm{exp}\left[J\delta (\sigma _0,\sigma _1)+h\delta (\sigma _1,0)\right]g_{n1}(\sigma _1).$$
(2)
where $`J=\text{}/kT`$ and $`h=\text{}/kT`$. The partition function may be written in the form

$$𝒵=\underset{\sigma }{}\mathrm{}^{/kT}=\underset{\sigma _0=0}{\overset{Q1}{}}\mathrm{exp}\left[h\delta (\sigma _0,0)\right]g_n^2(\sigma _0).$$
(3)
Introducing the notation

$$x_n=\frac{g_n(\sigma 0)}{g_n(\sigma =0)},$$
(4)
and using (2), the recurrence relation for $`x_n`$ may be found as

$$x_n=f(x_{n1}),f(x)=\frac{\mu +(z+Q2)x}{z\mu +(Q1)x},$$
(5)
where $`\mu =\mathrm{}^h`$, and $`z=\mathrm{}^J`$. $`x_n`$ does not have a direct physical meaning, but the thermodynamic functions, such as the magnetization, the specific heat and etc., may be expressed in terms of $`x_n`$. For example, the magnetization for our model in the thermodynamic limit has the form

$$m=(N𝒵)^1\underset{\{\sigma \}}{}\delta (\sigma ,0)\mathrm{}^{/kT}=\frac{\mu }{\mu +(Q1)x},$$
(6)
where $`x`$ is an attracting fixed point of the mapping (5). It will be shown below how the thermodynamic properties of the model may be defined from the dynamics of the recurrence relation (5).

The mapping (5) is a Möbius transformation, i.e. a rational map of the form

$$R(x)=\frac{ax+b}{cx+d},adbc0,$$

where

$$R(\mathrm{})=a/c,R(d/c)=\mathrm{},$$

if $`c0`$, while $`R(\mathrm{})=\mathrm{}`$ when $`c=0`$. The dynamics of such maps is well studied . The mapping (5) has only two fixed points which are solutions to the equation $`f(x)=x`$. According to the theory of complex dynamical systems, fixed points are classified as follows: A fixed point $`x^{}`$ is attracting if $`|\lambda |<1`$, repelling if $`|\lambda |>1`$, and neutral if $`|\lambda |=1`$, where $`\lambda =\frac{}{x}f(x)|_{x=x^{}}f^{}(x^{})`$ is called the eigenvalue of $`x^{}`$. It is easy to show that either both fixed points of the mapping (5) are neutral, or that one of them is attracting while the other is repelling. The correspondence between thermodynamical properties of the model and the dynamics of the mapping (5) is the following: If for a given temperature and magnetic field ($`z`$ and $`\mu `$) the mapping (5) has an attracting fixed point, then the system is in a stable state and its thermodynamical functions are defined by this fixed point (see for example (6)). The other fixed point is repelling and does not correspond to any phase. On the other hand, if the iterations of (5) do not converge to a fixed point, i.e. the mapping (5) has neutral fixed points only, the system undergoes a phase transition. Therefore, the existence of neutral fixed points of (5) corresponds to a phase transition in the model. It was mentioned in the Introduction that zeros of the partition function correspond to phase transitions in the model. Thus, zeros of the partition function are associated with neutral fixed points of the corresponding mapping. These may be found from the conditions of existence of neutral fixed points of the mapping (5). These conditions are the following

$$\{\begin{array}{c}f(x)=x\hfill \\ f^{}(x)=e^{i\varphi },\varphi [0,2\pi ],\hfill \end{array}$$
(7)
Excluding $`x`$ from the system (7) after some algebra the equation of phase transitions may be found

$$z^2\mu ^22\left[(z1)(z+Q1)\mathrm{cos}\varphi +1Q\right]\mu +\left(z+Q2\right)^2=0,$$
(8)
where $`\varphi [0,2\pi ]`$. Solutions to this equation for different values of $`\varphi `$ correspond to zeros of the partition function. Hence, the free energy of the model may be written in the form

$$F_0^{2\pi }\mathrm{ln}(z^2\mu ^2+2(Q1)\mu +\left(z+Q2\right)^22\mu (z1)(z+Q1)\mathrm{cos}\varphi )\varphi .$$
(9)
Formula (9) may also be derived from the transfer matrix method  using two non-degenerate eigenvalues of the transfer matrix

$$\lambda _{1,2}^{}=\frac{1}{2}\left[z\mu +z+Q2\pm \sqrt{(z\mu zQ+2)^2+4(Q1)\mu }\right],$$
(10)
and the following mathematical identity for any pair of scalars $`C`$, $`D`$

$$C^N+D^N=\underset{n}{}\left[C+\mathrm{exp}(2\pi \mathrm{}n/N)D\right],$$

where the product is from $`n=1,\mathrm{\hspace{0.17em}2},\mathrm{},N`$ if $`N`$ is odd; and from $`n=1/2,\mathrm{\hspace{0.17em}3}/2,\mathrm{},N1/2`$ for N even. Moreover, one can easily show that the eigenvalues of fixed points of (5) are related to the eigenvalues of transfer matrix (10) as follows

$$\lambda _{1,2}=\frac{4\mu (z1)(z+Q1)}{(\lambda _{1,2}^{})^2}.$$

Hence, the condition of a phase transition $`|\lambda _1|=|\lambda _2|`$ (both $`\lambda _1`$ and $`\lambda _2`$ are neutral fixed points) corresponds exactly to the $`|\lambda _1^{}|=|\lambda _2^{}|`$ condition for the eigenvalues of transfer matrix (10). It means that the phase transition point based on dynamical systems approach coincides exactly with the phase transition point based on free energy considerations. This correspondence seems to be general for spin models defined on recursive lattices .

The density of zeros of the partition function may be found from the equation (8) in a standard way and has the form

$$g(\xi )=\frac{B_\xi AA_\xi B}{2\pi B[(AB)(A+B)]^{\frac{1}{2}}},$$
(11)
where $`_\xi =\frac{}{\xi }`$,

$$A=z^2\mu ^2+2(Q1)\mu +\left(z+Q2\right)^2,B=2\mu (z1)(z+Q1)\mathrm{cos}\varphi ,$$

and $`\xi =\mu ,z`$ or $`Q`$ depending on whether the Yang-Lee, Fisher or Potts zeros are considered. It is interesting to note that in our case the numerator of (11) always contains the multiplier $`(A+B)`$. Hence, the density $`g(\xi )`$ is singular only when $`AB=0`$222Solutions to the equation B=0 should be neglected since these do not correspond to zeros of the partition function.. From (11) follows that $`g(\xi )`$ has a singular behavior $`g(\xi )|\xi \xi ^{}|^\sigma `$, where $`\xi ^{}`$ is a solution to the equation $`AB=0`$. We will see that $`\sigma =\frac{1}{2}`$ for all types of zeros of the partition function. Also, note that the equation $`AB=0`$ corresponds to the phase transitions equation (8) for $`\varphi =0`$. It defines the edge singularity points . In the subsequent Sections we will apply equations (8) and (11) to the study of the Yang-Lee, Potts and Fisher zeros.

## 3 The Yang-Lee zeros

According to the results of the previous section the Yang-Lee zeros of the $`Q`$-state Potts model may be found by solving the equation (8) with respect to $`\mu `$. Equation (8) is a quadratic equation of $`\mu `$ with real coefficients. Note that solutions to the equation (8) lie either on the real axis or in complex conjugate pairs on a circle with radius $`R=|z+Q2|/z`$ and have the form

$$\mu _{1,2}=E\left[2\mathrm{cos}^2\frac{\varphi }{2}F\pm 2\sqrt{\mathrm{cos}^2\frac{\varphi }{2}\left(\mathrm{cos}^2\frac{\varphi }{2}F\right)}\right],$$
(12)
where

$$E=\frac{(z1)(z+Q1)}{z^2}\text{and}F=\frac{z(z+Q2)}{(z1)(z+Q1)}.$$

A detailed study of (8) with respect to $`\mu `$ has already been performed . Here, we will present only the main results for the Yang-Lee zeros of both ferromagnetic and antiferromagnetic Potts models.

For the ferromagnetic Potts model ($`z>1`$) one can find that for $`Q>1`$ all solutions (12) are complex conjugate and lie on an arc of circle with radius $`R=(z+Q2)/z`$. Writing $`\mu `$ in the exponential form $`\mu =Re^{i\theta }`$ the angular distribution of the Yang-Lee zeros may be found in the form

$$\mathrm{cos}\frac{\theta }{2}=\sqrt{\frac{(z1)(z+Q1)}{z(z+Q2)}}\mathrm{cos}\frac{\varphi }{2}.$$
(13)
From (13) one can see a gap in the distribution of Yang-Lee zeros, i.e. there are no solutions to the equation (8) in the interval $`\theta _0<\theta <\theta _0`$, where $`\theta _0=2\mathrm{arccos}\sqrt{F^1}`$. This is the well known gap in the distribution of Yang-Lee zeros of ferromagnetic models above the critical temperature (formally $`T_c=0`$ for the one-dimensional case) first studied by Kortman, Griffiths and Fisher . The end points of the gap are called Yang-Lee edge singularity points. From (8) and (13) follows that the Yang-Lee edge singularity points correspond to $`\varphi =0`$. Hence, according to (11) the density of Yang-Lee zeros is singular in the Yang-Lee edge singularity points. Substituting $`\varphi =0`$ into (12) the formula for Yang-Lee edge singularity points we have

$$\mu _\pm =\frac{1}{z^2}\left\{\sqrt{(z1)(z+Q1)}\pm \sqrt{1Q}\right\}^2.$$
(14)
Substituting $`\xi =\mu R\mathrm{}^\mathrm{}\theta `$ in (11) after some algebra, the density of Yang-Lee zeros for $`Q>1`$ may be found in the form

$$g(\theta )=\frac{1}{2\pi }\frac{|\mathrm{sin}\frac{\theta }{2}|}{\sqrt{\mathrm{sin}^2\frac{\theta }{2}\mathrm{sin}^2\frac{\theta _0}{2}}},$$
(15)
where $`\theta _0=2\mathrm{arccos}\sqrt{F^1}`$. From the equation (15) follows that the density $`g(\theta )`$ diverges in the Yang-Lee edge singularity points $`\mu _\pm `$ with the critical exponent $`\sigma =\frac{1}{2}`$, i.e. $`g(\theta )|\theta \theta _0|^{\frac{1}{2}}`$ when $`\varphi 0`$ or $`\theta \theta _0`$.

For $`Q<1`$ the Yang-Lee edge singularity points $`\mu _\pm `$ are real and the density of Yang-Lee zeros has the form

$$g(\mu )=\frac{1}{2\pi \mu }\frac{|\mu \sqrt{\mu _+\mu _{}}|}{\sqrt{(\mu _+\mu )(\mu \mu _{})}}.$$
(16)
$`g(\mu )`$ diverges in the points $`\mu _\pm `$, i.e. $`g(\mu )|\mu \mu _\pm |^\sigma `$, with the critical exponent $`\sigma =\frac{1}{2}`$. The summary of results for the Yang-Lee zeros of the ferromagnetic Potts model is given in Figure 2.

The Yang-Lee zeros of the antiferromagnetic Potts model ($`z<1`$) may be studied in the same manner and the results are given in Figure 3. For complex $`\mu _\pm `$ the angular distribution of Yang-Lee zeros has the form

$$\mathrm{sin}\frac{\theta }{2}=\sqrt{\frac{(z1)(z+Q1)}{z(z+Q2)}}\mathrm{cos}\frac{\varphi }{2}.$$
(17)
In contrast to the ferromagnetic case now the Yang-Lee zeros lie in the interval $`\theta _0<\theta <\theta _0`$, where $`\theta _0=2\mathrm{arcsin}\sqrt{F^1}`$. The density of Yang-Lee zeros in this case has the form

$$g(\theta )=\frac{1}{2\pi }\frac{\mathrm{cos}\frac{\theta }{2}}{\sqrt{\mathrm{sin}^2\frac{\theta _0}{2}\mathrm{sin}^2\frac{\theta }{2}}}.$$
(18)
For real values of $`\mu _\pm `$ the density of Potts zeros has the same form as for the ferromagnetic case (16). Actually, the formula (16) is the most general from which formulae (15) and (18) may be derived (See also Section 4). Note, that here also the density of Yang-Lee zeros diverges at the Yang-Lee edge singularity points with the same index $`\sigma =\frac{1}{2}`$.

## 4 The Potts zeros

Recently, much attention has been attracted to the study of zeros of the partition function in the complex $`Q`$ plane for the Potts model on regular lattices and finite graphs . The partition function of the Potts model ($`H=0`$) on a finite graph $`G`$, $`𝒵_G(Q,v)`$, may be considered as a polynomial in the number $`Q`$ of Potts states and the temperature-like variable $`v=\mathrm{exp}(\text{}/kT)1`$. At zero-temperature the partition function of the antiferromagnetic Potts model ($`v=1`$) correspond to the chromatic polynomial $`P_G(Q)`$, which is closely related to the $`Q`$-coloring problem. By definition, for graph $`G`$ and positive $`Q`$, $`P_G(Q)`$ is the number of ways in which the vertices of $`G`$ can be assigned “colors” from the set $`1,2,\mathrm{},Q`$ in such a way that adjacent vertices always receive different colors. The original hope was that study of the real or complex zeros of $`P_G(Q)`$ might lead to an analytic proof of the Four-Color Conjecture, which states that $`P_G(4)>0`$ for all loop-less planar graphs. To date this hope has not been realized, although combinatoric proofs of the Four-Color Theorem have been found . Even so, the Potts zeros and zeros of $`P_G(Q)`$ are interesting in their own right and have been extensively studied in recent years .

In this section the Potts zeros of ferromagnetic and antiferromagnetic one-dimensional Potts models are studied. It was shown in Section 2 that zeros of the partition function correspond to solutions of the equation (8), which is a polynomial in $`Q`$. Hence, the Potts zeros may be found as solutions to (8) with respect to the parameter $`Q`$. It is convenient to rewrite (8) using the variable $`P=Q+z1`$

$$P^22P\left[1\mu +\mu (z1)\mathrm{cos}\varphi \right]+(1z\mu )^2=0.$$
(19)
Equation (19) is a quadratic equation in $`P`$ with real coefficients. Solutions of (19) lie either on the real axis or on the circle with radius $`R=|1z\mu |`$ in the complex plane $`P`$, and have the form

$$P_{1,2}=1\mu +\mu (z1)\mathrm{cos}\varphi \pm 2\mathrm{cos}\frac{\varphi }{2}\sqrt{\mu (z1)(1\mu \mu (z1)\mathrm{sin}^2\frac{\varphi }{2})}.$$
(20)
Since the analysis of (20) is standard, we will skip the details and give the results in Figures 4 and 5 for ferromagnetic and antiferromagnetic Potts models respectively. In analogy with the Yang-Lee zeros, the Potts edge singularity points are defined as solutions to the equation (19) for $`\varphi =0`$, and have the form

$$Q_\pm =P_\pm z+1=1z+\left(\sqrt{1\mu }\pm \sqrt{\mu (z1)}\right)^2$$
(21)
The general formula of the density of Potts zeros has the from

$$g(P)=\frac{1}{2\pi P}\frac{|P\sqrt{P_+P_{}}|}{\sqrt{(P_+P)(PP_{})}}.$$
(22)
From (22) it follows that the density of Potts zeros is singular at the Potts edge singularity points. For complex values of $`Q_\pm `$, (22) may be written in the form

$$g(\theta )=\frac{1}{2\pi }\frac{|\mathrm{sin}\frac{\theta }{2}|}{\sqrt{\mathrm{sin}^2\frac{\theta }{2}\mathrm{sin}^2\frac{\theta _0}{2}}},$$
(23)
for the ferromagnetic model, where $`|\theta |>\theta _0`$ and $`\theta _0=2\mathrm{arccos}(\mu (z1)/(\mu z1))^{\frac{1}{2}}`$. For the antiferromagentic model, $`g(\theta )`$ has the form

$$g(\theta )=\frac{1}{2\pi }\frac{\mathrm{cos}\frac{\theta }{2}}{\sqrt{\mathrm{sin}^2\frac{\theta _0}{2}\mathrm{sin}^2\frac{\theta }{2}}},$$
(24)
where $`|\theta |<\theta _0`$ and $`\theta _0=2\mathrm{arcsin}(\mu (1z)/(1\mu z))^{\frac{1}{2}}`$. Different formulae for the density of Potts zeros occur because the square root function in (22) is not a unique function in the complex plane $`P`$, i.e. one branch of it corresponds to the ferromagnetic model and the other to that of the antiferromagnetic.

It is noteworthy that the density of Potts zeros diverges at the Potts edge singularity points with the same critical exponent $`\sigma =\frac{1}{2}`$ as for the Yang-Lee zeros.

## 5 The Fisher zeros

In this section the Fisher zeros in a nonzero magnetic field are studied. Usually, the Fisher zeros are considered as zeros of the partition function with respect to a temperature-dependent parameter. In our case $`z`$ is such a parameter. Formally, the magnetic field in the Hamiltonian (1) may be presented in the form $`\text{}=H\text{}`$, where $`H`$ is the renormalized magnetic field. Later we will refer to $`H`$ as a magnetic field. Then, $`\mu `$ may be written as $`z^H`$ and the equation (8) will has the form

$$P_\varphi (z,H,Q)=0,$$
(25)
where

$$P_\varphi (z,H,Q)=z^{2H+2}2\mathrm{cos}\varphi z^{H+2}2(Q2)\mathrm{cos}\varphi z^{H+1}+2(Q1)z^H+\left(z+Q2\right)^2.$$

$`P_\varphi (z,H,Q)`$ is obviously a polynomial for integer values of $`H`$ and (25) may be solved numerically. For non-integer values of $`H`$ the equation (25) becomes a transcendental equation and the Fisher zeros may be found by numerically checking the condition of existence of neutral fixed points of the mapping (5).

The density of Fisher zeros may be found from (11) by substituting $`\xi =z`$ and has the form

$$g(z)=\frac{1}{2\pi }\frac{HG_1(z,H,Q)zG_2(z,H,Q)}{(z1)z(z+Q1)\sqrt{P_0(z,H,Q)}},$$
(26)
where

$`G_1(z,H,Q)=(z1)(z+Q1)(z^{H+1}zQ+2),`$
$`G_2(z,H,Q)=(z^H1)(\mathrm{\hspace{0.17em}2}(z+Q1)Qz)+Q^2.`$
From (26) it follows that the density of Fisher zeros is singular in the Fisher edge singularity points, which are defined as solutions of the equation (25) for $`\varphi =0`$.

Let us first study the Fisher zeros for the Ising model ($`Q=2`$). For $`Q=2`$ the equation (25) has the form

$$z^{2H+2}2\mathrm{cos}\varphi z^{H+2}+4\mathrm{cos}^2\frac{\varphi }{2}z^H+z^2=0.$$
(27)
This equation is symmetric under the $`HH`$ transformation because of the $`Z(2)`$ symmetry of the Ising model. Since all coefficients of (27) are real, its solutions will be either real or complex conjugate. For integer values of $`H`$ the Fisher zeros have the form shown in Figure 6.

From Figure 6 and (27) one can see that the Fisher zeros located on lines ending at the Fisher edge singularity points and the number of lines equals to the value of the magnetic field $`|H|+1`$. Fisher edge singularity points are divided on pairs and every pair corresponds to a line. Note that $`z=0`$ is a twice degenerate Fisher edge singularity point with the critical exponent $`\sigma =\frac{1}{2}`$. Formally, two degenerate $`z=0`$ edge singularity points are considered here as forming a “line”. One can see that for odd values of $`|H|`$ there are negative Fisher zeros and for $`|H|=4n`$, $`nN`$, the Fisher zeros are partially located on the imaginary axis. Since $`P_\varphi (z,H,Q)`$ is an analytic function of $`H`$ the location of Fisher zeros for non-integer values of $`H`$ are continuous deformations of Fisher zeros in the field $`[|H|]`$ or $`[|H|]+1`$, where $`[|H|]`$ is the integer part of $`|H|`$, i.e. the minimal integer less than $`|H|`$. The location of Fisher zeros for non-integer $`H`$ may be found numerically from the condition of existence of neutral fixed points of the mapping (5). As an illustration the Fisher zeros for non-integer values of $`0<H<2`$ are given in Figure 7.

The Fisher zeros for $`Q2`$ Potts model may be studied in the same way as for the Ising model. It is impossible to give all possible configurations of Fisher zeros for any $`H`$ and $`Q`$. Here we give only the summary of the main properties. First of all, there is no $`HH`$ symmetry for $`Q2`$ and the Fisher zeros are different for positive and negative values of $`H`$.

For $`Q=1`$ the Fisher zeros located on a closed curve and the density of Fisher zeros is not singular on this curve.

For $`Q1`$ Fisher zeros are located on several lines and the density of Fisher zeros is singular at the end points of these lines with the edge singularity exponent $`\sigma =\frac{1}{2}`$.

Numerical experiments show that for integer values of $`H`$ there are the following properties: for $`0<Q<1`$ there is an interval of real Fisher zeros only for even values of $`H`$; for $`1<Q<2`$, the Fisher zeros intersect the negative real semi-axis only for odd values of $`H`$; for $`Q>2`$ there is an interval of real Fisher zeros only for odd values of $`H`$. The number of the lines of Fisher zeros are defined as in the Ising model. In Figure 8 we give some plots that illustrate these properties.

## 6 Conclusions

In this paper, the Yang-Lee, Fisher and Potts zeros of the one-dimensional $`Q`$-state Potts model are studied using the dynamical systems theory. A recurrence relation for the partition function is derived. It is shown that zeros of the partition function may be associated with neutral fixed points of the corresponding recurrence relation. A general equation for zeros of the partition function is derived. It is shown that the density of zeros of the partition function is singular in the edge singularity points with the critical exponent $`\sigma =\frac{1}{2}`$.

The recursive method used in this paper is equivalent to the transfer matrix method, but in some sense it is more flexible. It gives a numerical algorithm to study the zeros of the partition function. The algorithm is based on testing the condition of existence of neutral fixed points of the corresponding recurrence relation: If neutral fixed points exist, it means that the partition function is zero for given values of $`z`$, $`\mu `$ and $`Q`$. In Section 5 we saw that this algorithm is very useful for studying the Fisher zeros in an arbitrary magnetic field. Moreover, it is the only method for studying the Fisher zeros for non-integer values of $`H`$, where $`H`$ is the renormalized magnetic field  and $`H=\text{}/\text{}`$.

In conclusion we would like to note that the results given in this paper show that the thermodynamic properties of the one-dimensional Potts model are completely defined by the recurrence relation (5) or equivalently, by two non-degenerate eigenvalues of the corresponding transfer matrix (10) (see also ). Moreover, it is proved that the phase transition point based on dynamical systems approach coincides exactly with the phase transition point based on free energy considerations. The dynamical systems approach used here gives a numerical method which may be used for studying other one-dimensional systems for which a one-dimensional recurrence relation may be derived.

One of the authors (R G) would like to acknowledge the Abdus Salam Centre for Theoretical Physics for hospitality extended during his visit where this work was partially performed.

## References


text:
# A quantum phase gate implementation for trapped ions in thermal motion

## Abstract

We propose a novel scheme to implement a quantum controlled phase gate for trapped ions in thermal motion with one standing wave laser pulse. Instead of applying the rotating wave approximation this scheme makes use of the counter-rotating terms of operators. We also demonstrate that the same scheme can be used to generate maximally entangled states of $`N`$ trapped ions by a single laser pulse. 03.67.Lx, 03.65.Ud ,42.50.Dv.

The existence of quantum algorithms for specific problems shows that a quantum computer can in principle provide a tremendous speed up compared to classical computers . This discovery motivated an intensive research into this mathematical concept which is based on quantum logic operations on multi-qubit systems . In order to implement this concept into a real physical system a quantum system is needed, which makes the storage and the read out of quantum information and the implementation of the required set of quantum gates possible. This system should be scalable and the isolation of the system from the environment should be very well in order to suppress decoherence processes. Several physical systems were suggested to implement the concept of quantum computing: cavity QED systems , trapped ion systems , nuclear magnetic resonance systems , solid state systems with nuclear spins , quantum dot systems  and Josephson-junction device systems . The ion trap quantum information processor, first proposed by Cirac and Zoller , was demonstrated to be a promising candidate for the realization of a small scale quantum computations . The original idea of Cirac and Zoller uses one mode of the collective ion vibration for a bus concept. The implementation of the Cirac/Zoller-gate requires to cool the mode of the collective vibration to the ground state. Any heating process on this mode, which is caused by the environment, will diminish the accuracy of the Cirac/Zoller-gate. Thus, an unreliable performance of the quantum computer as a whole will result, which makes the implementation of even simple algorithms impossible. The elimination of heating processes is a very demanding goal but not possible. Several gate proposal do not pose the requirement of ground state cooling . In this paper we propose a novel scheme to implement a quantum phase gate for trapped ions in thermal motion. In contrast to other schemes  only one standing wave laser pulse is needed to implement a quantum controlled phase gate. Our scheme does not require the rotating wave approximation. We also show that the same scheme can be used to generate maximally entangled states of $`N`$ trapped ions by using a single standing wave laser pulse.
In the case of one trapped ion interacting with a standing wave laser field, scheme have been proposed to generate and measure vibrational state of a trapped ion. Here we intend to consider $`N`$ trapped ions with ground state $`|g_i>`$ and excited state$`|e_i>`$ in a linear trap, which can interact with a monochromatic standing wave of frequency $`\omega _L`$. We assume that the ion-laser interaction strength characterized by Rabi frequency $`\mathrm{\Omega }`$ is the same for all ions and we also assume that the position of the center of the all ions relative to a node of the laser field is same, which is characterized by $`\theta `$. The total system is composed of the center of mass vibrational freedom of the ion string and the internal electronic freedom of the ions. The vibrational freedom is represented by a harmonic oscillator with ladder operators $`a^{}`$ and $`a`$ and frequency $`\omega _a`$ and the internal freedom is represented by Pauli matrices $`\sigma _\alpha ^{(i)}`$, $`\alpha =x,y,z`$. In a frame, which rotates with the standing wave frequency $`\omega _L`$, such a system can be described by the Hamiltonian

$`H=\omega _aa^{}a+{\displaystyle \underset{i=1}{\overset{N}{}}}\left[{\displaystyle \frac{\mathrm{\Delta }}{2}}\sigma _z^{(i)}+\mathrm{\Omega }\mathrm{sin}(kx+\theta )\sigma _x^{(i)}\right].`$ (1)
Here $`\sigma _z^{(i)}=|e_i><e_i||g_i><g_i|`$ and $`\sigma _x^{(i)}=|e_i><g_i|+|g_i><e_i|`$. The parameter $`\mathrm{\Delta }=\omega _0\omega _L`$ denotes the detuning between the standing wave laser field and the atomic frequency $`\omega _0`$. $`kx=\eta (a+a^{})`$ and $`\eta `$ is Lamb-Dicke parameter. In the Lamb-Dicke limit, the Hamiltonian(1) can be written as

$`H=\omega _aa^{}a+{\displaystyle \frac{\mathrm{\Delta }}{2}}J_z+\mathrm{\Omega }[\eta (a+a^{})\mathrm{cos}\theta +\mathrm{sin}\theta ]J_x.`$ (2)
Here we introduce collective spin operator $`J_\alpha =_{i=1}^N\sigma _\alpha ^{(i)}`$, $`\alpha =x,y,z`$. In the following, we consider a standing wave pulse, which is in resonance with the carrier transition: $`\mathrm{\Delta }=0`$. In this case, the Hamiltonian (2) simplifies to the form

$`H=\omega _aa^{}a+\mathrm{\Omega }[\eta (a+a^{})\mathrm{cos}\theta +\mathrm{sin}\theta ]J_x.`$ (3)
The exact time evolution operator of this Hamiltonian (3) is

$`U(t)`$ $`=`$ $`\mathrm{exp}\left[iJ_x^2\left({\displaystyle \frac{\eta ^2\mathrm{\Omega }^2\mathrm{cos}^2\theta t}{\omega _a}}{\displaystyle \frac{\eta ^2\mathrm{\Omega }^2\mathrm{cos}^2\theta \mathrm{sin}\omega _at}{\omega _a^2}}\right)\right]\mathrm{exp}(i\omega _aa^{}at)`$ (4)
$`\times \mathrm{exp}\left[{\displaystyle \frac{\eta \mathrm{\Omega }\mathrm{cos}\theta }{\omega _a}}J_x\left(a^{}(e^{i\omega _at}1)a(e^{i\omega _at}1)\right)\right]\mathrm{exp}(iJ_x\mathrm{\Omega }\mathrm{sin}\theta t).`$
Note, that the rotating wave approximation is not used. But our scheme uses a standing wave field of a particular time length $`\tau =2\pi /\omega _a`$. The time evolution operator reduces considerably

$`U(\tau )=\mathrm{exp}\left[{\displaystyle \frac{2i\pi \eta ^2\mathrm{\Omega }^2\mathrm{cos}^2\theta }{\omega _a^2}}J_x^2\right]\mathrm{exp}\left[{\displaystyle \frac{2i\pi \mathrm{\Omega }\mathrm{sin}\theta }{\omega _a}}J_x\right].`$ (5)
It is easy to see, at the time $`\tau `$, two subsystem are disentangled. The vibrational motion is returned to its original state, be it the ground state or any vibrational excited state, and we are left with an internal state evolution, which is independent of the external vibrational state. In contrast to other proposal, where the transformation (4) was realized by a sequence of laser pulses  our proposal needs only a single standing wave pulse.
One possible application of the scheme is implementation of a quantum phase gate of two ions by using a single pulse laser. For two trapped ions system, we choose the parameters $`\eta `$, $`\mathrm{\Omega }`$, $`\theta `$ and $`\omega _a`$ to satisfy

$`\mathrm{sin}\theta =\sqrt{{\displaystyle \frac{\eta ^2\pi }{\eta ^2\pi +4}}}\text{and}{\displaystyle \frac{\omega _a}{\mathrm{\Omega }}}=\sqrt{{\displaystyle \frac{64\eta ^2}{\eta ^2\pi ^2+4\pi }}}.`$ (6)
In this case, the time evolution operator becomes

$`U(\tau )=\mathrm{exp}\left[{\displaystyle \frac{i\pi }{8}}(J_x^2+2J_x)\right].`$ (7)
It is easy to check that this time evolution operator represents a quantum phase gate

$`|>_1|\pm >_2`$ $``$ $`|>_1|\pm >_2`$
$`|+>_1|>_2`$ $``$ $`|+>_1|>_2`$
$`|+>_1|+>_2`$ $``$ $`|+>_1|+>_2`$ (8)
in the basis states $`|\pm >_i=\frac{1}{\sqrt{2}}(|e>_i\pm |g>_i)`$.
In the following we will briefly demonstrate that it is possible to prepare a multiparticle entanglement of $`N`$ trapped ions

$`|\mathrm{\Psi }>={\displaystyle \frac{1}{\sqrt{2}}}[e^{i\phi _g}|gg\mathrm{}g>+e^{i\phi _e}|ee\mathrm{}e>].`$ (9)
by using a single pulse laser, irrespective of N even or odd. Quantum states of this kind were used to improve the frequency standard. Several schemes  were proposed to generate this kind of quantum states by a sequence of laser pulses. Our proposal shows an improvement over these schemes, since it requires only a single standing wave laser pulse. The quantum state (9) can be generated by this pulse, if the system was initially prepared in the ground state $`|ggg\mathrm{}>`$. If $`N`$ is even, the parameter $`\theta `$ has to be set to zero. If we choose this parameter to satisfy $`\frac{\mathrm{\Omega }}{\omega _a}=\frac{1}{4\eta }`$, the time evolution operator will at the time $`\tau =2\pi /\omega _a`$ transform the initial state with $`U(\tau )=\mathrm{exp}\left[\frac{i\pi }{8}J_x^2\right]`$. In Ref. such a kind of time evolution operator was used to generate quantum states of the form (9) with $`\phi _g=\frac{\pi }{4}`$ and $`\phi _e=\frac{\pi }{4}+\frac{N\pi }{2}`$. In addition, it is also pointed that for the time evolution operator (7) can be used to generate the maximally entangled state (9) of odd $`N`$ trapped ions.
In summary, we proposed a novel scheme to implement a quantum phase gate and entangle the electronic states of $`N`$ ions in a very efficient way by a single standing wave laser pulse. This is demonstrated in the case of the maximally entangled quantum state (9). The analysis of the system’s dynamics, which we present in this paper, doesn’t use the rotation wave approximation. In contrast to other schemes only the requirement of small Lamb-Dicke parameters is needed. We calculate the exact time evolution and set the time interval of the laser pulse exactly to the time of one period of the collective ion vibration. This requirement helps to fix the quantum number of the collective vibration. Since our scheme is not dependent on the value of the vibration quanta, it can be applied on ion systems in thermal motion. We pointed out how this scheme makes use of the counter-rotating operator terms. For the same purpose of implementing quantum phase gate and generating GHZ state of N ions, other schemes  need three laser pulses or more.


text:
# Unveiling the structure of the planetary nebula M 2-48:

## 1 Introduction

The morphology of M 2-48 has been previously discussed by Vázquez et al. (2000) \[hereafter Paper I\] who found the full extent of its bipolar structure and indications of the presence of high-velocity outflows from their monochromatic images. The occurrence of highly collimated, high velocity jet-like components in planetary nebulae is well established (e.g. Miranda & Solf 1992; López et al. 1995; Sahai & Trauger 1998; Miranda, Guerrero & Torrelles 1999; O’Connor et al. 2000; Miranda, Guerrero & Torrelles 2001 and Gonçalves, Corradi & Mampaso 2001, etc.). These bipolar outflows usually show indications of episodic outburst and rotation or precession of the symmetry axis (e.g. López et al. 1997). The origin of such outflows is intimately related to the formation and shaping of PNe and thus their kinematic characterisation is a key parameter for the models. Several models have been proposed for the generation of collimated outflows in PNe, such as pure hydrodynamical models (Icke et al. 1992, Mellema 1995, 1997), models involving stellar rotation and magnetized winds (Różyczka & Franco 1996; García-Segura 1997; García-Segura & López 2000), and the influence of binary system (Soker & Livio 1994; Reyes-Ruíz & López 1999, Soker 2002).

The bipolar planetary nebula M 2-48 has three main structural components: (1) an internal bipolar central region (CR); (2) an off-center semi-circular fragmented shell, containing a region with a bow-shock morphology (this has been interpreted in Paper I as the possible interaction of the bipolar, collimated outflow with the shell), and (3) a pair of symmetric knots with a seemingly bow-shock morphology, located at either side at $`2^{}`$ from the center (see Paper I).

In order to analyze the kinematical structure, excitation mechanisms and chemical abundances of M 2-48, spatially resolved high and low dispersion spectroscopy have been now obtained. We first describe the observations and the results extracted from the position-velocity (PV) diagrams in section 2. The kinematical structure and morphology are discussed in section 3, including a simple model to predict the orientation of the shocked structures with respect to the plane of the sky. In section 4, the excitation mechanism and chemical abundances at different regions of the nebula are derived. Finally, our main conclusions are summarized in section 5.

## 2 Observations and Results

Low-dispersion (LD) spectroscopy was obtained in 1999, July 18, 19 and 20, with the B&Ch spectrometer combined with the f/7.5 focus of the 2.1-m OAN telescope at the Sierra de San Pedro Mártir, B. C., Mexico. A 300 lines/mm grating was combined with a $`220\mu `$m ($`2.9`$ arcsec) slit width and a Tektronix CCD of $`1024\times 1024`$ pixel, $`24\mu `$m square pixels ($`0.84`$ arcsec) as detector. The resulting spectra cover from 3400 to 7500 Å, with a 12 Å spectral resolution. The slit was oriented in several position angles (PA) to cover different regions of the PN. Figure 1 shows the slit positions over a contour plot of the \[N II\] 6584 Å image (Paper I). The LD CCD frames were trimmed, bias-subtracted, flat-fielded and sky-subtracted with standard techniques using IRAF. Flux calibration was performed using sensitivity functions derived from the standard star HD 192281. Table 1 presents the log for each of the 1D spectra that have been extracted and analyzed and Figure 2 shows the one-dimensional spectra of CR and B1.

The high-dispersion (HD) spectroscopy was obtained in 1999, June 29 and 30, with the Manchester Echelle Spectrometer (MES; Meaburn et al. 1984) combined with the f/7.9 focus at the 2.1-m OAN telescope at the Sierra de San Pedro Mártir, B. C., Mexico. This spectrometer has no cross-dispersion. A 90 Å bandwith filter was used in order to isolate the 87th order, which contains the H$`\alpha `$ and \[N II\] $`6548,6584`$ Å emission lines.

A Tektronix CCD with $`1024\times 1024`$, 24 $`\mu m`$ ($`0.3`$ arcsec) square pixels was the detector. A 2 $`\times `$ 2 binning was employed in both spatial and spectral dimensions. The “seeing” varied between 1 and 2 arcsec during the observations. In Figure 3, the slit positions A-F are shown against a contour plot of \[N II\] $`6584`$ Å of M 2-48. The 150 $`\mu m(1.9`$ arcsec) wide slit was oriented at a position angle PA=75 for slits A – D and North-South for slits E –F. The spectral resolution is $`10`$ km s-1. Exposure times of 1800 s were carried out for each slit position. Images were reduced in the standard manner using IRAF. The spectra were wavelength calibrated with a Th-Ar lamp to an accuracy of $`\pm 1`$ km s-1.

The bidimensional arrays of \[N II\] 6584 Å position-velocity (PV) line profiles from the HD spectra (slits A-D), are shown in Figure 4. Figures 5, 6 and 7 show enlargements of the core region for slits B and C; the line profiles from slits E and F and the line profiles at the tip of north eastern lobe for slits B and C, respectively.

## 3 Morphology and Kinematics

M 2-48 has a bipolar morphology, as has been described in Paper I (see also the images in Manchado et al. 1996). From a bright bipolar core emerge relatively faint, broken filaments that delineate fairly tight collimated lobes extending over 2 arcmin on either side. Some of the filaments close to the core seem to trace an incomplete arc that gives the appearence of being the remnants of a spherical shell. The PV diagrams shown in Figures 4 – 7 reveal the complex kinematic structure of the main components of M 2-48. The labels B1, s1, s2, s4 and s5 in the PV arrays in Figure 4 match the location of the regions with the same labels in Figure 3, where we have followed the nomenclature of Paper I. The spectral element labeled as (f1) is probably related to the tip of the southwestern lobe.

### 3.1 Central region

Line profiles from the bright nebular core are shown in detail in Figure 5. The systemic velocity of the core, as defined by the central intensity maximum of the line profiles shown in Figure 5 lies around $`15`$ km s-1. The line profiles replicate a bipolar outflow emerging from a bright, compact core, extending bluewards to -50 km s-1 on the west side and redwards to $`+25`$ km s-1 on the east side. The line profile from slit C (bottom panel in Fig. 5) shows particularly extended wings, reaching remarkable values of $`\pm 85`$ km s-1 at FWZI within the central ($``$ 20 arcsec) region. The bright compact core itself is seen elongated and tilted in the velocity space, expanding towards negative velocities on the east side ($`25`$ km s-1) and positive velocities on the west side ($`+10`$ km s-1).

### 3.2 The arc-like structure

One of the curious features of M 2-48 is a group of several knots around the central bipolar region that seem to form an arc-like structure (Paper I). These knots (s1-s5 in Figure 1) were interpreted as possible segments of an expanding semi-circular shell. Some slits were obtained with the purpose of intersecting this possible expanding shell and clarify their nature. The line profiles shown in Figure 6, are found to follow the velocity pattern of the bipolar outflow emerging from the extended regions of the core, with redshifted velocities on the east side ($``$ $`+10`$ km s-1) and blueshifted velocities on the west side ($``$ $`25`$ km s-1). No indications of a velocity ellipse that could hint to the presence of an expanding shell are detected, at least with the present spectral resolution.

### 3.3 Bow-shock s1

The eastern lobe is the brighter one, with several high-speed knots with bow-shock looking morphologies being detected in the HD spectra, namely B1, s1 and s2. On the western lobe we detect and identify knots s4 and s5, and the feature f1 (see Figures 3 and 4).

A high-velocity element, labeled s1 in Figure 4, coincident with the position of a bright knot on the eastern side and near a portion of the arc-like structure, shows an apparent bow- shock morphology and kinematical structure (see Figures 3 and 8, bottom right and the image in Manchado et al. 1996). It is of interest to estimate the orientation angle of the bow-shock associated with knot s1 with respect to the plane of the sky, $`\varphi `$. To do that, we have used a “3/2-D” model to predict the position-velocity diagram of a bow-shock at different angles with respect to the plane of the sky. Raga & Noriega-Crespo (1993) found that the H$`\alpha `$ emission per unit time and area is approximately given by :

$$I_{\mathrm{H}\alpha }\left(v_{shock}\right)^{\gamma _{\mathrm{H}\alpha }},$$
(1)
where $`v_{shock}`$ is the shock velocity and $`\gamma _{\mathrm{H}\alpha }`$ is a constant dependent on the shock velocity regime.

If we assume that the jet is moving along the $`z`$-axis, the shape of the bow-shock can be parametrized by :

$$rz^p,$$
(2)
where $`p`$ determinates the aperture of the bow-shock. We have used a value of $`p=2`$.

Thus, the input parameter for the model essentially is the shock velocity. This velocity can be estimated from the PV diagrams obtained from the long-slit spectroscopy (Figure 4).

In addition, we use a value of $`\gamma _{\mathrm{H}\alpha }=3.49`$ corresponding to the 20.0–81.3 km s-1 range of velocities proposed by Raga & Noriega-Crespo (1993).

We can calculate the H$`\alpha `$ emission, and predict the position-velocity diagram at different orientations with respect to the sky plane, as shown in Figure 8. From this comparison we can establish that s1 is inclined at an angle $`10^{}`$ with respect to the plane of the sky.

This low inclination angle indicates that the radial velocity of the outflow producing the s1 bow-shock can be a factor of about 6 larger than the velocities derived from the PV diagrams.

We have used the radial velocity of the peak emission of s1, 15 km s-1, as the expansion velocity projected to the line of sight. Using the inclination estimated for s1, $`\varphi 10^{}`$, we can determine the tangential expansion velocity. Finally, taking into account the projected linear separation between s1 and the central region we can estimate a dynamical time of $`\tau _{\mathrm{s1}}3\left[d/kpc\right]\times 10^3`$ yr.

In addition, a peculiar feature labelled f2 can be distinguished in the PV diagram of slit C (Figure 4), with a constant heliocentric velocity of about $``$60 km s-1. This feature seems to be spatially coincident with the eastern arc, corresponding to the diffuse material between knots s1 and s2.

### 3.4 External bow-shock B1

In the \[N II\] 6584 Å direct images (see Paper I), it is possible to see two symmetrical structures with a bow-shock morphology. The eastern knot, B1, which is brighter than the western one, B2, is blueshifted and has a line width of $`80`$ km s-1. The PV diagrams of this bow-shock are shown in Figure 7, where we can see two different zones of B1 covered by the slits B and C. Slit B slightly covers the center of the bow-shock (larger line width), and C covers the wings (smaller line width).

It is possible to use the “3/2-D” model described in section 3.2 to predict the PV diagram of this bow-shock. From the comparison of the observations and the numerical model of the bow-shock B1 for different orientations (see Figure 9), we can estimate that B1 is inclined at an angle of $`10^{}`$ pointing inwards to the plane of the sky. The angle with respect to the plane of the sky suggests that the radial velocity of the bow-shock B1 can be a factor of $`6`$ greater than the velocities derived from the PV diagrams. In order to establish a sequence of events, using the same procedure as for s1 (see section 3.3), we estimate the dynamical age of this bow-shock to be $`\tau _{\mathrm{B1}}5.3\left[d/kpc\right]\times 10^3`$ yr. This is consistent with the idea that the bow-shocks B1 and s1 correspond to different bipolar ejections and probably with different axis orientation. This suggests a change in the orientation of the symmetry axis of the succesive outflows. In this way, a previous ejection of material would form the expanding arcs interacting with a bipolar outflow and the most external bow-shocks B1 and B2, and a later high-velocity bipolar outflow with different orientation would form the s1 bow-shock.

We have also noted that the symmetry axis of the most internal bipolar region is not coincident with either of these two high-velocity bipolar outflows.

In addition to the main features described in previous sections, we can see other interesting features in the PV diagrams obtained with the HD spectroscopy. In Figure 4 we can distinguish a very weak feature labelled f1 in slits B and C. This is a component with an heliocentric velocity of $``$ $`50`$ km s-1. This emission could correspond to the wing of the western bipolar lobe outflow. Observations at other slit positions or longer exposure times would be necessary to resolve in more detail the kinematical behavior of this feature.

## 4 Line intensity ratios and chemical abundances

As noted in section 2, five regions were selected from the slit positions observed with low resolution spectroscopy (see Figure 1), covering the four different morphological zones pointed out in Paper I. In Figure 2, the flux calibrated spectra of two different zones are shown. The top and middle panels correspond to different scales (high and low intensities, respectively) of the spectrum of the central region. The bottom panel corresponds to the spectrum of B1 (PA 47).

The extinction coefficient, C(H$`\beta `$), was derived for three regions: Central Region (CR), s1 and B1 (PA 47) from the comparison of observed and theoretical H$`\alpha `$/H$`\beta `$ ratios. We have used the calculations of Storey & Hummer (1995) assuming case B and an iteration procedure adopting finally the electron densities ($`N_e`$) derived for each particular zone and the electron temperature ($`T_e`$) obtained for CR ($`T_e`$=10800 K) for all of the regions. The reddening law of Seaton (1979) has been used. Dereddened line intensity ratios with respect to H$`\beta `$ as well as the values of C(H$`\beta `$) for these regions are presented in Table 2. The uncertainties in the line fluxes are of the order of 5% for the brighter lines and about 20-30% for the weaker ones. Emission line fluxes of blended lines such as the H$`\alpha `$+\[N II\] nebular lines and the \[S II\] doublet have been obtained from Gaussian fitting of the line profiles using the Starlink DIPSO package (Howarth & Murray 1990). The dereddened, integrated H$`\beta `$ flux of each region is also included in Table 2.

$`N_e`$ and $`T_e`$ (see Table 2) have been derived following standard methods of collisionally excited line ratios using the five-level program for the analysis of emission-line nebulae of Shaw & Dufour (1995). There seems to be a gradient in $`N_e`$, with parameter decreasing values towards the external zones. This is consistent with the expansion of stellar ejecta as the distance towards the nucleus increases. $`T_e`$ derived from \[O III\] and \[N II\] lines for CR are almost coincident, indicating that there is no strong temperature stratification in the bright bipolar nebular core of the nebula. The $`T_e`$ derived for s1 seems too high to be accounted for by pure photoionization, confirming contamination by shock emission in this zone, which morphologically and kinematically resembles a bow-shock (see section 3.3).

The spectra of the other two zones: s4 and B1 (PA 50), are so faint that H$`\beta `$ is not measured and therefore it is not possible to derive the extinction coefficient for these regions. The emission line ratios with respect H$`\beta `$ and the uncorrected H$`\alpha `$ flux for these regions are included in Table 2.

We have used several diagnostic diagrams in order to evaluate which is the dominant excitation mechanism prevalent in these regions. In Figure 10 we show two diagrams involving the locus of several line ratios: H$`\alpha `$/\[N II\], H$`\alpha `$/\[S II\], and \[S II\] 6717/6731 for H II regions, PNe and supernova remnants (adapted from Sabbadin et al. 1977). These diagrams clearly indicate that CR is mainly radiatively excited. However, the positions of s1 and B1 are not so definite, in the sense that some shock contribution is present. The fact that the positions of B1 and s1 are just below the locus of the SNRs in the upper pannel of Figure 10 could be also due to the high N enrichment of the nebula (see below).

In Figure 11 we show a diagram adapted from Phillips & Cuesta (1999) involving \[O III\]/H$`\alpha `$ and \[S II\]/H$`\alpha `$. This diagram includes the locus of the observed emission line ratios for a large number of radiatively excited PNe and predictions from plane-parallel and bow-shock models by Hartigan et al. (1987) and Shull & McKee (1979). The position of the different regions of M 2-48 in this diagram confirms that CR is radiatively excited and that shock excitation is contributing to the spectra of s1 and B1. Assuming purely shock excitation for these zones and comparing their line intensities with the models by Hartigan et al. (1987) we find the best agreement for models with shock velocities of 60 and 130 km s-1, respectively, in agreement with the values obtained in previous sections.

Ionic abundances of the CR are shown in Table 3 and have been derived using the emission line ratios given in Table 2, the code of Shaw & Dufour (1995), and assuming Te = 10 800 K and Ne = 1260 cm-3. To derive the He+/H+ and He++/H+ ratios we have used the effective recombination coefficients given by Péquignot et al. (1991) for all the lines involved. The He+ abundances obtained from each individual line have been corrected for collisional contributions following Kingdon & Ferland (1995). The final, adopted value of the He+ abundance is the average of the values obtained for the different individual lines excluding HeI 7065 Å, which suffers from the largest collisional effects.

The total O abundance has been determined from the derived abundance of O+ and O++ combined with the ionization correction factor (ICF) for the presence of O3+ proposed by Kingsburgh & Barlow (1994). The total He/H ratio has been obtained simply by adding the $`<`$He+/H$`{}_{}{}^{+}>`$ and He++/H+ ratios.

The log N+/O+=+0.30 obtained for the nebula clearly indicates that it is a nitrogen-enriched object. Taking into account the chemical classification proposed by Peimbert & Torres-Peimbert (1987), M 2-48 should be classified as a Type I PN (log N/O $``$ $``$0.30 and/or He/H $``$ 0.125), although it is not a helium-rich object. The O/H ratio is slightly lower than the average of the Galactic Type I PNe (12+log O/H = 8.65, Kingsburgh & Barlow 1994).

## 5 Conclusions

The kinematic mapping of M 2-48 has revealed a complex bipolar structure characterised by a series of shocked regions extending accross the major axis of the nebula and reaching the tips found in paper I. Indications that the bipolar outflow has suffered some changes in direction have also been detected. 3/2 D shock models of the line profiles are consistent with the main bipolar outflow having a small $`\pm `$ 10 angle with respect to the plane of the sky, thus making the real outflow velocity some six times larger than the observed one. The line fluxes derived from the low dispersion spectrum in the central region show that M 2-48 is nitrogen enriched; and should be classified as type I.

###### Acknowledgements.

We would like to thank the referee, You Hua Chu, for useful comments that improved our paper. López-Martín is in grateful receipt of a graduate scholarship from DGEP-UNAM (México). López-Martín and Raga acknowledge support from the CONACyT grant 34566-E. J.A. López ackowledges support from CONACYT and DGAPA (UNAM) through projects 32214-E and IN114199. Vázquez was partially supported by the CONACyT grant I32815E. Miranda and Torrelles acknowledge support from DGESIC grant PB98-0670-C02, and from Junta de Andalucía (Spain).Olguín is in grateful receipt of a graduate scholarship from DGAPA-UNAM (México). We thank R. Cook and the SPM staff for help during the observations. ——————————————————————–