text:
# The functional formulation of second-order ordinary differential equations

## 1. Introduction

The relationship between functional and ordinary differential equations was studied by many authors, e.g., M. A. Abdelkader , E. Castillo and R. Ruiz-Cobo , L. N. Eshukov . We will discuss here a functional formulation of second-order ordinary differential equations.

In Section 2, we present the functional equations (1), (2) for a continuous function of many variables. A solution for 1-dimensional linear case is found. This solution is related to Neuman’s result for linear second-order ordinary differential equations .

In Section 3, we investigate the relationship between the introduced functional equations and ordinary differential equations in general. A notion of so-called N-solvable second-order ordinary differential equation is presented. We prove that under given condition (12) the functional equations are equivalent to the family of all smooth N-solvable second-order ordinary differential equations.

Section 4 contains examples. The introduced functional equations are related to Jensen’s functional equation, to the functional equations for geodesics, and also to the functional equations for conic sections.

For the sake of simplicity, in this paper we consider smooth ordinary differential equations globally defined on the Cartesian product and globally N-solvable.

## 2. Functional equations

Let $`n`$ be a natural number, $`D=\{(\alpha ,\beta )^2\alpha =\beta \}`$ the diagonal of $`^2`$. Let $`F:\times (^2D)\times ^{2n}^n`$ be a continuous mapping. We consider functional equations

(1)
$$F(\tau ,\alpha ,\beta ,a,b)=F(\tau ,\gamma ,\delta ,F(\gamma ,\alpha ,\beta ,a,b),F(\delta ,\alpha ,\beta ,a,b)),$$

(2)
$$F(\alpha ,\alpha ,\beta ,a,b)=a,F(\beta ,\alpha ,\beta ,a,b)=b.$$

Let us start with an example.

Example 1 (linear case). Let $`n=1`$ and let $`F`$ be linear in the last two arguments. Then $`F`$ is a solution of the functional equations (1), (2) if and only if there exist functions $`x_1,x_2:`$ such that whenever $`\alpha \beta `$ the following statements hold

(3)
$$x_1(\alpha )x_2(\beta )x_1(\beta )x_2(\alpha )0,$$

(4)
$$\begin{array}{c}F(\tau ,\alpha ,\beta ,a,b)=\\ \\ \frac{(x_1(\tau )x_2(\beta )x_1(\beta )x_2(\tau ))a+(x_1(\alpha )x_2(\tau )x_1(\tau )x_2(\alpha ))b}{x_1(\alpha )x_2(\beta )x_1(\beta )x_2(\alpha )}.\end{array}$$

Let us check it. Conditions (1), (2) can be obtained from (3), (4) by direct computation. Conversely, put $`x_1(\tau )=F(\tau ,0,1,1,0)`$, $`x_2(\tau )=F(\tau ,0,1,0,1)`$. Then

(5)
$$\text{det}\left(\begin{array}{ccc}x_1(\alpha )& x_2(\alpha )& a\\ x_1(\beta )& x_2(\beta )& b\\ x_1(\tau )& x_2(\tau )& F(\tau ,\alpha ,\beta ,a,b)\end{array}\right)=0.$$

Indeed, from linearity and from (1) the third row is a linear combination of the first row and the second row with coefficients $`F(\tau ,\alpha ,\beta ,1,0)`$ and $`F(\tau ,\alpha ,\beta ,0,1)`$. Since from linearity and from (1), (2)

$$\left(\begin{array}{cc}F(0,\alpha ,\beta ,1,0)& F(0,\alpha ,\beta ,0,1)\\ F(1,\alpha ,\beta ,1,0)& F(1,\alpha ,\beta ,0,1)\end{array}\right)\left(\begin{array}{cc}x_1(\alpha )& x_2(\alpha )\\ x_1(\beta )& x_2(\beta )\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$

holds, we have the inequality (3). Finally, from (5) we can express $`F`$ in the form (4).

Remark $`1.`$ Consider a second-order linear differential equation with linearly independent solutions $`x_1,x_2`$. For this case the condition (5), where $`F`$ is the solution satisfying (2), was presented (in a more general form) by F. Neuman  in 2003. In the following chapter we will investigate the relationship between the functional equations (1), (2) and differential equations in general.

## 3. The relationship between functional and differential equations

Let $`f:\times ^n\times ^n^n`$ be a mapping. Consider a second-order equation

(6)
$$\ddot{x}=f(\tau ,x,\dot{x})$$

and Neumann conditions

(7)
$$x(\alpha )=a,x(\beta )=b,$$

where $`\alpha ,\beta `$, $`\alpha \beta `$, $`a,b^n`$.

The mapping $`F`$ will be called a dependence of the solution on Neumann conditions if and only if the following equality holds

(8)
$$F(\tau ,\alpha ,\beta ,x(\alpha ),x(\beta ))=x(\tau )$$

for any $`\alpha \beta `$, $`\tau `$ and for any complete solution $`x`$ of the equation (6). Notice that in this case $`\text{Dom}(x)=`$. We say that the dependence of the solution on Neumann conditions is smooth if and only if $`F`$ is a smooth mapping. The equation (6) is called N-solvable if and only if there exists a unique complete solution for any Neumann conditions (7). Consider the equation (6) with integral conditions

(9)
$$x(\alpha )=a,\underset{0}{\overset{1}{}}\dot{x}((1\gamma )\alpha +\gamma \beta )𝑑\gamma =v,$$

where $`\alpha ,\beta ,\gamma `$, $`a,v^n`$. The equation (6) is called I-solvable if and only if there exists a unique complete solution for any integral conditions.

Lemma 1. The smooth equation (6) is N-solvable if and only if it is I-solvable.

Proof. Let (6) be N-solvable and smooth. For $`\alpha \beta `$ we can rewrite (9) as

(10)
$$x(\alpha )=a,x(\beta )=a+v(\beta \alpha ).$$

For $`\alpha =\beta `$ we have (9) in the form

(11)
$$x(\alpha )=a,\dot{x}(\alpha )=v.$$

The existence and uniqueness of the solution of (6), (10) (respectively, (6), (11)) follows from the assumption of N-solvability (respectively, from the uniqueness of the solution of the Cauchy problem, see, e.g., ). Therefore (6) is I-solvable.

Conversely, let (6) be I-solvable. Then (6) is N-solvable. This follows from (10) immediately. $`\mathrm{}`$

Lemma 2. The dependence $`(\tau ,\alpha ,\beta ,a,v)x(\tau )`$ of the solution $`x`$ of the smooth I-solvable equation (6) on conditions (9) is smooth on an open set containing $`\times D\times ^{2n}`$, where $`D`$ is the diagonal of $`^2`$.

Proof. Let $`\tau x(\tau )`$ be the solution of (6) satisfying Cauchy conditions $`x(\alpha )=a`$, $`\dot{x}(\alpha )=u`$. The dependence of $`x`$ on $`\tau ,\alpha ,a,u`$ is smooth (see, e.g., \[6, chapter V\]). The mappings

$$h:(\alpha ,\beta ,a,v,u)\underset{0}{\overset{1}{}}\dot{x}((1\gamma )\alpha +\gamma \beta )𝑑\gamma v,$$

$$d:(\alpha ,\beta ,a,v,u)\text{det}\left(\frac{h(\alpha ,\beta ,a,v,u)}{u}\right)$$

are smooth too. It remains to check that the dependence of $`u`$ on $`\alpha ,\beta ,a,v`$ is smooth. If we put $`\beta =\alpha `$ in (9), using Cauchy conditions we get $`u=v`$. Since $`h(\alpha _0,\alpha _0,a_0,v_0,v_0)=0`$ and $`d(\alpha _0,\alpha _0,a_0,v_0,v_0)=1`$ for $`\alpha _0`$, $`a_0,v_0^n`$, we can apply \[9, Theorem 1.3.5. and Corollary 1.3.9.\] to the mapping $`h`$. From this there exists a neighbourhood $`U_1^{2n+2}\text{ of }(\alpha _0,\alpha _0,a_0,v_0)`$ and a neighbourhood $`U_2^n\text{ of }v_0`$ such that for any $`(\alpha ,\beta ,a,v)U_1`$ there is a unique $`u=u(\alpha ,\beta ,a,v)U_2`$, which satisfies $`h(\alpha ,\beta ,a,v,u)=0`$ and $`(\alpha ,\beta ,a,v)u(\alpha ,\beta ,a,v)`$ is a smooth mapping. Finally, the dependence of $`x`$ on $`\tau `$ and on parameters $`(\alpha ,\beta ,a,v)`$ from $`U_1`$ is smooth. $`\mathrm{}`$

Theorem. A mapping $`F`$ is the smooth dependence of the solution of some smooth N-solvable equation (6) on Neumann conditions if and only if the conditions (1), (2) hold and a mapping

(12)
$$\begin{array}{c}\times (^2D)\times ^{2n}(\tau ,\alpha ,\beta ,a,v)\\ F(\tau ,\alpha ,\beta ,a,a+v(\beta \alpha ))^n\end{array}$$

has a smooth extension $`S`$ to $`^{2n+3}`$.

Proof of Theorem. Assume that $`F`$ is the smooth dependence of the solution of some smooth N-solvable equation (6) on Neumann conditions (7). From (7), (8) we get $`F(\tau ,\alpha ,\beta ,a,b)=x(\tau )=F(\tau ,\gamma ,\delta ,x(\gamma ),x(\delta ))`$ for $`\alpha \beta `$, $`\gamma \delta `$, which implies (1). If we put $`\alpha =\tau `$ (respectively, $`\beta =\tau `$) in (8), using (7) we obtain first equality in (2) (respectively, second equality in (2)). We define a mapping $`S:(\tau ,\alpha ,\beta ,a,v)x(\tau )`$, where $`x`$ is the solution of the smooth N-solvable equation (6) satisfying conditions (9). The mapping $`S`$ exists for any $`(\tau ,\alpha ,\beta ,a,v)^{2n+3}`$, this follows from Lemma 1. Moreover, from (8), (10) we have $`S(\tau ,\alpha ,\beta ,a,v)=F(\tau ,\alpha ,\beta ,a,a+v(\beta \alpha ))`$ for $`\alpha \beta `$. Therefore $`S`$ is the extension of (12) to $`^{2n+3}`$. The smoothness of $`S`$ follows from the smoothness of $`F`$ and from Lemma 2.

Now we assume that (1), (2), (12) hold. With help of (12) we calculate the limit $`\delta \gamma `$ in (1) as

$$S(\tau ,\alpha ,\beta ,a,v)=S(\tau ,\gamma ,\gamma ,S(\gamma ,\alpha ,\beta ,a,v),\frac{S(\gamma ,\alpha ,\beta ,a,v)}{\gamma }).$$

Differentiating this with respect to $`\tau `$ two times and putting $`\beta =\alpha `$, $`\gamma =\tau `$ we see that $`\tau S(\tau ,\alpha ,\alpha ,a,v)`$ is the solution of (6), where

(13)
$$f(\tau ,x,\dot{x})=\frac{^2S(\tau ,\gamma ,\gamma ,x,\dot{x})}{\tau ^2}|_{\gamma =\tau }.$$

From (2), (12) we obtain $`S(\alpha ,\alpha ,\alpha ,a,v)=a`$ and

$$\frac{S(\tau ,\alpha ,\alpha ,a,v)}{\tau }|_{\tau =\alpha }=\frac{S(\tau ,\alpha ,\beta ,a,v)}{\tau }|_{\tau =\beta }|_{\beta =\alpha }=$$

$$=\left(\frac{}{\beta }\left(S(\tau ,\alpha ,\beta ,a,v)|_{\tau =\beta }\right)\frac{S(\tau ,\alpha ,\beta ,a,v)}{\beta }|_{\tau =\beta }\right)|_{\beta =\alpha }=$$

$$=\frac{}{\beta }\left(S(\tau ,\alpha ,\beta ,a,v)|_{\tau =\beta }S(\tau ,\alpha ,\beta ,a,v)|_{\tau =\alpha }\right)|_{\beta =\alpha }=$$

$$=\frac{}{\beta }\left(a+v(\beta \alpha )a\right)|_{\beta =\alpha }=v.$$

From (13) $`f`$ is smooth. Let $`x`$ be a complete solution of (6), (13). From above we can see that the solutions $`\tau S(\tau ,\alpha ,\alpha ,x(\alpha ),\dot{x}(\alpha ))\text{ and }\tau x(\tau )`$ satisfy the same Cauchy condition. Therefore

(14)
$$S(\tau ,\alpha ,\alpha ,x(\alpha ),\dot{x}(\alpha ))=x(\tau )$$

for any $`\tau `$ (see, e.g., \[6, chapter V\]). Finally, as $`\beta \alpha `$ in (1), using (7), (12), (14) we obtain $`F(\tau ,\gamma ,\delta ,x(\gamma ),x(\delta ))=x(\tau )`$, which is equivalent to (8). Hence $`F`$ is the dependence of the solution on Neumann conditions. The smoothness of $`F`$ follows from (12). The uniqueness of the complete solution of (6), (7), (13) follows from (8). Thus the equation (6), (13) is N-solvable. $`\mathrm{}`$

Remark $`2.`$ From Theorem we can see that under the condition (12) the functional equations (1), (2) are equivalent to the family of all smooth N-solvable ordinary differential equations (6). It follows from the proof of Lemma 2 that every smooth differential equation (6) is locally I-solvable, i.e., for any $`\alpha _0`$, $`a_0,u_0^n`$ there exists a neighbourhood $`U_0^{2n+1}\text{ of }(\alpha _0,a_0,u_0)`$ and a neighbourhood $`U_1^{2n+2}\text{ of }(\alpha _0,\alpha _0,a_0,u_0)`$ such that the equation (6) restricted on $`U_0`$ has a unique complete solution for any $`(\alpha ,\beta ,a,v)U_1`$. Moreover, from the proof of Lemma 1, every smooth equation (6) is locally N-solvable. Thus the functional equations (1), (2) hold locally for any smooth equation (6), but $`\text{Dom}(F)`$ is more complicated in such a case.

## 4. Examples

Example 2 (free fall). Let us consider the free fall equation $`\ddot{x}=g`$. It is easy to see, that for this equation we have

(15)
$$\begin{array}{c}F(\tau ,\alpha ,\beta ,a,b)=\frac{(\tau \beta )a+(\alpha \tau )b}{\alpha \beta }+\frac{1}{2}g(\tau \alpha )(\tau \beta )\end{array}$$

and a smooth extension (12) has the form

$$(\tau ,\alpha ,\beta ,a,v)a+v(\tau \alpha )+\frac{1}{2}g(\tau \alpha )(\tau \beta ).$$

Example 3 (geodesics). A geodesic can be viewed as a solution of a second-order differential equation. Thus this is the special case of the general situation presented above. Let $`\mathrm{\Gamma }`$ be a smooth linear connection on a manifold $`M`$. Let $`UM`$ be an open convex set and let $`a,bU`$. By $`G_{ab}`$ we denote the geodesic of $`\mathrm{\Gamma }`$ satisfying $`G_{ab}(0)=a`$, $`G_{ab}(1)=b,G_{ab}([0,1])U`$. The existence of $`U`$ and uniqueness of $`G_{ab}`$ is guaranteed by the Whitehead local convex lemma (see, e.g.,  or ). We define a mapping $`G`$ in the following way

$$G:(a,b,\rho )G_{ab}(\rho ).$$

Since there exist constant geodesics, we have

(16)
$$G(a,a,\rho )=a.$$

Let $`x`$ be a geodesic satisfying $`x(\alpha )=a\text{}x(\beta )=b`$. Since every affine reparametrization of a geodesic is a geodesic again, we have $`x(\tau )=G(a,b,\frac{\tau \alpha }{\beta \alpha })`$. If there exists a dependence of the solution on Neumann conditions, then from the previous equality and from (8) we obtain

(17)
$$F(\tau ,\alpha ,\beta ,a,b)=G(a,b,\frac{\tau \alpha }{\beta \alpha }),$$

where $`\text{Dom}(G)=^n\times ^n\times `$ and $`\text{Codom}(G)=^n`$. The following functional equations for geodesics were deduced by L. Klapka  in 2000

(18)
$$G(a,b,0)=a,G(a,b,1)=b,$$

(19)
$$G(a,b,(1\rho )\zeta +\rho \eta )=G(G(a,b,\zeta ),G(a,b,\eta ),\rho ).$$

Under the conditions (16), (17), the statements (18), (19) follow from (1), (2) for any $`\zeta ,\eta ,\rho `$.

Example 4 (Jensen’s functional equation). Consider a mapping $`G`$ from Example 3. Let $`G`$ be linear in the first two arguments. By putting $`\zeta =1`$, $`\eta =0`$ in (18), (19) we obtain $`G(a,b,1\rho )=G(b,a,\rho )`$. For any $`\rho `$ we define a mapping $`Q(\rho )`$ in the following way

$$Q(\rho ):^nbG(0,b,\rho )^n.$$

From this and from linearity we get $`G(a,b,\rho )=Q(1\rho )(a)+Q(\rho )(b)`$ and using (16) we obtain $`Q\left(\frac{1}{2}\right)(a)=\frac{1}{2}a`$. Finally, by substituting $`a=0`$, $`\rho =\frac{1}{2}`$ into (19) and using the previous equalities we get Jensen’s functional equation

$$Q\left(\frac{\zeta +\eta }{2}\right)=\frac{Q(\zeta )+Q(\eta )}{2}.$$

Example 5 (Conic sections). Let $`n=1`$. Let $`F`$ be continuous, linear inhomogeneous in the last two arguments, autonomous, i.e.,

(20)
$$F(\tau ,\alpha ,\beta ,a,b)=F(\tau +\lambda ,\alpha +\lambda ,\beta +\lambda ,a,b)$$

for any $`\lambda `$ and nondissipative, i.e., $`F(\tau ,\alpha ,\beta ,a,b)=F(\tau ,\alpha ,\beta ,a,b)`$. Under these assumptions we can find a solution of (1), (2). Let $`\alpha ,\beta ,a,b`$ be fixed. We denote $`x(\tau )=F(\tau ,\alpha ,\beta ,a,b)`$. From (1), (20) we get

$$x(\tau )=F(0,\mathrm{\Delta },2\mathrm{\Delta },x(\tau +\mathrm{\Delta }),x(\tau +2\mathrm{\Delta })),$$

$$x(\tau +3\mathrm{\Delta })=F(0,\mathrm{\Delta },2\mathrm{\Delta },x(\tau +2\mathrm{\Delta }),x(\tau +\mathrm{\Delta }))$$

for any $`\tau `$, $`\mathrm{\Delta }\{0\}`$. By using the inhomogeneous linearity and the nondissipativity we see that the difference $`x(\tau +3\mathrm{\Delta })x(\tau )`$ depends on $`x(\tau +2\mathrm{\Delta })x(\tau +\mathrm{\Delta })`$ linearly and a coefficient depends on $`\mathrm{\Delta }`$ only. Similarly, $`x(\tau +4\mathrm{\Delta })x(\tau +\mathrm{\Delta })`$ depends on $`x(\tau +3\mathrm{\Delta })x(\tau +2\mathrm{\Delta })`$ linearly and the coefficient is the same. Therefore

$$(x(\tau +4\mathrm{\Delta })x(\tau +\mathrm{\Delta }))(x(\tau +2\mathrm{\Delta })x(\tau +\mathrm{\Delta }))=$$

$$(x(\tau +3\mathrm{\Delta })x(\tau ))(x(\tau +3\mathrm{\Delta })x(\tau +2\mathrm{\Delta })).$$

This is the functional equation for characterization of conic sections presented by A. Angelesco (see \[2, section 2.5.2.\] or ) in 1922. All its continuous solutions are known and these solutions are real analytic. Differentiating this equation with respect to $`\mathrm{\Delta }`$ five times and using the limit $`\mathrm{\Delta }0`$ we get an equation

$$\dot{x}\stackrel{\dot{}\dot{}\dot{}\dot{}}{x}=\stackrel{\dot{}\dot{}\dot{}}{x}\ddot{x},$$

which is equivalent to

$$\frac{\text{d}}{\text{d}\tau }\left(\frac{\stackrel{\dot{}\dot{}\dot{}}{x}(\tau )}{\dot{x}(\tau )}\right)=0.$$

From this $`\stackrel{\dot{}\dot{}\dot{}}{x}=K\dot{x}`$ and by integrating we obtain $`\ddot{x}=Kx+g`$, where $`K,g`$. According to N-solvability $`K0`$. Putting $`K=k^2`$ we get

$$f(\tau ,x,\dot{x})=k^2x+g.$$

From (6), (8)

$$F(\tau ,\alpha ,\beta ,a,b)=$$

$$\frac{(k^2a+g)\text{sinh}(k(\tau \beta ))+(k^2b+g)\text{sinh}(k(\alpha \tau ))g\text{sinh}(k(\alpha \beta ))}{k^2\text{sinh}(k(\alpha \beta ))},$$

where $`k,g`$, $`k0`$. For $`k=0`$ we get the free fall (15).

Acknowledgment. The author would like to thank L. Klapka for his great help during the work. This research was supported by the grant MSM 192400002 from Czech Ministry of Education.


format: text

text:
# Magnetization steps in 𝐙𝐧_{𝟏-𝐱}𝐌𝐧_𝐱𝐎: Four largest exchange constants and single-ion anisotropy

## I INTRODUCTION

The most extensively studied diluted magnetic semiconductors (DMS’s) are II-VI materials ($`\mathrm{A}^{\mathrm{II}}\mathrm{B}^{\mathrm{VI}}`$, where $`\mathrm{A}=\mathrm{Zn},\mathrm{Cd}`$, $`\mathrm{B}=\mathrm{S},\mathrm{Se},\mathrm{Te}`$) in which some of the cations have been replaced by manganese.not (a) The magnetization-step (MST) method is one of the most effective techniques of measuring antiferromagnetic (AF) exchange constants in DMS’s.Shapira and Bindilatti (2002); not (b) This technique has been used to determine nearest-neighbor (NN) and distant neighbor (DN) exchange constants in several II-VI DMS’s with the zinc-blendeFoner et al. (1989); Bindilatti et al. (1998); Malarenko et al. (2000) and wurtziteShapira et al. (1989); Bindilatti et al. (1992); Cajacuri et al. (2001) structures. Relevant theoretical treatments of these exchange constants include those in Refs. Larson et al., 1988; Yu and Lee, 1995; Wei and Zunger, 1993. In addition to exchange constants, the MST method gives information about magnetic anisotropies, and about the distribution of the magnetic ions in the crystal, on a length scale of several atomic dimensions.

A new class of II-VI DMS’s based on ZnO, especially $`\mathrm{Zn}_{1x}\mathrm{Mn}_x\mathrm{O}`$, has attracted attention recently because theoretical calculations suggested the possibility of ferromagnetism above $`300\mathrm{K}`$ in p-type samples.Dietl et al. (2000) Experimental works on epitaxial thin films of $`\mathrm{Zn}_{1x}\mathrm{Mn}_x\mathrm{O}`$ gave different results: ferromagnetism was reported in Ref. Jung et al., 2002, but according to Ref. Fukumura et al., 2001 the largest exchange constant is antiferromagnetic, $`J/k_\mathrm{B}15\mathrm{K}`$. In the present work, MST’s from several $`\mathrm{Zn}_{1x}\mathrm{Mn}_x\mathrm{O}`$ single crystals, and from one powder sample, were used to determine the largest four exchange constants. The single-ion axial anisotropy parameter $`D`$ was also determined.

## II CRYSTAL STRUCTURE AND CLASSIFICATION OF NEIGHBORS

### II.1 Classification of Neighbors by Classes

The hexagonal wurtzite structure of ZnO (space group $`P6_3mc`$) is shown in Fig. 1. The cations, (open circles) form an hexagonal close-packed (hcp) structure. The cation marked by X is chosen as the “central cation.” The other cations are often classified by their distances from the central cation, i.e., nearest neighbors (NN’s), second-neighbors, etc.Larson et al. (1988) A major shortcoming of the classification by distance, is that in the (ideal) hcp structure some equidistant cations are not equivalent from symmetry point of view. Equidistant but symmetry-inequivalent cations have different isotropic exchange constants.

The classification of neighbors by symmetry, instead of distance, is discussed in Refs. Shapira and Bindilatti, 2002 and not, b. In this classification, neighbors are divided into “classes.” Neighbors of the same class have the following property: When cation sites are occupied by magnetic ions, neighbors of the same class have the same isotropic (Heisenberg) exchange interaction with the magnetic ion at the central site. The exchange constant $`J`$ is therefore the same for all neighbors of the same class. The underlying reason is that all pairs of cation sites consisting of the central site and a neighbor of a given class are related to each other by operations of the space group of the cation structure. Interactions other than isotropic exchange sometimes require distinctions between neighbors of the same class.Ladizhansky et al. (1997)

Properties of several classes of neighbors are given in Table 1. The number $`n`$, which is the same as in Ref. not, b, specifies the neighbor class. Note that $`n=1`$ and $`n=2`$ are two inequivalent classes of NN’s. As can be seen in Fig. 1, $`n=1`$ corresponds to “in plane” NN’s (i.e., NN’s which are in the same $`c`$-plane), whereas $`n=2`$ corresponds to “out of plane” NN’s. The distances in Table 1 are for the ideal hcp structure, with $`c/a=\sqrt{8/3}`$, where $`a`$ is the NN distance. A neighbor of the symmetry class $`n=3`$ corresponds to a second neighbor in the classification by distance. Neighbors of class $`n=4`$ are reached from the central cation by moving a distance $`c`$ along the $`c`$-axis. They are the closest neighbors along the hexagonal direction. Neighbors of classes $`n=5`$ and $`n=6`$ are equidistant but are inequivalent by symmetry. The remaining neighbor classes in Table 1, $`n=7`$ and $`n=8`$, are included in Fig. 1. The parameter $`z_n`$ in Table 1 is the “coordination number,” i.e., the number of neighbors of class $`n`$ surrounding the central cation.

The notation for the exchange constants $`J`$’s associated with different neighbors has been evolving, to accommodate newer classifications of these $`J`$’s. In early works, when neighbors were classified by their distance $`r`$, the notation for the $`J`$’s also was based on distance: $`J_1`$ for NN’s, $`J_2`$ for next-nearest(second)neighbors, $`J_3`$, for third neighbors, etc. For II-VI DMS’s with the zinc-blende structure (fcc cation lattice) this notation is still quite useful because each of the eight shortest distances $`r`$ is associated with a unique neighbor class $`n`$ (see footnote 115 in Ref. Shapira and Bindilatti, 2002). However, the distance-based notation is totally inadequate for DMS’s with the wurtzite structure (hcp cation structure). In the ideal hcp structure, the shortest distance $`r`$ already corresponds to two classes of NN’s, with different $`J`$’s.

An early apparent advantage of the distance-based notation followed from the predictionLarson et al. (1988) that the magnitudes (sizes) of the $`J`$’s decrease monotonically with increasing $`r`$. If true, this prediction would have made the classification by distance equivalent to a classification by size. However, later theories,Yu and Lee (1995); Wei and Zunger (1993) and recent experiments,Bindilatti et al. (1998); Hennion et al. (2002) indicate that there is no simple correspondence between size and distance.

The different notations for the $`J`$’s that are used in the present work serve different needs. The simplest notation, $`J(n)`$, associates the $`J`$’s with the neighbor classes $`n`$ listed in Table 1. For example, $`J(4)`$ is the exchange constant with a neighbor of the class $`n=4`$. The disadvantages of this notation is that neither the relevant distance $`r`$ nor the ranking by size are immediately obvious.

An alternative notation, similar to that in Ref. Larson et al., 1986, is given in the fourth row of Table 1. This notation too is based on division of neighbors into symmetry classes, but it also gives some information about the distance $`r`$. Instead of the number $`n`$, the neighbor class is specified by a combination of subscripts and superscripts. The information about distance is given by the numerical value of the subscript, which increases with increasing $`r`$. For the same subscript, the superscripts “in” and “out” are used to distinguish between equidistant but inequivalent neighbor classes. Thus, $`J_1^{\mathrm{in}}`$ and $`J_1^{\mathrm{out}}`$ are the $`J`$’s for the two classes of NN’s. A prime is added as a superscript to indicate that the distance $`r`$ is approximately, but not exactly, the same as for an unprimed exchange constant with the same subscript. For example, the exchange constants for the neighbor classes $`n=7`$ and $`n=8`$, whose distances differ by only 4%, are designated as $`J_4^{}`$ and $`J_4`$, respectively.

### II.2 Classification of Exchange Constants by Size

The preceding two notations for the $`J`$’s were both based on the division of neighbors into symmetry *classes*. This classification, however, has a serious practical drawback. Quite often the magnitude (size) of an exchange constant is measured *before* the neighbor class with which it is associated is determined. Prior to such a determination, any notation based on the neighbor class is not useful. It is then more practical to adopt a notation that is based primarily on the ranking of the $`J`$’s by size. In the present work, only the four largest $`J`$’s were measured. The chosen designations of these $`J`$’s in terms of their sizes are as follows:

1. The largest two exchange constants are labeled $`J_1`$ and $`J_1^{}`$, with $`J_1`$ chosen (arbitrarily) to be the smaller of the two. These two exchange constants are associated with the two inequivalent classes of NN’s. In the present work it has not been determined which of the two corresponds to $`J_1^{\mathrm{in}}`$ and which to $`J_1^{\mathrm{out}}`$.
2. By definition, any exchange constant $`J`$ which is not associated with either of the two classes of NN’s is a distant-neighbor (DN) exchange constant. The largest DN exchange constant (third-largest overall) is called $`J^{(2)}`$, the second-largest DN exchange constant (fourth-largest overall) is called $`J^{(3)}`$.

The assignment of $`J^{(2)}`$ and $`J^{(3)}`$ to specific neighbor classes is a major task class that will be discussed in detail.

## III EXPERIMENTAL

### III.1 Samples

Single crystals of $`\mathrm{Zn}_{1x}\mathrm{Mn}_x\mathrm{O}`$ ($`0.0056x0.030`$) were grown by chemical vapor transport using chlorine as the transporting agent. The growth temperature was $`900^\mathrm{o}\mathrm{C}`$. The Mn concentration $`x`$ was obtained using three methods:

1. From the Curie constant, obtained from a fit of the susceptibility between $`200`$ and $`300\mathrm{K}`$ to a sum of a Curie-Weiss susceptibility and a constant $`\chi _d`$ representing the lattice diamagnetism. The values $`S=5/2`$ and $`g=2.0016`$ (Ref. Dorain, 1958) for the $`\mathrm{Mn}^{2+}`$ ion were used.
2. From the apparent saturation value $`M_s`$ of the magnetization. The determination of $`M_s`$ (also known as the “technical saturation value”) is discussed in Sec. V. The relation between $`x`$ and $`M_s`$ was discussed in Refs. Shapira and Bindilatti, 2002 and Shapira, 1990, among others. This relation is based on the assumption of a random distribution of the Mn ions over the cation sites.
3. From atomic emission spectroscopy with inductively coupled plasma (ICP-AES).

Table 2 compares the results of the three methods. The good agreement indicates that the apparent saturation value $`M_s`$ is consistent with a random distribution of the Mn ions. The last column in Table 2 gives the chosen values of $`x`$ that will be used to label the various samples, namely, $`x=0.0056,0.021,0.029,\text{ and }0.030`$. These values will also be used in the data analysis.

X-ray powder diffraction data were obtained on small portions of the samples with $`x=0.0056,0.021\text{ and }0.030`$. These data were taken with a Bruker model “D8 Discover with GADDS” spectrometer, using Cu-$`\mathrm{K}_\alpha `$ radiation. All the powder diffraction patterns were in good agreement with the wurtzite structure (space group $`P6_3mc`$). No other crystallographic phase was detected.

The samples used in measurements of the magnetization $`M`$, and of the differential susceptibility $`dM/dH`$, had linear dimensions of $`2`$ to $`4\mathrm{mm}`$. The only exception was one set of pulsed-field data on a powder obtained by crushing the single crystal with $`x=0.029`$. Pulsed-field data on that single crystal were obtained before it was crushed.

### III.2 Magnetization measurements

Three types of magnetization measurements were performed:

1. The magnetization $`M`$ at $`T=20\mathrm{mK}`$ was measured in dc magnetic fields up to $`90\mathrm{kOe}`$. These data were taken with a force magnetometer operating in a plastic dilution refrigerator. The experimental techniques were described earlier.Bindilatti et al. (1996) The magnetic field $`𝐇`$ was either parallel to the $`c`$-axis, or parallel to the $`[10\overline{1}0]`$ direction (one of the directions perpendicular to the $`c`$-axis.) The dc magnetic field was produced by a NbTi superconducting magnet. Several traces of $`M\text{ vs. }H`$, in both increasing and decreasing $`H`$, were taken for each experimental configuration. All such traces were similar, and showed no hysteresis. They were averaged in order to improve the signal/noise ratio in the final result.
2. The magnetization $`M`$ was measured at $`0.65\mathrm{K}`$ in dc fields up to $`170\mathrm{kOe}`$. The samples were immersed in a liquid $`{}_{}{}^{3}\mathrm{He}`$ bath. A vibrating sample magnetometer (VSM) operating in an 18-tesla superconducting magnet ($`\mathrm{Nb}_3\mathrm{Sn}`$ wire) was used.Gratens et al. (2001) Again, no hysteresis was observed, and traces were averaged to improve the signal/noise ratio.
3. The differential susceptibility $`dM/dH`$ was measured in pulsed magnetic fields up to $`500\mathrm{kOe}`$. The experimental techniques were described earlier.Foner et al. (1989) The pulse duration was $`7.4\mathrm{ms}`$. The sample was in direct contact with a liquid $`{}_{}{}^{4}\mathrm{He}`$ bath maintained at a temperature $`T_{\mathrm{bath}}=1.5\mathrm{K}`$. However, the data showed that despite the direct contact, the sample was not in thermal equilibrium with the liquid-helium bath during the pulse. Such non-equilibrium effects were found in earlier pulsed-field experiments.Foner (1994); Shapira et al. (1999, 2001); Ajiro et al. (1998)

## IV NN EXCHANGE CONSTANTS FROM PULSED-FIELD DATA

### IV.1 NN Exchange Constants

Figure 2 shows $`dM/dH\text{ vs. }H`$ for the single crystal with $`x=0.029`$. This trace is from the field-down portion of a pulse with a maximum field of $`500\mathrm{kOe}`$. The two prominent peaks at high fields, labeled as $`1`$ and $`1^{}`$, correspond to two MST’s. The expanded view of these two peaks, shown in Fig. 3, indicates that both peaks have similar heights and widths. The maxima are at $`H_1=275\mathrm{kOe}`$ and $`H_1^{}=367\mathrm{kOe}`$. These two peaks were not resolved in the up portion of this field pulse.

Pulsed-field measurements were also performed on single crystals with $`x=0.021`$ and $`0.030`$, and on a powder obtained by crushing the single crystal with $`x=0.029`$ (after the data in Figs. 2 and 3 were obtained). In all cases the large peaks, $`1`$ and $`1^{}`$, were resolved in the down portion of the pulse. However, in the up portion of the pulse these peaks were well resolved only in the following situations: a) for the powder sample, in all field pulses, and b) for the single crystal with $`x=0.029`$ when the maximum field was $`420\mathrm{kOe}`$, as compared to $`500\mathrm{kOe}`$ for the pulse in Fig. 2. These results are explained later.

Based on the data in all samples, the two large MST’s observed in pulsed fields are at $`H_1=270\pm 8\mathrm{kOe}`$ and $`H_1^{}=362\pm 8\mathrm{kOe}`$. Because the uncertainties in these two values are correlated, the uncertainty in the difference is much smaller, i.e., $`(H_1^{}H_1)=92\pm 2\mathrm{kOe}`$. The average is $`\overline{H}_1=(H_1^{}+H_1)/2=316\pm 8\mathrm{kOe}`$.

Earlier data for other DMS’s with the wurtzite structureShapira et al. (1989); Bindilatti et al. (1992) showed that: a) each MST from NN pairs splits into a doublet, corresponding to the two inequivalent classes of NN’s, and b) the two corresponding exchange constants, $`J_1^{}`$ and $`J_1`$, are the largest. The two large MST’s at $`H_1^{}`$ and $`H_1`$ are therefore attributed to the two inequivalent NN pairs.not (c)

For $`\mathrm{Mn}^{2+}`$ pairs with intra-pair exchange constant $`J`$, the magnetic fields $`H_n`$ at the MST’s are given byShapira and Bindilatti (2002)

$$g\mu _\mathrm{B}H_n=2n|J|,$$
(1)
where $`n=1,2,\mathrm{},5`$. In the present case the calculated deviations from Eq. (1), caused by anisotropies and DN interactions, turn out to be smaller than the experimental uncertainties in $`H_1`$ and $`H_1^{}`$. Using $`n=1`$ and $`g=2.0016`$, we obtained $`J_1/k_\mathrm{B}=18.2\pm 0.5\mathrm{K}`$ and $`J_1^{}/k_\mathrm{B}=24.3\pm 0.5\mathrm{K}`$. Unfortunately, it was not possible to conclude which of the two exchange constants is $`J_1^{\mathrm{in}}`$ and which is $`J_1^{\mathrm{out}}`$. In the case of $`\mathrm{Cd}_{1x}\mathrm{Mn}_x\mathrm{Se}`$,Bindilatti et al. (1992) the smaller of the two NN exchange constants, defined here as $`J_1`$, was identified as $`J_1^{\mathrm{out}}`$, and the larger as $`J_1^{\mathrm{in}}`$. This identification was based on the effect of the Dzyaloshinkii-Moriya (DM) interaction on the widths of the MST’s, and it also agreed with Larson’s prediction.Larson (1990) In the present work the effect of the DM interaction was not apparent in the data, presumably because this interaction decreases rapidly as the atomic number of the *anion* decreases, i.e., much smaller for oxygen than for selenium.Larson and Ehrenreich (1989) Although direct evidence is lacking, based on the experimental result for $`\mathrm{Cd}_{1x}\mathrm{Mn}_x\mathrm{Se}`$ and on the theory, we speculate that in the present material too, $`J_1=J_1^{\mathrm{out}}`$ and $`J_1^{}=J_1^{\mathrm{in}}`$.

As discussed earlier,Bindilatti et al. (1992) the dominant superexchange path for both classes of NN’s, which is through the intervening anion, is the same. The difference between $`J_1^{\mathrm{in}}`$ and $`J_1^{\mathrm{out}}`$ is attributed to differences in the other exchange paths whose contribution is smaller. In the present material $`\mathrm{\Delta }\overline{J}_1=|J_1^{\mathrm{in}}J_1^{\mathrm{out}}|`$ is $`29\%`$ of the average $`\overline{J}_1/k_\mathrm{B}=21.2\pm 0.5\mathrm{K}`$. This percentage difference should be compared to $`13\%`$ in $`\mathrm{Cd}_{1x}\mathrm{Mn}_x\mathrm{S}`$, and $`15\%`$ in $`\mathrm{Cd}_{1x}\mathrm{Mn}_x\mathrm{Se}`$.Shapira et al. (1989); Bindilatti et al. (1992)

### IV.2 Other Features of the Pulsed-Field Data

Figure 2 also shows two small peaks, labeled as Q and F. Computer simulations indicate that these peaks are due in part to MST’s from quartets (tetramers). Each of these quartets consists of four spins that are coupled by some combination of $`J_1`$ and/or $`J_1^{}`$ exchange bonds. Other possible contributions to these small peaks may be due to cross-relaxation processes that can occur in the absence of equilibrium.Ajiro et al. (1998); Wernsdorfer et al. (2002); Paduan-Filho et al. (2003) The non-equilibrium behavior in pulsed fields, including cross-relaxation, will be discussed in a later publication. Very briefly, the cross relaxation processes that may contribute to Q are similar to those discussed in Ref. Paduan-Filho et al., 2003 in connection with the second-harmonic peak $`\mathrm{P}_{1/2}`$. The cross relaxation process that may contribute to F involves two NN pairs of different classes. The latter process is most rapid at $`\overline{H_1}=(H_1^{}+H_1)/2`$, where the energy separation between the two lowest levels for one class of NN pair matches that for the other class of NN pair.

The widths of peaks $`1`$ and $`1^{}`$ provide convincing evidence for other types of non-equilibrium processes, associated with a phonon-bottleneck, which restricts the sample-to-bath heat flow.Shapira et al. (1999, 2001); Nakano and Miyashita (2001); Waldmann et al. (2002); Chiorescu et al. (2000); Inagaki et al. (2003) From Fig. 3, the full width at half height of either of these two peaks is $`14\mathrm{kOe}`$. This value should be compared with a minimum equilibrium width of $`39\mathrm{kOe}`$ at $`T_{\mathrm{bath}}=1.5\mathrm{K}`$. In addition to narrow widths, peaks $`1`$ and $`1^{}`$ also show a pronounced asymmetry which is characteristic of non-equilibrium behavior resulting from a phonon bottleneck. The $`14\mathrm{kOe}`$ width at half height is the sum of $`5.5\mathrm{kOe}`$ from the rise and $`8.5\mathrm{kOe}`$ from the fall. (The data in Fig. 3 are for decreasing $`H`$, so that the rise corresponds to fields above the maximum of $`dM/dH`$, and the fall corresponds to fields below the maximum.) Both the narrowing and the asymmetry are predicted from models for the phonon bottleneck.

The phonon bottleneck also accounts for the difficulty of resolving the peaks $`1`$ and $`1^{}`$ during the up portion of the field pulse. At the beginning of the pulse the magnetocaloric effect associated with the alignment of the singles causes the sample to warm. This heating is basically the inverse of cooling by adiabatic demagnetization, except that it is not fully adiabatic. To observe the peaks $`1`$ and $`1^{}`$, much of this heat must be transferred to the bath before these peaks are reached. In the up portion of the pulse the peaks are reached earlier than in the down portion. Apparently, in the case of the single crystals an insufficient amount of heat was transferred before reaching these peaks on the way up. The peaks were therefore not resolved in the up portion of the pulse. Lowering the maximum field of the pulse leads to a slight delay of the time when the peaks are reached on the way up, which improves the chance of resolving these peaks. A more drastic change occurs when a single crystal is crushed into powder. The resulting much larger surface/volume ratio improves the sample-to-bath heat flow substantially, and the peaks are also resolved in the up portion of the pulse.

Another phenomenon seen in Fig. 2 is that (for decreasing $`H`$) a rapid rise of $`dM/dH`$ occurs below about $`50\mathrm{kOe}`$. This rise is due to MST’s from clusters in which the spins are coupled by DN exchange constants. The determination of these DN exchange constants is the main topic in the remainder of the paper.

## V DC MAGNETIZATION AT 0.65 K

Figure 4 shows magnetization curves at $`0.65\mathrm{K}`$, measured in dc magnetic fields up to $`170\mathrm{kOe}`$. As already noted, there was no hysteresis in any of the data taken in dc fields. The curves in Fig. 4 exhibit the expected behavior.Shapira and Bindilatti (2002) Above $`50\mathrm{kOe}`$ the magnetization $`M`$ shows an apparent saturation (“technical saturation”). The apparent saturation value $`M_s`$ is lower than the true saturation value $`M_0`$. The latter is expected to be reached in fields substantially above $`170\mathrm{kOe}`$. The relation between $`M_s`$ and $`x`$ was used in Sec. III as one of the three methods of determining $`x`$. An expanded view of the upper portion of each of the curves in Fig. 4 does not show any MST between $`50`$ and $`170\mathrm{kOe}`$. (The magnetization change associated with the small peak Q in Fig. 2 is estimated to be about $`0.2\%`$. This small change was not resolved in any of the dc data for $`x0.03`$). The absence of detectable MST’s in the dc data between $`50`$ and $`170\mathrm{kOe}`$ indicates that all MST’s from DN pairs occur below $`50\mathrm{kOe}`$.

A feature of Fig. 4 which is most obvious for $`x=0.030`$ is a magnetization “ramp” ending slightly above $`40\mathrm{kOe}`$. Magnetization ramps are produced by the coalescence of broadened MST’s.Shapira and Bindilatti (2002) The ramp ending just above $`40\mathrm{kOe}`$ is due to the coalescence of MST’s from clusters involving the largest DN exchange constant, defined earlier as $`J^{(2)}`$. These MST’s were not resolved at $`0.65\mathrm{K}`$, but were resolved at $`20\mathrm{mK}`$.

## VI DN EXCHANGE CONSTANTS AND SINGLE-ION ANISOTROPY FROM 0.02 K DATA

### VI.1 Overall View of the dc Magnetization at 20 mK

Figure 5 shows $`20\mathrm{mK}`$ data for $`x=0.029`$, taken both with $`𝐇𝐜`$ and $`𝐇𝐜`$. These data are normalized to the technical saturation value, $`M_s`$, and are corrected for the lattice diamagnetism. Expanded views of portions of these data are shown in Fig. 6. The main features are:

1. $`M`$ rises quickly at low fields. This initial fast rise is typical. It is mainly due to the alignment of singles ($`\mathrm{Mn}^{2+}`$ ions with no significant exchange coupling to other $`\mathrm{Mn}^{2+}`$ ions).Shapira and Bindilatti (2002)
2. The initial fast rise of $`M`$ is followed by a ramp which ends near $`5`$ or $`6\mathrm{kOe}`$, depending on field direction. This ramp is shown more clearly in the expanded view of Fig. 6(a).
3. A second ramp, smaller in height but spread over a larger field interval, follows the first ramp. The second ramp ends near $`40\mathrm{kOe}`$. Just above this field, $`M`$ reaches technical saturation.
4. The magnetization depends on field direction. This anisotropic behavior is more obvious at low $`H`$.

Each of the two ramps is due to the coalescence of broadened MST’s.Shapira and Bindilatti (2002) A well defined ramp usually corresponds to one series of MST’s. The high-field end of such a ramp is near the last MST from this series, and it can be used to estimate the relevant exchange constant $`J`$. A more accurate value for $`J`$ can be obtained if the MST’s on the ramp are well resolved.

The end of the second ramp, near $`40\mathrm{kOe}`$, leads to the estimate $`J^{(2)}/k_\mathrm{B}0.5\mathrm{K}`$ for the largest DN exchange constant. The end of the first ramp, near $`6\mathrm{kOe}`$, gives $`J^{(3)}/k_\mathrm{B}0.08\mathrm{K}`$. To improve on these rough estimates it is necessary to examine the MST’s which give rise to each ramp, taking into account the relevant weak anisotropic interactions.

### VI.2 Cluster Models and Anisotropies

Cluster models play a key role in analysis of MST’s involving DN exchange interactions.Shapira and Bindilatti (2002) Detailed information about cluster models and their statistical properties is given in Ref. not, b. In the present work, only the largest two DN exchange constants, $`J^{(2)}`$ and $`J^{(3)}`$, were determined. The cluster models that were used in the data analysis depended on the field range of the data:

1. For fields below $`8\mathrm{kOe}`$, both $`J^{(2)}`$ and $`J^{(3)}`$ are important. These exchange constants correspond to two classes of DN’s, but the identity of neither of these two classes was known at the beginning of the analysis. Therefore, the cluster models used for this field range included: a) the two classes of NN’s, and b) any possible two classes of DN’s selected from the six classes listed in Table 1.
2. Because the magnetization above $`10\mathrm{kOe}`$ is hardly affected by $`J^{(3)}`$, the cluster models for the field range from $`10`$ to $`50\mathrm{kOe}`$ involved only the one class of a DN associated with $`J^{(2)}`$, and the two NN classes associated with $`J_1`$ and $`J_1^{}`$. All six possible choices for the class for the DN were tried.
3. Analysis of the data between $`50\mathrm{kOe}`$ and $`90\mathrm{kOe}`$ was based primarily on a model which included only the two NN exchange constants. Simulations of the magnetization curves showed that in this field range the effects of the DN exchange constants $`J^{(2)}`$ and $`J^{(3)}`$ were very small.

In addition to the exchange interactions, two anisotropies were also included in the analysis: single-ion anisotropy, and dipole-dipole (dd) anisotropy. The effects of these anisotropies on the MST’s (or ramp) from pairs become more pronounced as the magnitude $`|J|`$ of the relevant intra-pair exchange constant decreases. That is, the effects caused by the anisotropies are very small for the MST’s originating from NN pairs of either class; small but easily detected for the MST’s (ramp) from pairs involving $`J^{(2)}`$; and very pronounced for the MST’s (ramp) from pairs involving $`J^{(3)}`$.

The single-ion anisotropy can be described by the hamiltonian

$$DS_z^2+\frac{a}{6}(S_x^4+S_y^4+S_z^4),$$

where the $`z`$ direction is along the $`c`$-axis. From Electronic Paramagnetic Resonance (EPR) at $`77\mathrm{K}`$, $`D/k_\mathrm{B}=31\mathrm{mK}`$ and $`a/k_\mathrm{B}=0.3\mathrm{mK}`$.Dorain (1958) Because the term involving $`a`$ is relatively small, it was neglected in the present work. Values of the dd anisotropy constant for different classes of pairs are given in Table 1. They are based on a simple model in which the two spins in a pair are represented by two points separated by $`r_n`$.

### VI.3 Results at 20 mK and their Analysis

#### VI.3.1 Objectives

The two objectives of the experiments at $`20\mathrm{mK}`$ were: 1) to determine the values of $`J^{(2)}`$, $`J^{(3)}`$, and $`D`$, and 2) to identify the DN classes associated with $`J^{(2)}`$ and $`J^{(3)}`$. The analysis consisted of a number of steps taken in sequence. In what follows, each step in this sequence is outlined, and the results are summarized. The main assumption in the analysis is that the Mn ions are randomly distributed over the cation sites.

#### VI.3.2 Procedures for Identifying DN classes

The general procedure of associating DN constants of known magnitudes with different possible classes of DN’s involves comparisons between experimental magnetization curves and computer simulations.Shapira and Bindilatti (2002); Bindilatti et al. (1998) Separate simulations are carried out for all competing possibilities for the DN classes. The six DN classes listed in Table 1 lead to $`6\times 5=30`$ possibilities for the two DN classes associated with $`J^{(2)}`$ and $`J^{(3)}`$. (Interchanging the order of the two classes of DN’s leads to a new possibility.)

In the present work these laborious simulations were postponed, because a preliminary identification of the DN classes was possible based on a simpler, and more physical, procedure. Two favorable circumstances permitted the simpler procedure: a large difference in the magnitudes of $`J^{(2)}`$ and $`J^{(3)}`$, and the availability of data for a sample with a very small $`x`$.

For the lowest Mn concentration, $`x=0.0056`$, the ramp which ends near $`40\mathrm{kOe}`$ is due primarily to $`J^{(2)}`$ pairs. To a good approximation, the total magnetization rise $`\mathrm{\Delta }M^{(2)}`$ associated with this ramp is therefore related to the fraction of Mn ions which are in such pairs. This fraction can be calculated for any possible choice of the DN class associated with $`J^{(2)}`$, so that $`\mathrm{\Delta }M^{(2)}`$ for any choice can be calculated and compared with experiment. Inclusion of $`J^{(2)}`$ triplets, in addition to $`J^{(2)}`$ pairs, in the calculation of $`\mathrm{\Delta }M^{(2)}`$ makes the comparison even more reliable.

Using the statistical tables in Ref. not, b, the contribution of pairs and triplets to $`\mathrm{\Delta }M^{(2)}`$ was calculated for each of the six possible choices of the DN associated with $`J^{(2)}`$. Only one of these choices agreed with the experimental value of $`\mathrm{\Delta }M^{(2)}`$. Thus, a unique identification of the cluster class associated with $`J^{(2)}`$ was achieved. A similar procedure was also used for a preliminary identification of the DN responsible for $`J^{(3)}`$, but in that case two possible DN classes gave good agreement with experiment, so that a unique choice could not be made. The full numerical simulations, *for all the samples*, carried out at the end of the analysis, confirmed the preliminary identifications of the DN classes on the basis of the data for $`x=0.0056`$.

#### VI.3.3 DN Class Responsible for $`J^{(2)}`$

The procedure of extracting the experimental ratio $`\mathrm{\Delta }m=\mathrm{\Delta }M/M_s`$ from the data, for the ramps associated with $`J^{(2)}`$ and $`J^{(3)}`$, is illustrated in Fig. 7. For the higher-field ramp, associated with $`J^{(2)}`$, it gives $`\mathrm{\Delta }M^{(2)}/M_s=1.3\%`$ for $`x=0.0056`$. The theoretical value of $`\mathrm{\Delta }M^{(2)}/M_s`$ is obtained by multiplying $`(\mathrm{\Delta }M^{(2)}/M_0)`$ by $`(M_0/M_s)`$. The first ratio depends on the DN class $`n`$ that corresponds to $`J^{(2)}`$. For each possible $`n`$, this ratio was obtained from probability tables for pairs and triplets that involve only the exchange constant for class $`n`$ \[not, b\]. The probabilities are based on a cluster model that includes only the two classes of NN’s and the DN class $`n`$. The second ratio, $`(M_0/M_s)`$, was calculated from the so-called NN cluster model,Shapira and Bindilatti (2002) in which only the two classes of NN’s are included.

For small $`x`$ the calculated ratio of $`\mathrm{\Delta }M^{(2)}/M_s`$ increases as the coordination number $`z_n`$ increases. This dependence on $`z_n`$ is the key for identifying the DN class that corresponds to $`J^{(2)}`$. For $`x=0.0056`$ the calculated values of $`\mathrm{\Delta }M^{(2)}/M_s`$ are: $`1.39\%`$ for the only DN class with $`z_n=2`$; between $`3.1\%`$ and $`4.2\%`$ for the three possible classes with $`z_n=6`$; and $`6.1\%`$ for the two possible classes with $`z_n=12`$. The experimental value $`1.3\%`$ agrees only with $`z_n=2`$, which corresponds to $`n=4`$ in Table 1. In the ideal wurtzite structure such a DN is reached from the “central” cation by moving a distance $`r_4=c=\sqrt{8/3}a`$ along the $`c`$-axis (see Fig. 1). In Table 1 the exchange constant $`J(4)`$ for $`n=4`$ is designated as $`J_3^{}`$. It is noteworthy that this largest DN exchange constant is not for the DN that is closest to the central cation. The closest DN is of class $`n=3`$, at a distance $`r_3=\sqrt{2}a`$. Thus, the magnitudes of the exchange constants do not decrease monotonically with distance. This non-monotonic dependence on distance was predicted by some theories,Yu and Lee (1995); Wei and Zunger (1993) and has been observed experimentally earlier.Bindilatti et al. (1998); Hennion et al. (2002)

#### VI.3.4 Values of $`J_3^{}`$ and $`D`$

The values for $`J_3^{}`$ and $`D`$ were obtained from analysis of the well resolved MST’s on the ramp associated with $`J^{(2)}=J_3^{}`$ (see part (b) of Fig. 6). These MST’s stand out more clearly in the derivative $`dm/dH`$ of the normalized magnetization $`mM/M_s`$, shown in Fig. 8. For three of the four samples in this figure the data are for both $`𝐇𝐜`$ and $`𝐇𝐜`$. The dependence of the fields at the MST’s on field direction is caused by the anisotropy. Because the $`c`$-axis is an easy axis ($`D<0`$), the spins in the pairs are aligned more quickly for $`𝐇𝐜`$ than for $`𝐇𝐜`$. The faster alignment for $`𝐇𝐜`$ is more apparent in part (b) of Fig. 6.

All curves in Fig. 8 show the last four MST’s from $`J_3^{}`$ pairs. For $`x=0.0021`$ the first MST is also seen, but only as a “shoulder” on the fast drop of $`dm/dH`$ at low fields. The fields at the last four MST’s are all above $`15\mathrm{kOe}`$. Values of both $`J_3^{}`$ and $`D`$ were obtained from an analysis of these fields, based on a pair hamiltonian which included the exchange interaction due to $`J_3^{}`$, the uniaxial anisotropy governed by $`D`$, and the dd-interaction. The latter was calculated using the value in Table 1 for $`n=4`$.

In the first step of the analysis, an approximate value for $`J_3^{}`$ was obtained from the field at the last (fifth) MST, and the EPR value of $`D`$ (Ref. Dorain, 1958) was adopted as the initial value. Later, the predicted fields at the last four MST’s, for both field directions, were calculated for many sets of ($`J_3^{}`$, $`D`$). The best match with the experimental values gave: $`J_3^{}/k_\mathrm{B}=0.543\pm 0.005\mathrm{K}`$, and $`D/k_\mathrm{B}=0.039\pm 0.008\mathrm{K}`$. The presumably more accurate EPR value for $`D/k_\mathrm{B}`$ is $`0.031\mathrm{K}`$. The difference may be related to our use of the dd-interaction constant given in Table 1. As noted, this constant was obtained from a simple model in which the spins in a pair are represented by two points separated by $`r_4`$.

#### VI.3.5 Possible DN Classes Associated with $`J^{(3)}`$

Possible assignments of the DN class which corresponds to $`J^{(3)}`$ were made on the basis of the magnitude $`\mathrm{\Delta }m^{(3)}`$ in Fig. 7. The experimental results give $`\mathrm{\Delta }m^{(3)}=\mathrm{\Delta }M^{(3)}/\mathrm{\Delta }M_s7.3\%`$ for $`x=0.0056`$. This $`\mathrm{\Delta }m^{(3)}`$ was attributed to the combined magnetization rise from pairs and triplets involving this exchange constant. The inclusion of triplets in the theoretical calculation of $`\mathrm{\Delta }m^{(3)}`$ was more important than in the calculation of $`\mathrm{\Delta }m^{(2)}`$, because the triplets/pairs population ratio was higher. The triplets/pairs ratio increases with $`z_n`$. The lowest possible coordination number, $`z_n=2`$, is for the neighbor class $`n=4`$ associated with $`J^{(2)}`$ and $`\mathrm{\Delta }m^{(2)}`$.

Theoretical values of $`\mathrm{\Delta }m^{(3)}`$ were obtained from cluster models which included the three largest exchange constants ($`J_1`$, $`J_1^{}`$, $`J_3^{}`$) and any one of the DN classes with either $`z_n=6`$ or $`z_n=12`$. The calculated values are: approximately $`3.1\%`$ for all three DN classes with $`z_n=6`$, and approximately $`6.0\%`$ for the two DN classes with $`z_n=12`$. On this basis, $`z_n`$ is equal to $`12`$, which leads to two possible neighbor classes: $`n=6`$ or $`n=7`$ (see Table 1). That is, $`J^{(3)}`$ is either $`J_3^{\mathrm{out}}`$ or $`J_4^{}`$. Another possibility (unlikely, but cannot be ruled out entirely) is that $`J^{(3)}`$ is associated with two different classes, both with $`z_n=6`$, that just happen to have very nearly equal exchange constants.

#### VI.3.6 Value of $`J^{(3)}`$

Figure 9 displays $`dm/dH`$ data for both $`𝐇𝐜`$ and $`𝐇𝐜`$, in the field range relevant for the analysis of $`J^{(3)}`$. Due to the smallness of $`J^{(3)}`$, the effects of the anisotropy are important, as can be judged from the strong dependence of the results on the direction of $`𝐇`$. The $`dm/dH`$ traces for $`𝐇𝐜`$, in Fig. 9(a), do not show the regular sequence of peaks observed in Fig. 8 for the MST’s from $`J_3^{}`$-pairs. However, all traces in Fig. 9(a) show a peak slightly above $`2\mathrm{kOe}`$. The best resolution is for the sample with $`x=0.021`$. In this sample, two additional peaks, just below and just above $`4\mathrm{kOe}`$ are also resolved.

The results in Fig. 9 were compared with simulations based on a pair hamiltonian which included $`J^{(3)}`$ and the two anisotropies. The simulations used $`D/k_\mathrm{B}=39\mathrm{mK}`$, as determined in Sec. VI.3.4, and the dd-interaction constant taken from Table 1. Because the latter constant depends on the DN class, simulations were carried out for the two possible neighbor classes with $`z_n=12`$, i.e., $`n=6`$ with $`J_3^{\mathrm{out}}`$, and $`n=7`$ with $`J_4^{}`$.

For $`𝐇𝐜`$, pairs involving neighbors of either class have different orientations relative to $`𝐇`$. The different orientations lead to different energies from the dd-anisotropy, and they give rise to MST’s at slightly different fields. Therefore, in the simulations for the perpendicular field direction, $`𝐇[10\overline{1}0]`$, pairs involving DN’s of either of the two possible classes were divided into groups Pairs in different groups had different orientations relative to $`𝐇`$. There were three such groups for the DN class $`n=6`$, with $`J_3^{\mathrm{out}}`$, and two groups for $`n=7`$, with $`J_4^{}`$. The simulated magnetization curve was obtained by adding the results from all the groups.

The simulations were carried out using the actual temperature, $`T=20\mathrm{mK}`$. Non-thermal broadening mechanism, such as local strains that give rise to a spread of the exchange and anisotropy interactions,Rubo et al. (1997) were ignored. In these simulations the values $`D/k_\mathrm{B}=39\mathrm{mK}`$, and the dd-interaction constant from Table 1, were kept fixed. The value of $`J^{(3)}`$ was adjusted to obtain the best match with experiment for both field orientations. The final $`dm/dH`$ simulations are shown in Figs. 9 (a) and (b) as the dashed and dotted lines. Obviously, the results of the simulations depend on field direction. For $`𝐇𝐜`$, the overall structure of the experimental $`dm/dH`$ traces is reproduced by the simulations. However, because some sources of line broadening were neglected in the simulations, the detailed structure is better resolved in the simulations than in the experimental curves. For $`𝐇𝐜`$ the individual MST’s are resolved in the simulations, even after the different orientations of the pairs relative to $`H`$ are included. Experimentally, however, the individual MST’s are not resolved for $`𝐇𝐜`$. This difference is attributed, again, to the neglect of some broadening mechanisms in the simulations. A crude way of accounting for the neglected broadening mechanisms is to replace the actual temperature $`T`$ in the simulations by a higher effective temperature $`T_{\mathrm{eff}}`$. The minimum $`T_{\mathrm{eff}}`$ which leads to unresolved MST’s for $`𝐇𝐜`$ is $`65\mathrm{mK}`$.

Some features of the experimental data in Figs. 9(a) and 9(b) are sensitive to the magnitude of $`J^{(3)}`$. For $`𝐇𝐜`$ these features include the field at the most prominent peak, and the field at the peak associated with the last MST. These peaks stand out most clearly in the trace for $`x=0.021`$. For $`𝐇𝐜`$ the field at the rapid drop of $`dm/dH`$, which is at the end of the ramp associated with this series of MST’s, is sensitive to the value of $`J^{(3)}`$. The value of $`J^{(3)}`$ was determined from comparisons of these experimentally observed features with simulations that used different values of $`J^{(3)}`$. Assuming that $`J^{(3)}`$ is $`J_3^{\mathrm{out}}`$, the results gave $`J^{(3)}/k_\mathrm{B}=0.074\pm 0.005\mathrm{K}`$. The alternative, $`J^{(3)}=J_4^{}`$ gave $`0.082\pm 0.005\mathrm{K}`$. The first choice gives a slightly better agreement with the data, but in our view the evidence is insufficient for concluding that the DN class is definitely $`J_3^{\mathrm{out}}`$.

#### VI.3.7 Simulations with both DN exchange constants

Simulations of the magnetization curves, in fields up to $`60\mathrm{kOe}`$, were carried out in order to confirm the preliminary identifications of the DN classes corresponding to $`J^{(2)}`$ and $`J^{(3)}`$. The simulations were for all samples, in contrast with the preliminary analysis that was carried out only for the sample with the lowest Mn concentration, $`x=0.0056`$. The simulations used cluster models which included the two NN exchange constants ($`J_1`$ and $`J_1^{}`$), $`J^{(2)}`$ and $`J^{(3)}`$, all having the values quoted above. Each cluster model was based on a specific choice of the two DN classes associated with $`J^{(2)}`$ and $`J^{(3)}`$. As noted earlier, there are $`30`$ possible such choices. The results of the simulations were not sensitive to any change of the DN class $`n`$ provided that the coordination number $`z_n`$ did not change.

Because all anisotropies were neglected in the simulations, the comparison was made with “isotropic” magnetization curves obtained from the experimental data using the relation

$$M_{\mathrm{isot}}=(M_{}+2M_{})/3,$$
(2)
where $`M_{}`$ is for $`𝐇𝐜`$, and $`M_{}`$ is for $`𝐇𝐜`$. The widths of the MST’s exhibited by $`M_{\mathrm{isot}}`$ are larger than the thermal width at $`20\mathrm{mK}`$ for several reasons: First, $`M_{\mathrm{isot}}`$ is an average over different groups of pairs with MST’s at slightly different fields. Second, DN exchange constants smaller than $`J^{(3)}`$ were neglected. Third, variations of the $`J`$’s, and of the anisotropies, caused by local strainsRubo et al. (1997) were not included. To match the widths of the MST’s on the ramp which ends at $`40\mathrm{kOe}`$, the simulations were carried out using an effective temperature $`T_{\mathrm{eff}}=100\mathrm{mK}`$ instead of $`20\mathrm{mK}`$. This change has no effect on the identification of the DN’s responsible for $`J^{(2)}`$ and $`J^{(3)}`$; the only effect is to smooth the curves. Simulations using the actual temperature, $`20\mathrm{mK}`$, lead to the same identifications of the DN classes.

Figure 10 shows simulations for $`x=0.029`$. The “experimental” curve represents $`M_{\mathrm{isot}}`$. Each simulation is labeled as $`n=(i,j,k,l)`$. This designation means that the four exchange constants $`J_1`$, $`J_1^{}`$, $`J^{(2)}`$, $`J^{(3)}`$ that were used in the simulation, correspond, respectively, to the neighbor classes $`n=i,j,k,l`$ in Table 1. All the simulations assumed that the first two neighbor classes $`(i,j)`$ were those of the NN’s, that is, $`(i,j)`$ = $`(1,2)`$ or $`(2,1)`$. (These two alternative choices are related to each other by an interchange of the NN classes assigned to $`J_1`$ and $`J_1^{}`$.) Both choices lead to nearly the same curves in this field range. The best agreement is with the simulation $`n=(1,2,4,6)`$ in which $`J^{(2)}=J_3^{}`$ and $`J^{(3)}=J_3^{\mathrm{out}}`$. This corresponds to one of the two possibilities which were identified earlier. The other possibility $`n=(1,2,4,7)`$, not shown in Fig. 10, leads to a very similar curve, so that a definite unique choice of the DN class for $`J^{(3)}`$ is not possible.

The same conclusions concerning the neighbor classes associated with $`J^{(2)}`$ and $`J^{(3)}`$ were reached from comparisons of the data for the other samples ($`x=0.0056`$, $`x=0.021`$ and $`x=0.030`$) with simulations. Thus, the earlier identifications of the DN classes are confirmed. The good agreement between the data and the simulations also lends support to a random Mn distribution in the studied samples, which is the main assumption in the simulations.

###### Acknowledgements.

The work in Brazil was supported by FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo, Brazil) under contract number 99/10359–7. The work in Poland was supported by grant PBZ-KBN-044/P03/2001. VB and NFOJ acknowledge support from CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brazil). Travel funds for YS were provided by FAPESP.


text:
# Entanglement interferometry for precision measurement of atomic scattering properties

## Abstract

We report on a two-particle matter wave interferometer realized with pairs of trapped 87Rb atoms. Each pair of atoms is confined at a single site of an optical lattice potential. The interferometer is realized by first creating a coherent spin-mixture of the two atoms and then tuning the inter-state scattering length via a Feshbach resonance. The selective change of the inter-state scattering length leads to an entanglement dynamics of the two-particle state that can be detected in a Ramsey interference experiment. This entanglement dynamics is employed for a precision measurement of atomic interaction parameters. Furthermore, the interferometer allows to separate lattice sites with one or two atoms in a non-destructive way.

The controlled creation of entanglement is one of the most subtle and challenging tasks in modern quantum mechanics, with both wide reaching practical and fundamental implications. In neutral atom based systems, significant progress has been made during recent years in the generation of large spin-squeezed samples of atomic gases Julsgaard et al. (2001) or the controlled creation of Greenberger-Horne-Zeilinger (GHZ) states Greenberger et al. (1990) in cavity QED systems Rauschenbeutel et al. (2000). In addition, it has been recognized early on that in binary spin-mixtures of Bose-Einstein condensed quantum gases (with spin states $`|0`$ and $`|1`$) a large amount of entanglement could be created by controlling the difference in interaction strengths $`\chi =1/2(U_{00}+U_{11}2U_{01})`$ between the particles in different spin states Sørenson et al. (2001); Jaksch et al. (1999); Helmerson and You (2001); You (2003); Sørenson and Mølmer (1999). Here $`U_{ij}`$ denotes the interaction matrix element between atoms in spin states $`i`$ and $`j`$. Such control can either be achieved by moving atoms on different sites in spin-dependent optical lattice potentials Jaksch et al. (1999); Sørenson and Mølmer (1999); Brennen et al. (1999); Mandel et al. (2003a, b) or by tuning the scattering lengths, such that $`\chi 0`$ Sørenson et al. (2001); You (2003); Micheli et al. (2003). In the latter case, the simple creation of a coherent spin-mixture, e.g. by an initial $`\pi /2`$-pulse, followed by a subsequent evolution of the spin system, would automatically lead to highly spin-squeezed or entangled $`N`$-Particle GHZ-like states. Here we demonstrate such entanglement dynamics with pairs of atoms trapped in the ground state of a potential well in an optical lattice. Such pairs form a unique and highly controllable model system to study interactions between two particles. By using a recently predicted inter-state Feshbach resonance in 87Rb van Kempen et al. (2002) we are able to control $`\chi `$ and observe the ensuing entanglement evolution in a Ramsey type experiment. We show that this dynamical evolution of the atom pairs into entangled and disentangled states can be used to obtain precise information on the scattering properties of the systems. The entanglement interferometer makes it furthermore possible to separate singly occupied lattice sites from doubly occupied sites in a non-destructive way.

The Ramsey interferometer sequence used in the experiment consists of two $`\pi /2`$-pulses with pulse separation $`t_{\text{hold}}`$, where the phase $`\alpha `$ of the last pulse can be varied. Between the two $`\pi /2`$-pulses the interaction behavior $`\chi `$ of the atoms is modified by using a Feshbach resonance occurring between atoms in hyperfine states $`|0|F=1,m_F=+1`$ and $`|1|F=2,m_F=1`$ van Kempen et al. (2002); Erhard et al. (2003). Here $`F`$ and $`m_F`$ denote the total angular momentum and its projection, respectively. In the case of a single isolated atom, initially in the state $`\psi _1^i=|0`$, the effect of a changing interaction behavior has no consequence for the particle since it does not interact with other particles. After applying the experimental sequence shown in Fig. 1 the probability of finding the atom in state $`|1`$ is simply given by $`P(\psi _1^f,\alpha )=\frac{1}{2}\left(1\mathrm{cos}\alpha \right)`$, which describes the usual Ramsey fringe without decoherence.

In the case of two particles, the change of interaction parameters leads to an additional entangling-disentangling dynamics which is markedly different. Let us consider two atoms in the internal state $`\psi _2^i=|0|0|00`$. First a $`\pi /2`$-pulse places the atoms in a coherent superposition of the two internal states $`\psi _2=\frac{1}{2}\left(|00|01|10+|11\right)`$. After a time $`t`$ each two-particle state obtains a phase factor $`e^{i\varphi _{ij}}`$ due to interactions, where $`\varphi _{ij}=\frac{U_{ij}}{\mathrm{}}t`$ is the collisional phase shift, with $`U_{ij}=(4\pi \mathrm{}^2a_{ij})/m\times d^3x|\psi _i|^2|\psi _j|^2`$ being the on-site interaction matrix element. Here $`a_{ij}`$ represents the elastic scattering length between particles in states $`i`$ and $`j`$, $`\psi _{i(j)}`$ is the ground state wave function of an atom in spin state $`i(j)`$ and $`m`$ the mass of a single atom. For $`t=t_{\text{hold}}/2`$ the two particle state then evolves into $`\psi _2=\frac{1}{2}(e^{i\frac{\varphi _{00}}{2}}|00e^{i\frac{\varphi _{01}}{2}}|01e^{i\frac{\varphi _{01}}{2}}|10+e^{i\frac{\varphi _{11}}{2}}|11)`$, where $`\varphi _{01}=\varphi _{10}`$. After a spin echo $`\pi `$-pulse, a further interaction time $`t_{\text{hold}}/2`$ and a last $`\pi /2`$-pulse the final state reads $`\psi _2^f=\frac{1}{2}\left\{c_c^+|00c_s(|01+|10)+c_c^{}|11\right\}`$, where $`c_c^\pm e^{i\alpha }(\mathrm{cos}\alpha \pm e^{i\varphi _\chi })`$, $`c_si\mathrm{sin}\alpha `$ and $`\varphi _\chi (\varphi _{00}+\varphi _{11}2\varphi _{10})/2`$. The probability of finding an atom in state $`|1`$ can then be expressed by $`P(\psi _2^f,\alpha ,\varphi _\chi )=\frac{1}{2}(1\mathrm{cos}\alpha \mathrm{cos}\varphi _\chi )`$, which is modulated in amplitude compared to the case of single atoms.

Four main cases illustrate the dynamics of the two-particle system: (i) For $`\varphi _\chi =0`$ the Ramsey fringe $`P(\psi _2^f,\alpha ,\varphi _\chi =0)=\frac{1}{2}\left(1\mathrm{cos}\alpha \right)`$ is identical to the fringe of an isolated particle as shown above. (ii) If the interactions lead to a phase difference of $`\varphi _\chi =\pi /2`$, the final state $`\psi _2^f(\varphi _\chi =\pi /2)`$ is a maximally entangled Bell-like state. For such a state the corresponding Ramsey fringe $`P(\psi _2^f,\alpha ,\varphi _\chi =\pi /2)=\frac{1}{2}`$ does not exhibit any modulation Mandel et al. (2003b). (iii) When the phase difference is increased to $`\varphi _\chi =\pi `$, the system is disentangled again, and the corresponding Ramsey fringe $`P(\psi _2^f,\alpha ,\varphi _\chi =\pi )=\frac{1}{2}\left(1+\mathrm{cos}\alpha \right)`$ is phase shifted by $`\pi `$ with respect to the case of a single particle. It should be noted that for this interaction phase the state vectors of isolated atoms and interacting atom pairs are orthogonal to each other. Therefore, by choosing a specific single particle phase, either species can be transferred into the $`|F=2`$ state and removed by a subsequent resonant laser pulse. The remaining atoms would form a pure lattice of either single atoms or atom pairs which again can evolve to Bell-like pairs. For even larger phase differences, the system entangles and disentangles again, until for (iv) $`\varphi _\chi =2\pi `$ the system exhibits a Ramsey fringe which is in phase with the fringe of a single atom. In a system containing $`N_1`$ isolated single atoms and $`N_2/2`$ isolated pairs of atoms, the total fringe will be a superposition of the two distinct fringes and will have a visibility according to:

$$V(\chi )=V_0e^{\frac{t}{\tau _1}}\left\{(1n_2)+n_2e^{\frac{t}{\tau _2}}\mathrm{cos}\varphi _\chi \right\},$$
(1)
where we have included decoherence and two-body losses with time constants $`\tau _1`$ and $`\tau _2`$, respectively, $`V_0`$ is a finite initial visibility 111Here the visibility is defined as $`V=(P_{max}^{tot}P_{min}^{tot})/(P_{max}^{tot}+P_{min}^{tot})`$. and $`n_i=N_i/(N_1+N_2)`$, $`i=1,2`$ with $`n_1+n_2=1`$. Whereas the contribution of isolated atoms to the fringe visibility $`V`$ remains unaffected under a change of $`\chi `$, the total fringe signal shows a dynamics with the same periodicity as $`P(\psi _2^f)`$. For zero phase difference $`\varphi _\chi =0`$ we expect to measure a Ramsey fringe with high visibility $`V`$. The visibility decreases for increasing $`\varphi _\chi `$ and reaches a minimum for $`\varphi _\chi =\pi `$ where the fringes from single atoms and from atom pairs are out of phase and partially compensate each other in the total signal. For larger interaction phase differences, the total visibility increases, until it shows a maximum for $`\varphi _\chi =2\pi `$, where the two Ramsey fringes are completely in phase again. It should be noted that the interaction time $`t_R`$ after which $`\varphi _\chi =2\pi `$ depends on the difference in the interaction matrix elements $`U_{ij}`$, and for constant overlap of the wave functions on the elastic scattering length difference $`\mathrm{\Delta }a_{s,\chi }\frac{1}{2}(a_{00}+a_{11}2a_{10})`$.

The experimental setup is similar to our previous work Greiner et al. (2002a, b). We start with a BEC of up to $`3\times 10^5`$ 87Rb atoms trapped in the hyperfine ground state $`|F=1,m_F=1`$. We load the BEC into a pure two or three dimensional optical lattice potential formed by three mutually orthogonal standing waves of far detuned light. The wavelengths used for the different standing waves are 829 nm ($`y`$ and $`z`$ axes) and 853 nm ($`x`$ direction), with trapping frequencies at each lattice site of $`\omega _x=2\pi \times 33`$ kHz, $`\omega _y=2\pi \times 43`$ kHz, and $`\omega _z=2\pi \times 41`$ kHz. In order to preserve spin polarization of the atoms in the optical trap, we maintain a 1 G magnetic offset field along the $`x`$-direction. The atoms are prepared in the Feshbach resonance sensitive spin mixture by transferring the population from the $`|F=1,m_F=1`$ state into the $`|0|F=1,m_F=+1`$ hyperfine level via a radio frequency (RF) Landau-Zener sweep. We then increase the magnetic field to 8.63 G. By applying a microwave field around 6.8 GHz and RF radiation around 6 MHz, we are able to coherently couple the two internal states $`|0`$ and $`|1`$ with a two-photon transition similar to Hall et al. (1998).

In order to locate the position of the Feshbach resonance through enhanced atomic losses, we load the BEC into a two dimensional optical lattice potential. The atomic density is thereby strongly increased compared to a simple dipole trap, thus loss processes occur with higher probability. The magnetic field is subsequently increased to different values within 10 $`\mu `$s. After holding the atoms for 1 ms at a specific magnetic field, we switch off all trapping potentials and magnetic fields to measure the remaining total atom number in a time-of-flight (TOF) measurement (see Fig. 2). The magnetic field has been calibrated by measuring the frequency of the $`|F=1,m_F=1|F=2,m_F=2`$ microwave transition at different magnetic field values and employing the Breit-Rabi formula to determine the actual field strength. Due to background magnetic field fluctuations, the magnetic field calibration has an uncertainty of 3 mG, and noise of the magnetic field creating current source introduces an additional uncertainty of 2 mG. The measured position of 9.121(5) G of the resonance agrees well with the predicted value of 9.123 G within our measurement uncertainty 222The inter-state Feshbach resonance has recently also been detected through atom loss measurements by Erhard et al. Erhard et al. (2003). There the resonance was found to be located at a magnetic field of 9.08(1) G, which disagrees with our result..

In order to determine the ratio of single to paired atoms in our three dimensional lattice potential, we monitor the loss of atoms when we hold the atomic sample at the resonance magnetic field for a variable time (see Fig. 3).

Lattice sites which are occupied by more than one atom are depleted within 3 ms due to the increased two- and three-body collision rates and the high density at single lattice sites. Isolated atoms, however, are protected from collisions and remain trapped. Assuming a negligible number of sites with three atoms, a fit with an exponential decay yields a time constant of $`1.3(2)`$ ms and a ratio $`N_1/N_21.1`$.

In order to determine the elastic scattering properties we apply the Ramsey interferometer sequence that has been described earlier (see Fig. 1). At each magnetic field value we record the Ramsey fringe visibility for different interaction times $`t_{\text{hold}}`$ (see e.g. Fig. 4 for a field of $`B`$=9.081 G).

The revival time $`t_R`$ at which the visibility shows its maximum is determined by a fit using eq. 1. This revival time depends on the difference in the interaction matrix element. A single revival of the fringe visibility could in principle also be caused by a complete loss of lattice sites with two atoms. We have however checked that even for the more pronounced losses at the Feshbach resonance the system exhibits dynamics due to interaction after the first revival. In this case losses would shift the revival time by less than 1%.

In order to extract information on the changes in the scattering length from the revival times, one can measure the the on-site matrix element $`U_{00}`$ through a collapse and revival experiment that we have demonstrated earlier Greiner et al. (2002b). For the same experimental parameters, we find $`U_{00}=h/396(11)\mu `$s. Using this information we can calculate $`\chi /U_{00}=\mathrm{\Delta }a_{s,\chi }/a_{00}`$, which expresses the change in scattering length measured in units of the scattering length $`a_{00}`$. In order to map out the change in the elastic scattering length on the Feshbach resonance, $`\mathrm{\Delta }a_{s,\chi }/a_{00}`$ has been measured for several magnetic fields and is shown in Fig. 5. Since the entanglement interferometer can only measure absolute values of $`\mathrm{\Delta }a_{s,\chi }/a_{00}`$, we perform a usual time of flight measurement to obtain information on the sign of the scattering length differences. For this, we leave the Feshbach field switched on for the first 3 ms of the TOF Inouye et al. (1998). During this time the altered interaction energy is converted into kinetic energy, and the size of the atom cloud is measured. For magnetic fields below the Feshbach resonance, we have found the cloud to be slightly larger (9% change in size) in the axial direction, whereas it is slightly smaller above the Feshbach resonance. From this we conclude that the interspecies scattering length grows and shrinks below and above the Feshbach resonance field respectively.

Since we ramp the magnetic field through the resonance in order to address the region above the resonance, we fit both branches of the scattering length change separately to a dispersive profile, with the constraint of a common center and width. The fit yields a center of 9.128(5) G and a width of 15(4) mG, in good agreement with our loss measurements. The two branches show an offset from the $`\chi =0`$ line which arises from the non-zero scattering length difference $`\delta a_s|a_{00}a_{11}|`$, which in principle can be extracted with high precision if decoherence permitted for longer interaction times.

The interferometric method presented allows for high precision measurements of relative changes of the scattering lengths. In order to demonstrate this, we assume an error-free scattering length $`a_{00}=100.4a_0`$ 333S. Kokkelmans, private communication., where $`a_0`$ is the Bohr radius. With this, we determine the change of elastic scattering length for $`B=`$ 9.081 G to be $`\mathrm{\Delta }a_{s,\chi }=4.2(1)a_0`$.

In conclusion, we have presented a novel interferometric method to create and investigate entanglement dynamics in an array of spin mixtures of neutral atoms. The observed entanglement oscillations allow the precise determination of interaction properties between atoms in different spin states. We have demonstrated the versatility of the interferometer by characterizing the elastic scattering properties of a newly predicted weak inter-state Feshbach resonance in 87Rb. We have found both the elastic and inelastic channel of the measured Feshbach resonance to be in good agreement with the theoretical prediction. The two-particle interferometer furthermore enables the direct creation of arrays of Bell-states together with the non-destructive separation of singly from doubly occupied sites.

We would like to thank Anton Scheich for support with the electronics and Tim Rom and Thorsten Best for assistance with the experiment. We also acknowledge financial support from the Bayerische Forschungsstiftung and AFOSR.


text:
# X-RAY BINARY POPULATIONS: THE LUMINOSITY FUNCTION OF NGC1569

## 1. INTRODUCTION

X-ray binaries (XRB) in nearby galaxies have been known since the Einstein era, but detections of significant samples were limited to just a few sources (e.g., Magellanic Clouds or M31) and interpretation of the source properties was often hampered by confusion problems. Chandra observations have revolutionized XRB studies with the discovery of large numbers of point X-ray sources in galaxies even beyond the Local Group (e.g., Fabbiano & White 2003). Short-term variability of many sources excludes the possibility of source confusion and strongly points toward accretion as the origin of the X-rays. The samples in most cases are large enough that we are now able to examine XRB populations in a wide range of galactic environments with different star formation histories. The samples are characterized by cumulative X-ray luminosity functions (XLFs) fitted by single or broken power-laws (e.g., Grimm et al. 2003; Zezas & Fabbiano 2002).

It has been noted that the XLF slopes follow a rather systematic behavior with population age and possibly star formation rate (SFR) (Grimm et al. 2003; Sarazin et al. 2003). Such behavior could constrain the SF history and properties of XRB populations in nearby galaxies. However, the development of reliable diagnostics requires a sound physical understanding of the observed correlations. With the exception of a recent study (Sipior et al. 2003, which does not focus on a specific set of observations), attempts to gain physical insight have been based so far on analytical models that assume the existence of power-law XLFs, and that X-ray lifetime is inversely proportional to X-ray luminosity and XRB properties do not depend on the SF history of the host galaxy (Wu 2001; Kilgard et al. 2002).

In order to develop a set of useful diagnostics, it is important to develop theoretical models for XRB formation and evolution that depend on the star-formation history of the host galaxies, and allow us to identify the main physical elements that determine the XLF slopes. Given that this study focuses on observed populations, it is clear that the employment of population synthesis models is necessary. First pre-Chandra XRB models for starbursts were developed and compared to ASCA observations of the total X-ray luminosity of WR galaxy He2-10 by Van Bever & Vanbeveren 2000. However, only now we may compare the specific theoretical models with the observed point source populations. A collaborative effort has led to the development of a detailed population synthesis code StarTrack (see § 2) that specifically focus on careful, self-consistent XRB calculations. Ultimately our goal is to construct a coherent picture of XRB formation and evolution based on theoretical models that have been first calibrated against observations of well-studied galaxies, and thus can be used in the interpretation of other XLF observations.

In this Letter we compare our model XLFs to Chandra observations of NGC 1569, a blue dwarf Irregular at a distance of 2.2Mpc (Israel 1988), characterized as a (post-)starburst galaxy. It has been selected as a good test case because its star formation history is relatively well constrained by HST optical and infrared observations (Aloisi et al. 2001; Greggio et al.1998; Vallenari & Bomans 1996) and has a long ($`80`$ ks) Chandra exposure providing a detection limit of $`10^{36}`$ erg s<sup>-1</sup> (Martin, Kobulnicky & Heckman 2002; hereafter M02). We are mainly concerned with two questions: (i) is it at all possible to theoretically reproduce the observed XLF? and (ii) do our models agree with the current constraints on the star-formation history of NGC 1569 derived by observations in other wavelenghts?

## 2. THEORETICAL MODELS

The StarTrack code was originally developed for the modeling of binaries with two compacts (BKB), but has recently undergone major revisions (Belczynski et al. 2003, in preparation) intended to treat in detail the formation and evolution of XRBs, for any choice of SF history and metallicity. The main revisions include a detailed treatment of tidal synchronization and circularization (Hut 1981), individual treatment of various mass-transfer (MT) episodes, full numerical orbit evolution with angular momentum losses due to magnetic breaking, gravitational radiation, mass transfer/loss, and tides.

We have calibrated the tidal implementation using observations of binary eccentricities and periods in stellar clusters (Mathieu et al. 1992) and of orbital decay in high-mass XRBs (Levine et al. 2000). Our MT implementation involves the detailed calculation of Roche-lobe overflow MT rates based on radius-mass exponents both for the donor stars and their Roche lobes. We have compared our results for a set of MT sequences to both published MT calculations (e.g., Beer & Podsiadlowski 2002) and results obtained within our group using an updated stellar evolution code (Ivanova et al. 2003) with very satisfactory agreement, much better than typical MT implementations. The modeled X-ray phases are identified as (i) stable, driven by nuclear evolution of the donor, magnetic braking or gravitational radiation (ii) thermally unstable with possibly anisotropic emission (King et al. 2001); (iii) Eddington-limited MT (cf. Kalogera et al. 2003, in preparation); and (iv) persistent or transient (critical MT rate from Dubus et al. 1999). We also account for wind accretion onto compact objects following Hurley, Tout & Pols (2002). The updated code has been already tested and used for a study of Galactic ultracompact binaries (Belczynski & Taam 2003).

In this paper we examine whether XRB models with rather standard binary-evolution parameters (i.e., not specifically selected for this study), but with star-formation history and metallicity consistent with what is known about NGC 1569 can produce an XLF shape in agreement with observations. We choose parameters from the reference model in BKB, with just a few differences: the maximum neutron star (NS) mass is set equal to $`2\mathrm{M}_{}`$, the most recently inferred natal NS kick distribution is incorporated (Arzoumanian, Chernoff & Cordes 2002), the primary masses are selected in the range $`4100\mathrm{M}_{}`$ (with an initial-mass-function slope of -2.7; see Kroupa, Tout & Gilmore 1993), and the secondary masses in the range $`0.08100\mathrm{M}_{}`$ with a flat mass ratio (secondary divided by primary mass) distribution.

We construct our XRB models with focus on the estimated metallicities and SFR properties of stellar populations in NGC 1569. A coherent picture has been developed (Aloisi et al. 2001; Greggio et al. 1998; Vallenari & Bomans 1996) for two populations: (i) one metal-rich ($`0.25`$ Z, although it could be comparable to 1 Z; M02), young population formed in a global burst of star formation that lasted about 100 Myr and seems to have stopped 5-10Myr ago, and (ii) one metal-poor ($`Z=0.0040.0004`$), old population formed by less vigorous (possibly by a factor of $`1020`$ in SFR), continuous star formation for the past 1-1.5 Gyr. For consistency with these estimates we adopt $`Z=0.005`$ and $`Z=0.0022`$, for the young and the old populations, respectively. We examine the consistency of our models by allowing the relative SFR weight and the population age and metallicity to vary.

## 3. THE X-RAY LUMINOSITY FUNCTION

We analyzed the archival Chandra data using the Ciao v3.0 data analysis suite. After screening for high background intervals we searched for sources in the full 0.3-7.0 keV band with the wavedetect source detection algorithm. Then, for each source within the D25 ellipse of the galaxy ($`\mathrm{D}_{\mathrm{max}}=3.6^{}`$, $`\mathrm{D}_{\mathrm{min}}=1.8^{}`$), we estimated the net number of counts using an aperture including all the emission from a source. The background was determined locally from a source-free annulus around each source. We detected 14 sources with a significance greater than 3$`\sigma `$ above the local background, in excellent agreement with M02, but we exclude one that has been identified as a supernova remnant by M02. Based on the LogN-LogS relation from the Champ survey (Kim et al. 2003) we estimate that at most two of our sources within the D25 area are associated with background or foreground objects.

The relatively high number of counts of the faintest sources ($`1520`$ counts) compared to their typical background ($`510`$ counts) and the lack of a turnover in the observed cumulative XLF (see Figure 2) indicate that incompleteness is not significant in the low luminosity end of the XLF. The higher level of diffuse emission in the central region of the galaxy (M02) may result to a higher detection threshold, but, given the very small number of sources in this region, we do not attempt to correct for this effect. In our XLF, apart from the standard Gehrels errors (Gehrels 1986), we include errors associated with the uncertainties in the observed count rate of each source, following Zezas & Fabbiano (2002). The luminosity of the sources is in the 0.1-10.0 keV band assuming a power-law spectrum ($`\mathrm{\Gamma }=1.7`$) with Galactic column density ($`2.3\times 10^{20}\mathrm{cm}^3`$) and is corrected for absorption.

In our theoretical models X-ray luminosities ($`L_X`$) are calculated based on the accretion rate $`\dot{M}`$ with appropriate efficiencies for NS and black holes (BH). At present no corrections specific to the Chandra energy band are taken into account, due to the lack of a reliable spectral model across a wide range of frequencies. We further apply the Eddington limit and calculate the effects of associated non-conservative mass transfer with the lost matter carrying away the specific angular momentum of the accretor.

To determine the contribution of transient systems to the XLF we need information on their X-ray duty cycle (DC), which however cannot be provided reliably by the disk instability theory. Empirically it is thought that $`DC1\%`$ (e.g., Taam, King & Ritter 2000). In StarTrack we adopt a DC value (in this study DC=1%) and, for each system we randomly sample the probability that the source is in outburst and then assign an X-ray luminosity equal to the Eddington limit, whereas for the systems in quiescence we assume that their X-ray luminosities are too low to be detectable. Donors more massive than the accretors can drive mass transfer on their thermal timescale (e.g., Kalogera & Webbink1996) and anisotropic emission is possible (King et al. 2001). We identify such systems based on the calculated radius-mass exponents. At present we do not account for the anisotropic-emission possibility, but do account for the possibility of super-Eddington transfer as described above. Both transient and thermal-timescale XRBs are Roche-lobe overflow (RLOF) systems, and typically dominate old ($`100`$ Myr) populations (see also Sipior et al. 2003).

Based on the current understanding of the SF history of NGC1569 (§ 2), we calculate combined XLFs for ages 1–1.5 Gyr and 70–110 Myr, for the old and young populations, respectively. The choice of the relative weight of the two populations also affects the XLF shape and we choose SFR weight factors in the range 10–50. Population models tend to have rather uncertain absolute normalizations and the absolute SFRs in NGC 1569 are not well constrained either (Vallenari & Bomans1996). Although, the total number of sources in a galaxy (corrected carefully for selection effects, partial galaxy coverage, etc) can provide additional model constraints, in this first step of our studies we restrict ourselves to calibrations and comparisons to observations based on the XLF shape. Therefore we used observed and model XLFs normalized to the total number of objects, and we take into account the small-number bias by randomly selecting from the simulated source populations samples of 13 sources at a time. The sampling of the underlaying model parent population is then repeated number of times ($`10^4`$) to yield the predicted XLFs and their associated errors.

Figure 1 indicates the behavior of XLF as a function of age of young XRB population of NGC1569, while the age of the old population is fixed at 1.5 Gyr. Unlike earlier suggestions (Wu 2001; Kilgard et al. 2002), we find that the dependence of the XLF slope on age is non-monotonic and this can be actually understood physically. At early times ($`10`$ Myr) the XLFs are by far dominated by MT-fed systems of the older population. The recent star formation episode barely begun producing wind-fed systems with the most massive compact objects (formed earlier) that accrete from strong winds from the next most massive donors. The combination of young and bright wind-fed sources with the MT-fed systems of older population tends to make XLFs rather flat. However, relatively quickly ($`70`$ Myr) the very massive stars of young population end their lives and are removed from the XRB sample. They are replaced by less massive, and fainter wind-fed systems and the XLF becomes steeper. As time passes ($`170`$ Myr) RLOF systems of younger population, that tend to be brighter than wind-fed sources of similar binary component masses, start becoming important in the total population and the XLF slope once again becomes flatter.

In Figure 2 we present XLFs calculated for NGC1569. Model A with parameters taken at face value from the HST observations: young population at 110 Myrs with SFR of 20 times larger than that of the old population at 1.5 Gyr. The calculated function tends to be flatter and extends to higher $`L_X`$ values compared to the observed XLF. At the age of 110 Myr since the onset of the burst (10 Myr after its end), the young population includes a significant number of a bright RLOF sources flattening the XLF. Model B corresponds to the young population at an earlier age of 70 Myrs (equal to the burst duration) when it is dominated by high mass XRBs, and we already see that the XLF becomes steeper, but not quite enough to match the observed one. There are still quite a few RLOF systems formed in the old star formation episode. However, the duration of that episode as well as its end time are not very precisely established. Model C corresponds to the young population at 70 Myrs (equal to the burst duration), but the old one at 1.3 Gyr with continuous star formation for 1 Gyr. The predicted XLF quite closely match the observed one within the relevant errors. We note the tendency for higher luminosities in the models, but we note that we adopt bolometric luminosities, which are higher than those in the Chandra band. In Table 1 we present the content of the X-ray binary population for our best model (C). The population is dominated by young high-mass XRBs, but with significant contribution of old RLOF systems. In the young population, where only the most massive stars have ended evolution, many accretors are BH, while in the older population sources with NS dominate.

## 4. DISCUSSION

We present our first results from XRB population models developed for comparison with current and future Chandra observations of nearby galaxies. We choose NGC 1569 as our first test case and we find good agreement between our models and the observations. This agreement is even more remarkable in view of the fact that we did not attempt to fine tune any of the model parameters related to X-ray binary evolution. However, we have explored other models with varying metallicities and IMF slopes, and found that both our quantitative and qualitative conclusions remain robust and the other models do not offer a better match to the observed XLF shape. This is true even when we account for the fraction of systems that can escape NGC 1569, due to systemic velocities acquired at supernova explosions.

Examination of various models with properties consistent with NGC 1569 constraints, lead us to conclusion that an age of 70 Myr for the young and 1.3 Gyr for the old population and a SFR relative weight of 40 are favored. In order to get agreement between the model and the observed XLF, we require a recent burst that is younger than inferred from the optical/infrared data. This slight discrepancy could be due to the fact that at this point we do not consider different black-hole binary spectral states and anisotropic emission from pulsar binaries. On the other hand the parameters of the older population (which is dominated by old non-magnetized neutron star binaries) are very consistent with the latest picture from the HST data (Angeretti et al. 2003, private communication).

We consider these encouraging results only a small, first step in our exploration of XRB models and their comparison to observations. As we gain experience with the study of specific galaxies, we expect to develop a reliable calibration system that will then allow us to extract information about origin of XRBs in other galaxies. A natural extension of this study will include two elements: the exploration of constraints on the absolute normalization of the XLF in addition to its shape, and the comparison models with a sample of starburst galaxies that form a time sequence with ages in a wide range of values to address the theoretical basis for correlations suggested by Grimm et al. (2003), for example. Detailed examination of degeneracies in the derived constraints is also important. Moreover, the modeling of supernovae remnants may prove to be necessary, since they may i) contribute significantly to and ii) be hard to remove from observed point source samples.

We thank the referee J. Irwin, and also A. King and T. Maccarone for useful comments, and the Aspen Center for Physics and the NU Visitors’ fund (to AZ) for support. We also thank M. Tosi and L. Angeretti for discussing their results prior to publication. This work is partially supported by a Packard Fellowship, and a Chandra theory grant to VK, NASA LTSA grant NAG5-13056 to AZ and VK and NASA grant NAS8-39073 to GF.